Assessing 20th century tidal range changes in the North Sea

In many places around the world, tide gauges have been measuring substantial non-astronomical changes. Here we document an exceptional large spatial scale case of changes in tidal range in the North Sea, featuring pronounced trends between -2.3 mm/yr in the UK and up to 7 mm/yr in the German Bight between 1958 and 2014. These changes are spatially heterogeneous, suggesting a superposition of local and large-scale processes at work within the basin. We use principal component analysis to separate large-scale signals appearing coherently over multiple stations from rather localized changes. We identify two leading principal components (PCs) that explain about 69% of tidal range changes in the entire North Sea including the divergent trend pattern along UK and German coastlines, which suggest movement of the region’s semidiurnal amphidromic areas. By applying numerical and statistical analyses, we can assign a baroclinic (PC1) and a barotropic large-scale signal (PC2), explaining a large part of the overall variance. A comparison between PC2 and tide gauge records along the European Atlantic coast, Iceland and Canada shows signiﬁcant correlations on time scales of less than 2 years, which suggests an external and basin-wide forcing mechanism. By contrast, PC1 dominates in the southern North Sea and originates, at least in part, from stratiﬁcation changes in nearby shallow waters. In particular, from an analysis of observed density proﬁles, we suggest that an increased strength and duration of the summer pycnocline has stabilized the water column against turbulent dissipation and allowed for higher tidal elevations at the coast.


Introduction
For thousands of years, tides have had a great influence on coastal areas globally and their residents.Today they play a critical role in influencing economic considerations, nautical safety, renewable energy schemes, assessments of land erosion, and the definition of geodetic datums (Haigh et al., 2020;Pugh & Woodworth, 2014).Tides not only control the navigability of some ports and sea routes, but also have a major influence on the intensity and timing of extreme sea levels during storm surges (e.g., Arns et al., 2020;Horsburgh & Wilson, 2007;Prandle & Wolf, 1978).Given their close connection to the periodic and predictable nature of astronomical variations, the amplitudes and phases of tidal constituents, and corresponding tidal water levels, are generally assumed to be constant on time scales over which basin geometry undergoes only minor changes (i.e., decades to centuries).However, Keller (1901) showed increased tidal amplitudes due to reflection and local resonance changes as a result of building measures such as weirs (e.g., in the Ems River).Similarly Doodson (1924) pointed to appreciable secular perturbations in the local tidal regimes of particular ports, weirs, and estuaries.More recently, the topic of changes in ocean tides has been revived and extended to the scales of shelves, basins and the global oceana development fueled by the digitization and publication of global data sets of tide gauge records, see Woodworth et al. (2017).In fact, statistically significant trends of tidal parameters of the order of a few percent (in relative terms) are now well documented around the world (e.g., Flick et al., 2003;Jay, 2009;Mawdsley et al., 2015;Ray, 2009;Talke & Jay, 2017;Woodworth et al., 1991;).Fluctuations of similar magnitude and regional extent have been observed on interannual time scales (e.g., Devlin et al., 2014;Feng et al., 2015;Müller, 2011;Ray & Talke, 2019).
Despite this ample evidence of changes in tides in water level series, the forcing factors and spatial extent of secular and short-term variability in tides remain uncertain.Woodworth (2010) succeeded in detecting coherent patterns of amplitude and phase trends in primary constituents along the North American coasts, but found less regional consistency in data from Asia, the Australian Seas or Europe.However, some spatially coherent changes could still be observed in smaller and well-instrumented areas.A major problem identified by Woodworth (2010) is that small-scale (often site-specific) and large-scale changes may occur simultaneously, thereby impeding research of the underlying physical processes.Over wider coastal sections, and at sites open to the sea, the effects of a rise in mean sea level (MSL) on tidal wave propagation explain only a fraction of the observed trends (Müller et al., 2011;Schindelegger et al., 2018).
Accordingly, the assumption persists that other mechanismssuch as changes in stratification, turbulent dissipation, and variations in shoreline position or bed roughnessplay major roles; see Haigh et al. (2020) for a review.The present consensus is that in many areas of the world a combination of different oceanographic processes may be at work.For instance, Ray & Talke (2019) suggest that the large secular changes of the lunar M 2 tide in the Gulf of Maine could be caused by both sea level rise and persistent stratification changes.Yet, as implied above, any contributing mechanism will act on its own characteristic spatial and temporal scales, overlaying and possibly reinforcing other processes.This particularly applies to anthropogenic construction measures (e.g., building of dykes and tidal barriers) that can cause transient perturbations to the local tidal regime and affect adjacent stretches of coastline (Talke & Jay, 2020).Therefore, a major challenge is the separation of local effects and large-scale changes and their subsequent attribution to certain forcing factors.
Exceptional changes in tidal range are found in the German Bight, documented in Führböter & Jensen (1985) and Jensen (2020).Between 1958 and 2014, different changes in tidal range were detected ranging from approximately 3% at some of the investigated tide gauges to more than 11% at others (Figure 1).The latter is equivalent to a trend of 5.7 mm/yr and outpaces the simultaneous local (Dangendorf et al., 2015) and global MSL rise (Dangendorf et al., 2019;Oppenheimer et al., 2019).To our knowledge, this magnitude of tidal range change is one of the highest in the world, only exceeded by developments in the Gulf of Maine (Ray & Talke, 2019).
It further seems that the overlap between local and large-scale effects in the North Sea is particularly pronounced, possibly nurtured by the region's character as a shelf sea with a tide generated in the Atlantic.Previous research (summarized in Jensen et al., 2014) has ruled out astronomical, large-scale morphological or tectonic causes (at least in the German Bight), but pointed to the generally non-linear and non-uniform behavior of water levels in the North Sea.
To improve our understanding of these puzzling tidal range changes, we aim to address the following questions through systematic data analysis: (1) Are these changes on different time scales detected within the German Bight a localized phenomenon, or are they part of a largerscale development spreading over larger areas within or even outside the North Sea region?(2) Is it possible to separate and quantify large-scale and small-scale effects from observed records?
(3) If (2) is the case; can we attribute physical causes to the observed changes?Below, we first discuss geographic and oceanographic characteristics that are fundamental to the understanding of the tidal regime in the North Sea, the available database, its limitations and major processing steps (Section 2).Section 3 introduces the analytical methods of Ordinary Kriging, which is here mainly used for gap-filling as the following PCA requires complete time series.The results of our analyses are described extensively in Section 4. To answer the abovementioned research questions, we start our analyses with the detection of observed changes in the tidal range at individual sites.In a second step, we apply a PCA to identify modes of variability common to all (or the majority of) sites and to distinguish them from local anomalies.
In a last step we analyze potential causes and drivers of the observed changes.The paper concludes with a summary and additional remarks in Section 5.

Study area
The North Sea is one of the largest shelf seas on Earth with a size of about 575.300 km² (Huthnance, 1991).Counted counter-clockwise, its margins comprise coastal sections of the United Kingdom, France, Belgium, the Netherlands, Germany, Denmark and the south of Norway (Figure 2).The North Sea is connected to the North Atlantic via a large inlet between Scotland and Norway in the north and a narrow opening through the English Channel in the southwest and it opens to the Baltic Sea in the east.Water depths in the North Sea are on average 90 m but vary greatly, generally increasing from south to north.While the southern parts are often shallower than 40 m with lowest depths in the German Bight, they increase to about 300 m at the continental shelf toward the Norwegian Trench and toward the entry into the Norwegian Sea in the northwest.There are also extensive shallow water regions off the south-eastern coast of the UK known as the Dogger Bank complex, with their western part extending to the coasts of Norfolk and Suffolk (Quante & Colijn, 2016).
The tidal regime in most parts of the North Sea is strongly influenced by the astronomical tides entering the basin from the Atlantic and therefore mainly semidiurnal.The greater part of these oscillations enters between the Shetlands and Scottish mainland and a smaller part through the English Channel.They travel counter-clockwise through the entire North Sea basin as Kelvin waves.The entry times of the tidal high and low waters are therefore shifted relative to each other according to the celerity of the tidal wave.This physical setting results in three amphidromic points, one close to the English Channel, one off the coast of Norway and one central in the North Sea basin (Proudman & Doodson, 1924).Since the North Sea's basin shape is close to the resonance frequency in the semidiurnal spectral band, the superposition of the principal lunar and solar tides M 2 and S 2 leads to a significant spring neap cycle.These two constituents cause a potential tidal range between 1 and 5 m (Quante & Colijn, 2016).
Accordingly, the tidal regime of the North Sea can be classified as macrotidal (>4 m), mesotidal (2-4 m) and microtidal (<2 m) (Haigh, 2017), with the actual tidal range being strongly influenced by local factors.For example, the mean spring tidal range at the east coast of the UK varies between 3.60 m (Aberdeen) and 6.20 m (Immingham) (Horsburgh & Wilson, 2007).The mean tidal range in the data set used below is about 3.40 m in the UK and the English Channel, 1.98 m at the Dutch west coast, 2.33 m at the Dutch north coast and 2.82 in the German Bight.
Figure 2: Bathymetry of the North Sea (Becker et al., 2009;Schrottke & Heyer, 2013).Also shown are the locations of tide gauges (black dots) used in this study including their respective numbering (see also Table 1).The black propellers indicate the location of the three semidiurnal amphidromic areas (including the amphidromic points for the M 2 and S 2 constituent) and the black dotted lines indicate contours of equal mean tidal range (Sündermann & Pohlmann, 2011).

Data
Time series of water level from 70 available tide gauges around the North Sea basin were collected from various sources.Data from GESLA (Global Extreme Sea Level Analysis, GESLA, Woodworth et al. 2017), Open Earth (Deltares) and the responsible German authorities (Wasser-und Schifffahrtsverwaltung des Bundes via the portals of the associated Central Data Management, ZDM) were used.The available time series vary considerably in length and completeness.The earliest measurements in the form of tidal high and low water readings are from 1843 (the tide gauge at Cuxhaven Steubenhöft, Germany), while on the Dutch coast data from some stations have only been digitally available since the 1980s.High-resolution data sets with an equidistant sampling between 1 and 60 minutes were used as well as time series of tidal high and low water.We excluded equidistant time series with a resolution lower than 60 minutes, as supplemental analyses have shown that they insufficiently describe the height and timing of individual tidal high and low waters.The tidal range was calculated as the difference between each tidal high water and the mean of the two surrounding tidal low waters, according to the German standard (DIN 4049-3, 1994).From those, we calculated monthly averages and removed the mean seasonal cycle, as we are mainly interested in longer-term changes.
Considering the 18.6-year nodal cycle and the end of numerous water level series in December 2014, we adopt an analysis period from January 1958 to December 2014; approximately 3 nodal cycles.Tide gauges known to be located near to weir installations or in rivers were excluded, as these are at least partially separated from the oscillation system of the North Sea.Seventy time series of tidal range remained in the data set, forming the basis for our investigations (Table 1, Figure 2 and Figure 3).Acknowledging the counter-clockwise propagation direction of the tidal wave, the tide gauges used in this study are counted by starting at Lerwick (Shetland Islands) and ending at Tregde (Norway).The average completeness of the stations is 64% in the UK, 65% in the Netherlands and around 88% in Germany.In addition to the procedures explained in the following sections, linear trend analysis, harmonic 199 analysis of tidal constituents and wavelet coherence analysis were carried out to characterize 200 multiple feature of the tide gauge records in the North Sea.Any significance statements made 201 throughout the manuscript are based on a 95% confidence level.We calculated linear trends 202 using ordinary least squares regression and assessed their significance by considering normally distributed but serially correlated residuals following an autoregressive process of the order 1 (e.g., Mawdsley & Haigh, 2016).Annual amplitudes for the leading constituents were determined by a harmonic analysis using the MATLAB toolbox U-Tide (Codiga, 2011) and the wavelet analyses were conducted with the MATLAB package of Grinsted et al. (2004).None of these methods are explained here in detail due to their general recognition and widespread use.Research in Environmental Sciences (CIRES) to describe the meteorologically induced effects on water levels (Compo et al., 2011).

Kriging
Kriging (also Gaussian process regression) is a geostatistical method to interpolate missing values based on information stemming from neighboring stations (i.e.their covariance matrix).It is here mainly used for gap-filling as the following Principal Component Analysis (PCA) requires complete time series.Originally developed in the 1950s for mining purposes (Krige, 1951), this method has been used increasingly in other areas including the analysis and interpretation of incomplete surface air temperature fields (Rigor et al., 2000;Rohde et al., 2013).In general, Kriging is a linear interpolation procedure.Missing values are determined according to a given covariance matrix, which is calculated from the existing observations (Cressie, 1990).Kriging provides some important advantages over other interpolation procedures.The interpolated values change smoothly and always pass through the observed values at the sample points.Problems related to the accretion of measurement points are avoided by considering the statistical distances between the neighbors used in the interpolation of a certain value, which means that the spatial variance is taken into account.If clustering occurs in a region, the weights of the affected sample points are reduced by including the density.In sparse regions, only the distance is considered.The procedure can be summarized with the formula where Z ̂ is the query value at the unobserved location x 0 and i = 1 … n represents a running index over n observations.Z ̂ is computed from a linear combination of all observed values z i = Z(x i ), which are weighted by the parameter w according to distance and density.A special property of the Kriging procedure is the convergence of interpolated values to the mean value of their region with increasing distance to the available samples.That is why Kriging estimates at query points tend to be conservative (Cowtan & Way, 2014).In keeping with this characteristic, the general tidal range behavior worked out later in Section 4.1 is also valid when the Kriging step is omitted.
We use Kriging for two different purposes.First, the temporal gaps in the tidal range data (Section 2.2) were closed for each monthly time step in the investigation period.Figure 3-a illustrates that this is a relevant issue in the Netherlands, in particular before 1970, while in the UK data gaps occur before 1990.Second, additional data points along the coastline of the North Sea were interpolated, allowing us not only an analysis of the temporal evolution of each station series in terms of a linear trend but also a spatial analysis of the different developments (Figure 4).For both applications, we use the Ordinary Kriging algorithm of Schwanghart (2020).Note also that in transitioning from

Principal Component Analysis
Principal Component Analysis (PCA), a method of multivariate statistics, is used to structure and simplify extensive data sets by approximating a large number of statistical variables with a smaller number of significant, non-correlated (orthogonal) linear combinations.If x is a vector with n random variables, first a linear function f 1 (x)dependent on constant coefficients c 1iis determined by calculating the eigenvector from the spatially weighted covariance matrix of x.
Then f 1 (x) represents the largest possible overall variance of all variables in x: This decomposition process is repeated for a function f 2(x) , which is uncorrelated with f 1(x) and describes the largest possible amount of the remaining variance.It is possible to find n such functions, but the purpose is usually to explain as much variance as possible with significantly fewer functions f i(x) , known as Principal Components (PCs) (Jolliffe, 2002).Therefore, the PC of a temporally and/or spatially varying physical process represents orthogonal spatial patterns, in which the data variance is concentrated.Using the leading PC, an approximate reconstruction of the observed variable can be generated.This type of analysis is often used in Earth system sciences to identify spatial and temporal patterns of climate oscillations.
In this study, we apply PCA to the entire monthly de-seasoned tidal range data set from the 70 sites (Figure 2), whose gaps were previously filled through Ordinary Kriging.If there are indeed large-scale signals affecting the tidal range in the North Sea, they should appear as a coherent pattern visible at multiple sites, and therefore be visible in the leading PCs.By contrast, spatially confined ("small-scale") anomalies in tidal range will be shifted into the higher PCs, as these can only be responsible for a small part of the overall variance.Such shifting includes not only the response of the local tidal system to, for instance, anthropogenic construction measures but also to changes in bathymetry or morphology.Local effects can explain more variance than largescale effects at individual sites or small subsets, but never for the entire data set.It is therefore important to consider the explained variance of the PCs at each tide gauge individually to ensure that large-scale effects with a very small influence on the overall variance are retained.With this approach, the PCA enables us not only to attribute tidal range changes to small-scale and largescale effects, but also to calculate the spatial extent and the temporal development of patterns that might reflect important environmental factors.

Trends of tidal range and tidal constituents
To address the three research questions defined in the introduction, we first map the spatial extent of the long-term changes in tidal range in the study area.We start our analysis by calculating linear trends for each individual record over a common period between 1958 and 2014 and map them in Figure 4.In this step of the analysis, the time series of Lerwick (Shetland Islands) and Tregde (Norway) were omitted, since both are the only available tide gauges within large areas and, therefore, there is not a sufficient data density for use by the Kriging algorithm.
We identify a variety of trends with a particularly pronounced spread in the southern parts of the basin.While there are no significant trends at the north-eastern coast of the UK, negative trends occur further south between Immingham and Dover.Here, six of eight stations show significant negative trends while the remaining two do not differ significantly from zero.In this area, Immingham shows the largest negative and statistically significant trend (-2.3 ± 0.5 mm/yr) of all sites, while the smallest negative trend of -0.7 ± 0.3 mm/yr is found in Felixstowe.The mean value for all tide gauges in this area is -1.0 mm/yr.In contrast, trends turn positive on the continental side of the English Channel and the European West Coast.Our assessment reveal increasing trends following the coastlines of France (trend at tide gauge Dunkerque 1.3 ± 0.4 mm/yr), Belgium and the western Netherlands up to the tide gauge at Huibertgat (0.8 ± 0.2 mm/yr), near to the German-Dutch border.On average, trends along the European West Coast are 0.8 mm/yr.Hereafter, sharp trend increases are found within a short distance, reaching values of more than 7 mm/yr in the German Bight area.Here, the average trend in tidal range amounts to 3.3 mm/yr (Table 2).Local changes affect some tide gauges like Den Oeverbuiten (Netherlands) or Büsum (Germany), which at first sight seem to contradict this spatial pattern.
We  The identified dipole-like trend pattern has its node approximately at the longitude of the English Channel (Figure 4) and suggests a westward displacement of the main low amplitude areas (including amphidromic points of M 2 and S 2 ) located in the central North Sea and near the English Channel (Figure 2).To obtain further indications of such a shift, we perform a harmonic analysis to determine the main semi-diurnal M 2 and S 2 tidal constituents, which make the largest contributions to the tides in the North Sea.Since high-resolution hourly time series with a coverage of at least 75% between 1958 and 2014 are required for a tidal analysis, only a subset of 28 tide gauge records is appropriate for our assessment.The available database is thus reduced and fewer stations show significant trends (20 for M 2 , 14 for S 2 ).Nevertheless, the overall finding (Figure 5) are similar to the assessment focusing on tidal ranges highlighted in Figure 4; that is for both constituents (though with larger magnitude for M 2 ), negative trends occur in the southeast of the UK and the highest positive trends are found in the German Bight area.A displacement of the M 2 and S 2 amphidromic point is, therefore, also inferred.

German Bight
).This argument is supported by several numerical modelling efforts (Idier et al., 2017;Pickering et al., 2012;Schindelegger et al., 2018), in which the impact of large (0-10 m) MSL increases on leading constituents (mainly M 2 ) were investigated.Complementary to our empirical assessment, they all detected (at least qualitatively) similar patterns as shown in Figure 4 and 5.However, closer examination also reveals some discrepancies and the model results do not correspond exactly to the measured data.For instance, both Pickering et al. (2012) and Schindelegger et al. (2018) predict an increase in M 2 amplitude in the southwestern part of the North Sea, between Suffolk/Essex and the Netherlands, while we detect negative trends in Suffolk/Essex and positive trends in the Netherlands.This disparity could be caused by two facts: (1) the assumed MSL rise projections used by these studies have not yet occurred and are therefore theoretical in character, and (2) a contribution of effects not yet considered by numerical models.

Principal Components and large-scale effects
Our results of the linear trend analysis point towards a distinct spatial pattern that is occasionally interrupted by diverging trends at individual locations.To further distinguish between the largeand small-scale effects of tidal range changescomprising both trends and short-term variability we apply PCA (Figure 6).The first two PCs, which are presented in Figure 6, explain about 69% of the total variance in the entire data set (PC1: 55%, PC2: 14%), while each of the remaining 68 PCs contributes between 0.01 and 4%.Additionally, no other PC represents significant parts of the variance at a larger number of tide gauges and is therefore rather local in character.This indeed suggests that the two leading PCs reflect coherent large-scale effects, represents the decrease in tidal range at the south-eastern coast of the UK.This contrast is also reflected in the correlation coefficients of the first two PCs with the measured tidal range changes (a metric that is mostly influenced by inter-and intra-annual variability).Figure 6-c show moderate but significant correlations of 0.3 -0.5 for PC1 at the south-western boundary of the North Sea and displays the highest values (~ 0.9) in the area of the German Bight.A contrasting picture emerges for PC2.In the area of the German Bight, correlations with tidal range changes are non-significant and close to zero but almost consistently above 0.7 and significant in the UK (Figure 6-d).
These patterns are also confirmed when considering the explained variance for particular clusters of tide gauges.Along southeastern UK coastlines, where negative trends are found, the explained variance of PC1 amounts to only 3%, while PC2 explains about 58% (Table 2).In the Netherlands, the mean explained variance for PC1 is 45% and only 10% for PC2.The contribution of the second mode drops to 3% in the German Bight, whereas PC1 explains 77% of the variance on average.This spatially reversing pattern is also detectable in the coefficients for PC1 and PC2 (Figure 6-b), just as in the linear trends of the tidal range observations.Apparently, PC1 with its positive slope is more pronounced in the area of the German Bight, whereas PC2 (negative slope) dominates in the southeast of the UK.This indicates different underlying physical mechanisms for these large-scale signals.

Impacts on local tidal range
After identifying two large-scale patterns relevant at the majority of tide gauge records in the North Sea, we next ask whether we also can identify small-scale effects using the residual signal after removing the linearly regressed PC1 and PC2 at individual sites.Figure 7-d shows the explained variances and indicates that alongside the described contrast between PC1 and PC2, local influences play a major role in some cases.Especially noticeable again are tide gauges Den Overbuiten (Netherlands, #33) and Büsum (Germany, #60) due to their high percentage of local effects.For example, PC3 (explained overall variance: 4%) captures more than 50% of the variance at tide gauge Büsum and around 30% at Cuxhaven (Germany, #59).This anomaly is reflected in the comparison of the measured trends with those from re-synthesizing PC1 and PC2 the linear trends from the reconstruction (significant trends outlined) and (d) explained variance of the two PCs as share of the total variance

Identifying physical causes
The PCA suggests two modes of variability that appear coherently at the investigated sites in the North Sea.Now the question naturally arises whether these signals are produced within or outside the basin.If the former is the case, then the corresponding PCs should show no correlations to tide gauge records from the adjacent North Atlantic, while an external forcing would possibly provide some sort of coherence with those records.No coherence is found for PC1 and we therefore conclude that it is produced within the basin, which will be explained later.
The opposite applies to PC2.A comparison between PC2 and available tide gauge records along the European Atlantic coast, Iceland and Canada is shown in Figure 8. Figure 8-c indeed documents high and significant correlations of about 0.7 on average between PC2 (calculated exclusively on the basis of North Sea data set) and Atlantic tide gauge records spanning the region from the English Channel southward to Spain.Moreover, there are significant correlations of 0.64 in the north (Reykjavik, Iceland), and even in the Northwest Atlantic (still reaching 0.46 in Port-aux-Basques, Newfoundland) (Figure 8-a/c).Further south towards the Gulf of Maine, these correlations disappear (not shown).A supplemental wavelet analysis (not shown) further reveals that the common oscillations between PC2 and the measured tidal range changes mainly occur on time scales from 6 to 24 months with particularly high coherence at around 12 months.
We interpret this finding as an indication for a common high-frequency signal in the North Atlantic, causing widespread changes in tidal range.
In order to narrow down the possible causes, outputs from an updated barotropic shallow-water model run by Arns et al. (2015a,b)  The effects of bottom friction are more involved, but some simple geometric considerations are instructive.As the tidal wave enters the extensive shallow water areas of the southern North Sea, energy losses due to friction become dominant, yet the influence of PC2 is increasingly attenuated in the direction of propagation (Figure 6-d).This discrepancy suggests that frictional effects do not represent the physical cause of PC2, although they might play a role in suppressing the magnitude of PC2 in the highly dissipative eastern North Sea region.As our simulations were performed with an invariant bathymetry and no changes to friction parameters, sea level rise and meteorological forcing remain as possible causes.We accordingly analyzed correlations between PC2 and these factors (MSL rise, atmospheric pressure loading, wind velocities and directions) but could not detect a clear and significant linear relationship.In this context, Arns et al. (2015a) already referred to the numerous non-linear relationships between the individual parameters in marginal seas.Specifically, the nonlinear interaction between tide and sea level rise as well as the dynamic response of the sea surface to meteorological forcing are important (see also Arns et al., 2020).Further analyses, in particular sensitivity studies taking into account altered tidal boundary conditions and time variable friction coefficients, will perhaps allow for a final identification of the ultimate driving factors (e.g., Rasquin et al., 2020).While the signal of PC2 is reproducible, PC1 cannot be detected in the simulated data, which means PC1 is absent in the barotropic model.At the beginning of this section we stated that there is no coherence to the Atlantic tide gauges for PC1, which suggests a origin of the signal within the basin.We thus conjecture that a baroclinic, density-related effect inside the North Sea is responsible for PC1 and attempt an explanation in terms of known relationships between tidal currents and turbulent energy losses in varying stratification conditions.This attribution primarily arises from considerations at seasonal time scales.Using hydrographic casts and baroclinic model simulations, Müller et al. (2014) linked M 2 elevation changes of 1-5 cm in the southern North Sea to the see-sawing of continental shelf stratification between statically stable summer and well-mixed winter conditions.Strong buoyancy gradients in mid-depths (20-30 m) of shallow waters arise during summer months (see e.g., van Haren et al., 1999) and stabilize the water column against energy losses to vertical mixing.The associated increase in barotropic tidal transport and surface elevations was found to be most pronounced in very shallow areas and for cyclonic rotation of strong tidal currents (Müller, 2012) conditions that are all present in the North Sea.
To relate at least parts of the PC1 content to this process, we analysed the temporal evolution of the North Sea's density structure based on gridded temperature and salinity profiles from the KLIWAS dataset (Bersch et al., 2016).These data are provided as annual values through to 2013 at comparatively high spatial resolution (0.25° x 0.5° latitude-longitude boxes, 2-5 m depth intervals).For consistency, the monthly PC1 series was binned to annual values  with respect to the length of the KLIWAS dataset) and cleaned from secular changes with periods longer than 30 years.Because it is unknown how well KLIWAS represents the smaller, more subtle secular trends of density across the water column, we limit our comparison between stratification and PC1 to changes on interannual time scales.To suppress noise in the climatology, vertical density profiles from a particular set of grid points around the German Bight were averaged to a mean water column structure per year (Figure 9).These query points, indicated by black dots in Figure 9-b, lie within 2° of 54.5°N/6.0°Eand have an exact depth of 35 m in the KLIWAS dataset.The sampled area is shallow, hosts strong tidal currents, and is not permanently mixed, thus favoring a potential effect of stratification on tides.The corresponding time-averaged density profile (Figure 9-a) indicates a pycnocline at 20-25 m, conforming in principle to modeling results (e.g., Guihou et al., 2018;van Leeuwen et al., 2015).While this agreement is reassuring, we also note that our crude spatial averaging ingests profiles in various states of stratification (i.e.homogeneous, seasonally or intermittently stratified conditions, see Leeuwen et al., 2015).Given the tendency for in situ measurements being taken in summer, the KLIWAS dataset may, however, mainly represent the seasonally stratified case.gradients amount to about 0.3 kg/m³ per 10 m of depth and thus correspond to the order of magnitude that maintains the seasonal cycle of M 2 in this region (Müller et al., 2014).Therefore, all indications are that changes to the intensity of summer stratification and/or the time spent in a stratified (or mixed) regime over the course of a year cause the variance in tidal range represented by PC1.A breakdown into different modes of stratification variability is tempting but beyond the scope of our study as it would call for consideration of several factors, including freshwater buoyancy input, variable local wind stirring, and the inflow of Atlantic water masses through the northern and southern boundaries (Mathis et al., 2015).Nevertheless, we have analysed long-term hydrographic data of the North Atlantic and detected high negative correlations (-0.8) between PC1 and temperature of the upper ocean off the Scottish (down to about 300 m) and Norwegian coasts (150 m).The anti-correlation is most pronounced in individual years prior to the 90s and still persists on decadal time scales.This preliminary finding suggests that a wider North Atlantic scope must be adopted to unravel the origin of the North Sea tidal range changes and the observed trends in particular.

Summary and conclusion
We have shown that the tidal range in the southwest and the southeast of the North Sea is characterized by a dipole-like pattern between 1958 and 2014, indicating that different forcing mechanisms of shelf-wide or larger spatial character may have been present.To separate these processes, and treat both trends and short-term variability in a unified framework, a PCA-based method was applied to 70 monthly time series of tidal range throughout the North Sea between 1958 and 2014.Data gaps were filled by the statistical method of Ordinary Kriging.A special property of the Kriging procedure is the conservative nature of its estimates at query points, resulting in under-rather than an over-estimation of the general system behavior with regard to trends and PCs.We were able to detect two large-scale signals and explain about 69% of the overall variability in the study area.We attribute the remaining variability of 31% to local effects, which vary widely; they may be absent or could well cause over 50% of variability at an individual tide gauge.In the overall variance, the maximum contribution of a single local effect is at 4%, the average is below 0.4%.
The second PC represents a large-scale barotropic signal and accounts for the negative trends in the UK area (up to -2.3 mm/yr).This mode of variability has a North Atlantic extent, as shown by supplementary analysis of tide gauges in Canada, Reykjavik, and the European Atlantic coast.
Correlations across the basin are high (0.5-0.7) and are caused by common oscillations on time scales between 6 and 24 months.By detecting the same barotropic signal in the shallow-water model of Arns et al. (2015a, b), and eliminating suspects that are not part of the model input or physics, we conclude that only sea level rise and meteorological forcing remain as possible causes.However, no linear correlations with these parameters were found, implying that nonlinear interactions must be present.A further indication for the presence of shallow water effects is the severe weakening of the signal as the tidal wave advances from the relative deep water at the UK into the shallow water areas at the southern and the eastern boundaries of the North Sea The absence of PC1 in the barotropic model and its confinement to the southern North Sea coast has prompted us to hypothesize that local stratification changes exert a strong influence on the tidal range in shallow water at various time scales.By analogy to the known seasonal tidal cycle in the area (Müller et al., 2014), we argue that a stronger pycnocline, possibly lasting over longer periods, stabilizes the water column against turbulent dissipation and allows for higher tidal elevations at the coast.The qualitative and quantitative agreement between inter-annual PC1 changes and an empirically derived stability index is certainly tentative, yet it provides an attractive first-order target for more systematic data analysis and numerical modeling.Further insight into the nature of large German Bight tidal range changesparticularly the underlying trendscould be furnished by a regional general circulation model with realistic background flow and open boundaries to the North Atlantic.

Figure 1 :
Figure 1: Time series of mean annual high and low tidal water levels for three exemplarily selected stations in the German Bight.For illustration purposes all records are shown with different artificial vertical offsets.The increase in the tidal range is illustrated for the three sites as grey shaded areas between high and low water level time series.
Furthermore, an existing barotropic tide and surge model of the North Sea and the adjacent Atlantic Ocean developed by Arns et al. (2015a, b) was updated and used to simulate total water levels from 1958 to 2014.At the open boundaries, we used the Technical University of Denmark DTU10 ocean tide model (Cheng & Andersen, 2010) as tidal input, and the MSL reconstructions of Wahl et al. (2013) were employed in order to incorporate the effects of rise in MSL.The entire model domain was forced with the 20th Century Reanalysis (20CR) data set of the US National Oceanic & Atmospheric Administration (NOAA) and the Cooperative Institute for Figure 3: Changes in tidal range before (a) and after (b) applying ordinary kriging and removing the nodal cycle.
suggest that these local exceptions are mainly caused by anthropogenic interventions such as the building of the Afsluitdijk at Den Oeverbuiten or dredging and dike constructions near to Büsum, as they provide different evidence to nearby stations and anomalies in their time series of tidal range coincide with known local anthropogenic interventions.From the aforementioned findings, we conclude that widespread and statistically significant secular changes in tidal range occurred around large parts of the southern North Sea between 1958 and 2014, although locally interrupted by opposing signals at individual sites.Furthermore, we note contrasting and dipolelike trends along south-western (significant negative values) and south-eastern margins of the North Sea (significant positive values).It remains to be critically noted that the changes in the tidal range at some individual tide gauges could also be instrumental.However, due to the largescale and the spatial homogeneity of the patterns, this cannot be causal for the overall picture.

Figure 4 :
Figure 4: Linear trends of tidal range between 1958 and 2014.Trends at measured sites are shown as dots with a black edge.Dots in between stations are based on Kriging.

Figure 5 :
Figure 5: (a) Linear trends of the M 2 and (b) S 2 tidal constituents between 1958 and 2014 (significant trends outlined).
while local effects through anthropogenic interventions are retained in the reminder of the lower PCs.The amount of these percentages depends to some extent on the spatial distribution of the tide gauges, which is why it is necessary to consider the PCA results at each tide gauge (Figure6-c/d, 7d).PC1 describes an increase in tidal range over time, as evident from its positive slope and the consistently positive values of the associated coefficients at all sites (Figure6-a).The magnitudes of the coefficients reveal that the signal represented by PC1 increases as one travels counterclockwise throughout the basin reaching its strongest expression in the German Bight.PC2 exhibits a negative trend and is most pronounced in the area of the southeastern coast of the UK.The coefficients of PC2 change sign from positive values along the UK coast to negative values in the area of the German Bight (Figure6-b).Similar to the trends of measured tidal range (Figure4), a dipole-like temporal evolution with a node in the area of the English Channel is detected.In general, PC1 accounts for the increase in tidal range in the German Bight and PC2

Figure 6 :
Figure 6: Results of the PCA.(a) Shown are time series of PC1, and PC2 and (b) their corresponding spatial patterns.Panels (c) and (d) map the correlations between observations and PC1 (c) and PC2 (d) for each site.
over the period 1958 to 2014 were used.To facilitate a rigorous comparison with our in situ data, simulated time series at the locations of the 70 tide gauge stations were extracted.A PCA revealed that the PC2 pattern is represented well in the simulated data.We find similarly high correlations between the model-based PC and the observations of the Atlantic tide gauges.While the mean correlation of the European tide gauge records (Figure 8-b) with North Sea PC2 from observations is 0.70 (p<0.05), it is only marginally lower with the barotropic model outputs (r=0.66).If the simulated signal is removed from the model, the correlation becomes insignificant and even disappears at most sites.In consequence, PC2 must be driven by a process initially included into the boundary conditions from the numerical model.Since we have used a barotropic formulation without buoyancy forcing and thermodynamic calculations, we can further infer a purely barotropic relationship.Amongst the possible relevant factors, the tidal input to the model can safely be neglected.The DTU10 tide model consists of ten tidal constituents, stationary in time and modulated only by the 18.6-year nodal cycle.The high correlations on the east coast of the UK and in the North Atlantic are unrelated to this forcing, since purely tide-induced changes would be periodic and present in the remaining parts of the North Sea.

Figure 8 :
Figure 8: (a) Extended network of tide gauges with additional stations shown in red, (b) correlations of all tide gauges (except 5-7) with PC2 and (c) comparison between measured and reconstructed values of tidal range at the newly added tide gauges 1-7.The reconstruction in c) is based on PC2, and the numbers in parentheses indicate the respective correlation.

Figure 9 :
Figure 9: (a) Vertical profiles of potential density as averaged over all query points in (b) at depths from 0 to 35 m for the years 1958 to 2013 (black), the year 1995 (blue) and the years 1998 (red).The two selected years feature the greatest deviation from the mean density profile.

Figure 10 :
Figure 10: (a) Spatially averaged density profiles (0-35 m) from the query area in Figure 9-b spanning the period 1958 to 2013.(b) Comparison between PC1 changes and the stability index (see main text), where both time series were scaled by their standard deviation and adjusted for long-term trends.

Table 2 :
Measured and reconstructed trends in tidal range and explained variance of the different regions.