Nathalie Giloy
email: nathalie.giloy@shom.fr

Alexa Latapy

Gaël André

Estimation of skew surge uncertainties when predicting tides in the past

Keywords:

The study of storm surges and more generally of extreme events is a subject of growing interest in the current context of coastal risk management in a changing climate. Recent studies show that the consideration of past events improves statistical models. In this context, important investigations have been carried out in the last decade to retrieve and quantify water levels and skew surges associated with historical events for the 18th or 19th century. The estimation of the historical surges depends on the calculated tidal prediction, and therefore directly on the tidal components and the mean sea level of the period under consideration.

Still, as tidal constituents and mean sea level change over time, it is difficult to estimate these parameters for historical events when sea level data are missing or incomplete. Two cases are considered for calculating tidal predictions. Either short time series are available contemporary to the event of interest, allowing a harmonic analysis and subsequent calculation of tidal constituents and mean sea levels. Or there are no or insufficient contemporary data available for the event of interest and current tidal characteristics are used in combination with a mean sea level adjustment. The present study aims to answer the following questions: depending on the data available for a historical event, which tidal constituents and mean sea level should be used when considering events that occurred several decades or even centuries ago? To what extent does the length of the observation series used in the harmonic analysis influence the tidal predictions? How does a mean sea level correction affect tidal predictions made with current tidal components? Finally, what are the uncertainties associated with the estimation of historical extreme skew surges, when there are not enough contemporary observations to estimate tidal components?

Two methods are considered to estimate the uncertainties associated with the calculation of tidal predictions. Either the use of contemporary tidal constituents and mean sea level or the use of current tidal constituents corrected for mean sea level. In this study, 14 sea level records of 100 years or more provided by the Global Extreme Sea Level Analysis (GESLA version 3) were analysed. First, the uncertainties of the astronomical tides are assessed. For this purpose, different lengths of observations are implemented in a harmonic analysis and the resulting tidal constituents and tidal predictions are compared. In a second step, current tidal constituents are used for past tidal predictions and a mean sea level trend is inferred to cope with sea level variations. Different time periods are used to estimate the mean sea level trend and their influence on the tidal predictions is analysed. Furthermore, these two methods are combined to study the effects on the estimated skew surges and to assess the associated uncertainties.

This study shows that uncertainties of less than 10 cm can be obtained for skew surges by using contemporary tidal components calculated from at least two years of observations. Larger differences are obtained for microtidal or river-influenced sites. When using current tidal constituents and adjusting a mean sea level trend, linear trend correction gives reliable results when more than 60 years of data are used for the trend estimation. In this case, the extreme skew surges obtained are within a range of more or less 5 cm. Finally, on the basis of the results presented, some suggestions are made on the possible applications of these methods according to the availability of data, stressing the importance of site-specific studies.

Introduction

In low-lying coastal areas, storms and their associated storm surges are natural hazards, that, depending on their severity, affect coastal communities, infrastructures and industries. Storm surges are the oceanic response to meteo-oceanic forcing associated with tropical or extratropical cyclones [START_REF] Pugh | Sea-Level Science: Understanding Tides, Surges, Tsunamis and Mean Sea-Level Changes[END_REF][START_REF] Pugh | Tides, Surges and Mean Sea Level[END_REF]. The height of the surge varies between geographical regions, with the highest surges reaching up 9 m in North America during Katrina in 2005 or in Southeast Asia, during Cyclone Bhola in the Gulf of Bengal region [START_REF] Needham | A review of tropical cyclone-generated storm surges: global data sources, observations, and impacts[END_REF]. On a smaller scale, Europe is also affected by storm surges. During the 1953 storm, the surge rose up to 3 m [START_REF] Rossiter | The North Sea Storm Surge of 31 January and 1 February 1953[END_REF] causing extensive flooding in the Netherlands, Belgium and the United Kingdom [START_REF] Wolf | Modelling waves and surges during the 1953 storm[END_REF].

As sea level rises due to climate change, the risk associated with storm surges will increase (IPCC AR6 [START_REF] Arias | Climate Change 2021: the Physical Science[END_REF]). In fact, smaller surges superimposed on a higher mean sea level will reach extreme levels more quickly. Storm surges can be defined as two different parameters 1) non-tidal residuals (NTR), i.e., the difference between the observed and predicted water level at each time step [START_REF] Pugh | Tides, Surges and Mean Sea Level[END_REF]; 2) skew surge as the difference between the highest observed and predicted water level per tidal cycle [START_REF] Haigh | Spatial and temporal analysis of extreme sea level and storm surge events around the coastline of the UK[END_REF]. Both surge definitions are used in the context of coastal flooding. The NTR provide information on the temporal evolution of the storm as it is continuous in time, but is highly dependent on the tides [START_REF] Horsburgh | Tide-surge interaction and its role in the distribution of surge residuals in the North Sea[END_REF]. On the other hand, the skew surges provide a more reliable measure of the meteorological effects on water levels [START_REF] Batstone | A UK best-practice approach for extreme sea-level analysis along complex topographic coastlines[END_REF][START_REF] Williams | Tide and skew surge independence: new insights for flood risk[END_REF] but only provide one surge component per tidal cycle. Nevertheless, [START_REF] Santamaria-Aguilar | Are extreme skew surges independent of high water levels in a mixed semidiurnal tidal regime?[END_REF] showed that there can be a dependence between skew surges and high-water levels in mixed semidiurnal tidal at sites on the shallow continental shelf. In the context of coastal flooding, the skew surge parameter appears to be more suitable than the NTR for use in risk assessment frameworks [START_REF] Feng | Comparison between the skew surge and residual water level along the coastline of China[END_REF].

Tide gauge observations, called systematic data, are used. Often, these data cover only short periods. For example, in the GESLA-3 dataset [START_REF] Haigh | GESLA Version 3: a major update to the global higher-frequency sea-level dataset[END_REF], more than 65% of the stations have less than 20 years of observations. These short durations are not sufficient to establish high return periods (e.g. 100 years) with acceptable uncertainties. Indeed, there is a general agreement that extrapolations should not be made for return periods longer than four times the length of the record [START_REF] Pugh | Sea-Level Science: Understanding Tides, Surges, Tsunamis and Mean Sea-Level Changes[END_REF]. In the assessment of flood risk, according to the European Directive 2007/60/EC, high return periods, reaching 100 or even 1 000 years, are required in public and industrial and coastal management policies. In order to anticipate future changes, it is important to better assess current risks and investigate past events. When historical events are included in extreme value analysis, they strengthen the upper tail of statistical distributions and reduce uncertainties especially when considering high return periods [START_REF] Gaál | Inclusion of historical information in flood frequency analysis using a Bayesian MCMC technique: a case study for the power dam Orlík, Czech Republic[END_REF][START_REF] Hosking | Paleoflood hydrology and flood frequency analysis[END_REF]). In the present paper, historical events are old events, for which no or very few systematic tide gauge data are available, or events that occurred before the tide gauge era.

For example, the Dunkirk tide gauge is located about 20 km away from the Gravelines nuclear power plant and its data are therefore used in the French nuclear safety assessment (IRSN, 2019). Systematic numerical tide gauge data are available since 1956. Historical events can be found either in systematic quantified sources (e.g., tidal ledgers or marigrams) or in written qualitative historical sources such as journals, local press, diaries etc. A data archaeology project (see [START_REF] Latapy | Data rescue process in the context of sea level reconstructions: an overview of the methodology, lessons learned, up-to-date best practices and recommendations[END_REF] for more information on the methodology of data archaeology) has been carried out and has enabled to add historical sea level observations made at Dunkirk in the second half of the 19th century and early 20th century to be added [START_REF] Latapy | Influence des modifications morphologiques de l'avant-côte sur l'hydrodynamisme et l'évolution du littoral des Hauts-de-France depuis le XIXe siècle. Littoral[END_REF]. However, these records are not continuous and contain gaps of several years or even decades. To complete the dataset, [START_REF] Giloy | Quantifying historic skew surges: an example for the Dunkirk Area, France[END_REF] suggest an approach that allows the reconstruction of historical water levels (usually during storms) based on written qualitative sources. Working with written qualitative sources is complex, as each data has to be considered individually and in its editorial context. [START_REF] Athimon | Quantification of historical skew surges: challenges and methods[END_REF] propose a new method to estimate the reliability of the sources. Thus, for major events such as the storm of 1953, or 1897 or for older events (1720,1763,1767,1778 or 1791, see [START_REF] Hamdi | Analysis of the risk associated with coastal flooding hazards: a new historical extreme storm surges dataset for Dunkirk, France[END_REF] and [START_REF] Giloy | Quantifying historic skew surges: an example for the Dunkirk Area, France[END_REF] for more details) information on water levels reached during extreme events can be found in archives [START_REF] Giloy | Quantifying historic skew surges: an example for the Dunkirk Area, France[END_REF][START_REF] Gorgeu | The rebuilding of the eastern dike of the harbor of dunkirk[END_REF]. Once these water levels have been estimated, tidal predictions must be calculated to deduce the skew surges. It is noteworthy, that when reconstructing historical water levels from archives, only the extreme high water level can be estimated which implies the use of skew surges. The question is how to compute reliable tidal predictions (with accurate tidal constituents (TC) and mean sea level (MSL)) for events where poor or no tide gauge data are available.

Tides are the ocean's response to astronomical forcing. They are a well-known phenomenon that has been studied since ancient times [START_REF] Egbert | Tidal prediction[END_REF]. While the physics behind the tides is well understood, several investigations have shown that the amplitude and phase of TC do not remain constant especially when focusing on secular or decadal timescales [START_REF] Bij De Vaate | Secular trends in global tides derived from satellite radar altimetry[END_REF]. Some of the many explanations include changes associated with large-scale changes in the ocean stratification [START_REF] Colosi | Tales of the venerable Honolulu tide gauge[END_REF][START_REF] Ray | Nineteenth-century tides in the Gulf of Maine and implications for secular trends[END_REF], sea level rise [START_REF] Idier | Sea-level rise impacts on the tides of the European Shelf[END_REF], changes in the coastal geomorphology [START_REF] Bowen | The tidal régime of the River Thames; long-term trends and their possible causes[END_REF], and also changes in the water depth near the harbours due to natural causes such as erosion or sediment accumulation or to human interventions (harbour expansion, channel deepening etc.) [START_REF] Talke | Changing tides: the role of natural and anthropogenic factors[END_REF].

Tides oscillate around a MSL. Still, like TC, these MSL are not stable over time, and changes are observed from seasonal to multi-decadal and secular time scales. These long-term variations follow linear and nonlinear long-term trends, that have been extensively studied [START_REF] Cazenave | Present-day sea level change: observations and causes[END_REF][START_REF] Church | A 20th century acceleration in global sea-level rise[END_REF][START_REF] Haigh | The tides they are A-changin': a comprehensive review of past and future nonastronomical changes in tides, their driving mechanisms, and future implications[END_REF][START_REF] Holgate | On the decadal rates of sea level change during the twentieth century[END_REF]. At the global scale, tide gauge observations show a rise of 21 cm from 1900 to 2020 [START_REF] Palmer | An ensemble approach to quantify global mean sea-level rise over the 20th century from tide gauge reconstructions[END_REF] while at the regional scale, several factors can affect the local MSL such as the contribution of meltwater, vertical land movement [START_REF] Wöppelmann | Vertical land motion as a key to understanding sea level change and variability[END_REF] short-term extreme events or long-term decadal (multidecadal) variations (e.g., ENSO, NAO, etc.) [START_REF] Zhang | Sea level trends, interannual and decadal variability in the Pacific Ocean[END_REF], and thus inducing a non-uniform rise across the globe. Furthermore, a significant acceleration of the global sea level rise (SLR) since the 1960s, fitting a quadratic curve, is now well established up to 0.07 mm. yr -2 in the subtropical North Atlantic [START_REF] Dangendorf | Persistent acceleration in global sea-level rise since the 1960s[END_REF][START_REF] Haigh | Timescales for detecting a significant acceleration in sea level rise[END_REF].

Variations in TC and MSL parameters have an impact on the tidal prediction at a given time, making predictions over historical periods more complex, especially if observations are not available for that period. The interannual variations of the MSL influence the tidal analysis and should therefore be taken into account to correct the estimated surges. Generally, tidal predictions and MSL components have been removed in studies dealing with the estimation of surges. Several methods can be used to estimate these MSL variations such as the annual median [START_REF] Calafat | The sources of seasea-level changes in the Mediterranean Sea since 1960[END_REF][START_REF] Calafat | Probabilistic reanalysis of storm surge extremes in Europe[END_REF][START_REF] Menéndez | Changes in extreme high water levels based on a quasi-global tide-gauge data set[END_REF], the annual mean sea level [START_REF] Haigh | Timescales for detecting a significant acceleration in sea level rise[END_REF][START_REF] Reinert | Seasonal shift in storm surges at brest revealed by extreme value analysis[END_REF] or a Butterworth low-pass filter of order 2 with a cut-off period of 1 or 2 years [START_REF] Marcos | Long-term variations in global sea level extremes[END_REF][START_REF] Marcos | Spatiotemporal changes in extreme sea levels along the coasts of the North Atlantic and the Gulf of Mexico[END_REF]. [START_REF] Wahl | Evidence for multidecadal variability in US extreme sea level records[END_REF] subtract the summer and winter medians (the median is preferred to the mean as it is more robust to outliers) for each year to account for MSL rise and intra/interannual variability.

When tide gauge data are available, one year is considered sufficient to compute stable TC and accurate tidal predictions for navigational purposes [START_REF] Simon | Coastal Tides[END_REF]. The use of contemporary tide gauge data allows the estimation of TC and MSL that represent the local historical hydrodynamics. This method is hereafter referred as method 1.

If contemporary tide gauge data are not available, current TC are used and an additional MSL correction is applied to the tidal predictions, to account for changes in MSL. This correction consists of estimating a linear regression on annual MSL. This method is an adaptation of the methods used by [START_REF] Bardet | Regional frequency analysis of extreme storm surges along the French coast[END_REF][START_REF] Frau | The use of historical information for regional frequency analysis of extreme skew surge[END_REF]. Both authors applied a linear MSL correction to observed tide gauge data in order to fit tidal predictions made using current TC and stable MSL for the whole period. This method is hereafter referred to as method 2.

The aim of this paper is to quantify the uncertainties associated with the method used to calculate tidal predictions for past events when little or no tide gauge data is available. It focuses, on the astronomical tide on the one hand and, the MSL on the other. The paper is structured as follows. First, the tide gauge records and the harmonic analysis software used for this study are described in section 2. This is followed by an estimation and comparison of uncertainties in tidal predictions associated with harmonic analysis and mean sea level trends. The magnitude of the estimated uncertainties is shown in section 3. These results are discussed in section 4 and finally, conclusions and recommendations are presented in section 5.

Data and methodology

GESLA tide gauge dataset

High-quality sea-level measurements are required for many oceanographic and coastal applications. This paper uses the recently released GESLA-3 data set, the extending global data set of higher-frequency (at least hourly) sea-level records from tide gauges [START_REF] Haigh | GESLA Version 3: a major update to the global higher-frequency sea-level dataset[END_REF]. The third version of the data set contains twice as many years of data (91, 021), and nearly four times as many records (5 119) as GESLA-2 [START_REF] Woodworth | Towards a global higher-frequency sea level dataset[END_REF].

As the aim of this study is to examine past surges, it was decided to focus on data sets with a long temporal coverage, i.e., approximately 100 years (or more) of data (Fig. 1, Table 1). 14 stations without abrupt or short time scale vertical land movements were selected. The selected stations for this study are located in different geographical areas (Atlantic Ocean, Pacific Ocean, North Sea, Mediterranean Sea, Baltic Sea), with different tidal types (diurnal, semidiurnal and mixed) and covering all tidal regimes ranging from macrotidal (>4 m), mesotidal (2-4 m) and microtidal (<2 m). Tidal ranges are calculated as the difference between the highest and lowest astronomical tides from tidal constituents estimated over the last 19 years.

An attempt has been made to cover the whole globe but unfortunately the longest series are mainly present in the northern hemisphere. As the data sets have variable sampling from one station to another and sometimes even for the same station over time, they have been interpolated to hourly water levels to simplify some further processes.

As this paper focuses on the study of skew surges, the observed highwater levels are calculated for each station over the entire observation period. To do this, water levels 3 h before/after the predicted high tide are interpolated with a 1-minute time step using a cubic spline, and the time and height of the highest water levels of each tidal cycle are extracted, corresponding to a local maximum of the predicted height.

Software for harmonic analysis

Prior to this study, a comparison was made between two software packages for harmonic analysis (HA) and tidal prediction. The first is MAS, which is the software used by Shom, the French National Hydrographic Office, to calculate official tidal predictions for maritime navigation in areas under French jurisdiction [START_REF] Simon | Coastal Tides[END_REF]. MAS is compared with UTide, an open-source toolbox developed for Matlab and Python to perform HA and tidal predictions [START_REF] Codiga | Unified Tidal Analysis and Prediction Using the UTide Matlab Functions 60[END_REF] which is widely used in different types of applications and has previously been validated by comparison with observations (Bij de [START_REF] Bij De Vaate | Secular trends in global tides derived from satellite radar altimetry[END_REF][START_REF] Tiggeloven | Exploring deep learning capabilities for surge predictions in coastal areas[END_REF].

MAS tidal software

This least-squares based software solves the harmonic equation of the tide by providing the amplitude and phase lag for each TC at the site of interest from a time series of measured water levels. The TC are obtained for each species analysed separately by solving a system of normal equations with reduced vectors using the least-squares method. The system of normal equations must be well conditioned to ensure that the answer will not be indeterminate, by applying a tidal harmonic constituent separation criterion. The number of TC solved by the software is 143, ranging from long periods (annual to bimonthly) to twelfthdiurnal constituents. The MSL is estimated by linear regression over the lunar daily mean, with the computed MSL being the regression value for the last observation. If the duration of sea level time series is too short, nodal TC will not be resolved. To cope for this, MAS applies nodal correction factors to the amplitude and phase lag function of the ratio of the tide generating potential and the phase differences between the reference and the inferred constituents when calculating tidal predictions.

UTide tidal software

Unified Tidal Analysis and Prediction (UTide) uses the uncertainty calculation principles of T_Tide [START_REF] Pawlowicz | Classical tidal harmonic analysis including error estimates in MATLAB using T_TIDE[END_REF] and r_T_Tide [START_REF] Leffler | Enhancing tidal harmonic analysis: robust (hybrid L1/L2) solutions[END_REF] based on the Monte Carlo uncertainty method. N. Giloy et al. This program has been specifically designed to deal with irregular data sets with many gaps. It has many options and allows to perform a robust harmonic analysis taking into account the nodal corrections for a data series longer than 2 years. UTide allows the estimation of 95% confidence intervals for harmonic constants. This measure is derived from linearised error propagation of the total residual power using the detided signal [START_REF] Bij De Vaate | Secular trends in global tides derived from satellite radar altimetry[END_REF][START_REF] Codiga | Unified Tidal Analysis and Prediction Using the UTide Matlab Functions 60[END_REF].

Comparison between the tidal software packages

The comparison between MAS and UTide software is made by confronting the TC issued from a one-year rolling HA over 19-year window. The UTide python package and the MAS tidal software are used selecting the same set of 143 TC. A threshold of 85% of available data is used to perform the analysis, i.e. HA is performed every year on 19-year time series with a maximum gap tolerance of 34.2 months. The resulting amplitudes and phase lags of eight TC (M2, S2, N2, K2, Q1, O1, P1, K1) computed with the MAS and UTide software are compared using statistical values as Root-Mean-Square Error (RMSE) and the correlation coefficient (R). The calculated validation values are given in Table 2 for the TC amplitudes. The results of the comparison between MAS and UTide show the ability of these two software to compute almost the same TC' characteristics.

The RMSE values for the amplitude for M2 do not exceed 0.07 cm, except at Saint-Nazaire and Delfzijl where it reaches 0.11 cm. The maximum differences are observed at Delfzijl with RMSE of 0.08 cm and 0.11 cm for the diurnal components O1 and P1 respectively. The correlation coefficients are generally significant with values higher than 90%, except for some diurnal constituents at Baltimore, Delfzijl, San Francisco and Seattle, where the amplitudes are weak. Similar results are obtained for the phases and are shown in Appendix Table A7.

In addition, for each station, hourly tidal predictions were calculated using both software packages. Using TC estimated over a recent 19-year period. These predictions were then compared to the available observations for each tide gauge. Both software show RMSE in the same order of magnitude with a difference of about 1 cm, except at Key West tide gauge where RMSE difference reaches 4 cm. These results can be found in Table A8.

As the comparison between MAS and UTide shows very good results, the following analysis is performed using the MAS software.

Methodology

Influence of the length of the observations for the harmonic analysis

HA can be carried out over very short periods such as two weeks or over long periods, such as several decades. The longer the time period, the more TC can be identified in the tidal signal. In general, in the literature, HA is performed over one or two years of observations [START_REF] Haigh | Spatial and temporal analysis of extreme sea level and storm surge events around the coastline of the UK[END_REF][START_REF] Menéndez | Changes in extreme high water levels based on a quasi-global tide-gauge data set[END_REF][START_REF] Pineau-Guillou | Large-scale changes of the semidiurnal tide along North Atlantic coasts from 1846 to 2018[END_REF], but if the year under consideration is stormier than others or if climatic conditions induce unusual MSL and TC, this could potentially bias the associated predictions [START_REF] Simon | Coastal Tides[END_REF].

The ideal case would be to use more than 18.6 years of observations for tidal analysis to properly integrate the nodal cycle [START_REF] Pugh | Sea-Level Science: Understanding Tides, Surges, Tsunamis and Mean Sea-Level Changes[END_REF]. Unfortunately, when going back in time, fewer tide gauge data are available. In order to see how timespan affects tidal prediction, 9 different tidal constituent data sets were used, with different lengths of tide gauge observations: 1, 2, 3, 4, 5, 6, 7, 9 and years. A summary of the analysis and tests carried out in this study is given in Table 3. For each station data set, the last year in which more than 85% of the data are available is considered as the 1-year data set. The previous year(s), all containing more than 85% of data, were then added for each range considered, to avoid seasonal modulations affecting the amplitudes. For example, for Baltimore, the most recent full year is 2020, so the 2020 water level observations are used to calculate the 1-year TC (TC1). For TC2, observations from 2019 to are used, for TC3 from 2018 to 2020 and so on up to 19 years of measurements (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020). Using these different data sets, HW and LW levels are estimated for the timespan of each time series. To be able to compare these HW and LW levels taking into account only the influence of TC, the MSL considered for each TC data set are normalised to the MSL of the TC19 data set. Finally, all predictions made with the different TC data sets are compared to the reference data set (TC19).

In addition to the variability of the TC obtained with different time periods, the long-term variability of the TC observed in previous studies and its influence on the resulting tidal predictions are highlighted. For this purpose, the one-year rolling datasets computed over a 19-year window are used, as for the comparison of MAS and UTide (cf. section 2.2.3).

Mean sea level trends

Calculating tidal predictions in the past involves taking into account changes in MSL, which are likely to be significantly different from the present the further back in time you go. In the absence of sea level records for a particular date, it can still be inferred from sea level trends.

Hourly sea level observations are used to calculate daily MSL, with a gaussian filter applied to reduce tidal effects. The Demerliac filter is used at Shom. It filters the tidal signal with a vector of 71 symmetric elements (see [START_REF] Simon | Coastal Tides[END_REF] for more details). Monthly MSL are then calculated by taking the arithmetic mean of at least 15 daily means for each month. Annual MSL correspond to the arithmetic means of monthly MSL weighted by the number of days actually observed during the month, if at least 11 monthly values are available [START_REF] Holgate | New data systems and products at the permanent service for mean Sea Level[END_REF]. Annual means were estimated for calendar years from January to December.

As the annual variability of the MSL affects the standard errors associated with the estimated trends, longer records are required to calculate trends. When analysing data sets from different parts of the world, [START_REF] Douglas | Global sea level rise[END_REF] recommends using 60-year records to smooth out decadal variability. [START_REF] Haigh | Mean sea level trends around the English Channel over the 20th century and their wider context[END_REF] and others have shown that and 50 years of records are required to obtain standard errors of the order of 0.5 and 0.3 mm. yr -1 , respectively.

To estimate the influence of mean sea level trends on tidal predictions, linear trends were estimated for the three different periods (Table 3). First, the period from 1990 to the present is chosen because it corresponds to the period for which satellite sea level observations are available. A second trend is estimated for a 60-year period starting in 1960, as this corresponds to the requirement of [START_REF] Douglas | Global sea level rise[END_REF]. Finally, a secular trend was estimated using observations from 1900 (or the beginning of the station) to the present.

In addition to the linear trend, quadratic coefficients are estimated. Here, a second order polynomial is estimated for the 14 sea level records for the three periods (1900 (or the beginning of the station) -present, 1960-present, 1990-present) to identify sea level acceleration [START_REF] Dangendorf | Persistent acceleration in global sea-level rise since the 1960s[END_REF][START_REF] Haigh | Timescales for detecting a significant acceleration in sea level rise[END_REF]. It should be noted that the estimation of linear or polynomial trends implies the smoothing of interannual variability.

Comparison of both methods

The aim of this study is to obtain uncertainties when calculating tidal predictions for the past and to compare method 1 and method 2. As a reminder, method 1 uses tide gauge data contemporary to a historic event for TC and MSL estimation to calculate tidal predictions. Method is used when no (or not enough) tide gauge data are available for an extreme historic event. It uses the TC from the HA on recent tide gauge data and applies an adjustment to the MSL to account for changes in MSL. To do this, the 10 highest annual skew surges dates are calculated for each station using the original tidal records. Two hypothetical cases are then considered (see Table 3), according to the two methods presented above.

Table 2

Statistical values of the comparison between the tidal constituents' (TC) amplitude results provided by MAS and UTide softwares. Respectively, root-mean-square error (RMSE) and the correlation coefficient (R). Comparison of surges estimated using method 1 and method 2

Surges estimated using method 1: HA on 1, 2, 5 and 19 years of observation contemporary to the event Fig. 6 Surges estimated using method 2: current TC and MSL adjustment using 60 years for the trend estimation + correction of annual MSL variations 1. In this method (method 1), one-year rolling HA are performed on the whole data sets for different following periods: 1, 2, 5 and 19 years. Tidal predictions are made for the 10 reference skew surges of each year using the different TC datasets, centred on the year of the event.

This enables taking into account the hydrodynamics of the event, e. g., contemporary MSL and TC. 2. In this method (method 2), the skew surges are inferred from predictions made with recent TC (estimated from observations from 2000 to the present). The tidal predictions are corrected for the sea level rise using a linear trend (calculated over the period of 1960-present. The resulting linear regression obtained is then applied to the MSL until 1890, before which the MSL is considered stable [START_REF] Wahl | Observed mean sea level changes around the North Sea coastline from 1800 to present[END_REF].

Correlation coefficients and RMSE are then calculated between the annual mean differences of the 10 annual maximum skew surges obtained by method 1 and 2, and the detrended MSL in order to assess the relationship between these annual mean differences and the interannual variability of the MSL. Next, to assess the uncertainties in estimating past skew surges, tidal predictions obtained with both methods 1 and 2 are used, and resulting skew surges are compared by calculating the RMSE and the 95th percentile (P95) of differences. Finally, skew surge differences between the two methods are corrected with the detrended MSL in order to assess whether incorporating the interannual variation of the MSL can reduce the uncertainties in estimating skew surges.

In the absence of bias and considering a Gaussian distribution of errors, the RMSE can be considered as equivalent to the 1-sigma uncertainty, representing 68% of the values.

Results

Uncertainties associated with harmonic analysis

Influence of length of tide gauge time series

Fig. 2 shows the synthesis for all RMSE between the TC19 data set (HA on 19 years) and the other TC data sets (HA on 1-9 years) for HW and LW levels respectively for all stations, calculated according to the method described in section 2.3.1. Detailed figures for each station can be found in the Supplementary Materials. For these 14 stations, with the exception of Delfzijl, a comparison of the TC1 and TC19 data sets shows a RMSE of less than 0.06 m and 0.08 m for HW and LW respectively.

Generally, using more than one year of observations for TC estimation reduces the RMSE. The largest differences are observed at the Delfzijl where the comparison of TC19 and TC1 data sets shows a RMSE of 0.09 m and 0.12 m, respectively, for HW and LW, which is reduced to 0.03 m for HW and LW when comparing TC19 and TC9. All other stations have smaller range differences. For example, the comparison of TC19-TC1 and TC19-TC9 gives differences of about 0.04-0.06 m for Baltimore, Newlyn and Saint-Nazaire to 0.02 m for Sydney. For most stations the RMSE stabilises between 4 and 6 years of observation at about 0.01-0.03 m, except for Delfzijl where the RSME is still at about 0.05 m.

Long-term changes in tidal constituents

The variability of the TC and thus of the associated tidal predictions is not only due to the observation period used for the HA. Temporal variations of amplitudes and phases are visible at all stations using a 19year periods. Fig. 3 shows the evolution for M2 amplitude and phase at the selected stations. The results for the eight selected TC for all stations are available in the Supplementary Materials. By using a 19-year moving HA, the nodal cycle is taken into account, and only the long-term variability of the TC can be analysed. While most of the stations show slow and long-term changes in amplitude (about 2 cm or less) and phase (less than 10 • ) for M2, some stations show larger variations.

For Delfzijl the amplitude of the M2 semi-diurnal constituent increases by more than 15 cm around the 1970s-1980s (Fig. 3). The same break is observed for the S2 and N2 time series with an increase of 3 cm and to a lesser extent for K2 (see Supplementary Materials). At Saint-Nazaire, a sudden change in the amplitudes is observed. For example, for M2, S2 and N2, the amplitudes between the end of the 19th and the beginning of the 20th century are about 170 cm, 60 cm and 34 cm respectively and oscillate around these mean values. An increase in amplitudes is observed, until the 1970s-1980s after which the amplitudes stabilise over time (Fig. 2 and Supplementary Materials). At San Francisco and Honolulu, a steady long-term increase in amplitudes is observed.

Although the amplitudes at Kungsholmsfort are very small, a change can be seen before the 1960s where all amplitudes show a decrease and start to increase until the 1960s. At Marseille, the tide gauge record contains many gaps but amplitudes and phases are stable at the beginning of the period, followed by a generalised decrease in amplitude and phase from the 1950s to the 1980s.

In terms of phase evolution, the largest changes are observed at Kungsholmsfort, where variations reach more than 50 • for M2 and S2 Fig. 2. RMSE in meters between High Waters (HW) and Low Waters (LW) estimated using TC19 data set and the other TC data sets.

with an abrupt change in the 1960s. Honolulu shows an increase for all 8 selected phases with a pronounced increase in the 1940-1960s for M2 and S2 of more than 10 • and for O1, P1 and K1 of more than 5 • . At Key West a decrease of each phase of about 10 • is observed for the semidiurnal TC and about 5 • for the diurnal TC.

Fig. 4 resumes the variations for the amplitudes and phases of the 8 main tidal constituents. Typically, two stations stand out in terms of amplitude: Delfzijl and Saint-Nazaire. Indeed, on the 8 waves plotted, significant variations are observed on N2, M2 and S2, for Delfzijl and Saint-Nazaire. For the phase, variations are observed for Gedser (on Q1, P1 and K2), for Honolulu (on semi-diurnal components) and for Kungsholmsfort for which a strong dispersion is seen for all diurnal and semidiurnal constituents. Like Kungsholmsfort, Marseille seems to stand out for almost all waves, with phase variations between 5 and 10 • over the whole period. In summary, of the 14 stations studied, Delfzijl, Gedser, Kungsholmsfort, Marseille, Saint-Nazaire and Honolulu show significant changes in the main tidal constituents over time.

Uncertainties associated to the mean sea level

Fig. 5 shows the annual MSL for the 14 stations. The MSL have been adjusted for graphical purposes. Table 4 shows the estimated relative linear sea level trends and accelerations for the three different periods. The smallest linear trends and associated uncertainties are observed for the secular trends. Kungsholmsfort is the only station showing no trend for the period 1900-present, all other stations show trends close to or greater than 1 mm. yr -1 reaching even more than 3 mm. yr -1 in Baltimore. In contrast, Boston, Honolulu, Marseille, San Francisco and Seattle show almost stable trends when comparing the 1900-present with the 1960-present trend. Significant increases for this latter period can be seen in Baltimore, Brest, Delfzijl, Key West, Kungsholmsfort Saint-Nazaire and Sydney. The largest increase is recorded at the Baltimore tide gauge with an increase of more than 1 mm. yr -1 .

Considering only 3 decades to estimate a trend, the highest linear increase and uncertainties with rates reaching 4 mm. yr -1 or more are observed at Baltimore, Boston, Key West, Marseille and Newlyn. The smallest increase is observed at Kungsholmsfort with about 1 mm. yr -1 . All stations have the highest uncertainties for the shortest time period with a maximum error of 0.62 mm. yr -1 , while the uncertainties for a secular trend are lower than 0.07 mm. yr -1 .

The acceleration estimated with a 2nd order polynomial trend over the three different time periods is shown on the right-hand part of Table 4. The acceleration corresponds to twice the quadratic coefficient. Baltimore, Delfzijl, Honolulu, Key West and Sydney show increasing accelerations as the timespan is reduced. Boston, Brest, Gedser, Marseille, Newlyn and Saint-Nazaire show an increase in acceleration when comparing the period 1900-present with the period 1960-present but then a decrease and a negative acceleration for the last period. Finally, Kungsholmsfort, San Francisco and Seattle show a decrease in acceleration between the secular period and the 60-year period and then an increase in acceleration that is even higher than the secular acceleration.

Uncertainties associated with the two methods

Skew surge derived from historical TC depending on the time period used to perform HA (method 1) are compared to the skew surge estimation using recent TC with MSL adjustment using a linear MSL trend over 60 years (method 2).

Where only 1 year of observations is used to estimate skew surges, the RMSE for the 14 stations are between 6 and 18 cm (10 and 35 cm when considering the P95 of the skew surges differences) (Fig. 6 a-b). Adding years of records for the HA results in a reduction of the RMSE and thus of the uncertainty. With 5 years of observations, the RMSE are lower than 6 cm for all stations except Delfzijl, Saint-Nazaire and Kungsholmsfort with RMSE of 8.7, 6.3 and 6.8 cm respectively. Considering the P95 of skew surge differences, the error in the estimated surge is in the order of tens of centimetres (7-13 cm) with a maximum of 17 cm for Delfzijl. With 19 years of records, the RMSE decrease further and are between 2 and 7 cm (3-13 cm for P95).

To assess the influence of the interannual variability of the MSL on the surge estimates, correlation coefficients are calculated between the mean annual surge differences (between surges derived from method 1 and 2) and the detrended MSL (Fig. 6 c). The mean annual surge differences and the detrended MSL are better correlated with 2 years of records and less correlated with 19 years of records. With 2 years of records, satisfying coefficients above 80% are obtained, for 10 stations (Fig. 6 c). The RMSE calculated between the mean annual surge differences and the normalised detrended annual MSL (normalised by centring the data to have a mean of zero) are shown in Fig. 6 d. These values are generally relatively low (<0.1 m). Delfzijl stands out from the others with a relatively low correlation coefficients remaining below 51% at most and the highest RMSE for the 4 time periods used. To a smaller degree, Gedser, Kungsholmsfort, Newlyn and Saint-Nazaire show coefficients below 80% with RMSE comparable to the others.

In a second step, the RMSE and the P95 are recalculated, this time correcting the surge differences with normalised interannual mean sea level variations. Generally, the MSL correction reduces the RMSE by 1-2 cm and the P95 by 1-4 cm. Even for stations with correlation coefficients below 60% (Delfzijl, Kungsholmsfort or Newlyn), the MSL corrections lead to a reduction of the RMSE and P95. However, for 19 years of records, it does not improve the surges differences, showing an increase in the RMSE and the P95.

Table 4

Linear mean sea level trends and accelerations for the period 1900-present, 1960-present and 1990- The use of method 1 to compute tidal predictions introduces uncertainties related to the TC and the observation period used for HA. The longer the observation period, the smaller the uncertainties and the higher the number of TC resolved. It should be noted that the MSL estimated in the HA may not be representative if only few years of contemporary records, as interannual variations may not be smoothed enough.

The use of method 2 introduces uncertainties related to the MSL trend applied. For example, using only 30 years of annual means to estimate a MSL trend may not be sufficient to estimate a MSL state for the beginning of the century.

Still, comparing method 1 and method 2, shows that the RMSE are inferior to 20 cm even when using just one year for HA (RMSE of about 0.05-0.18 cm when using 1 year for HA, decreasing to 0.03-0.08 cm when using 19 years for HA). Applying a correction for annual MSL variations further reduced the errors, when using less than 5 years for HA. Highest errors are seen for stations where TC are not constant over time (Delfzijl, Saint-Nazaire) or where the tidal range is small (Kungsholmsfort) (Fig. 6).

Discussion

The aim of this study is to estimate uncertainties when calculating tidal predictions for past events, especially when little or no tide gauge data is available. In most cases, the available data for a historical event do not allow a HA over a 19-year period. Thus, as shown in the results presented in section 3.1 the uncertainties in the tidal predictions increase as the period considered for the HA decreases. It is therefore important to take these uncertainties into account. On the other hand, punctual observations or reconstructions of water levels for a specific event may be available but do not allow the calculation of the TC and MSL contemporary with the event. Therefore, tidal predictions can be computed using TC calculated at another time. Here it is proposed to use TC estimated at the most recent period over a 19-year period and to apply a MSL correction to cope for sea level changes. The tidal predictions are sensitive to the TC variations (section 3.1) and the evolution of the MSL (section 3.2). These two aspects and the associated uncertainties are discussed in this section.

Tidal characteristics as a function of tide gauge observation length

Significant changes in the TC characteristics for a given coastal location can cause changes in the predicted tidal signal. The use of several years of tide gauge data to determine the TC and MSL induces an averaging of variabilities. For example, a stormy year, which may give more noisy results during the HA, will be smoothed out and have less impact than if only one year was considered.

One might ask, why tides in some locations and regions are particularly sensitive to changing boundary conditions and environmental factors, while tides in other locations have remained relatively stable. In their review, [START_REF] Haigh | The tides they are A-changin': a comprehensive review of past and future nonastronomical changes in tides, their driving mechanisms, and future implications[END_REF] highlight that analysis of long tidal series shows, that 1) any changes in the tide gauge location or instrumentation can lead to discontinuities in tides at the local scale and 2) changes in the hydrological regime, morphological changes (bathymetry, sedimentation rates, channel dredging) close to the tide gauge can also affect tides. For example, [START_REF] Ferret | Reconstruction de la série marégraphique de Saint-Nazaire[END_REF] showed that the changes in amplitude and phase lag of TC observed at Saint-Nazaire in the 1980s were due to a change in the tide gauge's stilling well. Most of the recent changes in tidal characteristics have been observed in estuaries and may be caused by changes in dissipation and turbulent mixing, changes in channel depth, changes in surface area, resonance and reflection, and river flow [START_REF] Haigh | The tides they are A-changin': a comprehensive review of past and future nonastronomical changes in tides, their driving mechanisms, and future implications[END_REF][START_REF] Talke | Changing tides: the role of natural and anthropogenic factors[END_REF]Winterwerp et al., 2013;Winterwerp and Wang, 2013). At coastal stations, shifting locations of amphidromic points may also play a role [START_REF] Haigh | The tides they are A-changin': a comprehensive review of past and future nonastronomical changes in tides, their driving mechanisms, and future implications[END_REF]. For each station, site-specific investigations can be carried out to determine the reasons for these changes.

To show the effect on extreme skew surges, the 99th percentile (P99) and the date and value of the maximum skew surges are calculated for the 14 stations for each dataset of skew surges obtained with 1, 2, 5 or 19 years of records for HA using method 1. Table 5 shows the direct impact of the different datasets used to calculate the skew surges on the maximum skew surge and the P99 threshold commonly used in extreme value analysis. For most stations, this does not imply significant changes in the P99, differences are less than 10 cm. Considering only the results obtained with more than 1 year of records for HA, the differences are less than 7 cm. In contrast, this has a much greater impact on the maximum recorded surge for each dataset, with differences of 18 cm and 16 cm for Brest and Seattle and up to more than 20 cm for Delfzijl and Kungsholmsfort, where a maximum of 27 cm is reached. In addition, for Baltimore, Delfzijl, Kungsholmsfort, Marseille, Newlyn, Saint-Nazaire and Sydney, the maximum skew surge does not occur on the same dates. The differences in the dates of the maximum skew surges may be related to the fact that depending on the number of years considered for the HA, some years could not be resolved. Thus, it is more difficult to have 19 years of continuous time series (even with the 85% tolerance threshold) than 2 or 5 years of records. Observations could have been discarded in this way.

It should be mentioned, that in this study the traditional HA was used (as it is the most common method), which is sensitive to instrumental and environmental noise. For example, for microtidal tides, with low amplitudes of the TC, the HA performed is less efficient in resolving phases, which can lead to higher variations such as those obtained at Gedser, Marseille and Kungsholmsfort presented in section 3.1 (Fig. 4). [START_REF] Dangendorf | North sea storminess from a novel storm surge record since AD 1843[END_REF] highlight that for tide gauges in shallow continental shelf sea locations (such as for meso-or macro tidal areas along the German Bight), non-linear effects can bias the harmonic representation of tides by generating additional interaction constituents. Delfzijl is located at the mouth of the Ems estuary and in the Wadden Sea, classical harmonic analysis may be less accurate for this station. More generally, for microtidal sites and/or stations located in estuaries or near estuaries, the use of traditional harmonic analysis could lead to greater uncertainties and variability in the skew surge calculation.

Except for microtidal or river-influenced sites, when more than one year of observations are used to calculate the skew surges, the maximum surges are of the same order of magnitude (<10 cm of difference).

Changes in mean sea level

Method 2, which implies the use of TC estimated from recent tide gauge data, requires a correction of the MSL trend to account for the evolution of the MSL for estimates of skew surges in the past.

Linear MSL trends are commonly used to quantify sea level changes at both global and regional scales, despite the possible presence of nonlinear trends [START_REF] Church | Sea-level rise from the late 19th to the early 21st century[END_REF][START_REF] Frederikse | The causes of sealevel rise since 1900[END_REF][START_REF] Jevrejeva | Nonlinear trends and multiyear cycles in sea level records[END_REF][START_REF] Steffelbauer | Evidence of regional sea-level rise acceleration for the North Sea[END_REF][START_REF] Woodworth | Trends in UK mean sea level revisited[END_REF]. As presented in section 2.3.2, three linear trends are derived from tidal records with different timespans (1900-present, 1960-present, 1990-present) to estimate their influence on the estimated surge and the associated uncertainties.

With the method 2, MSL are estimated in the past using a linear trend; any non-linear and abrupt changes (such as rapid ground movements) will introduce uncertainties in the historical MSL estimation. A site-specific analysis of the tide gauge data and further comparison, socalled "buddy checking" with neighbouring stations [START_REF] Pugh | Sea-Level Science: Understanding Tides, Surges, Tsunamis and Mean Sea-Level Changes[END_REF], is a preliminary step when working with tide gauge series to eliminate abrupt or temporary sea level variations. Quick releases in tensions, such as earthquakes, may cause rapid vertical movements which will subsequently have an incidence on mean sea levels, as for example in the Vanuatu Islands [START_REF] Ballu | Comparing the role of absolute sea-level rise and vertical tectonic motions in coastal flooding, Torres Islands (Vanuatu)[END_REF] or in the Mayotte Island [START_REF] Jeanson | Suivi des impacts de la subsidence, et de la montée du niveau marin relatif, à Mayotte[END_REF]. Anthropogenic patterns of MSL change have also been identified. Massive loading of buildings or massive groundwater extraction are two phenomena that cause land subsidence [START_REF] Esteban | Adaptation to sea level rise: learning from present examples of land subsidence[END_REF]. The latter has been studied for urban areas of cities such as Tokyo, Osaka, Jakarta, Manila, Venice or Trieste [START_REF] Abidin | Land subsidence of Jakarta (Indonesia) and its relation with urban development[END_REF][START_REF] Rodolfo | Global sea-level rise is recognised, but flooding from anthropogenic land subsidence is ignored around northern Manila Bay, Philippines[END_REF][START_REF] Sato | Land subsidence and groundwater management in Tokyo[END_REF]. In fact, the estimation of a MSL trend that includes anthropogenic or geologic temporal pattern may lead to an erroneous estimation of the mean sea level trend, which would then be reflected in tidal predictions.

It has been shown that uncertainties associated to linear trends are reduced when increasing the number of annual MSL used for its estimation [START_REF] Douglas | Global sea level rise[END_REF][START_REF] Tsimplis | Collection and analysis of monthly mean Sea Level data in the mediterranean and the black sea[END_REF]. However, it should be mentioned, that only 5% of the series in the GESLA-3 dataset are Table 5 99th Percentile (P99), Date and maximum skew surge estimated using a one year rolling HA on 1, 2, 5 or 19 years of observation. longer than 60 years [START_REF] Haigh | GESLA Version 3: a major update to the global higher-frequency sea-level dataset[END_REF]. Table 4 shows the estimated relative sea level trends for the three different time periods. For all stations, except Seattle, higher sea level rates and uncertainties are observed when reducing the timespan, with an increase of the calculated trends ranging from 1 mm. yr -1 to even more than 2 mm. yr -1 as a maximum at Saint-Nazaire.

It is clear that using only recent years for trend estimation may not be representative of a secular trend. At Saint-Nazaire, estimating trends using the last three decades gives a trend more than twice as large as the one estimated between 1900 and 2020 (Table 4). Using these trends to calculate skew surges for more remote periods may lead to important overestimation of surges especially when they are distant in time. Further analyses were carried out to estimate a 2nd order polynomial trend over a 30-, 60-and 120-years. The results are not as clear as for the linear trends, no clear acceleration can be seen when timespan is reduced. It should be noted that, as stated by [START_REF] Haigh | Timescales for detecting a significant acceleration in sea level rise[END_REF], the quadratic coefficient is highly dependent on the timespan to which it is fitted and is influenced by the interannual to multidecadal variability of the MSL.

As in section 4.1, the P99 and the maximum skew surge were estimated using the most recent TC (period of 2000-2022) and their MSL were adjusted using the 30-, 60-and 120-year linear trend as well as a 2nd order polynomial trend estimated over the period 1900-present (Table 6). At all stations, the P99 decreases as the number of annual means used to estimate a linear trend is increasing. Most stations show a decrease of about 10 cm or even less, while for Marseille, San Francisco and Seattle it remains stable. When analysing the P99 of each skew surge series, 4 stations have differences of more than 10 cm, with a maximum of 20 cm at Saint-Nazaire. However, if only 60-or 120-year linear trend is used, the difference is around or less than 5 cm. When considering the maximum skew surge, it is interesting to note that the differences at all stations and between all methods are less than 13 cm, except at Saint-Nazaire where the difference reaches almost 30 cm when comparing the surges estimated with the predictions made with a 30-or a 120-year trend.

For almost half of the stations, the date of the maximum skew surge changes when the different trends are applied. For Delfzijl, Newlyn and Saint-Nazaire the highest surge observed using a 30-year trend MSL correction is on a different date from the other three. Kungsholmsfort has three different dates. Effects are seen on the maximum skew surge estimated for the whole series. Simply applying a quadratic fit to the annual mean sea levels does not give clear results. The use of long linear mean sea level trends to correct tidal predictions may be sufficient for estimating skew surges in expertise, provided that the associated uncertainties are also taken into account in extreme value analysis (EVA).

Method 2 is required when historical records are not available. Still, as shown in Fig. 6, when the skew surges obtained by method 1 are corrected by the interannual variability of the MSL, the differences between the two methods are mostly reduced (except for 19 years). The question arises whether it is possible to correct and reconstruct the multiannual sea level variations.

Many studies show that MSL are a major driver of changes in extreme sea levels at interannual and longer timescales [START_REF] Lowe | Past and Future Changes in Extreme Sea Levels and Waves[END_REF][START_REF] Marcos | Spatiotemporal changes in extreme sea levels along the coasts of the North Atlantic and the Gulf of Mexico[END_REF][START_REF] Menéndez | Changes in extreme high water levels based on a quasi-global tide-gauge data set[END_REF][START_REF] Rohmer | Dominance of the mean sea level in the high-percentile sea levels time evolution with respect to large-scale climate variability: a Bayesian statistical approach[END_REF][START_REF] Woodworth | Evidence for centurytimescale acceleration in mean sea levels and for recent changes in Extreme Sea levels[END_REF][START_REF] Woodworth | Forcing factors affecting Sea Level changes at the coast[END_REF]. Climate indices have been widely studied in the recent decades to explain MSL variations. Regional sea levels are governed by climate modes such as the North Atlantic Oscillation (NAO), the Mediterranean Oscillation index in Europe, or the Pacific Oscillation and El Niño Southern Index [START_REF] Agha Karimi | Sea level variation around Australia and its relation to climate indices[END_REF][START_REF] Becker | Sea level variations at tropical Pacific islands since 1950[END_REF][START_REF] Tsimplis | The forcing of mean sea level variability around Europe[END_REF]. But the influence of the NAO on mean sea level variability has been shown to be variable in time and space [START_REF] Jevrejeva | Influence of large-scale atmospheric circulation on European sea level: results based on the wavelet transform method[END_REF][START_REF] Tsimplis | Towards a vulnerability assessment of the UK and northern European coasts: the role of regional climate variability[END_REF]. A possibility could be to use regional sea level indices [START_REF] Haigh | Mean sea level trends around the English Channel over the 20th century and their wider context[END_REF][START_REF] Wahl | Observed mean sea level changes around the North Sea coastline from 1800 to present[END_REF][START_REF] Woodworth | Trends in UK mean sea level revisited[END_REF] or long-term annual mean sea level reconstructions [START_REF] Ray | Experiments in reconstructing twentieth-century sea levels[END_REF], if they exist, to correct MSL using method 2 as they include the interannual MSL variability for specific regions. Still these indices cover only the timespan of the longest tide gauge series available. For extreme events that occurred before the tide gauge era, these interannual variations can not be taken into account. In addition, on a local scale, anthropogenic influences (harbour changes, urbanisation, development of coastal or estuarine areas) may also contribute to local sea level variations [START_REF] Talke | Changing tides: the role of natural and anthropogenic factors[END_REF] but these anthropogenic effects are still difficult to identify, to quantify and to apply without site-specific studies.

Influence of harmonic analysis and mean sea level on skew surge estimations

It has been shown that the estimation of tidal predictions is sensitive to the length of the observations used for the HA, as well as to an adjustment of MSL, when necessary. In the context of making tidal predictions for extreme historical events, in most cases tide gauge data are not available and method 2 must be used. In addition, even when sufficient data are available to make a HA (i.e. 1 or 2 years), there remains an uncertainty in the interannual variations in MSL, as these cannot be estimated.

To resume, for method 1, the consideration of using only one year of observations to estimate TC allows the calculation of tidal predictions but the uncertainties associated with the TC remain around 3-12 cm (Fig. 2). This is due to the fact that, long-period TC may not be detected or accurately determined, even with a threshold of 85% of available data used for HA. In addition, an annual MSL is estimated which may not be representative of the MSL influenced by interannual variabilities and meteo-oceanic conditions of the year analysed. Applying HA to 2 years of observations allows to reduce the uncertainties associated with TC, almost all tidal constituents can be resolved but still, the MSL may not be well estimated as it is influenced by interannual variations. With 5 years of observations used for HA, the skew surge uncertainty is less than 4 cm for the 14 stations studied in this paper (excluding Delfzijl) and is about 2 cm with 19 years of records. When the timeseries in HA are not long enough, the MAS software applies a parametric correction of the nodal factor. This is estimated using a bulk formulation based on the tidegenerating potential. Still, these parametric corrections are not suitable for some locations [START_REF] Feng | Nodal variations and long-term changes in the main tides on the coasts of China[END_REF][START_REF] Müller | Rapid change in semi-diurnal tides in the North Atlantic since 1980[END_REF][START_REF] Ray | Nineteenth-century tides in the Gulf of Maine and implications for secular trends[END_REF][START_REF] Woodworth | A survey of recent changes in the main components of the ocean tide[END_REF]. Therefore, the use of 19 years of records for HA allows to resolve all tidal constituents, including the nodal components, and provides a MSL contemporary to an event where interannual variations are smoothed. This explains why the comparison between method 1 and 2 with both 19 years of records for HA gives a lower RMSE and why the annual MSL adjustment does not improve the results. Table 5 shows that for extreme values, the effect on the 1% extreme surge remains less than 10 cm, but the maximum skew surge of each series can differ for more than 10 cm.

Unfortunately, long periods of tide gauge observations are not always available in the past, so the alternative is to use recent tide gauge data for HA over a 19-year period and apply a modified MSL following a trend (method 2). It should be noted that the application of a linear (or a quadratic) trend will only work if no short-term vertical land movements have been identified for the tide gauge station. According to [START_REF] Douglas | Global sea level rise[END_REF], the results show that 30 years may not be sufficient for a reliable assessment of a secular trend, especially in the current context where rapid changes in MSL are observed. Table 6 shows that the P99 surge is more sensitive to changes in MSL with differences of more than 10 cm being observed when using a trend estimated over only 30 years.

A comparison of method 2 with method 1 (using 19 years of observations for HA) shows good results with an RMSE lower than 6 cm and a P95 lower than or equal to 10 cm when Delfzijl and Kungsholmsfort are excluded. This shows that method 2 gives results within an acceptable range of uncertainties, which is important to know especially when using these values for further analysis. There may be a dependency between skew surges and high-water levels in mixed semi-diurnal tidal at sites located in shallow continental shelf [START_REF] Santamaria-Aguilar | Are extreme skew surges independent of high water levels in a mixed semidiurnal tidal regime?[END_REF]. As this study only focuses on the comparison of skew surges provided by both methods, this dependency does not affect the results presented, but should be considered when performing extreme value analysis (EVA).

Conclusion

In order to better assess and manage the risks associated with coastal flooding, there is a growing interest in historical events. The design levels of defences are often calculated using statistical analysis to estimate centennial or millennial return periods.

The aim of this study is to investigate the methodological uncertainties in the estimation of tidal predictions, especially for past events where no or poor tidal observations are available for the estimation of TC and MSL. Hourly observations of 14 long tide gauge series obtained from the GESLA-3 data set are used in this study. Selected stations were chosen according to their length and their relative stability, i.e. they are little or not affected by abrupt or discontinuous vertical movements. Firstly, a comparison is made between two HA software packages, MAS and UTide, with the Shom internal software MAS being used throughout the study. An analysis was then made of the number of years of observation used for HA, ranging from 1 year to 19 years. As the number of years used for HA increases, the differences with the TC obtained from 19 years of tide gauge observations decrease. The RMSE of the predicted water level decrease from 9 cm to 11 cm for HW and LW respectively, when 1 year is considered to 3 cm when 9 years are used. The evolution of 8 main TC over the observation period is then presented. Small changes (less than 5 cm in amplitudes and 10 • in phases) are observed at most of the stations, but can reach 15 cm and 50 • for a few TC. The estimation of MSL trends and their impact on past tidal predictions are then investigated, showing that 30 years may not be sufficient to represent secular sea level changes.

Focusing more on the extremes, the influence of the length of observations for HA (one-year rolling HA on 1, 2, 5 or 19 years of observations) is compared to tidal predictions made with current TC and an adjustment for different mean sea level trends (using 30, 60 or 120 years for estimating a linear trend and 120 years for estimating a 2nd order polynomial trend). Both methods show that the date of the maximum skew surge may change and differences in the maximum skew surge may reach more than 10 cm. Still, the use of the yearly rolling HA at different time periods also strongly influences P99 at some stations with maximum differences reaching more than 20 cm.

The two methods are then compared and show results that are within the same order of magnitude. Except for Delfzijl in an estuarine area and Kungsholmsfort with very small tidal ranges, the RMSE between the two methods are smaller than 15 cm with a maximum P95 of 25 cm, when using only one year of tidal observations for the TC estimation. Increasing the number of observational years for HA reduces the RMSE smaller than 10 cm for 2 years and about 6 cm for 5 years. It can also be noted that correcting skew surges by taking into account the interannual variations of the MSL reduces these values. However, the knowledge of the annual MSL is not always available, especially for historical events.

Based on these results, the following suggestions are made for the estimation of extreme skew surges in the past.

1. In the case where observations are available around the event, the authors recommend using the longest observation period available for a HA to estimate the TC and MSL. With at least 2 years of observations the average uncertainties reduces from 10 cm to 7 cm (compared to the used of 1 year). If 19 years are available the nodal components can be properly estimated, overwise they have to be inferred by parametric formulations. It should be noted, that uncertainties associated with interannual sea level variations remain, especially when the records used for HA are shorter than 5 years. Using 5 years of observations around the event gives skew surges with an uncertainty of about 5 cm.

2. If no or only few records are available, the current TC may be used (preferably estimated on a full nodal cycle (18.6 years)) for tidal predictions in the past and to apply a mean sea level trend correction to cope with sea level rise. Prior to this, an analysis of the station must be made, to check whether it is subject to short time-scale vertical land movements (e.g. earthquakes or subsidence) and whether it is subject to major anthropogenic changes (e.g. shoreline changes or major dredging). Short time-scale vertical land movements will affect the estimated mean sea level trend, and anthropogenic changes will affect the tidal characteristics. It is also suggested to use at least 60 years of observations for the trend estimations, which is more representative for a long-term change and allows to reduce uncertainties of the trends smaller than 0.2 mm. yr -1 . Linear trend estimates are sufficient as polynomial trends are highly dependent on the start and end dates of the analysis.

As many global stations may not have at least 60 years of sea level observations, mean sea level trends can be calculated for the available period and used to correct tidal predictions in the past. However, these trends may not be representative of a secular trend. Uncertainties associated with the skew surges might be higher (uncertainties of 9-10 cm) and could overestimate extreme skew surges (with variations of up to 30 cm).

These assumptions and associated uncertainties must be taken into account in EVA. Today, it is possible to add ranges of skew surges or lower bounds for the statistical estimation of high return periods [START_REF] Bulteau | How historical information can improve estimation and prediction of extreme coastal water levels: application to the Xynthia event at La Rochelle (France)[END_REF]. Another possibility would be to work with historical water levels and systematic skew surges such as proposed by Saint [START_REF] Saint Criq | Extreme Sea level estimation combining systematic observed skew surges and historical record sea levels[END_REF]. This study proposes a first approach to quantify the uncertainties related to tidal predictions and subsequently to skew surges. Further developments and analysis could be made considering more stations or using HA adapted to different coastal environments. In another sense, this study also raises questions about the calculation of an appropriate tidal prediction, i.e., whether the seasonal variations, modelled by the long period constituents as annual and semi-annual waves, should be taken into account in the tide predictions or at least in a method qualified as a reference. The methods tested here all have advantages and disadvantages, but above all they show the strong sensitivity of the surges to the method of calculating the tidal prediction. Furthermore, the uncertainty associated with the measurement of water levels and their instruments has not been considered in this article and would add to the estimated uncertainties.

Fig. 1 .

 1 Fig. 1. A) location of tide gauges used in this study; b) duration of the sea level records.

Fig. 3 .

 3 Fig. 3. Temporal changes of M2 amplitude and phase lag for the 14 stations studied.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Standard deviations of amplitudes (on the left panel) and phases (right panel) over time for all studied stations. The blue rectangle is a zoom of the upper figure between 0 and 20 • . (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 6 .

 6 Fig.6. Skew surge differences from method 1 (using 1, 2, 5 and 19 years of data) and method 2 (circles: adjustment of MSL using a linear MSL trend on 60 years; stars: by adding annual MSL variations to the linear adjustment): a) RMSE, b) 95th percentile (P95) of differences, and c) correlation and d) RMSE between the differences and detrended MSL variation.

Table 1

 1 Site characteristics and sea level data used in this study.

	Station Name	Location [ • ]		Range of years	Number of cumulated years	Type	Tidal range [m]
		Longitude	Latitude
	Baltimore (USA)	39.27	-76.58	1902-2021	117.8	Mixed	0.90
	Boston (USA)	42.35	-71.05	1921-2018	99.1	Semi-diurnal	4.58
	Brest (France)	48.38	-4.50	1846-2020	158.8	Semi-diurnal	7.68
	Delfzijl (The Netherlands)	53.33	6.93	1879-2018	139.7	Semi-diurnal	1.18
	Gedser (Denmark)	54.57	11.93	1891-2020	122.5	Mixed	0.23
	Honolulu (USA)	21.31	-157.87	1905-2018	112.4	Mixed	0.98
	Key West (USA)	24.55	-81.81	1913-2021	106.1	Mixed	1.06
	Kungsholmsfort (Sweden)	56.11	15.59	1886-2021	133.7	Diurnal	0.15
	Marseille (France)	43.28	5.35	1849-2018	103.3	Mixed	0.44
	Newlyn (UK)	50.10	-5.54	1915-2021	103.5	Semi-diurnal	5.95
	Saint-Nazaire (France)	47.27	-2.20	1821-2020	114.2	Semi-diurnal	6.55
	San Francisco (USA)	37.81	-122.47	1898-2021	121.1	Mixed	2.87
	Seattle (USA)	47.60	-122.34	1899-2021	122.3	Mixed	5.37
	Sydney (Australia)	-33.85	151.23	1914-2019	103.6	Mixed	2.08

Table 3

 3 List of analysis and tests performed on the datasets in this study to estimate the uncertainties associated with the method of tidal prediction in the past.

	Investigation	Test	Results
	Influence of length of	HA on 1, 2, 3, 4, 5, 6, 7, 9 and	Fig. 2
	observations for HA on	19 years of observation
	high and low waters
	Evolution of 8 main TC	1-year rolling HA on 19 years	Figs. 3 and 4
	through time	of observation	Supplementary
			Material
	Evolution of annual MSL	Estimation of annual MSL for	Fig. 5
		each station
	Evolution of MSL trend	Linear and 2nd order	Table 4
	relative to the length of	polynomial trend estimation
	observation	on 30, 60 and 120 years
	Influence of the HA	Use of 1, 2, 5 and 19 years of	Table 5
	contemporary to the event	observation for HA,
	on 1% extreme skew	contemporary to an event
	surges
	Influence of TC current and	Use of 1, 2, 5 and 19 years of
	MSL trend adjustment on	current TC and MSL
	1% extreme skew surges	adjustment using:
		-30 years linear trend
		-60 years linear trend
		-120 years linear trend
		-120 years 2nd order
		polynomial trend

Table 6

 6

  present.

	Station	Linear (mm.yr -1 )			Acceleration (mm.yr -2 )
		1900-present	1960-present	1990-present	1900-present	1960-present	1990-present
	Baltimore	3.23 ± 0.05	3.56 ± 0.13	4.56 ± 0.38	0.0083 ± 0.0016	0.0703 ± 0.0088	0.2191 ± 0.0495
	Boston	2.84 ± 0.04	2.96 ± 0.09	4.98 ± 0.25	-0.0036 ± 0.0016	0.1148 ± 0.0056	-0.0689 ± 0.0331
	Brest	1.51 ± 0.04	2.04 ± 0.10	2.92 ± 0.28	0.0083 ± 0.0012	0.0862 ± 0.0066	-0.0122 ± 0.0369
	Delfzijl	1.97 ± 0.07	2.30 ± 0.22	2.98 ± 0.62	0.0095 ± 0.0024	0.0191 ± 0.0142	0.0810 ± 0.0841
	Gedser	1.21 ± 0.05	1.69 ± 0.15	2.25 ± 0.43	0.0053 ± 0.0017	0.0374 ± 0.0096	-0.0148 ± 0.0540
	Honolulu	1.50 ± 0.02	1.55 ± 0.05	2.48 ± 0.14	-0.0045 ± 0.0006	0.0336 ± 0.0033	0.1747 ± 0.0190
	Key West	2.53 ± 0.03	3.19 ± 0.07	4.36 ± 0.22	0.0137 ± 0.0011	0.0627 ± 0.0047	0.4787 ± 0.0280
	Kungsholmsfort	0.01 ± 0.04	0.65 ± 0.14	0.95 ± 0.37	0.0194 ± 0.0016	0.0045 ± 0.0090	0.2990 ± 0.0498
	Marseille	1.26 ± 0.03	1.31 ± 0.08	3.89 ± 0.43	0.0005 ± 0.0010	0.1441 ± 0.0060	-0.1601 ± 0.0861
	Newlyn	1.92 ± 0.04	2.17 ± 0.10	3.96 ± 0.28	0.0115 ± 0.0016	0.0841 ± 0.0065	-0.0673 ± 0.0370
	Saint-Nazaire	0.95 ± 0.05	1.48 ± 0.12	3.23 ± 0.34	0.0191 ± 0.0016	0.0808 ± 0.0084	-0.3359 ± 0.0445
	San Franciso	1.94 ± 0.02	1.92 ± 0.07	2.16 ± 0.20	0.0050 ± 0.0008	0.0012 ± 0.0046	0.4888 ± 0.0251
	Seattle	2.08 ± 0.03	2.07 ± 0.09	2.13 ± 0.28	0.0090 ± 0.0011	0.0055 ± 0.0063	0.3599 ± 0.0359
	Sydney	1.07 ± 0.03	1.26 ± 0.07	2.50 ± 0.20	0.0156 ± 0.0011	0.0487 ± 0.0044	0.2078 ± 0.0251

Table 6 99th

 6 Percentile (P99), date and maximum skew surge estimated using a MSL adjustment applying a 30, 60, 120 year linear trend or a polynomial trend on 120 years.

	N. Giloy et al.
	Station	Linear trend on 30 years		Linear trend on 60 years		Linear trend on 120 years		Polynomial trend on 120 years
		P99	Date	Surge	P99	Date	Surge	P99	Date	Surge	P99	Date	Surge
		[m]		max [m]	[m]		max [m]	[m]		max [m]	[m]		max [m]
	Baltimore	0.52	24/08/1933	2.07	0.45	24/08/	1.98	0.43	24/08/	1.95	0.48	24/08/	1.99
			04:00			1933 04:00			1933 04:00			1933 04:00
	Boston	0.49	07/02/1978	1.34	0.37	07/02/	1.25	0.37	07/02/	1.25	0.4	07/02/	1.28
			02:52			1978 02:52			1978 02:52			1978 02:52
	Brest	0.55	15/10/1987	1.44	0.44	15/10/	1.41	0.38	15/10/	1.39	0.4	15/10/	1.42
			09:29			1987 09:29			1987 09:29			1987 09:29
	Delfzijl	1.06	13/01/1916	3.37	0.99	30/10/	3.34	0.96	30/10/	3.33	0.98	30/10/	3.34
			05:29			2006 16:16			2006 16:16			2006 16:16
	Gedser	0.69	29/12/1913	1.96	0.65	29/12/	1.9	0.62	29/12/	1.84	0.64	29/12/	1.86
			10:00			1913 10:00			1913 10:00			1913 10:00
	Honolulu	0.21	22/05//	0.62	0.14	22/05//	0.57	0.14	22/05//	0.57	0.16	22/05//	0.59
			1960 00:32			1960 00:32			1960 00:32			1960 00:32
	Key West	0.25	23/10/2005	0.81	0.16	23/10/	0.79	0.13	23/10/	0.78	0.18	23/10/	0.84
			19:53			2005 19:53			2005 19:53			2005 19:53
	Kungsholmsfort	0.51	21/09/1912	1.34	0.49	11/09/	1.31	0.45	07/02/	1.24	0.49	11/09/	1.28
			12:23			1912 13:23			1913 11:30			1912 02:33
	Marseille	0.54	07/03/	1.13	0.3	07/03/	0.87	0.3	07/03/	0.87	0.33	07/03/	0.9
			1917 16:52			1917 16:52			1917 16:52			1917 16:52
	Newlyn	0.49	05/01/1936	1.03	0.37	16/12/	0.91	0.36	16/12/	0.9	0.38	16/12/	0.93
			13:54			1989 19:18			1989 19:18			1989 19:18
	Saint-Nazaire	0.66	11/11/1875	1.67	0.49	16/11/	1.43	0.45	16/11/	1.38	0.46	16/11/	1.42
			02:01			1940 16:27			1940 16:27			1940 16:27
	San Francisco	0.24	22/12/1982	0.7	0.23	22/12/	0.69	0.23	22/12/	0.69	0.25	22/12/	0.71
			23:42			1982 23:42			1982 23:42			1982 23:42
	Seattle	0.34	27/01/1983	0.82	0.34	27/01/	0.81	0.34	27/01/	0.82	0.36	27/01/	0.84
			12:36			12:36			1983 12:36			1983 12:36
	Sydney	0.32	06/07/1931	0.66	0.23	09/06/	0.56	0.22	09/06/	0.55	0.24	09/06/	0.58
			14:13			1974 12:57			1974 12:57			1974 12:57

Table A8

 A8 RMSE (in m) between tidal predictions made with MAS or UTide and tide gauge observations.

	M2 RMSE [• •] 0.08	0.03	0.03	0.05	0.31	0.06	0.07	1.1	0.48	0.04	0.07	0.07	0.03	0.04
	RMSE [%] 99.98	99.88	99.95 99.97	98.98	99.99	99.99	99.89	99.83	99.97	99.9	99.96	99.96	99.86
	S2 RMSE [• •] 0.22	0.02	0.01	0.09	1.03	0.04	0.03	3.24	0.17	0.02	0.07	0.03	0.02	0.02
	RMSE [%] 99.2	99.95	99.97 99.8	98.6	100	100	98.77	99.99	99.96	99.94	99.97	99.93	99.99
	N2 RMSE [• •] 0.24	0.08	0.06	0.12	0.85	0.12	0.11	110.92	0.39	0.08	0.08	0.07	0.06	0.08
	RMSE [%] 99.41	99.83	99.88 99.78	93.8	99.99	99.98	24.1	99.93	99.87	99.92	99.96	99.62	99.7
	K2 RMSE [• •] 0.88	0.2	0.2	0.4	4.32	0.39	0.44	45.33	0.89	0.25	0.31	0.22	0.15	0.13
	RMSE [%] 97.59	95.66	99.18 98.55	90.63	99.84	99.75	27.87	99.37	98.65	98.53	99.54	98.75	99.6
	Q1 RMSE [• •] 1.26	0.53	0.69	1.57	108.82 0.33	0.23	8.55	1.95	1.08	0.71	0.14	0.12	0.31
	RMSE [%] 96.14	77.93	80.5	73.43	54.89	99.41	98.79	82.8	94.89	75.49	85.54	96.99	94.17	81.3
	O1 RMSE [• •] 0.31	0.07	0.15	0.49	1.39	0.11	0.06	2.54	0.31	0.18	0.4	0.06	0.05	0.07
	RMSE [%] 99.14	98.66	95.91 72.65	95.88	99.93	99.95	99.02	99.6	98.38	96.58	99.71	99.39	99.11
	P1 RMSE [• •] 0.52	0.06	0.11	2.08	2.67	0.03	0.06	4.9	0.25	0.55	0.6	0.04	0.03	0.06
	RMSE [%] 97.97	99.67	99.82 63.29	94.96	100	99.92	96.59	99.81	83.9	93.75	99.8	99.26	99.61
	K1 RMSE [• •] 0.29	0.08	0.06	0.88	0.78	0.1	0.04	1.65	0.42	0.16	0.24	0.04	0.02	0.04
	RMSE [%] 99.11	97.3	99.65 72.73	94.76	99.96	99.95	99.16	99.46	91.06	94.39	99.83	99.41	99.22
		Station			RMSE MAS -Observations		RMSE UTide -Observations
		Baltimore			0.20				0.19
		Boston			0.15				0.14
		Brest				0.14				0.13
		Delfzijl			0.33				0.32
		Gedser			0.23				0.23
		Honolulu			0.07				0.06
		Key West			0.11				0.07
		Kungsholmsfort		0.18				0.18
		Marseille			0.11				0.11
		Newlyn			0.13				0.13
		Saint-Nazaire			0.16				0.16
		San Francisco			0.10				0.08
		Seattle			0.13				0.13
		Sydney			0.09				0.09

Acknowledgements

This work was funded by the "Direction générale de la Prévention des risques" (DGPR) of the French Ministry of Ecology and by French Institute for Radiological Protection and Nuclear Safety (IRSN).

The authors would like to acknowledge Haigh et al. (2022), Woodworth et al. (2016) and Caldwell et al. (2015) who compiled all the observations available from more than 5100 stations records and make it available online with the GESLA data base (https://gesla787883612. wordpress.com/).

Data availability

We used the GESLAv3 dataset available online https:// gesla787883612.wordpress.com/

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.wace.2023.100597.

Table A 7

Statistical values of the comparison between the tidal constituents' (TC) phase results provided by MAS and UTide softwares. Respectively, root-mean-square error (RMSE) and the correlation coefficient (R).