Spatial sea-level reconstruction in the Baltic Sea and in the Pacific Ocean from tide gauges observations

Exploiting the Delaunay interpolation, we present a newly implemented 2-D sea-level reconstruction from coastal sea-level observations to open seas, with the aim of characterizing the spatial variability of the rate of sea-level change. To test the strengths and weaknesses of this method and to determine its usefulness in sea-level interpolation, we consider the case studies of the Baltic Sea and of the Pacific Ocean. In the Baltic Sea, a small basin well sampled by tide gauges, our reconstructions are successfully compared with absolute sea-level observations from altimetry during 1993-2011. The regional variability of absolute sea level observed across the Pacific Ocean, however, cannot be reproduced. We interpret this result as the effect of the uneven and sparse tide gauge data set and of the composite vertical land movements in and around the region. Useful considerations arise that can serve as a basis for developing sophisticated


Introduction
Sea-level variations are a mark of climate change at different spatiotemporal scales.This stems from the sensitivity of sea level to temperature and salinity variations in the oceans and to mass fluctuations in consequence of the waning and waxing of the continental ice sheets [Cazenave and Nerem 2004, Milne et al. 2009, Cazenave and Remy 2011].Since 1993, altimetry has provided detailed 2-D maps of sea level over time in a geocentric reference frame, which have been used to assess a rate of global mean sea-level rise (GMSLR) of 3.3±0.4mm yr −1 during 1993-2015 (see Nerem et al. [2010] and updates; http://sealevel.colorado.edu/content/global-mean-sea-level-time-series-seasonal-signalsremoved).However, this rate is not representative of the long-term (century) time scale GMSLR.In fact, sealevel change does not occur steadily due to the effect of high-energy quasi-periodic decadal fluctuations [Sturges and Hong 2001], of the time-varying contribution of the melting of ice sheets and glaciers [Bam-ber and Riva 2010] and of thermosteric and halosteric regional effects [Levitus et al. 2005].
Long records from tide gauges (TGs) deployed along the world's coastlines constitute the only observational evidence of long-term sea-level change (see Spada and Galassi [2012], and references therein).However, these records are scarce and the TGs are unevenly distributed geographically, with a dominant concentration in the northern hemisphere [Holgate et al. 2013].Furthermore, TG observations provide sea level relative to the coasts, which is not directly comparable to the absolute sea level detected by satellite altimetry.By using TG data properly corrected for the effects of the glacial isostatic adjustment (GIA), several sea-level reconstructions have been obtained so far, which have provided various estimates of GMSLR (see the review of Spada and Galassi [2012]) and of its acceleration during last century (see Spada et al. [2015] for a summary).The approaches followed so far to estimate the rate of sea-level change from TGs broadly fall into three categories: i) 1-D sea-level reconstructions based on various linear combinations of records as for the "virtual station" method [Jevrejeva et al. 2008] or standard stacking [e.g., Olivieri and Spada 2013]; ii) 2-D reconstructions by means of the empirical orthogonal functions (EOF) method derived from altimetry [Church et al. 2004] or ocean reanalyses [Calafat et al. 2014], and iii) spatial interpolation of TG observations [Choblet et al. 2014].Each of them suffers from some limitations.Moreover, since the "altimetry era" started in 1992, a comparison between TG-related models and altimetry for a time span longer than 20 years is not possible.
Methods i) and ii) above, with different implementations, have been extensively applied [see Spada et al. 2015] although the problem of extending the TG signals to open oceans, and the relationship with the record from satellite altimetry has not been fully solved yet [Calafat et al. 2014].Method ii) however, appears to be the more promising approach to model sea level and to extend to the pre-altimetry era the 2-D sea-level reconstruction.The validity of EOF-based models is commonly assessed comparing the resulting sea level to TG time-series.This is the case for the so-called CSEOF (cyclo-stationary empirical orthogonal functions; see Hamlington et al. [2011]) that evolved the standard EOF approach.The CSEOF method was applied to the complex case of the Southeastern Asian Sea [Strassburg et al. 2014] showing that long-term sea-level change is dominated by the Pacific decadal oscillation (PDO Index; Mantua and Hare [2002]).Further methodologies, in particular the objective analysis [Bretherton et al. 1976, Davis 1985, Robinson and Leslie 1985, Carter and Robinson 1987, Robinson et al. 1991] and the datainterpolating variational analysis (DIVA) [Brasseur et al. 1996, Troupin et al. 2012] where rarely applied, as discussed below, to sea-level reconstructions from coastal observations.
Here the focus is on method iii), since it has been only seldom exploited.The purpose is that of testing the viability of one of the existing spatial interpolation methods to delineate a 2-D map of sea level from local observations at coastal TGs, and of illustrating its main strengths and weaknesses.Spatial interpolation roots on the need of predicting a spatial variable in un-sampled locations and, broadly speaking, it uses a set of observations available at points located on a plane or across a portion of the sphere to reconstruct the entire 2-D field.A priori assumptions are generally adopted, such as the continuity of the field and of its derivatives, in order to constrain the resulting 2-D model.The details of the different implementations and applications to different physical and socio-economical spatial variables are not discussed in this work.The reader is referred to the works of Lam [1983], Robeson [1997], and Li and Heap [2014] for comprehensive reviews.Despite the fact that spatial reconstruction is a common practice in the field of ocean modeling since decades (see, e.g., Bretherton et al. [1976], Brasseur et al. [1996], and subsequent applications), the reconstruction of the sea-level variability using TG time-series has rarely been attempted (see, e.g., the local-scale case of the Estonian Sea; Raudsepp et al. [1999]).Recently, Choblet et al. [2014] implemented a Bayesian inference method based on a Voronoi [1908] tessellation of the oceans surface to reconstruct sea level during the twentieth century.The validity of the approach is manifest, but the choice of merging nonoverlapping TG observations and of not correcting for vertical land movement (VLM) hinders the comparison of the results with other reconstructions.
In this work, we take a fresh approach to the problem of reconstructing the spatial variability of sea level from TG data.Following Choblet et al. [2014], the idea is based on the work of Bodin et al. [2012], who reconstructed the 2-D structure of the Moho discontinuity beneath Australia.However, differently from Choblet et al. [2014], we do not impose an a priori tessellation of the reconstructed domain but we use the available TGs as nodes for a standard Delaunay triangulation [Delaunay 1934] at which the value of the rate of sea-level change is observed.Then, assuming continuity, smoothness and tension over the surface, we interpolate sea level by means of a curve-fitting method [Renka 1987[Renka , 1997b]].Similar approaches were used in the past to spatially interpolate ocean models [Ringler et al. 2013] and GPS data [Kass et al. 2009], but also for modelling wireless networks [Ghosh et al. 2007] or astronomical observations [Cardiel et al. 2011].The choice of Delaunay triangulation is motivated by its simplicity and by the ease of implementation.The tests we have conducted have been useful for learning the goods and odds of spatial interpolation applied to sea level from TGs.For these reasons, comparisons with other spatial reconstruction methods are not in the scope of our work.
In Section 2, we first apply the interpolation method to TG data for the Baltic Sea, during the time period 1993-2011.Then, we extend the results to the case study of the Pacific Ocean.The conclusions and future directions are drawn in Section 3.

The Baltic Sea
As pointed by Robeson [1997], the low sampling density is a major factor causing the inability of interpolation methods to reproduce detailed information.For this reason, our main focus will be on the Baltic Sea, since this small basin is well sampled by TGs.The Baltic Sea region is also of particular geodynamical interest because Fennoscandia is still experiencing a significant uplift rate (several millimeters per year) in consequence of the Late-Pleistocene deglaciation (see Steffen and Wu [2011], and references therein).Due to the large rates of crustal uplift associated with GIA, in this region a large discrepancy is expected to exist between the relative sea level observed at TGs and the absolute sea level revealed by altimetry.At the time scale of decades (and even centuries), Fennoscandia can be considered tectonically stable [Steffen and Wu 2011] and other potential sources of VLM as sedimentation and soil compaction, of importance in large sedimentary river deltas [Ericson et al. 2006], can be neglected in this region [Kuo et al. 2004].
The Baltic Sea is a small semi-enclosed basin that covers an area of ~415,000 km 2 , with a maximum width not exceeding 200 km and a total coastline length of ~70,000 km.More than one hundred TG stations [Holgate et al. 2013] provide the densest regional coverage worldwide.For this study, we have selected all the RLR (revised local reference; see Woodworth and Player [2003]) records available from the PSMSL (permanent service for mean sea level [PSMSL 2013]) for TGs stations in the Baltic Sea.In addition, to ensure a sufficient coverage of the western boundary of the Baltic Sea, we have also considered those along the coasts of the Kattegat and Skagerrak straits.The PSMSL database contains 111 time series from RLR TGs within the region but only 64 of them were operational during the time period of concern here (1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011), providing annual mean values.Gaps in some of the time-series decrease the effective amount of records available for a given year to a number varying between 57 and 63, with an average of 59.The TG time series are time differentiated by applying of a two-points two-sided derivative [e.g., Wood 1982] so that the ensuing analyses will be performed in the velocity domain.This removes possible biases due to the use of different TG benchmarks but it leads to a significant error on the annual rate of sea-level change (~0.3 mm yr −1 ; F. Raicich, personal communication, 2014).
Among the different interpolation methods [Robeson 1997], we choose the Delaunay triangulation [Brassel and Reif 1979] since this produces the simplest network connecting a set of points over a given domain [Manni et al. 2004].This choice reflects the idea of keeping the approach as simple as possible, without imposing constraints on the density and distribution of nodes that are common in other methods.We use a modification of the original method by Brassel and Reif [1979], implemented into the package sphinterpolate of GMT [Wessel et al. 2013].This program performs gridding under tension on a sphere with the assistance of the two packages SSRFPACK and STRIPACK developed by Renka [1997a,b].In our implementation, the input is a set of points of colatitude i and longitude m, at which the rate of sea-level change ṡ is known.The ṡ values are interpolated on the sphere with the possibility of including or excluding tension and of choosing between different smoothing options (local vs global gradient).In the following, sphinterpolate is used with the option -Q2 that enables tension in order to preserve local monotonicity and convexity, and smoothing is performed by means of global gradient.The Delaunay triangulation, and consequent interpolation, is performed yearly in the time span 1993-2011.Gaps in the time series modify, from year to year, the mesh ob-tained from triangulation.The final output of this process is a suite of yearly 2-D models for the rate of sea-level change.As an example, Figure 1 shows the model obtained for year 2002, in which the nodes are the locations of the TGs available for that year, the triangulation is displayed, and the color scale map represents the 2-D reconstruction of the rate of sea-level change.We remark that the simplicity of the chosen method has a major drawback that could harm the interpretation of the results.In fact, it does not provide an estimate of the uncertainties, and thus it is unable to propagate errors through the interpolated fields.
The large variability of the rate of sea-level change for year 2002, apparent in Figure 1, reflects the wide range of rates observed at the TGs for a single year (for 2002, the mean value is −7.4 mm yr −1 and the standard deviation is 9.4 mm yr −1 ). Figure 1 shows, not unexpectedly, that the overall shape and size of the triangles is determined by the distribution of nodes, which, in turn, sensibly affects the reconstruction.We note that general criteria for assessing the quality of the triangulation have not been defined.However, common sense suggests that small (with respect to the scale of the SEA-LEVEL RECONSTRUCTIONS sampled region) and equilateral triangles are more efficient than wide, acute or obtuse, ones [Chen and Xu 2004].However, in some specific cases, thin and acute triangles can provide better results [Rippa 1992].In our case study, the triangles are mostly small and nearly equilateral, but there are some exceptions as in the Gulf of Finland.Here, the absence of data along the Estonian coast results in one long and very acute triangle, marked in orange in Figure 1, whose southwestern vertex is not located at the coast of the gulf.Then, considerable portions of the sides of this triangle are on land, thus locally leading to a poor mesh.
To test the methodology outlined above, we have considered two TG stations, one located along the continental coast (Stockholm) and one in the Gotland Island (Visby).The two stations are both marked by triangles in Figure 1.This test is performed twice.The first time we removed one of the two stations from the input data set, the second time we repeated it by removing the other.Then, the rate of sea level at each location is extracted by the suite of 2-D models obtained by interpolation.This provides a simulated time-series for the rate of sea-level change at the two sites.In Figure 2, we compare the observed rate of sea-level change (data) with the reconstructed one (synthetic) for the two TGs.The agreement is apparent and this is confirmed by the value of the Pearson [1931] correlation coefficient, which is found to be 0.99 for Stockholm and 0.98 for Visby.This example indicates that, on the spatial scale of the Baltic Sea, the Delaunay triangulation can satisfactorily reproduce the 2-D sea-level field in spite of gaps in some of the time-series that cause a time variation of the mesh geometry.We remark that the time-series for TGs sited nearby Visby show a correlation with the Visby time-series comparable with that of its reconstruction.This confirms the low complexity of the basin, as originally observed by Wróblewski [1993] who proved that the first EOF accounts for more than 90% of the variance in the Baltic Sea over a time span of 40 years.
To obtain a model for the mean rate of sea-level change for the Baltic Sea over the two decades 1993-2011, we have stacked (i.e., averaged) the individual images obtained from the yearly 2-D spatial reconstructions following the method of Olivieri and Spada [2013].We did not extrapolate the mean sea-level rate in the basin to compare with mean sea-level rate from linear regression since the absence of associated errors to Delaunay triangulation would prevent a meaningful comparison.The map of the average rate of sea-level change, hereafter referred to as S -TG , is shown in Figure 3a.By a visual comparison with the map for year 2002, shown in Figure 1, the comparatively smaller spatial fluctuations in S -TG are apparent.This is an effect of the averaging over two decades, which smoothes local short lived relative sea-level fluctuations.The north to south gradient shown by S -TG is interpreted as the 2-D fingerprint of GIA, which dominates the relative sea-level change at these latitudes [Steffen and Wu 2011].Since TGs observe sea level with respect to the solid Earth, a straightforward comparison of S -TG with the altimetric rate of sea-level change is not possible, because altimetry observes sea level in a geocentric (or absolute) reference frame [Nerem et al. 2010].The relationship between relative and absolute sea-level change is known as the "sea-level equation" (SLE).In terms of rates, it reads: where N is the rate of absolute sea-level change, U is the rate of vertical displacement of the crust, i is colatitude, m is longitude, and t is time [Farrell andClark 1976, Spada andStocchi 2006].Taking the time average of the fields, the SLE simply reads The average rate of absolute sea-level change from altimetry, N -ALT , is shown in Figure 3b.This has been obtained from the AVISO re-analysis, which merges data from TOPEX/Poseidon, Jason-1 and Jason-2 missions, with no corrections for GIA nor for Inverted Barometer (data are available from: ftp://ftp.aviso.oceanobs.com/pub/oceano/AVISO/indicators/msl/MSL_Map_ MERGED_Global_IB_RWT_NoGIA_Adjust.nc).This re-analysis was preferred among others like CSIRO [Church and White 2011] or ESA CCI [Ablain et al. 2015] since it is widely employed and this facilitates the intercomparison of the results.
Using the SLE (2), in Figure 3c we show the field U -TG = N -ALT − S -TG , which can be compared with the rate of crustal motion observed by global positioning system (GPS), U -GPS , displayed in Figure 3d.This map, which is obtained in a totally independent way from U -TG , results from the spatial interpolation of rates observed at 85 GPS sites for the period 1993-2006 and extracted from Table 1 in Lidberg et al. [2010].Interpolation was performed using the same algorithm used for the sea-level reconstruction but, of course, using different nodes.We note that the time span covered by GPS data is not completely overlapping that of altimetry.As mentioned above, the VLM in this area is dominated by GIA and the other causes have a minor influence [Kuo et al. 2004, Steffen andWu 2011].This is of relevance here, since the occurrence of short-lived or local deformations in the vicinity of the TG or GPS sites, which in our case are not co-located, would compromise the comparison between U -TG and U -GPS .The juxtaposition of U -TG and U -GPS visually confirms that the dominating signal represented by the north-south decreasing pattern is properly reconstructed.Contours for 10 and 5 mm yr −1 occur almost at the same latitude with similar SW-NE trending.Two main differences are however apparent.First, the maximum at the northern edge of the basin is not reproduced although two TGs are providing data from that area.Second, the short-wavelength features shown by U -TG in the Arkona Basin (east of Denmark) are not visible in the U -GPS field, which can be possibly attributed to over-smoothing.

The Pacific Ocean
The encouraging results obtained for the Baltic Sea is posing the obvious question about the extensibility of the approach to other, larger regions.After several attempts, we have found that the Mediterranean Sea is not a viable candidate, mainly because its geometrically complex shape makes the triangulation inefficient and, most importantly, because reliable TG observations are missing from North Africa (PSMSL; Holgate et al. [2013]).In what follows, we extend the application of the interpolation method to the largest possible ocean basin, i.e., the Pacific Ocean.The motivation resides in the significant number of TG observations available from this region during the two decades of interest here, and its large-scale convexity that facilitates, at least in principle, the triangulation.The Pacific Ocean covers a surface ~1.65×10 8 km 2 wide and its maximum extension is ~12,000 and ~15,000 km in the east-west and in north-south directions, respectively.The data-set of RLR TGs operating during period 1993-2011 consists of 104 sites over the total of 206 archived at PSMSL.Comparing with the case of the Baltic Sea, the vastness of the basin combined with the reduced density of TGs suggests a lower resolution in the resulting model and, in general, a more cautious approach.
For the case study of the Pacific Ocean, we have interpolated the mean rate of sea-level change observed at each TG.This simplified approach is suggested by the gaps existing in individual time-series, which could change drastically the yearly distribution of nodes and, consequently, the mesh geometry.Compared with the case of the Baltic Sea, the gaps could have potentially more severe effects, in view of the longer sides of the triangles involved across this large ocean domain.We have selected all of the RLR TGs facing the Pacific Ocean with completeness ≥ 80% during 1993-2011.For each of them, we have determined the mean sea-level rate as the slope for the best fitting first order polynomial resulting from a standard linear regression.Triangulation and interpolation of this dataset provides the 2-D relative sea-level field drawn in Figure 4a.From a geometrical standpoint, we observe a large number of very acute triangles connecting regions separated by large distances, as for the case of the United States West Coast and Hawaii.This is consequence of the uneven distribution of TGs, dense on the West Coast and sparse in the North Pacific, and of the scarcity of TGs located on islands in the northern hemisphere.This shortcoming could be tapered by downsampling the United States West Coast, possibly my means of the "virtual station" method [Jevrejeva et al. 2008].
We first tested the capability of the reconstruction to reproduce the mean sea-level rate observed at some of the islands during 1993-2011.Results were not encouraging although the unavailability of error bars associated with the reconstruction does not facilitate the comparison.In Figure 4b, the map of the rate absolute sea-level change N -ALT observed from altimetry is shown, obtained from the same AVISO products exploited in Section 2. The map is characterized by short wavelength features (down to ~100 km) that reflect the complex spatial variability of contemporary sea-level rise [Cazenave andNerem 2004, Nerem et al. 2010].Clearly, these cannot be observed in the map reconstructed by TG observations in Figure 4a, mainly because of the poor spatial distribution of the observations.In spite of that, the two maps are showing a broad resemblance, with maxima and minima of N and S mostly located in the West and in the East Pacific, respectively.The large-scale similarity of the two fields suggests a multi-scale analysis that here has been accomplished using the orthonormal functions (ON) method of Hwang [1991Hwang [ , 1993]].ONs are preferable with respect to the traditional EOFs, because the former are uniquely determined by the shape of the oceans over which they are mutually orthogonal, and are independent from a specific geophysical dataset [see, e.g., Spada and Galassi 2016].To capture the very long-wavelength pattern of N -ALT , in Figures 4c and 4d, the altimetry maps have been reconstructed to harmonic degrees l = 3 and l = 10, corresponding to spatial wavelengths > ~4,000 km.After the ON reconstruction, similarities between the sea-level rates in Figure 4a and these maps are more apparent.
The characterization of the VLM for the case of the Pacific Ocean is more complex than for the case of the Baltic Sea, where GIA dominates and other components are small in comparison.Here GIA, a longwavelength phenomenon occurring at a constant rate at time scale of several decades and even centuries, acts only at high latitudes.On the contrary, other processes at shorter spatial and temporal scale provide large contributions to VLM.This is the case of tectonic deformation accompanying the seismicity and the volcanism along the "ring of fire".As an example, we focus on Japan whose TG records are crucial to cover with properly shaped triangles the northwestern sector of the Pacific Ocean (see Figure 4a).Here, however, tectonic motions are dominated by short-lived local deformations in consequence of large thrust earthquakes [Ekström et al. 2012] and of the inflation preceding volcanic eruptions [e.g., Miura et al. 2000].In this case, separate regional interpolations of GPS and TG data would provide two models that could not be compared because they are widely affected by local processes.In this context, the use of co-located GPS and TG instruments becomes crucial to convert relative sea level into absolute sealevel coatlines [Wöppelmann and Marcos 2016].

Conclusions
We have studied the usefulness of a spatial interpolation technique to generate realistic 2-D models for rates of relative sea-level change, when input data are rates observed along the coastline.The technique has been applied to the Baltic Sea and to the Pacific Ocean.For the case of the densely instrumented Baltic Sea, the comparison with rates of absolute of sea-level change observed by altimetry has shown the efficiency of the approach, although convex portions of the sea basins can be erroneously modeled when only a few TG stations are nearby.These results are in agreement with previous works [e.g., Wróblewski 1993] showing that sea-level variability is dominated by basin-scale wavelength.This goes along with the large coherence observed at different TGs in the Estonian Sea [Samuelsson and Stigebrandt 1996].We regret to mention that comparison with previous works is difficult because, with the exception of the work by Raudsepp et al. [1999], to our knowledge none of the available 2-D sea-level reconstruction have relied only upon TG time-series.The case study of the Pacific Ocean, at a much wider scale, has evidenced that the spatial resolution of the model should be calibrated in accordance with the density of the input data and that unevenly distributed TGs can lead to an erroneous model; in this case a downsampling should be preferred.To address some of the weaknesses of the method, improvements should be aimed to limit the interpolation to a specific domain and the inclusion of artificial nodes to force the triangulation out of land.However, the most remarkable limitation for this approach is the fact that it does not allow to estimate the uncertainties to be associated with the spatial reconstruction.This suggests that different approaches could be preferred.On the contrary, one of the goods of this method is that the tessellation of the reconstructed region depends only on the distribution of the nodes, which prevents any artifacts arising from possible arbitrary choices.Concerning its usability in different basins, we remark that this 2-D reconstruction from coastal observations necessarily requires input data all along the coast surrounding the basin, and an almost convex shape for the basin itself.In the current implementation this would strongly limit its usability.For example, the Mediterranean does not hold any long RLR TG time-series along a major portion of the African coast that hinders the application of Delaunay interpolation there.In this direction the recovery of ancient but not publicly available data, as for the case of La Valletta (MT) site (see discussion in Olivieri et al. [2013]) or Marseille (F) [Wöppelmann et al. 2014] would be very important for a proper reconstruction of sea level.
In conclusion, this work evidences that a realistic and reliable sea-level reconstruction by means of Delaunay triangulation from TGs would require i) a basin surrounded by coastlines, ii) homogeneous coverage of TGs, iii) TGs located on islands to ensure a dense distribution of small triangles, iv) no VLM with patterns at scales significantly smaller than the scale of the basin or proper spatial corrections.Finally, we remark that a suitable implementation of Delaunay triangulation would require the definition of a method for the reconstruction of the uncertainties assessment.This could be done by bootstrap application of the Delaunay triangulation randomly picking values from each TG estimated error-bar.Kriging [Matheron 1976] could be a viable alternative because it accommodates for the inhomogeneity in the spatial distribution of data and it also provides estimates for the uncertainties.The co-location of TG and GPS instruments would then facilitate the comparison between S and N and, indeed, ease the modelization of absolute sea-level change from TG observations.
Acknowledgements.Sea-level data have been downloaded from the Permanent Service for Mean Sea Level database on August 1st, 2012 (http://www.psmsl.org/data/obtaining/).The altimeter products were produced by Ssalto/Duacs and distributed by AVISO (Archiving, validation and interpretation of satellite oceanographic data) from http://www.aviso.altimetry.fr/duacs/ with support from CNES (Centre National d'Etudes Spatiales).Rates of vertical land motion for the Baltic region were extracted from Table 1 of Lidberg et al. [2010].GS is supported by a DiSPeA research grant (CUP H32I160000000005) and by Programma Nazionale di Ricerche in Antartide (PNRA 2013/B2.06,CUP D32I14000230005).All figures, except some free-hand insertions, have been drawn using the generic mapping tools (GMT) of Wessel et al. [2013].The use of "sea level" as a noun and "sea-level" as an adjective follows the editorial by Fairbridge [1997].

Figure 1 .
Figure 1.Map of the Baltic Sea region.Red triangles indicate the locations of the 57 TGs available for triangulation in the year 2002.Colors describe the 2-D reconstruction for the rate of relative sea-level change for the same year.Black triangles mark the sites (Stockholm and Visby) for which the observed sea level has been reconstructed.

Figure 2 .
Figure 2. Relative sea-level reconstruction (red) at the two target sites of Stockholm (top) and Visby (bottom), compared with the rate of relative sea-level change (black) obtained from the raw TG observations.The dashed line represents the difference between the two time series.

Figure 3 .
Figure 3. (a) Average rate of relative sea-level change reconstructed from the TG time series during 1993-2011; (b) average rate of absolute sea-level change from altimetry in the same period; (c) rate of vertical displacement given by (b)-(a); (d) spatial reconstruction of GPS observations of the rate of crustal uplift.Except in (b), all the color tables have the same ranges.

Figure 4 .
Figure 4. (a) Reconstructed rates of relative sea-level change across the Pacific Ocean; (b) rate of absolute sea-level change observed from altimetry; (c) altimetry map reconstructed from the ON coefficients to harmonic degree l = 3; (d) the same as in (c), but to degree l =10.