Oblique sand ridges in confined tidal channels due to Coriolis and frictional torques

The role of the Coriolis effect in the initial formation of bottom patterns in a tidal channel is studied by means of a linear stability analysis. The key finding is that the mechanism generating oblique tidal sand ridges on the continental shelf is also present in confined tidal channels. As a result, the Coriolis effect causes the fastest growing pattern to be a combination of tidal bars and oblique tidal sand ridges. Similar as on the continental shelf, the Coriolis-induced torques cause anticyclonic residual circulations around the ridges, which lead to the accumulation of sand above the ridges. Furthermore, an asymptotic analysis indicates that the maximum growth rate of the bottom perturbation is slightly increased by the Coriolis effect, while its preferred wavelength is hardly influenced.


Introduction
Tidal bars are rhythmic bottom patterns that occur in many tidal channels (e.g., the Western Scheldt in the Netherlands, the Exe Estuary in England, the Ord River Estuary in Australia, and the Venice Lagoon in Italy). These bars are several meters high and have wavelengths of 1-15 km. Their characteristics are determined by channel properties (depth, width, tidal amplitude, etc.), which may change due to, for example, dredging, sea level rise, and land reclamation. Tidal bars are invaluable for many organisms that feed on their rich grounds, but they also may hamper marine traffic. For proper management of tidal channels, it is therefore important to understand their behavior. Seminara and Tubino (2001), Schramkowski et al. (2002), and Hepkema et al. (2019), among others, studied the physical mechanism that causes tidal bars to form, as well as the sensitivity of their wavelength to channel properties. They explained that the initial formation of tidal bars can be understood by analyzing the residual currents generated by the topography (using arguments similar to those by Zimmerman (1981)). Hibma et al. (2004) showed that the results of the linear stability analysis of Schramkowski et al. (2002) compare well with results of a numerical morphodynamic model, Delft3D.
In these linear stability studies and in the study by Hibma et al. (2004), the Coriolis effect was neglected. However, tidal bars occur in natural systems (e.g., Western Scheldt) where the Coriolis force is a first-order term in the momentum balance. The importance of Coriolis on the hydro-morphodynamics in tidal channels is supported by several other studies, e.g., Valle-Levinson (2008), Winant (2008), Xie et al. (2017), and Olabarrieta et al. (2018). Furthermore, 2D morphological simulations similar (but now with and without the Coriolis effect) with those performed by Hibma et al. (2004) show clear differences between the initial formation of bottom patterns with and without the Coriolis effect taken into account (see Fig. 1 and supplementary information). In the left panel of Fig. 1, perturbations have a braided tidal bar pattern, whereas a ridge-like pattern emerges when the Coriolis effect is neglected (right panel).
This motivates an investigation into the physical mechanisms that explain how the initial formation of bottom patterns in a tidal channel is affected by the Coriolis effect. To this end, a semi-analytical model is developed in Section 2 and investigated by means of a linear stability analysis. In Section 3, the results of the linear stability analysis are presented and further analyzed with an asymptotic expansion considering a weak Coriolis force. Section 3 is followed by a discussion (Section 4) and the conclusions (Section 5).

Governing equations
The semi-analytical model developed in this section is similar to those used by Schramkowski et al. (2002) and Hepkema et al. (2019), which successfully explained the emergence of tidal bars in confined channels. The main extension here is that the Coriolis effect is taken into account.
The domain shape consists of an open section of a straight channel distant from the seaward and landward boundaries. The channel has a uniform width B * , depth H * , and a length which is small compared to the tidal wavelength and the length scale of channel width variations (Fig. 2). These choices imply that we consider a local model, ignoring sloping background topography, and apply the rigid-lid approximation (sea surface elevation only appears in the pressure gradient force). The hydrodynamics is governed by the depth-averaged shallow water equations, including the Coriolis effect. The flow is driven by a spatially uniform pressure gradient that oscillates with the principal tidal frequency σ * . The dimensional (henceforth denoted with an asterisk) continuity and momentum equations read: Here, t * and x * = (x * , y * ) are the time and space coordinate. Furthermore, ∇ * = (∂/∂x * , ∂/∂y * ), u * = (u * , v * ) is the depth-averaged current velocity, ζ * the free surface elevation, and h * the bottom elevation (with respect to the undisturbed bed). In Eq. 2, the bottom stress is linearized (see, e.g., Zimmerman 1982;Terra et al. 2005), i.e., it is modelled as ρ * r * u * , where ρ * is the water density and r * = c d U * 8/(3π), with c d the drag coefficient, and U * a typical amplitude of the tidal current along the channel direction. The parameter g * is the gravitational acceleration and with f * the Coriolis parameter (here assumed constant). The bed elevation h * evolves due to the divergence of volumetric sediment transport q * : Here, H * is the undisturbed depth, ζ * the free surface elevation, h * the bottom height, and B * the channel width. The x * , y * , and z * arrows denote the direction of the coordinate axes (x * is along-channel and y * cross-channel) Here, p is a porosity parameter and where u * 2 = u * 2 + v * 2 and s * 1 , s * 2 , b 1 , and b 2 are positive real numbers. The first term in Eq. 4 represents the advective transport of sediment and the second term accounts for bed slope effects. Equation 4 corresponds, with different choices for s * 1 , s * 2 , b 1 , and b 2 , to most bed load and total load sediment transport formulations (Soulsby 1997).
The boundary conditions imposed at the sides of the channel are: v * = 0 and ∂h * ∂y * = 0 at y * = 0, B * .
Equations 1-4 are made dimensionless (no asterisk) by first scaling time t * and space x * as: Subsequently, the depth-averaged velocity u * = (u * , v * ), the free surface elevation ζ * , and the bottom elevation h * are scaled as: The dimensionless continuity and momentum equations read: where ∇ = (∂/∂x, ∂/∂y), σ = σ * B * /U * , r = r * B * /(U * H * ) and The dimensionless bed elevation equation reads is a bed slope parameter and ε the ratio of the tidal time scale 1/σ * and the morphological time scale (1 − p )H * B * /Q * , with Q * = s * 1 U * b 1 +1 a typical volumetric sediment transport magnitude. Given that the bed evolves slowly and that the sediment transport varies almost periodically (with the tidal frequency), the bed evolution is approximated by the tidal average of the divergence of the sediment transport: where τ = εt and · denotes the average over one dimensionless tidal period 2π .

Linear stability analysis
The bottom pattern that initially forms when a tidal current flows over a horizontal sandy bed is analyzed with a linear stability analysis. An equilibrium solution (u 0 , v 0 , ζ 0 , h 0 ) to the system of Eqs. 6-9 is described by a spatially uniform symmetrical tidal current u 0 = (u 0 , v 0 ) = (cos(t), 0), driven by a spatially uniform pressure gradient: Small perturbations on the flat bed result in perturbations of the flow variables and vice versa. Let ξ be small and substitute: in Eqs. 6-8. At O(ξ ), this results in a linear system of partial differential equations: Here, Λ =Λ |u 0 | b 2 is a bed slope parameter and u 1 and v 1 are the components of u 1 = (u 1 , v 1 ). In the derivation of Eq. 13, the continuity (10) and |u 0 | b 1 u 0 = 0 are used to simplify the expression. Combining the continuity and momentum equations (substituting (10) in ∂/∂x (12)−∂/∂y (11)) eliminates ζ 1 and yields an equation for the vorticity Ω 1 : with The first term on the right-hand side of Eq. 14 is the torque due to bottom friction. The second term on the right-hand side is the Coriolis torque. The terms on the left-hand side of Eq. 14 describe the vorticity's inertia, its advection by the equilibrium current u 0 , and its dissipation due to bottom friction, respectively.
Equations 10, 13, 14, and 15, together with the boundary conditions, allow for solutions of the form: Here, c.c. stands for complex conjugate,û(t, y, k), v (t, y, k),Ω(t, y, k),ĥ(τ, y, k) are complex valued functions, and k is a dimensionless along-channel wavenumber, which relates to a dimensional wavenumber k * = k/B * . Substituting (16) in the continuity (10) and the vorticity (14) yields: The parameter γ = (rRo) −1 denotes the relative importance of the vorticity producing torque due to the Coriolis effect and the torque due to the bottom friction. Thus, when the Coriolis effect is neglected, γ = 0 and only the torque due to bottom friction produces vorticity.
Substituting (16) in the bottom evolution (13) results in: Inspired by the boundary conditions and the fact that without the Coriolis effect,ĥ ∼ cos(nπy) with n a natural number (Schramkowski et al. 2002), we writeĥ as a cosine series: withh n (τ, k) complex valued functions for every natural number n. In the Appendix, the currentû is expressed in terms of the bottom elevationĥ by solving (17)-(18). Furthermore, it is shown in the Appendix that substituting the cosine series (20) into the bottom evolution (19) and truncating the summation at N results in: Here,h = (h 0 , . . . ,h N ), D = diag(ω 0 0 , . . . , ω 0 N ) is a diagonal matrix and A is a matrix with elements a mn , where admits solutions of the form: with ω an eigenvalue of D + γ A and h the corresponding eigenvector. The superscript 0 of the elements in the diagonal matrix D represents the fact that these are the eigenvalues when γ = 0 (no Coriolis effect). For every wavenumber k, we calculate the eigenvector h j (k) corresponding to the j -th eigenvalue ω j (k). The eigenvalues and corresponding eigenvectors are sorted such that ω 0 ≥ · · · ≥ ω N . The wavenumber k, for which the largest growth rate ω 0 (k) is attained, is called the preferred wavenumber k pref . The corresponding eigenvector h 0 (k pref ) = (p 0 , . . . , p N ) sets the spatial structure of the fastest growing bottom pattern: The second and third largest eigenvalues are denoted by ω 1 and ω 2 , respectively. The dimensionless growth rates ω relate to dimensional ones by ω * = ωεσ * (since τ = εt = εσ * t * ).

Sensitivity to Rossby number
To investigate the role of the Coriolis effect on the initial formation of bars, the inverse Rossby number Ro −1 is varied. All other parameters are kept fixed such that varying Ro −1 corresponds to varying the Coriolis parameter f * . The results are summarized in Fig. 3. The top panels of the figure show the first three growth curves (i.e., ω 0 , ω 1 and ω 2 versus wavenumber k) for different values of the inverse Rossby number Ro −1 . The figure reveals that the growth rate of the fastest growing pattern slightly increases with increasing Ro −1 . The wavenumber k pref for which the growth rate attains its maximum value is hardly influenced by the Coriolis effect. The panels in the two middle rows in Fig. 3 illustrate the patterns and the residual current that correspond to the fastest growing patterns for different values of the inverse Rossby number. When the Coriolis effect is neglected (Ro −1 = 0), the patterns have a tidal bar structure of cos(nπy) cos(kx), with n a natural number and k the wavenumber. However, the pattern significantly changes when the Coriolis effect is considered. In that case, the fastest growing pattern seems a combination of tidal bars and oblique tidal ridges. Moreover, when the Coriolis effect is neglected or very weak, the cells of residual current are in between troughs and bars, whereas the residual currents go around the bars and troughs when the Coriolis effect is taken into account.
The bottom row of panels in Fig. 3 shows, for different values of the inverse Rossby number Ro −1 , the components p n of the eigenvector h 0 (k pref ). The figure reveals that the spectrum widens with increasing Ro −1 . Moreover, looking amplitudes and arguments of cosine series coefficients that correspond to the the bottom patterns in the second and third rows (see Eq. 25). A value Ro −1 = 0.5 is representative for the Western Scheldt. A dimensionless growth rate ω = 0.5 corresponds to a dimensional growth rate of approximately 0.16 yr −1 and a dimensionless wavelength 2π/k = 3 corresponds to a dimensional wavelength of 15 km at the phases of the components p n , it appears that, in all cases considered, the patterns are a good approximation of the form: h pref ≈ 2|p 1 | cos(πy) sin(k pref x) +2|p 2 | cos(2πy) cos(k pref x) −2|p 3 | cos(3πy) sin(k pref x).
The fact that only a narrow part of the spectrum is involved follows from expression (23). Since the coefficients a mn decrease when (m 2 − n 2 ) increases, the equation of ∂h m /∂τ is mostly dependent on the coefficientsh n for which n is close to m.

Asymptotic analysis for weak Coriolis force
Above, it was shown that when the Coriolis effect is considered, the maximum growth rate slightly increases and the fastest growing pattern changes significantly. An explanation for this is sought by analyzing the system for small values of γ (while fixing r), which corresponds to small values of the inverse Rossby number Ro −1 (i.e., weak Coriolis force). Let, as before, ω j and h j for j = 0, . . . N be eigenvalues and eigenvectors of the (perturbed) eigenvalue problem: and expand these eigenvalues and eigenvectors in powers of γ : , for j = 0, . . . , N. The eigenvectors of the O(1) eigenvalue problem are h 0 j = e j with e j the standard basis vectors of R N+1 (i.e., the j -th entry of e j equals 1 and the others 0). The eigenvalues ω 0 j are given in Eq. 22. At O(γ ), Eq. 26 reads: The O(γ ) correction to the eigenvalues is now computed by taking the standard inner product of Eq. 27 with e j , resulting in: Here, it is used that ω 0 j h 1 j · e j = Dh 1 j · e j . Given that a jj = 0, it follows that there is no O(γ ) correction to the eigenvalues. However, there is an O(γ ) correction to the eigenvectors. To see this, take the inner product of Eq. 27 with e m (m = j ) to obtain the m-th component of h 1 j : Figure 4 shows that the the superposition of the tidal bar pattern (colors) and the perturbation induced by the Coriolis Lastly, we show that the correction to the largest eigenvalue is positive by considering the O(γ 2 ) problem: Taking the inner product with e j and substituting ω 1 j = 0 results in: where again, it is used that ω 0 j h 2 j · e j = Dh 2 j · e j . For all m, a mj a jm > 0 and if ω 0 j is the largest eigenvalue of the unperturbed system, ω 0 j − ω 0 m > 0. Hence, a weak Coriolis effect increases the maximum growth rate of the bottom perturbations, which is consistent with the findings in Fig. 3.

Oblique tidal ridges versus tidal bars
The differences in patterns when the Coriolis effect is taken into account or not follow from the additional torque exerted by the Coriolis effect on water motion over the longitudinally sloping bed. When the Coriolis effect is neglected and only the frictional torque due to lateral bottom slopes is considered, the fastest growing bottom patterns consist of tidal bar patterns, identical to those obtained by Seminara and Tubino (2001), among others. An elaborate explanation of the physical mechanism of the formation of tidal bars is given by Hepkema et al. (2019). When the Coriolis torque is taken into account, the joint action of the two torques results in the formation of oblique tidal sand ridges. This mechanism is the same as the one responsible for the initial formation of offshore tidal sand ridges, as explained by Huthnance (1982). For completeness, we summarize the two mechanisms below.
Both the morphodynamic instability leading to tidal bars and the instability leading to tidal sand ridges result from the fact that perturbations of the flat bottom in the form of tidal bars or tidal sand ridges alter the tidal currents flowing over them (or vice versa). This tide-topography interaction results in residual currents, such that both during the ebb and flood phase, the currents become stronger upstream of the bars/ridges, whereas they are weakened downstream of the bars/ridges (and vice versa for the troughs) (Zimmerman 1981). This results in sediment transport converging at the bars or ridges and diverging at the troughs, hence the instability. Figure 5 illustrates the tide-topography interaction. In the top panels, a tidal bar pattern is considered with a tidal current flowing over it. Due to lateral gradients in the bottom, a frictional torque generates vorticity at the principal tidal frequency, between the bars and troughs (green round arrows in the figure). The M 2 vorticity is transported by the background M 2 current, resulting in a residual tidal vorticity as indicated by the light blue arrows in the figure. Adding the background current to the residual current results in higher velocities upstream of a bar and lower velocities downstream of the bar (and vice versa for the troughs).
In the bottom panels of Fig. 5, a tidal sand ridge pattern is shown with the same background current. Also here, due to lateral bottom slopes, a frictional torque generates vorticity at the M 2 frequency. Now, in addition, the longitudinal bottom slopes result in a Coriolis torque. On the Northern Hemisphere, the Coriolis and frictional torque are in the same direction when the ridges are rotated anti-clockwise with respect to the background current. Also, in this case, the background M 2 current transports the M 2 vorticity, resulting in residual vorticity as indicated by the light blue arrows in the figure. Again, when the residual current is added to the background current, the velocities are higher upstream of the ridge than downstream and vice versa for the troughs.
An essential difference between the residual currents due to tidal bars and those due to tidal ridges is that in the later case the residual currents are around the crests, whereas in the former case they are directed toward the crests. Note that pure tidal ridges as in the figure can not form in a confined channel because they violate the lateral boundary conditions, but a similar pattern is possible.

Sediment transport formulation
In the sediment transport formulation q in Eq. 8, we chose b 1 = b 2 = 2. This corresponds to advective bed load transport as in Bailard (1981) and a bed slope effect due to eddy diffusivity, which is larger than the bed load bed slope effect as in the Bailard formulation (Hepkema et al. 2019). Following the later study, the bed slope parameter s * 2 is calculated as: where μ * (≈ 5 m 2 s −1 ) is the horizontal eddy diffusion coefficient, w * s (≈ 0.013 m s −1 ) the sediment settling velocity, κ * v (≈ 0.01 m 2 s −1 ) the vertical eddy diffusivity coefficient, α * (≈ 4 · 10 −6 s m −1 ) an erosion parameter, and γ * (≈ 0.017 s −1 ) a deposition parameter.
When b 1 = 2 and b 2 = 3, the transport is the same as for the Bailard formulation for bed load sediment transport (Bailard 1981). When b 1 = 4, advective transport is similar to the Engelund Hansen parametrization for total load transport (Engelund and Hansen 1967). The magnitude of the dimensionless growth rate depends on the choice of b 1 and b 2 (see Fig. 6), but the shape of the fastest growing patterns is hardly affected by this choice (see supplementary information). When b 1 = 4, the growth rates increase. When b 2 = 3, the differences are minimal for k ≈ k pref . For larger wavenumbers, the growth rate is smaller when b 2 = 3 compared with b 2 = 2, due to an increased bed slope effect. However, for those wavenumbers, advection of suspended sediment becomes relevant and the total load sediment transport formulation (4) is no longer valid (Hepkema et al. 2019).

Comparison with a numerical model
The results of Section 3, obtained with the semi-analytical model, correspond qualitatively with those of the numerical model that has a similar setup as the one from Hibma et al. (2004), but includes the Coriolis effect (see SI). When the Coriolis effect is taken into account, oblique sand ridges also form in this numerical model. Moreover, their growth rate is larger than that of the tidal bars that form in the absence of Coriolis (see the color bars in Fig. 1), while their wavelength is hardly affected.
The wavelengths of the tidal bars in Fig. 1 produced with the numerical model are approximately 5-7 km. The average depth in the area considered is 6.5 m and the current velocity amplitudes are close to 1.0 m s −1 . In the semianalytical model, when H * = 6.5 m and U * = 1 m s −1 , the preferred wavenumber is k pref ≈ 3.6, corresponding to dimensional wavelengths of approximately 8.7 km. Hence, the wavelength of the semi-analytical model is slightly larger, but in the same order of magnitude as those simulated by the numerical model. The similarities between the two different models confirm robustness of the discussed mechanism.

Conclusions
To study the role of the Coriolis effect in the initial formation of bottom patterns in a tidal channel, a semianalytical model is extended to include the Coriolis effect. It was shown that the Coriolis effect breaks the (anti-) symmetry of the bottom pattern that initially forms. The fastest growing bed perturbation can be characterized as a combination of tidal bars and oblique ridges, unlike the case without Coriolis effect where only tidal bars form. The Coriolis effect also modifies the residual current so that it drives anticyclonic circulations around the ridges. The mechanism behind these modifications is the same as the one causing the formation of oblique tidal sand ridges on the continental shelf. Compared with the case where the Coriolis effect is neglected, the preferred wavenumber is similar, while the maximum growth rate slightly increases.
Funding Direct funding from Utrecht University.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creativecommonshorg/licenses/by/4.0/.
The inner products of P 1 n and P 2 n with 2 cos(mπy) read: Therefore, substituting the cosine series in the bed evolution equation and taking the inner product with 2 cos(mπy) yield: with ω 0 m and a mn as in Eqs. 22 and 23. Truncating the sums at N results in Eq. 21.