Paul Platzer
email: paul.platzer@ifremer.fr

Pierre Tandeo

Pierre Ailliot

Bertrand Chapron

Could old tide gauges help estimate past atmospheric variability ?

The storm surge is the non-tidal component of coastal sea-level. It responds to the atmosphere both through the direct effect of atmospheric pressure on the sea-surface, and through Ekman transport induced by wind-stress. Tide gauges have been used to measure the sea-level in coastal cities for centuries, with many records dating back to the 19th-century or even further, at times when direct pressure observations were scarce. Therefore, these old tide gauge records may be used as indirect observations of sub-seasonal atmospheric variability, complementary to other sensors such as barometers. To investigate this claim, the present work relies on tide gauge records of Brest and Saint-Nazaire, two portal cities in western France, and on the members of NOAA's 20th-century reanalysis (20CRv3) which only assimilates surface pressure observations and uses numerical weather prediction model. Using simple statistical relationships between storm surges and pressure maps, we show that the tide gauge records reveal part of the 19th-century atmospheric variability that was uncaught by the pressure-observations-based reanalysis. In particular, weighing the 80 reanalysis members based on tide gauge observations indicates that a large number of members are very unlikely, which induces corrections of several tens of Hectopascals in the Bay of Biscay. These findings support the use of early tide gauge records in sensor-scarce areas, both to validate old atmospheric reanalyses and to better probe old atmospheric sub-seasonal variability.

Introduction

Understanding the atmospheric system requires to understand all scales of variation, from daily to centennial. This cannot be done unless long observation records allow to disentangle these scales. The Twentieth Century Reanalysis Project, hereafter "20CR" [START_REF] Compo | The twentieth century reanalysis project[END_REF], which is now in its third version, hereafter "20CRv3" [START_REF] Slivinski | Towards a more reliable historical reanalysis: Improvements for version 3 of the Twentieth Century Reanalysis system[END_REF], is the only atmospheric reanalysis that runs through the 19th century. It relies on the International Surface Pressure Databank [START_REF] Compo | The international surface pressure databank version[END_REF], the largest historical global collection of surface pressure observations, and the NCEP Global Forecast System (GFS) coupled atmosphere-land model.

Because it is the longest atmospheric reanalysis available, the 20CR reanalysis is used to study possible long-term trends in atmospheric dynamics [START_REF] Rodrigues | Dynamical properties of the North Atlantic atmospheric circulation in the past 150 years in CMIP5 models and the 20CRv2c reanalysis[END_REF] or for extreme events [START_REF] Alvarez-Castro | Atmospheric dynamics leading to West European summer hot temperatures since 1851[END_REF]. However, although the 20th-century part of 20CR has been compared with other reanalysis [START_REF] Wohland | Inconsistent wind speed trends in current twentieth century reanalyses[END_REF] and observations (Krueger et al., 1 https://doi.org/10.5194/egusphere-2023-2997 Preprint. Discussion started: 2 January 2024 c Author(s) 2024. CC BY 4.0 License. 2013), comparisons with independent observations in the 19th century [START_REF] Brönnimann | Early ship-based upper-air data and comparison with the Twentieth Century Reanalysis[END_REF] are scarce. The present work is an effort to compare this reanalysis with tide gauge observations. More generally and to our best knowledge, this paper is the first attempt to use old tide gauges as indirect observations of the atmosphere. However, the opposite direction has been taken by [START_REF] Tadesse | A database of global storm surge reconstructions[END_REF], who extended storm surge reconstructions in the past using different atmospheric reanalysis products, in order to estimate past unobserved extreme storm surges.

Tide gauges are used primarily to measure the tide, which is the largest contributor to sea-level variations in many coastal cities. The astronomical tide is the result of gravitational attraction of the Sun and Moon on the ocean, combined with Earth's rotation. It results in periodic rise and fall of the water level [START_REF] Melchior | The tides of the planet Earth[END_REF], which have been predicted through harmonic decomposition for centuries. Other physical phenomena impact the water level: a low atmospheric pressure results in a high sea-level, a phenomenon called the "inverse barometer effect" [START_REF] Roden | Early Swedish contribution to oceanography: Nils Gissler (1715-71) and the inverted barometer effect[END_REF][START_REF] Woodworth | Forcing factors affecting sea level changes at the coast[END_REF], and wind blowing parallel to the coast results in Ekman transport either towards or away from the coast, increasing or decreasing the sealevel. These conditions are usually associated with storms, which is why the associated sea-level variations are called "storm surges". For instance, in Brest (France), the amplitude of tidal variations is close to 4m, and storm surges can amount to as much as 1.5m. Tide gauges are numerous, forming a dense global network in recent years, and a sparser one in the last centuries. As an example and from the GESLA-3 sea-level database [START_REF] Haigh | GESLA version 3: A major update to the global higher-frequency sea-level dataset[END_REF], 10 coastal tide gauge records start before 1907 in the North-American east coast, and 20 start before 1900 in Europe. Old tide gauges have varying observation frequencies, from hourly [START_REF] Wöppelmann | Brest sea level record: a time series construction back to the early eighteenth century[END_REF] to daily averages [START_REF] Marcos | Historical tide gauge sea-level observations in Alicante and Santander (Spain) since the 19th century[END_REF]. Although the sea-level measured by tide gauges is only an indirect tracer of atmospheric pressure variability, the scarcity of direct sea-level pressure measurements motivates the use of tide gauges to study past atmospheric fluctuations. Indeed, even when pressure measurements exist, they are often not yet digitized and even less available in global repositories [START_REF] Brönnimann | Unlocking pre-1850 instrumental meteorological records: A global inventory[END_REF].

It is possible to link sea-level variations with atmospheric phenomena using physical laws and models [START_REF] Lazure | An external-internal mode coupling for a 3D hydrodynamical model for applications at regional scale (MARS)[END_REF], or using statistical tools [START_REF] Quintana | Statistical forecast of the marine surge[END_REF]. This work adopts the second approach, but the underlying physical phenomena will often be used to motivate and interpret the statistical models. Multivariate linear regression (LR) will be used to relate the storm surge to local mean-sea-level pressure anomaly and pressure gradients. Hidden Markov Models (HMM) will allow to perform time-smoothing of probabilities given to members of 20CRv3, taking advantage of the time-continuity of each member.

We stress that the reader should be more concerned with the underlying ideas of this work than with the particular methodology and datasets used. Indeed, most of what is exposed here can be reiterated with other old tide gauge records, in particular in Europe and North-America, and using other statistical tools or physical models. The Brest and Saint-Nazaire tide gauges have the advantage of both being hourly sampled even in the 19th century, and of being close to each other. The statistical methods were chosen for their simplicity of use and interpretation.

The data and preprocessing are detailed in section 2. Section 3 studies changes in linear relationship between storm surges and 20CRv3 members SLP in the 19th century. Section 4 explores corrections in sea-level pressure induced by the conditional 2 https://doi.org/10.5194/egusphere-2023-2997 Preprint. Discussion started: 2 January 2024 c Author(s) 2024. CC BY 4.0 License. weighing of 20CRv3 members based on surge observations. Conclusions on the proposed methodology and experiments are drawn in section 5, along with potential applications of this work.

Data

The Twentieth Century Reanalysis version 3 (20CRv3)

The Twentieth Century Reanalysis Project [START_REF] Compo | The twentieth century reanalysis project[END_REF] aims at producing a global atmospheric reanalysis ending in 2015 and extending back to the 19th century. The present paper uses the latest version, 20CRv3 [START_REF] Slivinski | Towards a more reliable historical reanalysis: Improvements for version 3 of the Twentieth Century Reanalysis system[END_REF], which extends up to 1806. It is an atmospheric reanalysis with 80 members, using an Ensemble Kalman Filter data assimilation scheme [START_REF] Evensen | The ensemble Kalman filter: Theoretical formulation and practical implementation[END_REF]. It assimilates only surface pressure observation, from ships and fixed stations, as well as analysed cyclone-related IBTrACs data. These surface pressure observations are taken from the International Surface Pressure Databank (ISPD) which was created for the 20CR project but also exists as an independent product [START_REF] Compo | The international surface pressure databank version[END_REF]. In 20CR, the seasurface temperature and sea-ice cover are prescribed as boundary conditions. Sea-surface temperature and sea-ice cover both benefit from satellite observations from 1981 to 2015 (the end of the reanalysis), allowing more precise boundary conditions.

The surface pressure observation density is considerably lower in the 19th century then in the late 20th century. An online platform (https://psl.noaa.gov/data/20CRv3_ISPD_obscounts_bymonth) allows to consult the monthly observation count per 2 • ×2 • box. Fig. 1 shows yearly averages of the number of surface pressure observations per day, comparing years 1870 and 2000. The maximum value was set to 24 observations per day although in 2000 this value is mostly exceeded. In year 1870, approximately half of Europe's land surface has no observation at all, and only less than 10 points have more than 10 observations per day. Observations coming from ships allow to raise the number of observations to approximately one per day on dense traffic areas. Conversely, in year 2000, virtually all of western Europe's land has more than 24 observations per day.

Taking a spatial average over the whole map from Fig. 1 gives approximately one observation every three days in 1870, versus 44 observations per day in 2000. The number of available observations is also highly variable through time, especially in the 19th century. For instance, in the 2 • ×2 • box centered on 49 • -latitude, -5 • -longitude, the number of monthly observations in 1870 ranges from 2 (January, 1870) to 85 (May, 1870), while in 2000 it ranges from 2152 (June, 2000) to 3242 (May, 2000). This double constraint of data scarcity and high variations in sampling frequency in the 19th century legitimizes the search for other sources of observations to study centuries-old atmospheric variability.

Anomalies are considered with a reference climatology computed as an average over all members and over ±30 calendar days, ± 3 day hours (to include diurnal effects). In the following, we use two different periods, one is called hereafter "the satellite-era" because in this period the reanalysis is forced by satellite-derived sea-surface temperature and sea-ice cover, the second is 1870-1896. A different climatology is computed for each of these two periods to retrieve anomalies.

Tide gauges of Brest and Saint-Nazaire (France)

In this study, the tide gauges of Brest and Saint-Nazaire are used as indirect tracers of atmospheric circulation through storm surges. The storm surges of Brest and Saint-Nazaire are highly correlated (correlation coefficient of 0.92 using the whole overlapping period). The two stations are 206km apart, which is small compared to the typical size of storms and anticyclones.

The sea-levels in the two cities usually respond similarly to passing atmospheric systems. Therefore, one station can fill observational gaps of the other. However, when data from both stations is available, the distance between them is large enough to potentially get information on the direction, velocity, and size of passing systems, as the stations could provide a kind of gradient along the line between the two.

The availability of sea-level records in Brest and Saint-Nazaire in the GESLA-3 database is shown in Fig. 2. Brest's tide gauge with hourly sampling starts in 1846, while Saint-Nazaire starts in 1863. Apart from a few large gaps, both records are mostly continuous during periods 1863-1920 and 1953-present. This combination of historical and actual records is at the foundation of the methodology exposed in the next section.

Preprocessing of sea-level data

As mentioned earlier, the part of the sea-level which responds to atmospheric processes is the surge (also called "storm surge" or "skew surge"). To access the surge, one has to remove the tidal part of the signal, along with a possible time-linear trend in the 20th century due to sea-level rise [START_REF] Cazenave | Contemporary sea level rise[END_REF]. To remove the time-linear trend in Brest's and Saint-Nazaire's signals, we first compute the yearly average of the sea-level, excluding year 1952 in Brest because it contains less than 3 month of observations which leads to an abnormally high value.

Then, we compute the average of these levels before 1900, and perform a one-parameter, time-linear, least-squares fit starting in 1900 with the value in January 1900 fixed to the 19th-century average. This gives a two-part function for each tide gauge, constant before 1900 and time-linear after 1900. Removing these functions to the signal allows to access detrended sea-levels.

Then, the tidal constituants of the detrended sea-levels are computed using U-Tide [START_REF] Codiga | Unified tidal analysis and prediction using the UTide Matlab functions[END_REF], which performs harmonic (Fourier) decomposition with prescribed frequencies corresponding to planetary movements. This allows to compute and remove the tidal part of the sea-level, and access the surge. Then, we perform a 12h-moving average of the surge to filter out oscillations occurring at a frequency close to 12h. These oscillations are either due to tide-surge interactions [START_REF] Horsburgh | Tide-surge interaction and its role in the distribution of surge residuals in the North Sea[END_REF] or to measurement errors in the 19th century leading to phase shifts. Furthermore, tide-surge interactions lead to stronger surges in low-tide and weaker surges in high-tide [START_REF] Horsburgh | Tide-surge interaction and its role in the distribution of surge residuals in the North Sea[END_REF]. As these phenomenons are not linked to atmospheric processes, we chose to filter them out with a simple 12h-average. This also implies that these 12-hours-averaged surges will only respond to atmospheric events persisting for more than 12 hours. 3 Linear Regression (LR) between surges and sea-level pressure 3.1 LR in the satellite-era To estimate the statistical relationship between surges in Brest and Saint-Nazaire, and neighbouring fields of sea-level pressure, we first focus on the period 1981-2015 during which satellite data is used in 20CRv3 to constrain sea-surface temperature and sea-ice cover, and a large number of pressure observations gives high confidence in 20CRv3 fields of pressure reduced at mean sea-level (SLP in the following).

The filtered surges described in the previous section respond to sub-seasonal variations in atmospheric pressure and winds at time scales greater than 12 hours. First, the surges in Brest and Saint-Nazaire are highly sensitive to pressure anomalies, a phenomenon called the "inverse barometer effect" [START_REF] Roden | Early Swedish contribution to oceanography: Nils Gissler (1715-71) and the inverted barometer effect[END_REF]). An increase (respectively, decrease) of 1hPa in pressure at the mean sea-level leads to a decrease (respectively, increase) in sea-level of approximately 1cm, which can be expressed as a hydrostatic equilibrium between the atmospheric and water columns.

Second, wind imposes a stress on the sea-surface, creating fluid motion through the water column in a shape known as the Ekman spiral, governed by the balance between Coriolis and turbulent drag forces. At the top of this spiral, the ocean surface current deviates at 45 • to the right of the wind stress in the northern hemisphere. The resulting net fluid motion is called the "Ekman transport" [START_REF] Price | Wind-driven ocean currents and Ekman transport[END_REF], and is perpendicular to the wind stress direction, to the right in the northern hemisphere.

This transport has a non-linear influences on storm surges [START_REF] Bryant | An exploration of wind stress calculation techniques in hurricane storm surge modeling[END_REF][START_REF] Pineau-Guillou | Strong winds in a coupled wave-atmosphere model during a North Atlantic storm event: Evaluation against observations[END_REF].

To take advantage of these well-known physical processes, we use a simple linear regression as a statistical learning strategy, based on the 12h-averaged surges and data from 20CRv3. We take as regressors or "explanatory variables" the anomaly of pressure reduced at mean sea-level (SLP B and SLP SN in the following) interpolated at the coastal city of interest (Brest or Saint-Nazaire), and the two components of a local, horizontal gradient of SLP (∆ lon SLP B , ∆ lon SLP SN and ∆ lat SLP B , ∆ lat SLP SN ) expressed as longitudinal and latitudinal differences at ±2 • around the city of interest. In the satellite-era, we use the average SLP anomaly over all 80 members of the reanalysis, as in this period the ensemble spread in SLP between members is considerably low (less than 100Pa) due to the high observation density.

To summarize, the regressions in the satellite-era have the following form:

  s B s SN   ∼ N     α B SLP B + β B ∆ lat SLP B + γ B ∆ lon SLP B α SN SLP SN + β SN ∆ lat SLP SN + γ SN ∆ lon SLP SN   ,   σ 2 LR, B Cov B, SN Cov B, SN σ 2 LR, SN     (1)
where s B and s SN are the surges in Brest and Saint-Nazaire, considered as normal (i.e. Gaussian) random variables with aver- Also, we assume homoscedasticity: the variance of the surges are independent of the values of SLP, only the average depends on SLP values. These variances σ 2 sat,B and σ 2 sat,SN are estimated as the average of squared differences between the actual value of the surges and the SLP-based average with the LR coefficients:

σ 2 LR = {s(t) -(αSLP(t) + β∆ lat SLP(t) + γ∆ lon SLP(t))} 2 t∈[1981-2015] , (2)
where angle brackets mean average, and we have removed the indices B or SN to improve readability of the formula. This allows to define the coefficient of determination R 2 as follows:

R 2 = 1 - σ 2 LR (s -⟨s⟩) 2 , (3)
which can be defined for any city (Brest or Saint-Nazaire) and any period (here the averages are implicitly taken in the satellite-era). A value of 1 indicates perfect prediction, while a value of zero is obtained for a constant, average prediction.

The coefficient of determination R 2 can be understood as a percentage of explained variance.

The values of the LR coefficients and coefficients of determination are shown in Table 1. For comparison, a one-dimensional LR is performed with only α and imposing β = γ = 0. All three coefficients of these linear regressions can be understood from a physical point of view. The first coefficient, in front of SLP, is the "inverse barometer". A theoretical hydrostatic equilibrium between the atmosphere and ocean would give a coefficient of -1.02, which is very close to the two coefficients found in these regressions (see Table 1). The coefficients β and γ of the linear regression indicate the most favorable geostrophic wind directions [START_REF] Obukhov | On the question of geostrophic wind[END_REF], as well as the intensity of the interaction between geostrophic wind stress and surge. This allows to estimate the most surge-favorable wind directions according to the linear regressions, as shown in Fig. 4(c). Larger coefficients in Saint-Nazaire can be interpreted as a consequence of a shallower sea near Saint-Nazaire, leading to a stronger influence of winds, as can also be witnessed in other regions such as the North-Sea [START_REF] Pineau-Guillou | Impact of wave-dependent stress on storm surge simulations in the North Sea: Ocean model evaluation against in situ and satellite observations[END_REF]. For Brest, accounting for geostrophic wind only allows to explain an additional 6% of the surge variance (comparison of R 2 with one-dimensional linear regression, Table 1). For Saint-Nazaire, the geostrophic winds allow to explain an additional +12% of surge variance, which confirms the stronger influence of winds in Saint-Nazaire. Furthermore, although the linear coefficient related to the inverse barometer should approach -1.02 on physical grounds, it departs from this value when using the one-dimensional LR and neglecting the influence of wind for both portal cities (Table 1). This justifies the use of a threedimensional LR although most of the variability is explained by the first coefficient α.

A larger spread around the linear fit is observed in Saint-Nazaire (σ SN = 6.8cm) than in Brest (σ B = 5.4cm). This difference could reflect the difference in contributions of physical phenomenon such as surface currents and wave setup, which are not captured by the LR (by construction).

LR in the 19th century (1870-1896)

Is the statistical relationship between the surges in Brest and Saint-Nazaire and the SLP from 20CRv3 consistent through time ? The answer to this question is delicate, but could provide both an independent means of validation for 20CRv3 and a justification for the use of old tide gauges to probe 19th-century atmospheric variability. One way to try and answer this question is to compare the linear relationships established in section 3.1 using satellite-era 20CRv3 member's average with the linear relationships using 19th-century 20CRv3 members.

First, one must note that the spread in SLP between 20CRv3 members is a decreasing function of time, as it is anti-correlated with the observation density from the ISPD. The statistical relationship between surges and members-averaged SLP is therefore bound to evolve through time. More precisely, the 19th-century members-average is expected to be smoother, with weaker extreme values.

Fig. 5 provides a visual inspection as to whether the LR established in Fig. 4(a) between Brest surge and SLP fields holds in the 19th century, using either the members average (Fig. 5b) or individual members (Fig. 5c,d). In Fig. 5b,a bias is visible, the satellite-era LR consistently underestimating the amplitude moderately high and low values of surge. Indeed, for predictions around +20cm (respectively -20cm) there are more values above than below (respectively below than above) the x = y diagonal line. This is consistent with the expectation that the ensemble average underestimates SLP variability due to smoothing in the 19th century. However, the LR seems to hold for low values. A negative bias is observed towards positive extremes (i.e., the highest surges are underestimated) but this phenomenon also holds with the satellite-era data. This interpretation is further confirmed by comparing standard deviations in SLP, ∆ lon SLP and ∆ lat SLP, in the period 1864-1896 and in the period 1981-2015. Table 2 shows that, when taking the 20CRv3 ensemble average in Brest, the SLP and its gradients have a lower standard-deviation in 1864-1896 compared to 1981-2015, by more than 4%. and since this first term in the LR is the main driver of the response in surge to the atmosphere one finds a similar amplitude of linear surge-prediction based on individual members' SLP. Table 2 also indicates larger standard deviations in horizontal gradients of SLP in the 19th century, which could be interpreted as a tendency of the model used in 20CRv3 to create stronger gradients of SLP when not constrained by observations. In Brest, this has a minor impact on the LR between surges and SLP, as winds only mildly contribute to the surges. However, the situation might be different when using tide gauges of shallower seas which are more sensitive to wind conditions. These visual inspections can be refined by computing the LR as in section 3.1 but using the data from 1864-1896 (either ensemble average or individual members) rather than 1981-2015. The coefficients α, β, γ are reported in Table 3 for various choices of SLP in years 1864-1896, compared to the coefficients from the period 1981-2015. The first coefficient, α, is slightly larger with the 1864-1896 ensemble-average, and slightly smaller with the 1864-1896 ensemble members, which can be interpreted as a consequence of a smoothing in the average and an added variability in the individual members. The other coefficients β and γ are harder to interpret, because they contribute less to the explained variability of the surge (α sat alone explains 79% of the variability in 1981-2015, thus β sat and γ sat only explain the remaining 4%), and because of the differences in behaviour of the SLP gradients in 1864-1896 when considering ensemble averages or individual members as depicted in table 2. However, the angles of the corresponding geostrophic wind directions with respect to a purely westerly wind are very similar, respectively 80.3 • (1981-2015), 77.5 • (1864-1896, ensemble average), and 69.2 • (1864-1896, individual members).

Therefore, in any situation winds blowing from the south are associated with larger surges. The amplitude of the coefficient is similar between 1981-2015 and 1864-1896 ensemble average, while it is smaller for 1864-1896 individual members, which can be interpreted again as a consequence of stronger longitudinal gradients of SLP for less-constrained individual members of 20CRv3 (table 2).

Finally, coefficients of determination R 2 are given in table 3, using the same formula as in Eq. ( 3), but computin σ 2 LR either with "1864-1896" coefficients (i.e., the coefficients indicated in the second and third lines of the table) or with the satellite-era coefficients (i.e., coefficients in the first line of the table ). One can see that the LR explains 83% of the surge variability in 1981-2015, which is a quantitative view of the high correlation observed in Fig. 4(a). For the other time period 1864-1896 and taking the ensemble average, the percentage of explained variability with the optimal choice of coefficients is 79%, which is very close to the explained variability when using the coefficients computed in the satellite-era, indicating a good persistence in time of the statistical LR between Brest's surge and 20CRv3's ensemble average SLP. The slight decrease in explained variance with respect to 1981-2015 can be explained by the fact that the ensemble average SLP in 1870-1896 is less constrained by observatons, and therefore more likely to deviate from the truth. This added variability in SLP blurs the statistical relationship between SLP and surges. Finally, for individual members, another layer of variability is added, corresponding to the ensemble spread, so that individual members are even further from the true 19th-century atmospheric state, resulting in an explained variance of 69.9%. Again, however, a very similar score is explained when using the coefficients derived in the satellite-era LR, confirming the persistence through time of the statistical LR between the Brest surge and the local SLP and its gradients. 3. Coefficients of the linear regression between, on the one hand, the surge in Brest, and on the other hand, the SLP (coefficient α) and its horizontal gradients (coefficients β and γ). The satellite-era coefficient are compared with coefficients in 1870-1896, with values for the SLP taken either as the ensemble average of 20CRv3 members, or individual members. For the individual members, values given are in the format : average ± standard-deviation over the 80 members. Coefficients of determination R 2 are given in two case scenarios: the true surge value is compared to its predicted value either according to the LR with actual coefficients α, β, γ on the corresponding line, or according to the LR with satellite-era coefficients αsat, βsat, γsat (from the first line 1981-2015).

The same behaviour can be witnessed for Saint-Nazaire (not shown for conciseness). This investigation of the LR between observed surges and 20CRv3-SLP first shows that the statistical relationship between these two quantities appears to remain valid in the 19th century, which is by itself a validation result for 20CRv3 using an independent, indirect source of observation of atmospheric variability. Second, it appears that the added variability observed in the end of the 19th century is not due to a modified linear relationship, but rather to an increase of variability in 20CRv3. The fact that this increase is visible in the linear projection on the surges of Brest and Saint-Nazaire indicates that the surges of Brest and Saint-Nazaire contain valuable information on the atmospheric variability that could be used to reduce the uncertainties in 20CRv3. To reformulate, this analysis suggests that the surges contain information on the 19th-century atmospheric variability that is complementary with the observational information gather in 20CRv3 through surface pressure observations from the ISPD.

The next section explores the potential to use Brest and Saint-Nazaire surge observations to constrain the 20CRv3 through weighing of its members.

Probability of 20CR members based on surge observations

Hidden Markov Model (HMM)

In the 19th century, the spread between 20CRv3 members is much larger than in the period 1981-2015, and one of the aims of this work is to estimate conditional probabilities of each member of the reanalysis, based on surge values in Brest and Saint-Nazaire. Note that in the reanalysis, the members are assumed to have uniform probabilities, that is a probability of 1/80, since we have 80 members. Therefore, one aim of this work is to evaluate the change in probability distribution of the members when adding the information of the surges in Brest and/or Saint-Nazaire. One can estimate conditional probabilities of each member at time t based on the values of the surges at time t. To do that, we use the satellite-era-derived linear regression expressed in section 3.1. Using Baye's theorem, we know that the probability to observe a SLP map conditionally on observing a surge value is proportional to the probability of observing a surge given the SLP map, which we assume can be expressed simply by the satellite-era LR.

To differentiate these member probabilities from the ones we will derive later on using a hidden Markov model, we use the notation p HMM (i, t) for the probability of member i at time t.

     p HMM (i, t) ∝ p sat (s -[α sat SLP i (t) + β sat ∆ lon SLP i (t) + γ sat ∆ lat SLP i (t)]) , 80 i=1 p HMM (i, t) = 1 . ( 4
)
Although these probabilities already bear significant information, they have the undesirable property to be time-discontinuous. This is not coherent with the fact that the members of 20CRv3 are time-continuous: they are propagated in time using a NWP model. To remedy this issue, we compute smoothed (or reanalyzed) probabilities using a hidden Markov model (HMM) detailed below, which we write p HMM (i, t):

p HMM (i, t) := P      member(t) = i      surge(t = 1) . . . surge(t = T )      =      s 1 . . . s T           , (5)
where one uses an observational record of surges from time-index 1 to T . This quantity is a time-smoothed version of p HMM (i, t)

which takes into account past and future values of the surge. For this purpose, a simple Hidden Markov Model (HMM) is used.

The first ingredient of the HMM is the transition matrix T ij (t) from member i at time t -1 to member j at time t.

T ij (t) := P (member(t) = j | member(t -1) = i) . ( 6
)
To estimate the transition matrix, a strong hypothesis is made:

T ij (t) ∝ K θ (SLP map,j (t), SLP map,i (t)) , (7)
where SLP map,i (t) is the i-th member's map of mean-sea-level pressure in a squared box of 18

• W ≤ lon ≤ 18 • E, 28 • N ≤ lat ≤ 64 • N at time t, and K θ (•,
•) is a positive real-valued function that measures the similarity between SLP map,i (t) and SLP map,j (t) and depends on parameters θ. Note that in case of missing observations between t and t + q -1, where q ≥ 1, we assume the following form for the more general transition matrix T ′ :

T ′ ij (t, q) := P (member(t + q) = j | member(t -1) = i) , (8)
∝ {K θ (SLP map,j (t + q), SLP map,i (t + q))}

τ τ +q , (9)
where τ > 0 is a hyperparameter. We set τ to 16 in time-step unit, and since we use time-steps of 3 hours, our τ corresponds to a time of 2 days. When q is large compared to τ , the coefficients of the transition matrix become independent of i, j and transitions from any state to another are equally allowed (i.e., memory is lost).

Eq. ( 7) states that transitions from one member to another are more likely if the associated SLP map at time t are similar.

This prevents abrupt transitions to dissimilar atmospheric states. The size and location of the map was chosen to cover an area inside which storms and anticyclones which affect the surges in Brest and Saint-Nazaire would lie. Ideally, K θ (•, •) should be symmetric, semi-definite. Here, a simple Gaussian kernel of Euclidean distances is used, with normalization factor θ > 0, so that for two fields X and Y :

K θ (X, Y ) = exp - n ∈ lons l ∈ lats (X nl -Y nl ) 2 θ 2 , (10)
where the sum over n and l represents a sum over longitudes and latitudes. We then define Θ through the following equation:

θ Θ = 1 80 80 i=1 1 T T t=0 SLP map,i (t) 2 1/2 , ( 11
)
where the overbar denotes spatial average. This normalization will allow to optimize θ through grid search of Θ for a maximum of likelihood of the surge observations.

One can compute T ij (t) by setting a value of θ and using the hypothesis of Eq. ( 7) along with the fact that for all i, t, we have j T ij (t) = 1. This then allows to estimate p HMM (i, t) with the forward-backward algorithm [START_REF] Rabiner | A tutorial on hidden Markov models and selected applications in speech recognition[END_REF]. Let:

a i (t) :=P           surge(1) . . . surge(t)      =      s 1 . . . s t      , m(t) = i      , ( 12
) b i (t) :=P           surge(t + 1) . . . surge(T )      =      s t+1 . . . s T      m(t) = i      . (13)
These two quantities can be computed recursively, following the forward procedure:

a i (t + 1) = p HMM (i, t + 1) 80 j=1 a j (t)T ji (t) , (14)
and the backward procedure:

b i (T ) = 1 , ( 16
) b i (t) = 80 j=1 b j (t + 1)T ij (t)p HMM (j, t + 1) . (17)
Finally, this allows to estimate p HMM (i, t) by noting that:

p HMM (i, t) = P      member(t) = i ,      s(1) . . . s(T )      =      s 1 . . . s T           P           s(1) . . . s(T )      =      s 1 . . . s T           , (18)
which gives, in terms of a i (t) and b i (t):

p HMM (i, t) = a i (t)b i (t) 80 j=1 a j (t)b j (t) (19)
while keeping in mind that Eq (19) implicitly relies on hypothesis (7) and a fixed form of K θ .

Comparing p HMM (i, t) with the uniform distribution p(i, t) = 1 80 allows to see if the surge observations are coherent with the SLP fields of 20CRv3 (section 4.2) and to select the most relevant members given surge data (section 4.3).

To choose the parameter θ, we performed a grid-search of its normalized form, Θ, computed the log-likelihood of the surge observations as an output of the algorithm. Indeed, the log-likelihood l θ (0 . . . T ) is expressed as follows:

l θ (0 . . . T ) = log( 80 i=1 a i (T )) . (20)
Figure 6 shows variations of this quantity with Θ, for one year (1885) of surge observations in Brest (i.e. t = 0 is 01 January 1885 and T is 01 January 1886). The curve shows a distinct maximum around Θ ≈ 0.09, and plateaus for higher values. 11). For comparison, the loglikelihoood of the simple model without Hidden Markov Model (Θ → +∞) is also shown (dashed grey line). The log-likelihood was estimated using data from year 1885 for a first estimation of optimal parameter Θ. Values of Θ used for estimation are 0.02, 0.05, 0.07, 0.08, 0.09, 0.1, 0.11, 0.13, 0.15, 0.2, 1, and 10, with 0.09 giving the largest log-likelihood.

According to Figure 6, the difference of log-likelihood between the model without HMM (θ = +∞) and with HMM is close to 1000. The introduction of one extra parameter in the filtering model compared to the static one is thus clearly justified if the two models are compared using standard criteria usch as AIC, BIC or likelihood ratio tests [START_REF] Zucchini | Hidden Markov models for time series: an introduction using R[END_REF].

Note that in the limit Θ = +∞, we have a constant transition probability T ij (t) = 1 80 and p HMM (i, t) reduces to p HMM (i, t). Figure 6 thus supports the use of the HMM to estimate probabilities of SLP map conditioned by surge observations. The choice of restricting the estimation of log-likelihood to one arbitrary year (1885) is supported by the fact that estimation of T ij (t) is computationally expansive. We assume that the optimal value of θ generalizes well to other years. A better optimization of θ would necessitate further work that is out of the scope of this study. Setting Θ = 0.09 will already enable us to find interesting features of p HMM (i, t).

Modification of 20CRv3 ensemble when accounting for surges

This section is devoted to the study of δµ HMM (t), the difference between weighted and unweighted ensemble average, defined by:

δµ HMM (t) := 80 i=1 p HMM (i, t) - 1 80 SLP map,i (t) , (21)
where SLP map,i (t) is a short notation for the sea-level pressure field of 20CRv3's i-th member. δµ HMM (t) is defined equivalently using p HMM (i, t). This quantity shows how strong is the average deviation when taking into account surge observations. It will also be sometimes normalized by σ 20CR (t), the estimated standard deviation of the unweighted ensemble: SLP map,i (t) -1 80

80 i=1 SLP map,i (t) 2   1/2 . ( 22
)
Note that in this definition, σ 20CR (t) depends on time, latitude and longitude. Therefore at each grid point and for each time step the quantity δµ HMM (t) will be normalized by a different value, indicating the strength of the reanalysis ensemble spread at this location in time and space.

To further interpret the result of our HMM algorithm, we introduce the filtered effective ensemble size ν HMM (t) [START_REF] Liu | Nonparametric hierarchical Bayes via sequential imputations[END_REF]:

ν HMM (t) := 1 80 i=1 p HMM (i, t) 2 , ( 23
)
and we define equivalently ν HMM (t). These quantities are estimates of the number of ensemble members that can be retained according to surge observations, assuming one discards very unlikely members.

In Fig. 7, variables δµ HMM , δµ HMM /σ 20CR and ν HMM are shown as a function of time for the period 1846-1890. All these quantities show a strong seasonality. This is due to a much stronger SLP variability in winter, and a corresponding stronger response of the surges. The figure shows that the amount of correction δµ HMM and the decrease in ensemble size ν HMM are much stronger using smoothed probabilities with HMM rather than probabilities without HMM. Showing the deviation δµ HMM in the Bay of Biscay, where the standard deviation of δµ HMM /σ HMM is strongest (see Fig. 8), substantial absolute values of ∼600 Pascals are obtained in early 1850s winters, even after averaging over 3 months. These large deviations correspond to more than one standard deviations of the ensemble size. Using probabilities without HMM, deviations are weaker but still non-negligible (∼ 500P a, ∼ 0.7σ). The slow decrease in δµ HMM with time is coherent with slowly increasing observations used in 20CRv3, although with substantial decadal variations. However, δµ HMM /σ HMM and δµ HMM /σ HMM do not show a clear trend, indicating a persisting gain in information from surge observations throughout the 19th century.

In terms of effective size, Fig. 7 shows that the smoothing HMM algorithm imposes a strong member selection, with mostly only 1 member retained at each time step, in winter and before 1880. Probabilities without HMM mostly retain more than half of the members, although peak low values of ν HMM (t) show that even without the HMM sometimes more than half of the ensemble members are highly unlikely. Filtered effective ensemble size reaches very low yearly and seasonal average values, indicating that many 20CR members are highly unlikely with respect to surge estimates from tide gauge observations. A strong increase in ν HMM (t) is witnessed around year 1880. This can be explained by the availability of a large number of weather station data in Eastern Europe and Russia from 1880-on, and by an intensification of maritime traffic around 1880.

The spatial structure of δµ is examined in Fig. 8. The analysis of time-standard-deviation of δµ HMM and δµ HMM /σ 20CR

shows that the area of greatest influence of the corrections from surge-smoothing from Brest and Saint-Nazaire tide gauges is in the Bay of Biscay. This can be explained by the passage of strong storms in the Bay of Biscay, which can cause high surges in Brest and Saint-Nazaire, and by the sparsity of direct pressure measurements (ship logs) in this area in the 19th century.

Standard deviation of δµ HMM shows largest values to the north-west of the map, which is where strong storms travel. Indeed, These corrections also have a strong decadal variation, with non-trivial yearly averages persisting for several years, as shown in Fig. 9. The same behaviour can be witnessed for the surge, which is strongly anti-correlated to these deviations (Fig. 9). This can be explained by the fact that 20CRv3 smooths SLP values in areas of sparse measurements, and that surge-filtering corrections allow to retrieve more realistic intense values (either positive or negative). This interannual variability is related to the variability in storminess [START_REF] Bärring | Multi-indices analysis of southern Scandinavian storminess 1780-2005 and links to interdecadal variations in the NW Europe-North Sea region[END_REF].

Focus on two 19th-century storms

One of the aim of this study is to show that old tide gauge data can be used to better understand past severe storms. In this section, two known storms are studied for illustrative purposes. The "Moray Firth fishing disaster" of august 1848 was the result of a storm hitting the Northeastern coast of Scotland [START_REF] Coull | The Role of the Fishery Board in the Development of Scottish Fishing Harbours c. 1809-1939[END_REF]. The day before the storm, conditions seemed favorable for fishing, therefore the unexpected storm caused more than 100 casualties among Scottish fishermen, which later led to improvements in harbours and vessels in Scotland. This severe storm, although located hundreds of kilometers away from Brest, left a trace on the sea-level at the Brest tide gauge. Our estimation of surge shows a 12h-averaged positive surge of 15cm on the night between the 18th and 19th of august (Fig. 10). Probabilities without HMM give an effective ensemble size down to ν HMM ≈ 30 at 9am on august 19th, and smoothed probabilities give ν HMM ≲ 1.7 until 12am on that same day, indicating a strong selection among candidate SLP patterns given by members of the 20CRv3 reanalysis. The retained members give an average storm of greater intensity and spatial extent over the British Isles than using the uniform ensemble average from 20CRv3.

Although these results do not allow to assert firm statements about this storm, they show that even tide gauges located hundreds of kilometers away from the core of a pressure system can bring interesting pieces of information. Germany [START_REF] Zimmerli | Winter storms in Europe: messages from forgotten catastrophes[END_REF]. As shown in Fig. 11, this storm was more intense than the one from Fig. 10, and hit directly the Northern coast of France, leaving a strong footprint on Brest and Saint-Nazaire tide gauges. Peak positive 12haveraged surges were of more than 50cm for Brest and close to 80cm in Saint-Nazaire. At this moment even ν HMM shrinks down to ≈1.9, while ν HMM ≈ 1 more than two days before the extreme surge event. Therefore, the filtered averages on the left correspond to one single member from 20CRv3. Examining the difference µ HMM on the 11th (12am) and the following night (12th, 00am), it seems that the spatial extent of the storm off the Western coast of France was underestimated in most 20CRv3 members, according to surge observations in Brest and Saint-Nazaire. It shows again that tide gauges could be used to better estimate past atmospheric states, and in particular severe storms. This work has several potential applications. First, replicating this work with other tide gauges could help to validate reanalyses like 20CRv3 against independent data, and to potentially identify anomalous trends or wrong estimation of specific events. Second, tide gauge data could be directly used to constrain atmospheric reanalyses, for instance as "inverse barometers", only with higher uncertainties than direct pressure sensors. Tide gauges which are more sensitive to wind than pressure, such as the ones in the North Sea, would require a more complex treatment, but could still be used to constrain atmospheric reanalyses. Third, tide gauges could be used to constrain regional scale atmospheric simulations in order to better estimate the magnitude and spatial extent of known past severe storms. Fourth, tide gauge records could be combined with direct observations of atmospheric pressure to give statistical estimates of atmospheric fluctuations in the 19th century without the use of a NWP model, such as the optimal interpolation of [START_REF] Ansell | Daily mean sea level pressure reconstructions for the European-North Atlantic region for the period 1850-2003[END_REF] based on direct pressure observations only, or the analogue upscaling of [START_REF] Yiou | Ensemble meteorological reconstruction using circulation analogues of 1781-1785[END_REF] for the short period 1781-1785 of dense observations in western Europe. Finally, this work could be replicated in a more general context, using other types of variables and observations, learning the relationship between observations and large-scale features using recent observations and precise reanalyses, and applying these statistical relationship in the past to uncover past large-scale events.

Figure 1 .

 1 Figure 1. Number of surface pressure observations from the International Surface Pressure Databank (ISPD) assimilated in the Twentieth century Reanalysis version 3 (20CRv3). Top: Monthly in 2 • ×2 • box centered on Brest and Saint-Nazaire, for years 1870 (a) and 2000 (b). Bottom: Yearly average of daily number of observation in 1870 (c) and 2000 (d).

Figure 2 .

 2 Figure 2. Brest and Saint-Nazaire tide gauge record availability trough time, from the GELSA-3 database (Haigh et al., 2023).

Figure 3 .

 3 Figure 3. (a) Dots: Yearly averaged sea-level in Brest (blue) and Saint-Nazaire (orange). Dashed lines: constant (before 1900) and linear (after 1900) trends from least-squares fit. To improve readability, a constant +30cm was added to the Saint-Nazaire sea-level in this figure. (b) An example of detrended sea-level and associated storm surge, along with the 12-hours moving average of the surge (Brest tide gauge, March 1963).

5

  https://doi.org/10.5194/egusphere-2023-2997 Preprint. Discussion started: 2 January 2024 c Author(s) 2024. CC BY 4.0 License.

  ages and variances given by the linear regressions estimated on the satellite-era ensemble-average SLP anomaly and respective gradients at Brest and Saint-Nazaire. It is implicit in the above formulae that s B and s SN are univariate variables conditioned by the values of SLP and SLP gradients. For a given map of SLP, the surges in Brest and Saint-Nazaire are assumed to be normal random variables. 6 https://doi.org/10.5194/egusphere-2023-2997 Preprint. Discussion started: 2 January 2024 c Author(s) 2024. CC BY 4.0 License. Coefficients of the linear regressions with the surges of Brest and Saint-Nazaire as explained variables and the SLP (α) and its horizontal gradients (β and γ) as explanatory variables. The linear regressions are computed on the satellite-era (1981-2015) ensembleaverage from 20CrRv3. The results of one-dimensional linear-regressions (1DLR) are shown for comparison.

  Fig 4(a,b) shows scatter-plots of the data used to estimate the LR around the LR mean values.

Figure 4 .

 4 Figure 4. Two-dimensional histograms of surges in (a) Brest and (b) Saint-Nazaire, versus linear fit using ensemble-averaged sea-level pressure anomaly SLP, longitudinal and latitudinal gradients (at ±2 • ) of SLP at the city of interest. The predicted surge is αSLP(t) + β∆ lat SLP(t) + γ∆ lon SLP(t) with coefficients differing for Brest and Saint-Nazaire and given in Table 1. (c) Black arrows: direction of geostrophic wind projection associated with the coefficients facing ∆ lon SLP and ∆ lat SLP, leading to an increase in sea-level (surge).Arrows' lengths are proportional to the LR coefficients, thus large arrows indicate a stronger sensitivity of the sea-level to the winds.

  small systematic 8 https://doi.org/10.5194/egusphere-2023-2997 Preprint. Discussion started: 2 January 2024 c Author(s) 2024. CC BY 4.0 License.

Figure 5 .

 5 Figure 5. Histograms of Brest surge vs. linear combination of SLP and horizontal gradients, with the same coefficients for all sub-figures (ad). These coefficients are taken from the satellite-era LR (Fig. 4a and Table 1). (a) SLP is taken from the 20CRv3-members-average in the satellite-era as in Fig. 4a. (b) SLP is taken from the 20CRv3-members-average from years 1864 to 1896. (c) SLP is taken from the first member of 20CRv3 from years 1870 to 1896. (d) SLP is taken from the second member of 20CRv3 from years 1870 to 1896.

195

  The histograms of Fig 5(c,d) show a larger spread around the diagonal, and no clear sign of bias. This can be interpreted as the result of larger uncertainties in 1864-1896 regarding the value of the SLP in Brest due to a lower observation density. This added uncertainty leads to a higher spread in the linear relationship between SLP and surge. Table2 indicatesthat the standard 9 https://doi.org/10.5194/egusphere-2023-2997 Preprint. Discussion started: 2 January 2024 c Author(s) 2024. CC BY 4.0 License. deviation in time of individual members' SLP is consistent between 1864-1896 and 1981-2015 (less than 0.1% difference),

  11 https://doi.org/10.5194/egusphere-2023-2997 Preprint. Discussion started: 2 January 2024 c Author(s) 2024. CC BY 4.0 License.

  13 https://doi.org/10.5194/egusphere-2023-2997 Preprint. Discussion started: 2 January 2024 c Author(s) 2024. CC BY 4.0 License. a i (1) = p HMM (i, 1) ,

Figure 6 .

 6 Figure6. Log-likelihood of Brest surges (dotted blue line) as a function of parameter Θ defined in Eq. (11). For comparison, the log-

Figure 7 .

 7 Figure 7. (a) Average SLP deviation δµHMM in Pascals, in the center of the Bay of Biscay. (b) Same as (a) normalized by reanalysis ensemble standard deviation δµHMM/σ20CR. (c) Effective ensemble size νHMM.Orange: using smoothed probabilities with HMM pHMM(i, t). Blue: using probabilities without HMM pHMM(i, t). Bold: yearly average. Thin: 3-month average. All plots make simultaneous use of data from Brest and from Saint-Nazaire tide gauges.

Figure 8 .

 8 Figure 8. Time-standard deviation of δµHMM (a,d), δµHMM (b,e) and δµHMM/σ20CR (c,f) computed from October to March, for years 1848-1880 (a,b,c) and 1880-1890 (d,e,f). Probabilities pHMM and pHMM make use of data from both Brest and Saint-Nazaire tide gauges.

Figure 9 .

 9 Figure 9. Orange, bold: yearly-average of δµHMM at position 47.83 • lat, -7.57 • lon, in the Bay of Biscay, using data from both Brest and Saint-Nazaire tide gauges. Black, thin: Brest surge (inverted sign).

Figure 10 .

 10 Figure 10. Mean result of the smoothing procedure during the Moray Firth Fishing disaster that happened during the night between august 18th and 19th, 1848. Left column: average according to surge observations and HMM smoothing algorithm (probabilites pHMM(i, t)). Middle column: average using constant uniform weights on 20CRv3 members. Right column: difference HMM -uniform. Right panel: surge observations and effective ensemble sizes. The dates of the left-hand-side plots are indicated with dots on the right-hand-side figure.

Figure 11 .

 11 Figure 11. Same as Fig. 10 a storm that hit western Europe around the 12th of March 1876, later called Lothar's big brother. In this figure tide gauge observations from Saint-Nazaire are also available.

  doi.org/10.5194/egusphere-2023-2997 Preprint. Discussion started: 2 January 2024 c Author(s) 2024. CC BY 4.0 License. This study is a proof of concept for the use of century-old tide gauge data as a means of understanding past atmospheric subseasonal variability. Storm surges of Brest (starting in 1848) and Saint-Nazaire (starting in 1863) allow to assess part of the atmospheric variability that was uncaught in global 20CR reanalyses based on pressure observations only. Weighing 20CR members according to surge observations reduces the effective ensemble size, and implies significant deviations in membersaveraged sea-level pressure in the Bay of Biscay. Through the second-half of the 19th century, these deviations diminish and the effective ensemble size rises, however they remain non-negligible.

Table 2 .

 2 Standard deviations of SLP and horizontal gradients at the city of Brest, during the two different time periods used in Fig 5, expressed in Pascals. For the period 1981-2015, ensemble average and individual members show no significant difference. Relative differences with respect to 1981-2015 are expressed in percentages.

	1981-2015 (ensemble average)	949.1	259.9	486.1
	1864-1896 (ensemble average)	911.4 (-4.0%)	242.3 (-6.8%)	466.7 (-4.0%)
	1864-1896 (individual members) 948.7 (-0.04%)	292.1 (+12.4%)	517.2 (+6.4%)

std

(SLPBrest)

std(∆ lon SLPBrest) std(∆ lat SLPBrest)

https://doi.org/10.5194/egusphere-2023-2997 Preprint. Discussion started: 2 January 2024 c Author(s) 2024. CC BY 4.0 License.

Acknowledgements. This work was supported by European Research Council (ERC) Synergy grant STUOD -DLV-856408. Support for the Twentieth Century Reanalysis Project version 3 dataset is provided by the U.S. Department of Energy, Office of Science Biological and Environmental Research (BER), by the National Oceanic and Atmospheric Administration Climate Program Office, and by the NOAA Physical Sciences Laboratory. 21 https://doi.org/10.5194/egusphere-2023-2997 Preprint. Discussion started: 2 January 2024 c Author(s) 2024. CC BY 4.0 License.

Author contributions. Paul Platzer wrote and run the codes, created the figures, and wrote the original manuscript. Pierre Tandeo proposed major modifications to the manuscript. Pierre Ailliot helped design the Hidden Markov Model algorithm. All authors helped conceptualizing the work, reviewed the manuscript and approved its final form.

Competing interests. The authors declare that they have no competing interests.