*3.1. Model Description*

A finite volume community ocean model (FVCOM) was applied in this study. FVCOM is widely used for investigating coastal ocean hydrodynamics with complicated topography [52]. The governing equations are [52]:

$$\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial \mathbf{x}} + v \frac{\partial u}{\partial y} + w \frac{\partial u}{\partial z} - fv = -\frac{1}{\rho\_0} \frac{\partial P}{\partial \mathbf{x}} + \frac{\partial}{\partial z} \left( K\_m \frac{\partial u}{\partial z} \right) + F\_{u\prime} \tag{1}$$

$$\frac{\partial v}{\partial t} + u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y} + w\frac{\partial v}{\partial z} + fu = -\frac{1}{\rho\_0}\frac{\partial P}{\partial y} + \frac{\partial}{\partial z}\left(K\_m\frac{\partial v}{\partial z}\right) + F\_{\upsilon\_\tau} \tag{2}$$

$$\frac{\partial w}{\partial t} + u \frac{\partial w}{\partial x} + v \frac{\partial w}{\partial y} + w \frac{\partial w}{\partial z} = -\frac{1}{\rho\_0} \frac{\partial p}{\partial z} + \frac{\partial}{\partial z} \Big(K\_m \frac{\partial w}{\partial z}\Big) + F\_w - g\_\prime \tag{3}$$

$$
\frac{
\partial \mu
}{
\partial \mathbf{x}
} + \frac{
\partial v
}{
\partial y
} + \frac{
\partial w
}{
\partial z
} = 0,
\tag{4}
$$

where *x*, *y*, *z* are the three axes in the east, north, and vertical directions, respectively; *u*, *v*, *w* are the *x*, *y*, *z* velocities, respectively, ρ0 is the density; P is the total pressure of air and water; *f* is the Coriolis parameter; *g* is the gravitational acceleration; *Km* is the vertical eddy diffusion coefficient, determined by the Mellor and Yamada [53] level-2.5 (MY-2.5) turbulent closure scheme; *Fw* is the diffusion term of the vertical momentum; and *Fu*, *Fv* are the diffusion terms for the horizontal momentums.

The surface and bottom boundary conditions are:

$$K\_m \left(\frac{\partial u}{\partial z}, \frac{\partial v}{\partial z}\right) = \frac{1}{\rho\_0} (\tau\_{sx}, \tau\_{sy}), \ w = \frac{\partial \zeta}{\partial t} + \frac{\partial \zeta}{\partial x} + \frac{\partial \zeta}{\partial y'} \text{ at } z = \zeta(x, y, t), \tag{5}$$

$$K\_m \left(\frac{\partial \mu}{\partial z'}, \frac{\partial v}{\partial z}\right) = \frac{1}{\rho\_0} (\tau\_{\text{xx}}, \tau\_{\text{sy}}), \; w = -\mu \frac{\partial H}{\partial \mathbf{x}} - \mathbf{v} \frac{\partial H}{\partial y'}, \; \text{at } z = -H(\mathbf{x}, y), \tag{6}$$

where <sup>τ</sup>*sx*, <sup>τ</sup>*sy*and <sup>τ</sup>*bx*, <sup>τ</sup>*by*are the surface wind stress and bottom stress vectors, respectively. *H* is the water depth and ζ is the surface elevation. <sup>τ</sup>*sx*, <sup>τ</sup>*sy* is calculated by *Cd*ρ*a*|*U*10|*U*10, where *U*10 is the wind at 10 m height, ρ*a* is the air density (1.29 kg/m3), and *Cd* is the surface wind drag coefficient and is calculated by the following equations:

$$\mathbf{C}\_{d} = \begin{cases} (0.49 + 0.065 \times 11.0) \times 10^{-3}, \; \mathcal{U}l\_{10} < 11.0 \; m/s \\\ (0.49 + 0.065 \times |\mathcal{U}\_{10}|) \times 10^{-3}, \; 11.0 \frac{m}{s} \le \mathcal{U}\_{10} \le 25.0 \; m/s \\\ (0.49 + 0.065 \times 25.0) \times 10^{-3}, \; \mathcal{U}l\_{10} > 25.0 \; m/s \end{cases} \tag{7}$$

where <sup>τ</sup>*bx*, <sup>τ</sup>*by* is the bottom stress calculated by *Cd* √*u*<sup>2</sup> + *<sup>v</sup>*<sup>2</sup>(*<sup>u</sup>*, *<sup>v</sup>*), where *Cd* is the drag coefficient and is determined by the following equation:

$$C\_d = \max(\frac{k^2}{\ln\left(\frac{z\_{\text{eff}}}{z\_0}\right)^2}, 0.0025), \tag{8}$$

where *k* is the von Karman constant (0.4), *z*0 is the bottom roughness parameter, and *zab* is the height above the bottom.
