Effects of Sea Level Rise on Tidal Dynamics in Macrotidal Hangzhou Bay

: Sea level rise (SLR) due to climate change is expected to alter tidal processes and energy transport, disproportionately affecting coastal communities. Utilizing a nested hydrodynamics model, we provided an integrated investigation of tidal responses to SLR in the Hangzhou Bay (HZB). The scenarios of SLR in the next hundred years count for both non-uniform trends based on historical altimetry data and uniform trends from the latest IPCC projections. In a comparison of model results under different SLR scenarios, we found that the tidal range is ampliﬁed by SLR in HZB with stronger ampliﬁcation at the shallow southern coast. Tidal range change generally increases with the SLR scale; however, neglecting the heterogeneities in the spatial distribution of SLR tends to overestimate the SLR effects. The harmonic analysis illustrates that SLR exaggerates the dominated semidiurnal tides (M 2 and S 2 ) but dampens their overtides and compound tides (M 4 , M 6 , and MS 4 ), of which M 2 amplitude ampliﬁcation explains 71.2–90.0% of tidal range change. SLR tends to promote tidal energy entering HZB through the Zhoushan Archipelago (ZA) compared to the prototype, while dampened sea-bed roughness and reduced tidal velocity come with a less dissipative environment in HZB, resulting in 6–18% more tidal energy exported upstream. Numerical experiments indicate ZA has signiﬁcant effects on tidal responses and energy ﬂux generation, therefore, its quantitative inﬂuences and physical mechanism are also discussed in this paper.


Introduction
Sea level rise (SLR) is indeed one of the most threatening consequences of ongoing global warming [1]. Since the last century, global mean sea level from tide gauges and altimetry observations shows a virtually certain rising trend with a rate of 1.4 mm/a over the period ; however, an accelerating trend has been detected in the last fifty years as the SLR rate increased from 2.1 mm/a over the period 1970-2015 to 3.2 mm/a over the period 1993-2015 [2]. SLR leads these low-lying subsiding areas into the flood-prone coastal environment. For example, many areas of the lagoon from supratidal/intertidal to intertidal/subtidal, and the increase in mean sea level expected for the end of the century, could make the lagoon more vulnerable to the effect of frequent storm surges, harming mostly agricultural areas [3]. Comparisons of sea-level change and shore morphological change along the eastern coast of Bangladesh proved that the local SLR also has led to permanent land losses either by inundation or by erosion with approximately two years lag [4,5]. In recent years, a growing number of studies have highlighted that the long-term tidal amplifications occur alongside the rising sea level, the interaction of which can aggravate the impacts of coastal inundation and nuisance flooding [6][7][8]. In Apalachicola Bay, consideration of the interaction between tides and SLR outweighs static flooding (only SLR) by a factor of Tidal dynamics in the study area are controlled by the west-northwest branch of the tidal waves from the northwest Pacific Ocean. Semi-diurnal tides, dominated by the M 2 component, characterize HZB. Under this tidal regime, the tidal range in the ZA is 1.9-3.3 m. After propagating into the HZB, the tidal waves are seriously distorted upstream due to the width convergence and bed rising landwards. Tidal ranges are amplified from 3-4 m at the mouth to 4-6 m further upstream [25]. The tidal current near open water is usually rotary but is rectilinear inside the archipelago because of the restriction of coastlines [26]. Ebb and flood currents are asymmetric in the HZB, namely, the former dominates in the southern flats, while the latter dominates in the northern region.
The main river systems discharging runoff and sediment into HZB are Qiantang River and Yangtze River. The annual mean discharge of Qiantang River is 952 m 3 /s. Due to the monsoon climate, the river discharge shows a clear seasonal variation: low discharge occurs from August to the following March and high discharge occurs from April to July [27]. The coastal area of HZB is cyclone-prone. Of the 2157 typhoons that occurred in the Western North Pacific basin from 1949 to 2013, 202 affected HZB and its surroundings represent a mean annual frequency of 3.2 [28,29].

Model Configurations
The tidal model adopted in the present study was initially established as the Zhejiang Coast-China Sea nested model, based on the MIKE 21 flexible mesh model. The parent model ( Figure 1b) including a large portion of China coastal waters extends from 16.8 • N to 40.8 • N, and 105.8 • E to 133.8 • E. The bathymetric data, interpolated into the model mesh, was obtained from two sources, including high-resolution nautical charts in China shallow waters and ETOPO-1 data (National Geophysical Data Center) for offshore regions. The open boundary lies in the region where the potential tidal variability induced by SLR is limited (<0.08 m in amplitude and <1 • in phase), according to our previous study [15]. The parent model was driven by the astronomic tidal level derived from the NAOTIDE global model (NAO.99b) [30], which shows better performance in the China sea area [31].
The nested child model−Zhejiang Coast model, covering the entire Yangtze River Estuary (YRE) and Zhejiang Coast, is driven by the simulated tidal elevations extracted from the parent model. A higher resolution mesh was designed in HZB and ZA for the complex coastline and narrow channels. The bathymetric data, interpolated into the mesh of the coastal area, was extracted from a catalog of high-resolution nautical charts and publications. To better describe bed roughness, a spatially varied Manning number for the bottom friction is calculated by an empirical formula of water depth in shallow water. Specifically, the Manning number varies from 10 m 1/3 /s to 200 m 1/3 /s and 60 m 1/3 /s to 116 m 1/3 /s in the China Sea model and Zhejiang Coast model, respectively. Considering the computational expense and the water columns are well-mixed in most of the study area, the model is run in a two-dimensional mode and this simplification will not affect the tidal processes simulation [31].

SLR Trends
Both non-uniform and uniform SLR scenarios are considered in this study. For nonuniform SLR, we firstly obtained historical linear trends in MSL based on altimeter data from AVISO (MSL_Map_MERGED_Global_AVISO_NoGIA_Adjust) and projected it on the model domain ( Figure 2a). Then, this sea level trend was modified to include the influence of the Glacial isostatic adjustment GIA, using the ICE-7G_NA/ VM7 model (Figure 2b) [32]. The rates of relative sea level rise (RSLR) were finally obtained after this adjustment. For the Zhejiang coastal area, this RSLR rate ranges from 1.6 to 4.7 mm/a during 1993-2021, with a spatial mean of 3.3 mm/a ( Figure 2c). The RSLR scale by the end of 2100 is assumed to linearly increase to 0.16-0.47 m (Figure 2d), which was conducted in the non-uniform SLR case.
For uniform SLR scenarios, 0.33 m, the spatial mean of the non-uniform SLR is first considered in the study to assess the impact of SLR spatial distribution on tidal processes. The other SLR projections by the end of 2100 considered are based on the latest released Sixth Assessment Report of Intergovernmental Panel on Climate Change (IPCC) [2]. In detail, SLR projections of 0.56 m and 0.77 m respectively under SSP245 and SSP585 scenarios were formulated, where SSP585 is the scenario with the highest radiative forcing as a consequence of fossil-fuel-based development. These projections were modified to include the regional mean GIA (0.012 m). Finally, the SLR projections of 0.57 m and 0.78 m were imposed in the numerical simulation. This method involves increasing the water depth but not allowing new areas to flood because of the coastal defense works. Simulations with five SLR scenarios were run for 200 days with the first 26 days for spin-up and the last 174 days covering six spring-neap tidal cycles for data analysis. The flowchart of the methodological approach to evaluating SLR effects on tidal dynamics in this study is shown in Figure 3.

Responses of Tidal Range to SLR
The spatial pattern of the mean tidal range, modeled at the present sea level (PSL) scenario is mapped in Figure 4a. According to the results, the YRE and ZA are dominated by the meso-tide with a mean tidal range of 2-3 m, whereas the tidal range at the mouth of HZB is 3-4 m and further enlarged to 6 m at the head due to the width convergence and bed rising landwards. Figure 4b shows that the tidal range along the northern coast is generally larger than that along the southern coast; however, SLR produces more significant tidal range amplification along the southern coast. The largest increase in mean tidal range appears near the Andong tidal flat and the bay head, varying from 0.3 m to 0.8 m in four SLR scenarios. The ratio of tidal range change and the mean tidal range in each grid are calculated and presented in Figure 4c,d. In all the SLR scenarios, the SLR-induced tidal range changes along the northern and southern coasts are less than 10% and 18%, respectively. Dramatic changes mostly happen in shallow water, where SLR produces more inundation and reduces the shoaling effects by increasing the relative water depth, and consequently exaggerates the tidal ranges. Change ratios of tidal range are considered proportional to the depth change ratio increase in Figure 4e. A stronger linear correlation (0.7775) is captured along the southern coast. The spatial patterns of tidal range responses illustrate that tidal ranges in the HZB and ZA generally are enlarged by SLR, except for a slight reduction around the Luchaogang sea area in the RSLR scenario (shown in Figure 5). Besides overwhelming increases along the rim of the Andong tidal flat and the head of HZB, relative high tidal range amplification also occurs in the southeast part of HZB between Andong tidal flat and ZA, which induces the tidal ranges increases from NE to SW, perpendicular to the coastline between the Zhenhai and Andong tidal flat. Among all the scenarios, non-uniform SLR produces the smallest tidal range amplification. The mean tidal range change within the entire HZB is 0.037 m and only 4.8% of the area shows significant tidal range changes above >0.1 m (Figure 5a). Non-uniform SLR has low rates within HZB but high rates offshore. Neglecting the heterogeneity of SLR, its space averaged value overestimates the SLR rate in HZB. Therefore, the tidal range changes in HZB are exaggerated, where mean amplification is enlarged to 0.073 m and the area with tidal range changes above >0.1 m increases to 16.5%.
Larger SLR scales increase the mean amplification of tidal ranges within the entire HZB to 0.094 m and 0.116 m in SSP245 and SSP585 scenarios, respectively. Meanwhile, the 0.1 m isoline moves towards NE, with tidal range changes above >0.1 m increasing to 50.4% and 75.6% (Figure 5c,d).

Tidal Components Behavior
We conducted harmonic analysis on the simulated tidal elevation records at each mesh and drew the co-tidal charts for the four major tidal constituents (K 1 , O 1 , M 2 , and S 2 ) and three shallow-water components (M 4 , MS 4 , and M 6 ) in Figure 6a. The amplitudes demonstrate the semidiurnal tides dominate in this area, of which the M 2 tide has one order higher amplitudes than other tides. Specifically, the spatial averaged amplitudes of M 2 , S 2 , K 1 , and O 1 tides in the HZB-ZA system are 1.54 m, 0.42 m, 0.36 m, and 0.21 m, respectively, whereas the amplitudes of M 4 , MS 4 , and M 6 are only about 7.6%, 4.5%, and 3.2% of M 2 amplitude, respectively. SLR enlarges the water depth and accelerates wave propagation, resulting in an earlier tidal arrival. The phases of four major tidal constituents decrease in the domain and these reductions are enhanced from mouth to head of the HZB (Figure 6b). Comparisons of their spatial averaged phase changes in different SLR scenarios (Figure 7a) show the phase reductions increase with the SLR scale, however, the mean phase changes (also the spatial distribution, which is not shown) are comparable in RSLR and 0.33 m SLR scenarios. This is expected as follows: neglecting the bottom friction, the speed of tidal wave c can be approximated by c = gH. Because the tidal phases of major tidal constituents are controlled by far-field tidal wave systems in ECS and the Yellow Sea (YS) [15], regional characters of the water depth H as well as the SLR count for their phases changes rather than the local variations. Therefore, the offshore non-uniform SLR and its regional averaged uniform SLR are considered to result in similar tidal wave propagation acceleration and phase modification. For shallow water components, amphidromic points degenerate along the southern coast of HZB ( Figure 6a). The movement and anticlockwise rotation of these degenerate systems under the SLR effect produce complex phase decline (Figure 6b). Maximum phase reduction occurs near the amphidromic point, but the spatial mean is not closely related to the SLR scale ( Figure 7a).  SLR tends to enlarge the tidal amplitudes of semidiurnal tides rather than other constituents. Figure 6c shows the iso-lines of M 2 amplitude change increase from NE to SW with more amplitude increase on the southern coast, while the S 2 amplitude increases have an insignificant difference between two coasts. Although the amplitudes of shallow water components are generally smaller than those of diurnal tides, their changes are more significant. The movements of shallow water tidal systems trigger a larger amplitude reduction along the northern coast, while having little change around the Zhenhai sea area. Quantitative comparisons of mean amplitude changes in different SLR scenarios ( Figure 7b) illustrate that SLR-induced amplitude changes of M 2 and S 2 are similar, the changing ratios referring to their amplitudes of 1.0%, 1.9%, 3.0%, and 4.0% in RSLR, 0.33 m SLR, SSP245 and SSP585 scenarios, respectively, of which M 2 amplitude amplification explains 71.2-90.0% of tidal range increases in different SLR scenarios. The amplitude of M 2 and S 2 tides is more sensitive to the local depth changes. Because 0.33 m SLR diminishes more bottom fraction in HZB compared with non-uniform SLR, contributing to the tidal amplitude amplification, therefore M 2 and S 2 amplification in the RSLR scenario is only half of that in the 0.33 m SLR scenario. While the mean amplitude reductions of M 4 and MS 4 are comparable in RSLR and 0.33 m scenarios, which implies the movement of the tidal wave system is the trigger for tidal amplitude change for shallow water components within the HZB.

Tidal Energy Flux and Dissipation
The ZA is regarded as an ideal area where the tidal energy can be explored for power generation [24]; and the incoming tidal energy is a convergence at the head of the funnelshaped HZB and QE, resulting in the formation of the world-famous powerful tidal bore. Therefore, we present a detailed analysis of the tidal energy budget and the responses of tidal energy flux and dissipation to SLR in this area (details of the calculation refer to the Appendices B.1 and B.2).   Figure 9 illustrates the residual tidal energy flux compared to the prototype (Figure 8a) in vector and magnitude. It can be seen that the tidal energy flux generally increases with SLR in the study domain. Some of the residual energy flux propagating from the ECS flow toward the north, and some enter HZB through ZA, resulting in more than 10 KW/m of energy flux increase in the channels. After entering the HZB, the incoming flux is divided into two branches: one branch flows anti-clockwise toward the upper part of HZB and the other flows clockwise out of the bay through the northern mouth. As a consequence, the former branch enhances the energy flux around the head, leading to 6-18% more tidal energy transport into the QE through S5 (Table 1), while the latter diminishes the incoming flux and induces 12.6-50.1 MW energy declines through S4. Comparisons of energy flux trough S1 and S5 indicate that the residual energy is in proportion to the SLR scale; however, for S4, the non-uniform SLR tends to diminish more incoming flux than the 0.33 m SLR, which may be caused by the smaller local SLR and tidal range reduction (Figure 5a).  Bottom dissipation is a product of the bottom friction coefficient and the cub of the bottom tidal velocity, so the SLR-induced bottom dissipation change is expected to be more complicated. On one hand, the depth-related bottom friction coefficient becomes smaller, thus the tidal process tends to be less frictional. On the other hand, the SLR enhancing/diminishing tidal velocity can result in more/less bottom dissipation. In all SLR scenarios (Figure 10a-d), it can be seen that SLR reduces bottom dissipation in most areas of the shallow bay. This bottom dissipation reduction results from both diminished sea-bed roughness and reduced tidal velocity (Figure 10e-l). Significant bottom dissipation increase is found around the Andong tidal flat and Nanhui headland, which is dominated by a longer inundated period of these intertidal zones and the exaggerated tidal velocity induced by SLR. Similarly, the influence of increased velocity on the bottom dissipation overwhelms bottom friction coefficient effects in the ZA, so the bottom dissipation in this domain is gradually increased as the SLR becomes higher. From Table 1, it is evident that SLR damps influx and promotes outflux of the tidal energy, resulting in less energy dissipation within the HZB. It is notable that although the changes of dissipation are captured in most of the study areas, the changes in the tidal energy balance induced by the SLR are minimal, which is no more than 5.2% in the SSP585 scenario.

Role of ZA on the Tidal Modification
The tidal ranges and dominant M 2 amplitude are exaggerated along the southern coast compared to the tides along the northern coast. It indicates that ZA is a crucial geomorphic element that affects the tidal responses in the HZB. Therefore, we rerun the models eliminating the multiple islands under PSL and SSP585 scenarios to distinguish the moderating effects of ZA between tidal dynamics and SLR.
Semidiurnal tides and their shallow-water components in HZB are highly choked by constricted channels of ZA ( Figure 11). M 2 is the dominant tidal constituent, whose amplitude is sharply decreased when passing through ZA, and a 20% amplitude reduction is produced in HZB referring to the ZA eliminated reruns. In addition to amplitude attenuation, up to 22 min (11 • ) of a phase lag also occurred. This tidal choking effect is more significant in the higher frequency overtides, with a 32% larger amplitude attenuation rate and 36 min (36 • ) phase lag captured in M 4 . It is notable that though ZA produces amplitude attenuation and a phase lag, the imbalance of M 2 amplitude between two coasts remains unchanged. A faster and stronger M 2 tide, as well as its faster but smaller overtide M 4 , is also captured in the HZB under the SLR effect in the ZA eliminated cases (Figure 12). M 2 amplitude amplification is comparable at two side coasts in the ZA excluded case; however, ZA induces the asymmetry of amplitude amplification. It is because SLR comes with larger flow cross-sections of archipelago channels, which reduces the tidal chocking intensity [18]. Therefore, more M 2 tides are released, which enlarges the tidal amplitude on the southern coast of the bay.  Based on these numerical experiments, we infer that the ZA significantly dampens the mean tidal ranges of the study domain from 3.9 m to 3.2 m by the tidal choking effect. Though ZA plays a minor role in the formation of M 2 amplitude asymmetry between two coasts, it superposed on SLR effects, which produces a smaller tidal choking intensity and tends to narrow the tidal range gap between two coasts.

Standing Environment and Tidal Flux Variation
The model simulation showed that the deep channels between ZA are crucial waterways for the tidal flux, of which approximately 36% of tidal energy flux is transported into the HZB through them. This is consistent with previous findings of Wu et al. [24]. Results suggest that tidal ranges, tidal velocity, and tidal flux passing the archipelago are increased by the rising sea level ( Figure 9); however, the SLR has little effect on the incoming tidal energy of the HZB ( Table 1). The progressive character of the tidal wave is regarded as a useful local estimate of tidal energy dynamics since the Stokes transport is greatest for progressive waves and negligible for standing waves. The tides in the HZB are semidiurnal-dominated. Following Holleman and Stacey [33], we use the tidal phase ϕ c , (refer to the Appendix B.3 for calculation) to diagnose the standing wave or progressive wave dynamics. ϕ c = 90 • means peak flood velocity leads the peak sea surface elevation by 90 • , i.e., standing wave. ϕ c = 0 • means the free surface fluctuation and the tidal velocity are in the phase, i.e., progressive wave.
The progressive wave from the ESC sea becomes more progressive after entering the ZA with an enhanced tidal energy flux (Figure 13a). However, as the M 2 tide out of the archipelago, it sharply shifts toward a more standing wave (ϕ c > 60 • ) and develops into quasi-standing near the southern coast where vast tidal flats exist. In terms of the standing tidal wave in the bay, the energy flux is grossly damped throughout the archipelago and most of the energy is transported along the northern coast where the tidal wave is less standing. SLR shifts the mouth of HZB toward a more reflective, standing wave environment. It is apparent in Figure 13b that more standing wave blocks the residual incoming tidal energy flux flowing into the inner bay. In contrast, SLR damps the standing wave at the head of bay, leading to more exportation of the energy flux.
HZB is close to a standing wave resonance environment, with a mean velocity phase lead for the M 2 constituent of approximately ϕ c = 61.5 • . This resonance feature of the inner bay is dominated by the converging geometry and hardly affected by ZA. Even by eliminating all the islands, the mean velocity phase lead for the M 2 constituent remains at 61.3 • (Figure 13c). The SLR similarly produces a more standing/progressive wave environment at the mouth/head of HZB and the spatial changes of ϕ c and the tidal energy flux (Figure 13d) are comparable to Figure 13b. It infers that the inner bay resonance feature is less dependent on the archipelago scattered at the outer bay. Resonance is easily known as the superposition of an incident wave and a reflected wave. In a non-frictional basin where its length is near a quarter or three-quarters of wavelength, standing tidal waves occur as the coming wave is fully reflected [15,33]. The mean depth of HZB is 10.5 m at PSL, and the three-quarters of M 2 wavelength is 94.5 km. Increasing SLR moves threequarters of M 2 wavelength closer to the basin length (from Zhapu to Nanhui headland) of 118 km, thus causing a stronger resonance or standing wave. Notably, the HZB is not an idealized non-frictional basin. Friction is crucial in the upper branch of shallow water which dissipates incoming and reflected waves. SLR dampens the bottom dissipation as shown in Figure 10, which may result in the heterogeneity of ϕ c changes in Figure 13b,d.

Physical Mechanism for the Enhance of Progressive Waves around ZA
To gain insight into the mechanism for the enhancement of progressive waves around ZA, here we follow Hench and Luettich [34] to estimate the contribution of each term in the x-and y-momentum equations in ZA included and excluded runs. We choose five representative stations (S1-S5, Figure 1a) for analysis. Among them, S1-S4 help to describe the tidal dynamics around islands, while S5 helps to demonstrate the tidal dynamics of the inner bay. Nonlinear term variations in the momentum balance equations during one tidal cycle are shown in Figures 14 and 15, and their mean absolute values are listed in Tables 2 and 3. For S1-S4 where the progressive wave prevails, the momentum balance is primarily between advection and barotropic pressure gradient ( Figure 14, Table 2). This means once the strong currents form in the archipelago due to complex coastline and narrow channels between islands, water level tends to drop to balance the momentum, suggesting the momentum balance near the islands follows Bernoulli law [24,35]. While for the S5 in the shallow bay where standing wave prevails, momentum balance is mainly between local acceleration and barotropic pressure gradient, meanwhile bottom friction makes the third contribution. These momentum balance differences between S1-S4 and S5 imply that the Bernoulli-type balance is the substantial mechanism for the formation of local quasi-progressive waves in the ZA. Because the essence of the Bernoulli law is the transformation between tidal potential energy and tidal kinetic energy, the elevation drop and current velocity rise mean the tidal potential energy is transferred to tidal kinetic energy. Confined by energy conservation, the elevation drop and current velocity rise are always in phase, thus progressive environment forms in the archipelago. When the multiple islands are excluded, advection is not crucial in S1-S4, where the momentum balance is primarily between local acceleration and barotropic pressure gradient following linear wave dynamics. Meanwhile, the tidal phase ϕ c is increased by 17-30 • , so the tidal wave shifts more standing. However, the tidal phase ϕ c and momentum balance of S5 are comparable in Tables 2 and 3. These numerical experiments confirm that the archipelago is conducive to the formation of the local progressive wave by the Bernoulli-type momentum balance, while the standing environment in HZB is impervious to the ZA.     Table 3. Same as Table 2, but for the rerun of the models eliminating ZA.

Comparison with Previous Studies
Compared with previous numerical studies in HZB and ZA [24,36], our results show similar patterns of tidal properties under the present-day condition and SLR scenarios. Studies on a continental scale [15,37] have captured a maximum reduction of M 2 amplitude around the offshore YRE, however, differences in tidal response occur in the HZB, which could be due to the low resolution of bathymetry and coastlines in those studies. Regional-scale research with high resolution of bathymetry and coastlines [36] presents tidal amplification and tidal velocity reduction in the HZB under SLR. Similar to his results, we found the largest increase in mean tidal range appears near the Andong tidal flats and the bay head. Our results also confirmed the increasing trend of tidal amplification from the mouth to the head of HZB. However, in the study of Gu [36], the location of maximum reduction of the tidal range moves down towards the mouth of HZB, which induces a smaller area of the tidal amplification inner HZB. Different from numerical studies [15,36,38] that focus on the alarming SLR scales, we gain insight into the tidal response with plausible future SLR. In addition, the effects of non-uniform SLR, which are ignored in other studies, are offset in this study.
The covarying nature of the MSL increase and tidal amplification captured in our study is in line with the observation [22,39,40]. Based on our results, the mean rising rate of the tidal range is only 0.37-1.16 mm/a according to the different SLR rates, even the largest increasing rate on the Andong tidal flat does not exceed the SLR rate. However, the historical data presented an exaggerated rising rate of 0.85-1.48 cm/a in the tidal range covaried with a moderate rising rate of 3.2-5.2 mm/a in MSL. This difference provides insight that the historical rising in MSL played a minor role in tidal amplification. Anthropogenic activity, and other natural factors (i.e., storm events and discharge) may be the major contribution [22]. Especially, HZB is one of the coastal hotspot regions in China, with 1023.86 km 2 of coastal area being reclaimed during 1979-2014 [41]. These large-scale reclamation declined the tidal prism but aggravated the tidal amplification. Li et al. [42], using the numerical model, reviewed the tidal dynamics of HZB from 1974 to 2016 and confirmed a strong interplay between the rapid shrunk bay and tides with a maximum tidal range increased by >2 m. Since the year 2017, reclamation projects have been severely restricted; consequently, tidal amplification by artificial structures will gradually decline. In this view, our model results based on the current coastline are expected to provide more implications for the future tidal dynamics in HZB.
The seasonal or annual variation in runoff and temperature can induce short-term fluctuation in the water level, as well to the tidal dynamics [4,5]. However, the river discharge contribution is generally insignificant compared with tides, even in some riverdominated estuaries or coasts. For example, the daily-averaged and seasonal-cumulative river discharge contributes to 13% of total sea level variance in Kitimat Fjord System, Canada [43]. While in YRE, long-term variation of river discharge explained 6.8-8.9% of the total increase in the sea level. Two main rivers, the Yangtze River and the Qiantang River, impact the hydrodynamic characteristics of HZB. Kuang et al. [44] quantified that a 25% frequency of the Yangtze River discharge in the wet season could produce only 0.005-0.010 m SLR in most of HZB. Xie et al. [45] indicated though that tidal ranges in QE were amplified and damped in the wet and dry seasons, respectively; the high and low water levels in HZB were hardly influenced by the Qiantang River discharge because of the larger width and depth. These studies present that seasonal discharge variation plays a minor role in the sea level of HZB. In addition, our study focused on the tidal dynamic response to the SLR on a centennial time scale, therefore the seasonal perturbations can be filtered in this study. A typhoon event is a crucial factor in the rising sea level. Since the 1980s, the number and the intensity of typhoons that affected HZB and its surrounding represented an increasing trend [22]. However, a more detailed investigation of long-term typhoon effects on tidal dynamics is left for a future study.

Conclusions
Utilizing a nested hydrodynamic model with a variable future SLR, we provided an integrated investigation of tidal responses to SLR in HZB. The scenarios of SLR in the next hundred years count for the non-uniform trend based on historical altimetry data and uniform trends from the latest IPCC projections. In a comparison of different SLR trends imposed on the model, we found the evolution in tides on the regional scale is similar but different in magnitude.
The tidal range is amplified by SLR in HZB, but the changes along the northern and southern coasts are less than 10% and 18%, respectively. Dramatic changes mostly happen in shallow water, i.e., the Andong tidal flats and the bay head, where SLR produces more inundation and reduces the shoaling effects by increasing the relative water depth, and consequently exaggerating the tidal ranges by 0.3-0.8 m. Future SLR also can reduce the tidal choking intensity of ZA, contributing to the tidal ranges amplification at the southern coast of HZB, which induces the tidal ranges increases from NE to SW, perpendicular to the coastline between the Zhenhai and Andong tidal flat. Tidal range change generally increases with the SLR scale, whereas neglecting the heterogeneities in the spatial distribution of SLR, the simplified uniform 0.33 m SLR tends to overestimate the SLR effects.
Harmonic analysis results indicate SLR enlarges the water depth and accelerates wave propagation, resulting in an earlier tidal arrival. The phase reductions generally increase with the SLR scale. Because the tidal phases of major tidal constituents are controlled by farfield tidal wave systems, the offshore non-uniform SLR and its regional averaged uniform SLR are considered to result in a similar phase modification. For shallow water components, maximum phase reduction occurs near the amphidromic point, but the spatial mean is not closely related to the SLR scale. SLR tends to enlarge the tidal amplitudes of semidiurnal tides rather than other constituents. The changing ratios refer to their amplitudes of 1.0%, 1.9%, 3.0%, and 4.0% in RSLR, 0.33 m SLR, SSP245, and SSP585 scenarios, respectively, of which M 2 amplitude amplification explains 71.2-90.0% of tidal range increases in different SLR scenarios. Although the amplitudes of shallow water components are smaller than those of diurnal tides, their changes are more significant. The movements of shallow water tidal systems trigger a larger amplitude reduction along the northern coast, while having little change around the Zhenhai sea area.
The barotropic tidal energy budget in the HZB is also evaluated. The total incoming energy is 4974.0 MW, with 944.2 MW flowing out of the domain and 3805.8 MW dissipated due to bottom friction. SLR tends to promote more tidal energy entering HZB through the archipelago channels compared to the prototype. Then this residual tidal energy flux is divided into two branches. One branch flows anti-clockwise toward the bay head and the other flows clockwise outward the bay. Diminished sea-bed roughness and reduced tidal velocity come a less dissipative environment. Because of this energy dissipation reduction, the residual energy flux enhances towards the head of the bay, with 6-18% more tidal energy transported into the QE. This change in the tidal energy influx can also be explained by the regional reflection effect. Specifically, SLR shifts the mouth of HZB toward a more reflective, standing wave environment, which blocks the residual incoming tidal energy flux flowing into the inner bay. In contrast, SLR damps the standing wave at the head of the bay, leading to more exportation of the energy flux. A detailed momentum balance analysis suggests the archipelago is conducive to the formation of the local progressive waves by the Bernoulli-type momentum balance and harvesting substantial tidal energy; however, the inner bay resonance feature is less dependent on the archipelago scattered at the outer bay.   value less than 0.2.
Comparisons of the simulated and measured tidal processes at representative stations show a good agreement in both magnitude and phase, with all the Skill values of tidal levels exceeding 0.98 ( Figure A1). The Skill scores (0.804-0.895) for velocity are slightly less than those for tidal level ( Figure A2), however, our model has captured the second-order of current behavior during the hindcast periods and the model performance can be regarded as excellent. Therefore, this robust model gives us confidence in the SLR effects forecast.