Seasonal variation of the principal tidal constituents in the Bohai Sea

The seasonal variation of tides plays a significant role in water level changes in coastal regions. In this study, 10 seasonal variations of four principal tidal constituents, including M2, S2, K1, and O1, in the Bohai Sea, China, were studied by analysing one-year sea-level observations at E2 and 17-year sea-level observations at Dalian with an enhanced harmonic analysis. At E2, the M2 tidal amplitude and phase lag had annual frequencies, with large values in summer and small values in winter, while the frequencies of S2 and K1 tidal amplitudes are also nearly annual. In contrast, the O1 tidal amplitude increased constantly from winter to autumn. The maxima of phase lags appeared twice in one year for S2, K1 and O1, taking place near 15 winter and summer. The seasonal variation trends estimated by the enhanced harmonic analysis at Dalian are different from those at E2, except for the M2 phase lag. The M2 and S2 tidal amplitudes varied semi-annually and annually, respectively, and were relatively significant at Dalian. The results of numerical experiments indicate that the seasonality of vertical eddy viscosity induced seasonal variations of the principal tidal constituents at E2, while the variations at Dalian were due to the seasonality of stratification and vertical eddy viscosity. 20


Introduction
Although there is no primary seasonal cycle in the moon's orbit, a significant seasonal variation in the principal lunar tidal constituent has been observed and is dominant in coastal and polar regions . In particular, the seasonal variation in the major tidal constituent M2 has received considerable attention (Grä we et al., 2014). Corkan (1934) inferred a seasonal modulation of the M2 tide by analysing several sea-level records near the British coast. 25 Foreman et al. (1995) observed a seasonal cycle of the M2 amplitude at Victoria, which is on the southern tip of Vancouver Island off Canada's Pacific coast. Kang et al. (1995) revealed the seasonal variability of the M2 tidal harmonic constants in the seas adjacent to Korea. Huess and Andersen (2001) found a seasonal variation in the M2 constituent in the northwest European shelf. Kang et al. (2002) investigated the seasonal variability of the M2 tide in the Yellow and East China Seas. Georgas (2012) observed seasonal episodes of significant tidal damping and modulation in the Hudson River estuary. Müller et al. (2014) 30 studied the global seasonal cycle of the M2 tide and found significant seasonal variations in several coastal areas, including the North Sea, East China Sea and Yellow Sea, Sea of Okhotsk and regions of the Banda, Timor, and Arafura Seas north of Australia. Tazkia et al. (2017) found that the M2 amplitude changed markedly between winter and summer in the northern Bay of Bengal.
Several other studies have analysed the seasonal variability of the M2 tide in polar regions. Mofjeld (1986) observed seasonal fluctuations of the tidal harmonic parameters on the north-eastern Bering Sea shelf. Kagan and Sofina (2010) showed 5 that the seasonal variability of tidal constituents was widespread in the Arctic Ocean. Further, Müller et al. (2014) studied the global seasonal cycle of the M2 tide and likewise identified significant seasonal variations in the Arctic regions.
Previous studies have primarily focused on the seasonal variation in the M2 amplitude without considering the seasonality of other tidal constituents and their phase lags (Grä we et al., 2014). Indeed, several studies have investigated the seasonality of several constituents. For example,  studied the seasonal variations of M2, N2, O1 and M4 in the Bohai 10 Sea by introducing astro-meteorological constituents; Devlin et al. (2018) found that the diurnal (K1 and O1) and semidiurnal (M2 and S2) tidal amplitudes and phase lags exhibited strong seasonal variability in the seas of Southeast Asia.
In this study, sea-level observations at one mooring station (E2) and one tidal gauge station (Dalian) in the Bohai Sea were used to investigate the seasonal variability of the principal tidal constituents M2, S2, K1 and O1 using an enhanced harmonic analysis (EHA). The rest of the paper is organised as follows. In Section 2, the sea-level observations in the Bohai Sea are 15 reported and the analysis methods are described. In Section 3, the seasonal variability of the principal tidal constituents is estimated by analysing observations. The mechanisms underlying the seasonal variability are discussed by using numerical experiments in Section 4. Further discussions and conclusions occupy Sections 5 and 6, respectively.

Observations 20
From 0000 UTC 1 November 2013 to 0000 UTC 1 November 2014, total sea levels were observed hourly using a moored pressure gauge accurate to within 5 cm , at E2 station in the Bohai Bay, China (Figure 1). The time series of the total sea levels at E2 is shown in Figure 2a, demonstrating the continuous coverage of the observations. Hourly sea-level data at the Dalian tidal gauge station were obtained from the University of Hawaii Sea Level Center and used. After 1979, Dalian shared position with Laohutan (Feng et al., 2015), so the sea-level data at Dalian were comprised of 25 data from Dalian from 1980-1990and from Laohutan from 1991, as shown in Figure 2b.

Classical harmonic analysis
A sea level is composed of components from different sources (Godin, 1972;Foreman, 1977;Pawlowicz et al., 2002;Foreman et al., 2009) is the total sea level; 0  is the mean sea level; A and g are the amplitude and phase lag (UTC time, the same below), respectively; f and u are the nodal corrections to amplitude and phase lag, respectively; V is the astronomical argument; R is the nontidal component; K is the number of tidal constituents; NNR is the number of non-reference constituents; NR is the number of reference constituents; and NI is the number of constituents to be inferred from the j th reference constituent. 5 The mean sea level, amplitude and phase lag of each constituent can be solved by analyzing a time series of sea-level observations at a specific point using classical harmonic analysis (CHA). With different assumptions and conditions, CHA can be performed using the T_TIDE (Pawlowicz et al., 2002), U_TIDE (Codiga, 2011) or Institute of Ocean Sciences Tidal Package (Foreman et al., 2009). In this study, T_TIDE, in which the astronomical argument varies linearly and the nodal correction is performed after least squares fit, is used to realize CHA. 10

Segmented harmonic analysis
Following Foreman et al. (1995), Kang et al. (1995), Müller et al. (2014) and Devlin et al. (2018), sea-level observations are divided into monthly segments by calendar month and CHA with nodal and inference corrections is applied to each monthly segment to obtain the discrete tidal harmonic parameters (i.e., amplitude and phase lag). Then the discrete amplitude and phase lag at each monthly segment are interpolated using cubic spline interpolation to obtain the temporally varying amplitudes and 15 phase lags. This methodology is termed segmented harmonic analysis (SHA). Following Kang et al. (1995), the sea-level observations in every monthly segment are analyzed only when the duration of the observations is greater than 26 days.

Enhanced harmonic analysis
By combining CHA with independent point scheme and cubic spline interpolation, Jin et al. (2018) developed EHA to directly obtain temporally varying mean sea level and tidal harmonic parameters. In contrast, the harmonic parameters are 20 assumed to be constant in CHA and SHA. A MATLAB toolkit, S_TIDE, was released to support EHA by Pan et al. (2018b).
In this study, nodal and astronomical argument corrections are embedded into the least square fit, following Foreman et al. (2009); in addition, the harmonic parameters of the minor tidal constituents are assumed to be constant and computed together with the temporally varying harmonic parameters of the principal tidal constituents to resolve more constituents and retain computational stability. The sea level in EHA is as follows: 25 where I is the number of principal tidal constituents with temporally varying harmonic parameters; J is the number of minor tidal constituents with constant harmonic parameters; and the mean sea level and nontidal component are included in Similar to Jin et al. (2018) and Pan et al. (2018b), the independent point scheme and cubic spline interpolation are used to jointly determine solve the temporally varying and constant harmonic parameters, which are not shown here for brevity. As 5 mentioned in Pan et al. (2018b), the temporally varying harmonic parameters obtained using EHA with different numbers of independent points represent oscillations on different time scales. In this study, six independent points are used to obtain the seasonal variability of the principal tidal constituents.

Results
One-year sea-level observations at E2 were analyzed using CHA with the automated constituent selection algorithm 10 (Pawlowicz et al., 2002). According to the signal-to-noise ratio (Pawlowicz et al., 2002;Matte et al., 2013), M2, K1, S2 and O1 were selected as the principal tidal constituents to be investigated in this study.

Seasonal variability at E2
As shown in Figure 3, the significant constituent near K1 was P1, which was unresolved when analysing one-month observations , while that for S2 was K2 in the semidiurnal frequency band. Therefore, when the monthly 15 analysis was performed in SHA, the automated constituent selection algorithm in T_TIDE was used to determine the analysed constituents; in addition, the unresolved constituents P1 and K2 were inferred from K1 and S2, respectively, with the inference parameters taken from a yearly harmonic analysis of the one-year sea-level observations at E2 (Kang et al., 1995;Foreman et al., 2009). When EHA was used to directly analyse the sea-level observations at E2, the harmonic parameters of M2, K1, S2, O1, P1 and K2 were estimated together; in detail, the harmonic parameters of M2, K1, S2 and O1 were assumed to be temporally 20 varying and those of P1 and K2 were assumed to be constant. Figure 4, the estimated harmonic parameters obtained with SHA and EHA, including the temporally varying amplitudes and phase lags, were nearly equal and near those estimated using CHA, indicating that the temporal variations in the harmonic parameters of the principal tidal constituents at E2 can be reasonably estimated using both SHA and EHA. Based on Wei and Wang (2012) and Zhang et al. (2017), spring, summer, autumn and winter were defined as March to May, June to 25

As shown in
August, September to November and December to February of the following year, respectively. The temporally varying harmonic parameters of the principal tidal constituents showed seasonal variations ( Figure 4). For M2, the seasonal variations were significant; in detail, the amplitude reached maximum in summer and minimum in winter, just as the phase did, as in Müller et al. (2014). The seasonality of the S2 amplitude was not significant, but the estimated results using EHA increased significantly in summer. The temporal variation of the K1 amplitude spanned one period, with maxima and minima in summer and winter, respectively. In stark contrast, the O1 amplitude increased from winter to autumn. The phase lags of the S2, K1 and O1 components varied semi-annually, in which the phase lags were larger in summer and winter and smaller in spring and winter, respectively. In the four principal tidal constituents, only the M2 amplitude had the similar variation trend as the phase lag. 5 The seasonally averaged amplitudes and phase lags of the principal tidal constituents are listed in Table 1. The variation trends of the averaged harmonic parameters of these constituents were the same as those obtained using EHA in Figure 4.
Compared to the annual averages, the mean M2 amplitude increased by 6.90 cm (approximately 9.33%) in the summer and decreased by 6.68 cm (approximately 9.03%) in the winter, close to the estimated values in Foreman et al. (1995)  (7.72%) in the winter and increased by 7.93% (5.91%) in the summer, indicating a nearly annual variation as shown in Figure   4. The mean O1 amplitude in the summer increased by 3.45 cm compared to that in the winter. Only the M2 phase lag in winter was smaller than the annual average, and the phase lags in both winter and summer were larger than the corresponding annual average for the other three principal tidal constituents.

Seasonal variability at Dalian 15
The multiyear data at Dalian shown in Figure 2b were analysed year by year. In each year, one-year sea-level observations were analysed using CHA, SHA and EHA with settings similar to those used for data from E2. As shown in Figure 5, P1 and K2 were the significant constituents unresolved in the monthly analysis, just like at E2. Therefore, P1 and K2 were inferred from K1 and S2 in SHA and taken as minor constituents with constant harmonic parameters in EHA. The estimated harmonic parameters from various years were then averaged  and are shown in Figure 6. The estimated harmonic 20 parameters using both the SHA and EHA were near to those obtained using CHA, showing that the estimated results were reasonable. In addition, the estimated harmonic parameters obtained using EHA were much closer to those obtained using SHA for data from Dalian than those from E2.
The variation trends of the harmonic parameters estimated using EHA at Dalian were different from those at E2, except for the M2 phase lag ( Figure 6). The M2 amplitude at Dalian varied semi-annually, with large values in summer and winter and 25 small values in spring and autumn, respectively. The S2 amplitude had significant annual variation, with maximum in winter and minimum in summer, which is opposite of the variation trend of S2 amplitude at E2. The K1 amplitude was nearly constant from winter to spring and increased during the summer. The O1 amplitude reached the minimum in the winter and summer while increasing slightly in the spring and autumn. The estimated S2 phase lag reached the maximum in the spring with small variation. The K1 and O1 phase lags trended the same, increasing in the winter and summer while decreasing in the spring and 30 early autumn.
The averaged amplitudes and phase lags of the principal tidal constituents at Dalian, as listed in Table 2, show a seasonal variation that was generally smaller than that at E2. All of the seasonal changes of the principal tidal constituents were less than 1.80 cm, a variation that only appeared in the S2 amplitude at E2. In addition, all of the seasonal changes of the phase lags were less than 2.20°, while the S2 and K1 phase lags at E2 changed by at least 5.00° in the summer and winter, respectively.
The rate of change in the M2 amplitude at Dalian was less than 1%, which was significantly less than that at E2. The rates of changes in all phase lags were less than 1%, except for the M2 phase lag, which was larger than 2% in both summer and winter compared to the annual averages, and larger than those at E2. An increase in the S2 amplitude of 5.34% occurred in the winter, 5 larger than its decrease in the winter at E2.
In summary, the harmonic parameters of the principal tidal constituents at E2 and Dalian varied seasonally, and the seasonality of the tides at E2 was significantly different from that at Dalian. The amplitude of the principal tidal constituent M2 at E2 varied annually, while that at Dalian was semi-annual. The M2 phase lags at E2 and Dalian had the similar variation trend, with lager values in summer and small values in winter. The S2 amplitude in winter at E2 was less than that in summer, 10 which was opposite that at Dalian. The K1 amplitude at E2 had an annual frequency, with large values in summer and small values in winter, while the O1 amplitude increased steadily. In contrast, the variations of the K1 and O1 amplitudes at Dalian were small. The maxima of the S2, K1 and O1 phase lags at E2 appeared twice, like those of K1 and O1 and different from that of S2 at Dalian.

Mechanisms for the seasonal variability 15
Several previous studies have investigated the seasonal variability of the M2 amplitude. Three main mechanisms have been proposed: 1) Seasonal variations of the mean sea level. Corkan (1934) related the seasonal modulation of the M2 tide near the British coast to seasonal variations of sea level and atmospheric pressure. Tazkia et al. (2017) pointed out that the seasonal variability of the sea level generated by many processes can induce a seasonal variation of the M2 tide, as tidal wave propagation was 20 controlled by water depth on the first order.
2) Seasonally varying stratification. Foreman et al. (1995) presumed that the seasonal variability of the M2 amplitude at Victoria, Canada was induced by the changes in stratification due to seasonal variability in estuarine flow. Kang et al. (2002) used a two-layer numerical model to investigate the baroclinic response of the tide and tidal currents in the Yellow and East China Seas, and found that seasonal stratification had several noticeable effects on the tides, including varying degrees of 25 current shear, frictional dissipation, and barotropic energy flux. Müller (2012) indicated that in shallow seas, seasonal variations in stratification were a major factor for the observed seasonal modulation in tides. Müller et al. (2014) pointed out that the seasonal changes in stratification on the continental shelf affected the vertical profile of the eddy viscosity to further cause the seasonal variability of the M2 tide.
3) Seasonally varying ice coverage. St-Laurent et al. (2008) proposed that the significant seasonal variations of the M2 30 surface elevation in all regions of the Hudson Bay system were essentially caused by under-ice friction. Georgas (2012) pointed out that the seasonal episodes of significant tidal damping (reductions in tidal amplitudes by as much as 50%) observed in the Hudson River estuary were primarily caused by the under-ice friction as well. Müller et al. (2014) found that the frictional effect between the sea ice and ocean surface layer led to the seasonal variability of the M2 tide.
Other mechanisms, including long-term changes in the tidal potential (Molinas and Yang, 1986), interactions with other physical phenomena (Huess and Andersen, 2001;Pan et al., 2018a), changes in the internal tide with corresponding small changes in its surface expression (Ray and Mitchum, 1997;Colosi and Munk, 2006), as well as a number of other more 5 technical reasons, may also change the M2 tidal amplitude on various time scales. The above reasons have been presented or discussed in Woodworth (2010), Müller (2012), Müller et al. (2014) and Tazkia et al. (2017).

Design of numerical experiments
The Bohai Sea in north China freezes to varying degrees every winter for approximately 3-4 months (Su and Wang, 2012). According to the back-effect connection of the coastal shelf and open ocean via resonance mechanisms (Arbic et al., 10 2009;Arbic and Garrett, 2010), sea ice may be important to the seasonality of principal tidal constituents. However, Zhang et al. (2019) found that the damping effect of sea ice on the astronomical tides was almost negligible in the Bohai Sea employing numerical experiments with a three-dimensional ice-ocean coupled model. Therefore, ice coverage was not considered in this study. Several numerical experiments (Exp1-Exp4) were carried out to simulate the four principal tidal constituents in the Bohai Sea under different conditions using MITgcm (Marshall et al., 1997), testing the influence of seasonal variations of 15 mean sea level and stratification on the seasonal variability of the principal tidal constituents.
Identical model settings were used in all of the numerical experiments and described as follows. The simulation area was the Bohai Sea as shown in Figure 1b. The horizontal resolution was 2′×2′, with 16 layers in the vertical direction with thicknesses ranging from 2-5 m. The four principal tidal constituents M2, S2, K1, and O1 were implemented as tidal forcing at the east open boundary, whose data were predicted using the constant harmonic parameters extracted from the TPXO model 20 (Egbert and Erofeeva, 2002). Sea surface boundary conditions were not considered. The horizontal eddy viscosity coefficient was set to 1.0×10 3 m 2 /s, and the quadratic bottom drag coefficient was set to 1.3×10 -3 (Wang et al., 2014). The integral time step was 60 s and the total simulation time was 60 d. The results of the final 30 d were used to calculate the harmonic parameters using CHA. Table 3. In Exp1, the simulation started 25 from 0000 UTC 1 January 2014, while the simulation started from 0000 UTC 1 July 2014 in Exp2-Exp4. In Exp1, horizontally homogeneous profiles of the initial temperature and salinity (Figure 7) were extracted from the HYCOM global analysis results in winter, while those in summer were used in Exp2-Exp4. The vertical eddy viscosity coefficient was specified directly and no turbulence closure schemes were used. In Exp1, the vertical eddy viscosity coefficient was set to 2.0×10 -3 m 2 /s through a trial and error procedure. According to Müller et al. (2014), the eddy viscosity in summer was reduced by orders of magnitude 30 compared to well-mixed conditions in winter, as the stratification stabilized the water column. Therefore, the vertical eddy viscosity coefficient was decreased by one-half in Exp3 to test the influence of the vertical eddy viscosity caused by the stratification. As shown in Figure 8, monthly means of the low-pass sea levels, filtered using a cosine-Lanczos filter with a high frequency cut-off of 0.8 cpd, were nearly equal to the estimated mean sea level using SHA. They exhibited the same variation trend as those obtained using EHA, with large values in summer and small values in winter. Given the difference between the averaged mean sea level in summer and that in winter, Exp4 included 1 0.2-m increase of water depth to test the influence of mean sea level.

Results 5
The simulated harmonic parameters of the four principal tidal constituents in the numerical experiments and those obtained from observations at E2 and Dalian are shown in Figure 9. The simulated harmonic parameters were a little far from the observed results, except the M2 amplitude at E2 simulated in Exp1 and that simulated in Exp2, possibly because the constant bottom drag coefficient was used (Wang et al., 2014). Still, the differences between the simulated results in the different numerical experiments can be used to test the influence of potential factors on the seasonal variability of the principal tidal 10 constituents.
The observed amplitudes at E2 in summer were larger than those in winter for all four principal tidal constituents, as shown in Figure 9, but the simulated amplitudes in Exp2 were nearly equal to those in Exp1. In contrast, both the decreased vertical eddy viscosity coefficient in Exp3 and the increased mean sea level in Exp4 increased the amplitudes for all principal tidal constituents. The increases of the observed M2 and S2 amplitudes at E2 from winter to summer were 13.58 cm and 2.62 15 cm, respectively, while those were 13.34 cm (2.08 cm) and 2.75 cm (0.56 cm) for simulated results in Exp3 (Exp4) compared to those in Exp1. In addition, the increases of the observed K1 and O1 amplitudes were also captured better by the simulated results in Exp3 than those in Exp 4, as shown in Figure 9. Therefore, the seasonally varying amplitudes of all principal tidal constituents were primarily impacted by the seasonal variation of vertical eddy viscosity. Differences between the simulated phase lags in Exp3 and those in Exp1 indicated that the seasonal variation of vertical eddy viscosity caused the same trends as 20 the changes between the observed phase lags in summer and those in winter for all principal tidal constituents except K1. In contrast, the changes in stratification in Exp2 and mean sea level in Exp4 only captured the variation trend of the S2 and O1 phase lags, respectively. The aforementioned results demonstrated that seasonal variation in the vertical eddy viscosity was the most important mechanism influencing the seasonal variability of principal tidal constituents at E2.
The seasonal variation of the S2 amplitude was the most significant at Dalian, but the decrease in the simulated result 25 from winter (Exp1) to summer (Exp2) was less than 1 cm and the simulated S2 amplitudes in Exp3 and Exp4 were larger than that in Exp1, indicating the seasonality of stratification as a possible reason. However, the simulated seasonal variation was too weak, possibly because the simple horizontally homogeneous temperature and salinity profiles could not reflect reality.
The water depth is large in the eastern part of Bohai Sea (Figure 1b), so the stratification and ocean circulation were noteworthy and had significant effects on the tides. The increases of the M2 and O1 amplitudes were only captured by the changes in the 30 vertical eddy viscosity coefficient in Exp3 and stratification in Exp2, respectively, while the increases of K1 amplitudes in Exp3 and Exp4 agreed with the observed results at Dalian. The variation lags of the M2 and S2 phase lags were not captured well and the simulated results in Exp2 were the best results, while those for K1 and O1 were best captured by the decrease of the simulated results in Exp3 compared to those in Exp1. On the whole, the seasonal variations of the principal tidal constituents at Dalian were determined by the seasonality of stratification and vertical eddy viscosity.
The variations of the simulated amplitudes from winter (Exp1) to summer (Exp3) in the entire Bohai Sea are shown in Figure 10. The spatial distribution of the variations in M2 tidal amplitude had a strong positive correlation (R=0.98) with that in the S2 tidal amplitude, similar to that for the diurnal tides (R=0.98). Furthermore, the distributions were possibly related to 5 tidal wave propagation as their patterns were similar to the co-phase lines, as shown in Figure 10. For the semi-diurnal tides M2 and S2, the simulated amplitudes in summer were larger than those in winter in Bohai Bay, Laizhou Bay, Liaodong Bay, which was the same as that obtained by analysing the sea level data at several tidal gauge stations in the Bohai Sea in  and that simulated by numerical model in Kang et al. (2002), and smaller than those in winter in the middle of the Bohai Sea. The spatial distribution of the absolute differences between the M2 tidal amplitude in summer and that in 10 winter was similar to that in Müller et al. (2014). For the diurnal tides K1 and O1, the simulated amplitudes in summer were larger than those in winter in Bohai Bay, Laizhou Bay, Liaodong Bay and the middle areas, while smaller than those in winter in the northeast part of the Bohai Strait.

Discussions
In this study, the EHA developed in Jin et al. (2018) and Pan et al. (2018b) was further improved in order to resolve more 15 tidal constituents by adding the minor constituents whose harmonic parameters were assumed to be constant and computed together with the temporally varying harmonic parameters of the principal tidal constituents. The nodal and astronomical argument corrections were embedded into the least square fit to eliminate the influences of nodal cycle and linearly varying astronomical argument. In fact, there have been multiple improvements to T_TIDE in the past decades, such as R_T_TIDE (Leffler and Jay, 2009), versatile tidal analysis (Foreman et al., 2009), U_TIDE (Codiga, 2011) and NS_TIDE (Matte et al., 20 2013). In R_T_TIDE, versatile tidal analysis and U_TIDE, the harmonic parameters (i.e., amplitude and phase lag) are assumed to be constant, although they have some improvements to T_TIDE. However, harmonic parameters are not constant and have multiscale temporal variations, as shown in Corkan (1934), Kang et al. (1995), Müller et al. (2014), Devlin et al. (2018), and so on. Neglecting seasonal variation of tides will introduce significant error in sea-level prediction . Therefore, EHA, in which the harmonic parameters of the principal tidal constituents are assumed to be temporally varying 25 and computed directly within the least squares fit, improved T_TIDE in other ways. In NS_TIDE, the harmonic parameters are also assumed to be temporally varying. However, the temporally varying harmonic parameters are assumed to be functions of river flow and greater diurnal tidal range at the reference station, so it can be applied only to river tides, while EHA can be applied in analyzing any time series. On the whole, EHA used in this study is indeed enhanced than other methods.
The duration of hourly sea-level observations at E2 used in this study was only one year, which is a limitation. Although 30 the rare event might slightly skew seasonal pattern of tides, the seasonal variations of the principal tidal constituents M2, S2, K1 and O1 obtained using EHA were the same as those using traditional EHA, indicating that the results are not unreasonable and can reflect the seasonal variations of tides in the analysis period. The strong seasonal variations of the principal tidal constituents can be captured by the results of numerical experiments. The multi-annually averaged results at Dalian also showed the seasonal variations of the principal tidal constituents, but the results of numerical experiments were not in accordance with the observed results, which may be because only the horizontally homogeneous profiles of the initial temperature and salinity were used and the temporally varying ocean circulation were not considered. 5 The seasonality of the principal tidal constituents has been investigated widely. As shown in Müller et al. (2014), there were significantly seasonal variations in M2 tide in several coastal regions and the maximum annual tide was in July (±1 month) in most of the ocean. However, the spatial and temporal inhomogeneities were also existed and the summer amplitudes of M2 tide were less than those in winter in several areas, as shown in Kang et al. (1995), Müller et al. (2014), Devlin et al. (2018) and Figure 10 in this study. In general, the M2 tidal amplitude in summer was larger than that in winter in many areas, 10 including the Bohai Sea , the North Sea (Huess and Andersen, 2001;Grä we et al., 2014;Müller et al., 2014), most of the Ganges-Brahmaputra-Meghna delta (Tazkia et al., 2017), the seas of Southeast Asian (Devlin et al., 2018), Liverpool (Corkan, 1934, Victoria (Foreman et al., 1995), the western part of the Yellow and East China Seas (Kang et al., 2002), the Hudson Bay and Foxe Basin (St-Laurent et al., 2008), and so on, which was the same as that obtained by analysing the observations at E2 and Dalian in this study. Devlin et al. (2018) found that the diurnal and semidiurnal tidal amplitudes 15 and phases exhibited a high degree of seasonality in the seas of Southeast Asia, and  indicated that the O1 tidal amplitude in summer was also larger than that in winter in the Bohai Sea, which were similar to those concluded in this study. The seasonality of the principal tidal constituents obtained in the study was mainly similar to those in previous studies, but the novel EHA was firstly used to estimate the seasonal variations of the principal tidal constituents and the numerical experiments using three-dimensional MITgcm were performed to explore the physical mechanisms. 20 The seasonal variations of stratification and vertical eddy viscosity and their influences on the tidal amplitudes may be as follows. In winter, the strong northwest Asia monsoon develops a vertically will-mixed condition (Yanagi et al., 2001;Jeon et al., 2014), which will not stabilize the tidal currents and will lose more energy, leading to smaller tidal amplitudes. As the surface heating rate and freshwater discharge increase, the mixing is insufficient to homogenize the input potential energy and cause stratified conditions in summer (Huang et al., 1999;van Haren, 2000), during which the reduced vertical eddy viscosity 25 will increase the tidal amplitudes. However, the S2 amplitude at Dalian was larger in winter and smaller in summer, which is inconsistent with the other principal tidal constituents and should be further investigated in future studies.

Conclusions
In this study, based on one-year sea-level observations at E2 and 17-year sea-level observations in the Bohai Sea, the seasonal variability of the principal tidal constituents was investigated using different methods. When the sea-level 30 observations at E2 and Dalian were analysed, the seasonal variations of all principal tidal constituents obtained using EHA were nearly equal to those obtained using SHA (Figures 4 and 6), indicating that the seasonal variations were not related to the applied methods. At both E2 and Dalian, the principal tidal constituents M2, S2, K1 and O1 exhibited seasonal variations (Figures 4 and 6). The M2 amplitude at E2 varied annually, while that at Dalian was semi-annual. The M2 phase lags at E2 and Dalian had the similar variation trend, with large values in summer and small values in winter. The S2 amplitude in winter at E2 was less than that in summer, which was opposite that at Dalian. The K1 amplitude at E2 had an annual variation, with large values in summer and small values in winter, while the O1 amplitude increased steadily. On the contrary, the variations of the 5 K1 and O1 amplitudes at Dalian were small. The maxima of the S2, K1 and O1 phase lags at E2 appeared twice, which was the same as those of K1 and O1 and different from that of S2 at Dalian.
Through several numerical experiments, the mechanisms of the seasonal variability of the principal tidal constituents were investigated. The seasonal variations of the principal tidal constituents at E2 were caused by the seasonality of the vertical eddy viscosity, while the seasonal variations at Dalian were mainly induced by the seasonality of stratification and vertical 10 eddy viscosity, although the simulated results were not consistent well with the observed results. Therefore, the synchronous simulation of circulation and tides, and a reasonable parameterization scheme to convert the variations in stratification to those in vertical eddy viscosity were essential for precise simulation of the tides when considering the temporally varying harmonic parameters.

Data availability 15
The HYCOM global analysis data is available at http://hycom.org. New version of S_TIDE package can be downloaded from https://www.researchgate.net/project/Adaptation-of-tidal-harmonic-analysis-to-nonstationary-tides. The hourly sea level observations at Dalian are available at https://uhslc.soest.hawaii.edu/datainfo/. The hourly sea level observations at E2 used in this work are available from the authors upon request (xqinglv@ouc.edu.cn).