CC-BY 4.0
1. Introduction
Sea level rise is an important consequence of anthropogenic emissions of greenhouse gases (Fox-Kemper et al., 2021; Meyssignac et al., 2023). The question of sea level rise projections is crucial for coastal planning (Nicholls et al., 2021) and uncertainties on sea level projections are a major subject of concern (Hirschfeld et al., 2023). The main way to explore sea level in the future relies actually on numerical simulations; this approach is useful for projections until the end of the century and further (Le Cozannet et al., 2023). However, sea level rise observations are currently poorly used for future evaluations of sea level rise or for models simulations constraints such has been done with observations for global mean or local temperatures (Qasmi and Ribes, 2022; Ribes et al., 2021), in spite of the availability of sea level rise reconstructions and quadratic fits from tide gauges and/or space altimetry (e.g. Dangendorf et al., 2019; Nerem, Beckley, et al., 2018; Nerem, Frederikse, et al., 2022; Guérou et al., 2023, see also https://sealevel.colorado.edu/). Recently the United States Sea Level Rise Task Force proposed to quadratically extrapolate combined US tide gauges observations until 2050 (Sweet et al., 2022), based on the fact that the main causes of sea level rise are still on the road and that IPCC SSP scenarios can poorly be distinguished in terms of sea level rise before the second half of 21st century (Fox-Kemper et al., 2021). This operation has two justifications: first, sea level trends and accelerations estimates, detection and attribution is of major interest for climate change monitoring, climate services and coastal adaptation planning for next decades; second, combination of observations has several advantages: it diminishes the internal variability of sea level observed at only one site, it fills gaps in data, and it allows to characterize uncertainties of sea level observations.
Nevertheless, these two hypothesis pose some questions. Trends and accelerations have until recently been explored by different ways. For example, Visser et al. (2015) reviewed 30 trends estimates methods. They underlines the importance of the method chosen for the estimate (e.g. the time length), the way for uncertainties handling, the possibility to add external information (such as climate indices), and the ability to extrapolate sea level. In particular, past studies have estimated trends and/or accelerations by fixing different time lengths, e.g. 10 years for individual tide gauges trends (Holgate and Woodworth, 2004), 30 years for GMSL (Frederikse et al., 2020). However, systematic estimate of sea level trends and accelerations requires continuous datasets to explore all possible time spans, which is easy to obtain at global scales, either from long reconstructions based on tide gauges data (Jevrejeva, Moore, Grinsted and Woodworth, 2008; Haigh et al., 2014) or from shorter, but continuous, time series from space altimetry (Ablain et al., 2019; Guérou et al., 2023). At local scale, sea level variations can be measured either with space altimetry (Prandi et al., 2021) or with tide gauges. However, individual tide gauges have data gaps which require to be filled. Data gaps filling is a very wide scientific field, with various methods which have been developed, which this paper won’t review. Nevertheless, methods inspired by statistical interpolation (Asgharinia and Petroselli, 2020; Kandasamy et al., 2013) or machine learning Arriagada et al. (2021), and Kumar et al. (2021) have been developed; the former has the drawback to filter high frequencies, although the latter needs large learning database. Remote sensing has also its own methods, including a “shift and scale” approach for the study of leaf area index from satellites (Verger et al., 2013), which is fundamentally a close and simplified model to the Helmert seven parameters widely used in the geodesy subdomain of reference frames (Altamimi, Sillard, et al., 2002; Altamimi, Rebischung, et al., 2023).
In the case of tide gauges, the method to chose for gap filling must take the uses of long term mean sea level observations into consideration. A large part of sea level studies deal with the separation of internal variability and long term signal (Douglas, 1991; Jevrejeva, Moore, Grinsted and Woodworth, 2008; Jevrejeva, Moore, Grinsted, A. P. Matthews, et al., 2014; Church and White, 2011; Woodworth et al., 2009). Another, less known, application of long term sea level data at tide gauges is hydrodynamic leveling, which aims to assess slopes of vertical datum based on spirit leveling (Penna et al., 2013; Featherstone and Filmer, 2012), and to connect vertical datums of territories separated by the sea (Filmer and Featherstone, 2012), by the combination of long tide gauges time series and mean dynamic topography models.
However, mean sea level observed at tide gauges is affected by several effects which have to be assessed on long term mean, trends and accelerations estimates. A first effect comes from vertical land motion, which has long term consequences on relative sea level variations (Wöppelmann and Marcos, 2016). Another effect comes from the seasonal atmospheric influence on sea level, named “inverse barometer” (Wunsch and Stammer, 1997; Woodworth et al., 2009), which some studies do not correct for in hydrodynamic leveling Featherstone and Filmer (2012), although some others correct it either from atmospheric reanalysis (Penna et al., 2013) or by only selecting summer data (Afrasteh et al., 2023). Finally, astronomical cycles affect mean sea level observations on seasonal to multi-decadal time scales (Woodworth, 2012), and have also to be corrected for in hydrodynamic levelling, for instance if required to transform raw observations into mean tide observations (Mäkinen and Ihde, 2009). However, to our knowledge, no systematic study to assess and compare the consequences of these, eventually combined, corrections on long term mean, trends and accelerations has been carried out.
Finally, tide gauges observations have been used in many studies to discuss the acceleration of sea level rise, notably at global scale (Visser et al., 2015; Merrifield et al., 2009). Douglas (1991) recommended for tide gauges time series to be at least 60 years long to assess accelerations. Jevrejeva, Moore, Grinsted and Woodworth (2008) and Jevrejeva, Moore, Grinsted, A. P. Matthews, et al. (2014) reconstructed global then regional mean sea level back to 200 years to systematically explore global mean sea level trends and accelerations, and found a quasi-periodic cycle on such parameters with a typical timescale of 60 years. Then Houston and Dean (2011) did not conclude to any acceleration on the basis of 20th century tide gauges observations, but Rahmstorf and Vermeer (2011) demonstrated that this conclusion relied only on short time spans or regional analysis by systematically explore all time series based on Church and White (2011) reconstruction. In addition, Dangendorf et al. (2019) claimed that GMSL accelerates since the 1960s. Finally, tide gauges observation-based trends and accelerations have been used to extrapolate local mean sea level and assess uncertainties at timescales of decades, for instance in the USA (Ezer, 2013; Sweet et al., 2022), in Europe (Melet et al., 2023), or in single countries, such as Portugal (Antunes, 2019; Antunes and Taborda, 2009).
This paper addresses three science questions. The first one is individual tide gauges data gaps filling, taking into account the two constraints of the long term signal preservation and of the inclusion, as good as possible, of local internal variability. The second science question is the evaluation of the effects of all possible corrections on tide gauges mean sea level observations on trends and accelerations estimates. The third science question is the use of continuous tide gauges observations from gap filling step, possibly corrected for diverse effects, to systematically explore all time spans estimates of trends and accelerations to select the best fit to extrapolate sea level for next decades and constraint climate projections based on simulations. The paper is structured as follows: data and methods are described in Section 2, then results and discussion are presented in Section 3.
2. Data and methods
2.1. Data
We used monthly sea level RLR data from the Permanent Service for Mean Sea Level (Holgate, A. Matthews, et al., 2013) for tide gauges of Brest (PSMSL id. 1), Marseille (id. 61), Genova (id. 59), Sète (id. 958), Toulon (id. 980), Nice (id. 1468), Port-Vendres (id. 1469), l’Estartit (id. 1764), Saint-Jean-de-Luz (id. 469), La Rochelle—La Pallice (id. 466), Saint-Nazaire (id. 457), Brest (id. 1), Le Havre (id. 453), Boulogne-sur-Mer (id. 471), Calais (id. 455), Dunkerque (id. 468), and Newlyn (id. 202). We selected all data from epoch 1885.0 to 2023.0, covering 138 years of data. We also used the following climate indices: Atlantic meridional mode (AMM) (Chiang and Vimont, 2004), Atlantic multidecadal oscillation (AMO) (Enfield et al., 2001), Arctic oscillation (AO) (Higgins et al., 2002), North Atlantic oscillation (NAO) (Hurrell, 1995; Jones et al., 1997).
In addition, we downloaded IPCC AR6 WG1 regional sea level rise projections (Fox-Kemper et al., 2021) based on emulators FAIR (Smith et al., 2018) for the thermosteric component of sea level rise, and FACTS (Kopp et al., 2023) for the mass component of sea level rise, as computed by NASA (Garner et al., 2021), for all proposed socio-economic scenario (i.e. scenarios with medium confidence processes: SP1-1.9, SSP1-2.6, SSP2-4.5, SSP3-7.0, SSP5-8.5, and scenarios with low confidence processes: SSP1-2.6-low and SSP5-8.5-low). These projections are provided at a 2-degrees resolution, and are referred to mean sea level over 1995–2014; as proposed on the download interface, we selected projections for the selected tide gauge id. 1 (Brest) and id. 61 (Marseille).
2.2. Gap-filling combination model
In this subsection we propose a gap-filling method, closed to the “shift and scale” approach mentioned in the introduction, albeit based on an adaptation of the seven-parameters Helmert transformation model widely used in geometrical geodesy. We reduce it to two parameters needed to combine scalar datasets. We will apply it to our test-cases of tide gauges data.
Let Y be an observable variable; let S be the number of datasets (individually denoted s) each observing Ns epochs denoted tis. Let Ys be the vector of individual values yis. There are a total of N epochs (ti).
Our combination model relies on a transformation model where we consider two datasets Y1 and Y2 of the same observable variable. We write the transformation model as an adaptation of the well-known seven-parameters similitude, here adapted and reduced to only one translation T and one scale factor F:
| \begin {equation}\label {eq1} Y_1 = T + F\cdot Y_2 \end {equation} | (1) |
| \begin {equation}\label {eq2} Y_1 = T + Y_2 + D\cdot Y_2 \end {equation} | (2) |
Here, our aim is to estimate combined variables (Yc ≡ Y1) and transformation parameters Ts and Ds between a combined dataset c and individual datasets s. First we linearize the unknown vector Yc = Y0 +𝛿 Y . It comes for the transformation model:
| \begin {eqnarray} &\displaystyle Y_1 = Y_0 + \delta Y_2 + T + D\cdot Y_0 + \underbrace {D\cdot \delta Y}_{\approx 0} \label {eq3}\end {eqnarray} | (3) |
| \begin {eqnarray} &\displaystyle \iff \delta Y_1 = \left [ \begin {array}{ccc} I_N & \mathbb {1} & Y_0 \end {array}\right ] \cdot \left [ \begin {array}{c} \delta Y_2\\ T_2\\ D_2 \end {array} \right ] \label {eq4} \end {eqnarray} | (4) |
Generalizing this model to S datasets requires to build an observation matrix instead of the above-noted identity matrix. The observation matrix of dataset s is noted Hs and has N columns and Ns rows. Hs has zeros everywhere unless unknown variable yic referring to epoch ti of the combined solution is observed by solution s; in this case, the corresponding element of Hs is equal to 1. Therefore:
| \begin {eqnarray} \left [ \begin {array}{c} \delta Y_1\\ \vdots \\ \delta Y_S \end {array} \right ] & = & \left [ \begin {array}{cccc} H_1 & A_1 & & {0}\\ \vdots & & \ddots & \\ H_S & {0} & & A_S \end {array} \right ] \cdot \left [ \begin {array}{c} \delta Y_c\\ \theta _1\\ \vdots \\ \theta _S \end {array} \right ] \label {eq:combinationmodel}\end {eqnarray} | (5) |
| \begin {eqnarray} \mbox {with}\quad A_s & = & \left [ \begin {array}{cc} \mathbb {1}_s & Y_{0s} \end {array} \right ] \end {eqnarray} | (6) |
| \begin {eqnarray} \mbox {and}\quad \theta _s & = & \left [ \begin {array}{c} T_s\\ D_s \end {array} \right ] \end {eqnarray} | (7) |
| \begin {eqnarray} \Sigma = \left [ \begin {array}{ccc} \Sigma _1 & & {0} \\ & \ddots & \\ {0} & & \Sigma _S \end {array} \right ] \end {eqnarray} | (8) |
Under this form, the problem is however poorly conditioned, and needs to be constrained. We have to write constraints conditions equations, and to complete the above observation equations system to force the solution to respect some properties. Here we have a reference tide gauge dataset r and we want its properties T and D not to be altered. The constraints equations system is written:
| \begin {eqnarray} \left [ \begin {array}{c} K_1\\ \vdots \\ K_S \end {array} \right ] = \left [ \begin {array}{c} I_{2S}\\ \end {array} \right ] \cdot \left [ \begin {array}{c} \theta _1 \\ \vdots \\ \theta _S \end {array} \right ] \end {eqnarray} | (9) |
| \begin {eqnarray} \Sigma _K = \left [ \begin {array}{ccc} \sigma _{K_1}^2 & & {0} \\ & \ddots & \\ {0} & & \sigma _{K_S}^2 \end {array} \right ] \end {eqnarray} | (10) |
We now have a well-conditioned inverse problem that can be solved by a classical weighted least-square method; if we write the problem B = M ⋅ X, then the solution is simply $\hat{X} = \left (M^t\cdot P\cdot M\right )^{-1}\cdot M^t\cdot P\cdot B$, where P, the weighting matrix, is the inverse of the variance matrix (i.e. $P=\Sigma _K^{-1}$).
This combination model is applied in two different cases. The first one aims to fill gaps in Marseille tide gauge time series with neighboring tide gauges data from Genova, Sète, Toulon, Nice, Port-Vendres, l’Estartit. The second one aims to fill gaps in Brest tide gauge time series with neighboring tide gauges data from Saint-Jean-de-Luz, La Rochelle—La Pallice, Saint-Nazaire, Le Havre, Boulogne-sur-Mer, Calais, Dunkerque, and Newlyn. Brest and Marseille transformation parameters are, each on its own gap-filling process, fixed to zero in order not to alter their long-term behavior, so that short-term internal variability come from the combination of their data with other tide gauges. On their data gaps, internal variability only come from other tide gauges combined signals.
2.3. Corrections applied to gap-filled tide gauges time series
After the combination process, we have two sea level time series, one in Brest, the other in Marseille. The post-processing for both is the same, relying first on a quality check of data, then on subtracting several additional non-climate effects. First we reject epochs for which data are missing as indicated in the PSMSL file. Then, we estimated a degree-2 polynomial and reject observations for which the residual of the fit exceeds 20 cm. For instance, for Marseille tide gauge, only 7% of the data are rejected after this two-steps quality check. Data are also expressed relative to the 1995–2014 mean state, which is also the IPCC AR6 sea level reference level (Fox-Kemper et al., 2021). The a-posteriori variance of this quadratic fit is recorded.
The second step of tide gauges post-processing consists in correcting observations from diverse effects. We have tested several configurations of corrections. The first effect to correct is the inverse barometer effect (Wunsch and Stammer, 1997) by using NOAA 20CRv3 (Compo et al., 2011) reanalysis until 2015 (included), and ECMWF ERA5 reanalysis (Hersbach et al., 2020) afterwards (2016 and later). Another effect is the long-period astronomical cycles of period 18.6 years and 8.8 years (Woodworth, 2012). The last effect to correct are seasonal signals, which can be corrected for by estimating periodic signals of period 1.0 and 0.5 year with sine and cosine functions and subtracted from the signal.
In this study, the vertical land motion (Wöppelmann, Marcos, Coulomb, et al., 2014; Wöppelmann and Marcos, 2016) is not corrected, because we are interested in relative sea level rise projections comparisons between observations extrapolations and simulated projections which already incorporate vertical land motion in Fox-Kemper et al. (2021). Moreover, as we will see, our combination model allows to estimate any bias in term of vertical reference level as well as relative secular trends between datasets (such as it could be with different vertical land motions between tide gauges), so that the data of Marseille and Brest tide gauge will not be “polluted” by the data of other tide gauges with which they will be combined.
2.4. Time series fits exploration, extrapolation and milestones timing
With the Brest and Marseille combined, eventually corrected, sea level time series which are built in order to have their respective translation and scale factor null relative to their initial sea heights datasets, we build all possible periods of length longer than 20 years, by shifting consecutive time spans by one year, and consecutive periods lengths by one year. For each of these 7021 (equation to (138 − 20) × (138 − 20 + 1)/2) time series, we estimate first a linear time-regression, and second a quadratic time-regression (note that, for the quadratic regression, the acceleration is the double of the degree-2 estimated parameter). Each fit is calculated by the use of weighted least squares. The weighting matrix is taken as the inverse of the a posteriori variance matrix of observations of the pre-processing, rescaled by the unitary variance factor estimated from the pre-processing initial whole time series fit. In order to take into account the effect of internal variability, we multiply this variance matrix by the weighted RMS squared of the pre-processing fit.
Then, for both linear and quadratic regression models and for periods which ends at the last epoch of the full time span, sea level is extrapolated for each decade until 2100. Uncertainties are rigorously propagated by perturbing polynomial coefficients with a Monte-Carlo random draw of 1000 members, considering, for each, a gaussian distribution for which the best estimate of the coefficient comes from the estimate of the least square fit system, and the standard deviation is equal to the root squared diagonal terms of the a posteriori variance matrix of the least square fit system. We express extrapolated uncertainties as a 66% confidence interval, between 17% and 83% percentiles, as in Fox-Kemper et al. (2021). Note also that the extrapolated sea level is still expressed with respect to the 1995–2014 mean sea level, also as in Fox-Kemper et al. (ibid.).
Finally we look for roots of extrapolation polynomials in order to estimated the timing of some given sea level rise milestones, and we express timing uncertainties as a 66% confidence interval, between 17% and 83% percentiles.
2.5. Correlation of sea level trends and accelerations with climate indices
In order to clarify the physical links of sea level trends and accelerations variability at all timescales explored in previous methodological description subsection, for each period length, we low-pass filter climate indices listed above by considering the corresponding period length as a period cut. For each climate indice, we obtain as many time series as the number of period lengths explored; however, in order to not have only flat time series for longer periods, we limit our study to period length between 30 years and 80 years. Then, we calculate for each filtered climate indice time series the Pearson’s correlation coefficient with the corresponding time series of trends and accelerations extracted from the systematic calculation of these parameters for all possible period lengths and all possible time series central dates.
3. Results and discussion
3.1. Tide gauges datasets gap fillings
The combination of datasets is performed by dividing tide gauges data into two pools; the first one comprises Mediterranean tide gauges data (Marseille, Genova, Sète, Toulon, Nice, Port-Vendres, l’Estartit), and the second one comprises Atlantic tide gauges data (for the french coast: Saint-Jean-de-Luz, La Rochelle, Saint-Nazaire, Brest, Le Havre, Boulogne-sur-Mer, Calais, Dunkerque) to which we add Newlyn tide gauge, in England, to fill common french tide gauges gaps at the end of the 1940s and at the beginning of the 1950s. Our combination model allows “gap-fillers” data to also have data gaps. In this study, we maximize neighboring tide gauges of Brest and Marseille, respectively, to ensure that, altogether, combined solutions won’t have any data gap at the yearly scale.
The combination is performed by initially considering individual datasets diagonal identity variances matrices. Results of the combination are obtained after one re-weighting of the individual datasets variance matrices, but without any outlier rejection. The final unitary variance factor are equal to 0.972 for the Mediterranean pool, and 0.949 for the Atlantic pool. The standard deviation of transformation parameters are equal to the a priori values applied for these standard deviation of transformation parameters. In addition, initial time series of Marseille and Brest data ratio are 93.66% and 93.18%, respectively, before the gap filling step, then 100.0% for both after the gap filling step.
A striking feature of Figure 1 which display the individual tide gauges observations and the combined solution is that the Marseille inter-epoch variability considerably increases after the 1990s. This could be caused either from the increase of the number of tide gauges included in the combination, or to the change of tide gauge device which happened at the end of the 1990s, from mechanical water height observing system to electronic measurement system which are not subject to the mechanical friction of the historical device.
Top, left: Atlantic pool tide gauges sea level time series; top, right: combined Brest gap-filled sea level time series. Bottom, left: Mediterranean pool tide gauges sea level time series; bottom, right: combined Marseille gap-filled sea level time series.
We can see on Table 1 and Table 2 that tide gauges translations with respect to Marseille and Brest, respectively, have various values that reflect the conventional vertical references of different tide gauges. Since tide gauges origins are conventional, these values are of poor interest; however, our combination model could easily be used for the unification of PSMSL Revised local references (RLR). Note that translations are not corrected for on time series displayed on Figure 1, but that combined time series are, by construction, corrected for diverse contributions of individual translations from individual tide gauges origins.
Number of observations, transformation parameters and associated standard deviation (mm), RMS (mm), then final variance factor, estimated at the gap filling step for Marseille tide gauge time series
| Tide gauge | Time span | Nb. obs. | T | D | RMS | $\sigma _0^2$ |
|---|---|---|---|---|---|---|
| std T (1𝜎) | std D (1 𝜎) | |||||
| Genova | 1885–2022 | 1267 | −37.5465 | 0.0098 | 24.1580 | 1.09258 |
| 0.9999 | 0.0010 | |||||
| Marseille | 1885–2022 | 1551 | 0.0000 | 0.0000 | 17.2447 | 0.97086 |
| 0.0010 | 0.0010 | |||||
| Sète | 1992–2022 | 321 | −8.6675 | −0.0097 | 25.4999 | 0.92233 |
| 1.0000 | 0.0010 | |||||
| Toulon | 1961–2022 | 435 | 18.9891 | 0.0001 | 16.2592 | 0.87035 |
| 0.9999 | 0.0010 | |||||
| Nice | 1978–2022 | 474 | 14.2194 | 0.0038 | 27.9493 | 1.06432 |
| 1.0000 | 0.0010 | |||||
| Port-Vendres | 1984–2022 | 363 | −2.2627 | −0.0018 | 20.0580 | 0.85746 |
| 1.0000 | 0.0010 | |||||
| L’Estartit | 1990–2022 | 396 | 10.5637 | 0.0048 | 11.8920 | 0.77326 |
| 0.9999 | 0.0010 |
Number of observations, transformation parameters and associated standard deviation (mm), RMS (mm), then final variance factor, estimated at the gap filling step for Brest tide gauge time series
| Tide gauge | Time span | Nb. obs. | T | D | RMS | $\sigma _0^2$ |
|---|---|---|---|---|---|---|
| std T (1𝜎) | std D (1 𝜎) | |||||
| Brest | 1885–2022 | 1543 | 0.0002 | 0.0000 | 22.8671 | 0.71811 |
| 0.0010 | 0.0010 | |||||
| Newlyn | 1915–2022 | 1265 | −104.4754 | 0.0133 | 24.7139 | 0.74431 |
| 0.9999 | 0.0010 | |||||
| Saint-Nazaire | 1941–2022 | 465 | −71.6149 | 0.0093 | 39.3662 | 0.80143 |
| 1.0000 | 0.0011 | |||||
| La Rochelle | 1941–2022 | 513 | −51.9616 | −0.0103 | 36.2872 | 0.81203 |
| 1.0000 | 0.0011 | |||||
| Saint-Jean-de-Luz | 1942–2022 | 541 | 78.8907 | −0.0267 | 47.0729 | 0.91168 |
| 1.0000 | 0.0011 | |||||
| Le Havre | 1938–2022 | 695 | 15.7387 | −0.0083 | 39.4882 | 1.15245 |
| 1.0000 | 0.0010 | |||||
| Boulogne-sur-Mer | 1941–2022 | 427 | −15.6859 | −0.0067 | 56.1212 | 1.20251 |
| 1.0000 | 0.0011 | |||||
| Calais | 1941–2022 | 540 | 153.9639 | −0.0289 | 67.8163 | 1.20985 |
| 1.0000 | 0.0011 | |||||
| Dunkerque | 1942–2022 | 666 | 188.6161 | −0.0449 | 63.4725 | 1.19975 |
| 1.0000 | 0.0011 |
Scale parameters are of much greater interest, since they are a representation of relative sea level rise trends differences between tide gauges, reflecting different marine processes (such as higher sea level rise rates), land subsidence, or instrumental drifts. Here we don’t provide analysis on these estimates, which will be done in future studies, because our combination model in only used in a gap-filling purpose. We only note that, only Newlyn and Saint-Nazaire have lower (positive scale parameter) trends relative to Brest, based on their common epochs on the full reference time span which is Brest time series time span. In the Mediterranean pool, only Sète and Port-Vendres, in the Gulf of Lion, have higher (negative scale parameter) trends relative to Marseille, based on their common epochs on the full reference time span which is Marseille time series time span.
From now up to the end of the article, we refer to “Brest” and “Marseille” as to the gap-filled combined tide gauges datasets aligned on initial Brest and Marseille, respectively, tide gauges datasets. In this sentence, “aligned” means that the final combined time series have null transformations parameters with respect to the initial Brest and Marseille time series.
3.2. Sea level rise corrections
We have combined all possible combinations of corrections (see Supplementary Material Tables 2, 3, and 4), it appears that some cases are of most interest, namely the uncorrected time series, the inverse barometer (IB) only corrected time series, and the inverse barometer with all residual cycles corrected time series (see Figure 2 for Brest and Figure 3 for Marseille). Indeed, as expected, due to their periodic nature, cycles have little effect on trends and accelerations best estimates. Always as expected, these corrections are more important for shorter time series than for longer time series. Finally, periodic components corrections lower uncertainties on trends and accelerations, for about a quarter to a third of the uncorrected coefficients standard deviations. Periodic components amplitudes and phase estimates are robust among different configurations (see Supplementary Material Table 1). The main difference when combined with or without the inverse barometer effect is for periods of 0.5 year, 1.0 year and 8.8 years. In Brest, the 8.8 year component amplitude is about 6.5 mm without IB correction, but only 3.0 mm when estimated after IB correction; the annual component amplitude is 47.7 mm without IB correction, then 39.7 mm after IB correction. This suggests that the inverse barometer has statistically an annual to inter-annual influence on sea level in Brest. In Marseille, the semi-annual and the annual component seems to be affected, the former having an amplitude of 22.2 mm, then 24.0 mm without then with IB correction, respectively, the latter having 34.4 mm then 44.2 mm amplitude, without and with IB correction, respectively. This suggests a shorter typical timescale of IB influence on Marseille sea level than in Brest, typically at seasonal to annual timescales. IB effect can be positive (see Brest or Marseille 1900–2020 trends in Table 3) or negative (see 1971–2020 trends in Table 3) to linear trends, up to 0.3 mm/yr in these examples.
Triangle diagrams for Brest. Top: sea level trends estimated from linear fits of sea level time series. Middle: sea level trends estimated from quadratic fits of sea level time series. Bottom: sea level accelerations estimated from quadratic fits of sea level time series. Left: uncorrected time series. Center: inverse barometer correction only. Right: inverse barometer and all residual cycles correction. For all nine figures, the minimum length of time series is twenty years, and values exceding the colorbar boundaries are colored as the color boundaries; hatches indicates non-significant regressions (p-value higher than 0.01).
Triangle diagrams for Marseille. Top: sea level trends estimated from linear fits of sea level time series. Middle: sea level trends estimated from quadratic fits of sea level time series. Bottom: sea level accelerations estimated from quadratic fits of sea level time series. Left: uncorrected time series. Center: inverse barometer correction only. Right: inverse barometer and all residual cycles correction. For all nine figures, the minimum length of time series is twenty years, and values exceding the colorbar boundaries are colored as the color boundaries; hatches indicates non-significant regressions (p-value higher than 0.01).
Comparisons of our trends (mm⋅yr−1; mean, standard-deviation) estimate with linear fits for Brest with estimates of Melet et al. (2023)
| This study | Melet et al. (2023) | |||
|---|---|---|---|---|
| Correction | Estimate | |||
| Brest | 1900–2020 | Unfilled, not corrected | 1.5576 | 1.68 ± 0.24 |
| None | 1.5784 ± 0.1749 | |||
| IB only | 1.8019 ± 0.1231 | |||
| IB + all res. cycles | 1.7934 ± 0.0922 | |||
| 1971–2020 | Unfilled, not corrected | 2.9444 | 3.06 ± 0.72 | |
| None | 2.4281 ± 0.6517 | |||
| IB only | 2.1518 ± 0.4515 | |||
| IB + all res. cycles | 2.0879 ± 0.3381 | |||
| 1991–2020 | Unfilled, not corrected | 3.2124 | 3.73 ± 0.63 | |
| None | 3.5014 ± 1.6022 | |||
| IB only | 2.9408 ± 1.1362 | |||
| IB + all res. cycles | 2.8125 ± 0.8508 | |||
| Marseille | 1900–2020 | Unfilled, not corrected | 1.3271 | 1.3–1.4 |
| None | 1.2590 ± 0.1252 | |||
| IB only | 1.3833 ± 0.1068 | |||
| IB + all res. cycles | 1.3863 ± 0.0659 | |||
| 1971–2020 | Unfilled, not corrected | 2.4998 | ||
| None | 2.1540 ± 0.4615 | |||
| IB only | 2.0053 ± 0.3937 | |||
| IB + all res. cycles | 2.0298 ± 0.2429 | |||
| 1991–2020 | Unfilled, not corrected | 3.5198 | 2.5 | |
| None | 3.1552 ± 1.1561 | |||
| IB only | 2.5111 ± 0.9862 | |||
| IB + all res. cycles | 2.5590 ± 0.6084 | |||
For Marseille, Melet et al. (2023) do not provide estimate comparable to the same periods as Brest, but give some more general and less precise values: the value 1.3–1.4 mm⋅yr−1 is given as an estimate since the 19th century and the value of 2.5 mm⋅yr−1 is given since 1993 as an estimate of the whole Mediterranean basin. The “Unfilled, not corrected” trends are estimated from the raw PSMSL data; no uncertainty is provided since it is only formal uncertainties, although we have taken variance from de gap-filling combination step to weight data for the other trends.
3.3. Sea level rise fits
General aspect of all possible linear or quadratic fits seem very similar, when given the degree of fit, with or without any of the above mentioned effects corrected. As expected, long length periods appear to be the most similar, and short length periods reveal more differences (see Figures 2 and 3). The general pattern clearly shows that linear trends are almost every time positive, oscillating in Brest between 1 mm/yr and 3 mm/yr for 30 years period lengths (Figure 4). In Marseille, sea level rise rate almost continuously grows from 0.5 mm/yr at the beginning of 20th century, to about 2 mm/yr at the end of the 1950s. Then sea level trends decreases abruptly to −1 mm/yr around 1975, then recovers to more than 2 mm/yr as of the 1990s. Almost all linear fits are significant at more than 99%; the only exception arises from Marseille, where the visual sea level decrease around the 1970s corresponds to negative linear trends for periods length shorter than about 45 years, and significant at 99% confidence level on a short window centered on year 1970 for periods shorter than 35 years; IB correction provides more significant estimates only when combined with periodic components estimates.
Triangle period lengths-based sections for Brest (left) and Marseille (right) sea level. Top: linear trends for periods lengths of 30 years. Middle: quadratic degree-1 estimates for period lengths of 60 years. Bottom: quadratic degree-2 estimates for period lengths of 60 years.
Quadratic fits provide much more detailed features on trends and accelerations. Nonsignificant fits appear less frequent. In Brest, periods with alternate sea level rise and decrease are only present before the 1960s for period lengths longer than 80 years; afterwards, only period lengths shorter than 40 years, between 1950 and 1965 central dates, then 30 years (around 1980) can sometimes have negative trends; anytime outside these periods of possible negative trends, Brest sea level trends are positive. In Marseille, the two visual sea level decrease events are well represented by negative trends from quadratic fits. The first decrease event is only visible for period lengths shorter than 50 years. The second decrease event appear to have been so strong that it causes negatives trends to be estimated until one century period length around the 1960s. Until 1950, sea level rise appear to have been strong, between 2 mm/yr and 6 mm/yr depending on correction scheme and associated estimated uncertainties, for period lengths longer than 60 years, which corresponds to the beginning of the influence of the large sea level decrease event, and to be again very positive for central dates posterior to 1970 (around 2 mm/yr), albeit less positive than before the large sea level decrease event. At its maximum, the large sea level decrease event reached a negative trend of about −2 mm/yr for periods lengths of 60 years (Figure 4). These finding are in contrast with linear trends estimates, which display more positive trends after the large sea level decrease event than before, although linear trends fit display lower negative trends during the event. Another interesting fact is the time lag of minimum sea level trends when estimated by linear fits or quadratic fits; linear fits minimum is located at the end of the 1970s, while it appears ten years earlier with quadratic fits. This finding would be of a great help for operational monitoring of local sea level rise linked to coastal adaptation policies. As an illustration, the 30 years period length trend minimum during the decrease event is estimated to be equal to −1.4 ± 0.7 mm/yr between 1957 and 1987, while it is equal to −18.6 ± 6.8 mm/yr between 1947 and 1977 with quadratic fit. Our linear trend estimate appear slightly higher in absolute value than those from Calafat et al. (2022) (see also references therein), who find a decrease of −0.3 ± 0.5 mm/yr (−0.6 mm/yr if corrected from GIA to obtain relative sea level, versus −0.8 ± 0.7 mm/yr in our study with linear fit), mainly attributed to a positive anomaly in atmospheric pressure over the Mediterranean. With 60 years period lengths, linear fits no longer see negative trend, with a minimum equal to 0.7 ± 0.3 mm/yr, while quadratic fit estimates this minimum to be equal to −2.4 ± 1.4 mm/an between 1929 and 1988. These elements give confidence in quadratic fits to better display sea level trends and variability than linear fits.
Acceleration patterns appear to be more variable in Brest than in Marseille. Brest sea level acceleration oscillates around zero, but has a positive maximum around 1940, for period lengths of 60 years, then slowly increases (Figure 4). In Marseille, periods posterior to the large sea level decrease event display more positive accelerations than periods prior to this event (between 1 mm/yr/decade and 1.5 mm/yr/decade, versus 0.5 mm/yr/decade); during the decrease event, maximum negative acceleration reached between −0.5 mm/yr/decade and −1 mm/yr/decade. The large sea level decrease event even strongly affects sea level acceleration estimate on the full time span, and is the cause of the negative acceleration estimated on the full time span. On the contrary, strong acceleration estimated on posterior periods following this event could be interpreted as a recovery after an anomalous decrease. All in all, Marseille sea level quadratic fits are almost every time significant ate the 90% level, except periods lengths shorter than 30 years around 1950 in Brest, period lengths shorter than 40 years around 1910 in Marseille, and then period lengths between 35 and 55 years around 1970 corresponding to the end of the Marseille sea level decrease event.
Comparisons of our estimate of Brest accelerations (mm⋅yr−2; mean, standard-deviation, double of degree-2 estimate) estimates with quadratic fits with estimates of Haigh et al. (2014)
| This study | Haigh et al. (2014) | ||
|---|---|---|---|
| Correction | Estimate | ||
| 1900–2009 | Unfilled, not corrected | −0.0016 | −0.0003 ± 0.0096 |
| None | −0.0026 ± 0.0126 | ||
| IB only | −0.0080 ± 0.0090 | ||
| IB + all res. cycles | −0.0078 ± 0.0068 | ||
| 1915–2009 | Unfilled, not corrected | 0.0186 | 0.0200 ± 0.0124 |
| None | 0.0048 ± 0.0170 | ||
| IB only | −0.0080 ± 0.0120 | ||
| IB + all res. cycles | −0.0082 ± 0.0090 | ||
| 1930–2009 | Unfilled, not corrected | 0.0318 | 0.0340 ± 0.0165 |
| None | 0.0160 ± 0.0248 | ||
| IB only | 0.0002 ± 0.0170 | ||
| IB + all res. cycles | 0.0004 ± 0.0132 | ||
| 1960–2009 | Unfilled, not corrected | 0.1178 | 0.1208 ± 0.0520 |
| None | 0.0510 ± 0.0838 | ||
| IB only | 0.0036 ± 0.0594 | ||
| IB + all res. cycles | 0.0066 ± 0.0444 | ||
The “Unfilled, not corrected” accelerations are estimated from the raw PSMSL data; no uncertainty is provided since it is only formal uncertainties, although we have taken variance from de gap-filling combination step to weight data for the other accelerations.
Having explored all possible time series central dates and length, we can chose similar periods to compare our estimates to other studies (see Table 3 for comparison with linear trends from Melet et al. (2023), and Table 4 for comparison with accelerations from Haigh et al. (2014)). Our uncorrected trends estimates are consistent with Melet et al. (2023) for all periods, but appear lower, when corrections are applied. Another comparison comes from SONEL website, where Brest 1960–2018 trend is estimated to 1.80 mm/yr, while we estimated on the same time span 1.71 ± 0.47 mm/yr, making both estimates fully consistent. A greater difference occurs for Marseille, where SONEL 1960–2018 estimate is 1.39 mm/yr, while our uncorrected estimate is 0.84 ± 0.33 mm/yr. The only difference between both datasets is that ours has no data gaps, thanks to the combination model we used to include neighboring tide gauges time series. Since the combination model does not alter trends by construction on the full time span, this result suggests that data gaps in Marseille could be the cause of this difference of 0.56 mm/yr on 1960–2018 period. The sign and the order of magnitude of the difference is close to those from the comparison of our 1971–2020 estimate with Melet et al. (ibid.). Finally we compared our 1885–2022 Marseille uncorrected estimate with NOAA “tide and currents” website estimates, which are fully consistent, with values of 1.27 ± 0.11 mm/yr and 1.31 ± 0.12 mm/yr, respectively. By contrast, our uncorrected Brest sea level accelerations estimates reveal that the signs and and the variations of our estimates are consistent with estimates of Haigh et al. (2014), but that the amplitude is lower and the uncertainty greater (see Table 4). Times series corrected from inverse barometer only of in combination with periodic component have much lower absolute accelerations and delayed emergence of positive acceleration, although still uncertain.
3.4. Correlation between fits parameters and climate indices
The inter decadal variability of trends and accelerations adresses the question of physical mechanism possibly acting as additional processes to climate change driven by anthropogenic greenhouse emissions. Following the method exposed in Section 2, we have been able to find at the same time the climate index as well as the typical timescale maximizing the correlation with sea level trends and accelerations. Since most variability is displayed with quadratic fits, we focus here on degree-1 and degree-2 estimates from quadratic fits. In Brest and Marseille, as well as for trends and accelerations, we find that the Northern Atlantic Oscillation (NAO) correlates the most with these parameters. The typical time scales are between 70 years to one century. An interesting fact arises for Marseille trends and accelerations. At time scales maximizing the correlation with NAO, the large sea level decrease event remains visible. However, if we look at the maximizing time scale for inverse barometer corrected time series trends and accelerations, it appears that both trends and accelerations are non-negligibly less negative during the decrease event, as well as less positive before the decrease event (see Figure 5). This suggests that atmospheric pressure due to NAO multidecadal variability could be at the origin of both strongly positive trends before the decrease event and strongly negative trends during the decrease event, as proposed by several studies (Gomis et al., 2008; Calafat et al., 2022). As an illustration, maximized NAO correlation for uncorrected sea level time series is obtained with 79-yr period length, for which quadratic trend is equal to −1.3 ± 0.9 mm/yr between 1914 and 1992 (central epoch 1953.0); by comparison, with the same central-epoch, the IB corrected time series quadratic fit maximizing NAO correlation with 97-yr period length provides trend equal to −0.1 ± 0.5 mm/yr; for equal (but not maximizing NAO correlation) period lengths, the difference between uncorrected and IB-corrected time series with central epoch 1953.0 is about 0.3 mm/yr, suggesting the same order of magnitude as proposed by Calafat et al. (2022) or Gomis et al. (2008). No such effect is seen in Brest.
Marseille sea level trends from quadratic fits period lengths that maximize correlation with lowpass filtered NAO index.
3.5. Sea level rise extrapolations
After having estimated linear and quadratic fits on all possible time series of Brest and Marseille combined times series, we have selected all periods which ends at the end of the full time span, corresponding to top right inclined side of triangle diagrams, in order to extrapolate sea level at each decade until 2100. First we can see on Figure 6 that linear fits systematically provide sea level extrapolations that are almost every time lower than the IPCC projections, unless for shorter time series lengths which witness of recent higher, albeit more uncertain, sea level rise trends; it is also striking that linear fit extrapolations are even not consistent with the lowest emissions IPCC scenario (SSP1-1.9), and that, for period length longer than 50 or 60 years, their likely confidence level do not even have common confidence interval values. These inconsistencies probably come from the fact that linear fits, by construction, are unable to catch sea level rise acceleration due to recent ocean warming and fresh water contributions to sea level rise.
Brest (left) and Marseille (right) 2050 sea level rise extrapolation with respect to 1995–2014 mean sea level, depending on the interpolation period length for linear and quadratic fits of uncorrected, IB corrected and IB + periodic components corrected time series. IPCC model-based projected sea level projections are also displayed for comparison (horizontal colored lines). Associated likely confidence interval is 17%–83%.
Second, we can see that the feature is different for quadratic fits. Despite substantial uncertainties, we can see that sea level in Brest based on period lengths shorter than about 60 years are all consistent with IPCC projections, which do not diverge at this time horizon; for longer time series, the weight of long linear behavior becomes important and we retrieve the result found with linear fits, that long time series provide sea level extrapolations lower than the lowest emissions IPCC scenario SSP1-1.9. In Marseille, only sea level projections based on IB corrected time series extrapolation with period lengths between 50 and 70 years are consistent with IPCC projections and corresponds also to the highest sea level rise extrapolation: 23.6 [236.8; 194.0] cm (median, 17%–83% CL) in Marseille, with an interpolation length of 64 years, and 191.8 [125.5; 264.2] cm in Brest, with a quadratic interpolation length of 57 years. For periods shorter than this time window, sea level extrapolations are too uncertain to be considered robust. For longer period lengths, sea level extrapolations are lower than IPCC SSP1-1.9 scenario and cannot be considered as representative of future physical processes from models, because future sea level rise dynamics will be driven by too different processes than those which caused past observed elevation. An average period length of 60 years appears to be optimal for both Brest and Marseille for several reasons: first, 60 years corresponds to the period length identified by Douglas (1991) as a minimum for robust trends estimates, but also to the period of sea level cycle identified by Jevrejeva, Moore, Grinsted, A. P. Matthews, et al. (2014); second it is closed to the period length of 50 years used in the USA to extrapolate combined tide gauges data (Sweet et al., 2022); third, 60 years is a compromise between shorter periods, which display high uncertainties but are representative of recent sea level rise acceleration, and longer periods which propagate influences of far past processes in the future which could be not representative of what will really happen. In addition, it has to be noted the difficulty to use older data than the epoch of the large sea level decrease event for long-term estimates, because this drop acted as a temporary counter-forcing cause of climate change (due to dam buildings after Second World War, see Frederikse et al., 2020; Dangendorf et al., 2019) which still influences long-term sea level estimates. The recovery found after sea level drop observed in the 1960s is enhanced by the recent sea level rise acceleration driven by water warming and fresh water contributions.
3.6. Sea level rise milestones timing
As a consequence, the time to which some given levels are reached by sea level depends on both the degree of the interpolation polynomial and the period length chosen for extrapolation (see Supplementary Material Figure 1). Because of their interest for short-term coastal adaptation planning, we have chosen thresholds of 20 cm and 50 cm with respect to 1995–2014 mean sea level as milestones, and calculated epochs to which these levels will be reached by mean sea level (see Figure 7). If we consider the optimal period of 57 years for Brest, the IB-corrected linear fit model predicts that 20 cm sea level rise milestone will be reached in 2123 [2101; 2157], 17%–83%, although the quadratic fit model gives 2050 [2042; 2066]; level 50 cm will be respectively reached in 2309 [2253; 2396] and 2088 [2072; 2121]. These inconsistencies warns about intrinsic uncertainties due to the choice of the model and the interpolation period on which observations can extrapolated; for instance, it is likely that observations linear extrapolation is no longer valid to estimate a 50 cm sea level rise timing. For Marseille, if we retain an optimal period length of 64 years, the 20 cm level is reached in 2234 [2184; 2318] with the linear fit model, but in 2045 [2040; 2050] with the quadratic fit model; level 50 cm is reached during the 26th century with the linear fit model, and in 2077 [2069; 2086] with the quadratic fit model.
Extrapolated timing of 200 mm (left) and 500 mm (right) sea level rise, with respect to 1995–2014, for Brest (top) and Marseille (bottom), based on different interpolation period lengths. Associated likely confidence interval is 17%–83%.
4. Conclusion
On the methodological plan, we have proposed a combination model that proves to be efficient, in the context of tide gauge data, to fill gaps without altering neither the vertical datum of the chosen reference tide gauge nor its time-behaviour. In addition, transformation parameters allow to express observations of an individual tide gauge as virtual observations of the reference tide gauges, here chosen to be Brest and Marseille. Inversely, transformation parameters allow Brest and Marseille tide gauges to be considered as virtual stations of the other tide gauges used to the combination processes. If parameterized differently, the same model could be used to combine or compare other kinds of data, in different contexts. This will be the topic of future studies. Different corrections schemes have been tested, which do not affect the principal features of sea level trends and accelerations that we estimated. However, these corrections reduce uncertainties of estimated parameters, thus can be useful to constraint observations-based sea level extrapolations. The last methodological tool used here is the “triangle diagram”, which had been already used previously for other sea level analysis studies, and has the advantage to avoid any arbitrary choice in terms of time span selection to estimate variations of data. This tool is particularly interesting for a systematic analysis of variables displaying some variability, but it is conditioned to the obligation to have continuous time series; here these continuous times series were built with the help of the above-described combination model used to fill data gaps.
Concerning results, we have shown that quadratic fit are better adapted to represent the large sea level decrease event which happened in Marseille during the second part of 20th century, as well as the strong sea level rise trends displayed before and after this event. Now, Marseille sea level accelerates at rates never seen on the time span of our study, adding now 1 mm/yr every decade to trends already reaching almost 3 mm/yr. Brest sea level rises at smoother paces, with acceleration of about 0.5 mm/yr per decade and trends at about 2.5 mm/yr. We have added suggestion about the role of the Northern Atlantic Oscillation on past major sea level rise or decrease, notably in Marseille.
However, we have also shown that a structural uncertainty on sea level extrapolations comes from the choice of the time span used to interpolate the data and from the degree of the polynomial to fit them. Nevertheless it is also shown that a quadratic fit with the last sixty years of data provide extrapolations up to 2050 consistent with IPCC scenarios which are not well discriminated at this epoch, and that a too-short length of interpolation is too noisy. It is also shown that a too-long period length of interpolation provide underestimated extrapolated projections, because the recent period of strong acceleration seen both in Brest and Marseille comes after a non-negligible period of sea level drop after World War II.
Quadratically extrapolating sea level time series until 2050 shows, first, that time series longer than about 60 years cannot see current sea level rise acceleration and provide underestimated sea level rise projections compared to IPCC emulator projections based on physical processes; it is obviously the same with linear extrapolation whatever the length of the interpolated period. If quadratically extrapolated until 2100, sea level time series can display, despite uncertainties, some consistency with rather high IPCC scenarios for interpolated period length of about 60 years; shorter time series are too subject to variability and uncertainties and they provide very variable and noisy sea level projections, although longer time series do not see sea level rise acceleration and provide sea level rise projections that are underestimated compared to IPCC projections. This could mean first that the tide gauge data extrapolation model needs to be quadratic in order to model properly sea level rise acceleration, second that an optimal length for tide gauge sea level interpolation is about 60 years if extrapolation is needed (close to Douglas, 1991; Jevrejeva, Moore, Grinsted, A. P. Matthews, et al., 2014; Sweet et al., 2022), and third that observed tide gauge sea level rise time series are consistent with high emissions IPCC scenarios (confirming Fox-Kemper et al., 2021, (see their Section 9.6.3.3)).
Finally, the methods presented here could contribute to a sea level rise climate services aiming at monitoring long-term sea level rise trends, accelerations, milestones timing and associated uncertainties, constraining projections with observations, and identifying the most probable scenario based on available observations. We hope this could help practitioners for coastal adaptation planning.
Further studies are however needed to not rely only on observations extrapolations or on models simulations, but to rigorously combine both kinds of data, and to make simulations conditional to observations. Recent published papers developed such an approach for temperature and air humidity with Bayesian formalism (Ribes et al., 2021; Douville et al., 2022; Qasmi and Ribes, 2022) which, for instance, cannot be adapted to sea level rise because of the absence of full coupled sea level rise simulations of past sea level rise. We hope that this study will support efforts in this direction.
Data availability
PSMSL tide gauges data can be freely downloaded from https://psmsl.org. NASA sea level projections for individual tide gauges can be freely downloaded from https://sealevel.nasa.gov/ipcc-ar6-sea-level-projection-tool. Climate indices can be downloaded from the NOAA website: https://psl.noaa.gov/data/ climateindices/list/.
Acknowledgments
All figures were made using Python matplotlib module; all calculations were made using diverse Python Numpy module functions.
Declaration of interests
The author does not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and has declared no affiliations other than his research institution.
Supplementary materials
Supporting information for this article is available on the journal’s website under https://doi.org/10.5802/crgeos.302 or from the author.
Supplementary material file contains:
- periodic components parameters estimated for Brest and Marseille tide gauges over different corrected time series
- trends and accelerations obtained from other corrections schemes than presented in the main article.