Large-scale changes of the semidiurnal tide along North Atlantic coasts from 1846 to 2018

Abstract. We investigated the long-term changes of the principal tidal component M2 along North Atlantic coasts, from 1846 to 2018. We analysed 18 tide gauges with time series starting no later than 1940. The longest is Brest with 165 years of observations. We carefully processed the data, particularly to remove the 18.6-year nodal modulation. We found that M2 variations are consistent at all the stations in the North-East Atlantic (Cuxhaven, Delfzijl, Hoek van Holland, Newlyn, Brest), whereas some discrepancies appear in the North-West Atlantic. The changes started long before the 20th century and are not linear. The secular trends in M2 amplitude vary from one station to another; most of them are positive, up to 2.5 mm/yr at Wilmington since 1910. Since 1990, the trends switch from positive to negative values in the North-East Atlantic. Concerning the possible causes of the observed changes, the similarity between the North Atlantic Oscillation and M2 variations in the North-East Atlantic suggests a possible influence of the large-scale atmospheric circulation on the tide. Our statistical analysis confirms large correlations at all the stations in the North-East Atlantic. We discuss a possible underlying mechanism. A different spatial distribution of mean sea level (corresponding to water depth) from one year to another, depending on the low-frequency sea-level pressure patterns, could impact the propagation of the tide in the North Atlantic basin. However, the hypothesis is at present unproven.



Introduction
Tides have been changing due to non-astronomical factors since the 19th century (Haigh et al., 2019;Talke and Jay, 2020). In the North Atlantic, secular variations have been observed at individual tide gauge stations, e.g. Brest (Cartwright, 1972;Pouvreau, 2008), Newlyn (Araújo and Pugh, 2008;Bradshaw et al., 2016), New York (Talke et al., 2014), and Boston , but also at regional scale, e.g. Gulf of Maine (Doodson, 1924;Godin, 1995;Ray, 2006;Ray and Talke, 2019), at the North Atlantic basin scale (Müller, 2011), and at a quasi-global scale (Woodworth, 2010;Müller et al., 2011;Mawdsley et al., 2015). Long-term changes in tidal constituents are rather small at some coastal stations but tend to be statistically significant. The order of magnitude of these changes varies spatially and may reach a few centimetres per century for M 2 amplitude. For example, Ray and Talke (2019) found trends varying from −1 to 8 cm per century in the Gulf of Maine over the last century. Woodworth et al. (2010) and Müller et al. (2011) found trends of a few percent per century in the Atlantic. The changes can be larger in many estuaries and rivers (Talke and Jay, 2020).
The physical causes of these changes can be multiple and difficult to disentangle. In particular, the complexity comes from the possible interaction between local and largescale causes. Changes may have a local-scale origin, such as changes in the nearby environment (e.g. harbour development, deepening of channels, dredging, siltation) or changes in the instrumentation (e.g. tide gauge technology, observatory location, instrumental errors). For example, Familkhalili and Talke (2016) show that mean tidal range at Wilmington has doubled since the 1880s, due to channel deepening in the Cape Fear River estuary. Changes may also have a largescale origin, i.e. regional or global. Haigh et al. (2019) reported several possible large-scale mechanisms: (1) tectonics and continental drift, (2) water depth changes due to mean sea level rise or geological processes such as the Earth's surface glacial isostatic adjustment Pickering et al., 2017;Schindelegger et al., 2018), (3) shoreline position, (4) extent of sea-ice cover , (5) sea-bed roughness, (6) ocean stratification which may modify the internal tides and bottom friction over continental shelves (Müller, 2012), (7) non-linear interactions, and (8) radiational forcing (Ray, 2009).
Several authors have explored mean sea level (MSL) rise as a potential mechanism to explain M 2 changes. For example, simulations by Pickering et al. (2012) show that a 2 m sea level rise could modify M 2 from −20 to 20 cm around the whole ocean. Idier et al. (2017) show that depending on the location, the changes can account for ±15 % of the regional sea level rise. Schindelegger et al. (2018) find changes of about 1 %-5 % of the sea level rise. Beyond MSL rise, other mechanisms have been explored to explain M 2 changes. For example, Colosi and Munk (2006) attribute the changes of M 2 amplitude at Honolulu, Hawaii, to a 28 • rotation of the internal tide vector in response to ocean warming. Ray and Talke (2019) suggest that long-term changes in stratification could play a role in the Gulf of Maine. Müller (2011) suggests a possible link between M 2 changes and atmospheric dynamics in the North Atlantic; he reported that the time series of the North Atlantic Oscillation (NAO) show similar characteristics to those of the tidal amplitudes and phases. In the Gulf of Maine, Pan et al. (2019) suggest that changes in the response of the nodal modulation of the M 2 tide from 1970s to 2013 may be linked with the NAO. In Southeast Asian waters, Devlin et al. (2018) show that the impact of atmospheric circulation (via the wind stress, through Ekman current) on the M 2 seasonal cycle may be significant and comparable to the effect of permanent (geostrophic) currents. In the North Sea, Huess and Andersen (2001) explain a large part of M 2 seasonal cycle by the role of atmospheric dynamics, whereas Müller et al. (2014) and Gräwe et al. (2014) suggest a major role of the thermal stratification. These examples show the diversity of mechanisms that play a role in tide changes. In the present paper, we focus on the role of MSL and atmospheric dynamics. This paper has two main objectives. The first is to characterize the secular changes of the M 2 tide over the North Atlantic. We focus on the longest time series, i.e. starting no later than 1940. This approach is complementary to previous studies investigating M 2 changes focusing on smaller spatial scales, e.g. Brest Pouvreau, 2008), Gulf of Maine (Ray, 2006;Ray and Talke, 2019), or focusing on shorter temporal scales, i.e. recent decades (Woodworth, 2010;Müller, 2011). The second objective is to detect if there is any large-scale coherence in the observed changes in the North Atlantic and investigate the possible link with the at-mospheric circulation, already mentioned by Müller et al. (2011), on the basis of qualitative criteria. Here, we further provide quantitative insights into the possible influence of the NAO and discuss a possible NAO-related climate mechanism that can partly explain the observed changes.
The paper is organized as follows. The first section below describes the data: the sea level data (i.e. tide gauges and their processing) and the atmospheric data (i.e. climate indices and sea level pressure data). The following section presents the results (i.e. M 2 variations and trends). We then discuss a possible link between the observed tidal changes and MSL, as well as climate indices. The tide gauge data were retrieved from the University of Hawaii Sea Level Center (UHSLC, website accessed April 2020). The dataset consists of 249 stations in the Atlantic Ocean, with hourly sea level observations. Two additional long-term stations -Delfzijl and Hoek van Holland -were provided by Rijkswaterstaat (RWS) in the Netherlands.
We selected the stations following three criteria: time series (1) starting before 1940, (2) with at least 80 years of data, and (3) with tidal amplitude significant enough to detect trends, i.e. M 2 amplitude larger than 10 cm. Note that we selected only years with at least 75 % of data (see Sect. 2.1.2). Only 24 stations among the 249 from UHSLC fulfilled the two first criteria (Fig. 1). They are all located in the Northern Hemisphere. On the east side of the North Atlantic, Stockholm, Gedser, Hornbaek, Tregde and Marseille were discarded due to too small of an M 2 amplitude (i.e. lower than 10 cm). These stations are located in the Baltic Sea (Stockholm, Gedser), in the strait separating the Baltic and the North Sea (Hornbaek), in the North Sea (Tregde), and in the Mediterranean Sea (Marseille). On the west side of the North Atlantic, Galveston, Pensacola and Cristobal were also discarded due to too small of a tidal amplitude (i.e. lower than 10 cm). These stations are located in the Gulf of Mexico (Galveston, Pensacola) and the Caribbean Sea (Cristobal).
Finally, 18 stations followed the three criteria detailed above and were selected for this study (see stations in bold in Fig. 1, 16 stations are from UHSLC, and 2 from RWS). Among them, 5 are located on the North-East Atlantic coasts (Newlyn, Brest, Hoek van Holland, Delfzijl and Cuxhavennote that Hoek van Holland, Delfzijl and Cuxhaven are located in the North Sea) and 13 are located on the North-West Atlantic coasts (Halifax, Eastport, Portland, Boston, Newport, New London, New York, Atlantic City, Lewes, Wilmington, Charleston, Fort Pulaski and Key West). The main characteristics of the 18 selected stations are summarized in Table 1

Data processing
Harmonic analysis was performed in order to compute the M 2 amplitude. We used the MAS program (Simon, 2007(Simon, , 2013, developed by the French Hydrographic Office (SHOM). This program gives results similar to the T_Tide harmonic analysis toolbox (Pawlowicz et al., 2002). For instance,  found no differences in the yearly amplitudes of M 2 at Brest over the period 1846 to 2005 using either T_Tide or MAS. Hourly time series were analysed yearly. Note that at Delfzijl and Hoek van Holland, data had to be interpolated every hour before 1970, as the temporal sampling was 3 h. (We checked with hourly time series from recent years (1971-2018) that 3-hourly sampling did not result in a significant reduction of M 2 amplitude in a tidal analysis compared to hourly sampling.) We processed only years with at least 75 % of data, to avoid seasonal modulation affecting the computed amplitudes. In the North Atlantic, M 2 is affected by a seasonal variation of a few percent (Pugh and Vassie, 1976;Huess and Andersen, 2001;Müller et al., 2014;Gräwe et al., 2014). Considering only years with at least 75 % of data resulted in excluding up to 15 years for a given station (Table 1, columns 3 and 4). We care-fully removed the nodal modulation of M 2 amplitude (Simon, 2007(Simon, , 2013, as described briefly in Appendix A. Finally, 3 station-years were discarded due to problems in the record (1953( and 1962( at Delfzijl, 1953 Hoek van Holland), and 2 more station-years due to doubtful M 2 values (1972at Eastport, 1978.
At all the stations, we computed the normalized M 2 amplitude, removing the average and dividing by the standard deviation over the period 1910-2010: 1910,2010] σ M 2 [1910,2010] . (1) The average M 2 and standard deviation σ M 2 over the 1910-2010 period are given in Table 1 (column 5). The idea is to scale the data in order to compare all the stations together.

Climate indices
We investigated the correlation between secular changes in the tide and climate indices, such as the North Atlantic Oscillation (NAO) or the Arctic Oscillation (AO) -also called Northern Annular Mode (NAM) (Hurrell, 1995;Hurrell and Deser, 2009;. These climate indices are related to the distribution of atmospheric masses. They are based on the difference of average sea-level pressure between two centres of actions (i.e. stations) over long periods (e.g. monthly, seasonal, annual). The NAO is the major pattern of weather and climate variability over the Northern Hemisphere (Hurrell, (Hurrell et al., 2003). We used the wintertime (December to March) Hurrell station-based NAO Index (retrieved from https://climatedataguide.ucar.edu/climate-data/hurrellnorth-atlantic-oscillation-nao-index-station-based, last access: April 2020). It is based on the difference of normalized average winter sea-level pressure between Lisbon (Portugal) and Stykkishólmur/Reykjavik (Iceland). The normalization involves removing the mean  and dividing by the long-term standard deviation. The NAO index covers the period 1864-2019. The Arctic Oscillation (AO) is another index which resembles the NAO index. It is defined as the first EOF of Northern Hemisphere winter sea-level pressure data Wallace, 1998, 2000;. The AO index is highly correlated with the NAO. We used the wintertime Hurrell AO index (retrieved from https://climatedataguide.ucar.edu/climate-data/hurrellwintertime-slp-based-northern-annular-mode-nam-index, last access: April 2020). The AO index covers the period 1899-2019.
To remove the interannual variability and estimate lowfrequency variations, climate indices were low-pass filtered with a 9-year mean filter.

Sea level pressure
We employed the Twentieth Century Reanalysis (20CR version 3 dataset) (Compo et al., 2011;Slivinski et al., 2019), a historic weather reconstruction from 1836 to 2015, with a 1 • gridded global coverage. However, we made use only of data from 1850 to be more consistent with the temporal coverage of the tide gauge measurements. This will be discussed in Sect. 4.

M 2 variations
For the North-East Atlantic, the variations of normalized M 2 amplitude are presented in Fig. 2a.
The first result is that since 1910, the variations show similar patterns at all the stations; M 2 amplitude decreases up to the 1960s, then increases, and decreases again since the 1990s. This suggests that these changes are probably due to large-scale processes, rather than local effects due to changes in the environment (e.g. harbour development, dredging, siltation) or instrumentation errors. The similar patterns between Brest and Cuxhaven may be surprising, as Cuxhaven is located in the North Sea, and not in the open Atlantic Ocean, and far away from Brest, around 1300 km. This indicates that the spatial scale of the processes responsible for these changes must be at least as large as the North-East Atlantic.
. Different authors have noticed the increase of tidal range from 1960 to 1990 in the southern North Sea. Hollebrandse (2005) found a gradual increase during the period 1955-1980 at all the stations of the Dutch coast (five stations including Hoek van Holland) and the German coast (seven stations). Mudersbach et al. (2013) found a significant increase in M 2 amplitude at Cuxhaven since around the mid-1950s. Note that Cuxhaven is located in the German Bight; shallow depths and the shape of the coastline may induce some amplification. Variations in M 2 at Cuxhaven could therefore be sensitive to local effects, such as the migration of the underwater channels and the evolution of the tidal flats (Jacob et al., 2016). Moreover, Cuxhaven is located in the Elbe estuary, and some river engineering works, such as narrowing and deepening, may induce tidal amplification . Before 1910, normalized M 2 values are higher at Brest than at Delfzijl. The construction of dykes that have gradually closed the harbour of Brest since the end of the 19th century may have altered the tide at Brest. The high values before 1910 may be due to local changes, in addition to large-scale changes. To go further, the potential role of these successive constructions needs to be investigated (Wikipedia contributors, 2020). Cartwright (1972) made a first attempt to evaluate the influence of reducing the width of access to the harbour but did not take into account a potential role of dredging, for which we have no information. This example underlines the complexity of interpretation of the variations when changes of local and large-scale origin occur at the same time. Note that in the following, we focus mainly on the 20th century, as most of the stations start after 1900 (15 out of 18 stations).
The second result is that there is no obvious linear trend in M 2 variations, but rather break or change points, M 2 increasing and then decreasing, depending on the periods considered. Overall, M 2 decreases from 1910 until 1960, increases again until 1980-1990, to finally decrease since 1990; note that the curve flattens between 1920 and 1940.  already noticed these variations at Brest and Newlyn and suggested a long-period oscillation of around 140 years, rather than a steady secular trend. A careful analysis of the harmonic development of the tidal potential showed that no tidal component could explain this oscillation. Similarly, no linear combination of tidal harmonic components could explain it . This indicates that these variations are not due to an astronomical component. However, in contrast to Brest, M 2 at Delfzijl stays flat between 1880 and 1920. The decrease observed at Brest between 1880 and 1920 may be due to harbour development and/or dredging (see above). This underlines the importance of sea level data archaeology, for research studies related to long-term changes (Pouvreau, 2008;Woodworth et al., 2010;Marcos et al., 2011;Jay, 2013, 2017;Ray and Talke, 2019;Bradshaw et al., 2015Bradshaw et al., , 2020. The third result is that changes in M 2 have not the same order of magnitude at each station (see Fig. B1 in Appendix B for time series of M 2 ). Note that Fig. 2 represents normalized M 2 , i.e. removing the average and dividing by the standard deviation. The order of magnitude of unnormalized M 2 changes is roughly the same at Brest and Newlyn (standard deviations of 0.9 and 0.8 cm respectively, Table 1, column 5), but more than three times larger at Cuxhaven (standard deviation of 3.7 cm), and even larger at Delfzijl (standard deviation of 7 cm). This suggests that the North Sea may be more sensitive to the processes responsible for these changes. Note also that the environmental setting of Cuxhaven and Delfzijl in the Elbe and Ems estuaries, respectively, could introduce some amplification .
For the North-West Atlantic, the variations of normalized M 2 amplitude are presented in Fig. 2b and c. The first feature is that M 2 amplitude varies differently in the North-West and in the North-East Atlantic. The second is that there are discrepancies between stations, even when close to each other (e.g. Atlantic City and Lewes). We split the stations into two groups, in order to facilitate the detection of patterns, each being consistent in terms of trends: one with positive trend (group 1 in Fig. 2b), the other one with negative or no trend (group 2 in Fig. 2c).
The first group (with positive trends) consists of nine stations (Fig. 2b). Three outcomes can be highlighted. The first is that M 2 amplitude has increased overall since 1900. However, between 1980 and 1990, all the stations slightly decrease, and since 1990 they have increased again. The second outcome is that the rate of increase is very different from one station to another (keeping in mind that M 2 is normalized by standard deviation in Fig. 2). Portland is increasing 1.4 times faster than Charleston (standard deviations being respectively of 1.82 and 1.33 cm) and 28 times faster than Key West (standard deviation being only 0.36 cm at Key West). The large increase in Portland may be explained by some amplification in the Gulf of Maine. In many semienclosed basins, resonance leads to tidal amplification (Talke and Jay, 2020;Haigh et al., 2019). In the Gulf of Maine, Ray and Talke (2019) reported that the tides in the Gulf are in resonance, with a natural resonance frequency close to the N 2 tide (Garrett, 1972;Godin, 1993). Tides may be then very sensitive to any changes in the environment (e.g. basin configuration -shape, depth -but also external forcing). The third outcome, and probably the most interesting one, is related to the values of M 2 at Eastport, Portland and New York in the 1860s, estimated from Ray and Talke (2019) and Talke et al. (2014), and represented (after normalization) as stars in Fig. 2b. These values are not consistent with the positive linear trends observed since 1900, which provides some consistency with the hypothesis formulated from the analyses of the data prior to the 20th century in Fig. 2a: long-term variations introduce some breaks or change points, M 2 increasing and then decreasing, depending on the periods considered. The decrease observed between the 1870s and 1920s at the four stations (Brest, Eastport, Portland, New York) suggests a possible large-scale signal, in addition to local processes.
The second group (with negative or no trend) consists of four stations (Fig. 2c). Two points can be highlighted. The first is that M 2 decreases overall for Halifax, Newport and Lewes. This is less clear for Atlantic City, which is quite noisy and shows no significant trend. The second point is that at Halifax, M 2 values in 1896-1897 are higher than those after 1920. This suggests that the decrease may have started before the 20th century. However note that at Halifax, there is a long gap in the data recording (1898-1919), which raises the possibility of an instrumentation origin in the observed decrease of the M 2 amplitude.

Estimated trends
We estimated the trends for M 2 amplitude at each station, using linear regression. We computed the trends over two periods: 1910-2018, which corresponds roughly to the whole period of data (only five stations start before 1910), and 1990-2018, which corresponds to recent decades. Some tests showed that the later results were not very sensitive to the start date (moving 1990 to 1985 or 1995). The trend uncertainties were estimated considering the noise content in the time series using SARI software (Santamaría-Gómez, 2019). The noise was modelled as a white plus power-law noise, whose spectral index was found to be close to −1 (flicker noise). The results are summarized in Table 1 (columns 7 and 8) and Figs. 3 and 4.
The trends estimated since 1910 vary significantly from one station to another (Fig. 3). They are positive overall (up to 2.5 mm/yr at Wilmington), which is consistent with previous findings (Araújo and Pugh, 2008;Ray, 2009;Woodworth, 2010;Müller et al., 2011;Ray and Talke, 2019). They are slightly negative at three stations (Halifax, Newport, Lewes), and one station shows no trend (Atlantic City). The estimates are statistically consistent with those found previously by different authors (e.g. 0.14 ± 0.09 mm/yr at Newlyn compared to 0.19 ± 0.03 mm/yr in Araújo and Pugh (2008), 0.56 ± 0.06 mm/yr in Portland, compared to 0.59 ± 0.04 mm/yr in Ray and Talke, 2019). Note that our error bars are larger, because we considered the noise content in the time series as a white noise plus power law noise (we obtained the same error bars considering white noise only). In the North-East Atlantic, the trends are consistent with each other (in terms of sign), which is not surprising as the stations vary similarly (Fig. 2a).
The largest trends since 1910 are mainly observed in semiclosed basins: Wilmington in the Cape Fear River estuary, Delfzijl in Ems estuary, Cuxhaven in Elbe estuary, and Eastport and Portland in the Gulf of Maine. This suggests a possible amplification due to resonance effects (e.g. Gulf of Maine) and/or propagation in shallow waters (e.g. Cuxhaven), in addition to local effects. The stations located in estuaries or in a harbour with a channel may have been subject to dredging. Channel deepening increases the water depths, which reduces the effective drag and leads to tidal range amplification. This effect may be particularly large in estuaries (Ralston et al., 2019;Talke and Jay, 2020) and may explain the larger trends at Wilmington (Familkhalili and Talke, 2016) and Delfzijl. Finally, the shifting locations of amphidromic points could also play a role (Haigh et al., 2019). In the North Sea, different authors show a possible migration of the present-day amphidromes, under a 2 m sea-level rise scenario (Pickering et al., 2012;Idier et al., 2017).
The trends estimated since 1990 are quite different from those estimated since 1910 (Figs. 3 and 4), with more stations with negative trends: 9 stations out of 18 have post-1990 negative trends, whereas only 3 stations out of 18 have post-1910 negative trends (Table 1, columns 7 and 8). In the North-East Atlantic, they all switch from positive to negative trends. This underlines (1) some spatially coherent changes in recent decades (Müller, 2011;Ray and Talke, 2019) and (2) the difficulty in estimating long-term trends from short records (i.e. less than 30 years), especially if the data are noisy (interannual variability) and the underlying processes non-linear (change points).
The trends have to be interpreted very carefully as the M 2 variations are not linear and may increase or decrease depending on the years; as a consequence, the estimated trends depend strongly on the period considered to estimate it. The interannual variability also plays an important role, and when substantial, trends can vary depending on the computational period. For example, at Cuxhaven, the large interannual variability leads to a large uncertainty on the trend computed since 1990 (−0.47 ± 0.78 mm/yr).

Possible link with mean sea level rise
The MSL rise could partly explain M 2 changes. Simulations show that MSL rise can result in an change of M 2 up to ±10 % of the rise (Pickering et al., 2017;Idier et al., 2017;Schindelegger et al., 2018). Schindelegger et al. (2018) show that the sign of the observed M 2 trend is correctly reproduced at 80 % of the tide gauges on a global scale, but their simulated trends tend to differ from observations by a factor of 3 to 5; i.e. their simulations underestimate the M 2 response to MSL rise in terms of magnitude. Schindelegger et al. (2018) conclude that "magnitudes of observed and modeled M 2 trends are within a factor of 4 (or less) from each other in nearly 50 % of the considered cases". The large discrepancies between the simulations and the observations strongly suggest that MSL rise is not the only process that may explain M 2 changes -other large-scale processes, in addition to local processes, may also play a role. Figure 5 shows the annual MSL, after removing the average over the period 1910-2010 and filtering with 9-year time windows. The correlations between M 2 and MSL indicate that M 2 varies strongly with MSL (see Sect. 4.2). However, M 2 variations show some variability in the North-East Atlantic (Fig. 2a), which may not be explained with MSL rise alone.

Possible link with MSL and climates indices
Processes other than MSL rise may impact the tide (see Sect. 1), such as the atmospheric circulation and the ocean stratification. Ocean and atmosphere are fully coupled, and air-sea fluxes are responsible for the exchange of momentum, water (evaporation and precipitation budget) and heat at their interface. Among the wide range of possible interactions, two mechanisms have been explored for their ability   to modify the tide: (1) the momentum flux (wind stress) and the gradient of sea level pressure which act on the barotropic tide and (2) the water and heat fluxes which induce changes in both temperature and salinity distribution in the ocean. The latter effect acts on the stratification, which in turn could impact the tide in two different ways. The first way is the internal tide generation which transfers energy from barotropic and baroclinic motion and modifies surface tidal expression (Colosi and Munk, 2006). However, in the present study, most of the observations come from coastal stations sheltered by wide continental shelves which dampen internal waves. More important is the second way: the stratification acts on the eddy viscosity profile by modifying current profiles and bottom drag over continental shelves, which in turn modifies the M 2 surface expression (Kang et al., 2002;Müller, 2012;Katavouta et al., 2016).
Here, we focus on the effect of the atmospheric circulation on the tide. We used pressure indices (NAO and AO) that are relevant to represent atmospheric circulation. The NAO index represents the difference of normalized sea level pressure between the Azores high pressure system and the Iceland low pressure one (Hurrell, 1995). It indicates the redistribution of atmospheric masses between the subtropical Atlantic and the Arctic (Hurrell and Deser, 2009). In the North-East Atlantic, the similarity between the variations of the low-frequency winter NAO index and those of M 2 (Fig. 6) suggests a possible impact of large-scale atmospheric circulation on the tide.
The NAO index varies from positive to negative phases. Filtering the interannual variability, the NAO index tends overall to decrease between 1910 and 1970, then increase until 1990, and once again decrease. In the same way, M 2 amplitude tends to decrease up to 1960, then increase until 1990, and once again decrease. These similar patterns raise a possible connection between NAO and M 2 variation, already mentioned by Müller (2011) on the basis of qualitative criteria. In the following, we provide quantitative insights into the possible influence of NAO.
We computed the correlations (r value) between normalized M 2 and climate indices, NAO and AO (Fig. 7). M 2 , NAO and AO are filtered using the same time window (9 years). The correlations are computed since 1910, to have similar periods for all the stations. The correlations are considered as significant only if the p value is lower than 0.05 (95 % significance level). (Note that other statistics to measure the degree of association between the M 2 and NAO (AO) quantities would be worth exploring, for instance, nonlinear association using Spearman's correlation coefficient. In this respect, our study should be regarded as a first step that identifies sites worth considering in future investigations, especially investigating causal relationships with physics-based modelling.) The results are the following: (1) for NAO, 14 stations out of 18 show significant correlation. Note that at Brest, the correlation is significant since 1910, but not since 1864 (the NAO index used in this study starts only in 1864). This can be ex-  (4) For AO, we found similar, but overall larger, r values. This is not surprising as these two indices are closely related.
To go further in the relative contribution of MSL and NAO in M 2 variability, we fitted two linear regression models on M 2 variations. In the following, M 2 , MSL and NAO are filtered over 9-year time windows and normalized. At all the stations, we fitted M 2 variations with a MSL linear regression model (model 1) and a MSL and NAO multiple linear regres-sion model (model 2). Models 1 and 2 may be expressed as The correlations between M 2 and model 1 (MSL) and model 2 (NAO and MSL) are presented in Fig. 8. We checked if there was correlation between NAO and MSL at the stations: there is no correlation at six stations, and r value is between 0.2 and 0.6 at eight stations; see Fig. 8 and discussion below. The results are the following: (1) M 2 varies at first order with MSL (Fig. 8). (2) The introduction of the NAO (model 2) allows increasing the predictive performance of the model, beyond the inherent effect of adding an additional regression parameter. Indeed, on average, the Akaike information criterion (AIC) is 99.9 for model 2, instead of 112.7 for model 1. On average, the r 2 value is 0.67 for model 2 instead of 0.61 for model 1. At some stations, the increase  is quite large. For example at Cuxhaven, the r 2 value jumps from 0.42 to 0.64 between model 1 and 2. (3) The ratio β α+β represents roughly the relative contribution of the NAO compared to the total effect of MSL and NAO (Fig. 9), as MSL and NAO are normalized. We found a significant contribution at some stations (e.g. more than 30 % at Cuxhaven and Halifax), whereas it is negligible at others (e.g. only 5 % at Portland). A total of 8 stations out of 18 show large NAO contribution (> 20 %). The North-East Atlantic seems to be more sensitive to the NAO. Note that the interpretation of the results is tricky when MSL-NAO correlation is significant (orange bars in Fig. 8). For example, at Hoek van Holland, the relative NAO contribution is very small, mainly because MSL and NAO are highly correlated (r = 0.59). Figure 10 shows M 2 variations along with the predictions from the two models, at all four stations where the NAO contribution is significant ( β α+β > 0.25), and the correlation between M 2 and model 2 is large enough (r > 0.3). At Cuxhaven, Halifax and Key West, model 2 (MSL-and NAO-dependent) naturally captures the M 2 variations better than model 1 (MSLdependent); at Brest, the improvement is less significant. The trend switch observed since 1990 in the North-East Atlantic could be partly explained by the influence of the NAO on the tide.
These results suggest that a NAO-related mechanism may explain part of the variability of M 2 . As mentioned by Müller (2011), "it is shown that sea-level, sea surface temperature and Arctic ice thickness are correlated with the NAO index. Thus, changes in the dynamics of the atmosphere could affect both M 2 and S 2 tides by processes discussed under (1), (2) and (3)." An underlying mechanism linked with (2) -sea surface temperature -could be changes in the ocean stratification. This is one of main possible hypotheses invoked in Ray and Talke (2019) to explain secular changes in M 2 amplitude in the Gulf of Maine; this is also the main hypothesis in Müller et al. (2014) and Gräwe et al. (2014) to explain seasonal modulation of M 2 in the North Sea. The relationship between the NAO index and stratification is complex and spatially variable across the North-East Atlantic (Fromentin and Planque, 1996). In the North Sea, the sea surface temperatures are positively correlated with NAO (Becker and Pauly, 1996), while subsurface temperatures show no significant correlation with NAO (Tian et al., 2016). Stratification  could therefore be (positively) correlated with NAO, but a dedicated study, outside the scope of this paper, would be necessary. Another underlying mechanism linked this time with (1) -sea level -could be the difference of spatial distribution of water level, due to different sea-level pressure and wind stress patterns. This is the hypothesis invoked in Huess and Andersen (2011) to explain the seasonal modulation of M 2 in the North Sea. They ran a barotropic model, forced with tides only and with both tides and meteorological fields; their results show that the M 2 seasonal modulation is better captured when the model is forced with both tides and meteorological fields rather than with tides only. Figure 11a shows the average sea-level pressure during the period 1850-2015, derived from the Twentieth Century Reanalysis (20CR) (Compo et al., 2011;Slivinski et al., 2019). A positive NAO winter (e.g. 1989) corresponds to a situation with a stronger pressure gradient than average, between the two pressure systems of Azores and Iceland (Fig. 11c). By contrast, a negative NAO winter (e.g. 1969) corresponds to a weaker gradient pressure than usual (Fig. 11b). (We define winter here as December-February.) This way, from one year to another, the large-scale atmospheric masses are dis-L. Pineau-Guillou et al.: Large-scale changes of the semidiurnal tide tributed differently, and as a consequence, the water volumes are also distributed differently in the North Atlantic. In a situation of NAO + , the surface waters are pushed onshore by westerly winds, moving from Iceland to the European coasts of France, Spain and Portugal. Figure 12 shows the redistribution of the sea-level pressure, between two years with high and low NAO indices (here 1989 and 1969). Note that this is an extreme situation, as these years have strong positive and negative indices. Assuming an inverse barometer response of sea level, the changes in terms of water level may vary from more than 24 cm in the northwestern part of the area to around −12 cm in the region that includes most of the North-East Atlantic tide gauges considered in this study. This variation of a few tens of centimetres is probably negligible offshore but may have some impact on tide propagation along the continental shelves and in shallow waters. It could also shift slightly the amphidromic points. Assuming that these changes have a similar impact (in terms of magnitude) on M 2 as MSL changes, that is, ±10 % in shallow waters according to recent simulations (Pickering et al., 2017;Idier et al., 2017), we find that they can yield centimetric changes in M 2 amplitude. In other words, their order of magnitude is roughly in agreement with the changes observed in M 2 (Table 1). However, it is difficult to disentangle the effects of stratification and meteorological forcing (sea-level pressure and wind stress) in M 2 changes, and possibly both mechanisms coexist. Dedicated simulations should be conducted to assess the effects of atmospheric forcing on M 2 variability.

Conclusions
We investigated the long-term changes of the principal tidal component M 2 over the North Atlantic coasts. We analysed 18 tide gauges with time series starting no later than 1940. The longest is Brest with 165 years of data. We carefully processed the data, particularly to remove the 18.6-year nodal modulation.
We found that M 2 variations were consistent at all the stations in the North-East Atlantic (Cuxhaven, Delfzijl, Hoek van Holland, Newlyn, Brest), whereas variations appear between stations in the North-West Atlantic. The changes started long before the 20th century and are not linear. The trends vary significantly from one station to another; they are overall positive, up to 2.5 mm/yr, or slightly negative. Since 1990, in many stations, the trends switch from positive to negative values. The significant differences between the trends since 1910 and 1990 indicate caution when interpreting trends based on short records, i.e. less than 30 years, especially if the data are noisy (interannual variability) and the underlying processes non-linear (change points).
Concerning the causes of the observed changes, M 2 varies primarily with the MSL, but MSL rise is not sufficient to explain the variations alone. The similarity between the North Atlantic Oscillation and M 2 variations in the North-East At-lantic suggests a possible influence of the large-scale atmospheric circulation on the tide. Our statistical analysis confirms large correlations at all the stations in the North-East Atlantic. The trend switch observed since 1990 could be the signature of the large-scale atmospheric circulation on the M 2 tide. The underlying mechanism would be a different spatial distribution of water level from one year to another, depending on the low-frequency sea-level pressure patterns, and impacting the propagation of the tide in the North Atlantic basin. In the future, dedicated modelling studies should be undertaken to confirm or discard this hypothesis. These simulations should also allow estimating the effect of the wind (through the Ekman current) and currents on M 2 changes (Devlin et al., 2018).
In this study, we focused only on M 2 amplitude. A similar analysis on the phase lag would draw a more complete picture of the M 2 variations (Müller, 2011;Woodworth, 2010;Ray and Talke, 2019). Other constituents are also affected. Results show that S 2 amplitude decreases at all the stations located in the North-West Atlantic and, in contrast, tends to increase in the North-East Atlantic (not shown). The largescale decrease of S 2 observed in the North-West Atlantic is consistent with previous studies (e.g. Ray, 2006, in the Gulf of Maine). Further investigations should be definitely conducted to extend this work to more constituents.
The historic data show that the changes started long before the 20th century. This conclusion would not have been possible without the huge work of data rescue undertaken over the past decades (e.g. Pouvreau, 2008;Bradshaw et al., 2016). This underlines the great importance of sea level data archaeology, which allows extending and improving historical datasets (Pouvreau, 2008;Woodworth et al., 2010;Marcos et al., 2011;Jay, 2013, 2017;Ray and Talke, 2019;Bradshaw et al., 2015Bradshaw et al., , 2020Haigh et al., 2019). This is essential for studies related to climate change.
Finally, we should mention several additional limitations and perspectives in this study. (1) We processed the time series considering that they were quality controlled. A fuller analysis of the data quality before processing would probably be valuable. (2) We did not investigate the history of each station. There are probably some local changes (e.g. environment or instrumentation) that may explain a part of the variability of M 2 amplitude and some discrepancies between stations. (3) The tide gauges are located mainly in harbours. They are affected at the same time by local-and regional/global-scale changes, which are difficult to separate. Moreover, they may not be representative of changes offshore. A similar study based on satellite altimetry data would probably be of great interest, even if temporal scale for satellite data is still rather short (i.e. < 30 years) compared to climate-scale processes. (4) We focused mainly on the UHSLC dataset, which consists of 249 stations in the Atlantic Ocean. Other relevant stations (which are not in this dataset) may be considered in future studies. (5) We did not  investigate the impact of storminess on the tide. Dedicated studies are necessary to estimate if changes in storminess could affect significantly tidal constituents. (6) We used only winter AO and NAO indices, which show more variability than annual indices. A similar analysis with annual indices shows similar results for the correlation with AO or NAO (positive correlation on the North-East Atlantic). With annual rather than monthly indices, the difference of pressure fields will decrease, and as a consequence, the magnitude of the sea-level response will also decrease. Further investigations should be conducted on this point. The M 2 component is subject to a 18.6-year modulation, separated from a neighbouring line in the tidal potential (m 2 ) whose Doodson number differs in its fifth frequency (255 555 and 255 545 for M 2 and m 2 , respectively) (Doodson and Warburg, 1941;Pugh and Woodworth, 2014). This fifth frequency corresponds to N , the negative of the mean longitude of the Moon ascending node -hence the "nodal" term -whose period is 18.6 years. Note that there is also another component close to M 2 , whose Doodson number differs only from the fifth frequency (255 565), but it is negligible, its amplitude in the tidal potential being only 0.05 % of M 2 , whereas m 2 amplitude is 3.7 % of M 2 (Simon, 2007(Simon, , 2013. With one year of hourly data, the two components M 2 and m 2 cannot be separated by a yearly harmonic analysis (at least 18.6 years are necessary). As a consequence, M 2 amplitude is modulated by m 2 . However, we can estimate this modulation and remove it. The harmonic formulation is expressed schematically as a sum of harmonic components: where h(t) is the sea level height at time t, V i (t) is the astronomical argument (computed from Doodson number) and a i , κ i the amplitude and phase lag of each component. Considering that M 2 and m 2 are very close in terms of frequency, we can assume that their phase lags are similar (κ M 2 κ m 2 ). As their difference of astronomical arguments is V m 2 − V M 2 = N + π , the M 2 and m 2 contributions to the total water level may be expressed as h M 2 (t) + h m 2 (t) = h M 2 (t)[1 + f nod cos(N + π )], where f nod , the nodal modulation, is the ratio of the amplitude of m 2 and M 2 . As M 2 and m 2 are very close in terms of frequency, f nod is generally considered as close to the ratio of their amplitude in the tidal potential, A m 2 and A M 2 : f nod = a m 2 a M 2 A m 2 A M 2 0.037. (A3) Figure A1. (a) Estimation of the nodal modulation of M 2 amplitude (mean removed) at Newlyn. (b) Impact on M 2 amplitude of the nodal modulation correction at Newlyn. M 2 is detrended in (a) to better fit the nodal modulation.
The negative of the mean longitude of the Moon ascending node is expressed simply as a function of time (p. 116 in Simon, 2007, p. 112 in Simon, 2013: with N in degrees, and T the time elapsed since 1 January 2000 at 12:00, expressed in Julian centuries (36 525 d).
The tidal program we used (MAS) corrected M 2 applying the usual 3.7 % nodal modulation (Eq. A3). However, this value may vary significantly from one station to another; Ray (2006) reported values ranging from 2.3 % to 3.6 % in the Gulf of Maine. Here, we computed directly f nod from the observed data, proceeding as follows.
(2) We detrended the obtained signal removing the last intrinsic mode function (IMF) of an empirical mode decomposition (EMD) (Huang et al., 1998); note that the EMD is an analysis tool which partitions a series into "modes" (i.e. IMFs), the last one being the trend of the signal. (3) We fitted a function a m 2 cos(N + π ) to this detrended signal to estimate a m 2 , N being expressed as in Eq. (A4). (4) We finally computed f nod as the ratio between m 2 and M 2 amplitudes (Eq. A3). Figure A1a shows an example of estimate of M 2 modulation at Newlyn: the fit leads to a nodal modulation of 3.3 %. Note that this value is consistent with Woodworth (2010) (3.2 %), whereas Woodworth et al. (1991) gave a slightly different value (2.8 %). Figure A1b shows the impact of this value rather than the default one: oscillations of 18.6 years are clearly reduced. Note that in this study, the m 2 amplitude -and then the nodal correction -could have been computed from the full time series harmonic analysis, as records are longer than 18.6 years. However, the method presented here to compute the nodal correction can be applied even for time series shorter than 18.6 years.
The computed nodal modulations are summarized in Table 1 (column 6). They vary from 0.8 % to 4.1 %. Note that these values are consistent with those obtained by previous authors (Ray, 2006;Müller, 2011;Woodworth, 2010;Ray and Talke, 2019). Only the value at Charleston differs significantly: 3.0 % in our study compared to 3.7 % in Müller (2011).

31
Appendix B: Time series of annual M 2 amplitude at all the stations Figure B1. Annual M 2 amplitude at the 18 selected tide gauges.