1 Introduction

In tidal-dominated shallow seas, the tidal residual current controls the transport of the conservative substances beyond the time scale of tidal periods (Abbott, 1960; Nihoul and Ronday, 1975). Furthermore, the use of the Lagrangian velocity (u_LR), instead of the Eulerian velocity (u_E), is more appropriate in representing the mass transport. In weakly nonlinear systems, the first-order approximation of u_LR is the mass transport velocity u_L, first introduced by Longuet-Higgins (1969) in the study of wave-generated circulation in the open ocean. The governing equations for u_L in a three-dimensional (3D) case were first derived by Feng (1987). The equations include the 'tidal body force' terms resulting from the nonlinear interaction of the zeroth-order tides. The governing equations reveal the dynamic factors contributing to u_L; however, in practice, particle tracking is commonly used to obtain u_LR (e.g., Zimmerman, 1979; Cheng and Casulli, 1982; Feng et al., 2008; Muller et al., 2009; Liu et al., 2012; Liang et al., 2014; Quan et al., 2014; Xu et al., 2016; Deng et al., 2017, 2019).

The analytic solution of u_L is available for a few special cases (e.g., Ianniello, 1977, 1979; Winant, 2008; Jiang and Feng, 2011, 2014). In most other cases without available analytic solutions, numerical models are applied to solve u_L (e.g., Wang et al., 1993; Sun et al., 2000, 2001; Liu et al., 2002; Lei et al., 2004; Wang et al., 2006). In early studies, the accuracy of the modeled u_L was limited by the coarse spatial resolution of the models. Recently, Cui et al. (2019) assessed the accuracy of modeled u_Lfor the special case with analytic solution obtained previously by Jiang and Feng(2011, 2014). However, while the model and analytic solution show overall good agreement, the discrepancy remains noticeable.

In Cui et al. (2019), the tidal body force is computed by tidal-averaging the terms involving the time series of the tidal flow. In this process, the accuracy of the computed tidal body force is affected by errors related to discretization, numerical integration, and the open boundary condition, among others. In this study, we propose to compute the tidal body force using the harmonic constants of the tidal current obtained by performing a harmonic analysis of the tidal solution. After introducing the new method in Section 2, the solution of u_L is first obtained for an elongated rectangular bay for comparison with the analytic solution of Jiang and Feng (2014) in Section 3. In Section 4, u_L is calculated for a wide bay with different geometrical and physical parameters. The impacts of these parameters on u_L are discussed. The conclusions from this study are summarized in Section 5.

2 Improved Computing Method for Solving u_L 2.1 Governing Equations

The governing equations of u_L are the same as those in Feng (1987) and Cui et al. (2019). A summary of the key aspects is provided below. First, the equations are in nondimension form involving a parameter κ, which is the ratio of the tidal amplitude to the average water depth. In the weakly nonlinear case, i.e., O(κ) < 1, the tidal velocity and water elevation are presented as

$u = {u_0} + \kappa {u_1} + O({\kappa ^2}), $

$\zeta = {\zeta _0} + \kappa {\zeta _1} + O({\kappa ^2}), $

where u₀ and ${\zeta _0}$ are the zeroth-order tidal velocity and elevation, respectively, and u₁ and ζ₁ are the first-order velocity and elevation, respectively.

The residual water elevation ζ_E is defined as ζ_E = < ζ₁ > , and the zeroth-order tidal displacement ξ₀ is expressed as

${\xi _0} = (\xi, {\kern 1pt} {\kern 1pt} {\kern 1pt} \eta, {\kern 1pt} {\kern 1pt} {\kern 1pt} \varsigma) = \int_{{t_0}}^t {{u_0}(x, {\rm{ }}t')} {\rm{d}}t'.$

Meanwhile, the Eulerian residual velocity u_E is defined as

${u_E} = ({u_E}, {v_E}, {w_E}) = (< {u_1} >, < {v_1} >, < {w_1} >) = < {u_1} > {\kern 1pt}, $

where the tidal-averaging operator is expressed as follows

$ < \cdot > \; = \frac{1}{{nT}}\int_{{t_0}}^{{t_0} + nT} \cdot \;{\rm{d}}t.$

The mass transport velocity u_L named by Longuet-Higgins (1969) is the sum of the Eulerian residual velocity and the Stokes' drift, i.e.

${u_L} = {u_E} + {u_S}, $

where the Stokes' drift u_S is defined as

${u_S} = ({u_S}, \;{v_S}, \;{w_S}) = < {\xi _0} \cdot \nabla {u_0} >. $

The governing equations for the zeroth-order tidal motions are as follows:

$\nabla \cdot {u_0} = 0, $

(1)

$\frac{{\partial {u_0}}}{{\partial t}} - f{v_0} = - g\frac{{\partial {\zeta _0}}}{{\partial x}} + \frac{\partial }{{\partial z}}\left( {\upsilon \frac{{\partial {u_0}}}{{\partial z}}} \right){\kern 1pt}, $

(2)

$\frac{{\partial {v_0}}}{{\partial t}} - f{u_0} = - g\frac{{\partial {\zeta _0}}}{{\partial y}} + \frac{\partial }{{\partial z}}\left( {\upsilon \frac{{\partial {v_0}}}{{\partial z}}} \right){\kern 1pt}, $

(3)

$ z = 0, {w_0} = \frac{{\partial {\zeta _0}}}{{\partial t}};\frac{{\partial \left({{u_0}, {\kern 1pt} {\kern 1pt} {\kern 1pt} {v_0}} \right)}}{{\partial z}} = 0, $

(4)

$ z = - h, {u_0} = 0;{v_0} = 0;{w_0} = 0. $

(5)

The governing equations for u_L are as follows:

$\nabla \cdot {u_L} = 0{\kern 1pt}, $

(6)

$ -f{v_L} = - g\frac{{\partial {\zeta _E}}}{{\partial x}} + \frac{\partial }{{\partial z}}\left({\upsilon \frac{{\partial {u_L}}}{{\partial z}}} \right) + {\pi _x}, $

(7)

$f{u_L} = - g\frac{{\partial {\zeta _E}}}{{\partial y}} + \frac{\partial }{{\partial z}}\left({\upsilon \frac{{\partial {v_L}}}{{\partial z}}} \right) + {\pi _y}, $

(8)

$ z = 0, {w_L} = 0;\frac{{\partial \left({{u_L}, {\kern 1pt} {\kern 1pt} {\kern 1pt} {v_L}} \right)}}{{\partial z}} = 0, $

(9)

$ z = - h, {u_L} = 0;{v_L} = 0;{w_L} = 0, $

(10)

${\pi _x} = - \frac{\partial }{{\partial z}}\left( {\upsilon \frac{{\partial < {\xi _0} \cdot \nabla {u_0} > }}{{\partial z}}} \right) - < {u_0} \cdot \nabla {u_0} > - f < {\xi _0} \cdot \nabla {v_0} > , $

(11)

${\pi _y} = - \frac{\partial }{{\partial z}}\left( {\upsilon \frac{{\partial < {\xi _0} \cdot \nabla {v_0} > }}{{\partial z}}} \right) - < {u_0} \cdot \nabla {v_0} > + f < {\xi _0} \cdot \nabla {u_0} >, $

(12)

where f, g, υ, and h are the non-dimensional Coriolis parameter, gravitational acceleration, vertical eddy viscosity coefficient, and water depth, respectively.

The tidal body force π = (π_x, π_y) (Eq. (11) and Eq. (12)) consists of three terms (Cui et al., 2019): π₁ = (π_x1, π_y1) the vertical mixing of the Stokes' drift, π₂ = (π_x2, π_y2) the contribution of the advection, and π₃ = (π_x3, π_y3) the Coriolis effect of the local acceleration.

Eqs. (1)–(5) are commonly solved by using timestepping tidal models. Eqs. (6)–(12) can also be solved using the time-stepping method to reach a steady state (Cui et al., 2019). The key step is to compute the tidal body force terms (π_x, π_y).

2.2 Improving the Tidal Body Force Computation

In Cui et al. (2019), the tidal body force (Eqs. (11) and (12)) was calculated by applying tidal-averaging over the nonlinearity coupled tidal solutions. In the current work, we propose an alternative approach. First, the HAMSOM model (Backhaus, 1985; Pohlmann, 1996, 2006) was applied to obtain the time series of the zeroth-order tidal variables at half hourly intervals over a period of 8 days. Then, the tidal harmonic analysis package T_TIDE (Pawlowicz et al., 2002) was applied to obtain the tidal amplitudes and phases from the time series. This procedure can effectively filter out the clutter other than the single-frequency tide input from the open boundary and suppress the noise that can introduce errors into the tidal body force.

The time series of the three components of tidal current are represented as

${u_0} = {A_u}\cos (\frac{{2{\rm{ \mathsf{ π} }}}}{T}t -{\varphi _u}), $

${v_0} = {A_v}\cos (\frac{{2{\rm{ \mathsf{ π} }}}}{T}t -{\varphi _v}), $

${w_0} = {A_w}\cos (\frac{{2{\rm{ \mathsf{ π} }}}}{T}t- {\varphi _w}), $

where A_u, A_v, A_w and φ_u, φ_v, φ_w are the amplitudes and phases of u₀, v₀, w₀ respectively. The corresponding zeroth-order tidal displacements are expressed as

${\xi _0} = \frac{T}{{2{\rm{ \mathsf{ π} }}}}{A_u}[\sin (\frac{{2{\rm{ \mathsf{ π} }}}}{T}t- {\varphi _u}) + \sin ({\varphi _u})], $

${\kern 1pt} {\eta _0} = \frac{T}{{2{\rm{ \mathsf{ π} }}}}{A_v}[\sin (\frac{{2{\rm{ \mathsf{ π} }}}}{T}t- {\varphi _v}) + \sin ({\varphi _v})], $

${\varsigma _0} = \frac{T}{{2{\rm{ \mathsf{ π} }}}}{A_w}[\sin (\frac{{2{\rm{ \mathsf{ π} }}}}{T}t -{\varphi _w}) + \sin ({\varphi _w})].$

By substituting the above expressions into Eqs. (11) – (12), the tidal body force can be rewritten as

${\pi _x} = \\ -\frac{T}{{4{\rm{ \mathsf{ π} }}}}\frac{\partial }{{\partial z}}\left({\upsilon \frac{\partial }{{\partial z}}\left({A_u^2\frac{{\partial {\varphi _u}}}{{\partial x}} + A_u^{}A_v^{}\frac{{\partial {\varphi _u}}}{{\partial y}}\cos ({\varphi _u} -{\varphi _v}) + A_v^{}\frac{{\partial {A_u}}}{{\partial y}}\sin ({\varphi _u}- {\varphi _v})} \right.} \right. \\ + A_u^{}A_w^{}\frac{{\partial {\varphi _u}}}{{\partial z}}\left. {\left. {\cos ({\varphi _u} -{\varphi _w}) + A_w^{}\frac{{\partial {A_u}}}{{\partial z}}\sin ({\varphi _u} -{\varphi _w})} \right)} \right)\\ -\frac{1}{2}\left({A_u^{}\frac{{\partial {A_u}}}{{\partial x}} + A_u^{}A_v^{}\frac{{\partial {\varphi _u}}}{{\partial y}}\sin ({\varphi _v}- {\varphi _u}) + A_v^{}} \right.\frac{{\partial {A_u}}}{{\partial y}}\cos ({\varphi _u}- {\varphi _v}) \\ + A_u^{}A_w^{}\frac{{\partial {\varphi _u}}}{{\partial z}}\sin ({\varphi _w} -{\varphi _u})\left. { + A_w^{}\frac{{\partial {A_u}}}{{\partial z}}\cos ({\varphi _u}- {\varphi _w})} \right)\\ -\frac{{fT}}{{4{\rm{ \mathsf{ π} }}}}\left({A_u^{}\frac{{\partial {A_v}}}{{\partial x}}} \right.\sin ({\varphi _v} -{\varphi _u}) + A_u^{}A_v^{}\frac{{\partial {\varphi _v}}}{{\partial x}}\cos ({\varphi _v} -{\varphi _u}) + A_v^2\frac{{\partial {\varphi _v}}}{{\partial y}} \\+ A_w^{}\frac{{\partial {A_v}}}{{\partial z}}\sin ({\varphi _v}- {\varphi _w}) + A_w^{}\left. {{A_v}\frac{{\partial {\varphi _v}}}{{\partial z}}\cos ({\varphi _v} -{\varphi _w})} \right), $

(13)

${\pi _y} = \\ -\frac{T}{{4{\rm{ \mathsf{ π} }}}}\frac{\partial }{{\partial z}}\left({\upsilon \frac{\partial }{{\partial z}}\left({A_u^{}\frac{{\partial {A_v}}}{{\partial x}}\sin ({\varphi _v} -{\varphi _u}) + A_u^{}A_v^{}\frac{{\partial {\varphi _v}}}{{\partial x}}\cos ({\varphi _v} -{\varphi _u}) + A_v^2\frac{{\partial {\varphi _v}}}{{\partial y}}} \right.} \right. \\ +A_w^{}\frac{{\partial {A_v}}}{{\partial z}}\sin ({\varphi _v}\left. -{\left. { {\varphi _w}) + A_w^{}{A_v}\frac{{\partial {\varphi _v}}}{{\partial z}}\cos ({\varphi _v} -{\varphi _w})} \right)} \right)\\ -\frac{1}{2}\left({A_v^{}\frac{{\partial {A_v}}}{{\partial y}} + A_u^{}\frac{{\partial {A_v}}}{{\partial x}}\cos ({\varphi _u} -{\varphi _v}) + A_v^{}} \right.\frac{{\partial {\varphi _v}}}{{\partial x}}\sin ({\varphi _u}- {\varphi _v}) + \\ A_w^{}\frac{{\partial {A_v}}}{{\partial z}}\cos ({\varphi _v} -{\varphi _w})\left. { + A_w^{}A_v^{}\frac{{\partial {\varphi _v}}}{{\partial z}}\sin ({\varphi _w} -{\varphi _v})} \right)\\ {\rm{ + }}\frac{{fT}}{{4{\rm{ \mathsf{ π} }}}}\left({A_u^2\frac{{\partial {\varphi _u}}}{{\partial x}}} \right. + A_u^{}A_v^{}\frac{{\partial {\varphi _u}}}{{\partial y}}\cos ({\varphi _u} -{\varphi _v}) + A_v^{}\frac{{\partial {A_u}}}{{\partial y}}\sin ({\varphi _u} -{\varphi _v}) \\ + A_u^{}A_w^{}\frac{{\partial {\varphi _u}}}{{\partial z}}\cos ({\varphi _u} -{\varphi _w}) + A_w^{}\left. {\frac{{\partial {A_u}}}{{\partial z}}\sin ({\varphi _u}- {\varphi _w})} \right).$

(14)

Thus, the tidal body force can be computed using Eqs. (13)–(14) with the inputs of tidal current harmonic constants and without the need to perform tidal-averaging.

In Cui et al. (2019), the calculation of the tidal body force involves layer-integration and finite-differencing. In this study, this procedure involving the computation with Eqs. (13)–(14) is optimized. Specifically, the calculation of the last term related to the Coriolis parameter at its compute node takes the multi-point average of the finitedifference values of u_S.

3 Testing the Improved Computing Method

In this section, we tested the performance of the improved computing method for an elongated rectangular bay with varying topography across the bay. The configurations of the model are identical to those of Jiang and Feng (2014) and of Cui et al. (2019).

The distributions of the non-dimensional u_L at four cross-sections and the horizontal maps of depthintegrated u_L are shown in Fig.1 and Fig.2, respectively. The flow is directed toward the bay head (x = 0) in the deep central part and flows toward the bay mouth (x = L) over the shallow water and the slope. Down-welling occurs in the deep central part and upwelling occurs over the slope (Fig.1). For the depth-integrated flow, inflow occurs in the deep central part and outflow occurs in the shallow water (Fig.2). These characteristics are the same as those reported by Jiang and Feng (2014).

There are obvious improvements in the present result compared with that of Cui et al. (2019) (Fig.3 and Fig.4). The inflow in Fig.1 occupies a thinner width in the ydirection and reaches a deeper depth compared with the wider and shallower inflow displayed in Fig.3 in Cui et al. (2019). Fig.1 also shows less variation of vertical velocity at different crosssections compared with that shown in Fig.3 in Cui et al. (2019). Overall, the distribution of u_L along the four cross-sections (Fig.1) and depth-integrated u_L (Fig.2) from the present solution are closer to the analytical solution of Jiang and Feng (2014) than that shown in Cui et al. (2019) (Fig.3 and Fig.4).

4 Application to a Bay with Longitudinal Variation of Depth

Up to now, the governing equation of u_L can only be analytically solved in a narrow bay with lateral topography under the no-slip bottom boundary condition. In this section, the numerical solution of u_L is obtained without these limitations.

In the following, we apply the model to a wide bay with topography varying in the longitudinal direction according to the equation

$g(x) = 5 + 5{\rm{exp}}\{ {((x/180 -1000)/350)^2}\}, $

with an average depth of 6.5 m. The bay width is set as 54 km to represent the width-to-wavelength ratio δ > 0.1. The bay length is 180 km; hence, the corresponding nondimensional bay length is 0.5. The model grid has 400 grids in the x-direction and 100 in the y-direction and a grid spacing of 450 m by 540 m. The quadratic nonlinear bottom friction is adopted, i.e.

$\upsilon \frac{{\partial \left({{u_0}, {\kern 1pt} {\kern 1pt} {\kern 1pt} {v_0}} \right)}}{{\partial z}} = {C_d}\sqrt {u_0^2 + v_0^2} \left({{u_0}, \;{v_0}} \right){\kern 1pt} $

and

$\upsilon \frac{{\partial \left({{u_L}, {\kern 1pt} {\kern 1pt} {\kern 1pt} {v_L}} \right)}}{{\partial z}} = {C_d}\sqrt {u_L^2 + v_L^2} \left({{u_L}, \;{v_L}} \right).$

The bottom friction coefficient C_d is 0.0025. The vertical eddy viscosity coefficient v is set as 0.006 m² s⁻¹ to ensure $\beta = \nu \;{T_c}/h_c^2 = 1$. Hence, the vertical eddy viscosity term and the local acceleration term are both important in the governing equations of u_L. With the time-dependent method, the model is run with a time step of 12 s for 72000 time steps when the steady state of u_L is obtained.

4.1 Contribution of the Tidal Body Force to u_L and ζ_E

We now examine the contributions of different components of the tidal body force to u_L and ζ_E. First, Fig.5 and Fig.6 respectively compare the contributions of π₁ and π₂ with f = 0 (i.e., without considering the earth's rotation) in Eqs. (7), (8), (11) and (12). As can be seen, π₁ makes much more contribution than π₂ for both u_L and ζ_E. Although the topography of the bay in this study is different from that of Cui et al. (2019), the relative contributions of π₁ and π₂ are similar.

Next, we consider the solution with the Coriolis force not being zero. In this case, Fig.7 shows the distributions of u_L due to π, π₁ + π₂, and π₃. Note that the solution for π₁ + π₂ (Fig.7b) differs from that shown in Fig.5a. Set f ≠ 0 in the left hand side of Eqs. (7) and (8) induces asymmetry in the distribution of u_L relative to the center of the bay, in contrast to the case of f = 0. This contribution of the Coriolis force term to u_L is similar to that in Cui et al. (2019). Overall, π₃ makes a similar important contribution to u_L as π₁.

4.2 Influence of the Geometrical Shape on u_L

The geometry of the bay is characterized by two nondimensional numbers, δ = B/λ_c (bay width to tidal wavelength) and L = L^*/λ_c (bay length to tidal wavelength).

Fig.8 compares u_L with different values of δ. Overall, the distribution and intensity of u_L do not change significantly for these three cases with δ > 0.1. With increasing δ, the intensity of u_L is slightly reduced, and the width of the inflow extends in the bay. Moreover, the depth of inflow becomes shallower in the outer half of the bay.

The similarities in u_L for different values of δ can be attributed to the unchanged topographic gradient. Cui et al. (2019) showed that u_L is intensified significantly if the topographic gradient increases with different δ values. Without changing the topographic gradient, the different values of δ associated with changes in bay width only cause slight changes in u_L.

Fig.9 compares u_L with three values of L: 0.25, 0.5, and 1. The corresponding depth-integrated u_L is displayed in Fig.10. At x = L/4 near the bay head, u_L is similar for different values of L. The differences become more significant at the other three sections. Meanwhile, u_L becomes weaker as L increases. Three pairs of gyres appear for L = 0.25 and L = 0.5. The two pairs of gyres merge into one for L = 1.

With the presence of topography, the decreasing bay length results in the topographic uplift and the enhancement of the average incoming energy. Fig.11 shows u_L with different values of L, but with the topography set to be flat with a depth of 10 m. In this case, u_L decreases slightly with increasing L. This suggests that topography is more important than bay length on influencing u_L.

4.3 Influence of Incoming Tidal Strength on u_L

Fig.12 and Fig.13 display u_L for different tidal amplitudes (ζ_open) at the bay mouth: 50 cm, 100 cm, and 150 cm, respectively. As can be clearly seen, the magnitude of u_L increases with increasing ζ_open. With ζ_open increasing from 50 cm to 150 cm, the maximum of |u_L| increases by about nine times at the bay mouth and about three times at x = L/4. For the case of the depth-integrated u_L, three pairs of gyres can be found in the bay with ζ_open = 50 cm or 100 cm, and the two pairs of gyres near the bay head disappear with ζ_open = 150 cm. Furthermore, the region occupied by the inflow becomes larger with the increase of ζ_open.

The zeroth-order tidal velocity changes linearly with the tidal amplitude at the open boundary without considering the bottom friction dissipation. The tidal body force results from the nonlinear coupling of the zeroth-order tide. With ζ_open increased by three times from ζ_open = 50 cm to ζ_open = 150 cm, the magnitudes of the tidal body force and u_L should also increase by nine times. As the bottom friction also increases with the increasing tidal velocity, the magnitude of u_L increases by less than nine times. Furthermore, the intensification of u_L magnitude becomes lesser from the bay mouth to bay head due to the relative influence of bottom friction increasing toward the bay head.

4.4 Influence of the Momentum Dissipation on u_L

Momentum dissipation consists of two parts. The first part is the eddy viscosity, related mainly to the vertical eddy viscosity, and the second part is related to the bottom friction. The two parts are coupled together.

The vertical eddy viscosity coefficient affects u_L indirectly from changing the zeroth-order tide, but also directly by changing the intensity of the tidal body force and the vertical momentum transfer in the governing equations of u_L. Fig.14 shows u_L obtained using three values of the vertical eddy viscosity (υ): 0.002 m² s⁻¹, 0.004 m² s⁻¹, and 0.006 m² s⁻¹. Both the magnitudes and patterns of u_L are changed with different values of υ. The magnitude of u_L increases by less than two times as υ increases from 0.002 m² s⁻¹ to 0.006 m² s⁻¹. With the increasing υ, } the magnitude of π₁ increases. At the same time, the vertical transfer of u_L also increases. This results in smaller increases in the magnitude of u_L relative to the increase in υ. As υ decreases, the vertical transfer of u_L also gets smaller. This can explain the subtle differences in the vertical distribution of u_L shown in Fig.14.

Fig.15 and Fig.16 present the distributions of u_L obtained using four values of quadratic bottom friction coefficient, i.e., C_d = 0.00125, 0.0025, 0.005, and 0.01. As the value of C_d increases, the magnitude of u_L decreases and the flow pattern becomes simpler. The changes are significant in the inner half than in the outer half of the bay. The changes are minimal near the bay mouth and manifold in inner zone of the bay. With C_d increasing from 0.00125 to 0.01, the maximum value of u_L decreases by 30% at x = 3L/4 and by 60% at x = L/4. Obvious changes in the patterns of u_L with different values of C_d can also be observed. In particular, with the increasing C_d, the near bottom inflow extends further into the interior of the bay.

5 Conclusions

In this study, an improved computing method for the tidal body force is established to obtain a more accurate 3D Lagrangian residual current u_L. First, the solution of u_L obtained with the new method for an elongated rectangular bay shows better agreement with the analytic solution of Jiang and Feng (2014) and is more accurate than the results of Cui et al. (2019) using the tidalaveraging method. Next, the new method is applied to a wide bay with longitudinal topography. The impacts of various geometrical and physical parameters on the solution of u_L are then analyzed.

For the case of the tidal body force without including the Coriolis effect, it is clear that the intensity and pattern of u_L are governed by π₁, and π₂ has a much smaller contribution. This is because π₂ is one order of magnitude smaller than π₁. The rotation of the earth has two effects: the rotation effect within the residual current equations itself, and the rotation effect in π₃ transferred from the zeroth-order tide. The former one creates asymmetry of u_L along the center line of the bay. The latter one has the same order of magnitude as π₁. Therefore, π₃ and π₁ jointly determine the solution of u_L.

Different geometrical parameters of the bay have varying impacts on u_L. The changes of δ (i.e., ratio of bay width to tidal wavelength) do not change the topographic gradient of the bay. Therefore, the solution of u_L are similar for different values of δ. The uplift of the longitudinal topography with the decrease of L (ratio of bay length to tidal wavelength) leads to an increase of u_L and different circulation patterns. The enhancement of the residual ve locity can be attributed to the larger topographic uplift. Hence, topographic uplift has significant impacts on the Lagrangian residual current.

The incoming tidal strength can significantly affect the intensity of u_L. The magnitude of u_L increases nearly as much as nine times at the bay mouth and only three times at the bay head when the incoming tidal strength increases by three times. These results are attributed to the quadratic increase of the tidal body force and the enhancement of the bottom friction dissipation.

Two kinds of the dissipation effect in the bay affect u_L by different mechanisms. The intensity of the tidal body force rises by about three times when the vertical eddy viscosity coefficient υ also increases by three times, and the intensity of u_L increases by less than two times. The dissipation effect of the bottom friction and the eddy viscosity are enhanced correspondingly as the tidal body force rises by three times. This results in less increase of u_L than that of the tidal body force.

The intensity of u_L reduces with the increase of the bottom friction. The decrease of u_L with C_d increasing from 0.00125 to 0.01 is more signified near the bay mouth than near the bay head. This is because the nonlinear effect of the bottom friction dissipation enhances from the bay mouth toward head. The enhancement of the nonlinear bottom friction dissipation inhibits the seawater from flowing into the bay, delays the inflow area reaching the bottom, and finally simplifies the flow state of u_L.

Acknowledgements

This study was supported by the National Natural Science Foundation of China (No. 41676003) and the NSFC Shandong Joint Fund for Marine Science Research Centers (No. U1606402).