The model describing the surf zone hydrodynamics is based on the depth-averaged and time-averaged momentum and continuity equations coupled to the wave energy and phase equations. The momentum balance and water mass conservation equations read

In this notation, the Einstein convention is adopted; i.e. if an index appears twice in a term we assume a summation over that index. Here, x_i indicates the horizontal spatial coordinates (x₁∕x and x₂∕y are the cross-shore and alongshore directions, respectively), t is time, the vector $\mathit{v}(x_{\mathrm{1}},x_{\mathrm{2}},t)$ is the wave-averaged and depth-averaged mass flux current, $z_{\mathrm{s}}(x_{\mathrm{1}},x_{\mathrm{2}},t)$ represents the mean sea level and $D=z_{\mathrm{s}}-z_{\mathrm{b}}$ is the total mean depth, by which $z_{\mathrm{b}}(x_{\mathrm{1}},x_{\mathrm{2}},t)$ is the mean sea bottom bed level. Furthermore, $S_{ij}^{\mathrm{W}}$ and $S_{ij}^{\mathrm{R}}$ are the radiation stresses due to waves and rollers, while $S_{ij}^{\mathrm{t}}$ represents the turbulent Reynolds stresses, τ_bi indicates the bed shear stress, g is gravity and ρ is the water density.

The wave energy balance equation reads

where $E=\frac{\mathrm{1}}{\mathrm{8}}\mathit{\rho }g{H}_{\mathrm{rms}}^{\mathrm{2}}$ is the wave energy density, with H_rms being the wave height, c_gj are the components of the group velocity and 𝒟^W represents the wave energy dissipation rate due to wave breaking. The roller energy balance equation reads

R is the energy density of the rollers, c_j are the components of the phase velocity and 𝒟^R represents the wave energy dissipation rate due to the rollers. The wavenumber $\mathit{K}(x_{\mathrm{1}},x_{\mathrm{2}},t)$ ($\mathit{K}=(K_{\mathrm{1}},K_{\mathrm{2}})$) of the waves obeys the equation

where σ and ω are the intrinsic and the absolute wave frequencies, respectively. The wave energy dissipation rate is parameterized using the formulation by Church and Thornton (1993):

in which B (B³=2.2) is a parameter describing the type of breaking, and γ_b (=0.42) is the expected saturation value of H_rms∕D. The roller energy dissipation rate is modelled following Ruessink et al. (2001):

where β_rol (≤0.1) is the angle of the wave–roller interface. Wave radiation stresses, stresses due to roller propagation and turbulent Reynolds stresses (Svendsen, 2006) are expressed as

where δ_ij is the Kronecker delta symbol and the turbulent kinetic diffusivity is

with M (M=1.0) being a parameter that characterizes the turbulence. With respect to shear stresses, we use a linear friction law τ_bi=ρμv_i ($i=\mathrm{1},\mathrm{2}$) with $\mathit{\mu }=\left(\frac{\mathrm{2}}{\mathit{\pi }}\right)c_{\mathrm{D}}{u}_{\mathrm{rms}}$. We model the drag coefficient as

where z₀ (z₀=1.0 cm) is the bed roughness and u_rms is the wave orbital velocity at the edge of the wave-induced boundary layer:

To simulate the morphological evolution, the hydrodynamic field must be coupled to a sediment transport formulation and to the conservation of sediment mass. Bed evolution is described as

where z_b represents the mean sea level, p (=0.4) is the porosity of the seabed and q_j are the components of the volumetric sediment transport, whose parameterization is given by the Soulsby–van Rijn formula (see Soulsby, 1997), expressed in the form

where A_s depends on the sediment properties and γ is a bed slope coefficient. The term A_s(u_stir)^2.4 is the depth-integrated sediment concentration (C_di). Following Ribas et al. (2012), u_stir is a stirring velocity that takes into account the depth-averaged currents, the wave orbital velocity and the roller-induced turbulence velocity. The last term in the equation takes into account the tendency of the sand to move downslope towards an equilibrium profile, where $h(x_{\mathrm{1}},x_{\mathrm{2}},t)$ stands for the perturbation of the sea bottom with respect to an alongshore uniform background (equilibrium) bathymetry.

The system of equations, when alongshore uniformity is assumed, allows for a state of morphodynamic equilibrium (steady state) for the hydrodynamic forcing conditions. By following a standard linear stability analysis (e.g. Dodd et al., 2003; Calvete et al., 2005), the system of equations is linearized with respect to alongshore periodic perturbations of the form $\mathrm{\Phi }(x_{\mathrm{1}},x_{\mathrm{2}},t)=\mathit{\text{Re}}\left\{\mathit{\varphi }\left(x_{\mathrm{1}}\right)e^{\mathit{\sigma }t+ik{x}_{\mathrm{2}}}\right\}$, where Re stands for the real part of a complex variable. For each alongshore wavenumber k, an eigenvalue problem is solved and a number of modes ϕ_i with eigenvalues σ_i are found. For a given set of forcing conditions (wave height and period; normal incidence is assumed throughout this study) and a cross-shore profile, the output of the analysis consists of the cross-shore pattern ϕ, the alongshore wavelength of the pattern $\mathit{\lambda }=\mathrm{2}\mathit{\pi }/k$, the growth rate Re(σ), and the characteristic growing time or e-folding time, $\mathit{\tau }=\mathrm{1}/\mathit{\text{Re}}\left(\mathit{\sigma }\right)$. Boundary conditions and more details about the model are given in Ribas et al. (2012) and in Calvete et al. (2005).

The model has been applied to a series of different bathymetries to study the effect of the distance, Δx, and difference in water depth, ΔD, between the two sandbars. Figure 2 shows the series of cross-shore profiles that will be considered in this study. We have tried to isolate individual effects and, for example, profiles in red will specifically address the sensitivity to the difference in water depth between the two sandbars (notice that ΔD is varied by changing the depth over the outer sandbar). Similarly, profiles in blue will directly assess the role of the distance between sandbar crests.

Figure 2Geometry of cross-shore beach profiles used in this study. Different colours are used to highlight that some of the profiles were specifically designed to analyse the effect of variations in the distance or in the relative depth between sandbar crests.

Download

We initially present an example of the model analysis for a specific bathymetry (alongshore uniformity of the initial cross-shore profile is considered) and offshore wave conditions. For this case, we use a significant wave height H_rms=1.5 m and a wave period T=10 s with normally incident waves. The first step of a linear stability analysis is evaluating the equilibrium state, which represents the morphodynamic equilibrium previously discussed, of the equations presented in the previous section considering a fixed seabed. We assume that the bathymetry of the equilibrium state results from a morphodynamic evolution that occurs over a long temporal scale compared with the growth of the emerging morphological pattern.

Figure 3 shows the bottom cross-shore profile that is characterized by the presence of two sandbars with crests at about 200 and 480 m in the cross-shore direction, with the distance between the sandbar crests Δx=280 m and a difference of about ΔD=2.5 m between the bar depths. The other panels show other characteristics of the hydrodynamic and sediment transport (for example, notice the effect of the sandbar on wave transformation). The basic state, different for different cross-shore beach profiles, is then perturbed and possible emerging modes are analysed in terms of their growth rate. Figure 4 shows the growth rates for the example being analysed. Three different modes are present, with the fastest growing one, mode 1, characterized by an alongshore spacing close to 420 m (the wavenumber is about 0.015 m⁻¹). The second and third modes are characterized by slower growth rates and an alongshore spacing close to 170 and 500 m, respectively.

Figure 3Basic state variables for a typical cross-shore beach profile (shoreline is at x=0). From top to bottom: (a) water depth, (b) wave height, (c) wave energy, (d) roller energy and (e) depth-averaged sediment concentration.

Download

Figure 4Growth rates as a function of the wavenumber for H_rms=1.5 m and T=10 s. The bathymetry considered in this case is the same as presented in Fig. 3.

Download

Figure 5Fastest growing modes for the peaks in the growth rates shown in Fig. 4. The water depth (left) and depth-averaged currents (right) of each of the three modes are shown for bottom perturbations of 0.5 m of amplitude. The maximum velocities are indicated for each mode.

Download

The water depth and circulation patterns associated with the fastest growing value of each of the three modes are shown in Fig. 5. The patterns display some evident differences with respect to which of the two sandbars develops a crescentic shape. The mode 1 represents the classic emerging crescentic sandbar, and only the inner sandbar becomes crescentic. Circulation over the inner sandbar consists of onshore flow over the shoals and offshore flow in the lower and/or channel areas consistent with the traditional mechanism of crescentic sandbar or rip channel formation (Falqués et al., 2000; Calvete et al., 2005). The pattern of mode 2 is the result of morphodynamic feedbacks mainly acting in the zone between the inner sandbar and the shoreline. Circulation and morphology also develop close to the shoreline in the form of transverse bars aligned to the lower and/or channel areas of the inner sandbar crescents. Finally, mode 3 shows an instability of the outer sandbar with small in-phase signatures on the inner sandbar. The growth rate of the different modes can be understood following Ribas et al. (2015). The pattern related to the fastest growing mode, mode 1, arises in the areas of more intense dissipation of wave energy (both in wave and roller energy, Fig. 2) and where the gradients in depth-averaged sediment concentration are large (Fig. 2). Similarly, mode 2 is associated with an instability extending close to the shoreline, where the gradient in depth-averaged concentration leads to the development of transverse sandbars associated with an offshore flow (Ribas et al., 2015). Mode 3 is characterized by less intense circulation and depth-averaged concentration that extends to the inner sandbar.

Applying linear stability analysis to the beach profiles shown in Fig. 2 results in a variety of beach responses, each identified by a specific mode. The patterns that are predicted to emerge vary largely, and we have attempted to group them using a criterion based on the difference in the maximum amplitude of the perturbations over the inner and outer sandbars. If the amplitude of the perturbation of one of the sandbars is over 80 % larger than the amplitude of the other sandbar, we consider that only the sandbar with the largest perturbation amplitude will develop into a crescentic sandbar. If the amplitude of either the inner or outer sandbars is between 40 % and 80 % larger, that sandbar will dominate the coupling. If the difference in the maximum amplitude perturbation is below 40 %, the two sandbars are considered to be fully coupled. Also, just as shown in Fig. 5 for mode 2, several emerging configurations also involve changes close to the shoreline. If the amplitude of the perturbation close to the shoreline is at least 20 % of the largest amplitude, we consider that a pattern emerges close to the shoreline. It turns out that such shoreline patterns never dominate the dynamics in the present simulations, but they always appear coupled to a perturbation in the inner or outer sandbar. In Fig. 6 we show the different patterns obtained and use the letters I and O to indicate patterns that are associated only with the inner or outer sandbar, respectively. The symbols + and – are used to indicate possible in-phase or out-of-phase coupling so that, overall, a pattern indicated with the symbols O+ refers to a configuration where the dominant effect of the instability is over the outer sandbar (letter O), while the inner sandbar shows some limited in-phase coupling (symbol +). When the coupling between sandbars is obvious, we denote the patterns with the letters IO and add the symbol + or –, depending on whether the sandbars show in-phase or out-of-phase coupling. The possible effect on shoreline and/or inner surf zone morphology that is observed in some modes is indicated by the subscript s.

Figure 6Morphological patterns obtained in this study. The letters I and O indicate a dominance of the inner and outer sandbars, respectively. Patterns indicated with the code IO represent modes where both sandbars are simultaneously unstable. The symbols + and – are used to indicate possible in-phase or out-of-phase coupling. Modes affecting shoreline and/or inner surf zone morphology have been indicated by the subscript s.

Download

Figure 7Morphological patterns for three unstable modes as a function of ΔD and Δx. Each code represents a different pattern (detailed in Fig. 6). The top (a–c), centre (d–f) and bottom (g–i) panels represent results obtained for a wave height equal to 1.0, 1.5 and 2.0 m.

Download

Figure 8Growth time (a–c) and rip channel spacing (d–f) as a function of ΔD. The left (a, d), centre (b, e) and right (c, f) panels represent results obtained for Δx=200, 270 and 320 m, respectively. Blue, red and green symbols refer to wave height equal to 1.0, 1.5 and 2.0 m.

Download

We have run simulations over the bathymetries presented in Fig. 2 using three different values of wave height (equal to 1.0, 1.5 and 2.0 m) and keeping the wave period fixed (equal to 10 s). Results are presented in Figs. 7 and 8. Figure 7 shows the emerging modes as a function of wave height, sandbar distance (Δx) and sandbar depth difference (ΔD). Figure 8 shows the corresponding growth rates and spacings. Three unstable modes are usually present, but when wave height is smallest (H_rms=1.0 m), only two modes are unstable. The first mode, the fastest growing one, displays a similar pattern for the three values of the wave height considered. When the difference between the sandbar crests, ΔD, is large, the fastest growing mode is of type I, which implies that the inner sandbar develops into a crescentic shape. Because of the large difference in water depth between sandbar crests, the outer sandbar is essentially inactive, while when ΔD is small most of the wave breaking is concentrated on the outer sandbar, which develops crescents (type O). For intermediate differences in the water depth between the sandbar crests, a transition from type I to O can be observed. In most cases, the transition occurs through the development of an I+ pattern (the amplitude of the pattern is larger at the inner sandbar, and the outer sandbar reflects limited in-phase coupling). As ΔD decreases, an instability of type O– is more likely to develop (the amplitude of the pattern is larger at the outer sandbar, and the inner sandbar reflects limited out-of-phase coupling). For H_rms=1.5 m, the transition also results in the development of fully coupled patterns, namely type IO. While the patterns show an evident dependence on ΔD, the role of Δx on the emergent patterns is limited (Fig. 7). The pattern of the second mode, characterized by lower growth rates, can be characterized as follows. Given a specific combination of ΔD and Δx, if the mode 1 instability is of type I, then for the mode 2 the instability is type O. No mode 1 configuration affects shoreline morphology, while modes associated with changes at the shoreline appear more frequently as mode 2 and 3, especially if ΔD is large.

In order to understand the differences between the IO modes and the modes I or O, additional experiments have been carried out. For example, simulations for which modes I+ or I– were originally found were repeated, but without sediment transport in the outer bar. As a result, modes with similar growth rates and spacing were found, but with no extension to the outer sandbar. The same results were found in the equivalent experiments for O+ and O–. For conditions leading to modes IO, cancelling the sediment transport over any of the two sandbars resulted in pattern emergence focused over a single bar (the one where sediment transport was not cancelled), but now with large differences in the growth rate and spacing compared to the IO mode. IO should then be considered as a mode that develops as a result of the simultaneous interaction between the two sandbars.

Given that the sensitivity of the characteristics of the growing sandbars to Δx is limited, we fixed its value (equal to 200, 270 and 320 m) and specifically looked at the growth rates and spacings (Fig. 8). Results can be interpreted by looking at the dependencies of the individual types of patterns. For example, independent of the value of wave height, patterns of type I consistently show a marked decrease in the growth time, with increasing ΔD (see top panels in Fig. 8). The decrease in growth time is accompanied by a moderate increase in the spacing of the rip channels (the spacing of mode I never exceeds 400 m). The same is observed for all distances between the sandbar crests considered in this study. Mode O shows almost the opposite because the growth time increases with ΔD, while the spacing of the rip channels diminishes (for Δx=320 m a slight increase in spacing is observed for very large ΔD). The largest rip channel spacing observed for mode O is in excess of 1200 m, which is about twice the largest spacing observed for mode I. Finally, the modes that have a shoreline signature are all characterized by large values of ΔD, large growth times and short spacing (about 200 m).

We focused on the morphodynamics of double sandbar systems and investigated under which conditions the systems are unstable to perturbations that ultimately result in the development of surf zone patterns like rip channels and/or crescentic sandbars. We used linear stability analysis to discover the morphological configurations that can arise as a result of the feedbacks between hydrodynamics, sediment transport and morphological change. We primarily focused on the sensitivity to the initial seabed cross-shore profile, by varying the distance between sandbar crests or varying the difference between the water depth over the two sandbar crests. We generally observe large variability in the response of the system to changes in bathymetric details. This is not entirely unexpected, since the amount of wave breaking induced by the 3D morphology of the outer sandbar is critical to determine if the two sandbars are coupled or if the pattern emerges in correspondence of only one of the two sandbars. This is in agreement with the findings by Castelle et al. (2010a) showing that the type of horizontal flow circulation over the outer sandbar (driven either by refraction or by wave breaking) is ultimately responsible for the possible coupling between sandbars. For this reason, while the distance between the sandbar crests is unimportant in determining which pattern emerges (Fig. 7), the difference in water depth is a critical parameter to determine the shape and characteristics of the fastest growing mode. In our model, a large difference in the water depths over the two sandbar crests (e.g. ΔD=2.5 m) implies that limited wave set-up and breaking occur over the outer sandbar, which is essentially inactive. In this case, the fastest growing mode is always related to the inner sandbar, which becomes unstable (Fig. 9) following the physical mechanisms described by Calvete et al. (2005). When the difference between the sandbar water depths is small (e.g. ΔD=0.5 m), strong wave breaking occurs over the outer sandbar and the fastest growing mode is related to the outer sandbar (Fig. 9). Coupling between the two sandbars occurs for intermediate differences in the water depth of the sandbar crests, while the presence of emerging configurations that involve the shoreline only occurs for the largest water depth difference (Fig. 9). This behaviour is also evident when looking in detail at the spacing of the emerging rip channel pattern and at the growth time of the mode (Fig. 10). As synthesized in Figs. 9 and 10, results bear little dependency on Δx and the overall behaviour of the system is governed by ΔD and H_rms. When ΔD is small, the presence and characteristics of an emerging mode depend on the value of H_rms. As shown in Fig. 10, for small ΔD the outer sandbar spacing depends on H_rms, but tends to be large, while for large ΔD the dependency of the spacing to H_rms is smaller. Inner sandbar modes dominate instead for large ΔD and small H_rms.

Figure 9Sketch of the most likely fastest growing modes as a function of the geometry of the cross-shore profile. ΔD represents the difference between the water depth over the two sandbar crests while Δx is the distance between the two sandbar crests.

Download

Figure 10Sketch that summarizes our findings in terms of (a) growth time and (b) crescentic sandbar spacings. The colours refer to the modes in Fig. 9.

Download

Regarding the morphological coupling discussed by other authors (Castelle et al., 2010a; Price et al., 2014), our results derived from linear stability analysis can distinguish between modes that develop in one of the bars only and modes that induce the emergence of a pattern over the other sandbar. At the same time, we obtain modes that develop simultaneously over the two sandbars. In the first case, we interpret that there is a primary mode affecting one of the sandbars, with the other sandbar evolution being passively slaved to its morphodynamics. In the second case, the pattern developing over the two bars is related to the same mode and, therefore, the emerging pattern shows full sandbar coupling. This full morphodynamic coupling occurs for intermediate differences of sandbar depth. For small differences of depth, pattern emergence over the outer bar dominate, whilst for larger differences of bars' depths the main pattern is located at the inner bar (although the wavelength of the crescentic bars on the inner and on the outer bar appear to be very similar). For cross-shore profiles that allow for large waves to reach the shoreline, the model predicts the formation and coupling of shoreline patterns, even though the model does not include swash dynamics, and we considered a fixed shoreline. The transition from forced to fully coupled occurs smoothly in the parameter space that has been examined. Since our analysis of the model dynamics is linear, the concept of coupling is limited to the initial morphological formation and, since linear stability analysis focuses on the fastest growing wavelength, coupling at half of the outer bar wavelength cannot occur. Also, we do not simulate the nonlinear interactions between competing wavelengths, which might lead to coupling over longer timescales (days to weeks) or the final equilibrium configuration. Both important aspects can be studied using analysis that includes nonlinear mode interactions and that are suited to studying the long-term evolution and possibly the equilibrium of these systems.

Despite our attempts to provide a detailed description of hydrodynamics and morphodynamics, the model remains simplified and does not include a number of physical processes that in the context of surf zone morphodynamics can be relevant. As for the case of many surf zone morphodynamic studies, hydrodynamic forcing is simplified, and the effects of directional and frequency spread in the wave field as well as tidal variations are neglected. One could expect that the primary effect related to these processes was a decrease in the growth time of the features without necessarily affecting the type of morphodynamic patterns that grow. We also neglected the role of wave angle (we only considered normally approaching waves) that has been shown to be relevant for the coupling of sandbar systems (Price and Ruessink, 2011; Price et al., 2013). On the other hand, we include a detailed modelling of the effect of wave-induced rollers that has been shown to be important for the development of surf zone features (Ribas et al., 2011; Calvete et al., 2012), but whose effect on double sandbar systems had not been considered before. Finally, the study does not address some of the possible effects on sediment transport associated with undertow and wave asymmetry and, particularly for varying cross-shore beach profiles, it could quantitatively affect the results. Despite these shortcomings, the model reproduces morphodynamic patterns, which are consistent with the presence of coupled sandbar patterns. Although the objective of this contribution is limited to a numerical analysis of the possible emerging patterns arising in double sandbar configurations, model results are in qualitative agreement with observations of the Truc Vert (France) double sandbar system (Castelle et al., 2015), where transverse bars are coupled to inner bars during moderate conditions, and inner–outer bar coupling is observed for more energetic conditions (we stress that parameter settings are not necessarily representative of Truc Vert). Lack of detailed and systematic measurements of bathymetric evolution of coupled sandbar systems remains the biggest obstacle to model testing in this area of research. We envisage that future development in the extraction of bathymetry from video images will be hugely beneficial to this area of research (Van Dongeren et al., 2008).

Our findings have clear implications for the understanding of observed coupled sandbar patterns. Coupled sandbar systems are usually considered as the result of one sandbar affecting another. Our results indicate that coupling can also emerge as a result of a single mode. The apparent differential growth of each sandbar might lead one to think that one sandbar is forcing the coupling over the other sandbar. Our results indicate that a coupled pattern, with perturbations over each sandbar of different amplitude, can also arise without invoking one sandbar as a forcing mechanism. In addition, our results indicate that a variety of modes can grow for similar conditions. Although we do not deal with the nonlinear behaviour of the patterns, one can envisage that growth and interaction between multiple modes can become a source of spatial variability in the observed pattern.

For details about the numerical model, feel free to contact the authors.

GC and DC designed the study from inception to dealing with the execution, the analysis of the numerical experiments and the writing of the first draft of the manuscript. DC built the initial version of the model and FR contributed (a lot) to subsequent versions. FR, HEdS and AF contributed to subsequent manuscript drafts.

The authors declare that they have no conflict of interest.

This research has been supported by the Ministerio de Economía y Competitividad, Secretaría de Estado de Investigación, Desarrollo e Innovación (grant no. CTM2015-66225-C2-1-P), the European Regional Development Fund (grant no. CTM2015-66225-C2-1-P), and the Ministry of Business, Innovation and Employment (grant no. C05X0907).

This paper was edited by Orencio Duran Vinent and reviewed by two anonymous referees.