Viscoelastic Earth Model Research Articles

Time‐domain approach for the transient responses in stratified viscoelastic Earth models

We have developed the numerical algorithm for the computation of transient viscoelastic responses in the time domain for a radially stratified Earth model. Stratifications in both the elastic parameters and the viscosity profile have been considered. The particular viscosity profile employed has a viscosity maximum with a contrast of O(10²) in the mid lower mantle. The distribution of relaxation times reveals the presence of a continuous spectrum situated between O(10²) and O(104) years. The principal mode is embedded within this continuous spectrum. From this initial‐value approach we have found that for the low degree harmonics the non‐modal contributions are comparable to the modal contributions. For this viscosity model the differences between the time‐domain and normal‐mode results are found to decrease strongly with increasing angular order. These calculations also show that a time‐dependent effective relaxation time can be defined, which can be bounded by the relaxation times of the principal modes.

Pleistocene glaciation and the Earth's precession constant

We predict time-dependent variations in the Earth's precession constant arising from the ice and ocean mass fluctuations associated with the Late Pleistocene glacial cycles. Our predictions incorporate contributions from both the surface mass load redistribution and the adjustment of the solid earth. In this regard, we adopt spherically symmetrical, self-gravitating, Maxwell viscoelastic earth models and obtain results for a large suite of radial viscosity profiles. These profiles include a set obtained from published inferences based on post-glacial relative sea-level and uplift histories, as well as a set intended to sample the sensitivity of the results to variations in the viscosity within a number of major subregions of the mantle (e.g. The transition zone, the upper mantle, and the lower mantle). A more detailed measure of this sensitivity is also obtained by computing Frechet kernels for the predictions. We construct an ice model which incorporates the ICE-3G model for the final deglaciation event and which is constrained to yield a eustatic sea-level variation which matches observed fluctuations in oxygen isotope records over the last 800 kyr. In all cases, the ocean mass redistribution is constrained to be gravitationally self-consistent. The computed Frechet kernels indicate that the predictions are most sensitive to variations in viscosity in the deepest regions of the mantle; indeed, in some cases the sensitivity peaks at the core-mantle boundary. Both positive and negative perturbations to the precession constant are predicted, with the maximum peak-to-peak (relative) variation being ∼0.20 per cent for the published viscosity models and ∼0.32 per cent for all other models. Furthermore, the mean relative perturbation in the precession constant, with respect to the present-day value, is found to reach ∼ -0.08 per cent for the published viscosity models, and ∼ -0.20 per cent for other models. Recent solutions for the Earth's precession, obliquity and insolation parameters (Laskar, Joutel & Boudin 1993), indicate a passage through resonance, associated with a perturbation of Jupiter and Saturn, in the case when the mean relative perturbation in the precession constant is ∼-0.15 per cent. We find that this threshold is not achieved for any of the published viscosity models; however, it is reached, for example, for earth models with lower mantle viscosities which exceed 30–50 × 1021 Pa s, or models characterized by large (∼ two orders of magnitude) jumps in viscosity at mid-lower mantle depths.

Some comments on the 3-D impulse response of a Maxwell viscoelastic earth

SUMMARY We consider the 3-D normal mode response of a Maxwell viscoelastic earth model. We find that the partitioning of modal strengths, amongst the various modal branches, is significantly different for the case of the l Love nunber, which governs tangential deformations, than it is for the h and k Love number, which are associated, respectively, with radial deformations and gravitational potential perturbations. The partitioning and the differences are a strong function of the spherical harmonic degree.

Late Pleistocene and Holocene ice mass fluctuations and the Earth's precession constant

We computed, using a set of previously published Maxwell visco-elastic Earth models, time-dependent variations in the Earth's precession constant, H, associated with glacial cycles of duration ⋍ 10 5 yr. We found peak perturbations which are approximately an order of magnitude less than previous upper bound estimates, which assumed that variations in H arose solely from the redistribution of the ice age surface mass (ice plus ocean) load; that is, these studies neglected the visco-elastic isostatic adjustment of the solid planet and assumed a rigid Earth. A number of subsequent studies have cited and/or adopted the ‘rigid Earth’ variation in H to argue for a series of potentially important effects on the Earth's orbital dynamics, rotation and palaeoclimate. The results presented herein suggest that these arguments require significant re-appraisal. We also compute variations in H associated with the documented retreat of small ice sheets and glaciers over the last century and potential (recent) mass variations in the Antarctic and Greenland ice complexes. We find (in the case of a net mass loss from the Antarctic and Greenland ice sheets) that these effects may contribute a perturbation to H which is approximately the same amplitude as, but has the opposite sign to, the contribution from ongoing glacial isostatic adjustment during the last century.

A singularity free approach to post glacial rebound calculations

Calculating the post glacial response of a viscoelastic Earth model using the exponential decay normal mode technique leads to intrinsic singularities if viscosity varies continuously as a function of radius. We develop a numerical implementation of the Complex Real Fourier transform (CRFT) method as an accurate and stable procedure to avoid these singularities. Using CRFT, we investigate the response of a set of Maxwell Earth models to surface loading. We find that the effect of expanding a layered viscosity structure into a continuously varying structure is to destroy the modes associated with the boundary between layers. Horizontal motion is more sensitive than vertical motion to the viscosity structure just below the lithosphere. Horizontal motion is less sensitive to the viscosity of the lower mantle than the vertical motion is. When the viscosity increases at 670 km depth by a factor of about 60, the response of the lower mantle is close to its elastic limit. Any further increase of the viscosity contrast at 670 km depth or further increase of viscosity as a continuous function of depth starting from 670 km depth is unlikely to be resolved.

A spectral formalism for computing three‐dimensional deformations due to surface loads: 2. Present‐day glacial isostatic adjustment

Using a spherically symmetric, self‐gravitating, linear viscoelastic Earth model, we predict present‐day three‐dimensional surface deformation rates and baseline evolutions arising as a consequence of the late Pleistocene glacial cycles. In general, we use realistic models for the space‐time geometry of the final late Pleistocene deglaciation event and incorporate a gravitationally self‐consistent ocean meltwater redistribution. The predictions of horizontal velocity presented herein differ significantly, in both their amplitude and their spatial variation, from those presented in an earlier analysis of others which adopted simplified models of both the late Pleistocene ice history and the Earth rheology. An important characteristic of our predicted velocity fields is that the melting of the Laurentide ice sheet over Canada is capable of contributing appreciably to the adjustment in Europe. The sensitivity of the predictions to variations in mantle rheology is investigated by considering a number of different Earth models, and by computing appropriate Frechet kernels. These calculations suggest that the sensitivity of the deformations to the Earth's rheology is significant and strongly dependent on the location of the site relative to the ancient ice sheet. The effects on the predictions of three‐dimensional deformation rates of altering the ice history or adopting approximate models for the ocean meltwater redistribution have also been considered and found to be important (the former especially so). Finally, for a suite of Earth models we provide predictions of the velocity of a number of baselines in North America and Europe. We find that, in general, both radial and tangential motions contribute significantly to baseline length changes, and that these contributions are a strong function of the Earth model. We have, furthermore, found a set of Earth models which, together with the ICE‐3G deglaciation chronology, produce predictions of baseline length changes that are consistent with very long baseline interferometry measurements of baselines within Europe.

A comparison of methods for the inversion of viscoelastic relaxation spectra

SUMMARY The majority of investigations of the glacial isostatic adjustment problem have proceeded by invoking the correspondence principle and solving for the (Laplace) transformed impulse response of the viscoelastic Earth model that is represented in terms of Love number spectra (Peltier 1974). This formulation requires a final inversion of the solution into the time domain, and the present paper is concerned with a comparison and assessment of the three techniques (pure collocation, full normal mode analysis using residue theory, and a hybrid technique which we term mixed collocation) that have been developed to perform it. On the basis of the analysis presented here we conclude that both the full normal mode analysis and mixed collocation can generate accurate inversions of the Love number spectra. We also derive clear guidelines on the choice of collocation points that ensure that the same accuracy is achieved using pure collocation. As a final point we stress anew that, regardless of the technique employed, the accuracy of the inversions can and should be checked by comparing the predicted infinite time-scale response for a Heaviside loading history with an independent calculation of the response for an inviscid Earth with a lithosphere of appropriate thickness.

Lower-mantle viscosity constrained by seismicity around deglaciated regions

KNOWLEDGE of the viscosity structure of the Earth's mantle is important for constraining models of mantle convection and isostatic rebound. Here we show that seismicity around the margins of deglaciated areas provides a constraint on the viscosity of the lower mantle, in addition to those previously proposed1,2. Calculations using a spherical, viscoelastic Earth model show that the present-day magnitude of the stress fields induced in the lithosphere beneath the (now-disappeared) Laurentide and Fennoscandian ice sheets is very sensitive to the value of the lower-mantle viscosity. Stresses of ∼100 bar, sufficient to cause seismicity, can still remain in the lithosphere for lower-mantle viscosities greater than ∼1022 Pa s; for lower-mantle viscosities of ∼1021 Pa s, only a few tens of bars of stress persist in the lithosphere today. This influence of lower-mantle viscosity on the state of stress in the lithosphere also has implications for the migration of stress from earthquakes, and hence for earthquake recurrence times.

Computation of synthetic seismograms for stratified azimuthally anisotropic media

We outline a method of computing synthetic seismograms for stratified, azimuthally anisotropic, viscoelastic earth models. This method is an extended form of the Kennett algorithm that is efficient for multioffset vertical seismic profiling. The model consists of a stack of homogeneous plane layers, and the response is computed iteratively by successive inclusion of deeper layers. In each layer, the 6×6 system matrix A is diagonalized numerically; this permits treatment of triclinic materials, i.e., those with the lowest possible symmetry. Jacobi iteration is an efficient way to diagonalize A because the entries of A change little from one wavenumber to the next. When the material properties are frequency dependent, the wavenumber loops are inside the frequency loop, and the computation is slow even on a supercomputer. When the material parameters are frequency independent, it is better to make frequency the deepest loop, with diagonalization of A outside the loop, in which case vectorization gives a relatively rapid computation. Temporal wraparound is avoided by making use of complex frequencies, and spatial aliasing is avoided by using a generalized Filon's method to evaluate both the wavenumber integrals. Various methods of generating anisotropic elastic constants from microlayers, cracks, and fractures and joints are discussed. Example computations are given for azimuthally isotropic and azimuthally anisotropic (AA) earth models. Comparison of computations using single and double wavenumber integrations for a realistic AA model shows that single wave‐number integration often gives incorrect answers especially at near offsets. Errors due to use of a single wavenumber integration are explained heuristically by use of wave front diagrams for point and line sources.

Deformation of the Viscoelastic Earth Due to the Secular Change In the Earth's Axis of Rotation

The cause of the secular change in the Earth's rotation axis relative to the surface geography is regarded as the viscoelastic response of the earth to the mass displacement generated by late Pleistocene deglaciation and the sea level changes related to it. On the other hand, the secular drift of the pole will also bring about viscoelastic change in the Earth's shape because the Earth is not a sphere but a spheroid with an equatorial bulge. In the case of polar drift, we estimated the deformation of the Earth's surface and the Earth's core, as well as the change in the Earth's gravity field. We used a compressible, radially stratified, Maxwell viscoelastic earth model for the inverse calculations and also investigated if the viscous flow in the mantle can compensate for the mass redistribution caused by polar drift. As a result of these calculations, the viscous flow in the mantle seems not to be able to compensate for the mass redistribution caused by polar drift if the mantle viscosity is higher than 1022 poise. the agreement or disagreement between the rotation axis of the Earth and the principal axis of inertia determined by satellites can add a new constraint to the viscosity of the Earth's mantle.

Quasi-static deformation due to a dislocation source in a Maxwellian viscoelastic earth model.

In this paper, we treat quasi-static deformation field due to a dislocation source in a composite medium which consists of elastic layers overlying a strat-ified viscoelastic half-space. The rheologies of the viscoelastic layers are as-sumed to be elastic dilatational and Maxwell deviatoric. Integral representa-tions of the surface displacements are derived from those in the associated elastic problem by applying the correspondence principle of linear viscoelastici-ty. The effect of gravity is taken into account by modifying the boundary con-ditions with respect to normal stress. Features of the viscoelastic displacements are investigated for a three-layered structure model composed of an elastic surface layer, an intervening low viscosity layer, and a viscoelastic substratum with relatively high viscosity. A dislocation source is supposed to be located in the elastic surface layer. The effect of the stress relaxation in the substratum is examined for both a point source and a dimensional fault. Difference in the viscosity of the substratum affects the amount and extent of the viscoelastic deformation. In a model with a low viscous substratum, the amount of displacement is large and the deformation spreads in a broad area. The effect of gravity is investigated assuming a pure dip-slip fault. The viscoelastic deformation for the gravitating case is restrained in amount and extent as compared with the non-gravitating case. Such an effect of gravity becomes notable with time. In the gravitating case, the stress relaxation progresses fast, and the apparent relaxation time of deformation is much shorter than that for the non-gravitating case.

The periodic and eigen solutions of the equation governing a nonrotating viscoelastic earth model under transient force

The problem of existence of the periodic solution of the equation governing a nonrotating viscoelastic earth model under transient force is examined. By first formulating the governing equations, using the methods of Coleman and Noll ( Rev. Modern Physics 33 (2) (1961) , 239–249), Dahlen and Smith ( Philos. Trans. Roy. Soc. London A 279 (1975) , 583–624), and Biot (“Mechanics of Incremental Deformations,” Wiley, New York, 1965 ), these equations are subjected to oscillatory displacement resulting in an eigenvalue problem whose solutions are the viscoelastic-gravitational displacement eigenfunctions U ( x ) with associated eigenfrequencies ω. A theorem is then proved to show the existence of a periodic solution.

The impulse response of a Maxwell Earth

An extended form of the correspondence principle is employed to determine directly the quasi‐static deformation of viscoelastic earth models by mass loads applied to the surface. The stress‐strain relation employed is that appropriate to a Maxwell medium. Most emphasis is placed on the discussion of spherically stratified self‐gravitating earth models, although some consideration is given to the uniform elastic half space and to the uniform viscous sphere, since they determine certain limiting behaviors that are useful for interpretation and proper normalization of the general problem. Laplace transform domain solutions are obtained in the form of ‘s spectra’ of a set of viscoelastic Love numbers. These Love numbers are defined in analogy with the equivalent elastic problem. An efficient technique is described for the inversion of these s spectra, and this technique is employed to produce sets of time dependent Love numbers for a series of illustrative earth models. These sets of time dependent Love numbers are combined to produce Green functions for the surface mass load boundary value problem. Through these impulse response functions, which are obtained for radial displacement, gravity anomaly, and tilt, a brief discussion is given of the approach to isostatic equilibrium. The response of the earth to an arbitrary quasi‐static surface loading may be determined by evaluating a space‐time convolution integral over the loaded region using these response functions.

Stress relaxation in a semi-infinite viscoelastic earth model : ROSENMAN, M HARVARD UNIV. CAMBRIDGE, MASS. USA SINGH, SJ MIT, CAMBRIDGE, MASSACHUSETTS, USA 6F, 4R. SEISMOL. SOC. AMER. BULL. V63, N6, 1973, P2145–2148

Stress relaxation in a semi-infinite viscoelastic earth model : ROSENMAN, M HARVARD UNIV. CAMBRIDGE, MASS. USA SINGH, SJ MIT, CAMBRIDGE, MASSACHUSETTS, USA 6F, 4R. SEISMOL. SOC. AMER. BULL. V63, N6, 1973, P2145–2148

Stress relaxation in a semi-infinite viscoelastic earth model

Abstract Expressions for quasi-static surface stresses resulting from a finite, rectangular, vertical, strike-slip fault in a Maxwellian viscoelastic half-space are derived. Variation of the stresses with time and epicentral distance is studied. Contour maps are obtained in some representative cases. It is found that all nonvanishing stress components at the free surface die exponentially with time. This is in contrast to the behavior of the displacements and strains which, in general, do not vanish for large times.