Sensitivity of Tidal Bar Wavelength to Channel Width

Tidal bars are repetitive estuarine bedforms with heights of several meters and wavelengths in the order 1–15 km. Understanding their formation and sensitivity to changes in channel characteristics is important as they hamper marine traffic and play a crucial role in the local ecosystem. Recent observations suggest that the local width of the channel is dominant in determining the tidal bar wavelength. However, theoretical studies could not reproduce the sensitivity of the tidal bar wavelength to channel width. This discrepancy between theory and observations suggests that a mechanism is missing. In this study, one of the theoretical models is extended and results in tidal bar wavelengths, lateral mode numbers, and growth rates that agree fairly well with those of natural tidal bars, including the wavelengths dependency on estuary width. An important extension of the model concerns the bed slope induced diffusive suspended load transport of sediments. With this, it is explained why previously, the modelled tidal bar wavelengths depend only weakly on estuary width and why in the extended model it does. This has, from a modelling point of view, general implications for morphological models using a total load sediment transport formulation with a so‐called bed slope parameter.


Introduction
Many tidal channels have a rhythmic bottom topography consisting of so-called tidal bars. Examples are the Ord River estuary in Australia, the Exe estuary in England, the Netarts bay in the United States, and the Western Scheldt in the Netherlands (see Figure 1 and Leuven et al., 2016 for many more examples). These bars have wavelengths of 1-15 km and heights of several meters. As tidal bars hamper marine traffic and provide rich feeding grounds for many living organisms, it is important to understand their formation and sensitivity to tidal channel characteristics (i.e., channel width, depth, and currents). These characteristics may change due to, for example, dredging, sea level rise, or land reclamation. The sensitivity of the tidal bars to changes in channel characteristics was previously investigated by measurements, laboratory experiments, idealized models, and complex numerical models.
By analyzing data of 25 tidal bars in creeks in South Carolina and in the Salmon River estuary, Dalrymple and Rhodes (1995) found that the tidal bar wavelength is typically six times the channel width. Also, Leuven et al. (2016) found that the tidal bar wavelength correlates the strongest with channel width, compared to other channel characteristics such as depth and current velocity. They came to this conclusion by measuring the wavelength of 190 tidal bars in Google Earth. Furthermore, Tambroni et al. (2005) investigated morphodynamic equilibria in tidal channels and inlets in a laboratory. Their tidal bars formed after 50 tidal cycles with wavelengths of approximately three times the channel width.
To gain fundamental understanding about tidal bars, Seminara and Tubino (2001) analyzed a 3D-idealized model. In their model, the section of the channel, in which the tidal bars are studied, is assumed to be short with respect to the tidal wavelength, width variations, and depth variations. This setting differs from, for example, Schuttelaars and de Swart (1999) and ter Brake and Schuttelaars (2011), where the formation of tidal bars in a short semienclosed channel was investigated. Seminara and Tubino (2001) demonstrated that tidal bars can form as a free instability of an equilibrium state that describes a horizontally uniform tidal current over a horizontal bed. They found that the wavelengths of their modelled bars are comparable with those of observed bars. It further follows from their results that the ratio of channel width over wavelength of the bars scales close to linearly with the width-to-depth ratio. This implies that for a constant depth, the wavelength hardly depends on channel width. Schramkowski et al. (2002) showed that a depth-averaged model is able to model characteristics of tidal bars and that a 3D model is not needed. They included the local inertia terms in the momentum balance and showed that the current velocity and the water depth are also important in determining the tidal bar wavelength. This model also yields a tidal bar wavelength that only depends weakly on channel width (Leuven et al., 2016). Garotta et al. (2006) confirmed that a 2D model is sufficient and moreover showed that forcing by a tidal current that is either ebb or flood dominant leads to a net migration of bars.
Thus, the idealized models are not able to reproduce the observed sensitivity of tidal bar length on channel width. This could be a consequence of the many assumptions made in these models. To clarify this, numerical models, which contain less constraints, are helpful. Such a study was carried out by Hibma et al. (2004). Their numerical experiments were designed such that comparison with the idealized studies was possible. They found tidal bar wavelengths that are comparable with those found by Schramkowski et al. (2002). However, despite expectations after the work of Schramkowski et al. (2002), the wavelength was, besides sensitive to current velocity and channel depth, also sensitive to channel width. Yet they did not investigate this further. Also, van der Wegen and Roelvink (2008) found that the wavelength depends on both the estuary width, depth, and the current velocity. In addition, van de Lageweg and Feldman (2018) found that, in a strongly convergent estuary, an increase in basin width results in longer bars.
The considerations above show in particular that the idealized models are not able to explain the observed (and numerically reproduced) dependence of tidal bar length on channel width. A likely explanation for this discrepancy is that one or more essential physical processes are missing in these models. Compared with numerical model formulations, it turns out that the idealized models do not account for the effects of turbulent exchange of momentum and sediment in the horizontal plane. The intensities of these processes are measured by a horizontal eddy viscosity coefficient and a horizontal eddy diffusion coefficient, respectively. Therefore, in the present study, the model of Schramkowski et al. (2002) is extended with terms that describe horizontal turbulent exchange processes. Moreover, a critical current velocity for sediment erosion is added to achieve consistency with all other models.
The first aim of this work is to show that the extended model can mimic the pattern of observed tidal bars. The second aim is to show that the model can reproduce the sensitivity of tidal bar wavelength to channel width. The first and second aim are treated in section 3. Finally, we aim at explaining why the incorporation of the new processes results in a dependence of the tidal bar wavelength to channel width (section 4).

Governing Equations
The formation of tidal bars is modelled in a section of a tidal channel which is short compared to both the tidal wavelength and the e-folding length scale of channel width convergence. The model domain under consideration therefore consists of a 2D rectangular channel of fixed width B, mean depth H, and infinite length ( Figure 2). Note that by channel width B, the total width of the channel is meant (see also Figure 1). The hydrodynamics is modelled by the depth-averaged shallow water equations including horizontal eddy viscosity and assuming the rigid-lid approximation. The latter constitutes to only retaining the free surface elevation in the pressure gradient. Coriolis effects are not considered, and the bed shear stress and internal stresses are linearized. The depth-integrated sediment concentration is modelled by an advection-diffusion equation based on ter Brake and Schuttelaars (2010) and extended by including a critical velocity for erosion. The bed evolution (Exner) equation follows from the conservation of sediment mass. Since the convergence of sediment transport is almost periodic (with a period equal to the principal tidal constituent), the subtidal convergence is small compared to the instantaneous convergence of sediment transport. Therefore, the patterns evolve on a much larger time scale than the tidal period and is determined by the tidally averaged convergence of sediment transport (see Krol, 1991 for a discussion on the errors introduced by this approximation). The differences with the model of Schramkowski et al. (2002) are the critical velocity for erosion, the horizontal eddy viscosity, and horizontal eddy diffusivity. The resulting governing equations read with  the Heaviside function, which is zero (one) when its argument is (not) negative and || · || denotes the Euclidean norm of a vector. In these expressions, [s m −1 ] is an erosion parameter, u c [m s −1 ] the critical velocity for erosion, and [s −1 ] a deposition parameter. Lastly, q s [m 2 s −1 ] is the transport of suspended sediment and q b [m 2 s −1 ] the transport of sediment as bed load. In Schramkowski et al. (2002), it is argued that the bed load transport is much smaller than the suspended load transport but that the bed slope contribution to the bed load transport should be maintained. Here the suspended load transport is extended by a diffusivity term based on ter Brake and Schuttelaars (2010). The sediment transport formulations read whereŝ [s 2 m −1 ] is a bed load transport constant and ⋆ a bed slope parameter.
The bottom friction coefficient r = 8∕(3 )c d  is chosen to scale with a typical current velocity amplitude. The horizontal eddy diffusivity coefficients = c u h  B and = c c h  B scale with current velocity and the channel width. The vertical eddy diffusivity v = c v  H scales with current velocity and the undisturbed channel depth. Here c c h , c u h , c v are constants and  is a typically velocity, which will be specified in the next section. All parameters are summarized in Table 1. The bed load transport formulation and constantŝ s and ⋆ are taken from Schramkowski et al. (2002) and based on Bailard (1981). The value of the erosion parameter is also taken from Schramkowski et al. (2002) and based on Smith and McLean (1977), the deposition parameters and b are based on the derivation in ter Brake and Schuttelaars (2010), and the settling velocity w s corresponds to sediment with a median grain diameter d 50 = 0.15 mm. For the reasoning behind the values of the sediment parameters, we refer to the books of Dyer (1986) and Soulsby (1997).

The lateral boundary conditions imposed are
This means that there is neither transport of water nor sediment through the solid boundaries of the channel. Furthermore, all variables are assumed periodic in the longitudinal direction. Further details and underlying assumptions of the model are given in Schramkowski et al. (2002).

Linear Stability Analysis
Similar to the system in Schramkowski et al. (2002), the system (1)-(5) admits a spatially uniform equilibrium solution ((F p ) eq , u eq , C eq , h eq ), which is here chosen to be an M tidal flow 2 u eq = (u eq , v eq ) = (U cos( t), 0) over a flat bottom h eq = 0, with U [m s −1 ] the constant amplitude of the longitudinal equilibrium velocity u eq . The corresponding pressure gradient (F p ) = − eq g∇ eq and sediment concentration C eq are given in Appendix A. The typical current velocity scale  is chosen to be the amplitude of the equilibrium velocity U. Other harmonic compositions can be added to u eq , but are not considered here.
Tidal bars form as an instability on this spatially uniform equilibrium, ((F p ) eq , u eq , C eq , h eq ). Perturbations in the bottom h eq + h ′ result in perturbations in the flow u eq + u ′ , the pressure gradient −g(∇ eq + ∇ ′ ), and the sediment concentration C eq + C ′ . These perturbations, in turn, result in positive or negative feedbacks resulting in growth or decay of the bottom perturbations h ′ . The pattern with the largest growth overtakes the others and yields the length and time scale of the pattern that is formed initially (Dodd et al., 2003). (1)-(4) and linearizing the resulting equations, yields a system of linear partial differential equations in the primed variables, which are perturbations on the equilibrium variables. The linearized equations are given in Appendix Table 1 Parameters and Their Value (Range) After ter Brake and Schuttelaars (2010) and Schramkowski et al. (2002) Parameter Value(s) Units Name for every natural number n and real number k > 0. Here l n = n ∕B and̂, û,v,Ĉ, andĥ are complex valued functions of time. The number n is henceforth referred to as the mode number and is related to the number of bars and troughs in the lateral direction of the channel. For example, a pattern with mode number n = 1 corresponds to an alternating bar pattern and the pattern with n = 2 has one bar or a trough in the middle of the channel.
Since the perturbations in the bottom vary on a time scale that is much longer than the time scale on which perturbations of the other variables vary, the bottom is assumed to be constant during one tidal cycle. That is, h only depends on a long morphodynamic time scale whilê, û,v, depend on the short hydrodynamical time scale. This allows for calculating the fast variableŝ, û,v, andĈ, assuming a fixed bottom. The so-called "flow over topography problem" is solved during one tidal cycle.
Substituting equation (6) in the linearized continuity, momentum, and concentration equation for a constantĥ results in a system of ordinary differential equations for̂, û,v, andĈ. The solution to this system of equations is approximated by a truncated Fourier series. This Fourier series contains the residual component and overtides with frequencies p , with p an integer. Because the equations are linear and inhomogeneous, the resultinĝ, û,v, andĈ will be linear in the amplitude of the bed perturbationsĥ. Taking this into account, the amplitudes are approximated by choosing a large enough natural number N and substituting in the ordinary differential equations. Herẽp,ũ p ,ṽ p , andC p for p = −N, … , N are complex amplitudes, independent ofĥ. This results in a system of 4(2N + 1) algebraic equations iñp,ũ p ,ṽ p , andC p , which can be solved numerically. Recall that, since the resulting velocity and concentration field is solved for a fixed h ′ , the coefficients̃p,ũ p ,ṽ p , andC p depend on n and k. That is to say, the obtained ( ′ , u ′ , v ′ , C ′ ) describe how the flow, the pressure gradient, and the concentration react to the bottom perturbations h ′ = Re The next step is to determine the bed evolution. The linearized Exner equation reads with where Λ =ŝ ⋆ ⟨|u 3 eq |⟩. Since ′ , u ′ , v ′ , and C ′ are explicitly calculated in terms of h ′ , substituting these expressions in equation (8) results in a single equation for h ′ and hence forĥ, with complex growth rate in which the variables with the tildes are defined in equation (7) andC eq,0 in Appendix A. The real part of the complex growth rate represents the actual growth rate of the bottom perturbation and −Im{ }∕k its migration speed. Here turns out to be real because u eq consists of one tidal constituent.
The preferred wavenumber and mode number k pref and n pref are defined as those for which the growth rate is maximal, The preferred growth rate pref = (k pref , n pref ) and the preferred wavelength pref = 2 ∕k pref is the wavelength that can be compared to the distance between natural tidal bars. The preferred wavenumber, mode number, and growth rate are found using the optimization method "Brent" (Press et al., 2007).

Verification
The first objective is to verify that the model produces preferred wavelengths pref , mode numbers n pref , and e-folding growth time 1∕ pref that resemble those of bars that are observed in nature. For this, we selected four natural tidal channels. These four are chosen because they differ in width, depth, and current velocities and roughly satisfy the model assumptions. Table 2 summarizes the geometric and current characteristics of these channels, as well as the observed wavelength M pref and mode number n M pref of the tidal bars. The wavelength, mode number, and channel width are obtained from Google Earth (Figure 1). This is done by identifying the pattern of the form cos(n ∕B) cos(kx) (as in equation (6)), which represent the observed channels in Google Earth most accurately. The wavelength M pref is measured as the distance between successive crests. The uncertainty in the wavelength results from the uncertainty of the location of the crests of the bars. An alternating bar pattern (n M pref = 1) represents the patterns in the Exe estuary and the Netarts bay best. The patterns in the Ord and the Western Scheldt are described by patterns with a mode number close to three. The width of the pattern that fits best is as illustrated in Figure 1. The currents and undisturbed depths of the Western Scheldt, Ord River estuary, Exe estuary, and Netarts bay are based on Cancino and Neves (1999), Wolanski et al. (2001), Herrmann and Hübner (1982), and Glanzmann et al. (1971), respectively. Furthermore, these papers report that the considered channels are vertically well mixed (i.e., small vertical salinity gradients) and the tides predominantly semidiurnal with free surface elevation amplitudes of 1-2 m. The sediment properties are assumed to be similar in all channels (following Leuven et al., 2016). Typical time scales at which these bars evolve are of the order of decades.
The left column of Figure 3 shows the growth curves (i.e., graphs of (k, n)) for the four different channels. The red vertical line represents the measured wavenumber and the red area the measurement uncertainty.
The right column shows the bottom pattern corresponding to the preferred wavenumber and mode number on which the residual current is shown with arrows. The figure reveals that, using realistic parameter values, the model yields wavelengths of tidal bars that are of the same order of magnitude as those of bars observed in nature. Also, the preferred mode numbers n pref are in the right order of magnitude. They agree for the Exe estuary and the Netarts bay and for the Western Scheldt and the Ord River estuary the difference is one. Note that for the Ord River estuary, the maximal growth rates for the different mode numbers n = 1, n = 2, and n = 3 are similar. Lastly, the e-folding growth time is in the order of decades.

Sensitivity of Tidal Bar Wavelength to Channel Width
The second objective is to show that the extended model reproduces the sensitivity of the preferred wavelength to channel width. The lower lines in the top panels of Figure 4 (indicated by c c h = 0) show the relation between the preferred wavelength pref and the channel width B for the four channels when the horizontal eddy diffusivity , horizontal eddy viscosity , and the critical velocity for erosion u c are all zero. This case corresponds to the original model of Schramkowski et al. (2002), where the preferred wavelength, as Leuven et al. (2016) noted, remains of the same order of magnitude over the whole range of channel widths. Small increases in channel width result in slight increases of the wavelength pref . However, when the width B is increased even more, a pattern with a larger mode number n but with a smaller pref has the largest growth rate. That is, the preferred mode number jumps to a higher one.
The upper lines in the top panels of Figure 4 show the preferred wavelength versus channel width for the same parameter values, but with a nonzero horizontal eddy diffusivity. They reveal that, for all four channels, when horizontal eddy diffusivity = c c h UB, is nonzero, the relative change of the preferred wavelength pref over the considered range of channel widths is larger than when the horizontal eddy diffusivity is neglected. That is, when ≠ 0, the preferred wavelength pref is more dependent on channel width B than when = 0. The jumps in mode number occur later than in the case when horizontal eddy diffusivity is neglected. In fact, the middle (and top) row shows that the preferred mode number n pref decreases when the eddy diffusivity is taken into account. That is, the addition of horizontal eddy diffusivity reduces the growth of patterns with large gradients (high mode numbers and wavenumbers).
The addition of eddy viscosity or a critical velocity of erosion does not strongly alter the relation between the preferred wavelength and channel width. This is shown in Figure 5 for the Western Scheldt case. The right panels show the preferred wavelength versus channel width when eddy viscosity is taken into account. The left panels show the same when the critical velocity for erosion is nonzero. The eddy diffusivity is neglected in the top panels while it is taken into account in the lower panels. The difference between the lines in the top panels is minimal compared to the lower line in the top left panel of Figure 4, where both the eddy viscosity and the critical velocity for erosion are zero.

Different Contributions to the Growth Rate
The third objective is to explain why the incorporation of the horizontal eddy diffusivity results in a dependence of the preferred wavelength pref on channel width B. To facilitate the subsequent discussion in which the third objective is met, we analyze the different terms of the growth rate in more depth. To simplify expressions, the critical velocity of erosion u c and the horizontal eddy viscosity are henceforth neglected, since in the previous section, it is shown they do not significantly alter the relation between channel width and tidal bar wavelength.
The growth rate consists of four terms, where the factor −1∕(1 −p) is understood to be part of the different terms. The first term, adv , results from advective sediment transport. The second and third terms, diff depth-int and diff bed slope , are due to diffusive sediment transport and the fourth term bl bedslope originates from the effect of bed slopes on bed load sediment transport. Figure 6 shows the contribution of the different mechanisms to the growth rate versus wavenumber k for parameter values representative for the Western Scheldt. Next, we discuss the physics of the different terms of the growth rate . Figure 5. Preferred wavelength pref versus channel width B. In the left panels, a critical erosion velocity is included, and in the right panels, eddy viscosity is included. In the top panels, the eddy diffusivity is neglected while the lower panels, it is not. The remaining parameters are as in Tables 1 and 2 for the Western Scheldt. The first term in equation (12) is the advective part of the growth rate adv and relates to convergence of tidally averaged advective sediment transport, that is, −⟨∇ · (u ′ C eq + u eq C ′ )⟩.
To analyze adv , consider the p-th Fourier mode of the linearized concentration equation (3) with u c = 0, The dominant balance in the M 2 concentration equation for k ≈ k pref is between erosion and deposition in all four considered cases. For the Western Scheldt, this is shown in Figure 7. Then, neglecting the M 4 components in the perturbed current u ′ (as they are small compared to the M 2 components), the M 2 components of the perturbed concentration readsC  Tables 1 and 2 for the Western Scheldt case and c u h = 0, u c = 0 m s −1 and n = n pref = 2. The blue curve (inertia) lies below the yellow line (diffusion). When horizontal eddy diffusivity is neglected (as in Schramkowski et al., 2002), the main balance remains between erosion and deposition for k ≈ k pref , but n pref = 3 and k pref is larger. Figure 8. Generation of tidal residual vorticity cells for bottom pattern with mode number n = 1. (a) Situation during flood (u eq > 0). (b) Situation during ebb (u eq < 0). Friction (thin red arrows) experienced by water column is larger on the crests than on the troughs and generates tidal vorticity (red circles). Both during ebb and flood, positive vorticity is transported into the dashed box and the negative vorticity is transported out, resulting in the build-up of positive vorticity (black circle).
Substituting the latter in adv shows that the contribution of the advective transport scales with the magnitude of the residual currentũ 0 : The physical mechanism can now be understood by analyzing the residual current u 0 via vorticity arguments that follow from Zimmerman (1981). For this, consider a longitudinal tidal current that moves over a bottom that consists of bars and troughs as illustrated in Figure 8 (here the n = 1 pattern is assumed). The lateral depth variations result in frictional torques that generate tidal vorticity, as is indicated by the red circles. This vorticity is subsequently transported by the unperturbed tidal flow. As a result, both during ebb and flood, there is an influx of vorticity of the same sign into a region between crests and troughs. Thus, residual vorticity builds up in these areas (black circles in the figure), which is balanced by dissipation due to bottom friction. The resulting residual currents turn out to be directed from the troughs to the crests in the longitudinal direction and from the crest to the troughs in the lateral direction (see also right column in Figure 3). To understand why this leads to growth of bars and deepening of troughs, consider the convergence of the tidally averaged advective sediment transport (see equations (8) and (9)). This reads, using that the equilibrium velocity u eq and the equilibrium sediment concentration C eq are spatially uniform, where for the first equality sign, the perturbed continuity equation is used and for the second one, the fact that ⟨u eq C eq ⟩ = 0. Equation (16) shows that convergence of sediments occurs when u eq and C ′ ∕ x are negatively correlated. From equation (14), it follows that the perturbed concentration reads Since u 0 is directed towards the bars, this implies that both during ebb and flood, the perturbed concentration C ′ is positive upstream from a bar and negative downstream. On the other hand, C ′ is negative upstream a trough and positive downstream of it. Hence, the convergence of advective transport causes bars to grow and troughs to deepen for patterns with wavenumber k and mode number n close to the preferred ones. Figure 6 reveals that for larger k, adv is negative. The reason is that for these patterns, the dominant balance in the concentration equation is not between erosion and deposition but that the advection term also becomes significant (Figure 7).
The second term in equation (12) is diff depth-int . It originates from the fact that when at one location, the depth-integrated concentration C is larger than at another, a transport occurs from high to low C. Figure 6 reveals that diff depth-int is slightly positive for small k and negative for larger k and close to zero for wavenumbers close to the preferred wavenumber k pref . That is, for patterns with large wavenumbers, the depth-integrated concentration is larger above the troughs than above the crests, while the opposite is the case for patterns with small wavenumbers. Therefore, the depth-integrated diffusive transport results in a small contribution of the growth of patterns with large wavenumbers and slightly causes decay of patterns with small wavenumbers. Figure 9. Consider a uniform depth-integrated concentration over a sloping bottom. On the left, the water column is shallower than on the right. Hence, the concentration C 3 (t, x, y, z) is larger on the left than on the right. As a result, sediments diffuse from left to right even though the depth-integrated concentrations are equal. In this case, it is the gradient in bed level that induces diffusive sediment transport instead of a gradient in depth-integrated sediment concentration.
The third term in equation (12), diff bed slope , relates to a suspended load bed slope effect. It arises in the derivation of the depth-integrated concentration equation (3) from the 3D concentration equation when the diffusion term is integrated over depth, where C 3 is the "3D concentration" depending on x, y and z, and c b the concentration at the reference level above the bottom z = h(x). When a balance is assumed between vertical turbulent mixing (with vertical eddy diffusivity v constant over the vertical) and downward settling, the concentration at the bottom reads Schuttelaars, 2010). The first term in equation (17) results in diff depth-int and the second in diff bed slope . The physics behind the second term becomes apparent when considering, for example, a sloping bed and a uniform depth-integrated concentration. Then, the first term, ∇·( ∇C), is zero. However, in the shallow region, the concentration C 3 is higher than in the deep region, leading to sediment transport from the shallow to the deep region ( Figure 9). The second term models the convergence of this transport. As the magnitude of sediment transport from crests to troughs scales with the slope, diff bed slope is negative and becomes increasingly more negative if bottom gradients increase.
The last term of equation (12), bl bed slope , relates to bottom slope effects on bed load transport: sediment moves easier downslope than upslope. Therefore contributes bl bed slope , like the third term, to convergence of sediments in the troughs and divergence of sediments on the crests. Hence, it dampens the growth of the pattern. That is, bl bed slope < 0 and is stronger negative for larger bottom slopes. In contrast to diff bed slope results this contribution from bed load transport processes. Figure 6 reveals that bl bed slope is much smaller than diff bed slope , which, based on the values from Table 1, directly follows from the fact that The above gives a physical argument for using a larger bed load bottom slope parameter Λ than suggested by literature (Baar et al., 2018), when a total load formula is used. It then also represents the suspended load bed slope effect when this is not explicitly taken into account.

Mechanism Behind Sensitivity of Tidal Bar Wavelength to Channel Width
Here we address the third objective by discussing three important aspects of Figure 4, which shows tidal bar A 1. In the range of channel widths with constant n pref , pref increases with increasing B. 2. While increasing B, the preferred wavelength pref decreases when n pref jumps to a higher mode number.
A wavelength pref versus channel width B.

10.1029/2019JF005032
3.The range of channel widths for which n pref is the same is wider when horizontal eddy diffusivity in the sediment mass balance is taken into account.
Together, this shows that the wavelength depends on channel width when horizontal eddy diffusivity is taken into account and why the wavelength remains of the same order of magnitude for a large range of channel widths when the horizontal eddy diffusivity is neglected.
To analyze the three aspects above, the growth rate is first approximated as Here the two bed slope effects diff bed slope and bl bed slope are combined into an effective bed slope term and diff depth-int is neglected.
The first two aspects (A1 and A2) follow from the fact that, for a fixed mode number, the growth rate is sensitive to changes in channel width. For a fixed mode number n, the wavelength̃for which̃is maximal is approximately (see Appendix C)̃( where Here  f is a friction length scale and  t the tidal excursion length. When n = n pref in equation (20),̃(n pref ) is a first-order approximation for the preferred wavelength pref (a higher-order approximation is given in Appendix C). Equation (20) therefore yields an approximate relation between the preferred wavelength and three internal length scales of the system: the channel width B, the friction length scale  f , and the tidal excursion length  t . In its derivation, the approximate balance between erosion and deposition, see equation (14), is used and and u c are assumed to be zero, since including horizontal eddy viscosity in the momentum balance or critical erosion velocity u c was shown to be unimportant for the dependence on channel width. The first aspect (A1) is revealed by equation (20); in wider channels with the same mode number, is maximal for larger wavelengths̃. The second aspect (A2) also becomes apparent from equation (20). It shows that the preferred wavelength of patterns with higher mode number is smaller and hence pref decreases when the preferred mode number n pref jumps to a higher mode number.
The third aspect (A3) boils down to the statement that for a fixed channel width B, the preferred mode number n pref is lower when ≠ 0 than when = 0. This follows from a competition between the two mechanisms related to the first and second terms in equation (18). The first term increases while the second term decreases with increasing mode number n. That is, the advective term "promotes" strong lateral gradients in the bed, while the bed slope term does the opposite. As a result, if Λ eff increases, n pref decreases. In fact, as shown in Schramkowski et al. (2002), if Λ eff = 0, the preferred mode number n pref tends to infinity, but when Λ eff ≠ 0, the preferred mode number is finite. The essential difference between the situation where ≠ 0 compared to the case where = 0 is that Λ eff is larger in the former case and thereby a lower preferred mode n pref is selected. Hence, the range of channel widths for which n pref is the same is wider when horizontal eddy diffusivity in the sediment mass balance is taken into account. Van der Wegen and Roelvink (2008) found a wavelength of bar patterns that, besides its dependence on current velocity and depth, scales with the square root of the channel width. In their model, the sediment transport is governed by a total load formulation (after Engelund & Hansen, 1967), where the total transport is adjusted for bed slope effects. They state that the bed slope effect is overestimated in order to obtain realistic bottom patterns. Hence, their findings correspond to the ones above, where the dependence of the wavelength to channel width is only found if the bed slope effects are relatively strong. Furthermore,equation (20) suggests that the preferred wavelength scales with the square root of channel width. Hibma et al. (2004) used a concentration equation with the bed slope effect induced by the horizontal eddy diffusivity present in their formulation. They also found tidal bar wavelength to increase with channel width. Figure 10. The preferred wavelength pref , mode number n pref , and e-folding growth time 1∕ pref versus different parameter values for modelled tidal bars in four different tidal channels. The colors indicate the factor with which the default parameter value is multiplied, with the red squares representing the default setting for that channel. The red lines denote the measured wavelength M pref and mode number n M pref as in Table 2 and the red area the measurement uncertainty (as in Figure 3). Parameter affects the erosion, U the amplitude of the basic state tidal velocity, H the undisturbed depth, B the channel width, ⋆ the bed load bed slope term, w s the settling velocity, c c h the horizontal eddy diffusivity, c u h the horizontal eddy viscosity, and u c the critical erosion. The parameters values are as in Tables 1 and 2 Leuven et al. (2016) analyzed the fluvial bar model of Crosato and Mosselman (2009), based on Struiksma et al. (1985). They showed that, like in the tidal bar case, the number of fluvial bars and channels in the lateral direction of the channel increases with increasing river width. Furthermore, the model of Crosato and Mosselman (2009) reveals, for a constant mode number, the increase of bar wavelength with increasing river width and also the drop in bar wavelength when the mode number increases. Moreover, when their bed slope parameter is increased, the mode number decreases. This suggests that in rivers where suspended sediment transport is significant, the horizontal eddy diffusivity will increase the sensitivity of fluvial bar wavelength to river width by the same mechanism as described above. Note there are differences: fluvial bars migrate, while tidal bars can be steady (as is the case in this study with a symmetric equilibrium current).

Sensitivity Bar Wavelength to Other Parameters
In concurrence with findings of Hibma et al. (2004), van der Wegen and Roelvink (2008), and Leuven et al. (2016), it follows from equation (20) that the preferred wavelength is, besides the channel width, sensitive to changes in current velocity and channel depth. It is therefore relevant to assess the sensitivity of the modelled pref to these other parameters. Figure 10 shows the preferred wavelength pref , the preferred mode number n pref , and the e-folding time scale 1∕ pref for the four tidal channels that were introduced in section 3 (see also Figure 1). The red squares denote the results using the default parameter values for that channel (corresponding to Figure 3). In Figure 10, each default value is multiplied by factors between 0.5 and 1.5 to illustrate the sensitivity of pref , n pref , and 1∕ pref to these parameters. The figure reveals that, in particular for the Western Scheldt and Ord River estuary, pref is sensitive to both current velocity U, depth H, and channel width B. In the Exe estuary and Netarts bay, the current has less influence on the preferred wavelength than in the other two channels. The jumps in pref in the top panels correspond to jumps in mode number n pref . The growth rate is most sensitive to current velocity and settling velocity. To further analyze the sensitivity of the preferred wavelength to the current velocity and the depth, differentiate the wavelength̃with respect to the length scales B,  t , and  f to obtain, respectively, The first expression approximates the slopes of the lines in the upper panels of Figure 4 for a fixed n pref . It shows that when one of the length scales B,  t , and  f is very small (large), the sensitivity of̃to that length scale is large (small) compared to the others. The bottom lines in the top row of Figure 4 also show this behavior; the slope of the lines decreases with increasing width B. This was also observed by Leuven et al. (2016), who found that for estuaries where the current velocity is large, the bar length hardly depends on the current velocity. Also, it agrees with the results of Dalrymple and Rhodes (1995), who found that the wavelength correlates strongly with channel width when considering channels of only a few hundred meters wide.

Formulation of Horizontal Eddy Diffusivity Coefficient
The horizontal eddy diffusivity is chosen to scale with the product of a velocity scale (equilibrium current velocity U) and a length scale (channel width B) such that its order of magnitude is 10 m 2 s −1 (Deltares, 2019). One observation from Figure 4 was that the relative change of the tidal bar wavelength over the considered range of channel widths is larger when horizontal eddy diffusivity was taken into account. Figure 11 reveals that is also the case when the parameterization of is independent of channel width B.

Conclusions
The model of Schramkowski et al. (2002) for tidal bars in channels is extended to include horizontal eddy diffusivity, horizontal eddy viscosity, and a critical velocity for erosion. It is shown that this model is able to mimic tidal bar patterns with length scales and growth times that are of the same order of magnitude as those of observed bars, including the sensitivity of tidal bar wavelength to channel width. In particular, the observed sensitivity of tidal bar wavelength to channel width was reproduced by taking the eddy diffusivity into account. The inclusion of horizontal eddy viscosity and critical erosion velocity did not significantly alter the relation between the tidal bar wavelength and the channel width.
The reason the extended model shows a clear dependence of tidal bar wavelength on channel width, while in the original model of Schramkowski et al. (2002) this dependence was weak, is the following. First of all, in wider channels, the pattern with the largest growth rate (the preferred tidal bar pattern) has more bars 10.1029/2019JF005032 in the lateral direction. That is, it has a higher preferred mode number. Second, in the range of channel widths with constant preferred mode number, the preferred wavelength increases with increasing channel width. However, at the width where the preferred mode number changes (it increases by one), the preferred wavelength decreases. In the model of Schramkowski et al. (2002), the range of channel widths with constant preferred mode number is relatively short. This means that, with increasing channel width, before the preferred wavelength can increase significantly, the preferred mode number increases and with that, the wavelength decreases again. As a result, the tidal bar wavelength remains of the same order of magnitude for different values of channel width. In the present study, the effects of horizontal eddy diffusivity are added to the model. As a result, the effective bed slope effect is stronger. Since bed slope effects decrease the growth rate of bottom patterns with large gradients, the bottom pattern with the largest growth rate has a smaller mode number. This implies that the preferred mode number remains the same for a wider range of channel widths than in the model of Schramkowski et al. (2002). The width dependence now follows from the fact that, in this range, the tidal bar wavelength increases with increasing channel width and that this increase is stronger for smaller mode numbers.