Numerical Models Of Mantle Convection Research Articles

7th Workshop on “Numerical Modelling of Mantle Convection and Lithospheric Dynamics”, Aussois, France, June 2-6, 2001

7th Workshop on “Numerical Modelling of Mantle Convection and Lithospheric Dynamics”, Aussois, France, June 2-6, 2001

Degree-1 mantle convection and the crustal dichotomy on Mars

The surface of Mars consists of an old, heavily cratered and elevated southern hemisphere and younger, resurfaced and depressed northern hemisphere, a feature often termed the crustal dichotomy. The global crustal structure [Zuber et al., 2000] revealed by topography and gravity data from the Mars Global Surveyor spacecraft, and the possible late formation of the boundary zone between the hemispheres [McGill and Dimitriou, 1990], have been proposed to indicate an endogenic origin for the dichotomy. However, degree-1 mantle convection that is required for any endogenic process to be viable cannot be produced with conventional mantle convection models [Schubert et al., 1990]. We have studied the role of radially stratified viscosity on mantle deformation by using Rayleigh–Taylor instability analyses in a spherical shell geometry. Our analyses reveal that when mantle viscosity is stratified with a weak asthenosphere, deformation at long wavelengths is more efficient than that at short wavelengths. The weaker the asthenosphere, the longer the wavelength at which the deformation is the most efficient. A thicker asthenosphere also favors the deformation at long wavelengths. Both the Rayleigh–Taylor instability analyses and numerical modeling of mantle convection show that degree-1 convection can be produced within the Martian mantle provided that the mantle had a weak asthenosphere (∼500 km thick and ∼10 2 times weaker than the underlying mantle) early in planetary history. The degree-1 convection causes preferential heating of one hemisphere that may explain the primary features associated with the dichotomy in crustal structure.

Heat flow scaling for mantle convection below a conducting lid: Resolving seemingly inconsistent modeling results regarding continental heat flow

Numerical models of mantle convection below a conducting lid, meant to mimic a continent, have been used 1) to argue that heat flow variations in stable continental regions result principally from lithospheric thickness variations, and 2) to argue that the relationship between lithospheric thickness and continental heat flow is too weak to allow for this. We reconcile these results using a theoretical heat flow scaling which shows that the relationship between heat flow and lid thickness can take on two end‐member forms. As the ratio of lid thickness to convecting layer depth increases, the effects of lid thickness variations on heat flow move from weak to strong. Modeling studies that concluded that lithospheric thickness variations could lead to significant continental heat flow variations assumed a convecting layer depth appropriate to upper mantle convection and, thus, maximized the relative thickness of the continental lithosphere, i.e., the conducting lid. Studies that concluded that the relationship is weak assumed whole mantle convection which placed them on the flatter end of the heat flow versus relative lid thickness curve.

The Number of Hotspots in Mantle Convection: Effect of Depth-Dependent Viscosity and Internal Heating in Two-Dimensional Models

Two-dimensional numerical models of mantle convection have been calculated for a range of high Rayleigh numbers, depth-dependent viscosity and basal plus internal heating. Large aspect ratio boxes have been used in the calculations in order to estimate the expected areal density of upwellings in infinite fluid layers. The results are analyzed with regard to the number of the Earth’s hotspots which are assumed to be surface imprints of cylindrical upwellings in the mantle. For a pure whole-mantle situation, 6–7 upwellings can be expected. If the upper mantle convects separately above the 660 km discontinuity (allowing a second convective layer below 660 km depth), the theoretically estimated number of upper-mantle plumes can be as high as 250. Given the number of real hotspots (42 to 117 according to different compilations), it is suggested that the flow regime of the mantle is intermediate between the pure whole-mantle or pure two-layer circulation.

Numerical models of a dense layer at the base of the mantle and implications for the geodynamics of D″

To investigate the dynamics of the mantle's D″ layer, we explore numerical models of mantle convection which feature a dense basal boundary layer. We use a double‐diffusive finite element convection scheme and vary thermal and chemical Rayleigh numbers and properties including viscosity, thermal conductivity, and internal heating. For isoviscous models with heating only from below, the thermal Rayleigh number (Ra) is set at either 106 or 107. The negative chemical buoyancy of a dense layer produces a stable boundary layer when the ratio of chemical to thermal buoyancy (the buoyancy number B) is around 1. For B = 0.5 the layer is unstable, while B = 0.6 may produce a layer which remains stable for long periods, depending on other factors (such as Rayleigh number, layer thickness, and thermal diffusivity of the layer). At high enough Ra, small‐scale convection can occur within the layer. Using Ra = 107, we look at changes in layer thickness from 100 to 300 km at two values of B(0.6 and 1). Convection within the layer takes place most easily when the layer has a greater initial thickness. For B = 0.6, the layer is not always continuous along the lower boundary, so convection within the layer does not occur unless the initial thickness is fairly high (300 km). Increasing the thermal diffusivity of the layer (to simulate enrichment in metals) enhances heat conduction across the layer and can suppress internal convection within D″. Enhanced conduction also leads to higher plume temperatures. Including internal heating mainly increases the overall heat flow through the mantle, leading to higher surface heat flux. Finally, we examine models with temperature‐dependent viscosity and pressure‐dependent viscosity using Rayleigh numbers in the range of 107. Convection within the dense layer is enhanced by the relative reduction in viscosity. Our models exhibit only a modest decrease in viscosity (an order of magnitude) but illustrate how a dense, low‐viscosity layer interacts with cold, viscous downwellings.

The COSY Project: verification of global seismic modeling algorithms

Progress in determining the details of the global 3-D seismic velocity structure requires the ability to accurately model seismic wave propagation (e.g., travel times, waveforms, etc.) through heterogeneous 3-D Earth models. While for spherically symmetric models (quasi-) analytical solutions are available for the verification of numerical algorithms, this is not the case for general heterogeneous models. It is therefore desirable to establish global 3-D test models and verified reference seismograms, which allow us to assess the accuracy of numerical algorithms quantitatively. Prior to a workshop held on this issue at the 1997 IASPEI Meeting, a 3-D test model was handed out to various groups and long-period synthetic seismograms were returned. This workshop was the initiation of the COmparison of global SYnthetic seismogram techniques (COSY) Project, which aims at establishing a WWW page (http://www.geophysik.uni-muenchen.de/COSY), where the test models and seismograms as well as some of the algorithms can be accessed. In this paper, we study the accuracy of and compare solutions from different numerical methods for a spherically symmetric model and the 3-D test model. The algorithms compared use the normal-mode method, the Direct Solution Method (DSM), a direct evaluation of the Greens function for spherically symmetric media (GEMINI), and the finite-difference (FD) method. Our 3-D test model is a perturbation to the spherically symmetric background model (PREM) based on a (scaled) temperature field from numerical modeling of mantle convection. The model displays many features in common with recent seismic tomographic images. We suggest that in addition to (future) 3-D reference Earth models, verified reference synthetic seismograms should be established for use by the seismological community.

Plate tectonics and convection in the Earth's mantle: toward a numerical simulation

Numerical models of mantle convection are starting to reproduce many of the essential features of continental drift and plate tectonics. The authors show how such methods can integrate a wide variety of geophysical and geological observations. The goal is to combine the Stokes and energy equations with a realistic rheology, thereby letting us understand the complex dynamic coupling that occurs in the mantle and that gives rise to plate tectonics and other surface features. This approach holds great promise because it makes a tremendous amount of data relevant to understanding Earth's dynamics. The challenge is that the computational models must be inherently realistic, particularly when predicting observed geography or plate history, so that the models can be connected with observations. The authors view this as one of the most exciting future directions of computational geodynamics.

The evolution of continental roots in numerical thermo-chemical mantle convection models including differentiation by partial melting

Incorporating upper mantle differentiation through decompression melting in a numerical mantle convection model, we demonstrate that a compositionally distinct root consisting of depleted peridotite can grow and remain stable during a long period of secular cooling. Our modeling results show that in a hot convecting mantle partial melting will produce a compositional layering in a relatively short time of about 50 Ma. Due to secular cooling mantle differentiation finally stops before 1 Ga. The resulting continental root remains stable on a billion year time scale due to the combined effects of its intrinsically lower density and temperature-dependent rheology. Two different parameterizations of the melting phase-diagram are used in the models. The results indicate that during the Archaean melting occurred on a significant scale in the deep regions of the upper mantle, at pressures in excess of 15 GPa. The compositional depths of continental roots extend to 400 km depending on the potential temperature and the type of phase-diagram parameterization used in the model. The results reveal a strong correlation between lateral variations of temperature and the thickness of the continental root. This shows that cold regions in cratons are stabilized by a thick depleted root.

The influence of plate kinematics, convection intensity, and subduction geometry on the Earth’s upper-mantle dynamics in the vicinity of a subduction zone

SUMMARY This study reports the results of a 2-D numerical mantle convection model, featuring plate kinematics and non-vertical subduction. The main question we address is the combined eVects of plate velocities, convection intensity, and subduction dip angle on the Earth’s upper-mantle dynamics in the vicinity of a subduction zone. Through a systematic investigation using Rayleigh numbers in the range 1.2◊104‐6◊105 ,a motionless overriding plate, subduction dip angles of 70° ,4 5° and 30°, and subduction velocities in the range 0.3‐10 cm yr’1, we aim to study the interactions of the subducting and overriding plates with the underlying mantle. We focus on the possible transitions of flow structures, from a multicellular to a monocellular circulation, and from an unsteady to a steady global structure. Unsteady states can be either periodical or nonperiodical. The model predicts that unsteady structures are characterized by instabilities that aVect both horizontal thermal boundary layers under the moving plate, but only the upper thermal boundary layer under the motionless overriding plate. The movement of the subducting plate may be suYcient to drive the lower part of the circulation under the motionless overriding plate, characterized by an elongation of cells. This latter influence is favoured by a decreasing subduction dip angle. Thus diVerent transition relations for the flow structure are obtained according to whether we consider the moving plate or the motionless overriding one, and whether we consider a shallowor a steep-dipping subduction. By focusing first on the dynamics under the moving plate, we are better able to quantify the relationship between the buoyancy and boundary forces associated with the subducting lithosphere. For steep to intermediate subduction angles, we find, as in previous studies with no down-going slab boundary, that the greater the convection intensity, the greater the plate velocity that is required to achieve the dynamical mantle‐plate coupling. This is not, however, characteristic of a shallowdipping subduction, since we show that the down-going plate has a stabilizing influence, suggesting the existence of a critical subduction dip angle for the mantle‐plate coupling. Thus the investigation of internal dynamics under the motionless overriding plate reveals a surprising single-cell circulation under this plate for a suYciently high convection intensity. This single-cell circulation is obtained for smaller subducting plate velocities since the value of the subduction dip angle decreases. The diVerent lateral convection styles under this plate at moderate subduction velocities, namely multicellular and monocellular for dip angles of 70° and 30°, respectively, are consistent with stress states in the overriding lithosphere. Indeed, steep- and shallow-dipping subductions are likely to be associated with back-arc extension and compression, respectively. For a dip angle of 70° under the overriding plate, the model predicts an upwelling next to the trench, the location of which is consistent with observed back-arc basins.

The Role of Subduction‐Induced Subsidence in the Evolution of the Karoo Basin

Abstract The stratigraphic record indicates that during the Late Carboniferous to Early Triassic, South Africa experienced an episode of large‐scale subsidence that resulted in the development of the Karoo Basin. The southern margin of the region was contemporaneously dominated by an episode of compressional tectonics, and there is evidence for subduction of the paleo‐Pacific plate beneath the Karoo Basin. We propose that the long‐wavelength component of the basin subsidence resulted from the deflection of the lithosphere due to mantle flow coupled to the adjacent subduction. This was followed by a recovery phase during which the entire region was uplifted on cessation of subduction (and subsequent loss of dynamic support for subsidence), contributing to the present‐day high‐elevation topography of South Africa. Numerical models of mantle convection are presented and show that a shallow‐dipping slab is able to reconcile the observed long‐wavelength subsidence of the basin. A significant component of the n...

EMPLACEMENT OF DEEP MANTLE ROCKS INTO CRATONIC LITHOSPHERE BY CONVECTION AND DIAPIRIC UPWELLING

EMPLACEMENT OF DEEP MANTLE ROCKS INTO CRATONIC LITHOSPHERE BY CONVECTION AND DIAPIRIC UPWELLING

Mantle convection simulations with rheologies that generate plate-like behaviour

A long-standing problem in geodynamics is how to incorporate surface plates in numerical models of mantle convection. Plates have usually been inserted explicitly in convection models as rigidrafts1,2,3, as a separate rheological layer4,5 or as a high-viscosity region within weak zones6,7,8,9,10,11. Plates have also been generated intrinsically through the use of a more complex (non-newtonian) rheology for the entire model12,13 but with a prescribed mantle flow. However, previous attempts to generate plates intrinsically and in a self-consistent manner (without prescribed flow) have not produced surface motions that appear plate-like14,15,16. Here we present a three-dimensional convection model that generates plates in a self-consistent manner through the use of a rheology that is temperature and strain-rate dependent, and which incorporates the concept of a yield stress. This rheology induces a stiff layer on top of a convecting fluid, and we find that this layer breaks at sufficiently high stresses. The model produces a style of convection that contains some of the important features of plate tectonics, such as the subduction of the stiff layer and plate-like motion on the surface of the fluid mantle. However, the model also produces some non-Earth-like features, such as episodic subduction followed by the slow growth of a new stiff layer, which may be more consistent with the style of convection found on Venus.

Slab weakening by the exothermic olivine‐spinel phase change

Numerical modelling of mantle convection, with lateral variations of viscosity, and both the 670 km depth endothermic and the 400 km exothermic phase changes, reveals that downwellings are weakened upon crossing the exothermic phase change. A high temperature/low viscosity section is developed within the downgoing current, owing to both the olivine to spinel phase transition and the slowing down of the flow due to the 670 km endothermic transition. This mechanism may cause the creation of a weak zone at 400 km that results in a preferential location for horizontal or vertical slab tearing as revealed by seismic gaps in many slabs.

Evolution of starting mantle plumes: a comparison between numerical and laboratory models

Laboratory experiments on the initiation and development of mantle plumes suggest the formation of a bulbous plume head that grows in time by entrainment of surrounding material, followed by a thin feeder conduit. This type of plume is generally not seen in numerical models of mantle convection. To better determine the possible differences between numerical and laboratory experiments, we have numerically reproduced a published laboratory experiment on plume initiation. Theoretical predictions for the growth of the plume head have been verified and the experiments allowed for a new determination of a similarity constant. Models with a constant viscosity fluid show much less entrainment. Simulations using a plume rising from a thermal boundary layer with an olivine-type temperature-dependence show that most of the source material is contained in the plume head and that the conduit is largely composed of entrained material. Simulations with a more realistic non-Newtonian rheology show that plumes can rise much faster through the mantle, which will lead to higher temperatures of the plume head at the Earth's surface compared to what experiments with Newtonian rheology suggest.

Numerical models of mantle convection with secular cooling

SUMMARY A 2-D time-dependent finite-difference numerical model is used to investigate the thermal character and evolution of a convecting layer which is cooling as it convects. Two basic cooling modes are considered: in the first, both upper and lower boundaries are cooled at the same rate, while maintaining the same temperature difference across the layer; in the second, the lower boundary temperature decreases with time while the upper boundary temperature is fixed at 0 C. The first cooling mode simulates the effects of internal heating while the second simulates planetary cooling as mantle convection extracts heat from, and thereby cools, the Earth's core. The mathematical analogue between the effects of cooling and internal heating is verified for finiteamplitude convection. It is found that after an initial transient period the central core of a steady but vigorous convection cell cools at a constant rate which is governed by the rate of cooling of the boundaries and the viscosity structure of the layer. For uppermantle models the transient stage lasts for about 30 per cent of the age of the Earth, while for the whole mantle it lasts for longer than the age of the Earth. Consequently, in our models the bulk cooling of the mantle lags behind the cooling of the coremantle boundary. Models with temperature-dependent viscosity are found to cool in the same manner as models with depth-dependent viscosity; the rate of cooling is controlled primarily by the horizontally averaged variation of viscosity with depth. If the Earth's mantle cools in a similar fashion, secular cooling of the planet may be insensitive to lateral variations of viscosity.

A model for the layered upper mantle

Numerical modeling of mantle convection by Liu (1994, Science, 264: 1904–1907) favors a two-layer convection, if the results are reinterpreted for the correct phase relations in (Mg,Fe) 2SiO 4. The resulting chemical isolation of the upper and lower mantle suggests a highly differentiated and layered upper mantle to account for the discrepancy between the observed compositions of mantle xenoliths and the cosmic abundances of elements. It is shown that a layered upper mantle with a hidden reservoir can have a structure consistent with the observed seismic velocity profiles and an average bulk composition corresponding to the cosmic abundances. The evolution of the upper mantle and the origin of komatiites are discussed in the context of the proposed model.

Continental collisions in wide aspect ratio and high Rayleigh number two‐dimensional mantle convection models

Distinct rigidly moving oceanic and continental plates of finite thickness are incorporated into a two‐dimensional numerical model of mantle convection. We investigate upper mantle convection in models having aspect ratios as great as 24 and compare our findings with the results of earlier studies which were limited to aspect ratio 4 models. In addition, we implement models of whole mantle flow by specifying high Rayleigh number convection and thinner nondimensional plates. We are thus able to compare the results of continental collision models which include similarly sized continents in the cases of upper and whole mantle convection. For each case considered we model a pair of identical continents being carried toward a site of plate convergence by underlying counterrotating mantle convection cells. Upon collision, the continents form a motionless, rigid, model supercontinent, while oceanic plate material continues to recycle through the mantle. Following the continental collision, our models of upper mantle convection exhibit a reorganization of the convective planform below the model supercontinent into a smaller wavelength mode which is unable to generate the net stress needed to break apart the continent; alternating compressive and tensile deviatoric stress associated with the small scale flow results in a low integrated stress. In contrast the large scale of whole mantle convection enables flow reversals to produce shear stresses acting in a common direction over extensive areas of the base of a continent, the integrated effect of which is capable of causing continental rifting. The conventional view of the role of thermal blanketing in continental rifting does not apply in the whole mantle convection scenario.

Comparison between tomographic structures and models of convection in the upper mantle

SUMMARY Gruneisen's and third-order finite-strain theories are used to compute the density and seismic-wave velocities of minerals. Assuming the minealogical model of Ito & Takahashi (1987), seismic velocities of the upper mantle are calculated using the Hashin–Shtrikman averaging procedure. 1-D profiles are first obtained along adiabats, and compared to the IASP91 model. Different adiabats are considered in order to take into account the thermal effect of phase transitions. The best results are found with adiabats initiated at 1473, 1573 and 1613 K for α-olivine, β- and γ-spinel, respectively. The incorporation of thermal effects resulting from phase transitions gives velocity jumps at discontinuities close to those of IASP91. Next, a model of convection constructed by Dupeyrat, Sotin & Parmentier (1995), incorporating plate tectonics, is used to compute 1-D profiles and 2-D fields of seismic anomalies in the upper mantle. Averaged profiles show seismic-velocity gradients very close to those of IASP91, but individual values are much too high, suggesting that the mean temperature profile of the convection model is too cold by 400 K. When low-pass filtered to the resolution scale of presently available tomogrphic models, both the amplitude and shape of the computed seismic anomalies are consistent with the results of tomographic studies. The amplitude of the anomalies ranges between – 2.7 and 3.8% for P-wave slownesses, and from – 3.3 to 4.5% for S-wave slownesses. These anomalies correspond to lateral temperature variations of –465 to 520 K. These calculations are used (1) as an aid to the interpretation of global tomographic models, for instance by computing spectra of lateral heterogeneities, and (2) to test the adequacy of the basic assumptions used in the computation of numerical models of mantle convection, and to build a theoretical temperature profile that would give the best fit to IASP91. In the uppermost mantle this theoretical model has a shape close to both the convection model and the 1473 K adiabat, but in the transition zone the profile is highly subadiabatic. The spectra obtained for the synthetic seismic anomalies resemble that of tomographic studies, with most of the energy contained within gravest angular orders l, and a fast decrease of energy with increasing l. The spatial filtering has clearly different effects on heterogeneities, depending on their respective wavelengths. It is suggested that the change of decreasing rate observed in tomographic models at l= 7 is closely related to the filter wavelength and may correspond at a lesser extent to a characteristic wavelength to mantle heterogeneities.

The period of the free core nutation: towards a dynamical basis for an ‘extra-flattening’ of the core-mantle boundary

Indirect observations and theoretical predictions for the period of the free core nutation (FCN) differ by anywhere from 15 to 30 days, and various effects have been invoked in attempts to explain this difference. The favored explanation remains as much as 5% departure in the flattening of the core-mantle boundary (CMB) from that of its hydrostatic reference figure. This 5% ‘extra-flattening’ of the CMB is not seen at the Earth's surface, where the difference is only about 0.5%. In contrast to the a posteriori model adjustments used to determine this up to 5% value, and the kinematic results available from viscous flow modeling using the seismically determined lateral heterogeneity in density data, we consider this problem from the perspective of a forward-modeling dynamical study. More specifically, we investigate the related problem of flow-induced surface and CMB topography, arising from convection in the mantle. As such, we have completed a comparative and systematic study of relative surface and CMB topography resulting from numerical models of mantle convection. When effects resulting from boundary curvature are isolated, it appears that the magnitude of CMB topography produced is insufficient in producing a significant extra-flattening of the CMB. However, results concerning effects solely resulting from a depth-dependent mantle viscosity profile, indicate that this factor may indeed lead to enhanced topography at the CMB of the magnitude required to produce the extra-flattening there.

Penetration of plates and plumes through the mantle transition zone

Reliable evaluations of whether a phase transformation boundary with a negative Clapeyron slope will cause mantle convection to be stratified will require models to include good simulations of plates and plumes. Here numerical models of mantle convection are presented that incorporate reasonable simulations of subducting lithospheric plates and good simulations of rising low-viscosity mantle plumes, as well as being carefully scaled to other first-order parameters of the mantle. The models show that subducting plates and plume heads penetrate a phase transformation barrier more readily than does flow in constant viscosity convection models in which plates and plumes are poorly simulated. Plates can penetrate a larger magnitude of the Clapeyron slope than plume heads, with plume tails the least able to penetrate. These results, together with recent laboratory results on Clapeyron slopes and observational evidence that even weak plume tails reach the surface, suggest that subducting plates usually are not strongly resisted by this mechanism, and that putative episodes of mantle overturn would be much less dramatic than in some recent constant-viscosity models.