Influence of glacial isostatic adjustment upon current sea level variations in the Mediterranean

Tectonophysics 474 (2009) 56–68

Contents lists available at ScienceDirect

Tectonophysics j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / t e c t o

Influence of glacial isostatic adjustment upon current sea level variations in the Mediterranean P. Stocchi a, G. Spada b,⁎ a b

DEOS, Faculty of Aerospace Engineering, Delft University of Technology, The Netherlands Istituto di Fisica, Università di Urbino “Carlo Bo”, Urbino, Italy

a r t i c l e

i n f o

Article history: Received 15 January 2008 Received in revised form 18 August 2008 Accepted 2 January 2009 Available online 15 January 2009 Keywords: glacial isostatic adjustment sea level variations Mediterranean

a b s t r a c t Present-day sea level variations in the Mediterranean depend on various factors, including recent climatic forcing, tectonic activity, anthropogenic effects, and glacio-isostatic adjustment. The latter is governed by mantle rheology and the spatio–temporal distribution of the late-Pleistocene ice sheets and it is expected to produce a long-wavelength pattern of sea level variations across the Mediterranean, mostly determined by the response of the solid earth and of the geoid to loading effects of melt water since the end of deglaciation. Modeling glacio-isostatic effects in this region is necessary for a correct interpretation of tide gauge and GPS time-series, and thereby to constrain both the present-day climate-related sea level rise and regional or local geological, tectonic and human-driven displacements. By an exhaustive exploration of the parameter space of mantle rheology and ice sheet chronologies, in this work we outline upper and lower bounds on the current rate of sea level variation associated with glacial isostatic adjustment in the Mediterranean. This may contribute to a full assessment of coastal vulnerability by sea level rise on a regional and local scale. © 2009 Elsevier B.V. All rights reserved.

1. Introduction The global mean sea level rise is a direct response to an ongoing planetary climate change (Church et al., 2001). Various investigations motivated by a growing concern about an anthropogenic forcing on global warming show a widespread spatio–temporal variability of sea level change, which reflects the simultaneous action of a number of complex geophysical mechanisms (Douglas et al., 2001; Cazenave and Nerem, 2004). The most important long-term mechanisms, spanning thousand to hundred thousand years, are glacio- and hydro-isostatic readjustment and vertical deformations associated with active tectonics and marine sedimentation. Other natural mechanisms acting on a century scale are changes in winds (Fukumori et al., 2007), sediment compaction (Meckel et al., 2006), and sediment loading subsidence (Ivins et al., 2007). While the effects of glacial isostasy exhibits a smooth, long-wavelength character, the latter may be spatially heterogeneous and is generally characterized by short spatial scales. Furthermore, at a local scale, various human activities such as water, oil, and gas extraction and land reclamation may contribute to significant rates of subsidence that result in relative sea level variations. Decomposing the currently observed and spatially variable rates of sea level variation into individual components (see e. g., Carminati and Di Donato, 1999) thus demands information about local and regional crustal deformations and geoid height variation and ⁎ Corresponding author. Istituto di Fisica, Università di Urbino “Carlo Bo”, Via Santa Chiara n. 27, 61029 Urbino (PU), Italy. Tel./fax: +39 0722 303389 99. E-mail address: [email protected] (G. Spada). 0040-1951/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.tecto.2009.01.003

knowledge of the response of the Earth to the melting of the late Pleistocene ice sheets. Due to the presence of geological, geomorphological and archaeological indicators (Flemming and Webb, 1986; Lambeck et al., 2004a,b; Lambeck and Purcell, 2005; Pirazzoli, 2005; Antonioli et al., 2007) and the availability of instrumental observations (Zerbini et al., 1996; Serpelloni et al., 2006, 2007), the Mediterranean coasts offer a unique opportunity to investigate the various mechanisms that contribute to sea level variations on different time-scales. The whole Mediterranean basin is affected by a significant tectonic activity contributing to a widespread coastal instability. Both episodic and continuous vertical movements displaced, in the course of time, different coastal features above or below present-day mean sea level, preventing depositional or erosional features from the destructive wave action (Pirazzoli, 1991). Natural high rates of subsidence due to sediment compaction have been accelerated during the last decades by human activity causing, as in the case of NE Adriatic region, a significant beach retreat accompanied by a local increase of bathymetric slope and a reduction of the coastal sand reservoir (GNRAC, 2006). While it is now assessed that at least 70% of the world's beaches are in permanent retreat in response to global sea level rise (Day, 2004), recent studies of present-day coastal equilibrium in Italy indicate that 42% of the Italian beaches are eroding (GNRAC, 2006), a figure which is affected by considerable uncertainties because of the intrinsic difficulties in the assessment of the relationship between sea level rise and erosion (e. g., Day, 2004). Estimating present-day climate-related sea level rise (Fukumori et al., 2007) and constraining vertical deformations from geologically observed rates in the Mediterranean demands an accurate model of

P. Stocchi, G. Spada / Tectonophysics 474 (2009) 56–68

glacial isostatic adjustment (GIA) which, in turn, depends on assumptions about both the melting history of continental ice sheets since the Last Glacial Maximum (LGM), and the viscoelastic response of the solid Earth. Knowledge of the role of melting of late Pleistocene ice sheets on the relative sea level variations in the Mediterranean mainly comes from regional analyses aimed at the interpretation of various geomorphological and archaeological observations (Lambeck et al., 2004a,b; Lambeck and Purcell, 2005; Antonioli et al., in press). Due to evidence of significant sea level variations in the past, a number of studies have focused on the northern Mediterranean coasts (see e. g., Lambeck and Bard, 2000). These are, potentially, the most affected by the process of isostatic adjustment due to the proximity to the former Alpine and Fennoscandian ice sheets (Stocchi et al., 2005), since ice unloading shapes the overall pattern of land subsidence within a distance of ∼40° to ∼60° from the ice centers (Lambeck and Johnston, 1995; Sella et al., 2007). However, GIA-related deformation of the whole Mediterranean basin is mainly driven by water-loading, which contributes to a significant and widespread subsidence whose extension and strength, in turn, directly depend on ice sheet chronology and Earth viscosity. In the bulk of the basin, the average rate of sea level change has been close 1 mm yr− 1 during the last ∼ 6 kyr (Lambeck and Purcell, 2005; Stocchi and Spada, 2007), consistent with the analysis performed in this study. It is known that the viscosity contrast between upper and lower mantle and lithospheric thickness control significantly GIA in the Mediterranean (Spada and Stocchi, 2007). However, volume-averaged viscosity values inverted from global modeling may not be suitable on a regional scale, also considering the heterogeneity due to the complex tectonic setting of the region (e. g., Piromallo and Morelli, 2003). Since the parameters space of rheological layering and ice sheets chronology has not been fully sampled by sufficient data, our aim here is to estimate current rates of sea level change for a large number of forward models. The results will be expressed in terms of bounds instead of point estimates of relevant geophysical quantities. This contributes to a better understanding of the trade off between various factors presently contributing to sea level change in the Mediterranean. This goal will be accomplished by a straightforward solution of “Sea Level Equation” (hence after SLE) in the form originally proposed by Farrell and Clark (1976). Information about mantle viscosity and ice sheet chronology mainly come from inversion of different geological, geodetic and geophysical observables, often characterized by a high degree of nonuniqueness (Cianetti et al., 2002 presents a classical case study). For this reason, after a short description of the methods, in Section 3 we first assume ICE-5G chronology with the VM2 viscosity profile (Peltier, 2004) to predict the present-day rates of sea level change and geoid height variations across the Mediterranean. In a subsequent step, we analyze the effects of varying the Earth rheological parameters and ice sheet chronology to infer upper and lower bounds for rates of GIAinduced rates of sea level variation. The ice sheet time-histories considered here are based on the literature and differ by assumptions concerning the melting of remote ice sheets. In particular, we will consider ICE1 (Peltier and Andrews, 1976), ICE-3G (Tushingham and Peltier, 1991), ICE-5G (Peltier, 2004), and an ice model that has been recently calibrated using Holocene relative sea level observations from the Mediterranean by Stocchi and Spada (2007), here referred to as MDCR, which includes a catastrophic rise event of Antarctic origin (CRE3) 7 kyr BP (Blanchon and Shaw, 1995). Finally, we will compare model predictions with observations at various sites in which instrumental tide gauge observations are available. 2. Methods 2.1. Theory and numerical approach GIA is described by the SLE, the linear integral equation first proposed by Farrell and Clark (1976). Recently Spada and Stocchi

57

(2007) have implemented the SLE in a public-domain code (SELEN) based on three main approximations: (i) the Earth is radially stratified, incompressible and characterized by a linear Maxwell viscoelastic rheology, (ii) the shorelines do not migrate horizontally as the sea level changes, and (iii) the effects of Earth rotation on GIAinduced sea level variations can be neglected. In addition to providing a correct description of sea level variations, this zeroth order model also allows for a description of other relevant geophysical quantities related with GIA, such as vertical crustal deformations and geoid height variations. Sea level change S is determined by S = N − U;

ð1Þ

where N and U are the geoid height variation and vertical displacement of the solid Earth, respectively (Farrell and Clark, 1976) (to simplify notation, in the following we drop the dependence upon colatitude, longitude, and time of variables). Imposing the constraint of mass conservation gives     Φ Φ E S= −U −S − −U ; ð2Þ γ γ where Φ is incremental gravitational potential, γ is surface gravity, the overline indicates an average over the oceans, and E

S = −

mi ρw A o

ð3Þ

is the “eustatic term” of the SLE, where mi is the ice mass variation, ρw is water density, and Ao is the area of the surface of the oceans. In Eq. (2), both Φ and U depend on the surface load exerted by the ice sheets on land and by melt water load exerted on ocean floor. While the former can be reasonably assumed to be known from geological and glaciological evidence, the latter depends on S itself. This makes the SLE an integral equation S=

ρi ρ ρ ρ E G  I + w Gs o S + S − i Gs i I − w Gs o S; γ s i γ γ γ

ð4Þ

where I is ice sheet thickness variation, ρi is ice density, i and o denote spatio–temporal convolutions over the ice and ocean covered regions, and the last two ocean-averaged terms ensure mass conservation. The sea level Green's function Gs accounts for mantle viscoelasticity through the load-deformation coefficients for vertical displacement (h) and incremental potential (k) (Farrell and Clark, 1976; Spada and Stocchi, 2006, 2007). Following the original idea of Mitrovica and Peltier (1991) and Mitrovica et al. (1994), SELEN implements the “pseudo-spectral” method for the solution of the SLE (Spada and Stocchi, 2007). This method is based on a discretization of all the variables involved and, through the projection of the original unknown S onto the “ocean function”, it allows for a straightforward spectral approach. In our study spatial discretization takes advantage of the icosahedral pixelization of Tegmark (1996), which provides a set of equal-weight integration points that allow for a nearly optimal quadrature on the sphere. The temporal discretization is performed assuming that all variables vary stepwise (Δt = 1 kyr) in the time interval 0 ≤t ≤tp, where t = 0 denotes the beginning of melting and t =tp is present time. At the core of the algorithm is a recursive procedure in which as first guess it is assumed S =SE (see Eq. (3)). This test solution is substituted on the right-hand side of Eq. (4), and the process is iterated. The procedure quickly converges to a stable solution. In all the ensuing computations, a maximum degree lmax = 128 is employed, which represents an optimal compromise between accuracy and efficiency in computation. The resolution of the integration grid is R = 48, which corresponds to 40R(R − 1) + 12 = 90,252 pixels on the surface of the sphere (Tegmark, 1996). With these settings, our RSL predictions match those by

58

P. Stocchi, G. Spada / Tectonophysics 474 (2009) 56–68

Fig. 1. Equivalent sea level (ESL) for the ice sheets chronologies employed in this study. The total ESL variation between the LGM and the end of melting is shown next to each curve. The catastrophic rise event (CRE3) of Antarctic origin that is included in MDCR is marked in frame (b).

Lambeck et al., (2004a,b), provided that the same viscosity profile is employed (Antonioli et al., in press). Once S is obtained, a suite of geophysical quantities such as present-day rates of sea level change (Ṡ), rate of vertical deformation of the solid surface (U̇), and rate of geoid height variation (Ṅ) can be retrieved by discretizing the derivative  : dR  θ; λ; tp ; R= dt

ð5Þ

where θ is colatitude, λ is longitude, and R denotes any of the scalar fields in Eq. (1). Predicted rates can be directly compared with available observations, such as present-day trends of mean sea level variations at coastal tide gauges or vertical velocities recorded by continuous GPS stations on land (Stocchi et al., 2009b). 2.2. Ice sheets chronologies In the following, we will employ various ice sheet chronologies that differ for the history of melting and equivalent sea level (ESL), defined as ESL =

ρi V ; ρw A o

ð6Þ

where V is volume of melt water. The ESL of the ice sheets considered in this study is shown in Fig. 1 as a function of time. The recent ICE-5G chronology of Peltier (2004) constitutes our reference ice model (ICE-5G is available from the Special Bureau of Loading web page http://www.sbl.statkart.no/projects/pgs/). Storing an ESL of ∼127 m at the LGM (21 kyr BP), and consisting of 22,720 nonoverlapping quadrilateral ice elements, ICE-5G improves the fit between Holocene relative sea level observations and predictions with respect to previous low-resolution models, such as ICE-3G of Tushingham and Peltier (1991). In our implementation of ICE-5G, we assume isostatic equilibrium prior to the LGM, and deglaciation in steps of 1 kyr. With an Antarctic ESL larger by ∼50% relative to ICE-5G and a reduced ice thickness of North America and northern European ice sheets, ICE-3G also differs from ICE-5G for its spatial discretization, consisting of 808 axisymmetrical disk-shaped ice elements. ICE-3G assumes LGM to occur 18 kyr BP and stores a total ESL of ∼114 m. As discussed by Stocchi and

Spada (2007), ICE-3G does not fit satisfactorily the Holocene RSL curves in the Mediterranean, mainly because of the relatively large ESL variation of Antarctica which causes a general late-Holocene highstand that is not supported by the observations. In addition to ICE-5G and ICE-3G, here we also consider the two hybrid ice models ICE1+A3 and MDCR, which have been recently derived from the original ICE1 model of Peltier and Andrews (1976) to test the sensitivity of late to mid-Holocene relative sea level curves in the Mediterranean to the melting of remote ice aggregates (Stocchi et al., 2009a). Model ICE1 consists of 153 quadrilateral ice elements specified on a coarse 5° × 5° grid and describes the melting of northern hemisphere ice sheets since 18 kyr BP (ESL ∼ 78 m), assuming a constant ice thickness over Antarctica during the whole Holocene. ICE1+A3 combines the Antarctic component of ICE-3G with ICE1 (the latter converted into disk elements for computational convenience). Since A3 implies a gradual ESL variation of ∼28 m between 9 and 5 kyr BP (see Fig. 1), it produces a significant alteration of the melting history of ICE1, which ends 8 kyr BP. MDCR differs from ICE1+A3 for a reduced Antarctic ESL and assumes a sudden release of melt-water from Antarctica 7 kyr BP. This epoch approximately corresponds to the catastrophic rise event referred to as CRE-3 by Blanchon and Shaw (1995), whose magnitude in terms of ESL is estimated as 6.5 ± 2.5 m. From the analysis of Stocchi et al. (2009a), a CRE event of Antarctic origin improves the fit of model predictions with Holocene RSL observations in the Mediterranean. In MDCR, the Antarctic ice sheet is modeled by a single disk with radius of 20°, which is appropriate due to the large distance from the Mediterranean basin.

Table 1 Earth model parameters employed in this study. Layer

Range (km)

Density (kg m− 3)

Rigidity (GPa)

RVM2

RVM1 (×1021 Pa·s )

RVM3

LT SU TZ LM CO

6281–6371 5951–6281 5701–5951 3480–5701 0–3480

4120 4120 4220 4508 10,925

73 95 110 200 0

∞ 0.4 0.4 4.0 0

∞ 1 1 2 0

∞ 0.4 0.4 40 0

LT, SU, TZ, LM, and CO stand for (elastic) lithosphere, shallow upper mantle, transition zone, lower mantle, and (fluid) core, respectively. RVM columns show viscosity values.

P. Stocchi, G. Spada / Tectonophysics 474 (2009) 56–68

2.3. Earth model In all the computations that follow, we consider a five-layers Earth model characterized by an elastic lithosphere (LT), a three-layer Maxwell viscoelastic mantle including a shallow upper mantle (SM), a transition zone (TZ), a lower mantle (LM), and a homogeneous inviscid core. While density, rigidity and thickness of layers are kept fixed to values given in Table 1, lithospheric thickness HLT is varied, which implies a variable layer thickness UM. Upper and lower mantle viscosities ηUM and ηLM are varied as well, over a broad range of a priori possible values. Since we will assign the same viscosity to the shallow upper mantle and transition zone, hence after with ηUM we will denote mantle viscosity to a depth of 670 km.

59

The rheological profile of our reference Earth model is characterized by ηUM = 0.4, ηLM = 40 × 1021 Pa·s, and HLT = 90 km (see Table 1). Such profile, which here is referred to as RVM2, represents a rough discretization of profile VM2 of Peltier (2004), and differs significantly from previous VM1 profile of Tushingham and Peltier (1991) (here denoted by RVM1). The latter is in fact characterized by ηUM = 1 and ηLM = 2 × 1021 Pa·s that define the rheological profile implied in ice chronology ICE-3G. The need of a reduction of ηUM is also supported by a number of independent GIA investigations based on Holocene relative sea level or present-day geodetic data (Mitrovica, 1996; Spada, 2001). Since estimates of ηLM vary considerably in the range between 5 and 50 × 1021 Pa·s (Nakada and Lambeck, 1989; Mitrovica and Forte, 1997; Lambeck et al., 1998), in some of our computations below we

Fig. 2. (a) Predicted rates of present-day sea level change Ṡ and vertical deformation of the solid surface U̇ for our reference model ICE-5G(RVM2) (see Table 1). White circles show four relevant PSMSL tide gauge stations considered in the sensitivity analysis of Section 3.2: Cagliari (CA), Trieste (TR), Malaga (MA), and Hadera (HA). Black triangles show other four PSMSL sites considered in the comparison between observations and predictions (see Section 3.4): Tarifa (TA), Algeciras (ALG), Alicante (ALI) and Genova (GE). Inset shows the Po plain region. (b) Present-day rates of geoid height change Ṅ = Ṡ + U̇ (see Eq. (1)), for the same ice and Earth model considered in (a). Ṅ represents the absolute rate of sea level change measured with respect to the Earth's center of mass, which would be observed by an artificial satellite (Peltier, 2001).

60

P. Stocchi, G. Spada / Tectonophysics 474 (2009) 56–68

will also consider a rheological model (RVM3) characterized by a viscosity increase of two orders of magnitude across the lower and upper mantle. 3. Results 3.1. Reference model For our reference model ICE-5G(RVM2), Fig. 2 shows predicted values of Ṡ and U̇ (frame a), and Ṅ (b). These three quantities are not independent of one each other, being connected by Eq. (1). Despite the fact that both Ṡ and U̇ are defined both on land and on sea, in (a) these fields are shown separately to facilitate qualitative comparisons with observed rates, measured by tide gauges at sea and by GPS on land (Serpelloni et al., 2006; Bennett and Hreinsdòttir, 2007). In the bottom frame, Ṅ is shown both on land and on sea. All predictions have a characteristic amplitude of fractions of millimeter per year, roughly an order of magnitude smaller than those expected in the bulk of formerly glaciated areas (e. g., Di Donato et al., 2000; Spada and Stocchi, 2007). Fig. 2 shows a widespread subsidence of both the solid surface of the Earth and of the geoid across the Mediterranean region, and clearly indicate increasing rates moving from inland to the bulk of the basin. It is remarkable that |Ṅ| is not small compared to |U̇|, as it occurs in the depressed, formerly glaciated regions and in the surrounding uplifted forebulges, where present-day sea level changes are mainly driven by the vertical deformations of the solid surface (Di Donato et al., 2000). However, while Ṅ shows little variability across the study region, which reflects a large content of long-wavelength harmonics, U̇

displays a significant spatial variability and sharper spatial gradients. The shape of contour lines in Fig. 2 is similar to that obtained by Di Donato et al. (2000) using on model ICE-3G. Differences in numerical values of Ṡ are to be attributed, in addition to ice sheet chronology, to differences in mantle layering and spatial resolution. The pattern observed in Fig. 2 can be explained by the loading effect of melt water added into the oceans until 4 kyr BP, the epoch that marks the end of deglaciation according to model ICE-5G. The influx of water has caused a general subsidence accompanied by a collapse of the geoid (see frame b), which still results in a sea level rise reaching the largest amplitude close to the center of the basin. In particular, the largest rates of sea level rise (∼ 0.5 mm yr− 1) are expected in the bulk of the Tyrrhenian Sea (Sardinia), and southeast of Italy, between Sicily and Greece (Ionian Sea). Predicted Ṡ values decrease towards the coasts and vanishes between south Spain and Algeria (see thick black curve) and along the indented portions of the north African coasts between 10° and 35°, such as southeast Tunisia, south Lybia, and south Israel. Moving westward Ṡ changes its sign and a sea level fall is predicted between south Spain and northwest Africa. For the Po plain (northern Italy, see inset), predicted subsidence rates are close to 0.5 mm yr− 1, significantly larger than those obtained by Carminati and Di Donato (1999) (see their Fig. 3a), based on model ICE-3G (VM1). These subsidence rates are comparable with those driven by sediment loading in this region (Carminati and Di Donato, 1999), but significantly smaller than those predicted in other geophysical environments where sediment loading plays an important role as in coastal Louisiana (Ivins et al., 2007). The sensitivity of rates of vertical deformation to mantle viscosity further motivates our work in the following sections, which extends the sensitivity analysis

Fig. 3. Predicted Ṡ at the PSMSL sites of Cagliari (a) and Trieste (b) for upper and lower mantle viscosity values in the range 1020 ≤ η ≤ 1023 Pa·s. Shades of grey indicate sea level fall (Ṡ ≤ 0). Bottom frames show Ṡ values for two fixed values of ηUM, corresponding to the upper mantle viscosity of models RVM1 and RVM2, respectively (see solid and dashed lines in upper frames). In all these computations HLT = 90 km.

P. Stocchi, G. Spada / Tectonophysics 474 (2009) 56–68

61

Fig. 4. As in Fig. 3, but for the PSMSL sites of Malaga, Spain (a) and Hadera, Israel (b).

of Carminati and Di Donato (1999), limited to a few viscosity profiles and restricted to a single ice sheet chronology. The pattern of GIA-induced sea level change in the Mediterranean, which has been recently described by Stocchi and Spada (2007), is a peculiar feature of mid-latitude enclosed basins first noticed by Mitrovica and Milne (2002), and is broadly consistent with relative sea data from 2 kyr BP to present. Published paleo-sea level indicators from Sardinia suggest that 2 kyr BP sea level was lower than present of ∼ 1.5–2 m (Lambeck et al., 2004a). For the same time period, data from northern Adriatic indicate a sea level at ∼0.5 m below the present datum (Lambeck et al., 2004a), while archaeological indicators from Israel suggest a sea level stable at the present-day position (Sivan et al., 2004). Sedimentological evidence from the Almeria Delta suggests, for the south coast of Spain, a stable to smoothly falling sea level during since 5 kyr BP (Goy et al., 2003). This pattern of trends derived from RSL indicators is consistent with predictions in Fig. 2, indicating that while the bulk of the Mediterranean is still affected by the delayed viscoelastic relaxation of the Earth which results in a marked subsidence, coastal regions (and particularly north Africa), are close to being completely relaxed. The western margin of the Mediterranean, towards the Gibraltar Strait, is affected by the siphoning effect, which entails a long-wavelength sea level fall across the oceans (Mitrovica and Milne, 2002) that counterbalances the sea level rise within the basin. 3.2. Role of mantle viscosity and lithospheric thickness Since the Mediterranean Sea is in the intermediate field region with respect to the former ice sheets, both glacio- and hydro-isostasy may be expected to contribute to GIA-related sea level variations. However, as discussed in detail by Lambeck and Purcell (2005) and Stocchi and Spada (2007), water load effects are not dominated by the

ice-loading term, as it would be in the case in formerly glaciated regions. Since the horizontal extent of the basin exceeds the thickness of the upper mantle, we expect that present-day surface deformations (and the consequent sea level variations) may be sensitive to mantle viscosity beneath the 670 km depth transition. To test this conjecture, below we present predictions of Ṡ at the PSMSL1 stations of Trieste (Italy, TR), Cagliari (Italy, CA), Malaga (Spain, MA), and Hadera (Israel, HA), shown in Fig. 2a by white circles (see Section 3.4 for further details concerning observed rates of sea level change at these PSMSL tide gauges). Three of them (Trieste, Malaga, and Hadera) are located at the extreme margins of the basin towards North, West, and East, respectively, while the fourth (Cagliari) is closer to the bulk of the Mediterranean. As shown by predictions of reference model ICE-5G (RVM2) in Fig. 2a, amplitudes of Ṡ in response to GIA at Trieste and Cagliari differ significantly, due to the marked subsidence that characterizes the whole region of Sardinia. Among those considered, the Trieste tide gauge is expected to show the largest sensitivity to the melting of Northern European ice-sheets, since it is the closest to the region where fore-bulge surrounding the former Fennoscandian ice sheets is collapsing. Located approximately at the same latitudes, Malaga and Hadera may give information about the variance of Ṡ across the Mediterranean basin in the East–West direction. As in Section 3.1 above, we have implemented the ICE-5G model of (Peltier, 2004) but, keeping HLT = 90 km, we have generated a suite of Earth models varying both ηUM and ηLM in the range 1020–1023 Pa·s. Then, for each mantle viscosity profile, we have computed Ṡ at present time at the four tide gauges sites marked by white circles in Fig. 2.

1 Permanent Service for the Mean Sea Level. Data available on line from http:// www.pol.ac.uk/psmsl/.

62

P. Stocchi, G. Spada / Tectonophysics 474 (2009) 56–68

Results for Cagliari and Trieste are shown in Fig. 3a and b, respectively. Inspection of the diagrams clearly reveals that in these two sites Ṡ is indeed sensitive to lower mantle viscosity. Furthermore, the lack of symmetry with respect to the line ηUM = ηLM indicates a different role of upper and lower mantle viscosities. A sea level rise is always expected in Cagliari, where large-scale bending and subsidence of the ocean floor dominates and border effects are absent. Maximum values of Ṡ, close to 0.75 mm yr− 1, correspond to a stiff lower mantle with 1022 ≤ ηLM ≤ 1023 Pa·s and a “traditional” upper mantle with ηUM ∼ 1021 Pa·s. In the bottom frame of Fig. 3a, we analyze the Ṡ landscape keeping fixed ηUM to 1021 Pa·s (solid line) and 4 × 1020 Pa·s (dashed) and varying ηLM. These two upper mantle viscosity values, corresponding to those of models RVM1 and RVM2, have been obtained from global inversions of RSL observations,and are representative of volume-averaged mantle values. Ṡ predictions obtained using ηUM = 1021 Pa·s exceed those with 4 × 1020 Pa·s since a relatively soft sub-lithospheric mantle facilitates relaxation in the bulk of the basin during the early stages of loading. Fig. 3b shows that, in contrast to Cagliari, predicted GIA-induced sea level trends at Trieste may have both positive and negative sign (Ṡ ≤0 is denoted by shades of grey). However, average values are below those of Cagliari, consistent with the general pattern suggested by Fig. 3a. Maximum values of sea level rise of ∼0.5 mm yr− 1 are obtained for ηLM exceeding 1022 Pa·s and for 4×1020 ≤ηUM ≤1021 Pa·s, similarly to Cagliari. A maximum of comparable magnitude is found within the solution-space in which viscosity decreases with depth, a possibility that is generally ruled out in inverse approaches to the study of mantle rheology based upon global RSL observations, but that has been found to be compatible with observations on a regional scale (Cianetti et al., 2002). Maximum

rates of present-day sea level fall (close to −0.45 mm yr− 1) are expected for values of ηUM between 1020 and 4×1021 Pa·s and ηLM between 1020 and 1021 Pa·s. For ηUM below 3×1020 Pa·s, a sea level fall is predicted for almost all ηLM values, in the same region where a sea level rise is always expected in Cagliari (see left frame). By comparing predictions at Cagliari and Trieste, it is apparent that ηLM is a key parameter in shaping the pattern of the GIA-related present-day sea level change in these two Italian sites. At both sites, for values of ηLM between 1022 and 1023 Pa·s the Earth is far from isostatic equilibrium so that, for these viscosities, maximum values of Ṡ are predicted. As shown in Fig. 4, the peculiar sensitivity of the two Italian sites to mantle viscosity is not observed at Malaga and Hadera, where predicted rates of sea level change have generally a modest amplitude, well below 0.5 mm yr− 1, although they show a quite complex dependence upon upper and lower mantle viscosity values. For the two a priori reasonable ηUM values of RVM1 and RVM2 (solid and dashed lines, respectively) a sea level fall is basically always expected at Malaga, but it never exceeds 0.25 mm yr− 1. Remarkably, at the Israeli site of Hadera, a nearly stable present-day sea level is predicted for virtually any value of mantle viscosity, with maximum rates not exceeding 0.1 mm yr− 1. This is consistent with observations of relative sea level variations through the Holocene (Sivan et al., 2004). Comparing predictions for Malaga and Hadera (Fig. 4) with those of the Italian sites of Cagliari and Trieste (Fig. 3) shows that, despite the relatively small size of the basin, present-day sea level variations across the Mediterranean Sea in response to GIA are markedly variable. In addition, sensitivity to mantle viscosity varies considerably from place to place. At significant distances from the bulk of the basin and from the northern continental coasts (as in Fig. 4), the

Fig. 5. Present-day rates of sea level variation as a function of lower mantle viscosity, for values of lithospheric thickness in the range 60 ≤ HLT ≤ 120 km. Upper mantle viscosity is fixed to ηUM = 4 × 1020 Pa·s. Notice the different scales of plots.

P. Stocchi, G. Spada / Tectonophysics 474 (2009) 56–68

present-day sea level is nearly stable and the GIA-related signals do not vary appreciably with mantle viscosity profile. Conversely, in the bulk of the basin, as in the case of Cagliari (see Fig. 3a), the ongoing response of the mantle is still significantly operating for a broad range of a priori viscosity profiles, and results in a predicted maximum present-day sea level rise up to ∼0.8 mm yr− 1, a value comparable with current estimates of global mean sea level change due to recent climate forcing (Cazenave and Nerem, 2004). In the computations performed so far, we have only varied upper and lower mantle viscosity values, keeping the thickness of the lithosphere fixed to HLT = 90 km. However, it is known that the thickness of the elastic lithosphere may be potentially an important factor in determining the response of the Earth to unloading, especially in the presence of lateral thickness variations (Manga and O'Connell, 1995). Since spatial resolution of current 3D GIA models (e. g., Spada et al., 2006) is not sufficient to describe observed lateral variations of lithospheric thickness in the Mediterranean region (Piromallo and Morelli, 2003), here we assume a uniform thickness for the lithosphere. Furthermore, gradients in the lithospheric thickness may not be a substantial effect as shown by Davis et al. (2008) for the eastern seaboard of the United States and by Spada et al. (2006) on a global scale. Thus, using the ICE-5G chronology and keeping ηUM fixed to the RVM2 value of 4 × 1020 Pa·s, we have computed the present-day rate sea level change at the sites considered in Section 3.2 as a function of ηLM for HLT increasing from 60 to 120 km. The results, shown in Fig. 5, indicate that lower mantle viscosity is still the dominant factor, and that variations in the thickness of the lithosphere of the order of 20 km would imply small variations of Ṡ (of the order of 10− 2 mm yr− 1) in all the sites considered. These results stem from the relatively large lateral extent of the water load compared to lithospheric thickness, which enhances the sensitivity to the rheological parameters at larger depths.

63

3.3. Role of ice chronology To quantify the sensitivity of GIA-related signatures to assumptions about the ice sheets chronologies described in Section 2.2, here we perform the same analysis of Section 3.2 for chronologies ICE-3G, ICE1+ A3 and MDCR, and than we compare the results with those previously obtained using ICE-5G. Gray boxes in Fig. 6 represent, for each ice model, bounds of predicted values of Ṡ obtained by varying ηLM and ηUM in the range shown in Fig. 3 and keeping the thickness of the lithosphere fixed HLT = 90 km. Solid, dashed and dotted black horizontal lines superposed to the boxes represent Ṡ values obtained for viscosity profiles RVM1, RVM2 and RVM3, respectively (see Section 2.3). From inspection of Fig. 6a it is apparent that a present-day sea level rise is a stable feature at Cagliari, in spite of differences in the ESL and spatio–temporal features of the four ice models considered here. A maximum rate of sea level rise close to 0.8 mm yr− 1 is predicted for MDCR and positive rates are always obtained for the three RVM viscosity models considered. Negative rates, which do not exceed ∼− 0.2 mm yr− 1 are obtained only for ICE1+A3 and MDCR. In contrast to Cagliari, the site of Trieste (Fig. 6b) shows, for each ice model, predicted rates of sea level variation nearly symmetrically across the Ṡ axis. Furthermore, bounds for each ice model show a larger variance with respect to Cagliari, which indicates a sensitivity of Ṡ at Trieste to both mantle viscosity profile and ice sheet chronology. It should be observed that in our computations we do not account for the melting of the Alpine Würm glacier, which could significantly affect vertical deformations and sea level variations in southern France and in northern Italy according to previous simulations by Stocchi et al. (2005). Therefore, the results pertaining to Trieste in Fig. 6b could be subject to modifications that will be evaluated in forthcoming manuscript. Results pertaining to Malaga (c) and Hadera (d), are characterized by a general reduced sensitivity to ice sheets chronology

Fig. 6. Present-day rates of sea level change as a function of mantle viscosity, for four ice sheets models. Gray shaded bars show bounds on Ṡ with varying lower and upper viscosities in the range 1020 ≤ η ≤ 1023 Pa·s. Horizontal lines show predictions for viscosity models RVM1, RVM2, and RVM3.

64

P. Stocchi, G. Spada / Tectonophysics 474 (2009) 56–68

and small Ṡ amplitudes compared to (a) and (b). Predictions for these two sites are likely to be unaffected by the melting of the Alpine glacier, due to the relatively large distance from this ice sheet. For any of the ice models considered, results for RVM1, RVM2 and RVM3 are all positive at Cagliari, with maximum rates always obtained for the RVM3 mantle viscosity profile (dotted), i.e. the model characterized by the largest viscosity contrast across the 670 km depth interface, which implies a large isostatic disequilibrium at present time. This holds also for Trieste (b). Moving from viscosity model RVM1 to RVM2, which implies a viscosity increase of a factor of 5 between the upper and the lower mantle (see Table 1) does not alter significantly Ṡ at Cagliari, for both ICE-3G and ICE-5G. Conversely, in the case of Trieste, we observe that for the same ice models an increase of mantle viscosity contrast implies a shift from a marked sea level fall (solid) to low rates of sea level rise (dashed). When the two hybrid ice models ICE1+A3 and MDCR are employed, Cagliari and Trieste show opposite sensitivities to the viscosity increase that characterizes RVM2 relative to RVM1. At Malaga, all three mantle viscosity profiles generally imply a present-day sea level fall, with largest rates obtained for RVM1. Moving from RVM1 to RVM3 shifts the predictions towards

zero, and when MDCR is employed, a positive rate of sealevel change of ∼0.2 mm yr− 1 is predicted for RVM3. A similar trend is observed at Hadera where predicted rates are small and may attain both positive and negative values. As in the case of Malaga, the largest positive values are predicted for MDCR, the ice sheet chronology that includes a melt water pulse from Antarctica at 7 kyr BP. 3.4. Observed rates of sea level change and GIA In the sensitivity analysis above, we have been only concerned with predicted values of Ṡ as a function of mantle rheology and lithospheric thickness using previously published global models of ice sheets chronologies. To assess more quantitatively the effect of GIA upon ongoing sea level variations in the Mediterranean, we now consider Ṡ observations at the four tide gauges sites marked by white filled circles in Fig. 2a. In addition, we also consider the PSMSL sites of Tarifa (TA), Algeciras (ALG), Alicante (ALI), and Genova (GE) represented by black triangles in the same figure. Annual RLR (Revised Local Reference) mean values from these eight stations are shown in Figs. 7 and 8, respectively. Here we have

Fig. 7. Observed annual mean values of sea level (MSL) at the tide gauge of Cagliari, Trieste, Malaga, and Hadera. These sites are marked by white circles in Fig. 2a. Black and white dots show data effectively considered and those discarded because affected by missing monthly observations, respectively. Rate and standard error are shown next to each curve.

P. Stocchi, G. Spada / Tectonophysics 474 (2009) 56–68

65

Fig. 8. As in Fig. 7, but for the tide gauge of Tarifa, Algeciras, Alicante (I), and Genova. These sites are marked by black triangles in Fig. 2a.

computed long-term trends and uncertainties by straightforward least squares excluding annual values affected by missing observations, marked by white squares (details are given at http://www.nerc-liv.ac. uk/psmsl/datainfo/rlr.trends). We are aware that tide gauge time series shorter than 50 years cannot be considered reliable indicators of sea level rise or acceleration (Douglas, 1991, 1992). The inclusion of particularly short time series in our study, such as Hadera (see Fig. 7d) is motivated by the lack of other observations from the eastern coasts of the Mediterranean. The tide gauge of Trieste provides the longest time-series among those considered here, covering the period from 1905 to 2006 (see Fig. 7b). On the basis of ∼1 century of continuous observations, a secular sea level trend of 1.17 ± 0.12 mm yr− 1 is recorded at this station, representative for the stable northeastern Adriatic coast. Other long records for this region, such as those from Venice, are not considered in this study because of the large influence of anthropogenic effects on the observed trends (e. g., Carminati and Di Donato, 1999). The tide gauge of Cagliari records a sea level trend of 1.64 ± 0.40 mm yr− 1 on the basis of 30 years of observations between 1897 and 1934, while, for only 10 years of measurements (from 1995 to 2005), a trend of 1.43 ± 1.34 mm yr− 1 is derived for Hadera (Fig. 7d). The largest rate is observed at Malaga (57 years of observations,

Fig. 7c), where the sea level rose at a rate of 2.42 ± 0.43 mm yr− 1 between 1944 to 2001. This large rate of sea level change is hardly compatible with GIA, and may reflect a significant tectonic origin. Sea level trends from Tarifa (Fig. 8a) and Algeciras (8b) have not departed significantly from zero in the last ∼ 50 yr, while the low-quality time series of Alicante suggests a sea level fall in the same time span. The relatively long (but not continuous) time series of Genova (Fig. 8d) provides an average rate of sea level change (1.21 ± 0.08 mm yr− 1) that is fully consistent with the comparably long (but uninterrupted) record from Trieste (Fig. 7b). In Fig. 9a observed trends and their error bars are compared with Ṡ predictions (horizontal lines) obtained using ICE-5G and the viscosity profiles of models RVM1, RVM2 and RVM3. Predicted values of Ṡ and observations show a general coherence at the sites of Cagliari, Trieste, Genova and Hadera, while the large sea level rise observed in Malaga clearly contrasts with the GIA-related signal, which, for the RVM1 model, amounts to a sea level fall of ∼− 0.25 mm yr− 1. It is therefore worthwhile to consider other tide gauges close to Malaga. Time-series from Tarifa and Algeciras cover the same period as Malaga, but show a nearly stable present-day sea level, essentially consistent with the GIA signal. Alicante, to the northeast of Malaga, has recorded a sea level fall

66

P. Stocchi, G. Spada / Tectonophysics 474 (2009) 56–68

Fig. 9. Present-day PSMSL sea level trends compared with model predictions based on models ICE-5G(RVM1), (RVM2), and (RVM3). Site locations are shown in Fig. 2. In (a) trends are computed using all the data available from each site, while in (b), (c), and (d) we consider periods overlapping with the time-series of Trieste, the longest amongst those considered in this study.

of nearly 1 mm yr− 1 between 1952 and 1995 but the uncertainty is large due to the scatter of annual values (see Fig. 7c). This signal largely exceeds the magnitude of GIA at this site and has the opposite sign, which suggests significant local tectonic or local effects. Since the tide gauge of Tarifa, Algeciras, Malaga, Alicante, Cagliari, and Hadera cover different periods, in any case shorter than the nearly secular time window of Trieste, in frames (b), (c), and (d) of Fig. 9 we consider the latter as a reference site. Overall, this comparison indicates that, despite the sea level variations during the last century across the Mediterranean show a considerable scatter and relatively large error bars, at specific sites, the GIA-induced sea level variations may constitute a sizeable part of the observed trends. This is clearly noticeable in Cagliari and Trieste (b), where the GIA signal may amount to ∼40% of the observed trend for model RVM3 between 1905 and 1934. In the period between 1962 and 1994 (frame c) the tide gauges of Tarifa, Algeciras, Alicante, Genova and Trieste have all recorded negative trends generally below the GIA rates, but at Malaga a positive rate of 2.8 mm yr− 1 is still observed. The observed negative trends are consistent with the results of Tsimplis et al. (2005) suggesting, for the period between 1958 and 1993, a sea level fall of −0.4 to −0.7 mm yr− 1 driven by atmospheric pressure and wind effects over the Mediterranean, ultimately linked with the north Atlantic oscillation (a mean value for the meteorologically corrected trend representative of the Mediterranean for the period 1960–2000 is 0.7 ± 0.4 mm yr− 1 see Marcos and Tsimplis, 2007). Finally, between 1997 and 2005 (frame d) the tide gauges of Trieste and Hadera have recorded a nearly stable sea level and a positive rate of ∼1.5 ± 1.0 mm yr− 1, respectively. However, the large error bars reflects a significant scatter of the data on such a

short time period, whichprevents a reliable assessment of any sea level trend. 4. Conclusions The pattern of present-day GIA-induced rates of sea level change in the Mediterranean is characterized by positive values in the bulk of the basin, where the ocean load produces a marked sea floor subsidence. This pattern is similar to what is expected in other closed basins at mid-latitudes (Mitrovica and Milne, 2002; Stocchi and Spada, 2007). A similar effect is produced by sediment loading in the Gulf of Mexico, as recently observed by Ivins et al. (2007). The central portions of the Mediterranean appear to be the most sensitive to mantle viscosity profile, while the western and eastern margins of the basin show an approximately uniform response as a function of the radial viscosity variation. Among the rheological parameters considered in this study, GIA-related sea level change is particularly sensitive to lower mantle viscosity. Increasing lower mantle viscosity implies an increase of the rate of sea level rise as a consequence of the large isostatic disequilibrium attained in this case. Lithospheric thickness is a secondary factor. From our estimates, maximum rates of sea level rise due to GIA are expected along the eastern coasts of Calabria and Sicily, close to 0.8 mm yr− 1 assuming a lower mantle viscosity of ∼ 1023 Pa·s (if the traditional RVM1 viscosity profile is employed, this reduces to ∼0.5 mm yr− 1). The GIA rate of subsidence is thus significant portion of the GPS observed rates in these regions, which are in the range of 1 to 2 mm yr− 1 (Serpelloni et al., 2006). This supports previous modeling efforts of Di Donato et al. (2000), showing that peak values

P. Stocchi, G. Spada / Tectonophysics 474 (2009) 56–68

of subsidence in central Mediterranean are comparable with those due to active tectonics. The GIA induced rate of sea level rise diminishes towards the coastlines, and a stable present-day level is expected at the indented portions of the southern and eastern Mediterranean continental coasts. In the western margin of the Mediterranean, between Spain and Morocco, maximum negative values of the order of − 0.3 mm yr− 1 are expected for a nearly uniform mantle viscosity of the order of 1021 Pa·s. From our results it is apparent that the GIA parameter uncertainty overlaps with data uncertainties on the observed tide gauge trends, which makes the effect difficult to definitively detect at present-day. However, at specific tide gauge sites that are characterized by long time series, such as Trieste (northern Italy), the GIA contribution to sea level rise may be constitute a significant portion of the observed rate. At this site, depending upon assumptions about mantle viscosity, GIA may contribute up to ∼30% of the observed rate. In the bulk of the basin, where GIA produces its maximum effect regardless of mantle viscosity, this figure is close to 40% as we have verified in the specific case of the PSMSL site of Cagliari. For any a priori rheological profile and ice sheets chronology, the GIA induced rate of sea level variation in the Eastern Mediterranean is remarkably close to zero, as testified by our results pertaining to Hadera (Israel). Observed sea level variations in these regions may thus reflect the effect of tectonic movements or eustatic variations due to recent climate forcing. For all the PSMSL sites of western Mediterranean, GIA has the tendency to produce a sea level fall for all the ice chronologies considered. In this region the effect of GIA is negligible, unless the traditional viscosity profile of model RVM1 is adopted. The peculiar geographical position of Italy within the Mediterranean basin makes its coasts particularly exposed to at risk of sea level rise from GIA and consequently to erosion, which adds to the risk from the eustatic sea level variations due to ongoing climate change and anthropogenic causes. A full assessment of this issue is not our intention here. However, from the results obtained, it is apparent that the rate of GIA-induced sea level rise varies strongly with latitude along the Italian coasts, with maximum values towards the south. According to recent estimates (GNRAC, 2006), a similar pattern is broadly shared by the amount of beach erosion along the Italian coasts, which may represent a first direct evidence of the effect of GIA in the present-day coastal environment equilibrium of Italy (Stocchi et al., 2009b). Acknowledgments We thank Erik Ivins for very constructive suggestions which have greatly improved the manuscript and for the subsequent discussion. Florence Colleoni and an anonymous reviewer are acknowledged for their valuable comments. To facilitate reproducibility of the results, the numerical code employed in this study (SELEN, see http:// flocolleoni.free.fr/SELEN.html) is freely available and can be requested to GS (email: [email protected]). This work is funded by MIUR (Ministero dell' Università, dell'Istruzione, e della Ricerca) by the PRIN2006 grant “Il ruolo del riaggiustamento isostatico postglaciale nelle variazioni del livello marino globale e mediterraneo: nuovi vincoli geofisici, geologici, ed archeologici”. References Antonioli, F., Anzidei, M., Lambeck, K., Auriemma, R., Gaddi, D., Furlani, S., Orrù, P., Solinas, E., Gaspari, A., Karinja, S., Kovacic, V., Surace, L., 2007. Sea-level change during the Holocene in Sardinia and in the northeastern Adriatic (central Mediterranean Sea) from archaeological and geomorphological data. Quaternary Science Reviews 26, 2463–2486. Antonioli, F., Ferranti, L., Fontana, A., Amorosi, A., M., Bondesan, A., Braitenberg, C., Dutton, A., Fontolan, G., Furlani, S., Lambeck, K., Mastronuzzi, G., Monaco, C., Spada, G., Stocchi, P., in press. Holocene relative sea-level changes and vertical movements along the Italian and Istrian coastlines. Quaternary International. doi:10.1016/j. quaint.2008.11.008. Bennett, R.A., Hreinsdòttir, S., 2007. Constraints on relative vertical crustal motion for long baselines in the central Mediterranean region using continuous GPS. Earth Planeteray Sciences Letters 257, 419. doi:10.1016/j.epsl.2007.03.008.

67

Blanchon, P., Shaw, J., 1995. Reef drowning during the last deglaciation: evidence for catastrophic sea-level rise and ice sheet collapse. Geology 23 (1), 4–8. Carminati, E., Di Donato, G., 1999. Separating natural and anthropogenic vartical movements in fast subsiding areas: the Po plain (N. Italy) case. Geophysical Research Letters 26, 2291–2294. Cazenave, A., Nerem, R.S., 2004. Present-Day sea level change: observations and causes. Review of Geophysics 42, RG3001. doi:10.1029/2003RG000139. Church, J., Gregory, J.M., Huybrechts, P., Kuhn, M., Lambeck, K., Nhuan, M.T., Qin, D., Woodworth, P.L., 2001. Changes in sea level. In: Houghton, J.T., et al. (Ed.), Climate Change 2001: The Scientific Basis, Contribution of Working Group I to the Third Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge Univ. Press, New York, pp. 639–693. Cianetti, S., Giunchi, C., Spada, G., 2002. Mantle viscosity beneath the Hudson Bay: an inversion based on the Metropolis algorithm. Journal of Geophysical Research 107, 2352. doi:10.1029/2001JB000585. Davis, J.E., Latychev, K., Mitrovica, J.X., Kendall, R., Tamisiea, M.E., 2008. Glacial isostatic adjustment in 3–D earth models: implications for the analysis of tide gauge records along the U.S. east coast. Journal of Geodynamics 46, 90–94. Day, C., 2004. Sea-level rise exacerbates coastal erosion. Physics Today 24 (February, http://www.physicstoday.org/vol-57/iss-2/p24.html). Di Donato, G., Vermeersen, L.L.A., Sabadini, R., 2000. Sea-level changes, geoid and gravity anomalies due to Pleistocene deglaciation by means of multilayered, analytical Earth models. Tectonophysics 320, 409–418. Douglas, B.C., 1991. Global sea level rise. Journal of Geophysical Research 96, 6981–6992. Douglas, B.C., 1992. Global sea level acceleration. Journal of Geophysical Research 97, 12,699–12,706. Douglas, B.C., Kearney, M.S., Leatherman, S.P., 2001. Sea Level Rise: History and Consequences. Academic Press. 232 pp. Farrell, W.E., Clark, J.A., 1976. On postglacial sea level. Geophysical Journal of the Royal Astronomical Society 46, 647–667. Flemming, N.C., Webb, C.O., 1986. Tectonic and eustatic coastal changes during the last 10,000 years derived from archaeological data. Zeitschrift Für Geomorphologie 62, 1–29. Fukumori, I., Menemenlis, D., Lee, T., 2007. A near-uniform basin-wide sea level fluctuation of the Mediterranean Sea. Journal of Physical Oceanography 37, 2, 338–358. GNRAC (Gruppo Nazionale per la Ricerca sull'Ambiente Costiero), 2006. Studi Costieri, vol. 10. 174 pp. Goy, J.L., Zazo, C., Dabrio, C.J., 2003. A beach-ridge progradation complex reflecting periodical sea level and climate variability during the Holocene (Gulf of Almeria Western Mediterranean). Geomorphology 50, 251–268. Ivins, E.R., Dokka, R.K., Blom, R.G., 2007. Post-glacial sediment load and subsidence in coastal Louisiana. Geophysical Research Letters 34, L16303. doi:10.1029/2007GL030003. Lambeck, K., Johnston, P., 1995. Land subsidence and sea level change: contributions from the melting of the last great ice sheets and the isostatic adjustment of the Earth. In: Barends, F.B.J., et al. (Ed.), Land Subsidence. Balkema, Rotterdam, pp. 3–18. Lambeck, K., Bard, E., 2000. Sea level changes along the French Mediterranean coast from the past 30,000 years. Earth and Planetary Sciences Letters 175, 203–222. Lambeck, K., Purcell, A., 2005. Sea-level change in the Mediterranean Sea since th LGM: model predictions for tectonically stable areas. Quaternary Sciences Reviews 24, 1969–1988. Lambeck, K., Smither, C., Johnston, P., 1998. Sea-level change, glacial rebound and mantle viscosity for northern Europe. Geophysical Journal International 134, 102–144. Lambeck, K., Antonioli, F., Purcell, A., Silenzi, S., 2004a. Sea level change along the Italian coast from the past 10,000 yr. Quaternary Sciences Reviews 23, 1567–1598. Lambeck, K., Anzidei, M., Antonioli, F., Benini, A., Esposito, A., 2004b. Sea level in Roman time in the Central Mediterranean and implications for recent change. Earth and Planetary Science Letters 224, 563–575. Manga, M., O'Connell, R., 1995. The tectosphere and postglacial rebound. Geophysical Research Letters 22, 1949–1952. Marcos, M., Tsimplis, M.N., 2007. Forcing of coastal sea level rise patterns in the North Atlantic and the Mediterranean Sea. Geophysical Research Letters 34, L18604. doi:10.1029/2007GL030641. Meckel, T.A., ten Brink, U.S., Williams, S.J., 2006. Current subsidence rates due to compaction of Holocene sediments in southern Louisiana. Geophysical Research Letters 33, L11403. doi:10.1029/2006GL026300. Mitrovica, J.X., 1996. Haskell [1935] revisited. Journal of Geophysical Research 101, 555–569. Mitrovica, J.X., Peltier, W.R., 1991. On post-glacial geoid subsidence over the equatorial ocean. Journal of Geophysical Research 96, 20,053–20,071. Mitrovica, J.X., Forte, A.M., 1997. Radial profile of mantle viscosity: Results from the joint inversion of convection and convection and postglacial rebound observables. Journal of Geophysical Research 102, 2751–2769. Mitrovica, J.X., Milne, G.A., 2002. On the origin of late Holocene sea level highstand within equatorial ocean basins. Quaternary Sciences Reviews 21, 2179–2190. Mitrovica, J.X., Davis, J.L., Shapiro, I.I., 1994. A spectral formalism for computing threedimensional deformations due to surface loads. Journal of Geophysical Research 99, 7057–7073. Nakada, M., Lambeck, K., 1989. Late Pleistocene and Holocene sea level change in the Australian region and mantle rheology. Geohpysical Journal 96, 497–517. Peltier, W.R., 2001. Global glacial isostatic adjustment and modern instrumental record of relative sea level history. In: Douglas, B.C., Kearney, M., Leatherman, S.P. (Eds.), Sea Level Rise, History and Consequences. Academic Press, pp. 65–93. Peltier, W.R., 2004. Global Glacial Isostasy and the Surface of the Ice-Age Earth: The ICE5G(VM2) model and GRACE. Annual Reviews of Earth and Planetary Sciences 32, 111–149.

68

P. Stocchi, G. Spada / Tectonophysics 474 (2009) 56–68

Peltier, W.R., Andrews, J.T., 1976. Glacial isostatic adjustment, I, the forward problem. Geophysical Journal of the Royal Astronomical Society 46, 605–646. Pirazzoli, P.A., 1991. World atlas of Holocene sea-level changes. Elsevier, Amsterdam. 300 pp. Pirazzoli, P.A., 2005. A review of possible eustatic, isostatic, and tectonic contributions in eight late-Holocene relative sea-level histories from the Mediterranean area. Quaternary Sciences Reviews 24, 1989–2001. Piromallo, C., Morelli, A., 2003. PÐ-wave tomography of the mantle under the AlpineÐMediterranean area. Journal of Geophysical Research 108. doi:10.1029/ 2002JB001757. Sella, G.F., Stein, S., Dixon, T.H., Craymer, M., James, T.S., Mazzotti, S., Dokka, R.K., 2007. Observation of glacial isostatic adjustment in “stable” North America with GPS. Geophysical Research Letters 34, L02306. doi:10.1029/2006GL027081. Serpelloni, E., Casula, G., Galvani, A., Anzidei, M., Baldi, P., 2006. Data analysis of permanent GPS networks in Italy and surrounding regions: Application of a distributed processing approach. Annals of Geophysics, 49, 897–928. Serpelloni, E., Vannucci, G., Pondrelli, S., Argnani, A., Casula, G., Anzidei, M., Baldi, P., Gasperini, P., 2007. Kinematics of the Western Africa–Eurasia plate boundary from focal mechanisms and GPS data. Geophysical Journal International 169, 1180. doi:10.1111/j.1365-246X.2007.03367.x. Sivan, D., Lambeck, K., Toueg, R., Raban, A., Porath, Y., Shirman, B., 2004. Ancient coastal wells of Caesarea Maritima, Israel, an indicator for relative sea level changes during the last 2000 years. Earth and Planetary Sciences Letters 222, 315–330. Spada, G., 2001. Mantle viscosity from Monte Carlo inversion of very long baseline interferometry data. Journal of Geophysical Research 106, 16,375–16,385. Spada, G., Stocchi, P., 2006. The Sea Level Equation, Theory and Numerical Examples. Aracne, Roma. 96 pp. Spada, G., Stocchi, P., 2007. SELEN: a Fortran 90 program for solving the “Sea Level Equation”. Computers and Geosciences 33, 538. doi:10.1016/j.cageo.2006.08.006. Spada, G., Antonioli, A., Cianetti, S., Giunchi, C., 2006. Glacial isostatic adjustment and relative sea-level changes: the role of lithospheric and upper mantle heterogene-

ities in a 3-D spherical Earth. Geophysical Journal International 165. doi:10.1111/ j.1365-246X.2006.02969.x. Stocchi, P., Spada, G., 2007. Glacio and hydro-isostasy in the Mediterranean Sea: Clark's zones and role of remote ice sheets. Annals of Geophysics, 50 (6). Stocchi, P., Spada, G., Cianetti, S., 2005. Isostatic rebound following the Alpine deglaciation: impact on the sea level variations and vertical movements in the Mediterranean region. Geophysical Journal International 162, 137. doi:10.1111/ j.1365-246X.2005.02653.x. Stocchi, P., Colleoni, F., Spada, G., 2009a. Bounds on the time-history and Holocene mass budget of Antarctica from sea level records in SE Tunisia. Pure and Applied Geophysics 50. doi:10.1007/S00024-009-0488-Z. Stocchi, P., Girometti, L., Spada, G., Anzidei, M., Colleoni, F., 2009b. Post glacial readjustment, sea level variations, subsidence and erosion along the Italian coasts. Bollettino di Geofisica Teorica ed Applicata 50, 129–144. Tegmark, M., 1996. An icosahedron-based method for pixelizing the celestial sphere. The Astrophysical Journal Letters 470, L81–L84. Tsimplis, M.N., Álvarez-Fanjul, E., Gomis, D., Fenoglio-Marc, L., Pérez, B., 2005. Mediterranean sea level trends: atmospheric pressure and wind contribution. Geophysical Research Letters 32, L20602. doi:10.1029/2005GL023867. Tushingham, A.M., Peltier, W.R., 1991. Ice-3G: a new global model of late Pleistocene deglaciation based upon geophysical prediction of post-glacial sea level change. Journal of Geophysical Research 96, 4497–4523. Zerbini, S., Plag, H.-P., Baker, T., Becker, M., Billiris, H., Burki, B., Kahle, H.-G., Marson, I., Pezzoli, L., Richter, B., Romagnoli, C., Sztobryn, M., Tomasi, P., Tsimplis, M., Veis, G., Verrone, G., 1996. Sea level in the Mediterranean: a first step towards separating crustal movements and absolute sea-level variations. Global and Planetary Change 14, 1–48.