Depth-dependent Rheology Research Articles

Diapirs of crystal-rich slurry explain granite emplacement temperature and duration

The mechanism, temperature, and timescale of granite intrusion remain controversial, with wide-ranging implications for understanding continental growth, differentiation, rheology, and deformation dynamics. In this paper we present a method for determining intrusion emplacement temperature and timescale using the characteristics of the surrounding metamorphic aureole, and apply it to the Skiddaw granite in northern England. The estimated emplacement timescale (0.1–2 Myr) implies magma transport velocities of 1–100 mm/year. At the absent or low melt fractions relevant to our estimated emplacement temperature (580–650 ^{circ }C), such velocities are incompatible with pluton formation by successive injections through dykes. Instead, our results indicate the intrusion of a diapir of crystal-rich slurry, solidifying before emplacement, with a rheology governed by the solid crystals. The emplacement depth is likely to be governed by the depth-dependent rheology of the surrounding rocks, occurring close to the brittle-ductile transition. The wider implications of our results relate to (1) the appreciation that much of the chemical and textural characteristics of plutons may relate to pre-emplacement crystallisation at depth, passively transported to higher crustal levels, and (2) an explanation of the difficulty of seismically imaging active plutonism.

Planetary thermal evolution models with tectonic transitions

Thermal history calculations provide important insights into the interior evolution of planets, but incorporate simplified dynamics from the systems they represent. Planetary interiors typical incorporate complex rheologies, viscous layering, lateral heterogeneities, and time delays in processes, which have not been traditionally represented by parameterised approaches. Here we develop numerical models for planetary evolution, incorporating the physical complexity of Earth’s interior, and use them to generate statistically-based Nu-Ra scalings. These encapsulate the main effects of tectonic transitions, geometry, and depth-dependent rheology, and time-sensitivity. We find an exponent β of ~0.26 best describes the Nu-Ra relationship for evolving mobile lid systems, and β ~0.12 for stagnant-lid systems. Systems with time dependent subduction have β varying between ~0.26 during the Hadean, when external factors such as impacts facilitate tectonics, to ~0.12 during the Archaean, when the system is dominated by long periods of quiescence, and systems driven by external forcings (eg. due to impacts in the first 100Myr of Earth’s history) may exhibit much higher exponents. We also find a time-lag between Ra (which primarily depends on mantle temperature) and Nu (normalised surface heat flow) of around 200–300Myr, suggesting a significant delay between mantle thermal configuration, and its surface manifestation. These results provide an approach for the rapid characterisation of tectonic, volcanic, and atmospheric evolution of terrestrial exoplanets.

Enceladus’s crust as a non-uniform thin shell: I tidal deformations

The geologic activity at Enceladus’s south pole remains unexplained, though tidal deformations are probably the ultimate cause. Recent gravity and libration data indicate that Enceladus’s icy crust floats on a global ocean, is rather thin, and has a strongly non-uniform thickness. Tidal effects are enhanced by crustal thinning at the south pole, so that realistic models of tidal tectonics and dissipation should take into account the lateral variations of shell structure. I construct here the theory of non-uniform viscoelastic thin shells, allowing for depth-dependent rheology and large lateral variations of shell thickness and rheology. Coupling to tides yields two 2D linear partial differential equations of the fourth order on the sphere which take into account self-gravity, density stratification below the shell, and core viscoelasticity. If the shell is laterally uniform, the solution agrees with analytical formulas for tidal Love numbers; errors on displacements and stresses are less than 5% and 15%, respectively, if the thickness is less than 10% of the radius. If the shell is non-uniform, the tidal thin shell equations are solved as a system of coupled linear equations in a spherical harmonic basis. Compared to finite element models, thin shell predictions are similar for the deformations due to Enceladus’s pressurized ocean, but differ for the tides of Ganymede. If Enceladus’s shell is conductive with isostatic thickness variations, surface stresses are approximately inversely proportional to the local shell thickness. The radial tide is only moderately enhanced at the south pole. The combination of crustal thinning and convection below the poles can amplify south polar stresses by a factor of 10, but it cannot explain the apparent time lag between the maximum plume brightness and the opening of tiger stripes. In a second paper, I will study the impact of a non-uniform crust on tidal dissipation.

Tides on Europa: The membrane paradigm

Jupiter’s moon Europa has a thin icy crust which is decoupled from the mantle by a subsurface ocean. The crust thus responds to tidal forcing as a deformed membrane, cold at the top and near melting point at the bottom. In this paper I develop the membrane theory of viscoelastic shells with depth-dependent rheology with the dual goal of predicting tidal tectonics and computing tidal dissipation. Two parameters characterize the tidal response of the membrane: the effective Poisson’s ratio ν¯ and the membrane spring constant Λ, the latter being proportional to the crust thickness and effective shear modulus. I solve membrane theory in terms of tidal Love numbers, for which I derive analytical formulas depending on Λ,ν¯, the ocean-to-bulk density ratio and the number k2∘ representing the influence of the deep interior. Membrane formulas predict h2 and k2 with an accuracy of a few tenths of percent if the crust thickness is less than one hundred kilometers, whereas the error on l2 is a few percents. Benchmarking with the thick-shell software SatStress leads to the discovery of an error in the original, uncorrected version of the code that changes stress components by up to 40%. Regarding tectonics, I show that different stress-free states account for the conflicting predictions of thin and thick shell models about the magnitude of tensile stresses due to nonsynchronous rotation. Regarding dissipation, I prove that tidal heating in the crust is proportional to Im(Λ) and that it is equal to the global heat flow (proportional to Im(k2)) minus the core-mantle heat flow (proportional to Im(k2∘)). As an illustration, I compute the equilibrium thickness of a convecting crust. More generally, membrane formulas are useful in any application involving tidal Love numbers such as crust thickness estimates, despinning tectonics or true polar wander.

Plume-induced geoid anomalies from 2D axi-symmetric temperature- and pressure-dependent mantle convection models

A collection of numerical simulations of 2D axi-symmetric thermal convection is presented here. The aim is to investigate the shape of geoid anomalies and dynamic topography above a plume. The simulation is based on the Boussinesq approximation and infinite Prandtl number and is carried out in the spherical shell with strongly temperature- and depth-dependent Arrhenius-type viscosity. According to the Arrhenius law, plume models with purely depth-dependent rheology are unphysical and should be taken with care. The strongly coupled temperature- and depth-dependent viscosity enables us to better understand the plume's behavior inside the Earth.The topography and geoid anomalies produced from plumes are sensitive to rheology of the mantle and rheology of the plume; both have effects on shape and amplitude of the geoid anomalies. We determined different categories of the geoid which are related to various rheology. Depth-dependent viscosity models show a geoid with a negative sign above the plume, and temperature-dependent viscosity models depict a bell-shaped geoid. We identified different behaviors in the combined model with temperature- and pressure-dependent viscosity.

Flow trajectories in analogue viscous orogenic wedges: Insights on natural orogens

Flow trajectories in viscous orogenic wedges with depth-dependent rheology are modelled with an experimental thermomechanical apparatus and paraffin as the analogue material. Two different foreland-dipping backstop configurations (backstop dips of 60° and 30°, respectively) are used under constant convergence velocity (bulk deformation rate 5 × 10 − 6 s − 1 ). Results show that shortening is accommodated at depth by diverging shear zones (pro- and retrowedge-directed) , compensated at shallower levels by prowedge-directed shearing occurring in the pro-side of the axial region, whose efficiency increases with backstop dip. The resulting shear sense distribution allows the deep parts of the wedge axial region to be progressively exhumed along an extrusion channel in the retrowedge region. Exhumation is favoured by increasing the dip angle of the backstop. The ductile flow path at depth constitutes the main factor controlling shear sense distributions and metamorphic recurrences in the axial region of mature orogenic wedges.

Towards self-consistent modeling of the martian dichotomy: The influence of one-ridge convection on crustal thickness distribution

In order to find an explanation for the origin of the martian crustal dichotomy, a number of recent papers have examined the effect of layered viscosity on the evolution of a degree-1 mantle convection, e.g. Roberts and Zhong [Roberts, J.H., Zhong, S., 2006. J. Geophys. Res. 111. E06013] and Yoshida and Kageyama [Yoshida, M., Kageyama, A., 2006. J. Geophys. Res. 111, doi:10.1029/2005JB003905. B03412]. It was found that a mid-mantle viscosity jump, combined with highly temperature- and depth-dependent rheology, are effective in developing a degree-1 convection within a short timescale. Such a layered viscosity profile could be justified by martian mineralogy. However, the effect of a degree-1 convective planform on the crustal thickness distribution has not yet been demonstrated. It is not obvious whether a thinner crust, due to sublithospheric erosion and crustal thinning, or a thicker crust, due to enhanced crustal production, would form above the hemisphere of mantle upwelling. Also, the general shape of the dichotomy, which is not strictly hemispherical, has not yet been fully investigated. Here we propose a model of the crustal patterns produced by numerical simulations of martian mantle convection, using the finite-volume multigrid code StagYY [Tackley, P.J., 2008. Phys. Earth Planet. Int. 107, 7–18, doi:10.1016/j.pepi.2008.08.005] A self-consistent treatment of melting, crustal formation and chemical differentiation has been added to models of three-dimensional thermal convection. This allows us to obtain global maps of the crustal thickness distribution as it evolves with time. The obtained results demonstrate that it is indeed possible to form a crustal dichotomy as a consequence of near degree-1 mantle convection early in Mars' history. We find that some of the observed patterns show intriguing first order similarities to the elliptical shape of the martian dichotomy. In all models, the region of thick crust is located over the region of mantle upwelling, which itself is a ridge-like structure spread over roughly one half of the planet, a planform we describe as “one-ridge convection.”

Ultrafast mantle plumes and implications for flood basalt volcanism in the Northern Atlantic Region

Recent geochemical and geochronological data on Paleocene flood basalts from West Greenland, SE Greenland and the British Isles show that volcanic activity at these widely separated locations commenced nearly simultaneously around 61–62 Ma ago, and that the duration of this initial phase of flood basalt magmatism was of the order of a few million years or less. A small, fast moving upper mantle plume which rapidly spreads out horizontally, on encountering the base of the lithosphere, appears to be a viable mechanism to explain these observations. However, in order to reconcile the idea of an ultrafast plume with the observed plate velocities, a physical mechanism is needed for inducing a separation of timescales between the plume speed and the surrounding mantle circulation. Thermal convection with a non-Newtonian temperature- and depth-dependent rheology provides such a mechanism wherein extremely fast plumes ascending at velocities between one to tens of meters per year can be produced in an otherwise slowly convecting mantle moving at cm/yr. This transport mechanism is capable of bringing up very hot matter from the transition zone to the lithosphere. The fast upward velocities lead to the production of high viscous heating rates surrounding the plume upon impinging the lithosphere. It is thus possible to explain the near simultaneous onset of magmatism in West Greenland, SE Greenland and the British Isles in the early Tertiary by a fast moving mantle plume spreading out horizontally with a velocity of around 0.5 m/yr. For resolving numerically the thermal-mechanical state of these strongly time-dependent mantle flows, extremely high spatial resolution, on the kilometer scale, is required. Finally we suggest the possibility that the source of the fast upper-mantle plume under Iceland may be rooted in the lower mantle. This is consistent with the recent findings by seismic tomography of a deep mantle plume under Iceland.

The role of subduction in the evolution of the Apennines foreland basin

Seismic tomographic images of the Mediterranean region lead to the formulation of an hypothesis involving detachment (or alternatively break-off) of subducted lithosphere, propagating horizontally in the along strike direction: lateral migration of slab detachment [Wortel, M.J.R., Spakman, W., (1992). Structure and dynamics of subducted lithosphere in the Mediterranean region. Proc. K. Ned. Akad. Wet. 95, 325-347]. On the basis of this hypothesis variations in vertical loads are expected along strike of the convergent plate margin involved. We investigate whether such variations in vertical loads can be inferred from the Plio-Pleistocene foreland basin along the eastern side of the Apennines. Our aim is to test the hypothesis of a laterally migrating slab detachment beneath Italy. We employ a finite difference method to obtain the flexural downbending of the Adriatic plate resulting from known loads, such as topography and sediments. The unknown slab pull is inferred by constraining the model to fit the horizon that is assumed to reflect flexural downbending in two-dimensional sections perpendicular to the strike of the Apennines. Our method forms a consistent treatment of all loads involved in two-dimensional lithospheric flexure. If we assume that the downgoing Adriatic plate behaves purely elastically, we conclude that there is a trade-off relation between slab pull and elastic plate thickness for sections in north Italy. For these sections gravity anomaly data cannot constrain a best elastic model. If we make the more realistic assumption of a continental lithosphere with a depth-dependent rheology, we obtain a good match for sections in the north. In the south, however, we find that a match cannot be obtained. This may imply that the data reflect more than purely flexural processes — probably thrusting. The conclusion is that a trend in vertical loads along strike of the Apennines cannot be resolved from Plio-Pleistocene deflection data. Therefore, the hypothesis of lateral migration of slab detachment can neither be shown to be consistent nor inconsistent with the observations.

Lateral variations in lithosphere strength in the Romanian Carpathians: constraints on basin evolution

Lithospheric strength profiles constructed through extrapolation of rock mechanics data and constrained by crustal geophysics data, demonstrate the existence of important spatial variations of thermomechanical properties in different domains of the Romanian Carpathians, the Pannonian basin and the Transylvanian depression. These models show important lateral variations in lithospheric rigidity along the Southern and Eastern Carpathian foreland, compatible with inferences from flexural modelling studies. The modelling predicts an absence of mantle strength and the presence of weak lithosphere in the Pannonian and Transylvanian basins. A NWSE-trending rheological transect is presented, connecting the Moesian platform, the Vrancea region of the Carpathians, the Transylvanian basin and the Apuseni Mountains with the Romanian part of the Pannonian basin. This transect demonstrates the existence of a close correlation of the inferred spatial variations in mechanical properties of the Romanian lithosphere with observed patterns in large-scale gravity anomalies and with tectonic units. A notable feature of the rheological models is the presence of an intra-Carpathian weak zone. The combination of a strong foreland lithosphere with a weak lithosphere in the inner part of the Carpathian/Transylvanian/Pannonian system sheds light on a number of key features observed in the evolution of sedimentary basins in these areas. Remarkable differences in rheology predictions between flexural models and depth-dependent rheology models indicate that the crustal lithosphere is detached from the mantle part of the lithosphere in the Carpathian foreland.

Subsidence analysis and quantitative basin modelling in the Styrian Basin (Pannonian Basin System, Austria)

We present the results of subsidence analysis and of quantitative basin modelling using isostatical and flexural models for basin evolution along four cross-sections in the Styrian Basin, the westernmost subbasin of the Pannonian Basin System. Subsidence analysis reveals a first Ottnangian-Karpatian synrift phase. Our local isostatic models predict crustal stretching values up to 1.3 and subcrustal stretching values of 1.6 for this event. Stretching factors of a minor Sarmatian extension phase are below 1.04. The termination of subsidence during the Pannonian and a rapid Quaternary uplift phase can be explained by major changes in the regional stress field. A W-E cross-section through the northern Fürstenfeld Subbasin provides a key for the understanding of the dynamics of basin formation. It crosses a narrow Karpatian rift basin, the metamorphic core complex of the Penninic Eisenberg Window and shows an eastward tilting of the easternmost part of the basin during Pannonian times. Uplift of the Penninic window can only be modelled with an extremely weak lithosphere (equivalent elastic thickness (EET) ⪡ 2 km), whereas a best fit between observed and modelled tilting is obtained with an EET value of 5 km. These results suggest that the lithosphere was extremely weak during the onset of basin evolution in Ottnangian-Karpatian times, probably caused by high extension rates and high heat flows associated with Karpatian to early Badenian magmatic activity. Subsequent cooling led to a pronounced increase in flexural rigidity. An EET of 5 km fits well with estimations in other parts of the Pannonian realm. Depth-dependent rheology models based on palaeo-heat flow estimates indicate a similar increase in lithospheric strength with time. The impact of Plio-Pleistocene volcanism on rheology appears to be relatively modest, which can be explained by a deep position of the magma chamber for this event.

Effect of depth-dependent viscosity on the planform of mantle convection

LITHOSPHERIC plate motions at the Earth's surface result from thermal convection in the mantle1. Understanding mantle convection is made difficult by variations in the material properties of rocks as pressure and temperature increase from the surface to the core. The plates themselves result from high rock strength and brittle failure at low temperature near the surface. In the deeper mantle, elevated pressure may increase the effective viscosity by orders of magnitude2–5. The influence of depth-dependent viscosity on convection has been explored in two-dimensional numerical experiments6–8, but planforms must be studied in three dimensions. Although three-dimensional plan-forms can be elucidated by laboratory fluid dynamic experiments9,10, such experiments cannot simulate depth-dependent rheology. Here we use a three-dimensional spherical convection model11,12 to show that a modest increase in mantle viscosity with depth has a marked effect on the planform of convection, resulting in long, linear downwellings from the upper surface boundary layer and a surprisingly 'red' thermal heterogeneity spectrum, as observed for the Earth's mantle13. These effects of depth-dependent viscosity may be comparable to the effects of the plates themselves.

Response of the lithosphere to in-plane force variations

The Earth's lithosphere deforms when subjected to vertical and horizontal tectonic forces because it is not perfectly rigid. The deformation has both elastic (flexural) and inelastic (brittle and ductile creep) components. We use the yield strength envelope formulation of a depth-dependent rheology to predict how the lithosphere will deform in response to applied in-plane force. Brittle deformation of the upper lithosphere is influenced by the distribution of existing faults and fractures in the crust and their orientation with respect to the applied tectonic force. The sense of fault reactivation indicates the sign of the applied in-plane force. The induced flexural deformation is governed by the pre-existing lithospheric deflection, the effective elastic thickness of the lithosphere when the force was applied, and the magnitude of the applied force. Using this framework, we compare the predicted response of the lithosphere to in-plane force variations with that observed in the Central Indian Ocean and the Jeanne d'Arc basin, offshore Newfoundland. The brittle deformation in both regions is readily observed in seismic reflection data. The upper oceanic crust in the Central Indian Ocean is broken into blocks bounded by high-angle reverse faults. The dip of these basement faults and their orientation parallel to the marine magnetic lineations suggest that they represent a reactivated crustal fabric imparted during the formation of oceanic crust. The reverse sense of fault reactivation is consistent with the north-south compression determined from earthquake focal mechanisms. In the Jeanne d'Arc basin, episodes of northwest-southeast extension reactivated the basin-bounding faults in a normal sense during late Barremian, early Aptian, and late Aptian time. In the late Albian, the same crustal blocks were reactivated and uplifted in response to north-northeast in-plane compression. Flexural deformation due to applied in-plane force is observed only in the Central Indian Ocean case, where the oceanic crust and most of the overlying Bengal Fan sediments were deformed into broad, east-west trending undulations with wavelengths of ≈ 100–300 km and amplitudes up to ≈ 2 km. In the Jeanne d'Arc basin, the flexural deformation which should accompany the fault reactivation is difficult to recognize from the preserved stratigraphy. Consequently, brittle deformation within extensional basins, which appears to be the most recognizable response of the lithosphere to in-plane force variations, obscures the effects of the flexural deformation. However, it is this flexural deformation that has been postulated as a tectonic alternative to glacial-eustasy in explaining the 3 rd-order sea level cycles. Because it is difficult to recognize the effects of the flexural deformation in the preserved stratigraphy, we surmise that it is not a viable alternative to glacial-eustasy in generating either regional or global unconformities. Instead, we propose that the flexural and brittle deformation due to applied in-plane force may locally enhance or diminish the effects of other processes responsible for unconformity generation.

3-D Convection With Variable Viscosity

Accepted 1990 September 6. Received 1990 July 16; in original form 1990 May 3 SUMMARY report numerical calculations for 3-D convection with variable viscosity. A hybrid spectral and finite difference method is used. The coupling of modes in the equation of motion, which is caused by lateral viscosity variations, is treated iteratively. Solutions for bimodal, hexagonal, square, triangular and spoke patterns are reported for bottom heated convection at infinite Prandtl number. The Rayleigh number, based on the viscosity at the mean of top and bottom temperature, is between critical and lo5, and temperature-induced viscosity contrasts up to 100 are considered (lo00 in one case). In agreement with results from laboratory experiments we find that at low Rayleigh number temperature-dependent viscosity favours flow patterns like squares or hexagons, where a columnar rising current is surrounded by sheet-like descending flow. The dichotomy in geometry between upwelling and sinking flow becomes more pronounced with increasing viscosity contrast. The temperature dependence of viscosity gives rise to a toroidal velocity component; however, it amounts only to a few per cent of the total velocity. In contrast, at the earth’s surface an approximate equipartitioning of poloidal and toroidal energy is found. We show that with non-Newtonian and depth-dependent rheology the toroidal component at the free surface can become significant, and a pattern reminiscent of plate motion can arise in a free convection model. Although these results are obtained in a parameter range which is not directly applicable to the earth, they support the conclusions that (i) upwelling flow in the mantle is unlikely to be sheet-like and will probably be in the form of columnar plumes, and that (ii) the toroidal motion found at the earth’s surface is due to the highly non-linear rheology which leads to the existence of mobile surface plates and is not caused by viscosity variations related to lateral temperature contrasts deeper in the mantle.

Post-seismic stress evolution in a boundary-zone model with depth-dependent rheology

Tectonic boundary zones are characterized by a shallow, brittle region (seismogenic region) where frictional slip occurs on fault planes, and by a lower, ductile and aseismic region. The nucleation depth and the extent of slip planes in great interplate earthquakes are determined by the relationships between frictional resistance, applied stress and stress drop. The post-seismic stress evolution on a fault plane is investigated by a model which reproduces a transcurrent boundary zone subject to a constant strain rate. Dislocation models with frictional stress threshold for slip arrest must be considered in order to study post-seismic stress evolution in the proximity of dislocation edges. It is found that dislocations which are nucleated in the seismogenic region can propagate downward into the aseismic region: in this case earthquakes are followed by relaxation of the stress concentration at the lower dislocation edge. This has effects on the space-time distribution of aftershocks.

Rheological modelling and deformation instability of lithosphere under extension--II. Depth-dependent rheology

SUMMARY Bassi & Bonnin (1988) have recently discussed the stability of continental lithosphere in extension. They examine the hypothesis that thickness inhomogeneities can be amplified during extension and be responsible for lithospheric necking. We carry our analysis one step further by accounting for the increase of strength with depth in brittle layers, namely the upper crust and, in some cases, the upper part of the mantle. This variation had been neglected in the previous study to simplify the mathematical formulation. Rather than a constant, the equivalent viscosity is now a linear function of depth, as predicted by Byerlee’s friction law which is postulated to control the mechanical behavior of brittle layers. This more realistic rheological model confirms the results of the original analysis, specifically that the interfaces between layers tend to flatten during extension unless they correspond to an inversion of density. Therefore, extension of continental lithosphere appears to be stable if there is no other initial defect susceptible of localizing the deformation.

Model geoid anomalies due to subduction of inextensible lithosphere

We compute geoid slopes from models of subduction in which the subducted lithosphere is much stronger than the surrounding mantle. Geoid slope contributions from both the lithospheric slab and mantle boundary deformations are computed from finite element analysis of mantle flow. The finite element model includes a slab of finite length and a depth dependent Newtonian rheology for the surrounding mantle. We find that observed geoid anomalies at subduction zones, which are positive, cannot be matched by models with uniform mantle viscosity. However, even with a strong subducted lithosphere, the ratio of driving load to boundary deformation is significantly increased by a ten‐fold increase of viscosity with depth, resulting in a geoid high. We find that the sign of the geoid slopes within 3000 km of the trench are independent of maximum depth of the slab for maximum depths from 700 km to 2800 km.