A Combined Modeling and Measurement Approach to Assess the Nodal Tide Modulation in the North Sea

The Abstract The correct representation of the 18.61-year nodal tide is essential for an interpretation of the evolution of mean sea level, as errors cause misleading bias. The nodal tide is currently estimated by applying correction factors in harmonic analysis, which are derived from the equilibrium tide. From the equilibrium tide, correction values f for amplitude and u for phase are determined, which alter lunar tidal constituents, depending on the nodal cycle. This approach has proven to be valid for many tide gauges, even though the impact of the nodal tide in shelf seas has been shown to differ from their theoretical correction value. Hence, tidal constituents from tide records in the North Atlantic shelf region were analyzed for their nodal amplitude and phase lag with a new multiple, nonlinear regression approach, which is able to approximate the nodal modulation quantitatively and its agreement to the theoretical equilibrium tide. Results show an overestimation of the lunar M 2 and N 2 constituents by the equilibrium of more than 2.7% in the Wadden Sea, while O 1 and K 2 are underestimated by 1–4.6%, which would produce an error of 2–5 cm for example, in the German Wadden Sea. Additionally, a process-based model of the North Sea was applied at the diurnal minimum and maximum of the nodal cycle to calculate a spatial distribution of f and u . Results confirm the spatially varying nodal satellite modulation in friction dominated, shallow water regions.

of the nodal cycle variations on the tidal potential of for example, M 2 is assumed to be 3.7% on global average (Haigh et al., 2011) or 2.2 cm in global mean amplitude (Baart et al., 2012). However, since the actual distribution of land masses does not correspond to the assumption of the equilibrium tide theory of being a purely oceanic earth with a constant water depth, local effects of the nodal cycle can vary significantly and the greatest values occur in diurnal regions with tidal ranges of ≥ 4 m (Haigh et al., 2011). Spatial variations also occur not only due to variable water depths but also due to nonlinear frictional, tidal-tide interactions (Ku et al., 1985), changing energy propagation or/dissipation, and other nonlinear effects (Jay et al., 2015). The magnitude and spatial dependence of these changes already make it evident that the nodal cycle cannot be neglected for water level analysis based on tide gauge data, as it would cause misleading bias. Pugh (1987) also emphasized that the nodal cycle is significant, but difficult to separate for MSL determinations. Therefore, the nodal component of a lunar constituent must either be eliminated or equalized by considering only full nodal periods. This applies not only to tidal analysis but to MSL studies and tidal high and low water level analyses as well, especially when trends are estimated (Jensen et al., 1988(Jensen et al., , 1992. The Dutch coast is an example for this, where a nonconsideration of the nodal cycle conceals the changes in MSL (Baart et al., 2012). This necessity occurs not only in reconstructions, but also in projections. In the case of reconstructions, a correction is necessary for both tidal analyses and MSL determinations from gauge data in order to obtain unbiased results.
Even though long tide records do exist today, a reliable detection of the nodal tide signal remains difficult, for example, due to data quality limitations, noise (Rossiter, 1967;Trupin & Wahr, 1990), or overlying trends (Woodworth & Blackman, 2004;Woodworth et al., 1991). The commonly applied correction for the nodal tide in harmonic analysis (e.g., UTide by Codiga, 2011) is based on assumptions from the equilibrium tide theory (Pugh, 1987). Tidal constituent amplitudes from harmonic analysis are modulated by a percentage-wise adjustment (f, i.e., nodal modulation) and a nodal phase lag (u), respectively. These correction parameters are considered to be globally constant (Trupin & Wahr, 1990;Woodworth, 2012). The consideration of the nodal modulation (also nodal satellite variation Haigh et al., 2011) in amplitude and phase is therefore performed by applying the stationary correction parameters f and u (also f − u correction following Pugh, 1987) to lunar tidal constituents. Nodal amplitudes and phase lags, as well as nodal satellite variation, have been investigated globally (Cherniawsky et al., 2010;Menéndez & Woodworth, 2010) and regionally for example, in the North Sea (Amin, 1985;Hansen et al., 2015;Jensen et al., 1992;Woodworth et al., 1991), the Gulf of Maine (Ray, 2006), the Mediterranean Sea (Shaw & Tsimplis, 2010), the Chinese Sea (Feng et al., 2015), and the western coast of Australia (Amin, 1993) to name only a few. Evaluations of the f − u correction parameters have been published by Cherniawsky et al. (2010), who carried out a harmonic analysis of satellite sea surface elevation data to determine nodal satellites. Their analyzed nodal amplitudes overestimated the equilibrium tide parameters f and u, especially when dealing with small amplitudes. Contrary to this, an analysis of tidal constituents in the Chinese seas has shown an underestimation of the nodal modulation parameters for M 2 and N 2 , while O 1 and K 1 agree well with the theoretical values (Feng et al., 2015). In the Mediterranean Sea, the nodal variation agrees well overall with the equilibrium assumption (Shaw & Tsimplis, 2010). Hence, many findings either challenge or agree with the statements from Trupin and Wahr (1990), who stated that the amplitude and phase are close to their equilibrium amplitude and phase. However, it must be noted, that regional exceptions, for example, for the North Sea, were made in their results.
Thus, as the nodal modulation does not follow the equilibrium tide theory consistently across all coastal waters, the first major aim of this study is to develop a method to quantify nodal correction parameters f and u, from tide records. We use a multiple, nonlinear regression approach to calculate the nodal modulation at gauges in the European North Sea from tidal constituents, before comparing them to the equilibrium value. The second goal of this study is the identification and description of the processes causing the deviation from equilibrium modulation. In order to assess the spatial distribution for f and u, numerical simulations are carried out. In fact, Woodworth (2012) already suggested a barotropic approach for 19 consecutive years without explicit loading as a way to advance the topic of the accuracy of the current nodal correction formulation. In this study, we deploy a three-dimensional, astronomically forced, numerical model of the North Sea at the diurnal minimum, and maximum of nodal modulation in order to obtain spatial distributions through differences. These results provide detailed, spatial information about the regional nodal modulation, and phase lag. Based on these results, a recommendation on how to consider nodal satellite variation in the future can be made, which is the third major aim of this study.
This study is organized as follows: The applied data and the preprocessing are discussed in the first part from Section 2.1-2.3. Afterward a brief description of the numerical model and its validation are given in Section 2.4. Section 3.1-3.4 describe the results from the multiple nonlinear regression fitting approach, before the outcome of the numerical simulations is described in Section 3.5. Results are interpreted and discussed in Section 4.

Nodal Tide Modulation
Modern harmonic tidal analysis extracts amplitude and phase of tidal constituents, defined by given astronomical frequencies, from tide records. The astronomical frequencies (i.e., tidal constituents) are defined by their lunar or solar origin and their reoccurance interval, hence the semidiurnal moon tide is called M 2 . The harmonic analysis (i.e., satellite method) follows the development of tidal potential theory (Doodson, 1921(Doodson, , 1928 and is well documented with its modern formulations (Foreman et al., 2009;Godin, 1972;Pugh, 1987). Tidal constituents underlie interannual, annual, and perennial fluctuations, which have an influence on the results of a harmonic analysis (Feng et al., 2015;Gräwe et al., 2014;Hansen et al., 2015;Müller et al., 2014). In this study, we will focus on fluctuations with an interannual frequency. These are most prominently represented by the nodal tide (18.61 years) and the lunar perigee (8.85 years), while the solar perigee is negligible for practical applications due to its long timescale of 20,392 years. For a full review on the scientific background to the nodal tide and the lunar perigee, we recommend Haigh et al. (2011). In order to consider the nodal tide variation, tidal constituents are modified in harmonic analysis based correction parameters from the equilibrium tide theory (Pugh, 1987). The modulation effect amounts to an amplitude adjustment of 3.7% or phase lag of −2.1° for the M 2 or 28.6% and −17.7° for K 2 constituent for example. This correction procedure is referred to as f − u correction. Since tidal records are typically analyzed for less than a full nodal cycle in practice, tidal constituents are corrected for the nodal variation within a harmonic analysis. However, the f − u satellite correction (i.e., Table 1 following Pugh, 1987) assumes that nodal modulation always follows the equilibrium tide theory, which is valid for most gauges (Trupin & Wahr, 1990), but has been shown to be inappropriate locally, especially in shelf seas. Table 1 lists the correction values f and u for lunar tidal constituents. Note, that amplitude and phase lag are maximal for the diurnal constituents at N = 0° for M f , Q 1 , O 1 , K 1 and K 2 , which would be in June 2006, and for M 2 , N 2 at N = 180° which would be for example, in October 2015. Thus, if one imagines an exemplary mean amplitude of 1 m for M 2 , the f − u correction procedure would correct M 2 in June 2006 (minimum) to 0.963 m and in October 2015 (maximum) to 1.037 m, giving M 2 a nodal amplitude of 0.037 m. If K 2 had a mean amplitude of 1 m, its amplitude would be 1.286 m in June 2006 and 0.714 m in October 2015.

Tide Records
This study utilizes 31 North Sea and North Atlantic tide records from Germany, Great Britain, France, the Netherlands, and Denmark as a data basis for harmonic analysis. Figure 1 lists all gauges and the period of valid data hereinafter. The longest tide records are located in Brest, France (1845-2019), Delfzijl, the Netherlands (1879-2019), Esbjerg, Denmark (1889, Newlyn, Great Britain (1915), and Cuxhaven, Germany (1918-2018. Most Dutch and British tidal records start in the 1970s, while the German measurements start predominately in the late 1990s. Data were carefully controlled by visual inspection, checked for anomalies and harmonized to equidistant hourly values, as required for harmonic analyses (Codiga, 2011). Otherwise, tide records remain unchanged, thus they were not cleared from MSL trend, surge or noise. In order to determine the nodal component, a time series of at least 19 years is required and we define, in HAGEN ET AL.

Nodal Tide Fitting
In literature, different techniques are applied to extract the nodal tide from water level signals. Most recently, the quantile fitting method (Woodworth & Blackman, 2004) was applied to gauges worldwide showing clear nodal signals in 371 of 527 gauges (Peng et al., 2019) in the 90% quantile. The harmonic analysis for the nodal constituents has also been carried out globally (Cherniawsky et al., 2010), using a 16 years long satellite altimetry data set. Additionally, multiple, nonlinear regressions of annual tidal characteristic values (Jensen et al., 1988(Jensen et al., , 1992Woodworth et al., 1991), mean sea level (Baart et al., 2012;Hansen et al., 2015), and tidal constituents (Feng et al., 2015;Shaw & Tsimplis, 2010) were conducted to obtain information about the amplitude and phase of the nodal tide.
The results of the quantile method for tidal characteristic values or their percentiles as well as water level quantiles in the western English Channel and the Wadden Sea were inconclusive in this study, as neither nodal nor lunar perigee signal could be detected. We suspect that frequent wind and storm surge events as well as strong shallow water effects deform the signals and limit the applicability of the quantile analysis method. Thus, the nodal modulation model is applied to tidal constituents. The model is a nonlinear, least square fitting approach shown in Equation 1 (Feng et al., 2015;Jensen et al., 1992;Peng et al., 2019) including an acceleration trend term (Baart et al., 2012). We chose an annual analysis period (January to December) for a reliable estimation of true (Pugh, 1987) tidal constituents and to minimize the effect of interannual variation of partial tides in shelf seas (Gräwe et al., 2014;Menéndez & Woodworth, 2010;Müller et al., 2014).
HAGEN ET AL.

10.1029/2020JC016364
4 of 17 In Equation 1, H(t) represents an annual value at time t, a 0 denotes the initial tidal constituent, a 1 accounts for the linear and a 2 for an accelerated trend with a time lag φ 1 . The nodal tide term is represented by the amplitude a 3 , the nodal frequency ω and a phase lag φ 2 . The variable t represents the time in Julian years (365.25 days).
To account for the rapid change in tidal constituents from the 1980s (Haigh et al., 2019;Müller, 2011), a time lag φ 1 has been added to the acceleration term a 2 (Baart et al., 2012;Houston & Dean, 2011). The goodness-of-fit between data and the nonlinear model estimation is described by the coefficient of determination R 2 . R 2 compares the dependence of two data sets by identifying the fraction of the signal variance, which is explained by the regression model, with one being a perfect estimation. Following Peng et al. (2019), we chose results with a value of R 2 ≤ 0.5 to be statistically insignificant in this study. Note, that fitting of the lunar perigee cycle was attempted by adding an additional term to Equation 1 analog to a 3 with a period of 4.4 years. However, results rarely produced statistical significance within the entire study area for the lunar perigee, which is consistent with previous results (Haigh et al., 2011;Menéndez & Woodworth, 2010). Whenever there was a satisfying agreement for the lunar perigee (R 2 ≥ 0.5), its amplitude was more than 10 times smaller than the nodal amplitude. Therefore, the lunar perigee has not been considered in the following.
Tidal constituents consist of a radiational and a gravitational component. The radiational contribution results from periodic annual, semi-annual and diurnal meteorological phenomena involving variations in temperature, atmospheric pressure and wind variation. Therefore, only the gravitational component is influenced by the nodal tide. Though there have been efforts to separate radiational and gravitational components (Feng et al., 2015;Zetler, 1971), our results do not show an obvious disturbance from radiational influence. Seasonal M 2 modulation (i.e., Section 3.2) has not revealed a noticeable intraannual variation, which is why we suggest, that the radiational component is compensated by the long time-scales. Nevertheless, this imposes a simplification and therefore possible limitation of our results.

Hydrodynamic Model
Even though the North Sea is monitored by one of the most closely meshed measurement networks worldwide (Sündermann & Pohlmann, 2011), information between tide gauge locations and far-off the coast is limited. To close these gaps, this study deploys a numerical North Sea model, which is located on the European continental shelf in the northeastern Atlantic ( Figure 2). The majority of the unstructured grid cells is located within the German Bight and the Dutch Wadden Sea, as complex bathymetry with steep slopes and shallow embankments requires high grid resolution to reproduce realistic energy dissipation (Rasquin et al., 2020). 203,000 horizontal elements have been used with an increasing grid resolution from 7.5 km in the open North Sea to 350 m in the Wadden Sea to less than 50 m in the estuaries of the German Bight. The UnTRIM 2 model (Casulli & Walters, 2000) with the novel subgrid approach (Casulli, 2009) has been applied in order to consider complex bathymetry details at low computational cost. By applying a finer subgrid (increasing 4-12 times of the horizontal grid resolution) within the computational grid, the bathymetric information and therefore volume can be estimated with less effort than in conventional grids (Sehili et al., 2014). The hydrodynamic model has been thoroughly validated with measurements (BAW Technische Berichte et al., 2020) and the validation results from for example, the year 2012 show a root mean square error (RMSE) of 3.3 cm/2.8° for the M 2 , 1.1 cm/4.1° for the S 2 and 0.7 cm/3.5° for the N 2 tidal constituent, respectively. Water levels display a RMSE between 8 and 15 cm, which compares well to similar modeling approaches in the North Sea region (Zijl et al., 2015). The model uses a spatially varying bottom roughness, which has been calibrated to optimally fit the M 2 constituent.
Other processes have not been considered (i.e., waves, salinity, temperature, sediments, surge, and meteorology) to isolate astronomical nodal modulation and eliminate meteorological variability. Tidal constituents from the global tide model FES (FES2014 was produced by Noveltis, Legos and CLS and distributed by Aviso+, with support from Cnes [https://www.aviso.altimetry.fr/]) were reconstructed with ut_reconstr (Codiga, 2011) at the open boundary (blue line in Figure 2), using the default f − u correction. The model bathymetry within the German Bight has been adapted from the EasyGSH-DB (http://easygsh-db.org/) project for the year 2006 (BAW Technische Berichte et al., 2019). Bathymetry for the Dutch Wadden Sea was obtained from Rijkwaterstaat, for the English coast from UKHO and for the French coast from SHOM. The remaining data gaps have been filled with data from the EMODnet database (EMODnet Bathymetry Consortium, 2016). Turbulence closure uses a conventional k-ɛ model.

Fitting Evaluation
The dominant tidal constituents in the North Sea are the M 2 , S 2 , N 2 , O 1 , K 1 , and K 2 components. In the following, the tidal constituents from measured tide gauge records given in Figure 1 are analyzed for annual amplitude and phase, before applying the regression model from Equation 1 to these annual analysis results. All harmonic analyses are completed with the UTide algorithm (Codiga, 2011), which is applied using default settings with all tidal constituents available and without nodal correction.
Results of the analysis are shown in Figure 3 for two example locations, which have been chosen because of their different, geographical setting. Brest is located westwards from the North Sea in the North Atlantic Ocean while Büsum is placed within the German Bight near the Elbe and Eider estuary in the Wadden Sea. A nodal signal is present for all constituents at the exemplary gauge stations with the exception of S 2 in Brest and N 2 in Büsum. Slow quadratic trends for the M 2 and S 2 amplitude and phase modulation show at Brest, which demonstrate accelerated growth between 1970 and 1990. This trend is, however, not present at Büsum. To quantify the goodness-of-fit, the R 2 has been calculated for each constituents nodal amplitude and phase lag in Brest and Büsum in Figure 4.
The R 2 values indicate for Brest that the nonlinear fitting method is reliable (R 2 > 0.8) for all constituents but S 2 . In Büsum, the method yields only reliable results for O 1 and K 2 (R 2 > 0.8), while K 1 , M 2 , and S 2 are on the verge of statistical insignificance (R 2 > 0.5). However, visual inspection of the amplitude fitting of for example, M 2 or K 1 in Figure 3 suggests agreement between the fitting approach and available data. In Brest, R 2 is usually between 0.05 and 0.25 higher, while phase and amplitude R 2 are comparable. The fitting of N 2 was unsuccessful in Büsum for phase and amplitude with a R 2 < 0.5.
The R 2 distribution from Figure 4 is extended to all stations and tidal constituents in Figure 5 below. Gauges have been sorted in the direction of the counterclockwise propagating Kelvin wave in the North Sea with a break in Lowestoft, where the tidal waves from Scotland (Lerwick to Lowestoft) and the English Channel (Brest to Ijmuiden) unite toward the Dutch and German Wadden Sea (Vlieland Haven to Esbjerg). In the following, as the agreement between phase and amplitude coincides at most gauges, their R 2 will be discussed together. The index of agreement R 2 is consistently ≥ 0.85 at all locations for O 1 and K 2 (not included), showing that these constituents are steadily extractable from all tide records.
M 2 also demonstrates high values for R 2 with several local exceptions such as the Dutch west coast in the northeastern English Channel (Hoek van Holland to Ijmuiden), West-Terschelling and Cuxhaven, whose R 2 is ≤ 0.5 for the amplitude modulation, even though the phase lag agrees well for all stations. For K 1 , R 2 is between 0.6 and 0.93 with the exception of Dunkerque (0.5). K 1 's R 2 diminishes in the English Channel from 0.83 in Cherbourg to 0.61 Dover before improving slightly in the north eastern English Channel near Hoek van Holland to a mean of 0.7. N 2 shows strong correlation (0.8-0.95) from Lerwick to Dover, before its R 2 decreases in the northern English Channel in Westkapelle near its amphidromical point, similar to M 2 . Further eastwards, N 2 fitting remains inconclusive due to low statistical significance with the exception of Esbjerg where its nodal amplitude is lower than 0.5 cm. S 2 shows sufficient R 2 from the English Channel to the Wadden Sea, contrary to N 2 . Nevertheless, low correlation (R 2 ≥ 0.5) for S 2 can be seen for the amplitude in Lerwick to Dunkerque excluding Immingham, North Shields and Lowestoft. The phase of S 2 could not be fitted consistently, as R 2 is mostly ≤ 0.5, though slight phase modulation takes place in the Dutch and German Wadden Sea. The values for R 2 vary from 0.75 to 0.85 between Vlieland Haven and Büsum.

Seasonal Variation
Since tidal constituents underlie constant interannual variation (Gräwe et al., 2014;Müller et al., 2014), a sensitivity study is performed to quantify the effect of meteorological forcing on the estimation of nodal amplitudes with Equation 1. The method described with Equation 1 is now applied to all tide records in Table 2 for the constituent M 2 in the summer s (May-October) and winter w (November to April) half year as well as a full hydrological year a. In the following chapters, n denotes the nodal satellite of a tidal constituent. The chosen gauges represent a region in the North Sea study area. The English east coast (Aberdeen, Lowestoft), the English Channel (Brest, Roscoff, Dunkerque, Lowestoft), the Dutch Wadden Sea (Harlingen, Huibertgat) and the German Wadden Sea (Alte Weser, Büsum) are investigated for seasonal effects.
Results from Table 2 show, that the seasonal dependence of M 2 's nodal modulation is always lower than ± 0.4 cm. The M 2n,s and the M 2n,w to M 2n,a ratios deviate more from 1 with a R 2 ≤0.75. The nodal satellite for M 2 at robust estimations, such as Aberdeen, Brest, Dunkerque or Roscoff, shows seasonal variation of less than 6%, while weaker R 2 at Harlingen, Huibertgat or Büsum produce an over-or underestimation in summer or winter by a maximum of 25%. As seasonal variation in M 2n 's amplitude shows small variation, no bias due to seasonally heterogeneously sampled data is to be expected, when applying Equation 1. The analysis also shows, that R 2 is usually larger, after the analysis has been performed for an entire year. It is not surprising that we find the approximation of true tidal constituents to become more robust on a longer time timescale as this has already been stated by Pugh (1987).

Nodal Modulation at Gauges
After the nodal signal has been extracted successfully through multiple, nonlinear regression fitting, its amplitude modulation f and phase lag u can be calculated from regression parameters. If R 2 is ≤ 0.5, results have been regarded as statistically insignificant and will not be included in the subsequent analysis. Differences of the calculated modulation from the theoretical f − u correction factors (i.e., Table 1) are given in Figure 6. Positive values represent an overestimation of the theoretical nodal correction factors f and u. For all constituents, the fitted nodal modulation and phase lag agree well with equilibrium expectations at the English east coast (UK-east) and the southern English Channel (until Cherbourg).   amplitude modulation is underestimated throughout all locations except for Cherbourg and Dover. At the east coast of the UK it is underestimated by −2.5% in Aberdeen to −1.5% in North Shields. The westerly stations Newlyn, Brest and Roscoff underestimate K 1 by less than −0.7%, while the maximum underestimation is reached in Westkapelle (−2.8%), Ijmuiden (−3.1%), and West Terschelling (−3.4%). Though N 2 must be modulated by the nodal tide due to the laws of astronomy, its variation could only be detected until Cherbourg in the southern English Channel and at the UK east coast. Successfully fitted N 2 amplitude variation is underestimated slightly on the English east coast from Aberdeen (−0.1%) to Immingham (−1.1%). Gauges in the southern English Channel underestimate N 2 's amplitude modulation by less by −0.4% in Newlyn until −1.8% in Dunkerque. The nodal component of the semidiurnal M 2 is underestimated consistently after the Dover-Calais narrowing in the English Channel. In Westkapelle, the M 2 modulation starts to be underestimated by −1.4% and 0.7% in Dover. When moving further northeastwards, the nodal amplitude modulation of M 2 diminishes into statistical insignificance until reaching the Dutch Wadden Sea. In Vlieland Haven, M 2 modulation is reduced by −2.6%, although within the Wadden Sea itself, M 2 modulation is mostly reduced by −2.0% (Vlieland Haven to Wyk). Harlingen, West-Terschelling and Cuxhaven did not meet the criterion for statistical significance. S 2 shows a similar behavior, as its amplitude is modulated between Dover (2.0%) and Esbjerg (2.0%), peaking near the S 2 amphidrome at Ijmuiden (6.2%). In the Wadden Sea, S 2 is modulated consistently between Vlieland Haven (4.1%) and Büsum (3.5%). Contrarily, the nodal amplitude variation of K 2 is overestimated after the Dover-Calais narrowing by 1.1% in Dover and 2.8% in Westkapelle. The K 2 amplitude variation remains exaggerated by 4.8% in Hoek van Holland, 1.9% in Vlieland Haven, 2.4% in Huibertgat, 4.4% in Alte Weser and 2.6% in Büsum within the Wadden Sea before reaching a smaller overestimation in Esbjerg of 0.9%. Although, the phase lags of all constituents do not deviate more than ± 2.0° at any location, they show consistent departure from their equilibrium modulation value u with an overestimation between 0.5 and 2.0° of O 1 , K 1 , S 2 and M 2 for all gauges from Dunkerque to Büsum.

The Impact of Shallow Water Constituents
When tides travel from the open sea toward coastal waters, embayments or estuaries, energy dissipation through friction leads to the generation of new tidal constituents, which are called shallow water tides. They are related to larger parent constituents, which are astronomically predefined, e.g. M 2 or S 2 . In other words, the generation of shallow water tides drains energy from parent constituents, which underlie nodal tide modulation. Hence, it is not prudent to imply a connection between bottom friction on nodal modulation, especially because this has already been suspected in Feng et al. (2015) and Godin (1986). However, the quantity and amplitude of shallow water tides vary regionally, which is why the shallow water tide hypothesis above may not hold true for all coastal waters equally. As the underestimation of f for K 2 and O 1 HAGEN ET AL. Note. R 2 (coefficient of determination) is given in order to describe the goodness-of-fit for each scenario.

Table 2 Comparison Between the Summer and Winter Half-Year to Quantify the Meteorological Impact on the Nonlinear, Multiple Regression Fitting of Nodal Modulation of M 2 From Section 2.3. s Denotes Summer, w Winter, a the Full Hydrological Year and n a Nodal Satellite
in Section 3.3 coincides with an overestimation of the nodal amplitude modulation of M 2 , N 2 , and S 2 , we suspect interactions between the parent constituents through their shallow water tides to be the dominant driver for these deviations from the equilibrium tide in the North Sea.
To support the shallow water tide hypothesis, a harmonic analysis at representative gauges near the British east coast and in the Wadden Sea was carried out for M 2 and S 2 's shallow water constituents. In Büsum (Wadden Sea), for example, 10 shallow water components with an amplitude larger than 3 cm were present in 2012 of which MSf,M 4 , M 6 , MU 2 , MS 4 , NU 2 , 2MS 6 , and MN 4 are related to a lunar origin. If the constituents, which are either solar or lunar influenced and the meteorological tide MSf are ruled out, MS 4 and 2MS 6 remain. Multiple, nonlinear regression fitting (Equation 1) of MS 4 and 2MS 6 has shown that a clear nodal signal could not be extracted from these components at any of the locations presented in 3. However, HAGEN ET AL.

10.1029/2020JC016364
10 of 17 Figure 6. Differences between the fitted amplitude modulation in percent (left y-axis, black) and phase lag in degree (right y-axis, gray) to the equilibrium nodal correction parameters (i.e., Table 1) for the constituents M 2 , S 2 , N 2 , O 1 , K 1 , and K 2 . Data points with an insufficient R 2 are marked with a red dot. significant amplitude differences of MS 4 and 2MS 6 between the gauges do exist. If the M 2 to MS 4 or the M 2 to 2MS 6 amplitude ratios, respectively, are calculated, large values coincide with S 2 modulation as demonstrated in Table 3.
While the ratios do not exceed 0.021 in Brest, Lerwick, North Shields at the UK east coast, larger ratios are observed in Lowestoft (0.061), Vlieland Haven (0.061), Delfzijl (0.09), or Büsum (0.04). These larger ratios coincide with significant index of agreement values of R 2 , showing a relationship between the modulation of S 2 and large shallow water tide amplitudes. An analog analysis for K 2 has shown additionally, that the overestimation of K 2 (see also Figure 6) coincides with large K 2 to either MK 4 , 2MK 6 , or MKS 2 ratios.

Spatial Nodal Modulation
The amplitude modulation f and phase lag u has been shown to differ significantly from expected equilibrium conditions in measured tide records (i.e., Figure 6). Furthermore, the fitted nodal modulation shows regional tendencies, for example when the modulation of the M 2 in the Dutch side of the English Channel is reviewed. The nodal modulation corresponds to equilibrium values in Cherbourg, is overestimated near the Dover-Calais narrowing at Dunkerque/Dover, before diminishing into statistical insignificance in Westkapelle to Ijmuiden. Another example is the nodal amplitude variation of N 2 , which has not produced any results between Westkapelle and Wyk due to statistical insignificance. Therefore, we computed a complete spatial distribution of the amplitude and phase correction parameters for the North Sea with the numerical model from Section 2.4 to quantify regional variation. As a consequence of the long time scale of 18.61 years, a simulation would require a time period of roughly 19 years. For this reason, the model has only been applied for 10 months with the diurnal amplitude modulation and phase lag minima and maxima reached 5 months into the computation, as other solutions would be computationally expensive. The diurnal phase lag minimum is in October 2001, its maximum in February 2011, while the diurnal amplitude modulation minimum is in June 2006 and its maximum in October 2015, respectively (Pugh, 1987). The model is forced astronomically at the open boundaries by f − u corrected, reanalyzed tidal constituents from the FES 2014b (Section 2.4). Meteorology has been deactivated and the bathymetry is not altered between each simulation. Model results have been interpolated from the unstructured computational grid on a regular 7.5 km grid before the harmonic analysis UTide without nodal correction has been applied. The spatial nodal modulation is then calculated through the absolute amplitude difference between 2006 and 2015, divided by their mean, and the phase lag is calculated by subtracting 2001 from 2011. The resulting values must be halved, as an amplitude is only half of the sinusoidal range.

Assessing Computed Nodal Modulation at Gauges
Before determining the spatial distribution of amplitude modulation f and phase lag u, the model results are compared to the analysis results from tide records in Section 3.3 by calculating f and u at every gauge from Figure 1. The index of agreement (R 2 ) is applied for the observed and modeled (predicted) nodal satellite variation. To quantify the quality of the computed nodal variation, a mean absolute error (MAE) and a root mean square error (RMSE) are given in Table 4. The actual model validation is referred to in Section 2.4.
The validation of the computed amplitude modulation shows R 2 > 0.5 for the semidiurnal constituents N 2 , M 2 , S 2 and K 2 , proving a linear dependency between observed and computed nodal modulation. The semidiurnal MAEs range between 0.5% and 1.5 % and the RMSEs between 0.7% and 1.5 %, respectively. The differences result from a weaker model response to amplitude modulation, especially for K 2 and S 2 , though the HAGEN ET AL.

Table 4 Validation Parameters for Modeled and Observed Nodal Satellite
Variation Correction Parameters f and u prediction of over-or underestimation, of nodal modulation remains correct. The amplitude modulation of diurnal components K 1 and O 1 , however demonstrates lower agreement to the predicted counterpart. The model produces an overestimation of K 1 by an MAE of 1.8% for O 1 and 4.0% for K 1 . The larger error residuals result from an overestimation of nodal amplitude variation in the model. For K 1 especially, the model calculates an overestimation of the amplitude modulation, while observations have shown an underestimation from the equilibrium value. This holds true for O 1 as well, though prediction and observation display an overestimation of O 1 s amplitude modulation with weak R 2 of 0.12. For the nodal phase lag, R 2 ranges between 0 for O 1 and 0.14 for K 2 , which results from the models inability to reproduce a nodal phase variation, which differs from the equilibrium tide. Opposing to the fitting results from measurements, the model computes phase lags, which correlate with the equilibrium value for all constituents. This results in low R 2 for the phase lag as can be seen in Table 4.

Modeled Nodal Variation
The amplitude of a nodal constituent (left), its nodal amplitude modulation (middle) and its phase lag (right), as computed in the numerical model are shown in Figure 7. We will focus on the nodal modulation of M 2 , N 2 , K 2 , and O 1 subsequently, as these constituents have shown strong regional tendencies. In the following, the index n denotes the nodal satellite of a constituent.
The largest nodal amplitude is observed for M 2n with more than 15 cm in the southern English Channel. The amplitude in the Wadden Sea ranges between 2 and 4 cm, while the Dutch west coast and the Danish north coast show amplitudes of ≤ 1 cm. M 2 's nodal amplitude modulation corresponds with the equilibrium value in the southern English Channel and the UK east coast (3.7 %). It is underestimated on the Dutch side of the English Channel (1.5%-2.0 %) and the eastern German Wadden Sea (2.0%-3.0 %). The phase lag u of M 2 shows low deviations between − 2.1° and − 2.4° in the southern North Sea. The amplitude distribution of N 2n is similar to M 2n and K 2n , though the amplitude of N 2n is significantly lower with a maximum of 4 cm in the southern English Channel and 0.5-1.2 cm in the Wadden Sea. N 2n 's amplitude is below 0.5 cm on the western Dutch coast, Northern Frisia in Germany and the Danish North Sea coast, indicating, why nodal amplitude modulation could not be derived from tide records at these locations. The deviation from the theoretical phase lag u tends to be underestimated between 1.1° and 2.5° for N 2 . The amplitude distribution and magnitude of K 2n is again similar to M 2n , with maximum amplitudes of 15 cm in the southern English Channel. K 2 's amplitude modulation is overestimated by 1.5% in the English Channel and the Dutch Wadden Sea and by 3% near the Dutch west coast. UK's east coast and the southern English Channel remain unaffected as seen previously for M 2n and N 2n . The phase lag of K 2n deviates less than 1.5° from its correction value u almost in the entire study area. O 1n has its maximum amplitude in the UK Moray Firth and Thames estuary at 3.8 cm and its amplitude is ≤ 1 cm at the Norwegian coast, the Danish northwest coast and the southern English Channel. The equilibrium tide nodal modulation is well represented near the model boundaries (≤0.5%), before the model reveals an amplification of O 1 amplitude modulation at the English east coast (1.0%), the eastern English Channel (1.5%), the Dutch west coast, the Wadden Sea (2.5%). The nodal phase lag u for O 1n shows small regional deviation with maximum differences in the English Channel near its amphidromy and in the Dutch and German Wadden Sea (1.0°).

Discussion
This study has shown, based on measured tide records, that the nodal correction parameters f (amplitude) and u (phase) for lunar tidal constituents in harmonic analysis do not always follow the equilibrium tide theory, opposing to the statements in Trupin and Wahr (1990). This analysis corroborates results from previous studies, which discussed the nodal correction parameters based on the equilibrium tide assumption (e.g., Cherniawsky et al., 2010;Feng et al., 2015;Godin, 1986;Shaw & Tsimplis, 2010 to name only a few), and extends the current analysis methods for nodal tide estimation by introducing a multiple, nonlinear fitting approach for tidal constituents. Unlike the commonly applied quantile method (Woodworth & Blackman, 2004), the application of multiple, nonlinear regression to annually calculated, tidal constituents accounts for the accelerated change in semidiurnal tides (Baart et al., 2012;Müller, 2011) and is applicable for tide records with a duration of more than 19 years in the North Sea. In contrast to previous research (Godin, 1986;Shaw & Tsimplis, 2010), the nodal modulation of constituents does not follow the equilibrium tide theory strictly, as O 1 and K 2 are underestimated in the English Channel and the Wadden Sea (see Figure 6). M 2 's nodal modulation agrees well with previous findings (Feng et al., 2015;Ku et al., 1985;Woodworth et al., 1991), as it is reduced in the friction dominated areas of the study site such as the Dutch east coast and the Wadden Sea. Nevertheless, the fitting of N 2 in the same areas remained inconclusive, which was associated with the friction induced generation of shallow water tides and third order disturbances (Godin, 1986). The amplitude modulation of the constituent K 1 differs from the other semidiurnal results, as nodal modulation is already reduced in the Northern Atlantic (i.e., Lerwick, Newly, Brest in Figure 6), though its amplitude modulation is consistently reduced at the Dutch east coast and the Wadden Sea. These findings disagree with Godin (1986), who stated that the nodal correction of K 1 , K 2 , and O 1 is appropriate in any case. Moreover, the agreement found in the Mediterranean Seas (Shaw & Tsimplis, 2010) can only be observed outside the English Channel and the Wadden Sea. However, the interpretation of the results of this study must consider, that the amplitude of the nodal satellites is often smaller than 5 cm in the North Sea. Therefore, small error margins resulting for example, from the 95% confidence interval in harmonic analysis may lead to under-or overestimation when reviewing small amplitudes. The amplitude of K 1n for example rarely exceeds 1 cm in the entire study area. The underestimation of K 2 and O 1 coincided with an overestimation of M 2 , N 2 and a nodal modulation of the solar constituent S 2 . For this reason, we suspect interactions within these constituents and their shallow water tides must be responsible. Furthermore, as an overestimated nodal amplitude modulation for K 2 has not yet been documented in literature, other nonlinear effects such as diffraction, reflection or refraction may be present. The deviation from the equilibrium for the nodal phase lag from tide records was considered negligible for all constituents, as the difference rarely exceeds 2.5°, which would correspond to approximately 5 min for the M 2 .
To fortify the hypothesis, that shallow water effects cause the deviation from the equilibrium tide, a spatial distribution of nodal amplitude modulation f and phase lag u was determined by numerical modeling to distinctively identify affected regions as suggested by Woodworth (2012). The modeling approach has shown considerable skill for the amplitude modulation of semidiurnal constituents, but is less suitable concerning the diurnal components K 1 and O 1 and the nodal phase lag of all constituents. We find certain limitations to arise from the accuracy of the numerical model itself, which has been validated to an order of centimeters and minutes (BAW Technische Berichte 2020) and the harmonic analyses timespan, which differs marginally by 2 months between the model and the tide records. The poor agreement concerning the phase lag agreement could be related to the natural bathymetry changes in the North Sea between the years 2001-2015, but these are not included in the model. The link between tidal constituents and bathymetry changes in the Wadden Sea is well established (Jacob et al., 2016;Rasquin et al., 2020), which would suggest, that the deviation of the nodal phase lag in tide records originates from morphodynamic changes.
Still, the question arises, why the diurnal constituent representation did not achieve the same quality as the semidiurnals in the numerical model. Nevertheless, the model results have revealed regionally deviating amplitude modulation in areas with strong friction such as the English Channel or the Wadden Sea. Since meteorological and thus seasonal influences were neglected in the modeling approach, friction remains as the only major possible cause for shifts in nodal modulation. Friction is induced by the shallow water of the continental shelf and complex basin geometry, which becomes increasingly relevant in the English Channel and in the Wadden Sea. Therefore, we conclude, similar to Feng et al. (2015), a friction induced generation of shallow water tides in the inner North Sea, which leads to a reduced amplitude modulation of M 2 and N 2 . However, the amplitude modulation of K 2 , S 2 , and O 1 is also affected, which partially disagrees with the results from the Chinese Seas (Feng et al., 2015), as K 2 s modulation is significantly underestimated while O 1 is overestimated in the North Sea. Since K 2 s and O 1 s deviation coincide with the underestimation of M 2 and N 2 , we again suspect the energy transfer toward shallow water constituents to be responsible. The shallow water tide hypothesis shows most obvious in the modulation of the solar constituent S 2 , which should not be affected by the nodal tide whatsoever. In fact, every time the modulation of M 2 was reduced, the variation of K 2 and S 2 became enhanced. We have shown, that the modulation of S 2 links to high amplitudes of the shallow water tides MS 4 and 2MS 6 (i.e., Table 3), which for MS 4 has also been shown in the Chinese Seas (Feng et al., 2015). An analog harmonic analysis for the parent constituent K 2 has shown a similar relationship to its shallow water components MK 4 , 2MK 6 , or MKS 2 . Even though the relationship between shallow water tides and the deviation of the nodal amplitude modulation has been proven in this study, more work is required in order to quantify the influence of shallow water tides on the nodal tide, specifically. Other nonlinear effects, such as shoaling, reflection, refraction, or diffraction may also disrupt the nodal modulation, especially in the English Channel.

Summary and Conclusion
In sea-level science, an accurate estimation of low frequency tides, such as the lunar 18.61 years nodal tidal cycle is crucial, as the nodal amplitude can be up to 30 cm (Peng et al., 2019). Even though different methods, such as the quantile method (Woodworth & Blackman, 2004), have been used to quantify a nodal amplitude and phase lag from tide records, an accurate approximation is not yet consistently applicable at any geographical location. Nodal modulation is defined as the correction of lunar tidal constituents from harmonic tide analysis through the f (amplitude) and u (phase) correction parameters, which are derived from the equilibrium tide theory (Pugh, 1987). Past studies have shown these correction parameters to be accurate (Godin, 1986;Shaw & Tsimplis, 2010;Trupin & Wahr, 1990), overestimated (Feng et al., 2015) or underestimated (Cherniawsky et al., 2010) depending on the geographical location and analysis method. For this reason, the correction parameters are calculated and compared to the equilibrium approach at various tide records in the North Sea and North Atlantic region to find an approach, which accurately extracts the nodal tide from gauge records. Furthermore, a numerical model is deployed at the diurnal minimum and maximum of nodal amplitude modulation and phase lag to provide spatial information in between the gauge network at the North Sea. The overall aim is to develop and validate a method to extract the nodal signal from tide records and to supply large-scale information on the nodal modulation in shelf seas, such as the North Sea, as sea level sciences rely on accurate assumptions of the nodal cycle.
A multiple, nonlinear regression analysis of annual tidal constituents at North Sea gauges was chosen to approximate nodal amplitude and phase modulation. Results have shown, that the amplitude correction parameter f is significantly overestimated for M 2 , N 2 , and underestimated for K 1 , K 2 , O 1 in shallow, friction dominated parts of the North Sea, although the calculated phase lag coincides well with deviations of less than 2.0°. Additionally, the solar constituent S 2 was shown to be modulated regionally in the Northern English Channel and the Wadden Sea. We support the hypothesis from literature (Feng et al., 2015;Pugh, 1987), who state that energy transfer from M 2 and S 2 toward shallow water tides such as MS 4 or 2MS 6 leads to S 2 's modulation. Shallow water effects also influence other diurnal and semidiurnal constituents such as O 1 , K 1 , K 2 , or N 2 , though more research is needed on this subject. We could not identify the dominant process behind the underestimation of amplitude modulation f for O 1 and K 2 , as well as the overall overestimation of K 1 as similar effects have not yet been documented in literature. A link between K 2 and the shallow water constituents MK 4 , 2MK 6 , or MKS 2 , however, was established. A larger scale could provide more insight about the processes at play, though we suspect shallow water dynamics to be responsible for the overestimation of O 1 or K 2 .
The numerical modeling studies have confirmed, that friction affected areas, such as the English Channel or the Wadden Sea, reduce the nodal amplitude modulation of M 2 , while the variation of O 1 , S 2 , and K 2 is enhanced. The distribution of f is hereby inhomogeneous, as regional and local differences are present due to the variety of generated shallow water tides. An analysis of shallow water constituents of the tide record at Büsum, for example, has revealed a larger M 2 to MS 4 or 2MS 6 ratio when the modulation of S 2 became significant. Thus, we concluded that the current f−u correction should only be applied, whenever the influence of shallow water tides is negligible, as they influence the nodal modulation of lunar constituents and S 2 . The application of u can hereby be regarded reasonable due to low deviations of less than ±2.5°. It must be noted that these recommendations do depend on the field of application and the user-desired degree of accuracy.
Despite the wide acceptance of the f − u nodal correction methodology, it may significantly deviate from the equilibrium in friction affected areas. Additionally, the f − u correction does not consider the modulation of S 2 . Therefore, the f − u correction parameters must be determined appropriately and if necessary, corrected. This process is simplified with the nonlinear, multiple regression of tidal constituents presented in this study, which enables the calculation of accurate f and u correction values. Future work is recommended toward the extension of these regional results to a global scale. This can be performed numerically by global modeling scenarios or analytically by nonlinear, multiple regression analysis of satellite altimetry data. The resulting product would be a globally varying data set providing appropriate f and u correction factors. A correction layer, which includes nonlunar, yet modulated constituents such as S 2 , could stand behalf of the general correction formulation for a more accurate, spatially varying nodal modulation correction. Future work is also recommended towards quantifying the effect of shallow water tides on the nodal satellite variation, especially as many gauges are located in complex coastal or estuarine environments in practice. This would aid sea level science, as the correction of tide records for nodal modulation would be more accurate and help to understand yet unexplained phenomena.