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.


INTRODUCTION
Density anomalies are widely recognized to be important sources of stress field in the earth for a large range of wavelengths (Fleitout & Froidevaux 1982;Hager 1983;Richards & Hager 1984;Ricard, Fleitout & Froidevaux 1984). Mathematical models consisting of elastic plates overlying an inviscid fluid made it possible to estimate the shear stress induced in the upper mantle by crustal and lithospheric density anomalies (Caputo, Milana & Rayhorn 1984;Mareschal & Kuang 1986). It has been found that high stress levels of the order of several hundred bars can be reached in the elastic portion of the lithosphere as a consequence of lateral density variations. This forcing mechanism has thus been proposed as a valuable contributor to intraplate seismicity.
There are, on the other hand, some problems that remain unsolved if a rheology for the lithosphere is not taken into account in the modelling. The laterally varying density structure, usually inverted by gravimetric and seismic methods, is a consequence of tectonic activity that took several million years to deform the crust and lithosphere to its actual shape. On this time-scale, the stress induced by density anomalies is probably released by recurring earthquakes or plastic processes in the deep lithosphere (Bonafede, Dragoni & Morelli 1986). This leads to the conclusion that stress values inferred within the framework of models that do not properly account for lithospheric rheology have to be interpreted as upper bounds, with the implicit assumption that no releasing mechanism has been active during the deformation of the density structures.
To avoid the ambiguities suffered by the interpretation of the stress field, we built a stratified viscoelastic lithosphere to concentrate our attention on the rate of accumulation of shear stress and vertical ground motions driven by density anomalies. A realistic time history for the source mechanism is considered. A layered lithosphere with a temperaturedependent Maxwell rheology was first modelled by Quinlan & Beaumont (1984) to intepret the Appalachian Basin as a foreland basin formed under the loads of successive overthrusts. These authors studied the deflection of the lithosphere under the loads on timescales of lo2 Myr.
The behaviour of a perfectly elastic plate for the evolution of a foreland basin was addressed in the paper by Jordan (1981). Horizontal tectonic forces are not considered in the model. Cloetingh, McQueen & Lambeck (1985) have shown that discontinuities in horizontal tension or compression in a purely elastic lithosphere can be responsible for a rate of vertical motions of 10-*-10-' mm yr-' in the periphery of vertical applied loads. This amount of deformation is at least an order of magnitude lower than our estimate of vertical motion forced by density anomalies. Cloetingh (1986) suggested, on the other hand, that the introduction of a depth dependent (Goetze & Evans 1979) or viscoelastic rheology in his modelling, would significantly enhance the values of vertical motions. In general, the effects of density anomalies studied in this paper have to be superimposed on those induced by variations in the horizontal stress.
If the interaction between recurring earthquakes is neglected, as seems quite reasonable to a first order approximation, rate of stress accumulation can be used to estimate the recurrence time of earthquakes in a seismic area to test the efficiency of the proposed mechanism in triggering the recorded seismicity. Plastic processes in the crust can now be neglected since we consider the stress that is accumulated in a few thousands of years and, on this time scale, it is an excellent approximation to assume an elastic top layer as a model for the crust. At the same time our model can be used to make predictions on the rates of vertical deformation forced by density anomalies whose lateral extension is a function of time. We show how our results could be compared with geodetic measurements to test different tectonic models.
There is another crucial reason that suggests the need of a rheology for the lithosphere. Geological information indicates that crustal deformation is not a steady-state process. Crustal shortening or thinning occurs in a sequence of tectonic phases that push the crust-lithosphere system out of isostatic equilibrium. This is certainly true for the Mediterranean region in which we are mainly interested. From the results obtained during the Oceanic Drilling Program (ODP Leg 107) and previous geological and geophysical data, it is now recognized that crustal thinning associated with the opening of the Tyrrhenian Sea is a process that occurred in several steps (Sartori 1988). Tectonic activity in the Apennines is, on the other hand, responsible for the development of strong laterally varying structures (Baldi et al. 1982). The discontinuity in time that we find in the formation of the Tyrrhenian basin and surrounding chains seems to be a common feature of collision zones (Wortel & Cloetingh 1986). If the lithospheric viscosity is sufficiently high such that the memory of the system is comparable with the time scales of the tectonic phases, it becomes necessary to take into account in a proper way the rheology of the lithosphere. We built a stratified, spherical, viscoelastic earth to estimate the effects of sphericity, gravitation and selfgravitation on the flexural properties of the lithosphere when the planet is forced by lithospheric density anomalies. We find indeed that sphericity plays a crucial role in the flexure of the lithosphere for elongated crustal density structures. The model is incompressible and the explicit expressions for the inverse of the fundamental matrix entering the solutions for both the gravitating and self-gravitating cases are provided in the paper. The formalism developed can be easily adapted to any kind of source mechanism that perturbs a spherical earth model. Density contrasts or dislocations embedded at varying depth both in the elastic and viscoelastic layers and tidal forcing can be accounted for as driving mechanisms in the theory.
We show several cases of general interest that can be thought to be appropriate also for the tectonic configuration of central Italy (Sections 3 and 4). In Section 6 we make use of the crustal density structures inverted from Bouguer gravity data, the observed topography and estimated sedimentation in the Adriatic foredeep to evaluate the expected rate of vertical deformation along a transect perpendicular to the Apennines. Geological information on the velocity of overthrusting or timing of the various tectonic phases is taken into account.

MATHEMATICAL MODEL: T H E EFFECTS OF SPHERICITY, G R A V I T A T I O N A N D SELF-GRAVITATION
We develop the mathematical formalism that is needed to study the deformation and stress field forced in a gravitating and self-gravitating viscoelastic earth by excitation sources buried at varying depths. Our spherical model is differentiated into an elastic outer crust and inner viscoelastic Maxwell solid that reproduces the ductile properties of the lithosphere on geological time scales. The term 'lithospheric viscosity' is used throughout for the parameter that characterizes the Maxwell rheology of the inner layer. In our two-layer model we do not account for the viscosity contrast between the deep lithosphere and underlying mantle. Our mathematical description of the mechanical properties of the lithosphere provides a simplified version of a somewhat more complicated depth-dependent rheology. Sophisticated non-linear constitutive relations, based on the mechanism of flow in olivine, could be considered in these studies (Goetze 1978;Goetze  & Evans 1979). We think, on the other hand, that a simple Maxwell model clearly enlightens the basic physics of the deformation induced by density anomalies, with the advantage of being mathematically convenient.
In Table 1 we provide the reference model in which the thickness of the crust is fixed at 30 km. The parameters defining this model are also appropriate for the reference configuration in the Mediterranean region (Baldi et 01. 1982). The thickness of the crust is modified in a sequence of calculations to show the consequences of variations in the flexural properties of the elastic layer.
We simulate a realistic history for the density anomalies whose lateral extension is allowed to vary as a function of time where u is the displacement vector and r is the stress tensor appropriate for Maxwell rheology (Fung 1965). p denotes the background density field, g ( r ) the associated gravitational acceleration and f is the body force due to density anomalies. The solution obtained for a gravitating model will be generalized to take into account the effects of self-gravitation. Equation (1) is linearized with respect to perturbations from a reference configuration in which the stress is due to hydrostatic pressure (Rayleigh 1906;Dahlen 1971 where p is the shear modulus and v is the long-term viscosity. For a spherically symmetric earth it is convenient to expand the displacement field and stress tensor components in spherical harmonics. We introduce the solution vector y where Ul, V, denote the spheroidal part of the radial and tangential components of the displacement field expansion (Sabadini, Yuen & Boschi 1984). T', and TLo are the components of the radial and tangential stress, respectively. We get from equation (1) the following fourth-order system of ordinary differential equations for each degree 1 where r is the radial distance from the centre of the Earth and F denotes the component of the equivalent body-force due to the source. For a surface density anomaly embedded at depth r, we have (Hager & O'Connell 1981) Volume density anomalies are changed into surface density anomalies after integration in the radial direction.
The elements of the A matrix for a gravitating Earth can be obtained from equations (11)-(16) in Sabadini, Yuen & Boschi (1984), dropping the terms that are 'due to self-gravitation. We can express the solution of equation (4) in terms of the fundamental matrix Y as where yb denotes the solution at the crust-lithosphere interface b and a is the Earth radius. The explicit expression of the fundamental matrix Y is given in equations (46) and (48) in Sabadini et al. (1982). The inverse of the fundamental matrix is also needed if we want an explicit, self-consistent and analytical expression of the vector solution. In the Appendix (equations A5-A8) we provide the inverse Y -' ( r ) , explicitly given as a function of the density, rigidity and radial distance r. This allows us to handle analytically any radial dependence in the forcing term. Convolution of the Green function (A12) with the appropriate surface density anomaly embedded at depth r,, allows us to get the radial component of the displacement ~( e , S) = 2na2 (7) where a(@) denotes the angular dependence of the surface density anomaly and f(s) provides the time dependence of the source in the Laplace transform domain. We summed up to I = 1000 to ensure convergence in the solutions with lateral dimensions of the source as small as 0.5". The time derivative of the radial displacement gives the rate of vertical deformation that is shown in the paper.
In this work the forcing is produced by rings of surface density anomalies buried at varying depth and parallel to the equator of the model. The spherical expansion of a ring of surface density anomaly can be derived from equation (16) in Sabadini & Peltier (1981), by simply adding the contributions of two caps of opposite sign.
We point out that our model accounts automatically for the negative mass anomaly that is induced by a load displaced horizontally without variations in the lateral dimension. This is due to the self-consistent way in which gravitational forces and rheology are introduced. On the time scale controlled by the horizontal velocity of tectonic transport and relaxation in the lithosphere, the induced anomaly is compensated.
We now study the role played by sphericity, gravitation and self-gravitation on lithospheric deformation forced by surface density anomalies. The effects of sphericity are estimated from the comparison between our model and a flat approximation. Gravitation and self-gravitation pertain to the realm of the spherical earth. In Fig. 1 we show a section of the model with portions of the density anomalies extending also in the lateral direction. We study the 2-D deformation problem in the dotted plane containing the axis of symmetry of the model and perpendicular to the density structure. The crust is bordered by solid lines. The deformation pattern is symmetric around the polar axis, with two maxima located in the regions where the plane intersects the ring. Due to the symmetry of the problem, we consider only one of these regions. The ring of surface density anomaly is the spherical counterpart of an infinite strip loading a flat earth. Our approximation is thus valid for elongated lithospheric density structures and topography. The proposed modelling can also handle realistic finite sources, but these are beyond the purposes of the paper.
Throughout this work the lithospheric viscosity is varied from 10" to 1023Pas, in agreement with commonly accepted values for this parameter. Hager & Clayton (1987) built a global dynamic earth model in which the lithosphere has a viscosity of 10" Pa s. The range of values that is used in our model agrees well with the estimated effective viscosity for olivine under upper mantle conditions, at depth greater than 20 km and stress of the order of some hundred bars (O'Connell & Hager 1980). A lithospheric viscosity of around loz3 Pa s has been suggested by Walcott (1970) from his analysis of the Earth response to surface loads. A value of 5 X loz2 P a s for the lithosphere has been employed by Bonafede et al. (1986) in their modelling of recurring earthquakes on long transcurrent faults.
We first consider a simple time history for the source corresponding to a load of 1OOOm height and 2" lateral extension. The model parameters are given in Table 1 and the lithospheric viscosity is ld2 Pa s. The density anomaly is 'turned on' abruptly according to a box shaped time history. The rate of deformation evaluated at the initial stage of the motion is plotted in Fig. 2 for varying distance of the strip from the north pole of the model, which is the intersection point of the axis of symmetry with the surface. The time history is unrealistic but it can be useful to put bounds on the effects of sphericity when comparison is made with predictions derived from flat models. We point out that this time history is only used in Figs 2 and 3. For the flat earth we used the incompressible limit of Landau & Lifshitz (1967) solution, generalized to a viscoelastic halfspace by means of the Correspondence Principle. In the top panel the dashed curve corresponds to the flat model and the solid curve to the spherical one. Gravitational restoring forces are not included in this panel where attention is concentrated on the effects of sphericity. The values of the rate are evaluated at the centre of the strip. The anomaly degenerates into a cap when the radius of the ring is reduced. In this limit the results of the spherical and flat models merge together. Enlarging the amplitude of the ring, the flat model predicts an increase in the rate of deformation. Long wavelength deformation properties of the spherical earth act to reduce the rate of deformation with respect to the flat model. When the ring is located at the equator of the sphere and corresponding distance in the flat model, we have a factor two reduction in the predictions of the spherical earth with respect to the half-space. This is an upper bound in the estimate of spherical effects due to the continuous structure of the source, but clearly demonstrates the potential role migration rate of the boundary between the crust and lithosphere of different density is higher. For longer times the difference is reduced to a few per cent. For higher viscosity the deviation between gravitating and selfgravitating models is maintained longer. The two predictions become essentially identical after a sufficiently long interval of time.
We find that self-gravitation plays a minor role in the ground motion forced in a crust-lithosphere system by density anomalies.
-3.0 played by sphericity in global tectonics. This can have an impact on deformation characteristics of subduction zones extending in the horizontal direction for at least several hundred kilometers.
Within the framework of the spherical model, we now study in Fig. 2 (bottom panel) the effects of gravitation and self-gravitation. In the gravitating model we provide a constant gravity field that is responsible for the restoring force. Self-gravitation accounts for the effects of perturbation in the gravitational potential. The dotted curve reproduces the spherical case of the top panel while the solid curve corresponds to the gravitating model. Restoring forces reduce by around 30 per cent the deformation in the early stage of the motion. Self-gravitation is responsible for a further 15 per cent reduction in comparison with the gravitating case.
The effects of perturbations in the gravitational potential are better understood if the time evolution of the rate of deformation is considered. In Fig. 3 we show the rate of subsidence when the strip of density anomaly is located at the equator of the model to minimize boundary effects. The rate u is plotted as a function of the distance from the centre of the strip. The deformation is symmetric around the vertical axis. Solid curves correspond to the gravitating model, while dashed ones include self-gravitation. A lithospheric viscosity of 10' ' Pa s has been used for the top panel. At the bottom this parameter has been increased to 1023Pas. The parameter A provides information on the interval of time elapsed since the beginning of the loading event. The largest deviation between the two models occurs during the early stage of the deformation ( A = 0) when the Comparison between the gravitating model (solid) and the self-gravitating one (dashed). The rate of subsidence, driven by the same source of Fig. 2, is plotted as a function of the distance from the centre of the anomaly. A gives the interval elapsed from the loading event at r = 0. Panel (a) corresponds to v = 10" Pas. Lithospheric viscosity is increased to loz3 Pas in the panel (b).

A REALISTIC TIME HISTORY FOR DENSITY ANOMALIES: CONSEQUENCES FOR VERTICAL MOTIONS FORCED BY SURFACE T O P O G R A P H Y
When a rheology for the lithosphere is assumed, it becomes necessary to specify correctly the time history of density anomalies. If the lithosphere is described as an inviscid fluid it is possible to study only a configuration of equilibrium in which the density anomalies are compensated. In this case, the evolution in time of the anomalies does not enter play.
The validity of the model results depends on the assumption that the loading conditions did not change in the time scale that is necessary for the lithosphere to deform. This extreme configuration is certainly unrealistic as a consequence of the discontinuity in time of tectonic phases (Wortel & Cloetingh, 1986). In active areas the crust lithosphere system can hardly reach a configuration of isostatic equilibrium, being perturbed on time scales comparable with its relaxation time. In order to get a better picture of the real configuration of the system, it thus becomes necessary to match the definition of a rheology with a realistic time history for the density structure.
In this section we study the consequences on vertical ground deformation produced by a surface load that is expanding in the horizontal direction. Our aim is to model the overthrusting of a mountain range on the crust. We show several cases of general interest, but the parameter values that we have chosen are also appropriate to study the overthrusting of the Apennines on the Adriatic platform. A more detailed description of the tectonic setting of the Apenninic region is given in Sections. 5 and 6. In what follows, the lateral expansion of the surface load is the continuous counterpart of the emplacement of overthrusts on the lithosphere as given in the paper by Quinlan & Beaumont (1984). Variations in time of the lateral dimensions of the anomalies are obtained by the interference of elementary contributions whose history is a simple box in time domain. The angular and s-dependent parts of surface density anomalies entering the convolution (7) are now given by N a(e, S) = C a(en+, -en) exp (-n A~s ) / s , where N denotes the number of contributions that need to be summed in order to ensure the convergence of the solution. Each contribution extending between colatitudes 8, and On+, is added after an interval of time given by n At, as portrayed by the s-dependent terms. The parameters Ar and N are chosen in concert to reproduce the desired variation in time of the lateral extension. The number N of contributions that are needed is viscosity dependent. The lower the viscosity, the higher the number of steps that are necessary for convergence. For realistic time histories, corresponding to an overthrusting velocity ranging from 1 to 10 cm yr-', at least 40 steps are found to be necessary when the lower bound ld2 P a s and a lateral extension of a few degrees are considered. We fixed N = 40 in our calculations to ensure the quality of the solution in the range of considered viscosities.
In Fig. 4 we show the rate of subsidence due to a surface load of 1OOOm height that is expanding to the right in the horizontal direction at different velocities. The box in the top panel indicates the position of the load and the arrow is used to emphasize the motion in the horizontal direction.
The position of the front of the overthrusting load, at the instant in which the deformation is evaluated, is provided by the right border of the dashed portion of the box. The lateral extension of the load is 2" when the deformation is recorded; this could be appropriate for the Apennines.
The deformation pattern is characterized by a maximum subsidence localized just beneath the migrating front of the load. Curves 1, 2 and 3 in both panels correspond to horizontal velocity V of 10, 5 and 1 cm yr-', respectively. A rate of subsidence of around 0.7mmyr-' is produced by V = 10 cm yrF2 and v = 10" Pa s. Lower horizontal velo- cities inhibit the rate of subsidence. The same effect is produced by an increased lithospheric viscosity (lower panel). With v = 5 X 10" Pa s a subsidence of around 0.4 mm yr-l corresponds to V = 10 cm yr-l, reduced to 0.1 mm yr-' when V = 1 cm yr-l. It has to be pointed out that an horizontal velocity of 10 cm yr-l is not unrealistic for the tectonic phases that shortened the crust in the Apenninic region (Castellarin & Vai 1986). A subsidence of around lmmyr-' with v = 10ZZPas is consistent with the deposition of lo00 m of sediments on a timescale of 1 Myr in the Adriatic foredeep, as derived from the stratigraphic record (Castellarin et al. 1985). Comparison of the two panels allows us to visualize the effect of an increased lithospheric viscosity on the deformation pattern. High viscosity curves are smoother due to the longer memory of the system. For these viscosities the rate of subsidence is substantially different from zero even in the region where low viscosity solutions are almost negligible, at around 200 km in our scale.
In Fig. 5 we show the effects on vertical motion due to discontinuities in time of tectonic phases. Solid curves correspond to the realistic case in which the load is expanding to the right in the lateral direction. Dashed curves are obtained when the model is loaded instantaneously at t = 0. We first discuss the realistic case. We assume that the overthrusting of the load terminates at 1 = 0 and we study the rate of vertical motion after 0.25 and We have chosen these intervals because they seem to correspond to the time elapsed since the end of the last tectonic transport in the contact region between the Apennines and Adriatic platform, as discussed more closely in Section 5. A time scale of around 0.25-0.5 Myr can thus be resolved from geological data. The parameter A corresponds to 0.25Myr and gives the interval of time elapsed since the end of overthursting. Curves corresponding to t = O are curves 1 of Fig. 4, here redrawn for comparison. We find that vertical motions strongly depend on the parameter A, especially for v = Pas. In the top panel there is almost an order of magnitude reduction between the situation in which horizontal motion is still ongoing (t = 0) and that in which the deformation is sampled after 0.25 Myr since the end of overthrusting. In the bottom panel with higher viscosity we get only a factor 2 reduction and even after 0.5Myr there remains a rate of subsidence of O.lmmyr-'. We thus find that the rheology of the lithosphere plays a crucial role in the rate of vertical motions in a situation in which the tectonics are discontinuous in time. Our results suggest that altimetric geodetic surveying can be used to discriminate between different hypotheses on the tectonic style of a region. Ongoing horizontal motions must be associated with higher vertical motions.
If we compare the dashed curves with the corresponding ones at t = A we can see the impact on vertical motions produced by a realistic time history for density anomalies.
First of all the deformation is symmetric around the density anomaly when the load is applied instantaneously (dashed curves), while the subsidence is higher at the front of the laterally expanding anomaly. We also get a huge reduction in the vertical motions when realistic loading history is considered. This demonstrates the importance of a realistic description in the problem.
The 0.5-1.0 mm yr-' rate of vertical deformation driven by density anomalies is larger than the value of to lo-' mm yr-I due to modifications in the horizontal stress field in a purely elastic lithosphere (Cloetingh et al. 1985).

ANALYSIS OF THE STRESS FIELD
We study the stress field in the crust forced by a surface load modelling a topography. We plot the contour of principal stress differences evaluated in the vertical crustal section represented in Fig. 1. Our attention is concentrated on the top layer where seismic events can occur due to its brittle nature. Ductile properties of the lower lithosphere make it easier to release stress through plastic processes instead of seismic faulting. Values are also given in bar (10 bar = 1 MPa). This unit is used to allow a direct comparison with the stress characteristic of seismic events. The formulation for a realistic time history is analogous to equation (7), except that the appropriate Green functions for the components of the stress field have to be considered. The results for a load that is applied instantaneously are plotted in Fig. 6 for comparison. The top panel represents the elastic contribution to the principal stress differences due to a surface load whose lateral extension is 2", as in the preceding figures of Section 3. We make use of the reference model (Table 1); the centre of the box corresponds to 300 km in the horizontal scale.
We decided to follow the evolution in time of the principal stress difference because it is twice the maximum shear, within the plane stress approximation. This quantity is thus a good indicator for seismicity forced by density anomalies.
The pattern of the elastic deformation is characterized by, two regions beneath the load edges in the deeper portion of the crust in which the principal stress differences are higher. After relaxation took place in the viscoelastic lithosphere, the density anomaly is isostatically compensated. This configuration is given in the bottom panel where the pattern is somehow modified, with the highest stress differences concentrated beneath the load. The largest values are of the order of lo3 bar (10' MPa). This equilibrium situation corresponds to the configuration that is obtained when the deep lithosphere beneath the elastic crust is assumed to be an inviscid fluid (Mareschal & Kuang 1986). The high level of stress that is reached in this configuration supported the claim that lateral density anomalies play a crucial role in intraplate seismicity (Caputo et al. 1984;Mareschal & Kuang 1986).
We now study the effects of a realistic time history with the front of the load migrating to the right to simulate the overthrusting process. In Fig. 7 (top panel) we provide the contour of the principal stress differences for a horizontal velocity V of 10 cm yr-' and a lithospheric viscosity v = Pa s. This panel is a picture of the steady state configuration that the system reaches when the front of the anomaly has migrated continuously for a sufficiently long time. We get very high stress differences, comparable with those corresponding to the completely compensated density anomalies given in the bottom panel of Fig. 6. The pattern is now asymmetric, with the highest values and gradients concentrated beneath the moving front of the load. This asymmetry, that arises when a realistic rheology for the lithosphere is considered, clearly suggests the limits of simpler models in a situation of ongoing tectonic activity. If we stop the lateral motion of the load and we leave the model unperturbed for a sufficiently long time, the stress differences increase in the left part of the panel that corresponds to the crust underlying the oldest portion of the load. This modification in the stress pattern continues until a symmetric configuration is reached. This final situation is drawn in the bottom panel of Fig. 7 corresponding to 1 Myr after the overthrusting process has been stopped. This figure is essentially identical to the bottom panel of Fig. 6 .
Except for the asymmetry in the stress pattern, there are no substantial deviations with respect to simpler approaches in the level of principal stress differences that can be reached in the crust. The new feature that is brought out by our modelling is the rate of accumulation of principal stress differences. This quantity is more significant than the value of the stress itself, not being affected by possible plastic processes in the crust and recurring earthquakes. These phenomena reduce the stress level on time scales of 10"yr. We are in fact interested in the possibility that an accumulation of stress differences of the order of 5-10 bar (0.5-1 MPa) occurs on time scales of 102-103 yr, comparable with the historical memory of human beings. On this time length, the description of the crust in terms of an elastic material is quite reasonable.
In Fig. 8 we plot the rate of accumulation of the principal stress differences, in units of bar/104 yr (0.1 MPa/104 yr), for a load that is expanding to the right. The solid triangles in this figure and following ones give the position of the load at the instant of time in which the rate of stress accumulation is evaluated. The triangle on the right denotes the position of the overthrusting front with respect to the crust. Thin contours correspond to a reduction of the principal stress difference, while thick ones stand for an increase. In this figure the lithospheric viscosity beneath the elastic top layer is fixed to 10'' Pa s and we study the effects of a reduction in crustal thickness (panel b) and horizontal velocity (panel c). The crustal region in which the rate is higher corresponds to the portion beneath the migrating front of the load at the crust-lithosphere boundary. In the top portion of the crust, and proximity of the crustlithosphere boundary beneath the centre of the anomaly, we have a reduction in the principal stress differences. This reduction shows how results derived within the framework of inviscid models for the lithosphere can be misleading.
These models predict that the largest stress difference concentration occurs beneath the centre of the anomaly. We find on the contrary that this region experiences a slight reduction in this quantity when a realistic rheology is considered. This portion of the crust must be seismically inactive with respect to the deep crust beneath the front. We also have to note a feature that is common to the three panels. Close to the crustal portions in which the rate of stress accumulation is higher, we find another region that experiences a consistent reduction in the principal stress differences. If the largest rates of stress accumulation are linked with seismic activity, we thus get that seismicity in the deep crust appears tightly concentrated in a well defined region beneath the active front of the topography.
The rates of stress accumulation that we have obtained indicate that in order to gain 5 bar (0.5 MPa) in the principal stress differences, we have to wait around 700yr for the most favourable case of panel (b), if the only active tectonic mechanism is the lateral displacement of density anomalies. We considered this value for the stress difference because this is consistent with a reasonable estimate of the apparent average stress and stress drop of earthquakes in the Mediterranean region (North 1977) and Appalachians (Mareschal & Kuang 1986) where our model can be applied.
The tectonic mechanism that we considered predicts a recurrence time for earthquakes of ld-103yr. This is of course a rough estimate due to the simplicity of the model that does not account for lateral variations in rheology.
Comparison of panels (a) and (b) shows that rates of principal stress difference accumulation are inversely related to the thickness of the crust. Panel (c) emphasizes the linear relation between the horizontal velocity and rate of stress accumulation within this range of viscosity values.
We may ask whether the lateral dimension of the varying density anomaly plays any role in the estimated rates of stress accumulation. To point this out, we consider in Fig. 9 a surface load of 0.5" of lateral extension that is expanding to the right at a horizontal velocity V = lOcm yr-'. If we compare panels (a) of Figs 8 and 9 that correspond to the same parameters, except for the extension of the load, we see that the largest values of the rate of stress accumulation remains almost unaltered. What is modified are the gradients, the smaller load predicting smoother variations in the rates. The region in which we get the largest rates is, as before, beneath the active front. Regions of strong reductions in the rates almost disappear for the considered load.
We found that if the viscosity is sufficiently low, say lower than lo2' to 5 x 10" Pas, the rates of stress accumulation are substantially independent of the rheology of the lithosphere for fixed horizontal velocities. If the viscosity is increased beyond this threshold, around Pa s, the rate of accumulation suffers a strong reduction, as deduced from panel (b). A thinner crust with this load has the same effects already derived for the 2" wide anomaly, except that the rate is slightly lower.
In our analysis we assumed that the overthrusting process is still active. We want to consider the possibility that this process came to an end, and study the modifications in the stress pattern. This picture is depicted in Fig. 10. In panel (a) we provide the rate of principal stress difference accumulation in the same units as the preceding figures, at the instant in which the horizontal motion is stopped. Panel (b) is a picture taken at lo5 yr after the end of overthrusting, with the same parameters as panel (a). If we compare these contours with Fig. 8, we now get that stress differences are essentially increasing to the left, beneath the stable edge of the load in the backarc region. Rates are generally lower in comparison with the configuration of active tectonics, but   Figure 9. Contours of the rate of principal stress difference accumulation corresponding to a 0.5" wide topography. Horizontal velocity is fixed at 10 cm yr-'; viscosity is I d ' Pas except in panel (b) where it has been increased to Pas. Crustal thickness is varied to 15 km in the bottom panel.
the highest values in panel (b) are quite close to those of panel (a), Fig. 8. Another peculiarity is that for longer times after the end of overthrusting, the region where stress differences are large becomes wider, as deduced from the comparison of panels (a) and (b). Rates of stress accumulation decrease until the final equilibrium configuration of Fig. 7, bottom panel, is reached. Panel (c) represents a limit situation in which the lithospheric viscosity has been lowered to Id' P a s and horizontal motion just turned off.
In this case, around 200 yr are sufficient to accumulate 5 bar (0.5 MPa).
From Fig. 10 we thus find that seismicity in the elastic crust tends to be concentrated in the backarc region close to the ductile region if horizontal tectonics comes to an end.
It is easy to make an estimate of the efficiency of this tectonic mechanism in triggering earthquakes if results are compared with the accumulation of shear stress at plate margins along transcurrent faults. In the shear zone we have for the shear stress rate t = 2 p K , where K = 4 X 10-'5s-1 is appropriate for the strain rate in the Mediterranean region (North 1974). With p = 3 X 10"N m-' we thus find that around 750 bar (75 MPa) could be accumulated in lo4 yr. This is an order of magnitude larger than the rates of stress accumulation obtained in this paper, except in Fig. 10 where 300 bar (30 MPa) are built in the same span of time.
When a viscoelastic Maxwell rheology is assumed for the lithosphere, we find that density anomalies associated with collision processes are less efficient in inducing shear stresses in the brittle layer with respect to transcurrent motions. Only for viscosity values like 10" Pas, that is certainly a lower bound for the lithosphere, does shear stress accumulation for the two mechanisms become comparable.
When comparison is made with seismicity in transcurrent shear zones (Dragoni, Bonafede & Boschi 1986), we thus expect a lower level of seismicity and longer recurrence times. This is consistent with the moderate seismicity in the Appalachians and northern and central Apennines (Mareschal & Kuang, 1985;Mulargia, Gasperini & Tinti 1987), where overthrust is the prominent tectonic mechanism (Quinlan & Beaumont 1984;Vai 1988). For the Apennines a more detailed description is given in the next secton.

A SYNTHETIC DESCRIPTION OF T H E TECTONIC SETTING OF THE ITALIAN REGION (CENTRAL A N D NORTHERN APENNINES)
The geological setting of the Italian peninsula is characterized, in its central part, by the different phases of activation of the Apenninic front. In Fig. ll(a) we provide a tentative chronologic classification of the active fronts of the chain. This figure is redrawn from Vai (1988). Starting from the inferred position of the front 15Mya (Burdigalian front), we can recognize the migration of tectonic activity, until the last pulses at the end of late Pliocene and mid-late Pleistocene. These tectonics can be interpreted as the overthrust of the Apennines on the Adriatic foreland (Castellarin & Vai 1986). Downwarping of the plate beneath the loads results in the development of the Adriatic foredeep ( Fig. 1 la). The horizontal velocities of overthrusting, that are needed in this paper, can be estimated from the migration of the reactivated fronts (Vai 1988).
The reconstruction of the timing of tectonic phases has been allowed by the structural interpretation of seismic surveying in the contact region between the chain and the Adriatic platform. Large, syntectonic, sedimentary wedges, characterizing the main pulses of thrusting, can be recognized in the seismic profiles (Pieri 1983; Bally et af. 1986). Tectonic transport was still active in the middleupper Pliocene, as documented by sedimentary wedges of this age. The flat, lenticular shape of late Quaternary strata seems to indicate the end of tectonic transport during the various extensional phases of the basin. A coherence and late Pleistocene, that is around 0.5 Mya (Castellarin & Vai synchronism is found between the main deformation events 1986).
in the basin and surrounding Apennines (Sartori 1988). Crustal deformation in the Apenninic chain can be Density anomalies induced by tectonic activity, both in the compared with the evolution of the Tyrrhenian basin.
chain and Tyrrhenian basin, are needed for the study of Results from the ODP Leg 107 and previous geological and vertical motions. In Sections 3 and 4 we limited our geophysical data, allowed reconstruction of the timing of the discussion to the process of overthrusting. In the next section we specialize our model to take into account the density anomalies forced by the main tectonic processes characterizing the Tyrrhenian basin and surrounding chains.
To obtain the crustal structure, we averaged both in the horizontal and vertical direction the density model obtained by Baldi ef al. (1982) from Bouguer gravity data. Geological and seismotectonic information has also been used t o better constrain the results of gravity inversion. In  Table 2.
In the crust we have two major features that are associated with the opening of the Tyrrhenian sea (positive anomaly to the left) and crustal thickening (negative anomaly to the right). Crustal thickening is probably due to overthrusting of the Apennines on the Adriatic plate. Crustal thickening and thinning are described in our approach by surface density anomalies buried at the crust-lithosphere interface.
In the near future geodetic and satellite data on horizontal crustal deformation will become available to better constrain the horizontal velocities for the model. The optic WegenedMedlas program (Wilson, Reinhart & Pearlman 1985) will provide important information for the Mediterranean region.

A DETAILED CRUSTAL MODEL FOR THE APENNINES: CONSEQUENCES FOR VERTICAL MOTIONS
In this section we study the vertical motions forced by the laterally varying crustal structure in the central part of the Italian peninsula, as shown in Fig. l l ( b ) . When a vertical plane normal t o the Apennines and deep crustal anomalies is considered, we basically deal with a 2-D deformation problem, as discussed in the present paper. 1 mm yr-' is the only forcing mechanism for the dashed curve of panel (a). This rate is well constrained by geological data (Castellarin & Vai 1986). We have to emphasize that sedimentation appears to be an important source, characterizing the deformation pattern of the whole set of curves along the Adriatic coast. Dotted curves account for the effects of sedimentaion and underthrusting of the Adriatic plate beneath the crust. In this process light material is injected into the lithosphere and this causes uplift of the crust, as portrayed by positive vertical velocities. The topographic maximum of the Apennines, at around 320 km in our scale, is uplifting with respect to the Adriatic coast. The relative rate of uplift and subsidence between the topograpy and Adriatic foredeep, of around 0.5 mm yr-', is only slightly sensitive to the velocity of underthrusting (panels a and c). A small reduction is expected if vertical motions are recorded after 0.25 Myr from the end of underthrusting (panel b).
Dash-dotted curves correspond to the situation in which crustal thinning due to the formation of the Tyrrhenian sea is also taken into account. The positive density anomaly is responsible for the subsidence in the western portion of the peninsula with respect to the topographic maximum. This For a better understanding of the interplay between the different contributors, we start by adding separately the density anomalies. In Fig. 12 solid triangles denote the positions along the transect of the Adriatic coast (right) and Tyrrhenian coast (left). Lithospheric viscosity is kept fixed to 10" Pa s. Horizontal velocity is reduced from 10 cm yr-' in panels (a) and (b) to 5 cm yr-' in panel (c). Velocities V have been intended as relative velocities when density anomalies grow in the opposite directions, as those corresponding to the opening of the Tyrrhenian sea and crustal thickening. The parameters N and Ar in equation (8) have been chosen in such a way as to constrain the velocity of lateral expansion and the observed position of the density anomaly at the time in which deformation is recorded.
Panel (b) corresponds to an interval of 0.25 Myr elapsed since the end of horizontal motions while in (a) and (b) tectonic transport is active. We have chosen this value for the period of quiescence in tectonic activity in agreement with geological data (Castellarin & Vai 1986). A termination of tectonic transport is, on the other hand, controversial (Patacca & Scandone 1988) and so we were obliged to consider the different possibilities in the model.  1,22.5 km for 2 and 30 km for curve 3. In the bottom panel we fix D = 25 km while the lithospheric viscosity is increased from ld2 Pas for curve 1 to ld3 Pa s for curve 3. Index 2 corresponds to 5 X lo2* Pas. As before, solid triangles give the position of the Tyrrhenian and Adriatic coasts. subsidence is strongly dependent on the rate of lateral expansion of the Tyrrhenian basin. The rates of 5-10 cm yr-can be estimated from the amount of crustal shortening in the Tyrrhenian chains that seems to be synchronous with the phases of basin formation. (Sartori 1988;Vai 1988). Strictly speaking, the opening of the basin is not a continuous process as we have assumed in the model. On the other hand, a more precise reconstruction for the history of basin formation will not substantially alter our results. Solid curves correspond to the most complete tectonic model in which overthrusting of the Apennines on the Adriatic plate is accounted for. Overthrusting of the chain enhances the subsidence of the whole peninsula. The deformation pattern is characterized in the whole set of panels by a large subsidence of the order of 0.4 mm yr-' in the Adriatic part due to the combined effects of the overthrusting Apennines, the concomitant development of a large positive anomaly connected with the Tyrrhenian sea and sedimentation in the Adriatic foredeep. In the western part, subsidence is reduced to around 0.1 mm yr-'.
We show some numerical experiments in which the thickness D of the elastic top layer and lithospheric viscosity are varied to allow for deviations from our reference model. In Fig. 13 we assumed V = 10 cm yr-'. In the top panel, we vary the crustal thickness, keeping the viscosity fixed at 5 x 10Z2Pas. Curve 3 is the solid one of Fig. 12 that corresponds to the complete tectonic model. We find, as expected, that the thickness of the elastic top layer plays a crucial role in controlling the rate of vertical deformation. In particular we find that a thin elastic crust, which is probably more realistic at least in the Tyrrhenian portion of the Italian peninsula, predicts a larger subsidence in the western part with respect to the topographic maximum. We thus find that topography is uplifting with respect to the eastern and western coasts. This occurs because the crust supporting the topography is subsiding at a lower rate. If the crustal thickness is reduced to 15 km, subsidence increases to around 0.25 and 0.5 mm yr-' (respectively) to the left and to the right of the topographic maximum at 320 km.
The only available data on vertical deformation along the Tyrrhenian coast have been derived from the tide gauge recordings (Caputo & Pieri 1976). The increase in sea-level obtained by these authors is consistent with the subsidence of the coast. An increase in global sea water of 0.46 mm yr-' is attributable to ongoing glaciers melting (Meier 1984). If we take into account present day glacial activity, the 1.4 mm yr-' sea-level increase derived by Caputo & Pieri (1976) is reduced to around 0.9 mmyr-I plus or minus a 50 per cent deviation. We thus find that at least a fraction of this value can be due to tectonic activity. The other contribution can be explained by thermal expansion of the oceans (Meier 1984).
In the bottom panel we analyse the impact of viscosity variations. We find that lower viscosities enhance the effects of sedimentation as a consequence of the shorter memory of the system that reduces the contributions of the other sources. A viscosity increase has the effect of enhancing the subsidence in the Tyrrhenian coast.
Variations in the parameter A, denoting the interval elapsed since the end of horizontal motions, are considered in Fig. 14 for v = loz3 Pa s and V = 10 cm yr-'. Subsidence in the western part of the peninsula is very sensitive to this parameter. These results demonstrate the impact of discontinuities in tectonic phases on vertical deformation.

CONCLUSIONS
Sphericity plays an important role in long wavelength tectonic processes when a viscoelastic rheology is assumed for the lithosphere. This is valid for time-dependent density anomalies and must be true for other global geophysical phenomena, like stress diffusion due to large earthquakes in subduction zones. The model given in the paper can be easily adapted to deal with this class of problems, due to the simplicity of the analytical approach.
Vertical ground deformation in collision zones is found to be mainly controlled by lateral displacement of crustal and lithospheric density anomalies. If overthrusting is considered, the rate of subsidence of the plate beneath the load has a very peculiar pattern, with the largest deformation localized beneath the active front of the topography. The rate of subsidence is around 1 mm yr-' for a chain with characteristics appropriate for the Apennines. This result is consistent with the amount of sediments that are deposed in the Adriatic foredeep (Castellarin et al. 1985). The amplitude of subsidence decreases in the stable portion of the chain. In this region, the amount of deformation is highly sensitive t o lithospheric viscosity. Other parameters that have an impact on deformation pattern are the horizontal velocity of overthrusting and the timing in the discontinuities of tectonic phases. On the time scale of 1 6 yr, vertical motions can decrease drastically if tectonic transport terminates.
From analyses of the principal stress differences and their variations in time, it seems that the role of density anomalies on seismicity has t o be reconsidered and conclusions drawn from unrealistic models have to be somehow modified. The rate of accumulation of shear stresses is relatively low when it is compared with that induced by more powerful mechanisms, such as transcurrent motions along plate margins. W e found that the rate of shear stress accumulation is around an order of magnitude smaller when comparison is made with active shear zones. This implies a moderate seismicity characterized by relatively long return times. We estimated the recursive times of earthquakes triggered by the considered source mechanism in multiples of a hundred years. This is consistent with the clusters of the largest earthquakes in the Northern Apennines (Tinti, Lenzi & Mulargia 1985).
Also for the stress field, the rate of accumulation is very sensitive to the velocity of overthrusting and viscosity in the deep layer. We revealed that rates of stress accumulation are larger in the deeper portion of the crust beneath the active front, unlike the results from simpler models.
Discontinuities in time of tectonic phases a r e also important in characterizing the stress pattern. When our model is used to study the vertical motions in the Italian peninsula, it suggests that the deformation depends on the tectonic style that is assumed. Global subsidence of the peninsula is characterized by large subsidence in the Adriatic coast. The amount of horizontal migration of density structures o r discontinuities of tectonic phases plays an important role in the pattern of ground motions. This suggests that, in the near future, accurate geodetic surveying can become a useful tool to better constrain tectonic models.