Glacio and hydro-isostasy in the Mediterranean Sea: Clark’s zones and role of remote ice sheets

Solving the sea-level equation for a spherically symmetric Earth we study the relative sea-level curves in the Mediterranean Sea in terms of Clark’s zones and we explore their sensitivity to the time-history of Late-Pleistocene ice aggregates. Since the Mediterranean is an intermediate field region with respect to the former ice sheets, glacioand hydro-isostasy both contribute to sea-level variations throughout the Holocene. In the bulk of the basin, subsidence of the sea floor results in a monotonous sea-level rise, whereas along continental margins water loading produces the effect of «continental levering», which locally originates marked highstands followed by a sea-level fall. To describe such peculiar pattern of relative sea-level in this and other mid-latitude closed basins we introduce a new Clark’s zone (namely, Clark’s zone VII). Using a suite of publicly available ice sheet chronologies, we identify for the first time a distinct sensitivity of predictions to the Antarctic ice sheet. In particular, we show that the history of mid to Late Holocene sea-level variations along the coasts of SE Tunisia may mainly reflect the melting of Antarctica, by a consequence of a mutual cancellation of the effects from the Northern Hemisphere ice-sheets at this specific site. Ice models incorporating a delayed melting of Antarctica may account for the observations across the Mediterranean, but fail to reproduce the SE Tunisia highstand. Mailing address: Dr. Giorgio Spada, Istituto di Fisica, Università degli Studi di Urbino «Carlo Bo», Via S. Chiara 27, 61029 Urbino (PU), Italy; e-mail: giorgio.spada@gmail.com


Introduction
The Mediterranean sea-level variations result from complex geodynamical, geological, and metereological processes that span a wide range of time-scales.It is nowadays accepted that the melting of the Holocene remote ice sheets and the subsequent global glacio-and hydro-isostatic readjustment profoundly modi-climate change, new geological investigations have been recently complemented by geomorphological and archaeological evidence on a regional scale (Lambeck et al., 2004a).The latter have allowed us to better constrain the age of the paleo-sea-level indicators available in the Tyrrhenian and to address the influence of tectonic motions upon the observed uplift.
While sea-level signals caused by tectonic forces and other local mechanisms may exhibit a complex spatial and temporal variability (see e.g., Carminati and Di Donato, 1999; for a regional case study), those associated with glacial isostasy are characterized by a smooth, longwavelength pattern that discloses various regions (or zones) sharing the same RSL signatures and named after Clark (Clark et al., 1978;Clark and Lingle, 1979).The mechanisms that determine the shape of the Clark's zones can be identified solving the sea-level equation that describes simultaneously the effects of ice and melt water loads, the gravitational interactions between the solid Earth and the oceans, and the time-dependent Earth's viscoelastic response (Farrell and Clark, 1976).Due to the limited computing resources, early determinations of the Clark's zones failed to resolve important details of their global pattern, including their contours and their possible fine structure in the Mediterranean Basin.The availability of highresolution solutions of the sea-level equation (e.g., Lambeck, 1993Lambeck, , 1995) ) now allows us to scrutinize details of the shapes of the Clark's zones; in particular, the recent high-resolution approach of Mitrovica and Milne (2002) has shed new light on the mechanisms that determine their shape on a global scale.Using the ice sheet chronology ICE3G (Tushingham and Peltier, 1991) and assuming a moderate viscosity increase across the mantle, Mitrovica and Milne have shown that the bulk of the Mediterranean and of other mid-latitude basins is presently subject to a sea-level rise of glacio-isostatic origin.Taking advantage of the recognized sensitivity of sea-level predictions in the Mediterranean to model parameters (e.g., Lambeck and Purcell, 2005), here we investigate the details of the Mediterranean Clark's zones for various available ice-sheets chronologies, focusing on the role of remote ice sheets.The chronology of the Antarctic ice sheet since the Last Glacial Maximum (LGM) is still affected by large uncertainties.During the last two decades evidence supported by glaciological and glacio-isostatic adjustment modeling and ice core analysis have provided sometimes divergent estimates of its contribution to global sea-level (see e.g., Kaufmann, 2002 and references therein).In the following, we will consider several plausible models for the Holocene time-history of Antarctica and we will evaluate the consequences on the RSL observations in the Mediterranean.This is motivated by the preliminary results of Stocchi et al. (2005b), who first noticed the sensitivity of RSL predictions for SE Tunisia and Gulf of Sirte to the time-history of remote ice sheets and particularly of Antarctica.
After a presentation of the methods, we will discuss the shape of Clark's zones in the Mediterranean for a suite of ice models and we will qualitatively compare predictions to observations obtained from a public domain global database.To gain insight into the results presented, we will separately consider the effects of glacio-and hydro-isostasy.In the results section, we will consider more recent and reliable RSL observations to address the role of remote ice sheets on firmer grounds.

Theory
To define the contours of the Clark's zones for the Mediterranean, we solve the sea-level equation (hereafter SLE), the linear integral equation that describes the space and time variation of the geoid due to the mass exchange between the hydro-and the cryo-sphere and to the delayed response of a visco-elastic Earth.Paraphrasing Farrell and Clark (1976), the SLE reads where O is the ocean function (O=1 over the oceans and O= 0 elsewhere), Z = OS is sea- The SLE is solved according the «pseudospectral» approach of Mitrovica and Peltier (1991), by publicly available numerical tools (Spada et al., 2004;Spada andStocchi, 2006, 2007) based on the spatial discretization scheme introduced by Tegmark (2002).Here we neglect the effects of rotation on GIA-induced sea-level variations and we assume a constant ocean function.However, for most regions the rigorous and accurate solution of the Australian National University, which incorporates time-dependent shorelines and accounts for floating ice is recommended (Lambeck et al., 2003).In the study region, the inclusion of a time-dependent ocean function may me particularly appropriate to model the sea-level rise in Northern Adriatic, as done by Lambeck et al. (2004a).
In our computations, the viscosity profile of the mantle is that implied in the «standard» deglaciation chronology ICE3G of Tushingham andPeltier (1991, 1992), with a shallow upper mantle with viscosity η SM = 10 21 Pa•s, a transition zone viscosity η TZ = η SM , and a lower mantle viscosity η LM =2×10 21 Pa•s.Following Peltier (2004), we will refer to this viscosity profile to as VM1.The shear modulus and density profiles are the same as in Cianetti et al. (2002).

Ice sheet models
In the following we will consider three distinct global ice sheets chronologies: ICE1 (Peltier andAndrews, 1976), ICE3G (Tushingham andPeltier, 1991), and ICE5G (Peltier, 2004).Frames (a), (b), and (c) of fig.1a-f show ice thickness for these models at time t=18 kyrs BP.Using ice sheets with increasing complexity and spatial resolutions will allow us to assess more robustly the sensitivity of the Mediterranean data to the Antarctic ice sheet and to review some previous results in the literature.Minor constituents of the global ice sheets distribution, such as the Alpine aggregate (Lambeck and Purcell, 2005;Stocchi et al., 2005a), will not be considered here.Frames (d), (e), and (f) of fig.1a-f show Equivalent Sea-Level (ESL) as a function of time for the ice distributions considered in this study.Here ESL is defined as level change restricted to oceanic regions, S= =S(θ, λ, t) is sea-level change at colatitude θ, longitude λ and time t, I=I(θ, λ, t) is ice thickness variation, ρ i and ρ ω are the density of ice and water, respectively, * denotes a convolution over time and over the whole surface of the Earth, G s =G s /γ where G s is the (rheologydependent) sea-level Green function, γ is mean surface gravity, and the over bar represents the average over the surface of the oceans.In eq.(2.1) S EUS represents the uniform sea-level change for a rigid, non-gravitating Earth (G s ≡ ≡ 0), with where m ICE and A o represent the ice mass variation and the area of the oceans, respectively.
In the analyses of Section 3.3 it will be instructive to consider separately ice-and oceanload induced sea-level variations.Thus, following e.g., Mitrovica and Milne (2002), we write (2.3) where S ICE , S EUS , and S OCE account for glacio-isostasy, eustasy, and hydro-isostasy, respectively, with (2.4) and (2.5) where Z represents the solution of eq.(2.1).Due to the integral nature of the SLE, however, S OCE depends from SICE through Z, so that the decoupling given by eq. ( 2.3) is artificial (Mitrovica and Milne, 2002;Lambeck and Purcell, 2005) and is used here for illustrative purposes.Once the SLE has been solved, the RSL variations at coordinates (θ, λ) are obtained by (2.6) where t BP and t p are time before present and present time, respectively.
where V is melt water volume.
Model ICE1 describes the retreat of the Holocene Northern Hemisphere ice sheets according to geological and geomorphological evidence, and stores ∼78 m of ESL at the LGM (see fig. 1a, solid).It is characterized by a relatively low spatial resolution, being specified on a grid of 5°×5°spherical quadrilaterals; the thickness of each element varies at steps of 2 kyrs.Details on how this chronology is incorporated within the pseudo-spectral method is given e.g., by Spada and Stocchi (2007).
Differently from ICE1, model ICE3G (fig.1b) has been built to improve the fit to the global RSL data in the near field assuming the VM1 rheological profile and a 120 km thick elastic lithosphere (Tushingham andPeltier, 1991, 1992).While ICE1 assumes stationary Southern Hemisphere ice masses, ICE3G accounts for the melting of Antarctica and South America, and consequently stores a larger water volume at the LGM, with an ESL of ∼113.5 m (see fig. 1c).ICE3G has an improved resolution compared to ICE1, being composed by a set of disc-shaped elements with half-amplitude of ∼1°, with thickness varying at steps of 1 kyr.Antarctica stores ∼27 m of ESL at the LGM and its volume is constant until 9 kyrs BP; its deglaciation occurs between 9 and 5 kyrs BP, when the Northern Hemisphere ice sheets had already lost most of their mass.
To fit the constraints of the glaciological studies of Huybrechts (2002) and Denton and Hughes (2002), in the recent ICE5G model of Peltier (2004) the total ESL of Antarctica has been reduced (by ∼ 50%) with respect to ICE3G (ICE5G can be downloaded from http://www.sbl.statkart.no/projects/pgs/).Such reduction is accompanied by an increased ESL for the North American ice sheet to fit constraints from VLBI and gravity observations (Argus et al., 1999).The spatial resolution of ICE5G, whose ESL is shown in fig.1f as a function of time, is improved relative to ICE35, being defined on a gaussian grid with a spacing of ∼ 0.7°.Differently from ICE3G, ICE5G describes the whole history of ice thickness during the last glacial phase (i.e., since 122.0 kyrs ago), and assumes that the LGM occurred 21 kyrs BP.Here ICE5G has been implemented assuming the VM1 viscosity profile described above.
With the purpose of evaluating the sensitivity of Mediterranean Clark zones to the timehistory of remote Holocene ice aggregates, in the computations of Section 4.1 we will also consider alternative chronologies for the Antarctic ice sheet, which differ from those implicit in ICE3G and ICE5G.

Clark zones
Here we will first qualitatively characterize the RSL variations in the Mediterranean Sea studying the shape of the Clark zones and then we will compare results based on ICE1, ICE3G, and ICE5G with observations available from the database of Tushingham and Peltier (1993).The pattern of sea-level change in the Mediterranean will also be discussed by considering separately the ocean and the ice-load components S OCE and S ICE (see eq. (2.3)), and evaluating the role of individual regional components of the two ice sheet distributions.An in-depth analysis of the sensitivity of the Mediterranean RSL variations to the ice sheets chronology will be presented in the results section with the aid of further sea-level indicators.

Pattern of Clark zones for the Mediterranean Sea
As shown in fig.2a-e, the solution of the SLE discloses complex patterns of RSL change across the Mediterranean that strongly depend upon the assumptions about the chronology of the far-field ice sheets.Model ICE3G (frame a) implies a Late-Holocene highstand that marks the end of deglaciation (5 kyrs BP) along most of the Mediterranean coasts, with the exception of Southern Italy, Greece and part of the coasts of Algeria, Lybia and of the Southern Levant, while submergence is predicted in the bulk of the basin.The regions of emergence and submergence are separated by narrow transition zones in which the sea-level nearly follows the eustatic curve.The pattern of sea-level change shown in fig.2a might be mistakenly interpreted as a match of Clark's zones VI and II, which characterize the far-field continental shorelines as a consequence of «continental levering» and the collapsing forebulge regions, respectively (Walcott, 1972;Clark et al., 1978).However, as shown by Mitrovica and Milne (2002), Lambeck and Purcell (2005) and Lambeck et al. (2004b) (and independently confirmed by our computations), submergence in the Central Mediterranean Sea mainly stems from hydro-isostasy, the direct glacio-isostatic effects of the Northern Hemisphere aggregates, characteristic of zone II, being confined to the coasts of France and Northern Italy.In general, Clark zone VI shows up as bands of offshore sea-level rise and onshore sea-level fall, with the size of the submerging areas that show significant spatial variability with the tendency to increase for concave coast lines (this is clearly visible in the map of fig.1a-f of Mitrovica and Milne 2002).When the shorelines close on themselves to define a relatively small basin, as in the Mediterranean, the area of submergence tends to cover the central portion of the basin possibly leaving narrow regions of highstand on shore (this is also observed for the Black Sea and in other mid-latitude basins in the global map of Mitrovica and Milne).To characterize the peculiar RSL pattern predicted for mid-latitude basins, where zone VI is manifest as a central submergence region contoured by a narrow highstand zone, we propose the name of Clark's zone VII.The dominance of the influence of hydro-isostasy in the Mediterranean region when all the ice aggregates of ICE3G are considered is apparent from inspection of fig.1a-f of Mitrovica and Milne (2002), in which the pattern of ṠOCE is broadly similar that of the total signal Ṡ, where the dot indicates the time-derivative computed at present time.Such pattern has remained qualitatively stable in the last 5 kyrs, after the end of deglaciation of major ice sheets (Mitrovica and Milne, 2002).
When model ICE1 is considered (fig.2b) the portion of zone VII characterized by highstands is narrowed significantly, being only present along the coasts of Spain, Northern Morocco, Tunisia and Lybia.Furthermore, the two transition regions of fig.2a are not observed.Differently from in ICE3G, in ICE1 it is assumed that the Antarctic ice sheet is stationary during the Holocene, but the two models also differ for details of the melting chronologies of the North American and Fennoscandian ice sheets (see fig. 1a-f  Up to now, we have discussed the shape of Clark zones in the Mediterranean using quite dated models for the global history of deglaciation.For instance, ICE1 is known to overestimate the extent and thickness of ice in the Baltic region (e.g., Nakada and Lambeck, 1988), while ICE3G does not match some well constrained RSL curves, such those relative to Southern France (see our discussion below).In addition, neither ICE1 or ICE3G satisfy the far-field sea-levels for the entire period since the Last Glacial Maximum (Nakada and Lambeck, 1988).In spite of their limitations, these two dated ice models have allowed us to appreciate how the shapes of the Clark zones evolve with increasing spatial resolution of the surface load.However, the question of the validity of our results with respect to more recent models remains open.To address this point, we implemented ICE5G (Peltier, 2004) within our code, and repeated all the above computations.The results, not shown here, indicate that Clark zones of ICE5G match those of ICE1+A3 very closely.A more detailed comparison of the outcomes of ICE5G with those of ICE3G is considered, in terms of RSL curves, in the following section.

Synthetic RSL curves for the Mediterranean
Figure 3a shows the location of Mediterranean sites pertaining to the publicly available database (see ftp://ftp.ncdc.noaa.gov/pub/data/paleo/paleocean/relative_sea_level/) of Tushinham and Peltier (1993) (henceafter TP), while RSL indicators are displayed in fig.3b.Recent and more reliable observations, for which information on the exact meaning of geomorphological indicators used to derive sea-level curves can be more easily obtained from the literature, will be considered in Section 4.2 below.Focussing on the TP observations at this stage is useful since they are available on a Mediterranean scale, and thus are useful for a global (although qualitative) comparison with predictions across the whole region.
The TP sites shown in fig.3a are quite evenly distributed across the Mediterranean, but since only one indicator is available for North Africa (i.e., Algiers, site 12), here we also consider the site of Djerba (Tunisia), based on the work of Jedoui et al. (1998).According to fig.3b, the only indications of a sea-level highstand in the Late-Holocene come from Beirut (9), Djerba (11) and Algiers (12).While RSL indications from Beirut and Algiers are possibly influenced by local tectonic deformations (Meghraoui et al., 2004), those from Tunisia mainly reflect glacio-isostatic RSL variations (Morhange and Pirazzoli, 2005).We observe that sites indicating an highstand are only situated along the southern and eastern continental Mediterranean coasts: this is qualitatively consistent with the pattern of Clark's zones obtained for models ICE1, ICE1+A3, and ICE5G while it appears to be at variance with predictions based upon ICE3G (see fig. 2a-e).
To assess the role of far field ice sheets, in fig.4a-i we compare individual RSL indicators from the TP database with model calculations based on ICE3G (solid lines), ICE1 (dashed), ICE1+A3 (dotted), A3 (dash-dotted), and ICE5G (grey solid).Consistently with the qualitative study of fig.2e, RSL variations driven by the melting of A3 show a clear highstand at 5 kyrs BP in all of the sites considered, with a peak amplitude varying between ∼2 m in Roussillon and ∼1 m in Messenia.For models ICE3G and ICE1+A3, the curves are characterized by a knee at 5 kyrs BP caused by their Antarctic components that with the exception of the Aegean sites of Messenia (d) and SW Turkey (e) corresponds to a sea-level highstand.When ICE1 is considered instead (dashed), RSL predictions clearly show a monotonous sea-level rise that broadly agrees with the trend of TP indicators from France, Italy and SW Turkey.In Roussillon (a), Marseilles (b) and Civitavecchia (c), the offset between ICE1+A3 and ICE3G results reflects differences in the time-history of the Northern Hemisphere ice sheets.At lower latitudes, with increasing distance from the former margins of Fennoscandia, the gap between solid and dotted curves is reduced.In the case of Djerba (h), ICE1 implies a highstand of a few centimeters at 5.0 kyrs BP, and the predictions by A3, ICE1+A3, and ICE3G almost overlap, to indicate that RSL observations at this site are virtually sensitive to the sole Antarctic component of the ice sheets distributions.Although evidence in favor of a Late-Holocene highstand at Djerba is weak due to the large (and statistically dubious) data uncertainties, the sensitivity of North Africa RSL observations to the chronology of Antarctica merits further investigations in Section 4.2.3.Predictions based on ICE5G (grey solid curves) are generally intermediate between those of ICE1 and ICE3G, as a direct consequence of the significant reduction of the total ESL of Antarctica relative to ICE3G.
Even a cursory inspection of fig.4a-i reveals that ICE3G provides globally a poor fit to the TP RSL indicators across the Mediterranean.Conversely, with the sole exception of Djerba (h) and of sites belonging to tectonically unstable regions such as Beirut (Mahamoud et al., 2005;Morhange and Pirazzoli, 2005) and Algiers (12) (Meghraoui et al., 2004), ICE1 and ICE5G broadly match the RSL trends, but this should be tested by a more rigorous misfit study.The results suggest that the Antarctic component is the main responsible of the failure of ICE3G, but it is also possible that the rheological profile VM1 is not fully suitable to describe the RSL variations in this region.We will return to these issues in Sections 4.2 and in our future investigations.

Ocean and ice-induced RSL variations
Since the Mediterranean is moderately distant from the former Late-Pleistocene ice sheets, Fig. 4a-i.Observed and synthetic RSL for some of the sites of fig.3a.Predictions are obtained by VM1 and models ICE3G, ICE1+A3, ICE1, A3, and ICE5G.RSL data and error bars are taken from the TP database with the exception of Djerba (Jedoui et al., 1998).Since the TP RSL indicators are not calibrated, this comparison has mainly a qualitative character.In Djerba the melting of Antarctica provides virtually the whole signal during the last 6 kyrs (compare solid with dash-dotted curve in frame h).
the ice-induced RSL variations in this region are not expected to dominate the ocean components, as it is the case in the near-field Clark's zones I and II (Lambeck and Purcell, 2005).To address this point, using the rheological model VM1, in fig. 5 we separately consider the RSL variations of glacio-and hydro-isostatic origin at Roussillon, Jaffa and Djerba, representative of the northern, eastern, and southern coasts of the Mediterranean Basin, respectively (see fig. 3a,b).Dashed and dotted lines show ocean-and ice-induced RSL components, obtained by eq.(2.6) with S =S OCE and S=S ICE , respectively.Since S EUS =0 after the end of deglaciation, in the time window considered RSL variations are only given by ocean-and ice-load effects.For both ICE1 and ICE3G, the trend of the ice-induced component of RSL (dotted) is constant across the Mediterranean, being mostly determined by long-wavelength deformations driven by distant sources.However, they are clearly sensitive to the assumed time-history, with a sea-level fall for ICE3G (top frames) and a sea-level rise for ICE1 (bottom).As discussed below, this diametrically opposite trend is to be attributed to the sea-level changes driven by the melting of Antarctica in ICE3G.Differently from the ice term, the ocean term in fig. 5 (dashed) shows a significant spatial variability, with a clear sea-level fall in Roussillon and Djerba, and a moderate sea-level rise in Jaffa.Such variability is caused by the sensitivity of this component of RSL to local effects related to the irregular shape of the shorelines (e.g., Mitrovica and Milne, 2002;Lambeck and Purcell, 2005).The trend of the hydro-isostatic term is the same for both ice sheets, but its amplitude tends to be larger for ICE3G, due to the effect of the extra ocean load driven by Antarctica.In fig.6 we have repeated the computations above for the site  of Djerba using the recently published deglaciation chronology ICE5G (Peltier, 2004).The reduced ESL of Antarctica in ICE5G reduces the highstand amplitude relative to ICE3G (compare with fig.5c); the emergence following the peak is now controlled by the ocean component of RSL, to indicate a dominating role of continental levering at this site and for this specific surface ice load.
In fig.7 we consider individual contributions from major Late-Pleistocene ice sheets (North America, Fennoscandia, and Antarctica), neglecting other minor constituents.The results illustrate how the distance from the former ice sheets affects RSL curves.For Roussillon, the sea-level rise produced by the relatively nearby Fennoscandian component of ICE3G is counteracted by Antarctica and North America to produce a marked Late-Holocene highstand.This is only barely visible in Jaffa, since during the last 5 kyrs the melting of the Northern Hemisphere aggregates almost exactly compensates the effects from the Southern Hemisphere.Finally, for Djerba, in the last 6  kyrs, the North American and Fennoscandian ice sheets have produced equal but opposite trends for both ICE1 and ICE3G that make the RSL observations from the coasts of Tunisia particularly sensitive to the deglaciation of Antarctica, as anticipated in Section 3.2.In fig.8 the computations of fig.7c are repeated using ICE5G.In Djerba we still observe a substantial cancellation of the effects of North America and Fennoscandia that shows the robustness of the ICE3G results and motivates further data acquisition from Tunisia, which is underway.

Results
Assuming ICE3G and VM1 as a priori models for the Late-Pleistocene ice sheets and for the Earth's rheological profile, respectively, in the previous sections we have confirmed that Antarctica significantly affects the Holocene sea-level variations in the Mediterranean.However, from the results obtained (see in particular fig.4a-i) it is apparent that ICE3G qualitatively provides a poor fit to the field observations in this area, which are better reproduced by ICE1, which assumes a stationary Antarctic ice sheet throughout the Holocene, and by ICE5G, which implies a sizeable reduction of the ESL of Antarctica relative to ICE3G.To evaluate to what extent recent Mediterranean data might be useful to constrain the history of far-field ice sheets, here we will first review some existing deglaciation histories of Antarctica and then we will implement them in model ICE3G.Our predictions will be compared with field observations from various locations digitalized from the compilation of Pirazzoli (1991Pirazzoli ( , 1996) ) and from more recent sources.With respect to the TB database, employed so far in our study to describe qualitatively the RSL trends across the Mediterranean, the sea-level indicators used in Section 4.2 below provide a better spatial coverage and are more reliable in a number of aspects.

Three ice models for Antarctica
In the following, we will modify the original ice model ICE3G by including the three dis-tinct melting chronologies for Antarctica that in fig. 9 are denoted by S, G, and D, respectively.All of them contribute 14 m of total ESL (i.e., one half of the ICE3G value).Considering the distance of the Mediterranean from Antarctica, this ice sheet is modeled by a disc load of rectangular cross-section having an half-amplitude of 20 and a thickness of ∼355 m at the LGM.To fit the ∼113 m lowstand predicted by ICE3G in the far-field sites at the LGM (see fig. 1d), we have increased the volume of the Northern Hemisphere ICE3G components at this epoch keeping their isochrons unaltered.The global chronologies so obtained are named ICE3G(S), ICE3G(G), and ICE3G(D), respectively.Evidence in support of the chosen ESL for Antarctica comes from a number of sources.In his review on the volume of Antarctica at the LGM, Bentley (1999) proposed an ESL in the range of 6.1-13.1 m that matches the 12 m estimate of Huybrechts (1992) based on glaciological modeling.Through a 3D thermo-mechanical modeling, recently Huybrechts (2002) has derived an ESL in the range of 14-18 m, consistent with the value of 14 m obtained from geologic constraints by Denton and Hughes (2002), whereas from an ice-dynamical approach Bintanja et al. (2002) estimated that Antarctica has contributed ∼5 m of ESL since the LGM.These estimates are consistently smaller (by at least a factor of 2) than former values based on the classical reconstruction of Denton and Hughes (1981), and of the ICE3G figure (∼27 m).As we will show in the following, a reduced ESL for Antarctica improves the fit with the RSL observations in the Mediterranean, mainly because the amplitude of the Late-Holocene highstand is significantly reduced.
The simplest of the three time-histories considered for Antarctica (solid curve labeled by S in fig.9) is characterized by a constant rate of melting between 12 and 5 kyrs BP, when geological evidence indicates that deglaciation was complete (Goodwin, 1996).However, it has been suggested that some additional melt water was still added to the ocean during the last 6 kyrs.Since the main source of this Late-Holocene additional global sea-level rise of ∼3 m is assumed to be the Antarctic ice sheet (Nakada and Lambeck, 1988), we have also implemented a delayed (D) melting phase (dotted).This second chronology approximately follows S until 7 kyrs BP but subsequently the rate of deglaciation decreases until the ESL vanishes 1 kyr BP.The D chronology provides a Late-Holocene water release sufficient to increase sea-level by about 3 m since 6 kyrs BP (Lambeck and Bard, 2000), while its contribution since 3 kyrs BP is less than 1 m (Fleming et al., 1998).Support for a delayed model of deglaciation also comes from the work of Stone et al. (2003), who analyzed surface exposure ages of glacial deposits in the West Antarctic ice sheet.
Various pieces of evidence indicate a midto Late-Holocene sea-level highstand in the South Pacific, Indian Ocean and in parts of the Northern Atlantic and Pacific Oceans.The subsequent sea-level fall is generally attributed to «ocean siphoning» (Mitrovica and Milne, 2002).However, based on glaciological and geological evidence, Goodwin (1998) suggested that a Late-Holocene increase of the Antarctic ice volume may be partly responsible for this sea-level fall.According to this hypothesis, the expansion of mountain glaciers, ice sheet margins and the thickening of the ice sheet in-terior could account for ∼1.0 m of the sea-level fall on mid-oceanic islands.Hence, differently from D, the G chronology, shown by a dashed line in fig.9, is characterized by a Late-Holocene re-advance causing a general eustatic sea-level fall of ∼1.0 m since 5 kyrs BP, with ∼0.7 m of sea-level fall between 4 and 2 kyrs BP (Goodwin, 1998).

French coasts
From the collection of sea-level field observations of Pirazzoli (1991), we have borrowed data from nine sites along the Golfe du Lion and the coasts of Corsica, ranging from Roussillon (1) to South Corsica (9); sites locations are shown by filled circles in fig.10.In the original compilation of Pirazzoli (1991), only two sites show some indication of a Late-Holocene highstand, with an emergence of ∼2.0 m near Cap Romarin between 5 and 4 kyrs BP (Aloisi et al., 1978), and a highstand in the range of 2 to 4 m ∼4 kyrs BP deduced from undated beach deposits in the Nice area and in the surroundings (Dubar, 1987).From a reexamination of the available record, Lambeck and Bard (2000) concluded that along the French coast Holocene sea-levels have never exceeded the present level, the highstand being marked in the erosional notch at cap Romarin of pleistocenic age (Laborel et al., 1998).Some further evidence against a Late-Holocene highstand comes from archaeological excavations of the ancient harbor of Marseilles, which have provided a new set of high-precision data for the past 4 kyrs (Morhange et al., 2001), and from the preservation of half-submerged Paleolithic paintings on a wall of the Cosquer cave near Marseilles (Vouve et al., 1996), showing that during the Holocene the sea-level never never exceeded its present-day level.
In our previous calculations of fig.2a-e, ice models ICE3G and ICE1+A3 predicted a Late-Holocene emergence along the Mediterranean coast of France.As shown in fig.11, different results are obtained when the three modified Antarctic chronologies of fig. 9 are implemented in ICE3G.The solution for ICE3G(S) (solid lines) show a sea-level highstand of ∼1 m at 5 kyrs BP both in Roussillon and in the Rhone Delta region.Due to the reduced ESL of ICE3G(S) and to the 2 kyrs anticipation of the melting inception of its Antarctic component (see fig. 9), the maximum transgression is diminished by about a factor of 2 with respect to ICE3G (compare with fig.4a-i).The amplitude of the highstand diminishes eastward until it vanishes at Port Cros (site 6); in Corsica (sites 7-9) the highstand is canceled by the submergence that characterizes Clark zone VII in the bulk of the Mediterranean.Model ICE3G(S) generally provides a poor fit to the data.The Late-Holocene re-advance of the Antarctic ice sheet, implied in the ICE3G(G) chronology (dashed lines), enhances the development of highstands.Those with largest amplitude (∼1.5 m) are observed at Roussillon (1) and in the Rhone Delta (2).The predictions for this hypothetical chronology of Antarctica systematically stand above those based upon ICE3G(S) (solid) and are in clear disagreement with observations from Southern France.ICE3G(D) (dotted) implies a smooth sea-level rise through the Late Holocene, similar to that predicted by ICE1 and ICE3G-A3 (see fig. 2c,d) and a subsequent barely visible sea-level fall ending between 2 and 1 kyrs BP.The general agreement of ICE3G(D) predictions with sea-level observations from this region is apparent.

Tyrrhenian and Northern Adriatic coasts of Italy
The pioneering investigations of Alfieri and Caputo (Schmiedt, 1972) have shown the importance of the archaeological remains as sealevel indicators of past sea-levels along the Tyrrhenian coasts of Italy.For this region, they estimated a sea-level rise of ∼1.7 mm yr −1 between 600 BC and 100 AD; in particular, dating Roman fish tanks and submerged harbors, they showed that between 100 years BC and 100 years AD the sea-level was ∼1.0 m lower than present.Subsequently, the radiocarbon-based RSL curve of Antonioli and Frezzotti (1989) evidenced, for the southern coasts of Lazio, a sea-level similar to the present one between 7 and 5.4 kyrs BP, followed by a slight oscillation below the present level (Pirazzoli, 1991), while the cumulative Tyrrhenian sea-level curve of Alessio et al. (1994) 11), chronologies ICE3G(S) and ICE3G(G) do not improve the fit with the observations.It should be remarked, however, that since the ice models scrutinized here do not account for the effects of melting of the Alpine ice-sheet, a bias may be present in our computations (see Stocchi et al., 2005a).In addition, in the case of Northern Adriatic, assuming a constant ocean function may provide a poor approximation due to the shallow water environment.More refined modelizations can account for time-evolving shorelines (e.g., Lambeck et al., 2004a).4.2.3.Gulf of Gabes, Tunisia Paskoff and Sanlaville (1983) proposed a tentative sea-level curve for South Tunisia in the last 8000 years (see also Pirazzoli 1991).The fluctuating RSL curve of Paskoff and Sanlaville, relative to Djerba (site 18) and reproduced in fig.13 by filled squares, shows a maximum emergence of reproduce ∼2.0 m between 6 and 5 kyrs BP and a second minor highst and ∼3 kyrs BP.The subsequent study of Jedoui et al. (1998) evidenced two fossilized bioclastic beaches at different elevations (see the diamonds in fig.13a).The older palaeobeach deposit (6.4-4.3 years BP) is found at an elevation of about 40 to 100 cm above the present sea-level, while the younger, dated 1850 years BP, at present sea-level.Uncertainties on the interpretation of past sealevel indicators from this area have been discussed by Pirazzoli (1987) and Lambeck et al. (2004a).Recently Morhange and Pirazzoli (2005)   In fig.13 we compare the RSL observations from Djerba with predictions based on the three modified ICE3G chronologies described in Section 4.1.While ICE3G(S) (solid lines) and ICE3G(G) (dashed) support the Late-Holocene submergence suggested by the data of Morhange and Pirazzoli (2005) and Paskoff and Sanlaville (1983) in the last ∼5 kyrs, ICE3G(D) shows a monotonous emergence in contrast with observations.However, since reliable confidence bands for dotted curves are lacking, a rigorous evaluation of misfit is not possible.Provided that uncertainties reported by Jedoui et al. (1998) are reliable, fig.13 shows that none of the three time-histories considered is indeed in contrast with the observations.While ICE5G (grey curve) is broadly consistent with Jedoui et al. (1998), it does not account for the neat highstand suggested by Paskoff and Sanlaville (1983) and Morhange and Pirazzoli (2005), nor with its tim-ing.In every instance, from the results obtained for this region it is clear that improved RSL observations for the coasts of Tunisia could significantly contribute to constrain the time-history of Antarctica in the last 6 kyrs.The cancellation of the effects from the Northern Hemisphere ice sheets that occur here once again is manifest observing that ICE3G(D) (dotted) basically reproduces the effect from the sole Antarctic component of this ice sheet (D, dash-dotted).

The Levant Sea, Israel
Evidence for Late-Holocene sea-level higher than the present in Israel has been reported by Sneh and Klein (1984) and by Raban and Galili (1985) who derived two similar sea-level curves for the site of Dor, showing sea-level fluctuations of over 2.0 m of amplitude (Pirazzoli, 1991) as portrayed in fig.14b.In a subsequent study, Nir and Eldar (1987) proposed a curve characterized by small oscillations for the last 2.5 kyrs.From the investigation of submerged archaeological remains along the continental shelf of Israel, between Haifa and Atlit, Galili et al. (1988) obtained the RSL curve shown by a dashed line in fig.14b, indicating a monotonous and smooth sea-level rise until ∼1 kyrs BP (Pirazzoli, 1991).As shown in frame (c), archaeological evidence from the sites of Tel Nami (20), Dor (21), Michmoret (23) and Yavne Yam (24) indicate that 6 kyrs BP sea-level was ∼4 m lower than today and that it reached the present-day level between 3 and 2 kyrs BP (Sivan et al., 2001).The observational limits derived from the coastal water wells in Caesarea Maritima (site 22) have extended the record of the Late-Holocene sea-level change to 1300 AD (Sivan et al., 2004) and suggest that during the Byzantine period sea-level was higher by ∼30 cm than today.The RSL curve obtained using ICE3G(G) (dashed curve in fig.14c) shows a highstand between 5 and 4 kyrs BP, followed by emergence and lastly by a negative oscillation.When the ICE3G(S) and the ICE3G(D) ice chronologies are considered (solid and dotted lines, respectively), the predicted Late-Holocene sea-level curves define a narrow band broadly consistent with the archaeological evidence until ∼2 kyrs BP.These latter predictions do not match the more recent sea-level data and thus appear to be inconsistent with the higher than present sea-level between 2 and 1 kyr BP observed at Caesarea Maritima according to Sivan et al. (2004).

Conclusions
The main results can be summarized as follows: i) The pattern of sea-level change in the Mediterranean can be described by means of two types of RSL curves.The first denotes a rising sea-level through the Late-Holocene, while the second exhibits a Late-Holocene highstand that implies a sea-level fall in the last ∼5 kyrs.The extent of the so-called Clark's regions characterized by these distinct sea-level patterns is strongly dependent upon the assumptions about the time-history of the Late-Pleistocene ice sheets surrounding the Mediterranean.When the ICE3G chronology is employed, a well developed highstand region is expected along the Mediterranean coasts, while the submergence region covers the bulk of the basin.Such peculiar pattern for closed basins has been named here as «Clark's zone VII».The assumption of a stationary Antarctic ice sheet and the enhanced effect of the melting of Fennoscandia disrupt zone VII when ICE1 is employed, leaving highstand zones along the indented coast lines of the Alboran Sea and SE Tunisia.This clearly shows the significant role played by the melting of Antarctica upon the Holocene sea-level variations in the Mediterranean.ii) Using models ICE3G and ICE5G, we have shown that along the coasts of SE Tunisia (i.e., at Djerba) the effects due to the melting of Fennoscandia are almost exactly counterbalanced by that of North America.Such fortuitous cancellation makes the highstand of ∼2 m predicted in this region only sensitive to the effects of the remote Antarctica ice sheet.A highstand is indeed suggested by both available observations and tentative RSL curves at Djerba, but their resolving power is not sufficient to constrain its amplitude unequivocally.Incorporating within ICE3G a suite of plausible models for the melting of Antarctica during the last 6 kyrs, all characterized by a sensibly reduced ESL at the LGM, it has been possible to fully enlighten the sensibility of RSL observations from Southern Tunisia to the details of the time-history.A Late-Holocene highstand of ∼ 1 m in this region is predicted when Antarctica is assumed to melt at a constant rate between 12 and 5 kyrs BP, and its amplitude is enhanced (nearly doubled) if a late ice re-advance is assumed.On the contrary, a delayed melting of Antarctica until 1 kyr BP is responsible for a regular sea-level rise since 6 kyrs BP.Only an improved spatial coverage and sampling frequency of the RSL data, made however difficult by the large tidal excursions in the Gulf of Gabes (Sammari et al., 2006), could help to put tight bounds on the details of the melting chronology of Antarctica.Model ICE5G(VM1) is unable to account for the amplitude and timing of the highstand suggested by the recent observations of Morhange and Pirazzoli (2005).This motivates further data acquisition in SE Tunisia and a refinement of existing ice sheets models.
iii) An assessment of the effects of mantle viscosity is beyond our purposes here.However, we have verified that an increase in η LM from 2×10 21 to 10 22 , a value more appropriate according to a number of authors (see e.g., Nakada and Lambeck, 1989) has a significant influence on the predicted RSL curves for the Mediterranean.In particular, this implies the disappearance of the Late-Holocene highstand from the coasts of Southern France, that would improve the agreement with RSL observations from this region.In general, an increase in lower mantle viscosity limits the variability of the predicted RSL curves corresponding to different assumptions about the chronology of Antarctica.Assuming a high-viscosity lower mantle, the confirmation of the existence of a Late-Holocene highstand in Djerba would indicate for Antarctica an ESL consistent with the value implicit in ICE3G, whereas an essentially stationary sealevel would be suggestive of a sensibly reduced ESL.

Fig
Fig. 1a-f.Global ice thickness distribution according to ICE1 (a), ICE3G (b) and ICE5G (VM1) (c) at 18 kyrs BP.This epoch corresponds to the LGM for ICE1 and ICE3G, while in ICE5G the LGM is reached 21 kyrs BP.The contribution to ESL of major constituents of the ice models are shown in the right frames as a function of time, where FEN, and NAM stand for Fennoscandia and Laurentide, respectively.With A3 and A5 we label the Antarctic components of ICE3G and ICE5G, respectively.
).To understand the origin of the distinct patterns in fig.2a,b, fig.2cwe consider the ICE3G-A3 chronology, in which ICE3G is deprived of its Antarctic component.The pattern obtained is strikingly similar to that of ICE1 (b), indicating that the existence of Late-Holocene highstand in the Mediterranean mainly result from the melting of the Antarctic ice sheet implemented in model ICE3G.When Antarctica is built into model ICE1 (fig.2d,model ICE1+ +A3), the configuration of Clark zones closely matches that of ICE3G (a).Significant difference between the results obtained for ICE3G and ICE1+A3, visible along the Northern Adriatic and Tyrrhenian coasts, can be attributed to differences in the time-histories of the Northern Hemisphere components of ICE3G and ICE1, with the relatively contiguous Fennoscandian ice sheet that is likely to play a major role.In terms of Clark zones, for model ICE1+A3 zone II counteracts the highstand of zone VII and merges with its core approximately North of the 42°N parallel.In

Fig
Fig. 2a-e.Shape of Mediterranean Clark zones for models ICE3G (a), ICE1 (b), ICE3G-A3 (c), ICE1+A3 (d), and A3 (e).RSL variations expected within each zone, qualitatively shown in the legend, include a monotonous submergence, a Late-Holocene emergence with a highstand of a few meters, and possibly two narrow transition zones.In all computations, we have adopted rheological profile VM1.

Fig. 5 .
Fig. 5. Ice-and ocean-load induced components of RSL for the sites of Roussillon, Jaffa, and Djerba (the sites location is shown in fig.3a), using ICE3G (top) and ICE1 (bottom).Since in ICE1 the end of deglaciation occurs 8 kyrs BP (see fig. 1a-f), its eustatic component vanishes in this time window (not shown in figure).

Fig. 6 .
Fig. 6.Ice-and ocean-load induced components of RSL for the site of Djerba (see map of fig.3a) obtained using model ICE5G (VM1).

Fig. 7 .
Fig. 7. Contributions to RSL due to individual components of the ICE3G (top) and ICE1 (bottom) global ice sheets distributions at the sites of Roussillon, Jaffa, and Djerba (see fig. 3a).

Fig. 9 .
Fig. 9. Time-histories of ESL for the Antarctic component of ICE3G according to the S, G, and D models described in the text.They are constituted by a simple disc-shaped ice element that stores an ESL of 14 m before the beginning of melting (12 kyrs BP).Predictions for modified ICE3G chronologies including these Antarctic components are shown in fig.11 and 12.

Fig. 11 .
Fig. 11.RSL observations from Southern France and Corsica according to the compilation of Pirazzoli (1991), to which the reader is referred for the original sources.Data uncertainties are taken from sources with the exception of Roussillon (1) and the Rhone Delta (2), for which a qualitative error band is given.Solid, dashed, and dotted curves pertain to models ICE3G(S), ICE3G(G), and ICE3G(D), respectively.Fig. 12. RSL observations and uncertainties from the Tyrrhenian coasts of Italy (sites 10-16) and from the Northern Adriatic (17) according to Lambeck et al. (2004a).Predictions are based on the three variants of model ICE3G already considered in fig.9.
confirmed that during the Holocene the sea-level was never above the present level.From the recently published high-quality data of Lambeck et al. (2004a) we selected those relative to the tectonically stable Tyrrhenian sites marked by open circles in fig.10.RSL observations in fig.12 clearly indicate a monotonous submergence since 8 kyrs BP with no evidence of highstands.With the exception of the Northern Adriatic (site 17), ICE3G(D) (dotted) systematically overestimates RSL observations, but correctly reproduce their trends.In qualitative agreement with the results obtained for Southern France (see fig.
published a tentative sea-level curve for Southern Tunisia based on new indicators collected between the Gulf of Gabes and the Libyan border.The curve, reproduced by open circles in

Fig. 13 .
Fig. 13.RSL observations for Djerba according to various sources, and predictions using ice models ICE3G(S) (solid line), ICE3G(G) (dashed), ICE3G(D) (dotted), and ICE5G (solid grey).A dash-dotted curve shows results obtained assuming that the Antarctic ice sheet model D of fig. 9 is the only active load since the LGM.

Fig
Fig.14b,c.In (b) we have reproduced tentative RSL curves from the existing literature, while in (c) we compare the data ofSivan et al. (2001Sivan et al. ( , 2004) ) with results based on the same models previously considered in fig.13.