On the deep-mantle origin of the Deccan Traps

See allHide authors and affiliations

Science  10 Feb 2017:
Vol. 355, Issue 6325, pp. 613-616
DOI: 10.1126/science.aah4390

Double trouble for the Deccan Traps

The continental flood basalts in India known as the Deccan Traps formed from a massive outpouring of lava around the time that dinosaurs went extinct. The event dramatically reshaped the landscape and altered the climate. Glišović and Forte used time-reversed convection modeling to reconstruct the origin of this giant magmatic event. They found that two different deep mantle hotspots joined forces about 65 million years ago to produce one of the largest volcanic features on Earth.

Science, this issue p. 613


The Deccan Traps in west-central India constitute one of Earth’s largest continental flood basalt provinces, whose eruption played a role in the Cretaceous-Paleogene extinction event. The unknown mantle structure under the Indian Ocean at the start of the Cenozoic presents a challenge for connecting the event to a deep mantle origin. We used a back-and-forth iterative method for time-reversed convection modeling, which incorporates tomography-based, present-day mantle heterogeneity to reconstruct mantle structure at the start of the Cenozoic. We show a very low-density, deep-seated upwelling that ascends beneath the Réunion hot spot at the time of the Deccan eruptions. We found a second active upwelling below the Comores hot spot that likely contributed to the region of partial melt feeding the massive eruption.

Previous researchers associated the Deccan continental flood basalt (CFB) with the Réunion hot spot (1), a center of long-term magmatism presently located under the volcanic edifices of Réunion Island that appears to be fed by anomalously hot underlying mantle (2). The interpretation of the mantle dynamic processes that resulted in the eruption of the Deccan Traps has been mainly based on paleomagnetic (1, 3, 4), geochemical-geochronological (1, 47), and petrological (47) analyses of the surface geological record. Correlations between rapid changes in Indian plate velocity and time of the Deccan eruptions (3) along with theoretical modeling of plume-related plate-driving forces (8) support the longstanding hypothesis of a mantle plume origin for the Deccan CFB (9). Geological reconstructions of Indian Ocean paleotectonic activity (3, 8) provide surface constraints on possible plume sources for the Deccan CFB and subsequent Indian Ocean volcanism, but corresponding paleo-maps of internal Earth structure [similar to seismic maps of present-day mantle structure (10)] do not show an upwelling of hot material feeding the Deccan eruptions.

We present an alternative approach that links surface geological and geophysical processes to changes in mantle buoyancy determined by a back-and-forth iterative method for backward convection (11). A key aspect of this approach is the incorporation of present-day, three-dimensional (3D) structure of the mantle derived from the joint inversion of seismic, geodynamic, and mineral physics data in our convection model (1214). Seismic tomographic models contain a record of past states of mantle heterogeneity, allowing backward convection calculations to track the past evolution of this heterogeneity from the present-day. In particular, the ability to reconstruct mantle plumes depends on the adequate resolution of their remnant thermal structure in the upper mantle, as recorded in available tomography models. Our back-and-forth iterative solution of the time-reversed conservation of energy equation uses a quasi-reversible correction term for thermal diffusion (11), which is important in the thermal evolution of lithospheric-mantle structure.

We verified the reliability and geodynamic consistency of our convection model by comparing its global predictions of present-day convection-related surface signals—such as tectonic plate motions, gravity anomalies, and dynamic surface topography—against corresponding observations (table S1) (13). As a further check on the reliability of the reconstructed states of the mantle in the geologic past, we used them as starting conditions in convection modeling that we carried forward to the present day. We found that these reconstructions, when integrated forward in time from 70 million years ago (Ma) to the present day, yield high global correlations to density anomalies derived from the joint seismic-geodynamic tomography model at all depths in the mantle (11).

An important feature of our mantle convection model is the inclusion of a plate-like boundary condition (14) that incorporates the reconstructed geometry of the major tectonic plates over the Cenozoic (fig. S1) (15). Although we prescribed the plate geometries, their velocities are not. We predicted plate velocities at each instant from the evolving distribution of buoyancy forces in the mantle (11, 16). The uncertainties in the global plate geometry reconstructions augment rapidly with increasing time in the past, and we cannot easily quantify them. The global reconstructions of plate geometry and velocity nevertheless provide the only available constraints on the evolution of mantle buoyancy in the past. We therefore developed a mathematical inverse procedure that determines the smallest perturbations in the reconstructed density fields that yield a match to the geological estimates of past plate motions at selected instants (11). This procedure allows an objective mapping of those regions in the mantle that require a minimum “nudge” to the reconstructed buoyancy field in order to match the geologically determined plate motions.

2D cross sections of the mantle (Fig. 1, AA') show a deeply rooted hot upwelling, or plume, almost directly below the present-day location of the Réunion hot spot, whose activity appears well established at ~68 Ma (Fig. 1A). In the asthenosphere (near 220 km depth), lateral spreading of low-density material below Réunion emerges within a horizontal flow field, with a component directed toward Eurasia.

Fig. 1 A convection reconstruction of the Cenozoic evolution of mantle heterogeneity below the Indian Ocean.

(A to C) The mantle cross sections show the reconstruction of lateral density variations (11) at three different instants in the Cenozoic: (A) 67.5 Ma, (B) 50 Ma, and (C) 30 Ma. (D) The starting initial condition for backward convection modeling. These cross sections show the “nudged” pattern of density anomalies that yield a match to the geological inference of surface plate velocities at each instant (M models). Inset maps show the geographical orientation of the two cross sections, AA' and BB'; the reconstructed position of coastlines (black line) and plate boundaries (dark-green line) (25); and the present-day plate boundaries (green line) and continents (gray shading). I (blue rectangle), R1 (gray oval), R2 (green oval), R (violet oval), and C (light-blue oval) are indicators for the present-day mantle structure resolved by the joint tomographic model (D) that is convected backward beneath the position of India, Réunion, and Comores hot spots at 67.5 Ma (A), respectively. The distance between yellow dots on cross sections is 10°.

The model also resolves two distinct plume-like anomalies, 20° and 30° south of the Réunion plume (Fig. 1A), respectively. The first lies below the mid-ocean ridge that separates Madagascar and southern India, and the second is situated beneath the present-day Marion hot spot. The latter may be a remnant of the earlier upwelling that generated the Morondava CFB at ~93 Ma (8, 17). We also detected a large volume of subducting material under the reconstructed northern margin of the Neo-Tethys Ocean (Fig. 1A). The dynamical importance of these two anomalies, from a plate-driving perspective, becomes evident by comparing the “nudged” mantle reconstruction that matches the geological estimates of past plate velocities, the “M” models (Fig. 1), versus the direct reconstruction in which no matching to plate velocities is applied, the “N” models (fig. S2). Overall, the differences between N and M models are minimal, although we found that upwelling and downwelling under the “residual” Morondava plume and the northern margin of the Neo-Tethys Ocean, respectively, are more strongly expressed in the M model (fig. S2). We thus found that large-scale subduction in the north and two mantle upwellings in the south (the Réunion and residual Morondava plumes) operated simultaneously, generating a large mantle “roll” that propelled India northward during the Deccan eruptions.

The West-East mantle cross sections (Fig. 1, BB′) show penetration of hot, low-density mantle at sublithospheric depth beneath the Réunion hot spot at ~68 Ma. The N (no plate-velocity matching) model shows that sublithospheric injection of hot buoyant material by the Réunion plume propagated toward the paleo-Carlsberg ridge and the paleo-Owen Fracture Zone, aiding the divergence of India from Africa (fig. S2A). At ~68 Ma, we found a second lesser (or less resolved) mantle upwelling ascending below the present-day location of the Comores hot spot (18, 19) that injects buoyant material into the same ridge and fracture zone. The M model indicates that lithospheric injection by the Comores upwelling is channeled eastward toward the Réunion plume (Fig. 1A). These mantle reconstructions thus reveal the possible contribution of an unrecognized plume below the Comores hot spot as an additional source for the Deccan CFB, particularly their northern portion. Previous hypotheses for a “non-plume” origin of the Deccan CFB (20) are not supported by our mantle reconstructions (Fig. 1A).

The time-dependent viscous drag generated by the Réunion upwelling has effectively lifted a slab-like mantle anomaly, which was previously under northern India at ~68 Ma, to the top of the mantle. This surprising vertical displacement occurs relatively rapidly, so that between 55 and 50 Ma (after India has completely passed overhead), the top of the slab is hoisted up toward the lithosphere, at which point it starts to sink back into the mantle (Fig. 1B). This quasi-subduction, aft of the advancing Indian plate, generates a local upper-mantle roll that further attenuates the push exerted by the remnant Réunion plume.

The reconstructed distribution of buoyancy at 30 Ma (Fig. 1C) reveals substantially weakened Comores and Réunion plumes and correspondingly reduced northward velocity of India (especially evident in comparing reconstructions at ~68 and 30 Ma).

Several factors contribute to uncertainties in the backward convection simulations (Fig. 1) and to the growth of these uncertainties with increasing time (11, 21). One major source of uncertainty is the starting [time (t) = 0] tomography model (12). A quantitative assessment of the uncertainties in the starting mantle heterogeneity model is difficult owing to the intrinsic nonuniqueness of the tomographic inversions, which is further compounded by the application of damping and smoothing. We therefore tested a second global tomography model, independently derived by inverting seismic data alone (fig. S4), in backward convection simulations to reconstruct the Cenozoic evolution of the mantle (fig. S5). Despite the differences between the tomography models (fig. S4), especially with regard to their ability to fit geodynamic data (table S1), the backward convection procedure again reconstructs a deep-seated, anomalously hot upwelling under the Réunion hot spot at ~68 Ma (fig. S5A).

The choice of mathematical parametrization of lateral heterogeneity (for example, discrete cells or truncated series of basis functions) also contributes to uncertainty. These uncertainties in the starting tomography model are amplified as the backward convection calculations proceed further into the past (21). This is evident in Fig. 1, in which the “bumpy” or staccato structure of the reconstructed mantle heterogeneity is a consequence of the backward propagation of the discrete cells and vertical layers used to parametrize the starting tomography model (12). Although this imperfect short-scale resolution results in somewhat irregular shapes for the reconstructed upwellings, the associated flow patterns are more robust because they are the result of buoyancy integrated over much larger length scales (13).

The mantle viscosity structure (13) used in the time-reversed convection simulations (Fig. 1) is another source of uncertainty. The predicted mantle flow speeds are directly dependent on the absolute value of viscosity at different depths. The time-reversed evolution of mantle heterogeneity is a nonlinear function of these flow speeds, in which the nonlinearity grows with time. We extensively tested the viscosity profile we used in our backward reconstructions (11) by verifying that it yields good fits to a suite of geodynamic data sets that are sensitive to relative changes in viscosity across the mantle (global gravity and topography anomalies) and to absolute values of viscosity (plate motions and glacial isostatic adjustment) (13). Another key test was to verify that the predicted timing of surface geological phenomena—such as dynamic topography (22), synchronicity between paleo-ridge location and mantle buoyancy (14), and melt production (23)—agrees with available data.

The dominant features in the viscous flow prediction of dynamic surface topography at ~68 Ma are the strong surface uplift along the western coastal margin of India, with highest topography located between the Comores and Réunion hot spots (~3 km), and a deep subduction-related depression along the northern margin of the Neo-Tethys Ocean (Fig. 2A). Furthermore, this predicted uplift is compatible with the height of the Seychelles Plateau relative to the surrounding seafloor, the geochemical evidence for subaerial exposure of the offshore reflector sequences, as well as the emplacement of the Silhouette and North islands during the Deccan eruptions (4).

The pattern and amplitude of uplifted topography remain approximately constant between ~68 and 60 Ma; however, its position relative to India changed owing to the northeastward motion of the Indian plate (fig. S6). The resulting southward progression of the eastward tilt of the Indian coast provides a favorable orientation for the southerly dip and age progression of the successive Deccan lava flows (24).

Fig. 2 Dynamic surface topography, lithospheric structure, and mantle melt production.

(A) Dynamic surface topography predicted at 67.5 Ma. (B) Reconstructed lateral variations in density at 80 km depth at 67.5 Ma. (C) Map of mantle melting at 80 km depth at 67.5 Ma. The solid red contour line represents the region in which the melting temperature is exceeded by 300 K. The pink shading represents the region in which the melting temperature is exceeded by 100 K. (Inset) The depth-variation of total melt area, where melting temperature is exceeded by 300 K. The proto-Réunion hot spot (violet circle) is merely a location within the region of maximum melting that has the same latitude as the Réunion hot spot (magenta circle). (D) Total volume of mantle melting, as a function of time, predicted below each of the hot spots. The curves are color-coded with reference to the colored disks representing hot spots in map (C). These volumes are calculated by integrating all areas that undergo melting (where the melting temperature is exceeded by 300 K) within a 10° radius disk centered on each hot spot. The black dashed line shows the melt volume sampled by a 10° disk that moves northward with India and is exactly coincident with the Réunion hot spot at 67.5 Ma. In the maps shown in (A) to (C), the dark-green line shows the reconstructed plate boundaries (25); the green line represents the present-day plate boundaries; the black line represents the reconstructed position of coastlines (25); and the gray line shows the present-day coastlines.

The structure at the base of the lithosphere at ~68 Ma (Fig. 2B) shows a very good match between our time-reversed reconstruction of high-temperature, buoyant mantle along the western and southern margin of India and the independent geologic reconstruction (25) of the paleo-Carlsberg ridge location. This fit or synchronicity (14) is a further test of the good timing provided by the mantle viscosity used in our reconstructions.

To determine locations of mantle melting, we superimposed the reconstructed lateral temperature variations on our reference temperature profile that satisfies geothermal and mineral physics constraints (11) and mapped out all regions that exceed the solidus temperature of dry peridotite (26). We predict that the highest-temperature regions that undergo melting (Fig. 2C) are most widespread between ~50 and ~80 km and are located beneath the west coast of India, where they underlie regions of uplift in dynamic topography. These findings agree well with petrological inferences that the Deccan magmas were produced by high-temperature melting in a depth range between 60 and 100 km (23). The volume of melted mantle that was potentially accessible to the Réunion and Comores hot spots was ~40 and ~35 million km3, respectively, at ~68 Ma, decreasing to near zero around 20 and 40 Ma, respectively (Fig. 2D). Such ample amounts of mantle melting would have been sufficient to supply the Deccan lava flows as well as the subsequent eruptions that created the Chagos-Laccadive ridge and Mascarene plateau.

We found that the potentially accessible mantle melt volume, sampled beneath a point on the Indian continent that is coincident with the Réunion hot spot at 67.5 Ma, should have reached a maximum at about 67 to 68 Ma and subsequently should have decreased to zero at ~60 Ma as India moved northward (Fig. 2D). Therefore, the rapid northward advance of India implies that the volume of melted mantle that is potentially accessible to the Deccan lava flows, from the source region under the Réunion hot spot or anywhere else within the predicted melting region (Fig. 2C), should rapidly decrease with time. These considerations clearly do not make allowance for the influence of the near-surface magmatic plumbing system that delivers the melt produced in the mantle to the surface, nor for external influences such as the recently proposed impact-induced enhancement of the Deccan flows (27). Thermomechanical modeling of partial melting, magma transport from depth to the surface, and pressure-driven fracturing of the brittle crust (28) exceeds the scope of our current mantle convection models. Although our model cannot exactly predict the peak eruption of the Deccan lava flows at 66 ± 0.5 Ma, as shown in (6, 27), it does accurately predict a ~10-million-year time window that strongly peaked near 68 Ma, in which mantle melt was accessible to lava flows that preceded, coincided, and then followed the main Deccan event (29).

Our back-and-forth iterative reconstructions of the evolution of mantle heterogeneity (11) show that the rapid northward advance of India during the Cenozoic was entrained by a mantle-wide convective roll driven by upwelling beneath the Réunion hot spot and a large-scale downwelling beneath the Tethyan convergence zone. In the late-Cenozoic, the positive buoyancy associated with the initial (~68 Ma) Réunion plume has been advected northward and to shallower depths and continues to drive a present-day “conveyor belt” (30) that pushes the Indian plate northward. The Réunion upwelling thus plays a major role in the evolution of Indian Ocean tectonics, providing the buoyancy and mantle melting that caused the breakup of the Seychelles microcontinent from India and the associated Deccan lava flows. The remnant buoyancy from this upwelling, although substantially diminished, continues to be active today in the upper mantle between Réunion and the southern margin of the Indian plate (Fig. 1D).

Correction (13 February 2017): The Acknowledgments were corrected to include “A.M.F. is now at the University of Florida (UF) and he acknowledges the additional support for this work provided by UF.”

Supplementary Materials

Materials and Methods

Table S1

Figs. S1 to S6

Reference (31)

References and Notes

  1. Materials and methods are available as supplementary materials.
  2. Acknowledgments: The authors acknowledge the support of this work provided by the Natural Sciences and Engineering Research Council of Canada. A.M.F. is now at the University of Florida (UF) and he acknowledges the additional support for this work provided by UF. We thank D. Rowley for his reconstructions of the Cenozoic plate boundaries and for the discussions regarding the plate kinematic constraints over this time period. We also thank S. Zahirovic for providing us with plate configuration files generated with Gplates software. The authors are grateful to N. Simmons and S. Grand for the collaboration that resulted in the global tomography model used in this study. Additional thanks to N. Simmons for the schematic in fig. S3. The convection simulations in this study were carried out thanks to supercomputing facilities of Calcul Québec consortia at Université de Montréal. P.G. and A.M.F. both developed the scientific concepts and modeling techniques presented here and jointly wrote the paper. P.G. created the figures. All data needed to evaluate the conclusions are presented in the paper and supplementary materials. Additional data related to this paper may be requested from the authors.
View Abstract

Stay Connected to Science

Navigate This Article