New estimates of secular sea level rise from tide gauge data and GIA modelling

Summary

During the last three decades, at least 30 independent estimates of the secular global mean sea level rise (GMSLR) have been published, based on sufficiently long tide gauge (TG) records. Despite its apparent simplicity, the problem of GMSLR is fraught with a number of difficulties, which make it one of the most challenging questions of climate change science. Not surprisingly, published estimates show considerable scatter, with rates ranging between 1 and 2 mm yr⁻¹ for observations on the century timescale. In previous work, the importance of Glacial Isostatic Adjustment (GIA) upon the assessment of the GMSLR has been clearly demonstrated. In particular, starting from the 1980s, GIA models have been routinely employed to decontaminate TG observations from the effects of melting of the late-Pleistocene ice sheets, to fully highlight the sea level variations driven by climate change. However, uncertainties associated with the Earth's rheological profile and the time history of the past continental ice sheets can propagate into the GIA corrections. After revisiting previous work and estimates, we suggest a significant modification of the criteria for the selection of the TGs which are most suitable for the robust assessment of the secular GMSLR. In particular, we seek a set of TGs for which GIA corrections are essentially independent of the parametrization of the rheological profile of the Earth's mantle and of the detailed time chronology of surface loading. This insensitivity is established by considering predictions based upon three GIA models widely employed in the recent literature (namely, ICE–3G, ICE–5G and the one developed at the Research School of Earth Sciences of the National Australian University). Applying this approach and selection criteria previously proposed in the literature, we identify a set of 22 sufficiently evenly distributed TGs. By simple statistical methods, these records yield a ‘preferred’, GIA-independent GMSLR estimate since 1880, namely 1.5 ± 0.1 mm yr⁻¹ (rms = 0.4 mm yr⁻¹, wrms = 0.3 mm yr⁻¹). This value is consistent with various previous estimates based on secular TG observations and with that proposed, for the 20th century, by the Intergovernmental Panel on Climate Change (IPCC) Fourth Assessment Report (1.7 ± 0.5 mm yr⁻¹).

Sea level change, Global change from geodesy, Rheology: mantle

1 Introduction

In a global change scenario, current sea level rise is perhaps one of the most debated environmental problems (see the IPCC Fourth Assessment Report, Bindoff et al. 2007). On the century timescale, sea level variations mainly result from ocean water volume changes associated with two major climate-related processes, namely land-ice melting and water density variations (Cazenave & Nerem 2003). These processes are expected to have a central role in governing sea level changes also during the coming decades (e.g. Milne et al. 2009). The first implies a net mass variation of the oceans and mainly involves glaciers, ice caps, ice sheets and other land water reservoirs, but also hydrological variations caused by anthropic activity (mass contribution). The second, which is associated with ocean-water density and salinity changes, does not imply any mass variations (steric contribution). It is now recognized that both the mass and the steric contributions, beside varying the global mean sea level rise (GMSLR), also imply significant geographically non-uniform sea level fluctuations, which have a direct effect on coastal hazard and society. For these reasons, a major challenge for the next IPCC assessment will be to incorporate the spatially complex sea level variability into regional scenarios for the future sea level rise (Gehrels 2010). In spite of the growing role of sea level ‘fingerprints’ as tools for characterizing regional sea level variations (Mitrovica et al. 2001; Tamisiea et al. 2001), the assessment of spatial averages, hence of GMSLR, has been and continues to be of some importance. In fact, ocean-averaged sea level changes directly provide information on the net mass variations of continental ice masses (see e.g. Spada et al. 2012) and ocean-water density variations, which cannot be inferred by regional sea level observations. Tide gauges (TGs) have been historically employed to monitor sea level rise. Although the sparse distribution of TGs and local movements of tectonic and glacio-isostatic origin may hinder a precise evaluation of spatial averages, starting from the first decades of last century considerable attention has been paid to the assessment of long-term GMSLR from instrumental records. As discussed in Section 2.2, this quantity is generally defined as the simple mean of trends obtained from long-period TG records.

Since the seminal work of Gutenberg (1941), various GMSLR estimates have been produced, essentially based on the global TG records now supplied by the Permanent Service for the Mean Sea Level (PSMSL, see http://www.psmsl.org/). The available estimates obtained from TG observations are presented in Table 1, which also shows the period of observation and basic information on the methods employed to estimate the GMSLR (hereafter abbreviated by μ). From Fig. 1(a), showing the chronological sequence of the estimates, it is apparent that an increased interest into the sea level rise problem took place in the early 1980s, motivated by the awareness of the possible effects of greenhouse gases on the Earth's climate (Gröger & Plag 1993). Fig. 1(b) indicates that most of the studies are referred to the period 1900–1990, apparently reflecting a decreased interest in TGs observations after the introduction of satellite altimetry methods. In contrast to TGs, which are placed along shorelines mostly in densely populated areas, sea level observations by satellite radar altimetry permit a nearly complete sampling of the sea surface. Although altimetry data require a GIA correction, they are not affected by local vertical crustal movements that contaminate the TG observations (see e.g. Nerem and Mitchum 2001). The altimetry record is short compared to that of TGs, which usually provide records of several decades and in some cases of more than one century (this is detailed in Section 2.1). In fact, high-quality altimetry-based sea level measurements (with an interannual precision better than 1 mm yr⁻¹) are only available since the beginning of the 1990s with the Topex/Poseidon and ERS-1 missions (Cazenave et al. 1998). Consequently, a substantial portion of the satellite-determined sea level trends are possibly still affected by short-term decadal signals, unrelated to the long-term variations that are expected from global warming (Nerem 1995).

Table 1

Previous GMSLR estimates (data are displayed in Fig. 1). See Pirazzoli (1993), Gröger & Plag (1993), Gornitz (1995) and Douglas (2001) for similar summary tables.

Figure 1

(a) GMSLR estimates from TG data according to previous authors (see list in Table 1). GIA-corrected and uncorrected estimates are shown by solid and open circles, respectively. Our preferred estimate (estimate n. 15 of Table 7, with μ = 1.5 ± 0.1 mm yr⁻¹, rms = 0.4 mm yr⁻¹, wrms = 0.3 mm yr⁻¹) is shown by a grey dot. (b) GMSLR estimates proposed in the literature based on the analysis of TG observations. The extent of each rectangle shows the time period considered for each estimate, while its width shows its uncertainty. Among the GMSLR estimates considered in (a), only those corrected by GIA are considered in (b).

Traditionally, GMSLR is determined by a mean rate of sea level change, which can be either based on trends at individual TGs extracted from a global database (e.g. Gornitz & Lebedeff 1987; Cabanes et al. 2001) or, alternatively, on regional averages (e.g. Gutenberg 1941; Holgate & Woodworth 2004). In the following, GMSLR will be estimated by eq. (6), using individual TG trends from the PSMSL database. Following the work of Barnett (1983), empirical orthogonal functions (EOF) have been widely employed to study the spatial pattern of global sea level change, while the use of spherical harmonics analysis (SHA) appears to be limited to the study of Nakiboglu & Lambeck (1991). As shown in the last column of Table 1, different approaches have been adopted to correct the TG observations from the effects of Glacial Isostatic Adjustment (GIA). Before global GIA models were introduced in this context by Peltier & Tushingham (1989), sea level trends obtained from TG observations were often corrected by GIA trends extrapolated by geological sea level records (Gornitz et al. 1982; Gornitz & Lebedeff 1987), or not corrected at all. Since the 1990s, due to continuous models refinement, GIA corrections have been routinely applied, usually based on the ICE–3G model of Tushingham & Peltier (1991) and subsequent versions. Recently, non-traditional techniques (neural networks) have been employed to determine GMSLR (Wenzel & Schröter 2010). The application of a GIA correction is, in this case, unnecessary.

Though all estimates in Table 1 neatly point to a globally averaged sea level ‘rise’ and their central values are generally in the range between 1 and 2 mm yr⁻¹ (see also Fig. 1b), their distribution shows a significant spread. It is remarkable that according to some authors (Pirazzoli 1986; Stewart 1989; Emery & Aubrey 1991; Gröger & Plag 1993) the significance of a global sea level mean is doubtful, mainly because of the large spatial variations of the TG signals, the poor geographic distribution of the stations and various contaminating effects such as crustal deformations associated with GIA. The introduction of GIA corrections, based either on geological observations or on direct modelling, and the adoption of rigorous selection criteria for the TGs (Douglas 1991) contributed significantly to enlighten the coherence of the TG observations and to better constrain GMSLR. Furthermore, since the early 1990s, all the proposed GMSLR values have been supplied with an uncertainty estimate, though only in some of the works listed in Table 1 is it possible to ascertain its statistical significance (this is a major issue when some fundamental problems as the closure of the sea level budget are addressed, see e.g. Peltier 2009). The spread of the μ values in Table 1 is the consequence of a number of factors that may prevent direct comparisons among the outcomes of these analyses, though all of them are ultimately based on TG data. The most important include: (i) the different criteria adopted to select the TGs which are supposed to be suitable for the determination of μ, (ii) the nature of the data considered (i.e. monthly or annual means) and possible quality constraints (e.g. the exclusion of data that do not meet the PSMSL Revised Local Reference (RLR) requirements, see Section 2.1), (iii) possible pre-processing of the TG time-series (e.g. the application of moving averages to attenuate the noise level or to smooth decadal fluctuations), (iv) the length of the time period of observation considered (in Table 1 it ranges between 16 and 189 yr) and the details of averaging procedure employed to obtain μ, and finally (v) the effective application of GIA corrections to the sea level trends and the particular GIA model employed to perform the corrections. All these factors motivate the statement of Pirazzoli (1993), that is, that a ‘simple’ assessment of GMSLR is not possible, though subsequent studies (see Table 1) have shown that improvements in the robustness of the assessments are certainly not precluded.

The fundamental issue of the selection of the most appropriate TGs for global sea level analysis has been addressed by virtually all the authors listed in Table 1. Though it is clear that the quality of the trend estimate will increase with the record length, the minimum length has never been unequivocally determined. As pointed by Douglas (2001), this is certainly one of the fundamental reasons of the disagreement between the various estimates of μ proposed so far. The rule of thumb proposed by Sturges & Hong (2001) suggests that to determine meaningful statistics about sea level trend, the length of the time-series should be comparable to ∼10 times the period of the characteristic decadal oscillation. Since these are of the order of ∼20 yr (Sturges & Hong 2001), a rigorous application of this rule would reduce the number of appropriate TGs to a few units. It is now accepted that records as short as a few decades are not appropriate for the estimate of GMSLR. However, in some of the studies of Table 1, short records have been effectively employed. Perhaps not coincidentally, when this is done (see e.g. Emery 1980, who employed records with a span as short as a few decades), the value of GMSLR obtained (3 mm yr⁻¹) clearly represents an outlier in the population of the estimates. Starting from the work of Peltier & Tushingham (1989), it has been realized that a satisfactory trade-off between the length and the number of records can be obtained for records of 50 yr or longer. In the two possibly most influential papers published during the last two decades, this minimum number has been set to 60 (Douglas 1991) and 70 (Douglas 1997), respectively (in the following this last study will be referred to as D97). The minimum record length is, however, only one of several important criteria useful for estimating GMSLR from TG data. According to Douglas (1991), the time-series must (1) be at least 60 yr in length, (2) be obtained from TGs sufficiently distant from collisional tectonic boundaries, (3) have a sufficient completeness (>80 per cent), (4) be in reasonable agreement with nearby records at low frequencies and (5) not belong to previously ice-covered areas during the Last Glacial Maximum (LGM, ∼21 kyr ago) to avoid GIA contamination. Criterion (5) has been subsequently extended by D97, to exclude the TG sites from the peripheral bulge adjacent to these areas (the D97 criteria are summarized in the left part of Table 4). Applying the above criteria, D97 has used trends from 23 PSMSL TGs (listed in Table 2) for the time period 1880–1980 to determine a GMSLR of μ′ = 1.8 ± 0.1 mm yr⁻¹, where the prime indicates that a GIA correction has been applied (based, in this case, on model ICE–3G of Tushingham & Peltier 1991) and the uncertainty denotes the standard deviation of the mean.

Table 4

The D97 selection criteria compared with those proposed in this work.

Table 2

Computed rates of sea level change for the D97 set, compared with those determined by D97. The average span of the time-series is 95 yr (in D97 it was 83 yr), while the average number of valid yearly RLR records in each series is 88 yr. Note that in D97, uncertainties on the individual trends r_k were not provided. The average completeness c_k of the time-series for these TGs (i.e. the average of the ratio ) is 92 per cent. According to our computations, for this set of TGs the GMSLR is μ = 1.6 ± 0.1 mm yr⁻¹ (rms = wrms = 0.4 mm yr⁻¹). For reference, GIA corrections corresponding to model ICE–3G are shown in the last column. The 16 TG sites marked by a star also belong to the D97R set (see Section 3.3).

Though those proposed by D97 probably constitute the ‘best possible’ criteria for the selection of TGs useful to the assessment of GMSLR, a couple of caveats are apparent. First, the requirement (2) (distance from collisional tectonic boundaries) is motivated by the possible contamination of the TG signals that could result from vertical movements (and geoid height variations) at subduction zones. However, from forward modelling we know that vertical movements are also possible in transcurrent tectonic environments (e.g. Nostro et al. 1999). Since the D97 set contains a significant number of TGs located along the North American West Coast, the assessment of GMSLR could be affected, by an unknown amount, by tectonic deformations. It is unclear whether tectonic contributions in this region have been accurately accounted for (or avoided) in previous studies. As pointed by Nakada & Inoue (2005), it is likely that the tectonic trends at the TG of San Francisco (traditionally included in most of previous studies for its considerable record length and completeness), have changed as a consequence of the coseismic deformation associated with the 1906 M = 8.3 earthquake. This issue could be resolved, in principle, by modelling efforts aimed at estimating the co- and post-seismic displacements associated with the seismic cycle. However, the episodic nature of the earthquake energy release and possible viscoelastic effects could potentially prevent the accurate assessment of a meaningful tectonic sea level trend at TGs (Nakada & Inoue 2005). The second caveat concerns criterion (5) of D97, that is, the exclusion of TG stations strongly affected by GIA in the regions previously covered by the ice during the LGM and at their immediate periphery. According to D97, the motivation for the exclusion essentially resides in the difficulties inherent in the calculation of the GIA correction. The main reason to avoid sites on the peripheral bulge was that the response of such sites is strongly dependent on lower mantle viscosity, and that with a change in this viscosity one could reduce the systematic trends in TG rates along the US east coast (Davis & Mitrovica 1996). Using model ICE–3G, D97 argued that since the correction at some sites can be comparable with the value of GMSLR a posteriori assessed, even a relatively small error in the evaluation of the correction could impact the assessment. A problem with this approach is that different GIA models will provide different spatial uplift patterns, according to the viscosity profiles assumed and to the spatiotemporal details of ice melting (e.g. Nakada & Lambeck 1987). As we will discuss in the body of the manuscript, since the geometry of the forebulge region is sensitive to the choice of the GIA model, the selection criterion (5) will be affected accordingly. The a priori exclusion of the TGs located in previously ice-covered areas certainly represents a prudent approach when a unique GIA model is considered, as it is the case for D97. However, this could entail a rejection of long and sufficiently complete records, like those located in the northern European regions.

Similar to a considerable number of previous investigations, listed in Table 1, this study aims to determine the secular GMSLR by means of sufficiently long and properly selected TG observations. Despite the poor spatial distribution and the contamination of crustal movements, on these timescales TGs represent the only technique for directly observing sealevel variations (see e.g. Douglas 2001). Our focus is on the time period from year 1880 to present, for which revised yearly PSMSL records are today available (relative to D97, the time-series employed here are ∼15 yr longer). The main purpose of this study is to propose a new approach to the problem of the GMSLR assessment, aimed to fill gaps in previous works and to establish more robust estimates. This is accomplished by a revision of the selection criteria introduced by D97. The most significant implies the substitution of criterion (5) of D97 by a new requirement based on the values of GIA corrections obtained by a set of independently developed models. The idea behind is that suitable TGs are those for which GIA corrections should be largely independent from the details of GIA modelling. Along these lines, new and revised estimates for the GMSLR obtained by this new set of preferred TGs will be presented and discussed.

The paper is organized as follows. In Section 2 we introduce the TG observations (Section 2.1), we describe the basic statistic used for the analysis of the trends of the TG time-series (2.2), and we provide a first assessment of GMSLR from uncorrected TG data (2.3). Section 3 is devoted to the study of GIA effects on secular sea level trend. Here GIA modelling is introduced (3.1), the corrections on sea level trends are evaluated (3.2) and the modelling uncertainties are discussed (3.3). Section 4 presents a broad discussion in which we introduce a new criterion for the selection of TGs suitable for the assessment of GMSLR (4.1) and we present a set of TGs resulting from the application of this criterion (4.2). This is followed by the presentation of a suite of new GMSLR estimates (4.3) and by a discussion of the possible role of tectonic deformations (4.4). Finally, we briefly discuss the problem of secular sea level acceleration and recent sea level rise (4.5), and we draw our conclusions (5).

2 TG Data and Uncorrected GMSLR Estimates

2.1 TG observations

A broad review on the availability and quality of TG observations during the last two centuries has been provided by Douglas (2001), who discussed the errors in the observations and the accuracy of the various GMSLR estimates previously obtained. Critical assessments of the quality of TG data, particularly in connection with the poor geographical distribution of the instruments, the limited length of the time-series and the contaminating local vertical movements have been given by Emery & Aubrey (1991) and Gröger & Plag (1993). These issues ultimately explain the variance of GMSLR estimates obtained in the literature (see Table 1), which are based on different criteria for consideration or exclusion of specific time-series. As pointed by Douglas (2001), in view of the large number of factors that affect the TG observations and of the very broad scientific and socio-economic implications of global and regional sea level variations, international cooperative efforts are necessary.

In this work, we employ observations obtained from the PSMSL, which collects TG sea level data from individual national authorities since 1933 (Woodworth & Player 2003; PSMSL 2012). Currently, PSMSL distributes Metric and the RLR data, which correspond to different quality standards. Metric data are monthly and annual observations directly received from the national authorities. To construct time-series of sea level measurements, monthly and annual records are reduced to a common benchmark based on the history of the TG datum; the data revised in this way constitute the RLR data set. Here we use RLR annual mean data, which constitute the most appropriate record for analyses of long-term sea level rise (Spencer & Woodworth 1991). A visualization of each individual time-series is available from page http://www.psmsl.org/data/obtaining/rlr.annual.plots.

The distribution of the 1213 sites provided TG time-series including at least three yearly records since year 1880, hereafter referred to as the ‘ALL set’, is shown in Fig. 2(a). Fig. 2(b) shows the location of the 23 D97 TGs, and the corresponding RLR annual records are shown in Figs 3 and 4, respectively. For consistency with Douglas (1991) and D97, here we will only employ observations available since year 1880. Furthermore, following the suggestions of PSMLS (see http://www.psmsl.org/data/obtaining/rlr.php), for our GMSLR estimates we will rigorously consider only RLR observations. Apparently, no distinction between RLR and metric data was made in D97.

Figure 2

(a) Location of the 1213 PSMSL TGs with at least three yearly RLR observations since 1880 (ALL set). (b) Distribution of the 23 sites considered in D97 (the names of the sites and other relevant data are listed in Table 2).

Figure 3

RLR annual time-series for tide gauges located in the regions of the English Channel (1–2), Atlantic (3–5), Mediterranean Sea (6–8), New Zeland (9–11) and the Pacific (12), according to the selection of D97. Thick segments at the tip of some of the time-series mark yearly observations that were not available when D97 was published.

Figure 4

The same as in Fig. 3, but for the North American West Coast (13–16), Central America (17–18), South America (19–20) and SE North America (21–23), according to the selection of D97.

In Figs 3 and 4, the average sea level is subtracted before plotting and a constant shift of +150 mm has been applied for visualization purposes. Only data resulting from at least 11 months are shown, and subsequently considered in Section 2.2 for the computation of sea level trends. Neither smoothing has been performed, nor attempts made to cure the gaps in the time-series (we note that for some of them the problem of gaps is particularly severe, as for Dunedin II, in New Zeland). For all sites, the current record length is reported in Table 2. Though obviously we do not expect that an extended record, when available, would dramatically modify previous results of D97, the extension is in some case significant—for example, Balboa, Central America—and sometimes characterized by steep slope variations, as in the case of Trieste. Thus, the overall effect of these recent observations upon the long-term sea level trends is, in our opinion, worthy of investigation.

Despite the large annual (Willis et al. 2008), decadal and interdecadal fluctuations (Sturges & Hong 2001), the D97 time-series shown in Figs 3 and 4 coherently denote a long-term sea level rise. The spatial coverage of this specific subset of TGs, shown in Fig. 2(b), strongly suggests that sea level rise is indeed global, though characterized by regional variations that can be partly attributed to the effects of GIA. As noted by Douglas (1991) and D97, separating short-term fluctuations from the low-frequency variations is a difficult task, and even simple visualization of the time-series may be problematic without the aid of low-pass filters (e.g. moving averages). The fundamental reason is that TG records are characterized by a red spectrum, in which power continues to increase out to the low frequencies (Sturges & Hong 2001). Hence, establishing the amplitude of long-term sea level variations from time-series with periods of a few decades can be particularly problematic. Noise reduction (see e.g. Rompelman & Ros 1986), can be realized by a technique often employed in seismology and known as stacking (Gilbert & Dziewonski 1975). Stacking has been used by Trupin & Wahr (1990b) to study the 18.6 yr lunar tide and the 14 month pole tide using TG records; here this method is adopted to visualize the low-frequency component of sea level rise.

In Fig. 5(a), the 23 TG time-series shown in Figs 3 and 4 have been stacked in an attempt to reduce the signal-to-noise ratio and to visualize qualitatively the overall coherence of the datum at different epochs. Before stacking, the average sea level change has been subtracted by each series, but no smoothing has been applied. The stack shows no apparent sea level rise between 1880 and 1900, when only a few stations were in operation (namely, Genova, Marseille, Brest, Cascais, Fernandina and San Francisco). This supports previous observations of Nakada & Inoue (2005), indicating that modern sea level rise probably corresponded to the onset of the industrial revolution, approximately 100 yr before present. Support of this hypothesis is provided by observations based on TG time-series not belonging to the D97 set (e.g. Swinoujscie, along the southern coast of the Baltic Sea), though this datum is not consistent with that of nearby sites. Starting from ∼1905, RLR records from the sites of Trieste, Honolulu, San Diego and Buenos Aires become available. After a few years, a coherent signal is neatly emerging from the stack of Fig. 5(a), and is maintained to the present time, apparently unaffected by the systematic decrease in the number of operating TGs (and/or the number of available RLR observations) since the early 1980s. Fig. 5(a) suggests that to observe a globally coherent sea level rise, a minimum number of ∼10 high-quality TG stations is required, providing a sufficient geographical coverage. This is consistent with previous studies (e.g. D97), in which sea level trends from the 23 stations have been grouped on a regional basis, practically reducing the number of time-series down to nine (see Table 2). We note that a very limited set of TGs (seven) has also been considered by Nakada & Inoue (2005), in view of their very long record. Though from their fig. 1 it is clear that the geographical coverage is poor, their estimate of the 20th century GMSLR (1.5 mm yr⁻¹), has been found to be largely consistent with other studies based on the D97 set. The effectiveness of stacking in reducing the signal-to-noise ratio is apparent when we compare Fig. 5(a) with individual TG series in Figs 3 or 4. In the stack, the peak-to-peak amplitude of decadal and interdecadal oscillations is reduced by a factor of ∼3, leaving a low-frequency signal that does not seem to indicate, at least visually, any acceleration since the beginning of last century (Woodworth 1990).

Figure 5

(a) Stack of the 23 time-series belonging to the D97 set (see Figs 3 and 4). The gray shaded curve shows the stack for the 22 time-series belonging to the SG01 set, introduced in Section 4.2. (b) Stack for the time-series of the ALL set (). Dotted curves show the number of TGs that, at a given epoch, are operating and provide RLR sea level observations according to the PSMSL record (see second y-axis).

In Fig. 5(b), all the TGs with a minimum number of three yearly records are considered (ALL set). Clearly, in this case we end up with a more noisy stack as a consequence of the variegated and largely incoherent time-series that we are averaging. However, as soon as the number of operating TGs sharply increases to ∼200 at the beginning of the 1930s, a coherent sea level rise can be discerned. The detection of similar ‘change points’ in the time-series of sea level change is a particularly challenging problem, which has also been recently addressed by Kemp et al. (2011) on the millennium timescales. On the secular scale, it is known that abrupt variations could be associated with transients that do not imply an effective change of the long-term sea level rise. This is illustrated by Sturges & Hong (2001) for the case of the records at New York and Charleston, which suggest an apparent increase in the rate of sea level rise beginning in the late 1920s. As previously pointed by Sturges (1987), it is very difficult to argue for any statistically significant change in sea level trend in the 1920s, since the energy in the multidecadal fluctuations at periods between 40 and 50 yr prevents robust assessments. We also note that Church & White (2006) have suggested 1930 as the point in which the slope of sea level trend abruptly changes, implying an acceleration, and that Donnelly et al. (2004) detected an increased sea level rise to modern values in the late 19th century, roughly coincident with the climate warming observed in both instrumental and proxy records. Clearly, the matter of possible change points in the instrumental sea level record is not settled. Beside this, the results of Fig. 5 show that an assessment is possible from the analysis of a limited number of coherent and sufficiently long TG time-series, also satisfying specific additional requirements. This is consistent with Douglas (2001), who observed that the detection of sea level variation by TGs is indeed feasible by means of a few long records, which do not demand the dense coverage of short ones to establish a globally coherent signal. However, contrary to the pessimistic point of view of previous investigators (see e.g. Pirazzoli 1986; Gröger & Plag 1993), Fig. 5(b) shows that it can also be qualitatively observed from consideration of a sufficiently large number of globally distributed instruments with relatively short observation periods.

2.2 Sea level trends at TGs

Previous estimates of GMSLR have relied upon rates of sea level change over a specific time period, computed for individual TGs or groups of TGs (see Table 1). Consistent with the guidelines of previous investigations, to extract the background sea level trend from the records, we will use a simple linear regression. Differently from Douglas (1991) and D97, however, we do not start from monthly data and we do not perform any filtering on the time-series to reduce the amplitude of decadal cycles. Rather, following Sturges & Hong (2001), we seek a good estimate of long-term trends adopting a ‘pencil and ruler’ approach directly to annual (RLR) observations. This approach is motivated, a posteriori, by the consistency of our (GIA-uncorrected) results with D97. Furthermore, since in this section we are not concerned with possible sea level accelerations (Woodworth 1990; Douglas 1992) or decelerations (Houston & Dean 2011), a constant sea level trend is assumed. This is motivated by our analysis in Section 4.5, where we find that the sea level acceleration, evaluated comparing rates of sea level change during different time periods, is negligible.

Using basic statistics (see e.g. Taylor 1997), the rate of sea level change for the kth TG is computed by

1

where k = 1, …, N_tg, N_tg is the number of TGs considered, y_j is sea level at the year x_j (), is the number of valid annual records, and ∑_j stands for . The values of x_j and y_j are directly obtained from the annual RLR PSMSL record. At this stage, no GIA corrections are performed on the annual data or on the trends computed by eq. (1). The formal uncertainty on the computed rate of sea level change is determined building a 95 per cent confidence interval for r_k. This is done by evaluating

2

where is the average of the x_j's and is the 0.975-th quartile of Student's t-distribution with degrees of freedom (henceforth, we will only consider TG series with ). In eq. (2), the standard error of the estimate is defined as the rms of the deviations

3

(e.g. Zar 2010). We note that the uncertainties on the trends evaluated by Douglas (2001) (see his eq. 1) do not define a confidence interval, though they are found to be basically comparable to those obtained by eq. (2). Hence, to properly account for uncertainties on the rates r_k, the rate of sea level change is expressed by

4

with r_k and σ_k given by eqs (1) and (2), respectively. We note that in previous relevant studies (e.g. Douglas 1991, and D97), the uncertainties on the computed trends σ_k have not been considered explicitly in the analysis. It is however clear that they have an important role, since they can be used as weighting factors in the assessment of the uncertainties of GMSLR (see e.g. Mitrovica et al. 2001, and our eq. 9). For the subsequent discussion, it is convenient to model the observed secular rate given by eq. (1) as the combination of several contributions:

5

where k = 1, …, N_tg and the terms on the right-hand side are associated with GIA, the current melting of continental ice masses, the steric effect, tectonic motions and other factors, respectively.

The basic statistics of the TG trends are shown in Fig. 6. Mainly because here we dispose of a comparatively larger number of observations (there are 1213 TGs with a minimum number of yearly records ), details of this figure differ somewhat from similar plots presented in previous studies (see e.g. figs 3.1 and 3.11 of Douglas 2001). Their main features, however, are largely reproduced. The distribution of , shown in Fig. 6(a), appears to be bimodal, with the maximum for and a secondary peak for . This last feature, unnoticed by Douglas (2001), is due to the increased rate of installation that immediately followed the end of the Second World War. The cumulative distribution, shown by a stairs-step histogram in Fig. 6(a), indicates that most of the time-series (90 per cent +) are less than ∼60 yr long, which makes the assessment of the secular sea level trend particularly problematic (the issue is discussed by Sturges 1987), also in view of their uneven spatial coverage (see Douglas 2001, and Fig. 8). The rates of sea level change obtained by least-squares fitting of the time-series, shown in Fig. 6(b), are largely scattered, though visually it can be appreciated that they are not symmetrically distributed across zero, even for relatively small values of . Neatly positive r_k values are shown for most of the longest series (), but the values are scattered significantly. Uncertainties σ_k, shown in Fig. 6(c), neatly decrease with . Relatively precise r_k values, with error levels σ_k≤ 0.1 mm yr⁻¹ can only be obtained from some of the time-series with , as marked by the dotted horizontal line.

Figure 6

Distribution of the number of tide gauges N_tg (a), of the observed rates of sea level change r_k (b) and of their uncertainties σ_k (c) as a function of the number of valid records in each time-series. Dashed segments in (b) show the ranges of values considered in Fig. 8. Only TGs belonging to the ALL set () are considered. Circles in (b) and (c) mark the D97 TGs.

Figure 8

Spatial distribution (left-hand side) and histogram of observed rates (right-hand side) for the TGs with a number of valid RLR annual records , 60, 90 and 110 (from top to bottom). For each subset, the right-hand frames also show the basic statistics (μ, rms and wrms). The effect of GIA corrections on these distributions is shown in Fig. 12.

For reference, the numerical values of the individual rates of sea level change ρ_k (given by eq. 4) for the D97 set are listed in Table 2. The locations of these sites are shown in Fig. 2(b). Following Mitrovica et al. (2001), the site of Lyttelton (New Zealand), is not considered because its rate is inconsistent with those at nearby sites. Rates are rounded to the first significant digit since it is our opinion that a more precise evaluation is not possible given the noise level of TG observations. It should be remarked, however, that some of the previous estimates of GMSLR are given with two significant digits, see Table 1. Rates of sea level change for the D97 set, computed since 1880, are all neatly positive and coherently concur to indicate a sea level rise of about 1.5 mm yr⁻¹. Our computations in Table 2, based on eqs (1) and (2), are generally consistent with those of D97. Some of the discrepancies can be attributed to our choice of not pre-processing the observations by a low-pass filter and possibly to different details of the statistical regression adopted to evaluate the trends r_k; some others to the significantly extended RLR record that is currently available for some specific TGs, as shown in Figs 3 and 4. In the case of Balboa (Central America), the period of the currently available RLR observations exceeds by more than 30 yr the one considered by D97, but the two rates do not differ significantly. However, a significant inconsistency with D97 can be observed for the nearby site of Cristobal. Inspection of the annual data for this site, shown in Fig. 4, reveals that the cause is a sea level rise of ∼10 cm that occurred between 1970 and 1972, which has altered significantly the slope of the best fitting line over the whole time period of observation.

In the following, we will often consider trends of sea level change obtained from subsets extracted from the whole set of available RLR time-series, containing a minimum number of yearly observations or obeying more specific criteria. In both circumstances, for a subset composed by a number N_tg of TGs, each providing a rate of sea level change ρ_k = r_k±σ_k, the ‘best estimate’ of GMSLR will be computed by the arithmetic mean

6

where the observed rate r_k is obtained by eq. (1), and ∑_k stands for . Others, like Gröger & Plag (1993) preferred to use the median of the distribution to minimize the possible influence of exceptionally large values of the trends. It is now widely recognized that evaluating the uncertainty associated with g is of utmost importance. However, as it appears from Table 1, in early studies published until the 1980s no uncertainty was generally provided. Subsequently, we will consider three possible error estimates for g. The first is the rms:

7

which characterizes the average uncertainty of the individual trends r_k. The second is the standard deviation of the mean

8

which represents the uncertainty of the best estimate g (see e.g. Taylor 1997). Since in the definition of rms and sdom above the individual uncertainties σ_k do not play any role, it is sometimes convenient (see e.g. Mitrovica et al. 2001) to consider a third estimate of the error to properly account for the large variability of the σ_k values (see Fig. 6). This is defined as the weighted rms

9

where the weights are defined as . To characterize the GMSLR estimates obtained in the following from a given set of TGs, these will be written as

10

and the corresponding numerical values of rms and wrms will be also provided for the sake of completeness.

2.3 GMSLR from uncorrected TG data

In Fig. 7(a), we show the histogram of the distribution of the rates of sea level change r_k for all TGs with a minimum number of yearly records (the ALL set, see Fig. 2a). For this set our GMSLR estimate is μ = 1.4 ± 0.2 mm yr⁻¹, where the small value of sdom compared to m indicates that the latter is, in this case, quite well constrained (the fractional uncertainty is sdom/|g|≈ 15 per cent). According to eq. (8), this does not result from a small rms value, which in fact amounts to 5.8 mm yr⁻¹; rather, it is due to the large number of TGs belonging to the ALL set (N_tg = 1213). In spite of the small sdom, the relatively large rms value clearly confirms that the one obtained using the whole set of available TGs is not a robust GMSLR estimate.

Figure 7

Distributions of the rates of sea level change obtained from eq. (1) for the ALL set of TGs (a), and for the D97 set (b). A vertical dashed segment marks the average value m. Other statistics are given in the insets.

The overall accuracy of the result shown in Fig. 7(a) is difficult to assess, since this would imply the knowledge of the ‘true value’ of g (see e.g. Taylor 1997), which here is not given a priori. It is certain, however, that because of its global character, GIA potentially constitutes a source of systematic errors which can deteriorate the accuracy of the estimate (Nakiboglu & Lambeck 1991; Mitrovica & Davis 1995; Peltier 1996). At individual TGs, these errors cannot be eliminated by increasing the record length, since the rates of sea level change associated with GIA can be safely assumed to be constant on a century timescale. Regarding the propagation of GIA-induced systematic errors into the estimate g, we observe that the GIA contribution to present day sea level change is globally vanishing across the oceans, on the average (see Section 3). However, since the TG distribution is certainly not uniform across the oceans (Fig. 2), we should not expect that the cumulative effect of GIA will vanish leaving g unaffected. For these reasons, the distribution of Fig. 7(a) is likely to be biased by GIA by an unknown amount. While a quantitative assessment will be given starting from Section 3, we note that the uncorrected estimate provided in Fig. 7(a), though possibly inaccurate, is not totally inconsistent with some of the previous, GIA-uncorrected estimates listed in Table 1, which in some cases are based on the thorough selection of a limited number of TGs meeting specific criteria. In Fig. 7(b), we illustrate the particularly important case of D97. The data are taken from Table 2. We observe that the GMSLR estimate obtained in this case (μ = 1.6 ± 0.1 mm yr⁻¹), is consistent with that obtained from the ALL set in Fig. 7(a), where μ = 1.4 ± 0.2 mm yr⁻¹. Hence, in spite of the sparse distribution of the ALL stations, decadal and interdecadal fluctuations and gaps in the time-series and spurious tectonic effects, when taken collectively, the TG observations appear to provide a reliable (though imprecise, because of the relatively large rms) GMSLR estimate.

In Fig. 8, sites with short records are progressively removed from the ALL set considered in Fig. 7(a); the corresponding distribution of rates r_k are shown in the right frames. In particular, in Figs 8(a)–(d), we consider the TG sites with a number of valid yearly observations , for , 60, 90 and 110, respectively. Ideally, assuming negligible systematic errors from, for example, GIA, with increasing we would expect that the computed rates are less and less influenced by the interdecadal components of oceanographic origin that affect individual time-series. Consistent with Fig. 6(a), the number N_tg of sites considered in each of the maps of Fig. 8 is strongly decreasing with increasing , which deteriorates the estimate of μ and the smoothness of the histograms of the r_k values. In this respect, we note that for (Fig. 8b), the histogram is similar to that in fig. 4 of Douglas (1991), corresponding to . A smooth best-fitting normal curve is clearly questionable for such a distribution, contrary to what we have seen in Fig. 7(a) for very large N_tg values (the ALL set). Even a cursory inspection of the histograms reveals that the rms values are virtually unaffected by the reduction of , and the same can be observed for wrms. As a consequence, the sdom value progressively increases until the average sea level μ becomes largely undetermined, and eventually falls to values close to zero. Figs 8(c) and (d) show that negative values of μ, corresponding to an average sea level ‘fall’ cannot be discounted by our analysis for and . This effect is clearly associated with the increasingly dominant role, in the assessment of μ, of strongly negative sea level trends from sites located in areas deeply covered by ice at the LGM. According to the maps of Fig. 8, these sites with long records are predominantly located in northern Europe, and particularly in the Baltic Sea. Gröger & Plag (1993) have reached similar conclusions, observing that the median of sea level trends decreases with increasing , due to the increasing weight of the Scandinavian records.

The dependence of μ and of the respective errors upon and N_tg is better illustrated in Fig. 9. The increase of sdom with increasing and the insensitivity of rms and wrms in a very broad range of values are now apparent. We note that g is monotonically decreasing for , while for larger values it follows a more erratic curve that reflects the progressively increasing spread of the histograms shown in Fig. 8. Also because of the relatively large sdom values for , the estimate μ=g± sdom becomes progressively unstable. However, all the μ values fall approximatively between 0 and 1 mm yr⁻¹. Incidentally, we note that the number of TGs for is N_tg = 23, the same as the specific subset of TGs selected by D97 according to rigorous requirements not uniquely based on . The geographical distribution of these sites, shown in Fig. 2(b), differs significantly from that in Fig. 8(d), at least in two respects. First, the D97 set does not include TGs from formerly ice-covered areas and from the peripheral bulge. Second, it contains a small but significant number of TGs (four) from the Southern Hemisphere, which are totally absent if the selection of sites is merely based on the length of the record, as in Fig. 8. By the histogram of Fig. 7(b), it is apparent that a careful selection of the sites implies a much better constrained value of μ (sdom = 0.1 mm yr⁻¹) and a considerably reduced rms compared to Fig. 8(d). This confirms that the estimate of μ involves a more thorough approach than the simple decimation of the ALL set and that, in general, it cannot be performed by a straightforward automatic selection. As discussed in Section 4, such assessment is prone to imperfect evaluations of the reliability of individual TGs as genuine indicators of climate-driven sea level variations.

Figure 9

Statistics for sets of TG stations containing records with a minimum number of RLR annual records, as a function of . The corresponding value of N_tg is shown in the second x-axis on top. No GIA correction is applied before the statistics are computed.

3 GIA Correction of TG Observations

3.1 GIA modelling

The GIA component of sea level change is evaluated solving the ‘Sea Level Equation’ (hereinafter referred to as SLE). It reads:

11

(Farrell & Clark 1976) where S is sea level change, ω = (θ, λ), where θ is colatitude and λ is longitude, t is time, and S^gi(ω, t), S^e(t) and S^hi(ω, t) are the glacio-isostatic, eustatic and hydro-isostatic components of sea level change (e.g. Spada & Stocchi 2006). All terms of the SLE depend on the history of ice thickness variation. In particular, terms S^gi and S^hi also depend on mantle rheology and are integral functions of S itself, which makes the SLE an implicit equation. The eustatic term reads

12

where m_i(t) is the ice mass variation (m_i(t) =−m_o(t), m_o(t) is the mass variation of the oceans), ρ_w is water density (assumed to be constant) and A_o represents the current area of the surface of the oceans. The eustatic term represents the spatially uniform sea level variation that would be observed for a rigid earth, also neglecting the gravitational interactions between the solid earth, the oceans and the ice sheets. Regardless of the assumptions about the rheological layering of the Earth, an important property of the SLE is

13

where the overline indicates the average over the oceans surface (e.g. Spada & Stocchi 2006).

To address the role of GIA corrections upon the estimate of GMSLR, in this work we will consider three previously published time histories of the late-Pleistocene ice sheets: ICE–3G (Tushingham & Peltier 1991), ICE–5G (Peltier 2004) and the ice model progressively developed at the Research School of Earth Sciences (RSES) of the National Australian University by Kurt Lambeck and coworkers (see Lambeck et al. 1998, and subsequent improvements). Since this latter ice model is valid as of 2005, we will refer to it as to KL05 in the following (Kurt Lambeck, personal communication, 2012). The main features of these models are shown in Fig. 10. Left frames show the time histories of equivalent sea level (ESL) for ICE–3G, ICE–5G and KL05. At a time t before present,

14

where ρ_i and ρ_w are the densities of ice and water, respectively, V_i(t) is the ice volume, V_i(t_p) is present day ice volume and A_o is the current area of the oceans surface. In Fig. 10(a), ESL is shown for the total ice aggregates, while in (b) and (c) we separately consider the major ones: Laurentide and Antarctica, respectively. In all cases, isostatic equilibrium (constant ice mass) is assumed before the LGM and no ice mass variations are considered during the last 1000 yr. For all models, we employ a time discretization with steps of 1 ka. Right frames show the distribution of ice complexes for the three models, at the times when the corresponding time histories initiate. The maps show T(θ, λ, t) −T(θ, λ, t_p), where T is ice thickness. The total ESL for the three models is ‘broadly’ following (to within ∼± 10 m) the same curve until ∼5 ka (a). This curve is constrained by evidence from coral reefs at Barbados, see for example, Fairbanks (1989) and subsequent refinements. However, the ice volumes are clearly partitioned in different manners and partial ESL curves differ considerably. Though barely visible in (c), in KL05 the melting of Antarctica continues essentially until 1000 yr ago (Nakada & Lambeck 1988), while in the other models it terminates between 5 and 4 kyr BP. In addition to the distinct viscosity profiles assumed, summarized in Table 3, these features have important consequences on the assessment of present day regional sea level variations. For the motivations and hypotheses behind the three reconstructions, the readers are referred to the original works quoted earlier.

Figure 10

Left-hand side: equivalent sea level (ESL, see eq. 14) for models ICE–3G, ICE–5G and KL05, relative to the total ice aggregates (a), the Laurentian (b) and Antarctica (c). Right-hand side: spatial distribution of the ice thickness for the three models, at the times when the corresponding time histories initiate.

Table 3

Lithospheric thickness and mantle viscosity for the GIA models employed in this study. ICE–3G data are those presented in Tushingham & Peltier (1991). KL05 data are taken from Lambeck et al. (1998), while ICE–5G values have been obtained by volume-averaging the finely layered, original VM2 viscosity profile of Peltier (2004).

All ice models are based on a respective 1-D viscoelastic stratification of the solid Earth, as shown in Table 3. Possible effects from 3-D lateral variation in mantle viscosity on final predictions, discussed by Kendall et al. (2006) and Spada et al. (2006), here are neglected. Eq. (11) is solved numerically by a pseudo-spectral iterative approach using an improved version of code SELEN (Spada & Stocchi 2007). In our reference runs, we have employed a maximum harmonic degree ℓ_max = 128 on a quasi-regular hicosahedral geodesic grid (Tegmark 1996) with a spacing of ∼45 km, assuming a fixed shorelines and neglecting the sea level variations associated to the Earth's rotational fluctuations (Milne & Mitrovica 1998). The effects of relaxing this latter assumption will be discussed in Section 4. We remark that details of our numerical implementation of the SLE, described by Spada & Stocchi (2007), may differ from those in the original works where the ice models have been presented. Therefore, the results may deviate somewhat from those obtained by other GIA modelers, although the ice models used are nominally the same. A recent GIA benchmark study has shown that there is a broad consensus in the community on the numerical implementation of some of the major components of the SLE (Spada et al. 2011b).

Maps of the GIA component of present day rate of sea level change are shown in Fig. 11. This has been computed by evaluating numerically the time derivative

15

where S represents the solution of the SLE (eq. 11), and t = t_p is present time. Left and right frames show results obtained neglecting and including the effects of the Earth's rotation on sea level variations, respectively (for the latter simulations, we have strictly followed the theory outlined by Milne & Mitrovica 1998). Since for the loads of the size of the major late-Pleistocene ice sheets the Maxwell relaxation time of the mantle is of the order of a few kilo-years (see e.g. Schubert et al. 2004), the rates of sea level change shown can be effectively considered as constant through the period of the instrumental sea level record, at least until the role of transient (i.e. non-Maxwellian) rheological components of deformation will be fully ascertained (Ivins & Sammis 1996; Spada et al. 2011a). Since the three GIA models considered in Fig. 11 assume no ice melting during the last thousand years (i.e. ), according to eq. (12) the ocean-averaged rate of GIA sea level change is

16

regardless of the effective implementation of rotational feedbacks into the SLE, if we do not consider moving coastlines, that is, A_o is a constant in eq. (12).

Figure 11

GIA component of present day rate of sea level change, evaluated according models ICE–3G (a), ICE–5G (c) and KL05 (e). Frames (b), (d) and (f) show results of computations based on the same ice models, but taking into account for the effects of the Earth's rotation on sea level change. The contour corresponding to = 0 is marked by thick curves; positive and negative values of are shown by black and grey thin lines, respectively. Contour lines are also plotted through the continents. TG sites belonging to the D97 are marked by dots.

While in Fig. 11 the global features of are quite similar for the three ice models considered, it is possible to note remarkable regional and local differences (sometimes at the 0.5 mm yr⁻¹ level), reflecting the different chronologies of the ice sheets and the distinct viscosity profiles adopted (see Table 3). Though a precise identification of the causes of the different patterns shown in Fig. 11 is not of concern here, some have a very simple explanation. For instance, it is apparent that the sea level change patterns suggested by ICE–5G (Fig. 11c) and KL05 (e) across north America and the north Atlantic differ significantly from those of ICE–3G (a), and we expect that this may have significant impact on the GIA correction at the numerous nearby TGs. These differences remain in the simulations that include rotational effects (right frames). Inspection of Fig. 10 shows that the cause is the larger mass of the Laurentian ice sheet in ICE–5G and KL05, which implies a more vigorous sea level fall across north America and a much broader peripheral region of sea level rise compared to ICE–3G. Another major difference between the models considered is the total amount of meltwater from Antarctica (see Fig. 10c) and the details of the history of melting in this region. This feature, along with the recognized role played by rotational effects (Milne & Mitrovica 1998), explains the different shape of contour lines in the Southern Hemisphere, and is expected to influence the rates of sea level change at south America TGs.

3.2 GIA corrections of sea level trends

GIA corrections at TGs are directly obtained from the maps of Fig. 11. Following eq. (15), they can be expressed as

17

where ω_k = (θ_k, λ_k) are the coordinates of kth TG and t_p is present time. Hence, the GIA-corrected rates of sea level change are

18

where and are given by eqs (17) and (1), respectively. Assuming ‘exact’ GIA corrections, the uncertainty on r_k will directly propagate on . Hence

19

where σ_k is given by eq. (2). For reference, our ICE–3G corrections are shown in the last column of Table 2 for the D97 set. According to our computations, at these sites their amplitudes range between −0.5 and +0.2 mm yr⁻¹ and show little regional variations as a consequence of the long-wavelength pattern of the GIA-induced sea level variations. Comparison of these estimates with those in table I of D97, corresponding to the same GIA model, shows significant differences. This confirms that the details of the numerical implementation of the SLE can play an important role in the assessment of GMSLR, and suggests that the assumption of fully reliable GIA corrections is questionable until a broad consensus within the community is met (Spada et al. 2011b). Statistics derived by corrected rates will be indicated with primed symbols (e.g. g′, wrms′, …). These are obtained by substituting r_k with in eqs (6)–(9). In particular, we note that

20

where

21

simply represents the average of the GIA corrections at TGs. According to eqs (16) and (17), this term would vanish, hence providing g′=g, only in the very idealized case of a regularly spaced, dense network of TGs located across the oceans. Obviously, in view of the spatially inhomogeneous distribution of TGs in the real world (see e.g. Fig. 2), we should not expect that g^gia effectively vanish. This is well illustrated in Table 2, showing that the GIA corrections to TG trends performed by means of model ICE–3G have, with a few exceptions, a negative sign.

With the aid of eq. (5), it is possible to discuss how much the GIA-corrected GMSLR g′ contributes to assess a volume change of global oceans, and to show the relevance of this quantity. Eq. (20) gives

22

where the terms on right-hand side, associated with the secular melting of continental ice masses, the steric effect, tectonic motions and other possible factors, have the same form of eq. (21), with r^gia substituted by , , and , respectively. A relationship between g′ and the rate of volume change of the global oceans can be established relying upon a series of assumptions. The first is that the GIA correction g^gia is computed correctly and is substantially independent from the GIA model adopted (this will be accomplished in Section 4 by a proper selection of TGs). Second, in the idealized case of a dense and regular network of TGs, or considering sufficiently long observations that can be assumed to be representative of a large-scale sea level change (see Douglas 2001), according to eqs (12) and (13) we have , where is the rate of change of the mass of the continental ice sheets. By mass conservation = , where is the rate of mass variation of the oceans in response to the melting of continental ice masses and is the corresponding rate of volume change. This gives

23

which states proportionality between g^cma and ocean volume variations associated with the mass contribution to secular sea level rise. Third, assuming a spatially uniform steric sea level rise on the secular timescale, we similarly obtain

24

where represents the rate of change of the steric component of the ocean volume variation. Fourth, in the hypothesis of a substantial tectonic stability of the sites where TGs are located, or vanishingly small globally averaged tectonic movements, and also neglecting other unknown causes of sea level rise, by substitution of eqs (23) and (24) into (22), for the GIA-corrected GMSLR we obtain

25

where is the rate of variation of the total volume of the oceans. Of course, because of the assumptions mentioned above, eq. (25) should be viewed as an idealized approximation, which however represents the only means for estimating secular variations of the ocean's volume from the available instrumental record.

In Fig. 12 we show the distributions of the GIA-corrected rates for various subsets of TGs. Here we have used GIA predictions by model ICE–3G, since it has been widely employed in the past (see Table 1); corrections for GIA models ICE–5G and KL05 will be discussed in Section 4. In Figs 12(a)–(e), the minimum number of valid yearly observations is increased from (the ALL set) to , while in Fig. 12(f), the D97 set is considered. Solid and stairs-step histograms, characterized by the same area, show the distribution of and of , respectively (these latter are reproduced from Figs 7 and 8). Vertical dashed segments show the GIA-corrected average g′ (eq. 20); other primed statistics obtained from eqs (7) and (9) are also shown. The effect of GIA corrections, which can be appreciated even visually, is that of reducing the spread of the distributions around their average value, or increasing the overall coherence of the results. We note, however, that looking for the geographic coherency of the GIA-corrected results may be counter productive, since here we are not considering the sea level variability associated with the melting of polar ice sheets and glaciers, as well as steric effects (Mitrovica et al. 2001). The weighted variance reduction, estimated by a Fisher F-test for the ratio wrms′/wrms, is significant at the 95 per cent confidence level for all the sets considered from Fig. 12(a)–(e). For the D97 set in Fig. 12(f), the variance reduction obtained by applying the GIA correction is not significant (at the same confidence level), and the average value of sea level change μ′ = 1.7 ± 0.1 mm yr⁻¹ overlaps the uncorrected estimate (μ = 1.6 ± 0.1, see Fig. 7b and Table 7). However, from the statistics in Fig. 12(f), it appears that the five requirements introduced by D97 produce a sizeable variance reduction with respect to the case when the only criterion is minimizing the number of yearly records (, see Fig. 12e).

Figure 12

Frames from (a) to (e) show the effects of the GIA correction (by ICE–3G) on the distributions of the sea level trends for increasing values of (corrected and uncorrected distributions are shown by lighter grey (green in online article) solid and stairs-step histograms, respectively). The D97 set, shown in darker grey (magenta in online article) is considered in frame (f). The statistics are summarized in (g) for more values of .

Table 7

Estimates of GMSLR (μ=g± sdom) obtained in this work, corresponding to different subsets of TGs and GIA corrections evaluated using code SELEN (Spada & Stocchi 2007). These trends have not been corrected for any tectonic effect. The statistics rms and wrms are also shown. Our preferred estimate is starred.

The statistics for the corrected rates, shown in frame Fig. 12(g) as a function of , confirm that with increasing the GIA correction plays an increasing role. Furthermore, they differ from the corresponding uncorrected statistics, shown in Fig. 9, in several respects. First, g′ shows a reduced sensitivity to compared to g. Second, the standard error of the mean (sdom) is smaller, especially for the longest records, as an effect of the reduced rms. Third, g′ and wrms′ broadly follow the same trend with varying . According to eq. (9), this indicates that for the time-series with the longest record (hence, from Fig. 12, those having generally the largest weight w_k in the expression of wrms′), the GIA correction removes a large fraction of the observed trend, that is, . As pointed by D97, since these series are from the regions deeply covered by ice at the LGM, even small errors in the GIA predictions for these regions can be comparable to the expected value of g′. This justifies, in D97, the exclusion of these sites from the analysis by the application of criterion (5).

3.3 Uncertainties in GIA modelling

A major drawback of the approach followed so far is the assumption that the GIA correction is perfectly determined. This is not the case, however. In fact, there are at least three sources of error which can affect . The first is associated with the physical ingredients implemented in the SLE (eq. 11). For instance, the sea level variations associated with fluctuations of the Earth's rotation vector (Milne & Mitrovica 1998) can be included or not in the SLE (indeed, this can be accomplished in various ways, as discussed by Spada et al. 2011b). Similarly, to minimize the computational burden, the SLE can be solved assuming fixed shorelines or, alternatively, allowing for the development of a ‘gravitationally self-consistent’ palaeotopography (Peltier 2004). Accounting for or neglecting these effects can certainly have some impact on the GIA correction, both on a global and a regional scale. Once the physics of the SLE is set, a second source of uncertainty is associated with its numerical implementation. For instance, different temporal discretization schemes (the SLE is normally solved assuming stepwise or, alternatively, a linear piecewise discretization), as well as different geometries and spatial resolutions of the spatial grids are factors that can potentially impact the numerical results. Third, there are uncertainties associated with the spatio-temporal features of the late-Pleistocene ice sheets, which derive from errors in the Relative Sea Level (RSL) data that constrain the ice models itself (see the review of Whitehouse 2009). In this respect, since the amount and the overall quality of the RSL observations is increasing with time, GIA models are continuously improved. Consequently, GIA predictions are not definitive. A rigorous evaluation of all the possible sources of errors described earlier is certainly not feasible in this work, since it would imply running a very large number of sensitivity experiments. Hence, we use our code SELEN to estimate the GIA uncertainties adopting the three ice models described in Section 3 and the viscosity profiles recommended by their authors, shown in Table 3. This is done assuming that they provide three reliable (though not completely error-free) descriptions of the sea level variations associated with the GIA phenomenon.

To illustrate the uncertainties inherent to the GIA correction, in Fig. 13 we compare the rates of sea level change at the D97 TGs (a) with the GIA corrections based on models ICE–3G, ICE–5G and KL05 (b). The three GIA models show comparable corrections (to within ±0.3 mm yr⁻¹) at most of the sites, but the ICE–5G and the KL05 predictions strongly disagree with those of ICE–3G along the North American West coast (sites 13–16) and in South East North America (21–23). This mismatch can be, at least partly, attributed to the larger lateral extent and thickness of the Laurentian ice sheet in models ICE–5G and KL05 compared to ICE–3G, as illustrated by the maps and ESL curves in Figs 10(a)–(c) for the three models. To ascertain possible effects of the Earth's rheological layering, keeping fixed the time history of loading to that of ICE–3G, we have varied the upper mantle (see Fig. 13c) and the lower mantle viscosity values (d; values are given in the caption). The results indicate that the GIA corrections for sites along the North American West coast (sites 13–16) and South East North America (21–23) are less sensitive to variations in the upper-mantle viscosity than to variations in the lower-mantle viscosity. This result is found to be in agreement with a number of previous studies in the context of TG records (see, in particular, Davis & Mitrovica 1996).

Figure 13

(a) Rates of sea level change and their uncertainties at the 23 sites of the D97 set (the shaded region shows μ=m± sdom = 1.6 ± 0.1 mm yr⁻¹, which represents our GIA-uncorrected estimate for the sites in D97). The numerical values of the rates are listed in Table 2. (b) GIA corrections for the rates shown in (a), computed for the GIA models ICE–3G, ICE–5G and KL05, respectively. In (c), ICE–3G is used, with upper-mantle viscosity having the values of 0.3 × 10²¹ Pa s, see curve ICE–3G(UM03) and 0.5 × 10²¹ Pa s for ICE–3G(UM05). In (d), the lower-mantle viscosity is varied, with values of 5 × 10²¹ Pa s for ICE–3G(LM5) and 10 × 10²¹ Pa s for ICE–3G(LM10).

The findings illustrated in Fig. 13 have important consequences on the assessment of the GMSLR. In fact, the lateral extent of the peripheral forebulge of ICE–3G is significantly smaller than in ICE–5G and KL05 (see Fig. 11), and TGs located in the North American West coast and South East North America can be safely considered outside the collapsing region. However, this is not the case when ICE–5G and KL05 are considered. This means that, by a rigorous application of the D97 criteria [in particular, of criterion (5), see Table 4], these sites should be excluded from the analysis, unless models ICE–5G and KL05 are considered fully unreliable compared to the previously released ICE–3G. In the following, we denote by D97R the set of 16 TGs obtained excluding these sites (these are marked by a star in Table 2). We anticipate that the GIA-corrected GMSLR estimate based on D97R is, contrary to D97, basically insensitive to the model adopted for the correction. This result suggests a modification of the D97 criteria for the selection of the TGs suitable for the assessment of GMSLR, which is discussed later.

4 Discussion

4.1 A new, GIA-based selection criterion: the SGX set

In Fig. 14(a), we show the rates of sea level change for 44 TGs for which at least yr of valid RLR annual records are available. In the following, this requirement will be referred to as criterion (I) (the new or revised criteria introduced in this study are summarized in Table 4 and compared with those by D97). Criterion (I) constitutes a stronger constraint compared to its counterpart (1) of D97, which is based on the record length (and not on the number of valid years) of possibly non-RLR records. The sites of Fig. 14(a) also satisfy another constraint, namely that the GIA corrections, shown in Fig. 14(b), are nearly the same for all stations, regardless of the model employed. This has been accomplished applying the following selection rule to the ALL TGs:

26

where M₁ and M₂ are any two of the three GIA models considered in this work, and Δr is a small tolerance parameter. Hereinafter, we will use Δr = 0.3 mm yr⁻¹, which roughly represents 20 per cent of the average magnitude of GMSLR (≈1.5 mm yr⁻¹), based on previous assessments (see Table 1). The resulting set, presented in Table 5, will be referred to as ‘SGX set’ in the following. Note that the constraint expressed by eq. (26) is applied to GIA simulations that do not include the effects from the rotational feedback. This is motivated, a posteriori, by the negligible role played by this modelling feature on the assessment of the GMSLR, as it will be discussed later. Since constraint (26) will substitute criterion (5) of D97 (see Table 4), it will be referred to as criterion (V). In the following, criterion (V) is adopted using the three GIA models described in Fig. 10, which have had a particularly important role in past GMSLR determinations (see Table 1). However it can be applied to future GIA models constrained by a specific set of Holocene RSL observations.

Figure 14

(a) Observed rates of sea level change at the 44 TGs belonging to the SGX set. Black symbols denote sites of the SG01 subset, obtained by decimation of SGX according to the discussion in Section 4. In (b) the GIA corrections are shown for all sites. The spatial distribution of TGs is shown in (c), where black dots denote the location of the SG01 TG sites.

Table 5

TGs belonging to the SGX set. The time-series from this set are characterized by at least 60 valid yearly records (criterion I) and obey the constraint expressed by eq. (26; criterion V). The rates with error bars are displayed in Fig. 14.

Imposing simultaneously constraints (I) and (V), our aim is to discard TGs with too short records and those at which different parametrizations of models ICE–3G, ICE–5G and KL05 produce conflicting GIA corrections. The SGX sites so selected should be only sensitive, in principle, to the gross features of the spatial pattern of GIA-induced sea level change, shared by any realistic model. An obvious hypothesis is that this feature is the amplitude of the S^hi term of the SLE (see eq. 11), which is mainly determined by the time history of the meltwater loading and follows broadly similar curves for the three models considered (see the ESL curves in Fig. 10a). This is confirmed by the geographical distribution of the TG sites belonging to SGX (see Fig. 14c), showing that the constraints (I) and (V) are very effective (with a few exceptions, discussed later) in eliminating the TGs located in the vicinity of the major formerly ice-covered areas and those across the collapsing forebulges, where the three GIA models provide broadly different predictions, as a consequence of the different time histories of ice melting at a regional scale. These are expected to be mostly affected by the S^gi term of the SLE and, only to a lesser extent, by hydro-isostatic effects. Since it is unlikely that the future GIA models will adopt global ESL curves that dramatically diverge from those displayed in Fig. 10(a), we believe that SGX constitutes the most natural set of TGs on which to rely upon for GMSLR assessments.

4.2 Decimating SGX: the SG01 set

The anomalous rates of sea level change observed for a significant number of the SGX sites (in some cases as large as several millimetres per year, see Fig. 14a), demands the application of selection criteria in addition to (I) and (V), introduced in Table 4. They are: (II) the TG station must be sufficiently distant from collisional tectonic boundaries, (III) they must have a sufficient completeness (>70 per cent, where completeness is the ratio where span is the difference between the newest and the oldest annual RLR record), and (IV) the low-frequency sea level trend from each TG must reasonably agree with that of nearby instruments. Criteria (II) and (IV) are directly borrowed from (2) and (4) of D97 (see Table 4). We note that sites showing ‘suspect accelerations’ or those where human-driven effects are ascertained, will be prudently avoided to prevent data contamination also following the suggestion of Pirazzoli (1993). This criterion, which has been often implicitly or explicitly applied in previous studies (see e.g. Hagedoorn et al. 2007), will be referred to as criterion (VI). In what follows, starting from set SGX, criteria (II)–(IV) and (VI) are employed to define a new of preferred TGs (set SG01), which will be adopted in Section 4.3 to produce new GMSLR estimates.

4.2.1 Tectonic contamination

Based on criterion (II), we have first expunged from the SGX set all TGs located close to active collisional tectonic boundaries. These include, in particular, the stations along the coast of Japan, whose low-frequency signals (periods longer than 50 yr) have long been recognized to be associated to the subduction of the Pacific and Philippine plates (Aubrey & Emery 1986). According to the current RLR PSMSL record, at the six Japanese SGX sites of Mera, Aburatsubo, Uchiura, Hosojuma, Tonoura and Wajima, the long-term rates of sea level change show a considerable spread (with values ranging between ∼−0.5 to ∼+4 mm yr⁻¹), which prohibits the assessment of a reliable regional average value. We note that other sites that could be potentially affected by tectonic deformations along the San Andreas Fault, like those along the North American West coast (San Francisco, Santa Monica, La Jolla and San Diego), which are included in the D97 set (see Table 2), have been automatically removed by the application of constraint (V).

4.2.2 Human-driven effects and anomalous accelerations

Applying criterion (VI), several sites have been removed from the SGX set for possible contaminations from human-driven subsidence, as in the case of Venice (see e.g. Hagedoorn et al. 2007; Carbognin et al. 2010), whose sea level trend exceeds significantly those at other Mediterranean TGs located at approximately the same latitude (e.g. Trieste). For some other SGX sites, suspect accelerations are apparent in the time-series. This is particularly evident for the sites of Batumi and Poti (Black Sea), Takoradi (Ghana), Mumbai and Calcutta (India), Ko Lak, Ko Taphao Noi and Fort Phrachula Chomklao (Thailand), Manila S. Harbor (Philippines) and Hilo (Hawaii Island). Possible origins of these anomalous trends, discussed in the geophysical signals page of PSMSL http://www.psmsl.org/train_and_info/geo_signals/, include increased groundwater extraction, recent deposit from river discharges, land reclamation works, various coastal processes and possibly volcanic deformations, as for example in Hilo (Hawaii).

Anomalous sea level trends at the sites above have been reported and discussed in a number of previous studies. The subsidence at Batumi and Poti is confirmed by Garcia et al. (2007) and by Stanev & Peneva (2002). Caccamise et al. (2005) have described an inconsistency of the Hilo TG record compared to the nearby Honolulu TG, whereas Yanagi & Akaki (1994) and Emery & Aubrey (1991) have explained the anomalous increase in sea level from 1965 to 1982 at Manila and Fort Phrachula Chomklao invoking a large withdrawal of groundwater; the anomalies shown by the Manila record are attributed by Douglas (1991) to the tectonic movements in that area and to harbour development. For the Indian TGs, in addition to the subsidence problems (Subrahmanya 1996), it is recognized that the records were influenced by the high frequency of tropical cyclones and storm surges (Das & Radhakrishna 1991). Finally, a large interannual variability is the cause of the unreliability of the trend observed at the Ko Taphao Noi TG (Unnikrishnan & Shankar 2007). Following criterion (VI), these sites have been prudently eliminated from set SGX.

4.2.3 Eurasian sites

The SGX set contains the Northern UK sites of Aberdeen and North Shields, which were discarded by D97 because of their proximity to the peripheral bulge adjacent of the former ice sheets. However, it appears that the observed rate of sea level change at these sites is not anomalously large to motivate their exclusion (Table 5). Furthermore, from our GIA computations (see Table 6), the GIA rates at these sites have a negative sign (contrary to what it is expected for sites located across the lateral forebulge region) and comparable with GIA rates of far field sites, such as Buenos Aires. Thus, also following Douglas (1991), these Northern UK sites are not discarded (though this is done for Aberdeen I, whose trend of sea level change duplicates that of the nearby location of Aberdeen II).

Table 6

Basic data and computed rates of sea level change for the SG01 set, determined according to the discussion in Section 4. The average number of valid yearly RLR records for SG01 is 87 yr. GIA corrections corresponding to models ICE–5G are also shown. According to eq. (26), those pertaining to ICE–3G and KL05 do not differ from those shown by more than 0.3 mm yr⁻¹. The last two columns show ICE–5G GIA corrections that include the effects of the rotational feedback and that of the horizontal migration of shorelines, respectively.

By similar arguments, the two North European sites of Smogen (Sweden) and Heimjso (Norway), approximately located along the margins of the former Fennoscandian ice sheet, have not been excluded by our analysis. Among the SGX sites, these are the only two that are not located in the far field of the former ice sheets, and are experiencing secular rates of sea level change of −1.5 ± 0.2 and −1.9 ± 0.1 mm yr⁻¹, respectively (see Table 5). At these sites, the GIA models considered in this study coherently indicate a sea level fall of ≈2.5 mm yr⁻¹ (see Fig. 14). We note that Smogen was indeed included in the analysis of Douglas (1991) but it was later discarded by D97 applying criterion (5). Hemisjo was not considered by Douglas (1991) probably because, at that time, the TG record was not long enough to merit attention (the Smogen TG record starts only in 1928). Presently, Smogen and Heimjso meet the requirements of minimum length and of completeness (77 per cent for Heimjso and 99 per cent for Smogen). We observe that the strongly negative values of the rate of sea level change observed for these TGs, if corrected for the GIA effect becomes largely coherent with the values found for other TGs.

Finally, contrary to D97, we did not have any reason to exclude the site of Tiksi (Tiksi Buktha, Siberia), which is remarkably complete (100 per cent). Since the Tiksi record has been employed in previous studies aimed to a regional assessment of sea level trend (see Pavlov 2001; Proshutinsky et al. 2004, 2007), it will not be expunged by the SGX set.

4.2.4 Australian sites

The SGX set contains five Australian TGs (namely, Sydney, Fort Denison, Fort Denison II, Newcastle and Fremantle). Since the two Sydney stations provide inconsistent trends (see Table 5) and it is impossible to determine which one is effectively correct, both have been excluded from our assessment, also in agreement with criterion (IV). In addition, we also note that the observed sea level rate from the nearby site of Newcastle is largely inconsistent with the two Sydney trends [as pointed by Douglas (1991), the major discrepancy between Newcastle and Sydney is observed during the period 1950–1980]. Consequently, we have also expunged the site of Newcastle from set SGX. However, in view of its considerable distance from Sydney and Newcastle and the fairly long and complete record (see Table 5), we did not discard the site of Fremantle. As pointed by Church et al. (2004), the variability related with the El Niño–Southern Oscillation (ENSO) does not appear to deteriorate the assessment of the secular sea level trend obtained from this site. Previous authors have assumed distinct viewpoints regarding the Australian TG sites. For instance, D97 excluded Sydney and Newcastle and did not consider the TG of Fremantle, apparently without any explicit motivation. We also note that more recently, Hagedoorn et al. (2007) considered various Australian sites in their study, including Sydney and Fremantle but excluding Newcastle.

4.3 New GMSLR estimates

The basic data pertaining to the SG01 sites, surviving the decimation of set SGX, are displayed in Table 6. The rates of sea level change for these sites are marked in Fig. 14(a) by black symbols, while those discarded from set SGX are shown by grey symbols. The corresponding GIA corrections are shown in Fig. 14(b) for the three models considered in this study. As shown in Fig. 14(c) by filled circles, the spatial coverage of set SG01 appears not too severely biased towards the Northern Hemisphere. As discussed in Section 3.2, a sufficiently uniform geographical distribution is one of the conditions that would make the GIA-corrected GMSLR a useful indicator of the volume change of the global oceans.

The various GMSLR estimates obtained in this work are summarized in Table 7. As anticipated in Section 3.3, for set D97R, the three GIA-corrected GMSLR values point to a value μ′ = 1.5 ± 0.1 mm yr⁻¹ and are broadly consistent one to each other (to within ±0.1 mm yr⁻¹), regardless of the model employed to compute the GIA corrections at individual TGs. This is evidently in contrast with the estimates based on D97, which show a significant sensitivity to the model employed, with μ′ in the range between 1.3 ± 0.1 and 1.7 ± 0.1 mm yr⁻¹. As discussed in Section 3.3 and Fig. 6, the reason for this sensitivity resides in the largely different predictions provided by the GIA models at TGs located along the North American West coast and in South East North America. Another appealing feature of D97R is the relatively small difference between the estimates with and without a GIA correction (0.1 mm yr⁻¹). This indicates that though the D97R TGs are sparsely distributed, the average GIA corrections for these sites is close to zero, as it would be for a set of TGs uniformly spaced across the oceans (see eq. 20). We also observe, for D97R, a systematic reduction of rms′ and wrms′ with respect to D97, which indicates that the GMSLR estimate obtained from the former is slightly more precise than the one based upon D97.

As an alternative to D97 and D97R, in Section 4.2 and in the subsequent discussion we have defined set SG01 (the 22 SG01 TGs, with individual rates and basic data, are listed in Table 6). According to eq. (26), TGs belonging to SG01 are characterized by a GIA correction that is essentially independent from the particular model employed to determine it. The application of this constraint makes the process of TG selection less prone to contamination from possible errors in GIA models and, consequently, the assessment of GMSLR (and of the total variation of the global oceans volume) more robust overall. The SG01 results of Table 7 (see estimates 11–14) reveal that for this set of TGs the application of the GIA correction implies a sizable modification of μ, contrary to what we observe for D97. The reason is the presence, in SG01, of the two sites of Heimsjo and Smogen, located close to the former Fennoscandian ice sheet (see Fig. 14), and characterized by negative sea level trends. Consistent with constraint of eq. (26), after GIA correction the GMSLR obtained using the SG01 TGs is largely independent on the model adopted for compute the correction, and point to a value of μ′ = 1.5 ± 0.1 (rms = 0.4 mm yr⁻¹, wrms = 0.3 mm yr⁻¹). This value represents our ‘preferred’ GMSLR estimate so far. By further computations, we have verified that excluding from SG01 the two sites of Heimsjo and Smogen, which clearly would appear as outliers with respect to the ρ_k values listed in Table 6, would not change the value of μ′ or that of rms and wrms at the 0.1 mm yr⁻¹ level. In this case, similarly to D97R, the condition μ′≈μ would be attained.

In all the GMSLR estimates of Table 7 discussed so far (n. 1–15), in our GIA modelling we have neglected the effects from the rotational feedback on sea level variations (Milne & Mitrovica 1998). Even a cursory inspection of Fig. 11 reveals that polar motion produces a long-wavelength sea level pattern of a few fractions of millimetres per year, particularly visible in the Southern Hemisphere, where it is not overwhelmed by the sea level changes directly associated with the melting of the major Northern Hemisphere ice sheets. This pattern is dominated by a degree 2 order 1 spherical harmonic term, with is associated with the variations of the Earth's centrifugal potential (Munk & MacDonald 1960). A site-by-site analysis (see last columns of Table 6), shows that that present day rates of sea level change at the SG01 TGs are perturbed by as much as 0.2 mm yr⁻¹, well above the typical sdom in Table 7. This fully agrees with the estimates by Milne & Mitrovica (1998), although these were totally based on rates of sea level change extrapolated from RSL curves in response to the melting of late-Pleistocene ice sheets. Of course, the significance of individual GIA corrections does not necessarily imply that the Earth's rotational response significantly impacts the GMSLR estimate, since according to eq. (20) the value of g′ depends on the TG-averaged values. When this average is computed, all the GIA-corrected GMSLR estimates indicate the same μ′ value (see estimate n. 15 in Table 7), which does not depart (at the 0.1 mm yr⁻¹ level), with our ‘preferred estimate’μ′ = 1.5 ± 0.1 mm yr⁻¹ (rms = 0.4 mm yr⁻¹, wrms = 0.3 mm yr⁻¹), obtained by neglecting the rotational effects on sea level change. As a further test on the effects of approximations in GIA modelling, we have abandoned the fixed-shorelines approximation. In particular, we have implemented, for model ICE–5G, the iterative method outlined by Peltier (2004), which allows for a reconstruction of the palaeo-topography since the LGM. However, as shown by estimate n. 16 of Table 7, adding the realistic feature of migrating shorelines does not modify, at the 0.1 mm yr⁻¹ level, our preferred value of GMSLR (individual GIA corrections are displayed in the last column of Table 6). It should be observed, however, that our numerical code does not implement the more refined Sea Level Equation theory proposed by Mitrovica (2003) and Mitrovica & Milne (2003) and therefore this last result should be considered as preliminary and possibly not fully accurate.

4.4 Role of tectonic deformations

By application of criterion (II), the SG01 set should not contain TGs that are ‘too close’ to active plate boundaries. However, tectonic vertical deformations can never be totally discounted (see http://www.psmsl.org/train_and_info/geo_signals/). Their possible effects on our GMSLR estimate are difficult to assess, since no global predictive model for tectonic deformations can be invoked on a secular timescale, as it is the case of GIA. Previous attempts to estimate the contribution of coseismic and post-seismic vertical deformation to sea level change at TGs (Melini et al. 2004; Melini & Piersanti 2006) using data from the Seismic Moment Tensor catalogue (see e.g. Dziewonski et al. 1981) have clearly shown that earthquakes have cumulatively the tendency to produce a positive sea level trend at locations of TGs. This is manifest globally, but also at the 23 TGs that are included in the D97 set. According to these studies (see, in particular, Melini et al. 2004), estimates of secular sea level variations that do not account for the effects of global seismicity are likely to be overestimated by ∼0.2 mm yr⁻¹. This will not change dramatically our estimates in Table 7, but nevertheless such effects exceed the sdom associated with the secular sea level trend and therefore call for further investigations, also in view of some recent advancements in the theoretical framework for such calculations (Melini et al. 2009).

Coseismic and post-seismic effects, however, are only one component of the total tectonic deformation field. While in a non-expanding Earth scenario it is reasonable to assume that tectonic vertical movements average out globally, it is clear that their cumulative effect upon μ depends upon the geographical distribution of TGs, which is far from being uniform. Though a direct modelization is not possible, by simple arguments we can estimate the potential effect that regional tectonic deformations could have in the GMSLR estimates of Table 7. For instance, let us assume that the four Mediterranean TGs belonging to the SG01 set (Marseille, Genova, Trieste and Bakar, see Table 6) are coherently affected by a rate of secular sea level change of tectonic origin of constant amplitude r^tec (see eq. 5), while it is assumed that random deformations acting on the remaining SG01 sites have a cumulatively negligible effect. Using eq. (6), it is easy to see that to affect the value of g by an amount ≥2 · sdom, the condition r^tect≥ 1 mm yr⁻¹ should be approximately met if one assumes sdom ≈0.1 mm yr⁻¹. This appears to be unrealistic as an order of magnitude of the typical tectonic trend in the Northern Mediterranean area (Tsimplis et al. 2011).

Though not supported by rigorous quantitative arguments, the above suggests that estimates of μ obtained by the D97, D97R and the SG01 sets are reasonably insensitive to tectonic effects, although individual sites can be certainly affected by local deformations associated with the global seismic activity (Melini et al. 2004).

4.5 Sea level accelerations and recent sea level rise

Up to now, the SG01 set has been solely employed to estimate the secular GMSLR. The question arises, however, as to whether this set of TGs is useful for the assessment of (i) possible sea level secular accelerations and (ii) the sea level variations during the last few decades. To provide a qualitative answer, we have first performed a quadratic regression on the stacked SG01 time-series (see grey curve in Fig. 5a). This has provided a secular sea level trend of ∼1.2 mm yr⁻¹ and no evidence of secular acceleration (the quadratic term of the regression is close to 10⁻⁵ mm yr⁻²). Second, we have carried out a linear regression of the SG01 time-series using the same methods illustrated in Section 2.2. The analysis has been performed on two time periods, both encompassing ∼50 per cent of the total lapse of time initially considered in our analysis, namely 1900–1955 and 1956–2010, respectively. With this choice, the length of the time-series (∼60 yr) is expected to ensure a meaningful evaluation of trends (see e.g. D97). For each of the two time periods, there are 13 time-series with completeness >75 per cent. The values of GMSLR so determined are μ = 1.2 ± 0.3 and μ = 1.1 ± 0.3 mm yr⁻¹, for the two time periods, respectively. Their consistency provides no indication of a sea level acceleration in the middle of last century. By further tests, we have verified that this result is robust with respect to variations of ±20 yr in the location of the midpoint (we note, however, that the number of time-series with completeness >75 per cent is, for these simulations, <10). The results of this simple analysis are not modified by the application of the GIA correction (by model ICE–5G) to the time-series considered, which in the two periods provides μ′ = 1.6 ± 0.1 and 1.6 ± 0.2 mm yr⁻¹, respectively.

Our findings, based on elementary methods, contrast with those by Church & White (2006), who analysed 1658 metric and RLR TG time-series. Based on the time integration of the first differences of the EOF amplitudes, they detected a significant slope variation in ∼1930, and particularly during the periods 1870–1935 and 1936–2001, with a sea level acceleration of +0.017 ± 0.007 mm yr⁻². The question of an acceleration (or deceleration) after ∼1930 is a matter of debate. While Woodworth et al. (2009) argue that an acceleration occurred from around 1870 to the end of 20th century, Houston & Dean (2011) suggest a ‘deceleration’ of −0.0014 ± 0.0161 mm yr⁻² beginning since ∼1930. However, this has been recently criticized by Rahmstorf & Vermeer (2011), who found a ‘acceleration’ since 1930, strongly correlated with global temperature. It is possible that the lack of acceleration suggested by our analysis is mainly due to the limited number (22) of TGs employed. As pointed by Bindoff et al. (2007), sparse TG networks (such as those constituting our SG01 set) are inadequate to detect possible sea level accelerations. As a matter of fact, with small sets, it is likely that the geographical variations of sea level change are not sufficiently sampled (Gregory et al. 2001), and possible accelerations unresolved.

Finally, using the time-series in the SG01 set, we have estimated the GMSLR for the period 1993–2009 (Church & White 2011). Starting from 1993, altimeter satellite data became available and a comparison with results based on the TG record is now possible. The sea level trend for the recent period 1993–2009 is broadly recognized to exceed the average secular value. Our ICE–5G-corrected result μ′ = 3.6 ± 0.5 mm yr⁻¹, which is based on the 12 TGs with completeness >75 per cent during this time period, is found to be consistent with the independent evaluation by Church & White (2011) for the same period, based on altimetry data (μ′ = 3.2 ± 0.4 mm yr⁻¹) and is comparable with that assessed by the same authors using only TG observations (μ′ = 2.8 ± 0.8 mm yr⁻¹). The not exact coincidence of these estimates is likely to reflect the different methods of analysis adopted in these studies. According to Bindoff et al. (2007) these large deviations from the secular trend may partly reflect decadal variability rather than a genuine acceleration, though starting from 1993 thermal expansion has contributed significantly more than during the whole 20th century and melting of continental ice has been responsible for a sea level rise of 1.2 ± 0.4 mm yr⁻¹.

5 Conclusions

After a revisitation of previous work and an overview of the currently available data, we have proposed a new solution to the long-standing problem of the assessment of GMSLR from secular TG observations. Although local sea level variations are of much larger importance for the respective community or society, the global changes detected from properly selected TGs with long records of observations can provide, in principle, information on the ocean volume variation in response to current climate change. GIA corrections to the sea level trends obtained from specific sets of TGs have long been recognized to have an important role in the assessment of GMSLR (see the review of Table 1). Up to now, these have been generally computed a posteriori, once a suitable set of TG has been identified according to specific criteria. In this work, GIA corrections have been employed in a non-traditional way, since their amplitude is quantitatively adopted a priori as one of the selection criteria imposed to constrain the TG set. In particular, we have proposed to substitute criterion (5) suggested by D97 (i.e. exclusion of the TGs belonging to regions which were deeply covered by ice at the LGM or to their surroundings), with a new requirement, that is, that GIA corrections to TG series should be largely independent from the particular GIA model employed to compute them. Our motivation essentially resides in the markedly different ice distributions at the LGM shown by current GIA models, which is the cause of distinct subsidence patterns, especially across the lateral forebulge regions.

In our study, the selection of TGs is based on numerical results obtained from three specific GIA models widely employed in the literature and characterized by distinct viscosity profiles and ice sheets chronologies (namely, ICE–3G, ICE–5G and KL05). However, the approach outlined is easily reproducible and could be directly applied to any other plausible GIA model, provided that it is consistent with a specific set of RSL observations since the LGM. Using the selection criteria (I)–(VI) summarized in Table 4, we have identified a restricted set of 22 TGs (SG01) which are useful for the assessments of the GMSLR. As the result of the application of criterion (V), set SG01 is mostly composed of TGs in the far field of previously ice-covered regions, where GIA predictions are essentially only sensitive to the history of meltwater loading. SG01 significantly differs from (but partly overlaps) the traditional set of 23 TGs proposed by D97. The imposed constraints make our ‘preferred’ GMSLR estimate, namely μ′ = 1.5 ± 0.1 mm yr⁻¹ (rms = 0.4 mm yr⁻¹, wrms = 0.3 mm yr⁻¹) independent on the GIA model employed to compute the corrections. Furthermore, we have found that it is stable (at the 0.1 mm yr⁻¹ level) with respect to the introduction of refined features in the numerical solution of the SLE such as the rotational component of the sea level variations and the horizontal migration of shorelines.

Our GMSLR estimate is consistent with other GIA-corrected estimates obtained from TG observations during comparable timescales (∼ one century), but based on different methods of analysis. Among those listed in Table 1, we find a consistency (i.e. the ranges of estimates are overlapping) with previous results by Mitrovica & Davis (1995), Church et al. (2001), Nakada & Inoue (2005), Church & White (2006) and Hagedoorn et al. (2007). Though the statistical meaning of the error bar on the GMSLR estimate by Bindoff et al. (2007) is not explicitly stated, our estimate is clearly consistent with the results of the last IPCC report. We remark that consistency with the results of D97 is only met for the GIA-uncorrected estimate. The corrected one ( mm yr⁻¹) is inconsistent with ours (and exceeds it), essentially because of the significantly different set of criteria imposed for the selection of the TGs (see Table 4). According to our results, the discrepancy is specifically to be attributed to the assumption that the NA West coast and SE North America TGs are located outside the peripheral bulge of the formerly ice-covered regions. While this can be certainly assumed for GIA computations based on model ICE–3G, our computations have shown that this is not the case when models ICE–5G or KL05 are employed. This serious inconsistency between GIA model predictions has ultimately prompted us to substitute criterion (5) of D97 with our requirement (V) (i.e. the GIA correction should be essentially model independent).

Since our GMSLR estimate is largely independent of assumptions regarding the GIA corrections, we have possibly improved its robustness with respect to previous studies. However, some problems still remain open. One of these concerns possible sea level accelerations (e.g. Douglas 1992) or steep variations of the rate of sea level rise, or change points, during the last century (Church & White 2006). Detecting these changes would be of paramount importance to fully assess the effects of global warming on the Earth's climate. According to our computations based on SG01, we have found no evidence of a secular acceleration or changes in the rate of sea level change in the mid of the last century. Of course, this does not rule out the possibility of an acceleration, and could merely reflect the inadequacy of the SG01 set as a tool to detect it. A second major problem is that of the possible influence of tectonic movements. Though it is now recognized that coseismic deformations may contribute significantly (i.e. at the 0.1–0.2 mm yr⁻¹ level) on the D97 set of TGs (Melini et al. 2004), the effect on the SG01 TGs, hence on our preferred GMSLR estimate, remains to be addressed.

Acknowledgments

Very constructive comments by Jerry Mitrovica and an anonymous reviewer have greatly helped to improve the manuscript. We have particularly benefited from stimulating discussions with Marco Olivieri, Valentina Barletta, Mark Tamisiea and Philippe Huybrechts. We also thank Florence Colleoni, Daniele Melini, Gian Luigi Alberghi, Federica Fabbri and Roberto Casadio for various comments and suggestions. GS also acknowledges Sofia and Francesco for their support, and Paola for her patience. To facilitate intercomparison of results, code SELEN (v. 2.9) is freely available from the web page http://www.fis.uniurb.it/spada/SELEN_minipage.html. All figures have been drawn using the GMT public domain software of Wessel & Smith (1998). Kurt Lambeck is acknowledged for kindly providing the ice sheets chronology of the KL05 model. This work is supported by COST Action ES0701 ‘Improved Constraints on Models of Glacial Isostatic Adjustment’. We acknowledge the ice2sea project, funded by the European Commission's 7th Framework Programme through grant number 226375 (ice2sea manuscript number 074). We have benefited from the warm atmosphere of the ‘Dopo lavoro ferroviario’ of Lugo di Romagna, where much of the revision work has been carried out.

1

Data are obtained from the archive http://www.psmsl.org/data/obtaining/rlr.annual.data/rlr_annual.zip, extracted from database 2011 February 9.