\makeatletter
\@ifundefined{HCode}
{\documentclass[screen,CRGEOS,Unicode,biblatex,published]{cedram}
\addbibresource{crgeos20240586.bib}
\newenvironment{noXML}{}{}
\let\citep\parencite
\let\citet\textcite
\def\xcitealp#1#2{\citeauthor{#1}, \citelink{#1}{#2}}
\def\citealt#1#2{\citeauthor{#1}, \citelink{#1}{#2}}
\newcommand*{\citelink}[2]{\hyperlink{cite.\therefsection @#1}{#2}}
\def\defcitealias#1#2{}
\let\citetalias\textcite
\let\citepalias\parencite
\def\thead{\noalign{\relax}\hline}
\def\endthead{\noalign{\relax}\hline}
\def\tbody{\noalign{\relax}}
\def\tabnote#1{\vskip4pt\parbox{.87\linewidth}{#1}}
\def\tsup#1{$^{{#1}}$}
\def\tsub#1{$_{{#1}}$}
\def\ndash{\text{--}}
\def\hyphen{\text{-}}
\def\0{\phantom{0}}
\def\mn{\phantom{$-$}}
\skip\footins=10pt
\makeatletter
\g@addto@macro{\UrlBreaks}{\UrlOrds}
\gappto{\UrlBreaks}{\UrlOrds}
\RequirePackage{etoolbox}
\usepackage{upgreek}
\def\jobid{crgeos20240586}
%\graphicspath{{/tmp/\jobid_figs/web/}}
\graphicspath{{./figures/}}
\def\xmorerows#1#2{#2}
\newcounter{runlevel}
\usepackage{multirow}
\def\nrow#1{\@tempcnta #1\relax%
\advance\@tempcnta by 1\relax%
\xdef\lenrow{\the\@tempcnta}}
\def\morerows#1#2{\nrow{#1}\multirow{\lenrow}{*}{#2}}
\let\MakeYrStrItalic\relax
\def\refinput#1{}
\let\ubreak\break
\csdef{Seqnsplit}{\\}
\def\botline{\\\hline}
\def\xbotline{\hline}
\def\back#1{}
\newcommand*{\hyperlinkcite}[1]{\hyper@link{cite}{cite.#1}}
\DOI{10.5802/crgeos.302}
\datereceived{2024-07-15}
\daterevised{2025-03-21}
\datererevised{2025-07-16}
\dateaccepted{2025-07-16}
\ItHasTeXPublished
}
{
\PassOptionsToPackage{authoryear}{natbib}
\documentclass[crgeos]{article}
\usepackage[T1]{fontenc}
\def\CDRdoi{10.5802/crgeos.302}
\def\citelink#1#2{\citeyear{#1}}
\def\xcitealp#1#2{\citealp{#1}}
\def\xmorerows#1#2{\morerows{#1}{#2}}
\let\ubreak\relax
\def\xbotline{\botline}
\makeatletter
\def\CDRsupplementaryTwotypes#1#2{}
\def\hyperlinkcite#1#2{\citeyear{#1}}
\def\href#1#2{\url[#1]{#2}}
\let\newline\break
}
\makeatother

\usepackage{upgreek}

\dateposted{2025-07-30}
\begin{document}

\begin{noXML}

\CDRsetmeta{articletype}{research-article}

\TopicFR{Sciences du climat}
\TopicEN{Climate Sciences}

\title{Observational uncertainties on past and future sea level rise
for Marseille and Brest tide gauges}

\alttitle{Incertitudes observationnelles sur l'\'el\'evation pass\'ee
et future du niveau de la mer \`a partir des mar\'egraphes de Marseille
et Brest}

\author{\firstname{Jonathan} \lastname{Chenal}\CDRorcid{0009-0006-5565-4046}}
\address{Institut national de l'information g\'eographique et foresti\`ere,
Service de g\'eod\'esie et m\'etrologie, Centre national d'\'etudes spatiales, 
IGN/B\^atiment Euler, 18 Avenue \'Edouard Belin, 31400 Toulouse, France}
\email{jonathan.chenal@ign.fr}

\keywords{\kwd{Sea level}
\kwd{Brest}
\kwd{Marseille}
\kwd{Tide gauge}
\kwd{Trend}
\kwd{Acceleration}}

\altkeywords{\kwd{Niveau de la mer}
\kwd{Brest}
\kwd{Marseille}
\kwd{Mar\'egraphe}
\kwd{Tendance}
\kwd{Acc\'el\'eration}}

\begin{abstract}
Brest and Marseille tide gauges time series between 1885 and 2022 data
gaps are filled with data from neighboring tide gauges thanks to a
combination model that does not alter original trends. Continuous
relative sea level time series obtained from this step are then
corrected from the inverse barometer effect and/or periodic components.
Time series are systematically explored to form all possible sub-time
series, from 20 years length to the full time span, and adjusted with
linear and quadratic fits. It is shown that corrections have various
effect on parameters best estimates, depending on the chosen period,
and that they reduce uncertainties on estimated parameters;
corrections schemes comparison and correlation with climate indices
suggest the influence of the Northern Atlantic Oscillation on sea level
trends and accelerations in both locations. It is shown that linear fit
model have less capabilities to represent sea level variability than
quadratic fit model. Quadratic fit model shows that Marseille sea level
strongly dropped after WWII, with trend of almost 4 mm/yr before this
large sea level decrease event, \tralicstex{--}{${-}$}2 mm/yr during the event, then 3 mm/yr
during the recovery following this one, and than Brest sea level did
not see such variability. It is shown that sea level extrapolation is
consistent with IPCC projections only for interpolation period lengths
around 60 years. On this basis, 20 cm additional sea level rise is
expected until 2050, and level 50 cm will be reached before the end of
21st century in both Brest and Marseille.
\end{abstract}

\begin{altabstract}
Les donn\'ees manquantes des s\'eries temporelles des mar\'egraphes de
Brest et Marseille entre 1885 et 2022 sont combl\'ees avec les
donn\'ees des mar\'egraphes voisins gr\^ace \`a un mod\`ele de
combinaison qui n'alt\`ere par leur tendance initiale. Les s\'eries
continues de niveau de la mer relatif obtenues \`a l'issue de cette
\'etape sont ensuite corrig\'ees de l'effet de barom\`etre inverse
et/ou des composants p\'eriodiques. Les s\'eries temporelles sont
ensuite syst\'ematiquement explor\'ees pour former toutes les
sous-s\'eries temporelles possibles, d'une longueur minimale de 20 ans
jusqu'\`a la dur\'ee enti\`ere des donn\'ees, puis ajust\'ees avec des
mod\`eles lin\'eaires et quadratiques. Il est montr\'e que les
corrections ont des effets vari\'es sur les param\`etres estim\'es,
selon la p\'eriode choisie, et qu'elles r\'eduisent les incertitudes
sur les param\`etres estim\'es ; la comparaison des sch\'emas de
comparaison et la corr\'elation avec les indices climatiques sugg\`ere
l'influence de l'Oscillation de l'Atlantique Nord sur les tendances et
acc\'el\'erations des deux sites. Il est aussi montr\'e que le mod\`ele
d'ajustement lin\'eaire est moins capable de repr\'esenter la
variabilit\'e du niveau de la mer que le mod\`ele quadratique. Le
mod\`ele quadratique montre que le niveau de la mer \`a Marseille a
fortement baiss\'e apr\`es la Deuxi\`eme Guerre Mondiale, avec une
tendance de 4 mm/an avant cet \'ev\'ement de baisse profonde du niveau
de la mer, de \tralicstex{--}{${-}$}2 mm/an pendant l'\'ev\'enement, puis de 3 mm/an
pendant la p\'eriode de rattrapage suivant celui-ci, et que Brest n'a
pas connu une telle variabilit\'e. Il est montr\'e que l'extrapolation
du niveau de la mer est coh\'erente avec les projections du GIEC
uniquement pour des p\'eriodes d'interpolation proches de 60 ans. Sur
cette base, 20 cm d'\'el\'evation du niveau de la mer est attendu d'ici
2050, et le niveau 50 cm sera atteint avant la fin du 21e si\`ecle \`a
Brest comme \`a Marseille. 
\end{altabstract}

\maketitle

\vspace*{9pt}

\twocolumngrid

\end{noXML}

\defcitealias{IPCC_AR6_WG1_ch9Oceancryospheresealevelchange_2021}{ibid.}
\defcitealias{Meletetal_2023}{ibid.}

\section{Introduction}\label{sec1}
Sea level rise is an important consequence of anthropogenic emissions
of greenhouse gases
\citep{IPCC_AR6_WG1_ch9Oceancryospheresealevelchange_2021,Meyssignacetal_2023_NCLIM}. 
The question of sea level rise projections is crucial for coastal
planning  \citep{Nichollsetal_2021} and uncertainties on sea level
projections are a major subject of concern \citep{Hirschfeldetal_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 \citep{LeCozannetetal_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 \citep{Qasmietal_2022,Ribesetal_2021}, in spite of
the availability of sea level rise reconstructions and quadratic fits
from tide gauges and/or space altimetry \citep[e.g.][see also 
\url{https://sealevel.colorado.edu/}]{Dangendorfetal_2019,Neremetal_2018,Neremetal_2022,Guerouetal_2023}. 
Recently the United States Sea
Level Rise Task Force proposed to quadratically extrapolate combined US
tide gauges observations until 2050 \citep{Sweetetal_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
\citep{IPCC_AR6_WG1_ch9Oceancryospheresealevelchange_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, \citet{Visseretal_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
\citep{HolgateWoodworth_2004}, 30 years for GMSL
\citep{Frederiksetal_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
\citep{Jevrejevaetal_2008,Haightetal_2014} or from shorter, but
continuous, time series from space altimetry
\citep{Ablainetal_2019,Guerouetal_2023}. At local scale, sea level
variations can be measured either with space altimetry
\citep{Prandietal_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
\citep{Asghariniaetal_2020,Kandasamyetal_2013} or machine learning
\citet{Arriagadaetal_2021},
and
\citet{Kumaretal_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 \citep{Vergeretal_2013}, which is fundamentally a
close and simplified model to the Helmert seven parameters widely used
in the geodesy subdomain of reference frames
\citep{Altamimietal_2002,Altamimietal_2023}.\looseness=-1

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
\citep{Douglas_1991,Jevrejevaetal_2008,Jevrejevaetal_2014,ChurchWhite_2011,Woodworthetal_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
\citep{Pennaetal_2013,FeatherstoneFilmer_2012}, and to connect vertical
datums of territories separated by the sea
\citep{FilmerFeatherstone_2012}, by the combination of long tide gauges
time series and mean dynamic topography\break 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 \citep{WoppelmannMarcos_2016}. Another effect comes from the
seasonal atmospheric influence on sea level, named  ``inverse
barometer'' \citep{WunschStammer_1997,Woodworthetal_2009}, which some
studies do not correct for in hydrodynamic leveling
\citet{FeatherstoneFilmer_2012}, although some others correct it either
from atmospheric reanalysis \citep{Pennaetal_2013} or by only selecting
summer data \citep{Afrastehetal_2023}. Finally, astronomical cycles
affect mean sea level observations on seasonal to multi-decadal time
scales \citep{Woodworth_2012}, and have also to be corrected for in
hydrodynamic levelling, for instance if required to transform raw
observations into mean tide observations \citep{MakinenIhde_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
\citep{Visseretal_2015,Merrifieldetal_2009}. \citet{Douglas_1991}
recommended for tide gauges time series to be at least 60 years long to
assess accelerations. \citet{Jevrejevaetal_2008}
and \citet{Jevrejevaetal_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 \citet{HoustonDean_2011} did not conclude
to any acceleration on the basis of 20th century tide gauges
observations, but \citet{RahmstorfVermeer_2011} demonstrated that this
conclusion relied only on short time spans or regional analysis by
systematically explore all time series based on
\citet{ChurchWhite_2011} reconstruction. In addition,
\citet{Dangendorfetal_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
\citep{Ezer_2013,Sweetetal_2022}, in Europe \citep{Meletetal_2023}, or
in single countries, such as Portugal
\citep{Antunes_2019,AntunesTaborda_2009}.\looseness=1

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~\ref{sec:datamethods}, then results and discussion
are presented in Section~\ref{sec:resultsdiscussion}.\looseness=-1

\section{Data and methods}\label{sec:datamethods}

\subsection{Data}\label{sec2.1}
We used monthly sea level RLR data from the Permanent Service for Mean
Sea Level \citep{Holgateetal_2013} for tide gauges of Brest (PSMSL 
id.~1), Marseille (id.~61), Genova (id.~59), S\`ete (id.~958), Toulon 
(id.~980), Nice (id.~1468), \mbox{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)  \citep{Chiangetal_2004}, Atlantic
multidecadal oscillation (AMO)  \citep{Enfield_2001}, Arctic
oscillation (AO)  \citep{Higginsetal_2002}, North Atlantic oscillation
(NAO) \citep{Hurrell_1995,Jonesetal_1997}.

In addition, we downloaded IPCC AR6 WG1 regional sea level rise
projections \citep{IPCC_AR6_WG1_ch9Oceancryospheresealevelchange_2021}
based on emulators FAIR \citep{Smithetal_2018} for the thermosteric
component of sea level rise, and FACTS  \citep{Koppetal_2023} for the
mass component of sea level rise, as computed by NASA
\citep{Garneretal_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). 

\subsection{Gap-filling combination model}\label{sec2.2}
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 $N_s$ epochs denoted
$t_{is}$. Let $Y_s$ be the vector of individual values $y_{is}$. There
are a total of $N$ epochs ($t_i$).

Our combination model relies on a transformation model where we
consider two datasets $Y_1$ and $Y_2$ of the same observable variable.
We write the transformation model as an adaptation of the well-known
seven-parameters similitude, here adapted and 
\mbox{reduced} to only one
translation $T$ and one scale factor~$F$:	
{\begin{equation}\label{eq1}
Y_1 = T + F\cdot Y_2
\end{equation}}\unskip
If $F$ is sufficiently close to 1, we can write $F = 1 + D$ and:
{\begin{equation}\label{eq2}
Y_1 = T + Y_2 + D\cdot Y_2
\end{equation}}\unskip
$T$ can be interpreted as a difference of origins between $Y_1$ and
$Y_2$, such as a different reference period for a climate variable
anomaly or an arbitrary origin. $F$ (or equivalently $D$) has a more
interesting physical meaning, since it can reveal drifts between $Y_1$
and $Y_2$, coming from different physical behaviors, but also from
fundamental different physical meanings between $Y_1$ and $Y_2$ (e.g.\ 
the difference between Global mean surface temperature denoted GMST and
Global surface air temperature GSAT
\citep{IPCC_AR6_WG1_ch3humaninfluenceclimate_2021}), different emergent
properties (e.g.\ consequences of different climate sensitivities on
surface temperatures if we compare two climate models outputs
\citep{IPCC_AR6_WG1_ch7energybudget_2021}), or more simply different
physical units or imperfect realizations of a SI unit. In the case of
tide gauges, the scale factor can reveal differences of sea level rise
between different tide gauges, due to oceanic processes or to different
vertical land motion affecting the instrument, or instrumental drifts
between devices.

Here, our aim is to estimate combined variables ($Y_c\equiv Y_1$) and
transformation parameters $T_s$ and $D_s$ between a combined dataset
$c$ and individual datasets $s$. First we linearize the unknown vector
$Y_c = Y_0 + \delta 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}\Seqnsplit
&\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}}\unskip
where $I_N$ is the $N$-dimension identity matrix if $Y_1$ and $Y_2$ do
have the same points (i.e.\ $N=N_1=N_2$), and $\mathbb{1}$ is a
$N$-vector of 1.

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 $H_s$ and has $N$ columns
and $N_s$ rows. $H_s$ has zeros everywhere unless unknown variable
$y_{ic}$ referring to epoch $t_i$ of the combined solution is observed
by solution $s$; in this case, the corresponding element of $H_s$ 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}\Seqnsplit
\mbox{with}\quad
A_s & = & \left[
\begin{array}{cc}
\mathbb{1}_s & Y_{0s}
\end{array}
\right]
\Seqnsplit
\mbox{and}\quad 
\theta_s
& = & 
\left[
\begin{array}{c}
T_s\\
D_s
\end{array}
\right]
\end{eqnarray}}\unskip
The left-hand side vector of  Equation~(\ref{eq:combinationmodel}) is
the observations vector; a variance matrix $\Sigma$ is associated to
this vector, and it is built with individual datasets variance matrices
$\Sigma_s$: 
{\begin{eqnarray}
\Sigma  = 
\left[
\begin{array}{ccc}
\Sigma_1 & & {0} \\
 & \ddots & \\
{0} & & \Sigma_S
\end{array}
\right]
\end{eqnarray}}\unskip
The right-hand side  of  Equation~(\ref{eq:combinationmodel}) is the
product of the design matrix and the unknowns vector.

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}}\unskip
where $K_r = [1,1]^t$ and $K_s = [0,0]^t$ if $s\neq r$; I$_{2S}$ is the
identity matrix of dimension $2S$. These pseudo-observations are
associated to a diagonal variance matrix $\Sigma_K$:
{\begin{eqnarray}
\Sigma_K  =  
\left[ 
\begin{array}{ccc}
\sigma_{K_1}^2 &  & {0} \\
 & \ddots & \\
{0} & & \sigma_{K_S}^2
\end{array}
\right]
\end{eqnarray}}\unskip
Here, $\sigma_{K_r}$ is chosen very small to constraint $\theta_r = 0$,
whereas $\sigma_{K_s}$ for $s\neq r$ is to be relaxed. The conditions
block $K$, the identity matrix $I_{2S}$, the pseudo-observations vector
$\Theta=[\theta_1,\hdots,\theta_S]^t$ and the variance matrix
$\Sigma_K$ adequately complete the above-written equations system.

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\cdot 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\`ete, 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. 

\subsection{Corrections applied to gap-filled tide gauges time
series}\label{sec2.3}
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
\citep{IPCC_AR6_WG1_ch9Oceancryospheresealevelchange_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 \citep{WunschStammer_1997} by using NOAA
20CRv3 \citep{Compoetal_2011} reanalysis until 2015 (included), and
ECMWF ERA5 reanalysis \citep{Hersbachetal_2020} afterwards (2016 and
later). Another effect is the long-period astronomical cycles of period
18.6 years and 8.8 years \citep{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
\citep{Woppelmannetal_2014,WoppelmannMarcos_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
\citet{IPCC_AR6_WG1_ch9Oceancryospheresealevelchange_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. 

\subsection{Time series fits exploration, extrapolation and milestones
timing}\label{sec2.4}
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)\times(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
\citet{IPCC_AR6_WG1_ch9Oceancryospheresealevelchange_2021}. Note also
that the extrapolated sea level is still expressed with respect to the
1995--2014 mean sea level, also as in
\citetalias{IPCC_AR6_WG1_ch9Oceancryospheresealevelchange_2021}.

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.

\subsection{Correlation of sea level trends and\newline accelerations with
climate indices}\label{sec2.5}
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.

\section{Results and discussion}\label{sec:resultsdiscussion}

\subsection{Tide gauges datasets gap fillings}\label{sec3.1}
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\`ete, 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~\ref{fig:TGtimeseries} 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.

\begin{figure*}
\vspace*{3pt}
\includegraphics{fig01}
\vspace*{4pt}
\caption{\label{fig:TGtimeseries}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.}
\vspace*{6pt}
\end{figure*}

We can see on  Table~\ref{tab:transfoMarseille} and 
Table~\ref{tab:transfoBrest} 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~\ref{fig:TGtimeseries}, but that combined time series are, by
construction, corrected for diverse contributions of individual
translations from individual tide gauges origins.

\begin{table*} %t1
\caption{\label{tab:transfoMarseille}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}
\tabcolsep9pt
\begin{tabular}{ccccccc}
\thead
Tide gauge & Time span & Nb. obs. & $T$ & $D$ & RMS & $\sigma_0^2$\\
 & & & std $T$ (1$\sigma$) & std $D$ (1 $\sigma$) && \\
\endthead
Genova & 1885--2022 & 1267 & ${-}$37.5465   &    \mn0.0098 & 24.1580 & 1.09258\\
 & & & \mn\00.9999   &    \mn0.0010\vspace*{4pt} &&\\
Marseille & 1885--2022 & 1551 & \mn\00.0000 & \mn0.0000 & 17.2447 & 0.97086\\
 & & & \mn\00.0010   &    \mn0.0010 \vspace*{4pt}&&\\
S\`ete & 1992--2022 & \0321 & \0${-}$8.6675   &   ${-}$0.0097 & 25.4999 & 0.92233\\
 & & & \mn\01.0000   &    \mn0.0010 \vspace*{4pt}&&\\
Toulon & 1961--2022 & \0435 & \mn18.9891   &    \mn0.0001 & 16.2592 & 0.87035\\
 & & & \mn\00.9999   &    \mn0.0010 \vspace*{4pt}&&\\
Nice & 1978--2022 & \0474 & \mn14.2194     &  \mn0.0038 & 27.9493 & 1.06432\\ 
 & & & \mn\01.0000   &    \mn0.0010 \vspace*{4pt}&&\\
Port-Vendres & 1984--2022 & \0363 & \0${-}$2.2627  &    ${-}$0.0018 & 20.0580 & 0.85746\\
 & & & \mn\01.0000   &    \mn0.0010 \vspace*{4pt}&&\\
L'Estartit & 1990--2022 & \0396 & \mn10.5637   &   \mn0.0048 & 11.8920 & 0.77326\\
 & & & \mn\00.9999   &    \mn0.0010 &&
\botline
\end{tabular}
\vspace*{4pt}
\end{table*}

\begin{table*} %t2
\caption{\label{tab:transfoBrest}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}
\tabcolsep9pt
\begin{tabular}{ccccccc}
\thead
Tide gauge & Time span & Nb. obs. & $T$ & $D$ & RMS & $\sigma_0^2$\\
 & & & std $T$ (1$\sigma$) & std $D$ (1 $\sigma$) & &\\
\endthead 
Brest & 1885--2022 & 1543 & \mn\0\00.0002   &    \mn0.0000 & 22.8671 & 0.71811\\
& & & \mn\0\00.0010    &   \mn0.0010\vspace*{4pt}&&\\
Newlyn & 1915--2022 & 1265 & ${-}$104.4754  &     \mn0.0133 & 24.7139 & 0.74431\\
& & & \mn\0\00.9999   &    \mn0.0010\vspace*{4pt}&&\\
Saint-Nazaire & 1941--2022 & \0465 & \0${-}$71.6149   &    \mn0.0093 & 39.3662 & 0.80143\\
& & & \mn\0\01.0000   &    \mn0.0011\vspace*{4pt}&&\\ 
La Rochelle & 1941--2022 & \0513 & \0${-}$51.9616  &    ${-}$0.0103 & 36.2872 & 0.81203\\
& & & \mn\0\01.0000   &    \mn0.0011\vspace*{4pt}&&\\
Saint-Jean-de-Luz & 1942--2022 & \0541 & \mn\078.8907   &   ${-}$0.0267 & 47.0729 & 0.91168\\
& & & \mn\0\01.0000   &    \mn0.0011\vspace*{4pt}&&\\
Le Havre & 1938--2022 & \0695 & \mn\015.7387   &   ${-}$0.0083 & 39.4882 & 1.15245\\
& & & \mn\0\01.0000  &     \mn0.0010\vspace*{4pt}&&\\
Boulogne-sur-Mer & 1941--2022 & \0427 & \0${-}$15.6859   &   ${-}$0.0067 & 56.1212 & 1.20251\\
& & & \mn\0\01.0000   &    \mn0.0011\vspace*{4pt}&&\\
Calais & 1941--2022 & \mn540 & \mn153.9639   &   ${-}$0.0289 & 67.8163 & 1.20985\\
& & & \mn\0\01.0000   &    \mn0.0011\vspace*{4pt}&&\\
Dunkerque & 1942--2022 & \mn666 & \mn188.6161  &    ${-}$0.0449 & 63.4725 & 1.19975\\
& & & \mn\0\01.0000   &    \mn0.0011 &&
\botline
\end{tabular}
\vspace*{2pt}
\end{table*}

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\`ete 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 
\mbox{sentence,} ``aligned'' means that the final 
combined time series have
null transformations parameters with 
respect to the initial Brest and
Marseille time series.

\subsection{Sea level rise corrections}\label{sec3.2}
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~\ref{fig:trianglesBrest} for Brest and 
Figure~\ref{fig:trianglesMarseille} 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. \mbox{Finally,} \mbox{periodic} \mbox{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 \mbox{Marseille}
1900--2020 trends in  Table~\ref{tab:SLtrendscompMELET}) or negative
(see 1971--2020 trends in  Table~\ref{tab:SLtrendscompMELET}) to linear
trends, up to 0.3~mm/yr in these examples.

\begin{figure*}
\includegraphics{fig02}
\vspace*{4pt}
\caption{\label{fig:trianglesBrest}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).}
\vspace*{4pt}
\end{figure*}

\begin{figure*}
\vspace*{-1pt}
\includegraphics{fig03}
\vspace*{-1pt}
\caption{\label{fig:trianglesMarseille}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).}
\vspace*{-2pt}
\end{figure*}

\begin{table*}%t3
\caption{\label{tab:SLtrendscompMELET}Comparisons of our trends
(mm${\cdot}$yr$^{-1}$; mean, standard-deviation) estimate with linear
fits for Brest with estimates of \citet{Meletetal_2023}\vspace*{-1pt}}
\tabcolsep9.5pt
\begin{tabular}{ccccc}
\thead
 & & \multicolumn{2}{c}{This study} & \xmorerows{1}{\citet{Meletetal_2023}}\\
\cline{3-4}
 & & Correction & Estimate & \\
\endthead
\morerows{11}{Brest} &
\morerows{3}{1900--2020} & 
Unfilled, not corrected & 1.5576 & $1.68\pm0.24$\\
 & & None & $1.5784 \pm 0.1749$ & \\
 & & IB only & $1.8019 \pm 0.1231$ & \\
 & & IB $+$ all res. cycles & $1.7934 \pm 0.0922$ & 
\\\cline{3-5}
 & \morerows{3}{1971--2020} &
Unfilled, not corrected & 2.9444 & $3.06\pm0.72$\\
 & & None & $2.4281 \pm 0.6517$ & \\
 & & IB only & $2.1518 \pm 0.4515$ & \\
 & & IB $+$ all res. cycles & $2.0879 \pm 0.3381$ &
\\\cline{3-5}
 & \morerows{3}{1991--2020} & 
 Unfilled, not corrected & 3.2124 & $3.73\pm0.63$\\
 & & None & $3.5014 \pm 1.6022$ & \\
 & & IB only & $2.9408 \pm 1.1362$ & \\
 & & IB $+$ all res. cycles & $2.8125 \pm 0.8508$ &
\\\xbotline
\morerows{11}{Marseille} & 
\morerows{3}{1900--2020} & 
Unfilled, not corrected & 1.3271 & 1.3--1.4\\
 & & None & $1.2590 \pm 0.1252$ & \\
 & & IB only & $1.3833 \pm 0.1068$ & \\
 & & IB + all res. cycles & $1.3863 \pm 0.0659$ & 
\\\cline{3-5}
 & \morerows{3}{1971--2020} & Unfilled, not corrected & 2.4998 & \\
 & & None & $2.1540 \pm 0.4615$ &  \\
 & & IB only & $2.0053 \pm 0.3937$ & \\
 & & IB + all res. cycles & $2.0298 \pm 0.2429$ & 
\\\cline{3-5}
 & \morerows{3}{1991--2020} & Unfilled, not corrected & 3.5198 & 2.5\\
 & & None & $3.1552 \pm 1.1561$ & \\
 & & IB only & $2.5111 \pm 0.9862$ & \\
 & & IB + all res. cycles & $2.5590 \pm 0.6084$ &
\botline
\end{tabular}
\tabnote{For Marseille, \citet{Meletetal_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${\cdot}$yr$^{-1}$ is given as
an estimate since the 19th century and the value of
2.5~mm${\cdot}$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.}
\vspace*{-3pt}
\end{table*}

\subsection{Sea level rise fits}\label{sec3.3}
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~\ref{fig:trianglesBrest}
and~\ref{fig:trianglesMarseille}). 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~\ref{fig:trianglessections}). 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 
\mbox{periods} length shorter than about 45 years, and 
\mbox{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.

\begin{figure*}
\vspace*{-3pt}
\includegraphics{fig04}
\vspace*{-3pt}
\caption{\label{fig:trianglessections}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.}
\vspace*{-3pt}
\end{figure*}

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 
\mbox{visible} for period lengths shorter than 50 years. 
The~\mbox{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~\ref{fig:trianglessections}). These finding are in
contrast with linear trends estimates, which display more positive
trends after the large sea level decrease event than before, although
\mbox{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 
\mbox{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 \pm 0.7$~mm/yr between 1957 and 1987, while it is equal to 
${-}18.6 \pm 6.8$~mm/yr between 1947 and 1977 with quadratic fit. Our
linear trend estimate appear slightly higher in absolute value than
those from \citet{Calafatetal_2022} (see also references therein), who
find a decrease of  ${-}0.3 \pm 0.5$~mm/yr (${-}$0.6~mm/yr if corrected
from GIA to obtain relative sea level, versus  ${-}0.8 \pm 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 $\pm$ 0.3~mm/yr, while quadratic fit estimates this minimum to
be equal to  ${-}2.4 \pm 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~\ref{fig:trianglessections}).
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 \mbox{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. \looseness=-1

\begin{table*}%t4
\caption{\label{tab:SLBrestaccelerationsHAIGH}Comparisons of our
estimate of Brest accelerations (mm${{\cdot}}$yr$^{-2}$; mean,
standard-deviation, double of degree-2 estimate) estimates with
quadratic fits with estimates of \citet{Haightetal_2014}\vspace*{-1pt}}
\tabcolsep15pt
\begin{tabular}{cccc}
\thead
 & \multicolumn{2}{c}{This study} & \xmorerows{1}{\citet{Haightetal_2014}} \\
\cline{2-3}
 & Correction & Estimate & \\
\endthead
\morerows{3}{1900--2009} & 
Unfilled, not corrected & ${-}$0.0016 & ${-}0.0003\pm0.0096$\\
 & None & ${-}0.0026 \pm 0.0126$ & \\
 & IB only & ${-}0.0080 \pm 0.0090$ & \\
 & IB $+$ all res. cycles & ${-}0.0078 \pm 0.0068$ & 
 \\\cline{2-4}
\morerows{3}{1915--2009} & 
Unfilled, not corrected & 0.0186 & $0.0200\pm0.0124$\\
 & None & $0.0048 \pm 0.0170$ & \\
 & IB only & ${-}0.0080 \pm 0.0120$ & \\
 & IB $+$ all res. cycles & ${-}0.0082 \pm 0.0090$ & 
 \\\cline{2-4}
\morerows{3}{1930--2009} & 
Unfilled, not corrected & 0.0318 & $0.0340\pm0.0165$\\
 & None & $0.0160 \pm 0.0248$ & \\
 & IB only & $0.0002 \pm 0.0170$ & \\
 & IB $+$ all res. cycles & $0.0004 \pm 0.0132$ & 
 \\\cline{2-4}
\morerows{3}{1960--2009} & 
Unfilled, not corrected & 0.1178 & $0.1208\pm0.0520$\\
 & None & $0.0510 \pm 0.0838$ & \\
 & IB only & $0.0036 \pm 0.0594$ & \\
 & IB $+$ all res. cycles & $0.0066 \pm 0.0444$ &
\botline
\end{tabular}
\tabnote{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.}
\vspace*{-2pt}
\end{table*}

Having explored all possible time series central dates and length, we
can chose similar periods to compare our estimates to other studies
(see  Table~\ref{tab:SLtrendscompMELET} for comparison with linear
trends from \citet{Meletetal_2023}, and 
Table~\ref{tab:SLBrestaccelerationsHAIGH} for comparison with
accelerations from  \citet{Haightetal_2014}). Our uncorrected trends
estimates are consistent with \citet{Meletetal_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 \pm
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 \pm 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  \citetalias{Meletetal_2023}. \mbox{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
\pm 0.11$~mm/yr and  $1.31 \pm 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  \citet{Haightetal_2014}, but that the 
\mbox{amplitude} is lower
and the uncertainty greater (see 
Table~\ref{tab:SLBrestaccelerationsHAIGH}). 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.

\subsection{Correlation between fits parameters and\newline 
climate indices}\label{sec3.4}
The inter decadal variability of trends and accelerations adresses the
question of physical mechanism possibly acting as additional processes
to 
\mbox{climate} change driven by anthropogenic greenhouse 
\mbox{emissions.}
Following the method exposed in Section~\ref{sec:datamethods}, 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~\ref{fig:MarseilleCORRMAX}). 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 \citep{Gomisetal_2008,Calafatetal_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 \pm 0.9$~mm/yr between 1914 and 1992 (central epoch
1953.0); by 
\mbox{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 \pm 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
\citet{Calafatetal_2022} or \citet{Gomisetal_2008}. No such effect is
seen in Brest.

\begin{figure}
\includegraphics{fig05}
\caption{\label{fig:MarseilleCORRMAX}Marseille sea level trends from
quadratic fits period lengths that maximize correlation with lowpass
filtered NAO index.}
\end{figure}

\subsection{Sea level rise extrapolations}\label{sec3.5}
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~\ref{fig:2050SL} 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\break level rise.

\begin{figure*}
\includegraphics{fig06}
\caption{\label{fig:2050SL}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\%.}
\vspace*{-3pt}
\end{figure*}

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 \mbox{(median,} \mbox{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 \citet{Douglas_1991} as a minimum for
robust trends estimates, but also to the period of sea level cycle
identified by \citet{Jevrejevaetal_2014}; second it is closed to the
period length of 50 years used in the USA to extrapolate combined tide
gauges data \citep{Sweetetal_2022}; third, 60 years is a compromise
between shorter periods, which display high \mbox{uncertainties} but are
\mbox{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 \citep[due to dam buildings after Second World War, 
see][]{Frederiksetal_2020,Dangendorfetal_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 
\mbox{contributions.}

\subsection{Sea level rise milestones timing}\label{sec3.6}
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~\ref{fig:niv}). 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. 

\begin{figure*}
\includegraphics{fig07}
\vspace*{-1pt}
\caption{\label{fig:niv}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\%.}
\vspace*{-2pt}
\end{figure*}

\section{Conclusion}\label{sec4}
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
\mbox{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
\mbox{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
\citep[close to][]{Douglas_1991,Jevrejevaetal_2014,Sweetetal_2022},
and third that observed tide gauge sea level rise time series are
consistent with high emissions IPCC scenarios 
\citep[confirming][(see their Section~9.6.3.3)]{IPCC_AR6_WG1_ch9Oceancryospheresealevelchange_2021}.

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
\citep{Ribesetal_2021,Douvilleetal_2022,Qasmietal_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.

\section*{Data availability}
PSMSL tide gauges data can be freely downloaded from
\url{https://psmsl.org}. NASA sea level projections for individual tide
gauges can be freely downloaded from
\url{https://sealevel.nasa.gov/ipcc-ar6-sea-level-projection-tool}.
Climate indices can be downloaded from the NOAA website:
\href{https://psl.noaa.gov/data/climateindices/list/}{https://psl.noaa.gov/data/\unskip\break climateindices/list/}.

\section*{Acknowledgments}
All figures were made using Python matplotlib module; all calculations
were made using diverse Python Numpy module functions.

\section*{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.

\section*{Supplementary materials}
Supporting information for this article is available on the journal's
website under \printDOI\ or from the author.

Supplementary material file contains:
\begin{itemize}
\item 
periodic components parameters estimated for Brest and Marseille tide
gauges over different corrected time series
\item 
trends and accelerations obtained from other corrections schemes than
presented in the main article.

\end{itemize}

\CDRsupplementaryTwotypes{supplementary-material}{\cdrattach{crgeos-302-suppl.pdf}}

\back{}

\printbibliography
\refinput{crgeos20240586-reference.tex}

\end{document}
