Analysis and Exploration of the Impact of Average Sea Level Change on Navigational Safety in Ports

Abstract

The primary clientele of a harbor is vessels, and vessels are primarily influenced by external forces such as wind (on the water surface), currents (underwater), and waves (affecting vessel stability). Therefore, it is necessary to comprehensively consider safety factors such as marine environmental forces and port characteristics. As ship sailing falls under applied science, acquiring marine meteorological information regarding ship routes can enhance port navigational safety. However, in the face of changes in the environmental conditions of harbor waters, it is essential to fully consider the impact of the external environment on ship maneuvering. One can effectively navigate complex operating environments by devising reasonable ship-handling plans. In the context of sea level rise caused by extreme climatic events, long-term variations, trends, and random factors are at play. Previous assessments of sea level rise have often relied on linear regression and the least squares method to determine coefficients. However, these methods fail to accurately capture the actual trend of sea level rise. Additionally, traditional harmonic analysis methods are unable to analyze sea level rise as well. Therefore, in this study, the techniques of simple moving average (SMA), empirical mode decomposition (EMD), and ensemble empirical mode decomposition (EEMD) were applied to analyze sea level rise. The obtained results of sea level rise under different analysis conditions were integrated with a hydrodynamic model that incorporates both wave and tidal characteristics to calculate the overall coastal dynamics parameters, which are crucial for ship navigation. The research findings contribute to the study of ship navigational safety issues by examining the distribution characteristics of port meteorology under climate change conditions. They offer valuable insights for mariners to assess navigational safety and devise maneuvering strategies based on the actual water flow conditions. Furthermore, the findings help identify and address potential risks and issues, ultimately ensuring the safety of navigation.

1. Introduction

Ports are the foundation for maritime engineering development and are the hub for sea and land transportation. More than 90% of global trade is conducted through maritime shipping. The effectiveness of port construction relies on ensuring the safe navigation of ships during entry, safe navigation within port waters, and secure mooring and loading/unloading operations. Over time, ship designs have evolved toward larger sizes to meet the advancements in shipbuilding technology and the demands of maritime transportation. The water area between the outer channel and the breakwater entrance of the harbor is primarily influenced by environmental factors such as wind and currents. Safety considerations in this area require a comprehensive assessment of ship handling, environmental forces, and port characteristics. Evaluating the safety of entering the port relies on considering these factors. Therefore, the requirements of a port are closely related to climate conditions.

Climate change issues related to global warming have become increasingly prominent in recent years. One of the critical concerns is the rise in sea levels, along with other environmental changes. Research on sea level rise [1,2,3,4,5,6] indicates that the global mean sea level has been rising at approximately 1.5 to 2.4 mm per year during the 20th century. According to the final draft of the Working Group I contribution to the Intergovernmental Panel on Climate Change (IPCC) Sixth Assessment Report (AR6) released on 9 August 2021 [7], observations of extreme weather and climate events since 1950 have shown an increase in the frequency and duration of warm periods, rainfall intensity, intense tropical cyclone activity, and the likelihood of extreme high sea levels. The report also states that from 1901 to 2018, the global mean sea level has risen by 200 mm, with a rate of 1.3 mm per year from 1901 to 1971, 1.9 mm per year from 1971 to 2006, and an increase of 3.7 mm per year from 2006 to 2018. Climate change has significantly impacted most terrestrial and marine areas worldwide. Changes in hydrological conditions, for example, directly affect the characteristics of coastal environments and can even impact the functionality of coastal protection structures and port operations. Due to global climate change, there has been a rise in average sea levels, leading to changes in local hydrodynamic features such as tides, waves, and currents in coastal areas. Considering the safety of vessel navigation, it is primarily influenced by environmental factors such as wind and currents. It is essential to effectively understand the dynamics of tides, waves, and currents in the vicinity of ports and the distribution characteristics of nearshore flow fields influenced by the water depth and topography near the ports. Therefore, studying tidal and sea level variations holds significant importance and practical value [8,9,10,11,12,13,14].

Keelung Port is an important harbor in northern Taiwan. During the era of economic takeoff, it played a crucial role as the country’s artery for international trade and cargo transportation. Container shipping is the primary mode of transportation, supplemented by bulk cargo. In December 2009, the largest cargo vessel to ever dock in the port, the MSC MARIA ELENA, arrived. The MSC MARIA ELENA had a capacity of 9100 TEUs (TEU: twenty-foot equivalent unit), with a length of 336 m, a width of 45.6 m, and a maximum draft of 14 m. In recent years, with the government actively promoting the transformation of Keelung Port into a tourism and recreational harbor, it has become Taiwan’s primary international cruise port. In May 2023, the port welcomed the MSC Bellissima cruise ship, with a total tonnage of 172,000 tons (length: 315.83 m, width: 43 m, draft: 8.75 m), setting a record as the largest cruise ship to dock at Keelung Port. However, due to their significant inertia and small length-to-width ratio, large vessels require the careful consideration of ship-handling characteristics, environmental forces, and port characteristics when maneuvering to assess their safety during port entry [15,16,17,18,19,20,21,22,23].

Due to the influence of coastal dynamics factors such as tides, waves, currents, and topography, the distribution of water flow near port channels is affected. Considering the rising sea levels caused by global warming, there is concern that it may alter the characteristics of tidal patterns, waves, and currents in the marine environment. Among these external conditions, the selection of basic flow field parameters is crucial for vessel navigation, as the flow distribution can be affected by changes in port structures and topography, resulting in more complex flow conditions in port areas [24,25,26,27]. In general, when the direction of water flow is not aligned with the ship’s heading, the flowing water exerts a hydrodynamic pressure on the vessel. The impact on the ship becomes more significant as the angle between the flow direction and the ship’s heading increases.

Furthermore, as the ship moves, the effects of the water flow on the vessel shift gradually from the bow to the stern [28,29,30]. To ensure navigation safety, the assessment of the safety impact on ship navigation during port entry and exit is often conducted using computer simulations. Therefore, this study employed a harmonization analysis framework to investigate coastal tidal characteristics and sea level trends. The analysis involved the application of simple moving average (SMA), empirical mode decomposition (EMD), and ensemble empirical mode decomposition (EEMD) to analyze sea level rise. Next, the obtained results of a sea-level-rise analysis will be validated using a numerical model considering wave–flow interactions, with field observations as background data. This study investigates the impact of port sea–meteorological distribution characteristics on ship navigation safety under climate change conditions. The analysis results will provide navigation personnel with a safety assessment and maneuvering strategy considering the actual distribution of water currents. Additionally, it can serve as a reference for maritime and relevant authorities in ensuring navigation safety in the fairway.

2. Sea-Level-Rise Analysis

2.1. Analysis of Water Level Data

This study selected the location of the tide gauge station at Keelung Port, as shown in Figure 1 (25.155° N, 121.7522° E). However, the instrument’s position may change due to human factors such as maintenance or natural disasters like typhoons, which can result in instrument damage or the cessation of recording. These factors can affect the integrity of the actual tide data. The statistical analysis of hourly water level data, as shown in Figure 2, indicates that the completeness of the tidal data at the Keelung Port tide gauge station from 2005 to 2022 is approximately 88.9%. To compensate for the data gaps, the measured data were subjected to harmonic analysis, and the missing astronomical tide data were filled in at the gaps in the original data. The missing tidal data in the original dataset were filled with predicted astronomical tide data, as shown in Figure 3. Subsequently, three different methods, namely SMA, EMD, and EEMD, were applied separately to analyze the sea level rise at Keelung Port.

2.2. Methods for Analyzing Sea Level

In this study, the harmonic analysis framework was utilized to analyze sea level rise using three different methods: SMA, EMD, and EEMD. Moving average (MA) aims to smooth out data fluctuations by taking the average, making the trend of the data more apparent. MAs can be categorized into three methods: simple moving average (SMA), weighted moving average (WMA), and exponential moving average (EMA), depending on the calculation method. WMA and EMA assign weights to data from different periods based on the SMA. The weights are higher for data closer to the current period, indicating a more significant impact on the average. In simple terms, WMA and EMA are more responsive to recent data changes. However, the SMA is still the most commonly used method, especially in various technical indicators and analyses. The formula for calculating the SMA is as follows:

$$F_{t+1}=\frac{D_t+D_{t-1}+D_{t-2}+\dots +D_{t-n+1}}{n}$$

(1)

where $F_{t+1}$ represents the forecast value for the next period, $n$ denotes the number of periods for the moving average, and $D_t$, $D_{t-1}$, $D_{t-2}$, and $D_{t-n+1}$ represent the data for the previous period, two periods ago, three periods ago, and up to $n$ periods ago, respectively.

This study applied the SMA to analyze the tidal water level data from the Keelung tide gauge station from 2005 to 2022. The moving average was calculated over one year, resulting in the Keelung tidal water level trend line shown in Figure 4. The estimated rise in sea level is approximately 17.68 mm, with an average sea-level-rise rate of 0.98 mm/yr.

Next, this study applies the Hilbert–Huang transform (HHT), which consists of empirical mode decomposition (EMD) and Hilbert transform (HT), for analysis. The HHT can analyze non-stationary and non-linear data, making it applicable in various fields such as seismic signal analysis, the detection of damage in bridges or highway structures, the processing of biomedical signals, and the measurement of tides or waves. The common characteristic of these applications is that they can be presented in frequency form, and the input signals are often non-stationary data [31]. HHT decomposes a signal into multiple components called intrinsic mode functions (IMFs). These IMF components are not necessarily complete and symmetrical sinusoids. Their periods and amplitudes are not fixed, allowing for meaningful instantaneous frequency and instantaneous amplitude to be directly obtained using the Hilbert transform (HT). HHT exhibits physical significance, high precision, and adaptability. The empirical mode decomposition (EMD) in HHT refers to decomposing a signal into multiple IMF components. This process requires multiple iterations to obtain each IMF. The procedure is as follows.

The original signal $s\left(t\right)$ is inputted, and all local maxima and local minima are identified. Then, curve fitting is performed on all the local maxima and local minima. Cubic spline interpolation is commonly used for curve fitting in HHT. Cubic spline interpolation provides results closest to the actual curve compared to other interpolation methods. After the interpolation is completed, upper and lower envelopes are obtained separately. Then, the average of the upper and lower envelopes is calculated to obtain the mean envelope $m_1\left(t\right)$. The original signal is then subtracted by the mean envelope $m_1\left(t\right)$, resulting in the first component $h_1\left(t\right)$, as shown in Equation (2):

$$h_1\left(t\right)=s\left(t\right)-m_1\left(t\right)$$

(2)

After identifying the first component $h_1\left(t\right)$, it is necessary to check if it meets the conditions for being an IMF component. There are two conditions for an IMF component. First, the sum of the total number of local maxima and local minima should equal the number of zero crossings or differ by one. Second, the average of the upper and lower envelopes should approach zero at any given time. If these conditions are not met, further filtering is required. The component $h_1\left(t\right)$ is treated as the original signal, and the steps of interpolation and obtaining the mean envelope are repeated. The relationship can be expressed as Equation (3):

$$h_2\left(t\right)=h_1\left(t\right)-m_2\left(t\right)$$

(3)

Repeat the above process $k$ times, as described in Equation (4), until the conditions for an IMF component are met:

$$h_k\left(t\right)=h_{k-1}\left(t\right)-m_k\left(t\right)$$

(4)

Once $h_k\left(t\right)$ satisfies the conditions for an IMF component, the first IMF component $c_1\left(t\right)$ can be obtained, as shown in Equation (5):

$$c_1\left(t\right)=h_k\left(t\right)$$

(5)

The residual signal $r_1\left(t\right)$ can be obtained by subtracting $c_1\left(t\right)$ from the original signal $s\left(t\right)$, as shown in Equation (6):

$$r_1\left(t\right)=s\left(t\right)-c_1\left(t\right)$$

(6)

Next, the process is repeated with $r_1\left(t\right)$ as the new original signal, resulting in sequentially obtaining new residual signals $r_2\left(t\right)$, $r_3\left(t\right)$,..., $r_n\left(t\right)$, as shown in Equations (7)–(9):

$$r_2\left(t\right)=r_1\left(t\right)-c_2\left(t\right)$$

(7)

$$r_3\left(t\right)=r_2\left(t\right)-c_3\left(t\right)$$

(8)

$$r_n\left(t\right)=r_{n-1}\left(t\right)-c_n\left(t\right)$$

(9)

When $r_n\left(t\right)$ becomes a constant or a monotonic signal, meaning it cannot be further decomposed into IMF components, the decomposition process of EMD is considered complete. Finally, the original signal can be represented as the sum of all IMF components and the average trend component $r_n\left(t\right)$, as shown in Equation (10):

$$s\left(t\right)={\displaystyle \sum _{k=1}^n{c}_k\left(t\right)+r_n\left(t\right)}$$

(10)

After applying EMD to the original signal, we obtain various IMF components. By using the HT, we can then obtain the Hilbert spectrum, which is represented by Equation (11):

$$X\left(t\right)=a_j\left(t\right)e^{i2\pi {\displaystyle \int f_j\left(t\right)dt}}$$

(11)

By defining time boundaries, we can obtain the marginal spectrum, which is expressed by Equation (12):

$$h\left(\omega \right)={\displaystyle {\int }_0^{T}H\left(\omega ,t\right)dt}$$

(12)

The entire process of HHT is illustrated in Figure 5.

However, the conventional EMD method suffers from mode mixing, where different scales (i.e., different amplitudes) can be observed within a single IMF component, or the same scale (i.e., the same amplitude) can occur in different IMF components. Mode mixing leads to the loss of physical meaning in IMF components, affecting subsequent analysis results. To address this issue, ref. [32] proposed an improved EEMD method. In EEMD, the original signal is combined with white noise, and the relationship between the random noise and the signal is given by Equation (13):

$${\epsilon }_n=\frac{\epsilon }{\sqrt{N}}$$

(13)

where $N$ represents the total number of ensembles, $\epsilon$ denotes the amplitude of the added white noise to the original signal, and ${\epsilon }_n$ represents the standard deviation of the absolute error. White noise is a signal with small amplitudes that behaves randomly, and its energy is uniformly distributed. It is recommended to use 0.2 times the original data’s standard deviation as the noise’s amplitude [32]. Ultimately, the added white noise can be eliminated, ensuring that the original signal is not affected by its presence.

In this study, the characteristics of water level data were analyzed using the EMD method developed in [31,33]. After applying EMD to analyze the water level data from 2005 to 2022, the water level variation trend is shown in Figure 6. The estimated rise in sea level is approximately 51.20 mm, with an average sea-level-rise rate of 2.84 mm/yr. The EEMD method developed in [32] was also employed to analyze water level data. The standard deviation of the original data was 0.3, and the amplitude of the white noise was set to 0.06. The water level variation trend obtained from this analysis is illustrated in Figure 7. The estimated rise in sea level is approximately 57.73 mm, with an average sea-level-rise rate of 3.21 mm/yr. It is observed that the water level shows an initial rise followed by a decline. This trend may be attributed to the deviation caused by the end-point effect resulting from missing data. Therefore, it is not possible to infer a continued decline in the water level trend at the Keelung tide gauge station. The results of analyzing the sea level rise using the SMA, EMD, and EEMD methods for 2005 to 2022 at the Keelung tide gauge station are summarized in

Table 1. Analysis results of sea level rise.

Analysis Method	Avg. SLR	Avg. SLR Rate
SMA	17.68 mm	0.98 mm/yr
EMD	51.20 mm	2.84 mm/yr
EEMD	57.73 mm	3.21 mm/yr

. The analysis results indicate that the Keelung tide gauge station’s average sea-level-rise rate ranges from 0.98 to 3.21 mm/yr.

According to the IPCC AR6 report, the global average sea-level-rise rate from 2006 to 2018 was 3.7 mm/yr. Based on the analysis of other researchers’ studies [34,35,36,37,38] in Taiwan, from 1975 to 2012, the average sea-level-rise rate ranged from 2 to 7 mm/yr. Comparing the global average sea-level-rise rate from the research report with the analyzed average sea-level-rise rate at the Keelung tide station can help assess the trend of sea level rise in the port of Keelung with the global trend, thereby increasing the credibility and significance of the study. After obtaining the sea-level-rise results, numerical models considering wave–current interactions were used, along with field observations as background data for model validation. This was performed to investigate the distribution characteristics of hydrodynamic environments near the port of Keelung under different sea-level-rise scenarios.

3. Hydrodynamic Characteristic Analysis

3.1. Model Setup

The present study aims to understand the impact of hydrodynamic characteristics on ship navigation in Keelung Harbor under different sea-level-rise scenarios. A hydrodynamic model that incorporates wave and tidal current characteristics is used for the calculations. The numerical model first analyzes the overall coastal dynamics and adjusts specific numerical parameters and input/output data formats and content for areas that may be affected. Subsequently, the model analyzes the influence of sea level rise on the wave and flow fields [39,40,41]. The calculation of the planar wave field in the hydrodynamic model utilizes the unsteady-type governing equations for water waves, including the effects of currents on inclined slopes [42]:

$$\frac{D^2\phi }{D{t}^2}+\left(\nabla \cdot \stackrel{\rightharpoonup }{U}\right)\frac{D\phi }{Dt}-\nabla \cdot \left(C{C}_g\nabla \phi \right)+\left({\sigma }^2-k^{2}C{C}_g\right)\phi =0$$

(14)

Under the assumption of non-rotational, single-frequency linear surface waves, the wave potential function can be expressed as follows:

$$\phi \left(\stackrel{\rightharpoonup }{x},\stackrel{\rightharpoonup }{y},z,t\right)=f\left(a,h\right)\phi \left(\stackrel{\rightharpoonup }{x},\stackrel{\rightharpoonup }{y},t\right)$$

(15)

where $f\left(z,h\right)=\frac{\cosh \left[k\left(h+z\right)\right]}{\cosh kh}$. For the harmonic motion of waves on a single periodic wave, Equation (15) can be rewritten as follows:

$$\phi \left(\stackrel{\rightharpoonup }{x},\stackrel{\rightharpoonup }{y},t\right)=\mathrm{Re}\left\{a{e}^{is}{e}^{-i\omega t}\right\}$$

(16)

The potential wave function in the equation only represents the forward-scattered wave component. It neglects the contribution from the reflected wave, which means it is not applicable when there is a significant reflection effect on the structure. By substituting Equation (16) with Equation (14) and considering the real and imaginary parts separately, the following expressions can be obtained:

Real part:

$$\frac{1}{aC{C}_g}\left\{\left(\stackrel{\rightharpoonup }{U}\cdot \nabla a\right)\left[\left(\stackrel{\rightharpoonup }{U}\cdot \nabla \right)+\left(\nabla \cdot \stackrel{\rightharpoonup }{U}\right)\right]\right\}-\frac{1}{a}\left[{\nabla }^{2}a+\frac{1}{C{C}_g}\left(\nabla C{C}_g\cdot \nabla a\right)\right]-k^2+{\left|\nabla s\right|}^2=0$$

(17)

Imaginary part:

$$\nabla \cdot \left[a^2\sigma \left(U+C_g\right)\right]=0$$

(18)

Equations (17) and (18) represent the motion equations interacting with breaking waves and currents. When the flow velocity $\stackrel{\rightharpoonup }{U}$ is known, these two simultaneous equations can be solved to obtain the amplitude function $a\left(x,y\right)$ and the wave number $\left|\nabla s\right|$. When $\stackrel{\rightharpoonup }{U}=0$, Equations (17) and (18) become:

Real part:

$$\frac{1}{a}\left\{\frac{{\partial }^{2}a}{\partial x^2}+\frac{{\partial }^{2}a}{\partial y^2}+\frac{1}{C{C}_g}\left[\nabla a\cdot \nabla \left(C{C}_g\right)\right]\right\}+k^2-\left|\nabla s\right|=0$$

(19)

Imaginary part:

$$\nabla \cdot \left[a^{2}C{C}_g\nabla s\right]=0$$

(20)

Furthermore, within the surf zone, energy dissipation occurs, and therefore, Equation (18) for the energy representation needs to be modified. Based on the principle of energy flux, the bed friction effects are neglected in this case:

$$\frac{d\left(E{C}_g\right)}{dx}=-\epsilon \begin{array}{cc},& \epsilon =\frac{1}{2}\rho V_e{\left(k{H}_B\right)}^2\end{array}\begin{array}{cc},& V_e=V_{eB}{\left(\frac{H_B/2-c\prime h_B}{\gamma \prime h_B}\right)}^m\end{array}\begin{array}{cc},& V_{eB}=\frac{5S_{B}g}{8k_B\rho }\end{array}\frac{1}{\sqrt{1-C_0}}\begin{array}{cc},& S_B=\frac{\tan \beta }{1+3r^2/2}\end{array}$$

(21)

where $c\prime$ represents the wave number-to-water depth ratio in the recovery zone; it is suggested that $c\prime$ should be set to 0.17 when the wave recovery zone is not significant in a sloping beach condition [43]. In the coexistence of wave and current interaction zone, the energy dissipation caused by nearshore currents within the surf zone is negligible and can be ignored. Therefore, according to Equation (21) for energy amplitude representation, it can be expressed as follows:

$$\nabla \cdot \left[\frac{E}{\sigma }\left(\stackrel{\rightharpoonup }{U}+C_g\right)\right]=-\frac{5}{16}\frac{\rho g^2{k}_B}{{\sigma }^2}\frac{\tan \beta }{1+3r{\prime }^2/2}\frac{1}{\sqrt{1-C_0}}\sqrt{\frac{H_B/2-c\prime h_B}{r\prime h_B}}{\left(H_B\right)}^2$$

(22)

Using Equation (19), the energy representation within the surf zone is modified as follows, in Equation (22):

$$\nabla \cdot \left[a^2\sigma \left(\stackrel{\rightharpoonup }{U}+C_g\right)\right]=\nabla \cdot \left[\frac{2g}{\rho }\frac{E}{\sigma }\left(\stackrel{\rightharpoonup }{U}+C_g\right)\right]=-\frac{5}{8}\frac{g^2{k}_B}{{\sigma }^2}\frac{\tan \beta }{1+3r{\prime }^2/2}\frac{1}{\sqrt{1-c\prime /r\prime }}\sqrt{\frac{H_B/2-c\prime h_B}{r\prime h_B}}{\left(H_B\right)}^2$$

(23)

In Equations (21)–(23), the subscript “$B$” indicates the values within the surf zone. Additionally, since the phase function of $\varphi$ is given as $x\left(\stackrel{\rightharpoonup }{x},t\right)=s\left(\stackrel{\rightharpoonup }{x}\right)-\omega t$, the wave number obtained from the modified shallow water wave equation can be expressed as follows:

$$\stackrel{\rightharpoonup }{k}=\nabla x=\nabla s$$

(24)

To determine the wave number $\left|\nabla s\right|$ from Equations (17), (18), or (23), the knowledge of the wave direction angle is required. In other words, with only two equations involving $a$, $\left|\nabla s\right|$, and $\theta$, there needs to be more information. However, assuming non-rotational gradients of the wave phase function based on linear wave theory, we can obtain the following expression:

$$\begin{array}{l}\nabla \times \left(\nabla s\right)=0\\ \nabla s=\left|\nabla s\right|\cos \theta \stackrel{\rightharpoonup }{i}+\left|\nabla s\right|\sin \theta \stackrel{\rightharpoonup }{j}\\ \frac{\partial }{\partial x}\left(\left|\nabla s\right|\sin \theta \right)-\frac{\partial }{\partial y}\left(\left|\nabla s\right|\cos \theta \right)=0\end{array}$$

(25)

The tidal effect on the nearshore flow field is incorporated into the computational model in this study. Additionally, the tides are treated as long waves, considering the Coriolis effect due to the Earth’s rotation. The governing equations for the flow field, including the continuity equation and the momentum equation, can be expressed as follows:

Continuity equation:

$$\frac{\partial \eta }{\partial t}+\frac{\partial }{\partial x}\left[U\left(h+\eta \right)\right]+\frac{\partial }{\partial y}\left[V\left(h+\eta \right)\right]=0$$

(26)

Momentum equation:

$$\frac{\partial U}{\partial t}+U\frac{\partial U}{\partial x}+V\frac{\partial U}{\partial y}=fV-g\frac{\partial \eta }{\partial x}+\frac{1}{\rho }\left(\frac{\partial {\tau }_{xx}}{\partial x}+\frac{\partial {\tau }_{yx}}{\partial y}\right)+\frac{1}{\rho \left(h+\eta \right)}\left({\tau }_{sx}-{\tau }_{bx}\right)-\frac{1}{\rho \left(h+\eta \right)}\left(\frac{\partial S_{xx}}{\partial x}+\frac{\partial S_{yx}}{\partial y}\right)$$

(27)

$$\frac{\partial V}{\partial t}+U\frac{\partial V}{\partial x}+V\frac{\partial V}{\partial y}=fV-g\frac{\partial \eta }{\partial y}+\frac{1}{\rho }\left(\frac{\partial {\tau }_{xy}}{\partial x}+\frac{\partial {\tau }_{yy}}{\partial y}\right)+\frac{1}{\rho \left(h+\eta \right)}\left({\tau }_{sy}-{\tau }_{by}\right)-\frac{1}{\rho \left(h+\eta \right)}\left(\frac{\partial S_{xy}}{\partial x}+\frac{\partial S_{yy}}{\partial y}\right)$$

(28)

The variables $U$ and $V$ represent the depth-averaged velocity components along the fixed coordinate axes $x$ and $y$, respectively. Specifically:

$$U=\frac{1}{\left(h+\eta \right)}{\displaystyle {\int }_{-h}^{\eta }udz}\begin{array}{cc},& V=\end{array}\frac{1}{\left(h+\eta \right)}{\displaystyle {\int }_{-h}^{\eta }vdz}$$

(29)

The shear stress component ${\tau }_{xx}$, ${\tau }_{xy}$, ${\tau }_{yx}$, and ${\tau }_{yy}$ includes the viscous stress caused by fluid viscosity and the Reynolds stress generated by turbulent flow effects. Comparatively, the value of viscous stress is minimal compared to Reynolds stress. Therefore, viscous stress is generally neglected, and only the Reynolds stress represents the momentum exchange between fluid particles.

$${\tau }_{xx}=\rho E_v\frac{\partial U}{\partial x}\begin{array}{cc},& {\tau }_{xy}=\rho E_v\frac{\partial U}{\partial y}\end{array}\begin{array}{cc},& {\tau }_{yx}=\rho E_v\frac{\partial V}{\partial x}\end{array}\begin{array}{cc},& {\tau }_{yy}=\rho E_v\frac{\partial V}{\partial y}\end{array}$$

(30)

where the vortex viscosity coefficient is obtained from the semi-empirical formula of Prandtl’s mixing length theory [44].

$$E_v=\frac{k_v\sqrt{g}\left(d+h\right)\sqrt{U^2+V^2}}{6C_c}$$

(31)

The components of the wind shear stress on the sea surface represent the shear stress in the x and y directions, respectively [45].

$${\tau }_{sx}=\rho k_w{W}^2\cos a\begin{array}{cc},& {\tau }_{sy}\end{array}=\rho k_w{W}^2\sin a$$

(32)

$$k_w=\left\{\begin{array}{lll}1.2\times {10}^{-5}\hfill & ,\hfill & W\le W_c\hfill \\ 1.2\times {10}^{-6}+2.25\times {10}^{-1}{\left[1-\frac{W_c}{W}\right]}^2\hfill & ,\hfill & W>W_c\hfill \end{array}\right.$$

(33)

The components ${\tau }_{bx}$ and ${\tau }_{by}$ of the bottom friction force represent the frictional shear stress in the x and y directions, respectively [46].

$${\tau }_{bx}=\rho E_{r}U\sqrt{U^2+V^2}\begin{array}{cc},& {\tau }_{by}=\rho E_{r}V\sqrt{U^2+V^2}\end{array}$$

(34)

The bottom friction coefficient is denoted as $F_r=g/C_c^2$. The radiative stress components $S_{xx}$, $S_{xy}$, $S_{yx}$, and $S_{yy}$ are the main factors causing longshore currents, represented by the linear wave theory, as follows [47]:

$$\left[\begin{array}{cc}{S}_{xx}& S_{xy}\\ S_{yx}& S_{yy}\end{array}\right]=\overline{E}\left[\begin{array}{cc}n\left(1+{\cos }^2\theta \right)-\frac{1}{2}& \left(\frac{n}{2}\right)\sin \left(2\theta \right)\\ \left(\frac{n}{2}\right)\sin \left(2\theta \right)& n\left(1+{\sin }^2\theta \right)-\frac{1}{2}\end{array}\right]$$

(35)

where $\overline{E}$ represents the total wave energy per unit time per unit area, in the theory of small-amplitude water waves, which can be expressed as:

$$\overline{E}=\frac{\rho g{H}^2}{8}$$

(36)

The boundary conditions for the flow field are shown in Figure 8, and the water level variation includes the elevation changes caused by waves and tides.

The tidal variation is represented by a sinusoidal function, with a time lag ($T_L$) indicating the time difference for the tide to reach the left and right boundaries. $T_{start}$ represents the phase difference:

Left boundary:

$${\eta }_L=A_L^t\sin \left[\frac{2\pi }{T_t}\left(t+T_t+T_{start}\right)\right]\begin{array}{cc},& T_t=\frac{L_y}{\sqrt{g{h}_{\max }}}\end{array}$$

(37)

Right boundary:

$${\eta }_R=A_R^t\sin \left[\frac{2\pi }{T}\left(t+T_{start}\right)\right]$$

(38)

Offshore boundary:

$${\eta }_0=\left[A_R^t+\left(A_L^t-A_R^t\right)\left(\frac{N_y-j}{N_y-1}\right)\right]\sin \left\{\frac{2\pi }{T_t}\left[t+T_t\left(\frac{N_y-j}{N_y-1}\right)+T_{start}\right]\right\}$$

(39)

The water level induced by a wave ($\overline{\xi }$) neglects the reflection effect [48]. Equation (40) represents the water level drop outside the surf zone, while Equation (41) represents the water level rise inside the surf zone:

$$\overline{{\xi }_d}=-\frac{H^2}{8}\frac{k}{\sinh \left(2kh\right)}{\left(\cos \theta \right)}^{2/3}$$

(40)

$$\frac{d\overline{{\xi }_u}}{dx}=-K\frac{dh}{dx}\begin{array}{cc},& K=\frac{1}{1+\left(8/3{\gamma }^2\right)}\end{array}$$

(41)

The left boundary of the flow velocity is given by Equation (42), the right boundary is given by Equation (43), the offshore boundary is given by Equation (44), and the alongshore boundary is given by Equation (45):

$$U_{j=1}=U_{j=2}\begin{array}{cc},& {\left(\frac{\partial V}{\partial y}\right)}_{j=1}\end{array}=0$$

(42)

$$U_{j=NY}=U_{j=NY-1}\begin{array}{cc},& {\left(\frac{\partial V}{\partial y}\right)}_{j=NY}\end{array}=0$$

(43)

$$U_{j=NX}=U_{j=NX-1}\begin{array}{cc},& {\left(\frac{\partial U}{\partial x}\right)}_{j=NX}\end{array}=0$$

(44)

$$U=0\begin{array}{cc},& V=0\end{array}$$

(45)

The stability condition for the flow field calculation requires $\mathsf{\Delta }t\le \frac{2\mathsf{\Delta }s}{\sqrt{g{h}_{\max }}}$, where $\mathsf{\Delta }s$ is the grid size. The maximum value of the two consecutive time steps must be smaller than the allowable error with Equation (46), and then proceed to the next time step of the calculation:

$$\begin{array}{l}Max\left({\eta }_{ij}^{k+1}-{\eta }_{ij}^k\right)\le {\epsilon }_{\eta }{\eta }_{ij}^k\begin{array}{cc},& {\epsilon }_{\eta }=0.0001\end{array}\\ Max\left(U_{ij}^{k+1}-U_{ij}^k\right)\le {\epsilon }_U{U}_{ij}^k\begin{array}{cc},& {\epsilon }_U=0.0001\end{array}\\ Max\left(V_{ij}^{k+1}-V_{ij}^k\right)\le {\epsilon }_V{V}_{ij}^k\begin{array}{cc},& {\epsilon }_V=0.0001\end{array}\end{array}$$

(46)

3.2. Model Calculation Range and Settings

The model calculation area in this study is shown in Figure 9. It is centered around Keelung Harbor and covers a coastal distance of 5 km in the alongshore direction and 6 km in the offshore direction. The numerical topography was established using interpolated bathymetric data with a grid resolution of 200 m obtained from the National Science and Technology Council—Taiwan. The overall computational conditions of the numerical model mainly refer to the statistical data obtained from on-site observations conducted by the acoustic doppler current profiler (ADCP) installed at the port of Keelung by the Institute of Transportation, MOTC, Harbor and Marine Technology Center, as shown in Figure 1. The model settings and conditions are summarized in

Table 2. Simulation conditions.

Item	Model Setup
Area	5 km × 6 km
Grid size	$\mathsf{\Delta }x=\mathsf{\Delta }y=$25 m
Number of grid points	200 × 240
Time step	$\mathsf{\Delta }t=2s$
Time scale	1.0
Wave condition	2.5 m; 8.0 s; NNE

. The selection of the model calculation time considers the need to effectively reflect the characteristics of water flow distribution near Keelung Harbor. Therefore, the calculation is performed during the spring tide period, with a simulation period starting from minor wave conditions and gradually increasing. The comparison of model calculation results with field data is shown in Figure 10. The figure demonstrates that the water level, flow velocity, flow direction, wave height, and wave direction obtained from the model calculations generally match the field data, indicating that the model can reasonably represent the hydrodynamic characteristics in the vicinity of Keelung Harbor. Subsequently, flow field simulations were conducted for Keelung Harbor under different sea-level-rise scenarios using the validated boundary conditions and parameters. The effects of sea level rise on the regional flow field were analyzed, particularly concerning ship navigation safety.

3.3. Hydraulic Simulation Results Analysis

The simulated results of the model, as verified in the previous section, are presented in Figure 11. The flow characteristics in the nearshore area primarily reflect the interaction between tidal motion and the features of the nearshore topography. In the vicinity of Keelung Port, the ebb flow generally moves from west to east, while the flood flow moves from east to west. It is observed that there is an increase in flow velocity between Keelung Island and Keelung Port due to the influence of the topography. Along the outer channel to the vicinity of the revetment, during the flood tide, the average flow velocity ranges from 0.3925 to 0.6945 m/s. After entering the breakwater, the average flow velocity ranges from 0.1209 to 0.2400 m/s. The distribution of flow directions indicates the formation of a nearshore circulation pattern due to the influence of the eastern breakwater head. During the ebb tide, the average flow velocity from the outer channel to the revetment ranges from 0.6094 to 0.6906 m/s. After entering the breakwater, the average flow velocity ranges from 0.1421 to 0.4722 m/s. The flow direction is influenced by the structures, transitioning from south to north in the channel before bypassing the breakwater and then flowing eastward.

Due to the influence of environmental factors such as wind and currents, the waters from the outer channel to the vicinity of the revetment are primarily affected. Ship maneuvering is a complex process involving human, vessel, and environmental factors. The ship’s operators must process and comprehend much information to maintain or alter the vessel’s motion to achieve the intended maneuvering objectives. Various ships may face threats to their navigation in the channels due to external forces such as crosswinds and currents, which can prevent them from staying in the center of the channel. The selection of basic flow parameters is crucial for safe ship navigation. Suppose the regional water flow direction differs from the vessel’s heading. In that case, it may create a transverse flow that alters the hydrodynamic forces, pushing the ship off its intended course and compromising its safety. Large vessels are particularly susceptible to external forces during inbound and outbound sailing, making maneuvering more challenging than other vessels. Therefore, it is necessary to consider a combination of environmental forces and harbor characteristics to ensure safety. The hydraulic simulation results are obtained using different analysis methods, as shown in Figure 12, Figure 13 and Figure 14, to analyze the effects of sea level rise. For better readability, reference is made to the existing harbor channel planning in Keelung Port. The simulated results are extracted near the “main fairway sector” indicated in Figure 1, and data extraction and discussion are conducted along the inbound fairway at 150 m intervals.

The simulation results, shown in Figure 15 and Figure 16 and summarized in

Table 3. Statistics table of flow velocity and direction within the outer fairway to the breakwater (a–i) under flood tide conditions.

	Original		SMA		EMD		EEMD
	V (m/s)	Dir.	V (m/s)	Dir.	V (m/s)	Dir.	V (m/s)	Dir.
a	0.6943	272.72	0.6945	272.72	0.6946	272.73	0.6947	272.73
b	0.6706	271.07	0.6707	271.08	0.6709	271.08	0.6709	271.08
c	0.6033	269.17	0.6034	269.18	0.6035	269.20	0.6036	269.20
d	0.4952	268.05	0.4953	268.08	0.4953	268.13	0.4954	268.13
e	0.3925	266.77	0.3926	266.82	0.3926	266.90	0.3927	266.90
f	0.2400	293.97	0.2401	294.03	0.2402	294.12	0.2402	294.11
g	0.1623	351.56	0.1624	351.59	0.1626	351.63	0.1626	351.62
h	0.1975	27.90	0.1976	27.87	0.1977	27.80	0.1977	27.80
i	0.1209	31.43	0.1210	31.33	0.1210	31.12	0.1211	31.13

and

Table 4. Statistics table of flow velocity and direction within the outer fairway to the breakwater (a–i) under ebb tide conditions.

	Original		SMA		EMD		EEMD
	V (m/s)	Dir.	V (m/s)	Dir.	V (m/s)	Dir.	V (m/s)	Dir.
a	0.6834	84.04	0.6866	83.98	0.6871	83.97	0.6901	83.91
b	0.6906	80.70	0.6940	80.65	0.6945	80.64	0.6977	80.58
c	0.6714	75.96	0.6749	75.90	0.6754	75.89	0.6786	75.84
d	0.6370	67.68	0.6405	67.63	0.6410	67.62	0.6441	67.57
e	0.6094	59.68	0.6129	59.64	0.6135	59.63	0.6166	59.59
f	0.4722	42.57	0.4752	42.54	0.4757	42.53	0.4783	42.51
g	0.3531	33.85	0.3555	33.83	0.3561	33.83	0.3580	33.80
h	0.2614	38.18	0.2634	38.11	0.2638	38.15	0.2654	38.07
i	0.1421	70.36	0.1429	70.21	0.1435	70.25	0.1440	70.11

, provide a visual representation and quantification of the current pattern along the inbound and outbound fairways at 150 m intervals. In the vicinity of the outer fairway to the entrance (a–e), the simulated results under different sea-level-rise conditions provide valuable insights when incorporating different methods such as SMA, EMD, and EEMD. During the rise tide, the average current velocity ranges from 0.3926 to 0.6945 m/s when utilizing the SMA method, from 0.3926 to 0.6946 m/s with the EMD method, and from 0.3927 to 0.6947 m/s with the EEMD method. On the ebb tide, the average current velocity ranges from 0.6129 to 0.6866 m/s (SMA), 0.6135 to 0.6871 m/s (EMD), and 0.6166 to 0.6901 m/s (EEMD). After entering the entrance (f–e), the current velocity generally decreases as it enters the harbor area. During the flood tide, the average current velocity ranges from 0.1210 to 0.2401 m/s (SMA), 0.1210 to 0.2402 m/s (EMD), and 0.1211 to 0.2402 m/s (EEMD). On the ebb tide, the average current velocity ranges from 0.1429 to 0.4752 m/s (SMA), 0.1435 to 0.4757 m/s (EMD), and 0.1440 to 0.4783 m/s (EEMD). It is observed that the current velocity generally increases from the outer fairway to inside the entrance (a–i) as sea levels rise, indicating the influence of hydrodynamics and the impact of changing water levels on current patterns. Conversely, the distribution of flow direction shows less significant changes, suggesting that the dominant direction of current remains relatively consistent despite variations in sea level.

The results of the simulations provide a basis for discussion and a theoretical analysis of the flow velocity patterns under different sea-level-rise conditions in the studied area, particularly in the vicinity of the outer fairway to the entrance (a–e). By utilizing the simulation results, potential hazards and issues can be identified, aiding mariners in conducting safety assessments and devising effective navigational strategies based on the actual distribution of current. Overall, this comprehensive analysis of the simulation results provides valuable insights into the current characteristics within the studied area, aiding in understanding the impact of sea level rise and enhancing navigational safety considerations.

4. Conclusions

In the past, sea level rise has been assessed using linear regression and the least squares method to determine coefficients. Additionally, traditional harmonic analysis methods have been unable to analyze sea level rise. Therefore, this study examines the impact of sea level rise on ship navigation safety by analyzing changes in current characteristics. Three analysis methods, SMA, EMD, and EEMD, were employed to investigate sea-level-rise trends from 2005 to 2022 at the Keelung tide gauge station. The SMA analysis indicated an average sea level rise of 0.98 mm/yr, totaling approximately 17.68 mm. The EMD analysis showed an average rise of 2.84 mm/yr, reaching around 51.20 mm. The EEMD analysis revealed the highest average rise of 3.21 mm/yr, approximately 57.73 mm. However, the EEMD analysis also indicated a declining trend, possibly influenced by data gaps and endpoint effects. It was challenging to ascertain a continuous downward trend in water levels at the Keelung tide gauge station.

Due to the external forces acting on a vessel, including wind (on the water surface), currents (underwater), and waves (affecting vessel stability), understanding the magnitude and direction of current is crucial for effective vessel maneuvering. When a vessel carries cargo, its draft increases, increasing the current impact area. The velocity and direction of the currents are influenced by harbor structures, leading to variations. Moreover, the current pressure exerted on the vessel by the current changes with the vessel’s heading angle during navigation. Therefore, understanding the magnitude and direction of current is significantly beneficial for subsequent vessel operations. Therefore, using the obtained sea-level-rise data, hydrodynamic simulations were conducted, focusing on the outer fairway to the entrance area. During the flood tide, current velocities ranged from 0.3926 to 0.6945 m/s (SMA), 0.3926 to 0.6946 m/s (EMD), and 0.3927 to 0.6947 m/s (EEMD). On the ebb tide, current velocities in the same area ranged from 0.6129 to 0.6866 m/s (SMA), 0.6135 to 0.6871 m/s (EMD), and 0.6166 to 0.6901 m/s (EEMD). Upon entering the harbor area, current velocities generally decreased. During the flood tide, velocities ranged from 0.1210 to 0.2401 m/s (SMA), 0.1210 to 0.2402 m/s (EMD), and 0.1211 to 0.2402 m/s (EEMD). On the ebb tide, velocities ranged from 0.1429 to 0.4752 m/s (SMA), 0.1435 to 0.4757 m/s (EMD), and 0.1440 to 0.4783 m/s (EEMD). The current velocity distribution from the outer fairway to the entrance generally increased with sea level rise, while changes in current direction distribution were less pronounced. Since selecting basic current parameters is crucial for vessel navigation, these analysis results offer valuable insights for mariners to assess navigation safety and devise appropriate strategies based on the actual flow conditions. Furthermore, they can serve as a reference for maritime and management authorities in ensuring safe navigation in harbor fairways.