Viscoelastic Lithosphere Research Articles

Feedback mechanisms controlling Antarctic glacial-cycle dynamics simulated with a coupled ice sheet–solid Earth model

Abstract. The dynamics of the ice sheets on glacial timescales are highly controlled by interactions with the solid Earth, i.e., the glacial isostatic adjustment (GIA). Particularly at marine ice sheets, competing feedback mechanisms govern the migration of the ice sheet's grounding line (GL) and hence the ice sheet stability. For this study, we developed a coupling scheme and performed a suite of coupled ice sheet–solid Earth simulations over the last two glacial cycles. To represent ice sheet dynamics we apply the Parallel Ice Sheet Model (PISM), and to represent the solid Earth response we apply the 3D VIscoelastic Lithosphere and MAntle model (VILMA), which, in addition to load deformation and rotation changes, considers the gravitationally consistent redistribution of water (the sea-level equation). We decided on an offline coupling between the two model components. By convergence of trajectories of the Antarctic Ice Sheet deglaciation we determine optimal coupling time step and spatial resolution of the GIA model and compare patterns of inferred relative sea-level change since the Last Glacial Maximum with the results from previous studies. With our coupling setup we evaluate the relevance of feedback mechanisms for the glaciation and deglaciation phases in Antarctica considering different 3D Earth structures resulting in a range of load-response timescales. For rather long timescales, in a glacial climate associated with the far-field sea-level low stand, we find GL advance up to the edge of the continental shelf mainly in West Antarctica, dominated by a self-amplifying GIA feedback, which we call the “forebulge feedback”. For the much shorter timescale of deglaciation, dominated by the marine ice sheet instability, our simulations suggest that the stabilizing sea-level feedback can significantly slow down GL retreat in the Ross sector, which is dominated by a very weak Earth structure (i.e., low mantle viscosity and thin lithosphere). This delaying effect prevents a Holocene GL retreat beyond its present-day position, which is discussed in the scientific community and supported by observational evidence at the Siple Coast and by previous model simulations. The applied coupled framework, PISM–VILMA, allows for defining restart states to run multiple sensitivity simulations from. It can be easily implemented in Earth system models (ESMs) and provides the tools to gain a better understanding of ice sheet stability on glacial timescales as well as in a warmer future climate.

Hydrogeologic and Thermal Effects of Glaciations on the Intracontinental Basins in Central and Northern Europe

We use a fully coupled hydro-thermal model (TH) to quantify changes in the pore pressure and temperature distribution following the Last Glacial Maximum (LGM) in the intracontinental basins in Central and Northern Europe. We demonstrate that even without considering a direct mechanical coupling from the visco-elastic lithosphere rebound, the system is, at present-day, in a state of hydrogeologic and thermal disequilibrium as a result of the past ice sheet dynamics. We find that the local geology exerts an additional control on the subsurface response to imposed glacial loading, as evidenced by a contrasting thermal and pore pressure configuration in time and space. Highest rates of pore pressure dissipation are restricted to crustal domains that underwent substantial glacial loading, while the majority of the sedimentary sub-basins show a prominent signature of hydraulic disequilibrium (overpressure) at present. Groundwater-driven convective cooling and heating during the advance and retreat of the ice cap occurred mainly within sedimentary rocks, domains where thermal equilibration is ongoing. The spatial correlation between modeled pore pressure dissipation rates and postglacial uplift rates is indicative of a complex and transient hydrogeological system structurally connected to the viscous tail of the ongoing isostatic adjustment after the LGM, with important implications for assessing the long-term mechanical stability of this intraplate setting.

Melt Percolation, Concentration and Dyking in the Hawaiian Mantle Plume and Overriding Lithosphere: Links to the Evolution of Lava Composition along the Volcanic Chain

Abstract Oceanic island basalts and related magmatic rocks from Hawaii are derived from a compositionally heterogeneous mantle plume. Here we describe how this heterogeneity results from the transport of filaments of a specific composition in the plume, representing a relatively small volume of rocks (~15 %) interbedded inside a dry peridotite mantle. Four types of filaments are considered: sub-primitive mantle, ultralow-velocity zone, fertilized-harzburgite and eclogite type filaments. We present a model that describes the flow within a plume and the stress field in the overriding viscoelastic lithosphere and that can determine, from depth to the surface, the melting rate, composition and trajectory of melts produced within each type of filament. Our model shows that (1) the filaments melt at a depth corresponding to >5 GPa, where the temperature gap between the solidus and liquidus is narrow (~40–80 °C), and (2) the volume of filaments is small relative to the total volume of mantle, which therefore allows the latent heat required for the partial melting to be provided via conduction inside the hot plume. The primitive melts produced inside the filaments, occasionally mixed with the melt derived from an eclogite filament, represent a volume comparable with that expected in a plume composed only of dry peridotite that partially melts to a degree of ~10 % at the interface between the spinel and garnet fields (60–70 km depth). In particular, in the centre of the plume, sub-primitive mantle filaments produce up to 30 % tholeiite–picrite melts, whereas in fertilized-harzburgite filaments, the mantle melts completely to produce a melt having a meimechite-like composition. A key finding is that the fractional crystallization of these melts probably forms the so-called ‘primary mantle-derived alkaline magmas’ along with dunites and olivine-rich cumulates. Our plume model shows that the mantle flow divides into two parts. The first corresponds to hot flowlines that originate at a depth of ~200 km and at a distance of less than 25 km from the plume axis. Along these flowlines, when the mantle reaches a pressure of 5 GPa, the partially molten horizon in filaments is sufficiently thick for the interstitial melt to be squeezed out via dykes. This melt eventually ponds as sills in a subrectangular zone that is located inside the overlying lithosphere, between 70 and 50 km depth and centred over a distance of less than 40 km on either side of the axis. This zone is designated as the shield magmatic reservoir. The volatile-rich melt inside the sills infiltrates the surrounding mantle lithosphere and partially melts it. After ~0·1 Myr, the melt resumes its vertical ascent via dykes and eventually ponds and differentiates within subcrustal magma chambers located below active shield volcanoes. This sequence of processes matches the expected volume, petrology and geochemistry recorded for shield volcanoes. The second part of the melt flow does not pond within the shield magmatic reservoir. Rather, the mantle cold flowlines, originating at ~200 km depth and at 25–35 km from the plume axis, discharge their interstitial melt through dykes that were initially generated deeper, at ~5 GPa. The melt reaches the Moho at 100–150 km from the plume axis, where it forms magmatic bodies within which the melt differentiates. This melt probably represents that observed in pre- and postshield volcanoes. Finally, at ~70 km from the plume axis and at a depth greater than 200 km, the flowlines are subvertical. They then deflect at ~180 km depth and rotate toward the horizontal and eventually transit at 10–20° to the horizontal across an ~200 km distance from the axis and reach ~140 km depth. The fertilized-harzburgite and sub-primitive mantle/ultralow-velocity zone filaments that flow along these elbows partially melt by a few to several per cent. The resulting interstitial melt has a kimberlite-like composition. Thereafter, the excess pressure at the top of the filament at ~200 km from the axis overcomes the threshold for dyking and thus allows the escape of the interstitial melt via dykes ponding in subcrustal magma chambers or emerging directly at the surface. These melts have a composition similar to that associated with rejuvenated volcanism. We use the nature and the composition of whole erupted magmas and the seismic structure along the Hawaiian chain to validate this model.

Time-dependent rotational stability of dynamic planets with viscoelastic lithospheres

We extend previous work to derive a non-linear rotational stability theory governing true polar wander (TPW) of terrestrial planets with viscoelastic lithospheres. We demonstrate, analytically and using numerical examples, that our expressions are consistent with previous results in the limiting cases of low and infinite (i.e., purely elastic) viscosity lithospheres. To illustrate the stability theory, we compute TPW on Mars driven by a simple, prescribed mass loading. Our calculations demonstrate that on short time scales relative to the relaxation time of the viscoelastic lithosphere, the rotation axis follows a constrained path that reflects stabilization by remnant strength in the lithosphere, but that on long times scales this stabilization disappears and the load ultimately reaches the equator. Earlier work based on the assumption of a permanent remnant bulge in the case of a purely elastic lithosphere has suggested that Martian TPW would not persist for any significant time after a load is emplaced, and thus an equilibrium stability theory is sufficient to model long-term (order 1 Myr or longer) polar motion of the planet. Our results suggest, in contrast, that TPW on Mars can continue over time scales on the order of the relaxation time of the lithosphere after load emplacement; for sufficiently high lithospheric viscosities, this time scale may be comparable to the age of the planet.

Dynamics of plate bending at the trench and slab‐plate coupling

The bending strength of subducting lithosphere plays a critical role in the Earth's plate tectonics and mantle convection, modulating the amount of slab pull transmitted to the surface and setting the boundary conditions under which plates move and deform. However, it is the subject of a lively debate how much of the potential energy of the downgoing plate is consumed in bending the plate and how the lithospheric strength is defined during this process. We model the subduction of a viscoelastic lithosphere, driven solely by the downgoing plate's buoyancy, freely sinking in a passive mantle, represented by drag forces. To investigate the dynamics of bending, (1) we vary the density and the viscosity profile within the plate from isoviscous, where strength is distributed, to strongly layered, where strength is concentrated in a thin core, and (2) we map the stress, strain, and dissipation along the downgoing plate. The effective plate strength during bending is not a simple function of average plate viscosity but is affected by rheological layering and plate thinning. Earth‐like layered plates allow for the transmission of large fractions of slab pull (∼75–80%) through the bend and yield a net slab pull of FSPnet = 1 to 6 × 1012 N m−1, which varies with the buoyancy of plates. In all models, only a minor portion of the energy is dissipated in the bending. Surprisingly, bending dissipation hardly varies with lithospheric viscosity because in our dynamic system, the plates minimize overall dissipation rate by adjusting their bending curvature. As a result, bending dissipation, ΦB, is mainly controlled by the bending moment work rate exerted by slab pull. We propose a new analytical formulation that includes a viscosity‐dependent bending radius, which allows for assessment of the relative bending dissipation in the Earth's subduction zones using parameters from a recent global compilation. This yields estimates of ΦB/ΦTOT < 20%. These results suggest that plates on Earth weakly resist bending, yet are able to propagate a large amount of slab pull.

The effect of a vertical temperature gradient in the lithosphere on seismic and tectonic waves

A linear analysis of the stability of the lithosphere considered as a viscoelastic layer with an equilibrium vertical gradient of temperature is carried out. The problem is solved with a complete system of linearized equations of a continuous medium represented in the dimensionless form and containing a set of dimensionless parameters that determine thermomechanical properties of the lithosphere. As a result of the stability analysis, decrements are found that give the time dependence of perturbations and correspond to high-frequency seismic waves and low-frequency tectonic waves. The frequency and velocity of seismic waves are determined by the elasticity and inertial properties of the lithosphere, and their attenuation, by viscous properties of the lithosphere. The temperature gradient existing in the lithosphere influences seismic waves very weakly. On the contrary, the pattern of tectonic waves is controlled by the temperature gradient and viscous properties, while the effect of elastic and inertial properties on these waves is negligibly small. The stability of a viscoelastic lithosphere is examined using such rheological models as the Maxwell, standard linear, and Andrade media (the frequency of tectonic waves is zero in the Maxwell medium).

Polar wander caused by the Quaternary glacial cycles and fluid Love number

Perturbations of the Earth’s rotation caused by the Quaternary glacial cycles provide an important constraint on the viscosity of the deep mantle because they represent a long-wavelength response of the Earth to surface load redistribution. The predicted present-day polar wander speed (PWS) is, however, sensitive to both the lower mantle viscosity ( η lm), the density jump at 670 km depth, and the lithospheric thickness and viscosity (e.g., Sabadini and Peltier, Geophys. J. R. Astron. Soc. 66 (1981) 553–578; Yuen et al., J. Geophys. Res. 87 (1982) 10745–10762; Peltier and Wu, Geophys. Res. Lett. 10 (1983) 181–184; Wu and Peltier, Geophys, J. R. Astron. Soc. 76 (1984) 753–791; Peltier, J. Geophys. Res. 89 (1984) 11303–11316; Vermeersen et al., J. Geophys. Res. 102 (1997) 27689–27702; Mitrovica and Milne, J. Geophys. Res. 103 (1998) 985–1005; Johnston and Lambeck, Geophys. J. Int. 136 (1999) 537–558; Nakada, Geophys. J. Int. 143 (2000) 230–238). For earth models with η lm<5×10 21 Pa s and an elastic lithosphere, the present-day PWS is very sensitive to the M1 mode (buoyancy mode) related to the density jump at 670 km depth [Mitrovica and Milne, J. Geophys. Res. 103 (1998) 985–1005]. The contribution of the M1 mode, however, is less significant for earth models with a viscoelastic lithosphere [Nakada, Geophys. J. Int. 143 (2000) 230–238]. This is due to the fact that this contribution depends on the relative strength of the M1 mode, Δ k 2 T(M1)/ k f T, where Δ k 2 T(M1) is the magnitude of tidal Love number ( k 2 T) of the M1 mode and k f T is the value of k 2 T in the fluid limit (fluid Love number). The magnitude of k f T for earth models with a viscoelastic lithosphere is larger than that for an elastic lithosphere, and it is smaller for a thicker elastic lithosphere than for a thinner one. Thus, for earth models with a viscoelastic lithosphere, the PWS is mainly sensitive to the lower mantle viscosity regardless of the behavior of the 670 km density discontinuity. This relation also explains why the predicted PWS increases with increasing thickness of an elastic lithosphere. That is, since the value of Δ k 2 T(M1)/ k f T with a thicker elastic lithosphere is larger than that with a thinner elastic lithosphere, the M1 mode will have a higher contribution in the case of a thicker elastic lithosphere.

Modelling of viscoelastic plume-lithosphere interaction using the adaptive multilevel wavelet collocation method

SUMMARY Modelling of mantle flows with sharp viscosity contrasts in a viscoelastic medium is a challenging computational problem in geodynamics because of its multiple-scale nature in space and time. We have employed a recently developed adaptive multilevel wavelet collocation algorithm to study the dynamics of a small rising diapir interacting with a stiff lithosphere in a Maxwell viscoelastic mantle. In this kinematic model we have prescribed the upward velocity of the diapir and then we need to integrate in time only the momentum equation governing the temporal evolution of the pressure, stress and velocity components, which together constitute a sixth-order system in time. The total number of collocation points did not exceed 10 4 , compared to more than 10 6 gridpoints using conventional evenly spaced grid methods. The viscosity of the diapir is 10 x4 times lower than that of the surrounding mantle, while the viscosity of the thin lithosphere, about 5‐10 per cent of the entire layer depth, is 10 4 ‐10 8 times stiffer than the ambient mantle. Our results demonstrate the efficacy of wavelets to capture the sharp gradients of the stress and pressure fields developed in the diapiric impingement process. The interaction of the viscoelastic lithosphere with the rising viscoelastic diapir results in the localization of stress within the lithosphere. The magnitude of the stress fields can reach around 100‐300 MPa. Our simple kinematic model shows clearly that viscoelasticity can potentially play an important role in the dynamics of the lithosphere, especially concerning the potential severage of the lithosphere by mantle upwellings.

Effects of a viscoelastic lithosphere on the isostatic bedrock response

The lithosphere responds to loading by elastic flexure which is followed by viscous relaxation, the amount of which depends on the stress duration. This study compares results of an Earth model in which the lithosphere is modelled as a purely elastic layer and as a more general viscoelastic solid overlain by a rigid crust. It shows the emergence of a noticeable difference in the short wavelength (1×10 2–5×10 2 km) component of the bedrock deformation after loading durations of the order of 10 5–10 6 yr assuming a lower lithosphere viscosity of the order of 10 23–10 24 Pa.s. In particular, for the long-term loading hypothesised to be imposed by the East Antarctic ice sheet, we find that a viscoelastic lithosphere yields a more local deformation pattern to which ice sheet dynamics are highly sensitive. It confirms that modelling of the Antarctic long-term evolution would benefit from a fully coupled ice/bedrock approach in which the lithosphere would be represented by a viscoelastic solid.

Effect of the viscoelastic lithosphere on polar wander speed caused by the Late Pleistocene glacial cycles

Previous studies of the wander of the rotation pole associated with the Late Pleistocene glacial cycles indicate that the predicted polar wander speed is sensitive to the density jump at the 670 km discontinuity, the thickness of the elastic lithosphere, and the lower mantle viscosity. In particular, the M1 mode related to the density jump at 670 km depth has been shown to contribute a dominant portion of predicted polar wander speed for sufficiently small lower mantle viscosities. In this study, we examine the sensitivity of polar wander to variations in the viscosity of the viscoelastic lithosphere using simplified compressible Maxwell viscoelastic earth models. Model calculations for earth models with a viscoelastic lithosphere of finite viscosity indicate that the contribution of the M1 mode is similar to those associated with the density discontinuity at the core–mantle boundary (C0 mode) and the lithosphere (L0 mode). We speculate that this is due to the interaction between the M1 mode and the transient mode associated with the viscoelastic lithosphere, which reduces the magnitude of polar wander rates. Therefore, the M1 mode does not contribute a dominant portion of the predicted polar wander speed for earth models with a viscoelastic lithosphere of finite viscosity. In this case, predictions of polar wander speed as a function of lower mantle viscosity exhibit the qualitative form of an ‘inverted parabola’, as predicted for the J˙2 curve. We caution, however, that these results are obtained for simplified earth models, and the results for seismological earth models such as PREM may be complicated by the interaction between the M1 mode and the large set of transient modes.

Stress fields and tectonic motions induced by time-varying density anomalies in the lithosphere

Abstract We explore the extent to which a viscoelastic rheology of the lithosphere could affect the deviatoric stress and vertical deformation driven by deep seated, time-dependent density anomalies within the lithosphere. These sources can be associated with fast density instabilities of thermal, chemical or tectonic origin. In particular, we focus on the stress transfer from the ductile lower lithosphere to the elastic crust. The model is radially stratified and employs the Maxwell rheology. Stress transfer takes place from a wide region within the lower lithosphere to the thinner elastic top layer where it is amplified. The time history of the source and the lithospheric viscosity play a crucial role in determining the time evolution of the stress and deformation patterns. Attention is drawn to the temporal rates of maximum shear accumulation, which are interpreted as earthquake recurrence times. For realistic values of the model and source parameters, the increase of the deviatoric stress within the crust can be of a few bars on time-scales of 10 3 years, for lithospheric viscosities of around 10 22 to 10 23 Pa s −1 . Lower levels of stress variations are found in the viscoelastic lithosphere. These findings are consistent with the low levels of stress release and long recurrence times that characterize the intraplate seismicity, commonly attributed to density instabilities. As a specific application, we considered the tectonic subsidence of the Tyrrhenian basin. From the rate of subsidence derived from the ODP Leg 107 drilling campaign, we can estimate an upper bound for the lithospheric viscosity of around 5 × 10 22 Pa s −1 , which is appropriate for active tectonic regions.

Stress accumulation in a visco-elastic lithosphere containing a buried creeping vertical strike-slip fault : Pal, B S Phys Earth Planet InterV54, N1–2, April 1989, P55–63

Stress accumulation in a visco-elastic lithosphere containing a buried creeping vertical strike-slip fault : Pal, B S Phys Earth Planet InterV54, N1–2, April 1989, P55–63

Stress accumulation in a visco-elastic lithosphere containing a buried creeping vertical strike-slip fault

A model of a long aseismically creeping vertical strike-slip fault buried in a visco-elastic lithosphere is considered. It is assumed that tectonic forces in the lithosphere maintain a steady shear stress far away from the fault, while the stresses near the fault may change due to fault creep. Exact solutions are found for the displacements and stresses in the system. The mathematical results are compared with surface deformations near strike-slip faults.

Ground motion and stress accumulation driven by density anomalies in a viscoelastic lithosphere. Some results for the Apennines

SUMMARY We provide the analytical formulation for calculating the displacement and stress field forced by internal sources in a stratified, self-gravitating, viscoelastic earth. This model is specialized to study the rate of vertical motion and shear stress accumulation produced by lithospheric density anomalies. These sources are allowed to vary in the lateral direction. We show that sphericity plays a crucial role for elongated lithospheric anomalies while self-gravitation produces minor deviations from a gravitating Earth. When the model is applied to the Apennines we get, for lithospheric viscosity ranging between 10 and loz3 Pas, the subsidence of the plate underlying the active front of the overthrusting load to be around 0.5-1.0 mm yr-'. This is consistent with the amount of sedimentation in the Adriatic foredeep. The deformation pattern is very peculiar, with the largest subsidence localized beneath the active front of the topography. Our model enlightens the impact of discontinuities of tectonic phases on vertical motions in collision zones. If lithospheric viscosity is around 10'' Pa s, vertical motions decay drastically on time scales of 105yr if lateral migration of density anomalies comes to an end. For higher viscosities, deformation rates are maintained longer. This correlation between horizontal and vertical motions suggests that altimetric geodetic surveying along levelling lines of a few hundred kilometers can be an important tool to constrain the tectonics of the studied region. Results are also shown for vertical motions along a transect perpendicular to the Apennines, when the crustal structure inverted from Bouguer gravity data is considered. Analysis of the stress field induced by an overthrusting load shows that principal stress differences of a few bar (or a few tenths of MPa) can be accumulated on time scales of 102-103yr. These low values agree with the average stress drop of earthquakes in the Appalachians and northern Apennines where our modelling can be applied. We find that lateral density variations certainly contribute to intraplate stresses, but they are less efficient in triggering earthquakes than other mechanisms, such as transcurrent motions along active plate margins. Seismicity induced by lateral variations of crustal and lithospheric density must be moderate, characterized by long return times. These results are in agreement with the recorded seismicity in the northern Apennines where the largest earthquakes have return times of lo2 yr. If shear stress is forced by an overthrusting load, we find that the largest rate of stress accumulation in the crust is concentrated beneath the active front, close to the boundary with the ductile lithosphere. Discontinuities of tectonic phases play an important role in controlling the amount of shear stress due to density anomalies.

A thermo-mechanical model of continental lithosphere

An understanding of the mechanical response of the continental lithosphere to vertical surface loading is crucial when considering geodynamic processes such as intraplate basin formation, postglacial rebound or mountain building. Previous models of the mechanical behaviour of the lithosphere, which are more successful in oceanic than in continental regions, have used specific isotherms to delineate an effectively elastic zone of the lithosphere and result in mathematically convenient descriptions of mechanical behaviour. Heat flow and lithosphere flexure data for continents, however, suggest that this relationship is more complex, which can be demonstrated by numerical modelling of the loading of a heterogeneous visco-elastic lithosphere with temperature-dependent viscosity. Viscous stress relaxation in the lower and hotter section of the lithosphere results in a region of stress concentration dependent on the thermal state of the lithosphere, its structure and the time since loading. Here we present initial numerical models, based on a relatively simple steady-state thermal model constrained by present-day heat flow, which explain observations of continental flexure except in regions which have experienced a long or complex thermal history since loading.

Stress and relaxation time in the viscoelastic lithosphere inferred from the outer topographic rise

The deformation of a thin viscoelastic plate is usually analyzed by use of an equation that is derived with an artificial assumption of the Poisson's ratio equal to 1/2. It is shown in this paper that a simple treatment applicable to an arbitrary value of Poisson's ratio is possible by introducing a new variable that represents the pressure gradient between the top and bottom of the plate. In this treatment a one‐dimensional stationary deformation of moving plate is given by a simple analytic equation. This equation is used to analyze the bathymetry of the outer topographic rise seaward of the deep‐sea trenches. The analysis reveals that the lithosphere in the western Pacific trenches is subject to an extensional horizontal force of 100 MPa or greater. The horizontal force may be much less or even compressive in the eastern Pacific trenches. To prevent horizontal force and bending stress from being extraordinarily large, the Maxwell relaxation time of the plate must be 1 m.y. or shorter, and the effective viscosity must be 1024 Pa s or less.

Flexural models of continental lithosphere based on the long-term erosional decay of topography

The most easily observed geophysical characteristic of continents is their surface topography. Topographic power in the spectral range 1000-100 km has been computed for North American geological provinces and is found to decrease systematically as the tectonic age of the underlying crust increases. These data have been quantitatively modelled in terms of topographic erosion and stress relaxation within the continental lithosphere through geological time spans. The lithosphere is treated as a thin (Maxwell) viscoelastic plate which deforms flexurally in response to the vertical surface loads of topography and its erosion. Effective rheological model parameters are flexural rigidity D, a measure of the lithosphere's elastic strength, and viscous relaxation time constant τ, a measure of its viscosity. Erosion is assumed to occur at a rate linearly proportional to the height of topography at any given time, the proportionality being determined by the erosion time constant ∂. The models are constructed in space and time frequency domains using simple linear systems theory. Two loading schemes are considered. Model 1 assumes that a flexurally competent lithosphere is suddenly loaded by topography which subsequently erodes whereas model 2 ignores the emplacement and erosion of topography during orogenesis when the lithosphere is thermally activated and considers only the load effected by erosion beginning some time later when the lithosphere has cooled and strengthened. Both models predict that the remnant topography of very old regions such as the Canadian Shield can be maintained, in fact is favoured, by a viscoelastic lithosphere with τ as small as 1-10Ma. The continued existence of such topography does not necessarily imply significant long-term lithospheric strength. Model 1 results indicate low values of the flexural rigidity D, ∼ 1021-1022Nm, which are interpreted to characterize not so much presentday D but that which was effective when topography was formed. Model 2 more successfully reproduces the data and constrains the product ∼D to values, ∼ 1024-1025Nm Ma, which compare very favourably to those found from studies of individual loads on continental lithosphere. The model results show that the erosion time constant τ is probably wavenumber-dependent and falls in the range 200–400 Ma for the topographic wavelengths under consideration.

Thermal contraction and flexure of intracratonal basins: a three-dimensional study of the Michigan basin

Cooling of continental lithosphere by thermal conduction to the surface following a heating event imposes a load on the lithosphere and causes subsidence as the basement rocks contract. Accumulation of sediments in the resulting depression forms a sedimentary basin. Studies of the geometry of sedimentary basins with horizontal dimensions of a few hundred kilometres suggests the lithosphere responds to loads by regional flexure of a strong elastic or viscoelastic upper lithosphere. Two- and three-dimensional models for flexure of the lithosphere owing to thermal contraction are applied to the Michigan basin. 2-D flexural models are compatible with previously estimated rheological parameters. For an elastic lithosphere, an effective flexural rigidity of about 5 × 1029 dyne cm can explain observed gravity anomalies and structural contours across the Michigan basin. A viscoelastic lithosphere requires an effective flexural rigidity of about 1032 dyne cm for a viscosity of 1025 poise. For a 3-D elastic or viscoelastic lithosphere, effective flexural rigidity must be approximately an order of magnitude smaller: 2 × 1028 and 1031 dyne cm, respectively. Larger values leave an excessive remaining load at the centre of the basin, which is not indicated by observed gravity anomalies. The decrease in effective flexural rigidity is attributed to greater difficulty in equi-dimensional bending of a 3-D plate. The magnitude and spatial distribution of the computed thermal contraction load are consistent with the presumed mechanism. As the thermal anomaly of the initial heating event was present mainly at great depths in the lithosphere, the lack of geological evidence for an initial heating event prior to subsidence is not a fatal objection to the thermal contraction mechanism of subsidence. A driving load, presumably from subsurface processes including crustal stretching and emplacement of dense rocks into the crust, is necessary to explain the net subsidence of the basin from the time before the heating event until the lithosphere beneath the basin cooled. Without this driving load, the thermal expansion would produce uplift, but the surface would subside only back to sea-level after the lithosphere cooled. Theoretical gravity results indicate that the driving load is centred at a depth of approximately 15 km. Deviations in subsidence curves from exponentials associated with thermal contraction can be explained by changes in sediment supply. Spatial variations in sediment load, caused by regional facies changes, produce migration of the centre of maximum deposition. Water and basement depths are determined for sequential time intervals during basin development. The predicted depositional environments are consistent with lithofacies maps of the Middle Devonian.

Intraplate lithosphere strength and heat flow

The magnitude levels of intraplate stresses and the strength of intraplate lithosphere are critical for understanding the mechanisms of intraplate deformation. The decay of stress in the ductile lower lithosphere and its resulting amplification in the brittle upper lithosphere constitute an important mechanism for providing stress levels in the upper brittle lithosphere sufficient for fracture. Such stress transfer and the subsequent lithosphere deformation are controlled by lithosphere rheology. A model is examined here for the deformation of viscoelastic lithosphere, with temperature- and stress-dependent rheology, subjected to laterally applied stress for different lithosphere thermal structures. The geotherm is found to be of critical importance in determining the strength of the lithosphere and the extent of the resulting lithosphere deformation. For lithosphere with heat flow ⩾75 mW m−2, 0.25 kbar tensile stress is sufficient to cause geologically significant deformation, while for compressional deformation, either greater stress levels or a hotter lithosphere are necessary.

Lithosphere response to externally and internally derived stresses: a viscoelastic stress guide with amplification

Summary Lithosphere subjected to an externally derived horizontal stress undergoes creep in the lower lithosphere resulting in the decay of lower lithosphere stress and the associated amplification of stress within the upper lithosphere. This stress response of lithosphere has been investigated for a lithosphere model with power-law stress and temperature dependent viscoelastic properties. The rate and extent of stress decay and associated stress amplification is greatly dependent on the lithosphere geotherm. Oceanic lithosphere subject to an applied stress of ± O.1 kb undergoes upper lithosphere stress amplification of × 1.5, × 1.8 and × 2.0 at 104, 106 and 108 yr respectively. At 106 yr the effective lithosphere thickness is reduced to approximately 40 km. The stress decay and amplification proceeds more rapidly for an applied stress of ± 1.0 kb and in the case of a tensile applied stress results in some upper lithosphere fracture. For continental and Basin and Range type lithosphere the comparable stress amplification at 106 yr, for a stress of ± 0.1 kb, is × 2.0 and × 6.5 with effective lithosphere thicknesses of 60 and 20 km respectively. The large values of stress amplification for the Basin and Range lithosphere result in complete upper lithosphere fracture which gives rise to a cyclic process of upper lithosphere faulting and lower lithosphere creep in which extensive lithosphere deformation can occur. The stress amplification process also occurs for stresses generated by lateral density contrasts. For an isostatically compensated plateau uplift structure, deviatoric stresses of the order of 1 kb can be generated in the upper lithosphere by this process and are sufficient to cause tensile fracture of the upper lithosphere. The response of viscoelastic lithosphere to constant geometry bending stresses has also been examined and results in substantial but not complete reduction of the bending stresses.