Nonlinear Processes in Geophysics Patch behaviour and predictability properties of modelled finite-amplitude sand ridges on the inner shelf

The long-term evolution of shoreface-connected sand ridges is investigated with a nonlinear spectral model which governs the dynamics of waves, currents, sediment transport and the bed level on the inner shelf. Wave variables are calculated with a shoaling-refraction model instead of using a parameterisation. The spectral model describes the time evolution of amplitudes of known eigenmodes of the linearised system. Bottom pattern formation occurs if the transverse bottom slope of the inner shelf, β, exceeds a critical valueβc. For fixed model parameters the sensitivity of the properties of modelled sand ridges to changes in the number(N−1) of resolved subharmonics (of the initially fastest growing mode) is investigated. For any N the model shows the growth and subsequent saturation of the height of the sand ridges. The saturation time scale is several thousands of years, which suggests that observed sand ridges have not reached their saturated stage yet. The migration speed of the ridges and the average longshore spacing between successive crests in the saturated state differ from those in the initial state. Analysis of the potential energy balance of the ridges reveals that bed slope-induced sediment transport is crucial for the saturation process. In the transient stage the shoreface-connected ridges occur in patches. The overall characteristics of the bedforms (saturation time, final maximum height, average longshore spacing, migration speed) hardly vary withN . However, individual time series of modal amplitudes and bottom patterns strongly depend onN , thereby implying that the detailed evolution of sand ridges can only be predicted over a limited time interval. Additional experiments show that the critical bed slope βc increases with larger offshore angles of wave incidence, Correspondence to: H. E. de Swart (h.e.deswart@uu.nl) larger offshore wave heights and longer wave periods, and that the corresponding maximum height of the ridges decreases whilst the saturation time increases.


Introduction
Field data collected at various storm-dominated inner shelves of coastal seas (depths of 5−20 m) reveal the presence of patches of large-scale shoreface-connected sand ridges (Duane et al., 1972;Swift and Field, 1981;Dyer and Huntley, 1999;Harrison et al., 2003, and references herein), hereafter abbreviated as sfcr.Typically, a patch consists of 4−8 sfcr, where the latter have heights of several meters and are spaced several kilometers apart.The ridges make an angle of 20-50 • with the storm-driven current along the coast, such that their seaward ends are shifted upstream with respect to their attachments to the shoreface.Furthermore, the ridges have asymmetrical profiles, with gentle slopes on the landward (upstream) sides and steep slopes on the seaward (downstream) sides.The estimated evolution time of sfcr is several thousands of years and they migrate several meters per year in the downstream direction.As these ridges seem to affect the stability of the beach (Van de Meene and Van Rijn, 2000), gaining more understanding about their behaviour is relevant for coastal zone management purposes.
Early studies focused on the initial growth of sfcr.Trowbridge (1995) demonstrated that these ridges can grow due to feedbacks between the waves (which stir sediment from the bottom), the storm-driven current (which transports the sediment) and the sandy bottom.Calvete et al. (2001a) extended this model by including Coriolis terms, bottom shearstresses, sediment transport due to bottom slopes, suspended load sediment transport and spatial variations in stirring of Published by Copernicus Publications on behalf of the European Geosciences Union and the American Geophysical Union.sediment by the waves.Improvements with respect to the Trowbridge model were i) the occurrence of a preferred mode, ii) the results are hardly sensitive to the cross-shore profile of the storm-driven current and iii) the predicted spacing between successive ridges, as well as their growth rate and migration speed are in fair agreement with field data if all these effects are accounted for.Calvete and De Swart (2003) investigated the long-term dynamics of sfcr by expanding the flow and bottom perturbations in a truncated series of eigenmodes.The result is a set of equations describing the time evolution of the amplitudes.They showed that after the initial growth stage ridge profiles become asymmetric and reach a finite height.They also found that adding subharmonic modes (eigenmodes with wavelengths being larger than that of the most preferred mode) causes the spacing between successive ridges in the saturated state to be larger than that during the initial stage.
The models cited above have two important drawbacks.First, stirring of sediment by waves is described by a severe parameterisation.Second, they do not provide an explanation for the patchiness of observed sfcr.Therefore, in this paper a new nonlinear model will be presented, in which processes like shoaling and refraction are explicitly accounted for.This was also done in recent studies by Lane and Restrepo (2007) and Vis-Star et al. (2007), but they only examined the initial formation of sfcr.Here, an analysis will be presented of the long-term behaviour of sfcr in dependence of offshore wave characteristics and the transverse bottom slope of the inner shelf.Specific emphasis is put on the role of subharmonics in generating patches of sfcr.This is because adding subharmonics implies a higher resolution in the spectral domain, thereby potentially allowing for the occurrence of group (or modulation) behaviour.The physical processes controlling the migration and saturation behaviour will be identified by analysing the energy balance of the bedforms.Our method extends the one used by Garnier et al. (2006) in the sense that also the instantaneous global longshore migration speed is investigated.
The model is presented in Sect.2, after which the method of analysis is discussed in Sect.3. In Sect. 4 the results are presented, followed by a discussion in Sect. 5 and the conclusions in Sect.6.

Model
We adopt the model formulation of Vis-Star et al. (2007).The coupled dynamics of waves, storm-driven currents and the sandy bed is considered on an idealised inner shelf (see Fig. 1).The latter is bounded on the landward side x=0 by the transition to the shoreface (depth H 0 ) and at the seaward side (x=L s ) by the outer shelf (depth H s ).The mean transverse bottom slope is β= (H s −H 0 ) /L s .The water motion is forced by imposed wave conditions at the outer shelf (offshore angle of wave incidence s , offshore root-meansquare wave height H rms,s and wave period T ) and a wind stress that forces a net current.Sand is only transported during severe weather conditions (which occur ±5% of the time) and is the result of the joint action of waves (stirring sand from the bottom) and a longshore storm-driven flow (causing net sand transport).Hence, the model is representative for storms.After application of the rigid-lid approximation (sea surface elevation mean depth) and the quasi-steady approach (hydrodynamics adjusts instantaneously to a new bed level), the wave equations are Here, ω, g, κ and D are the wave frequency, acceleration due to gravity, wavenumber (length of wave vector κ) and water depth, respectively.The wave vector has components κ x =−κ cos θ and κ y =κ sin θ in the x-and y-direction, respectively, with θ the angle of wave incidence with respect to the shore normal.The horizontal nabla operator ∇ has components ∂/∂x and ∂/∂y in the x-and y-direction, respectively.The wave energy density E is governed by the energy balance Eq. ( 3), where c g is the group velocity vector of the waves (magnitude c g ) with components c gx =−c g cos θ and c gy =c g sin θ.The source of energy is F and D denotes the energy dissipation by bottom friction.In the definition for the wave energy density ρ is the water density and H rms is the root-mean-square wave height.Finally, the root-meansquare amplitude of the near-bed wave orbital motion u w (hereafter called wave orbital velocity) depends on the other variables and is input in the modules for the current and the sediment transport.
The currents are described by the quasi-steady depthaveraged shallow water equations Here, v denotes the depth-and wave-averaged flow velocity with components u and v in the x-and y-direction, respectively, and z=z s is the level of the mean sea surface.Furthermore, f is the Coriolis parameter, e z is the unity vector in the z-direction, τ s is the wind stress vector, τ b the mean bed shear-stress vector and e y is a unit vector in the longshore direction.Note that τ b depends linearly on the current (friction coefficient r), as waves are strong compared to currents during stormy weather.
The bed evolution equation and the formulations for the sediment transport read In Eqs. ( 8)-(11) t is time, p is the porosity of the bed, z = z b is the bed level, and q b and q s denote the wave-averaged sediment transport as bedload and suspended load, respectively.The latter two are divided into advective parts q .,aand diffusive parts q .,d(related to bed slopes).In the expressions for bedload and suspended load transport ν b , α and γ are known coefficients, C is the depth-integrated volumetric concentration of sediment and λ b and λ s are the bed slope parameters for bedload and suspended load, respectively.Note that the expressions for q b,a and q s,a in Eqs. ( 9)-(10) do not contain any net advective transport due to waves.This is a consequence of the assumption that, although wave orbital motion is strong compared to currents, the wave profiles are symmetric.Moreover, it is assumed that u w is much larger than the critical velocity u c for erosion of sediment, so the influence of u c on sediment transport is not explicitly modelled.
As boundary conditions offshore wave properties (rootmean-square wave height, period and angle of wave incidence) are imposed.Furthermore, the cross-shore flow component u vanishes at x=0 and far offshore, the bed level z b is kept fixed at these two positions and the sediment concentration vanishes far from the coast.

Basic state and linear stability analysis
From here on primary wave variables are denoted by X =(κ, θ, E) and other dependent variables by The system of equations of motion of Sect. 2 allows for a basic state that is steady and longshore uniform.It is characterised by a shelf topography z=−H (x), an incoming wave field with wavenumber K(x), angle of wave incidence (x) and wave energy E(x), a storm-driven current v=(U (x), V (x)), a free surface elevation ξ(x) and a depth-integrated volumetric concentration of sediment C(x).Hence, X =X b (x)=(K, , E) and Both X b and b are defined in Vis-Star et al. (2007); in particular U (x)=0 and V (x) results from a balance between the longshore wind stress and bottom stress.The basic state describes shoaling and refracting waves and a storm-driven flow along the coast.
The stability properties of the basic state are considered by studying the dynamics of small perturbations evolving on this basic state.Hence, and (x, y, t)=(u , v , η , c , h ) denote the perturbed variables, which are assumed to have small values with respect to their basic state values.In principle, also perturbations in the wave variables have to be considered: X =X b +X .Similar as in previous studies, we assume X =0.Physically, this means that wave-topography interactions, i.e., wave refraction and shoaling and dissipation of wave energy due to the presence of bedforms are neglected.We will return to this aspect in Sect. 5. Substitution of Eq. ( 12) into the equations of motion of Sect. 2 and linearizing the results yields the system The 5×5 matrix S contains the temporal information of the perturbations and has one non-zero element: S(5, 5)=1−p.The linear matrix operator L contains all the linear terms; its elements are given in Appendix A. This system admits solutions = e ψ(x) exp (iky where k is a longshore wavenumber, σ the complex growth rate and ψ(x) denote the cross-shore dependence of the solutions.Substitution of Eq. ( 14) into system Eq.( 13) results www.nonlin-processes-geophys.net/15/943/2008/ Nonlin.Processes Geophys., 15,[943][944][945][946][947][948][949][950][951][952][953][954][955]2008 in an eigenvalue problem, where σ are the eigenvalues and ψ(x) the eigenfunctions.This problem can be solved by standard methods.The eigenmode that has the largest growth rate σ r ≡ e(σ ) is called the initially most preferred mode; it has a wavenumber k=k p , with a corresponding wavelength p =2π/k p , and a migration speed V m =− m(σ )/k p .

Nonlinear analysis
In order to investigate the long-term evolution of sfcr the full set of nonlinear equations for the perturbations has to be considered.Here, S and L are as in the linearised system Eq.( 13), whilst N ( ) is the nonlinear vector operator, which includes all nonlinear terms in the equations of motion for the perturbations.The expression for the nonlinear vector operator N is also given in Appendix A. Following Calvete and De Swart (2003) the perturbations are written as where the unknown contributions 0 have a longshore uniform structure and the contributions ψ are expanded into a truncated series of (known) eigenmodes of the linear system: ψjn j (t) ψjn j (x) exp(ik j y) For each alongshore wavenumber k j , n j refers to the crossshore modenumber, ψjn j (t) are the unknown modal amplitudes and ψjn j (x) the known cross-shore structures of the eigenfunctions of the linear problem.Furthermore, J and N j are the largest longshore and cross-shore mode number, respectively.Expansions Eqs. ( 16) and ( 17) are substituted in the nonlinear equations of motion.After averaging over the alongshore direction, equations are obtained for the longshore uniform flow modes and bottom mode, which are subsequently subtracted from the original equations.The results are projected onto the adjoint linear eigenmodes.This procedure yields a set of nonlinear algebraic equations for the flow amplitudes ǔjn j , vjn j , ηjn j and čjn j and a set of nonlinear differential equations for the amplitudes of the bottom modes ȟjn j .A third-order time integration scheme is used to solve the resulting system of equations (see Karniadakis et al., 1991).
Here, the choice is to include the most preferred mode with wavenumber k p (largest initial growth rate) in the nonlinear expansion, as we are interested whether it is still the dominant mode in the bottom pattern on the long term.Furthermore, several superharmonics and generally some subharmonics are used, such that a total number J of different longshore wavenumbers k j is included.Thus, if only superharmonics are included, the system is solved on a domain with longshore length L y =2π/k p .By including (N−1) subharmonic modes the domain becomes longer: L y =2π N/k p .In both cases periodic boundary conditions in the longshore direction are applied.Modes that fit into the domain have wavenumbers k j =2πj/L y = j N k p (j =1, 2, 3, . ..).An individual mode is denoted by (j/N, n j ).

Energy balance of the bedforms
In order to investigate the saturation behaviour of bedforms Garnier et al. (2006) developed a method to study the global properties of the bedforms on the whole domain.It boils down to deriving a potential energy balance of bottom perturbations by multiplying the linearised version of the bed evolution Eq. ( 8) with the bed elevation h and integrating over the whole domain.It reads where Note that 1 2 |h| 2 measures the potential energy density of the bedforms and the overbar indicates the average over the domain.Equation ( 21) results from application of Green's theorem and the definitions for q b,d and q s,d .Term P , which describes the production of potential energy due to the advective bedload and suspended load transport, is positive if the advective sediment transport and bed slopes are positively correlated.Term describes the loss of potential energy due to bedslope-induced sediment transport.
An instantaneous global growth rate of the sfcr σr is defined as The definition is such that if the bed pattern is represented by a single wave h = e h(x)e iky+σ t , then σr →σ r , i.e., the initial growth rate.A new variable that is considered in this study is the instantaneous global longshore migration speed, defined as Again, for a single wave Ṽm →V m , i.e., the initial migration speed.

Parameter setting
Runs were performed with parameter values that are representative for the micro-tidal inner shelf of Long Island, located at latitude ∼40 • N. Here, a patch of 8 sfcr is observed (Duane et al., 1972).The depth varies from H 0 =14 m to H s =20 m over an inner shelf width of L s =5.5 km (β=1.1×10−3 ).However, the default experiment is performed for an inner shelf slope which is approximately 25% of its observed value.The motivation for the latter choice is that the behaviour of the system for larger values of β is not principally different, but a much higher number of eigenmodes is required to obtain accurate solutions.This aspect will be further discussed in Sect. 5.The alongshore wind stress τ sy =−0.4 N m −2 , the offshore root-mean-square wave height H rms,s =1.5 m, the wave period T =11 s and the offshore angle of wave incidence s =−20 • (waves approach from the northeast).Values of the other parameters are: A motivation for these choices can be found in Vis-Star et al. (2007).Results of the simulations are shown for continuous storm conditions.
In the default case N =10.In the time integration a time step of about 1 full-storm year was used.Both adding more modes and decreasing the time step did not change the solutions.

Linear stability analysis and eigenmodes
Figure 2 shows the initial growth rate σ r as a function of the longshore wavenumber k for the default parameter setting.The largest growth rate is attained for k=k p ∼0.6 km −1 .Thus, the initially most preferred mode has a longshore wavelength λ p ∼10 km, an e-folding time scale T g ∼1100 yr and a migration speed V m ∼−26 m yr −1 .If N=10 this mode is labeled (10/10,1).Its bottom pattern, as well as that of the (5/10,1) subharmonic mode, the (20/10,1) superharmonic mode and the (10/10,2) secondary cross-shore mode is shown in Fig. 3. Further details about the results of the linear stability analysis are shown in Supplementary Note 1 (see http://www.nonlin-processes-geophys.net/15/943/2008/ npg-15-943-2008-supplement.zip).

Patch behaviour and sensitivity to number of subharmonics
The temporal evolution of the maximum height of the bedforms is shown in Fig. 4 for different N. At t=0 all resolved bottom modes have small amplitudes (in the order of 0.1 mm) and random phases.At first, scale selection takes place: all resolved modes initially have the same amplitude and it takes time before the fastest growing mode dominates over the other modes.Next, a stage occurs in which the (10/10,1) (20/10,1) (10/10,2) (5/10,1) height of the bedforms grows exponentially, followed by saturation towards a constant finite value.The saturation time scale is defined as the time at which the maximum height is 98% of its value in the saturated state.If N=1 (no subharmonics) saturation of the bedform height takes place after a period of ∼9000 yr at a value of 0.57 m.If subharmonics are included the final height is larger (∼0.66 m) and the saturation time is slightly longer.However, both the increase in the finite bedform height and saturation time are almost independent of the choice of number N. It turns out that the saturation time also hardly depends on the initial amplitude of the resolved modes.This is because for larger initial amplitudes the stage of scale selection becomes longer, whilst the stage of exponential growth becomes almost equally shorter.
The time evolution of the amplitude of the individual modes resolved yields information about their mutual interactions.Figure 5 shows for N=10 the time evolution of the five largest bottom modes at t/T m ∼8, t/T m ∼15, t/T m ∼50 and t/T m ∼100.For the present bottom slope (β=2.7×10−4 ) the morphodynamic time scale T m ∼1000 yr.
During the first stage, up to about t/T m ∼8, the initially most preferred (10/10, 1) mode has the largest amplitude and the growth is predominantly exponential.In the course of time nonlinear processes become important and subharmonic modes, i.e., the (9/10, 1) and (8/10, 1) mode, become more mode (5/10,1) mode (10/10,1) mode (20/10,1)  dominant than the initially most preferred mode.Note that, although the bedform height is saturated for t/T m 8, the amplitudes of the individual modes are not and thus the length scale and shape of the bedforms are still changing.Slightly after t/T m =15, the (9/10, 1) mode rapidly decreases in amplitude, whereas the (8/10, 1) mode increases in amplitude and becomes dominant.Just before t/T m ∼50 the initially most preferred mode disappears from the graph and is no longer one of the five largest bottom modes.It takes a very long time before the amplitudes of the individual modes are more or less saturated.The wavelength of the saturated bedforms is dominated by the (8/10, 1) mode and is approximately 10 8 ×10=12.5 km.Furthermore, a difference in the height of individual bars and depths of individual troughs can be observed.The lengthening of patterns is consistent with Calvete and De Swart (2003) and Garnier et al. (2006), although the shift in length scale is smaller than in these studies.
Figure 6 shows the temporal evolution of the instantaneous global growth rate and migration speed for N=10.The shaded areas mark specific time periods during which a sudden change in growth rate and migration speed is visible.These changes are linked to a re-ordering of modes and will be discussed in Sect. 5.Both the instantaneous global growth rate and migration speed seem to stabilise on the long term, where the latter is smaller than during the linear evolution.The decrease in migration speed with increasing height is typical for anormal dispersive waves.Indeed, for the present  y (km) x (km) system it turns out that the group velocity of its linear eigenmodes is larger than the phase velocity.
The time evolution of the amplitude of the five fastest growing modes depends strongly on the amount of subharmonics included in the calculations.In case of N =8, N=12 or N=16, a period of t/T m ∼60 is not sufficient to reach saturation of the amplitudes of individual modes, whereas for N=10 modal amplitudes are saturated after that period.This sensitivity to N is remarkable, as Fig. 4 indicates that the saturation time scale for the bedform height shows only a very weak dependence on N. If N=8, modal amplitudes become constant after t∼100T m .At that time the (7/8, 1) mode is dominates over all others, which indicates a wavelength of the final bedforms of approximately 8 7 ×10=11.4km.If N =12 saturation of amplitudes occurs only after t/T m ∼200 and at that time the (10/12, 1) mode is dominant, indicating a wavelength of the sfcr of approximately 12.0 km.If N=16, mode saturation is found after t/T m ∼300, the (13/16, 1) mode is dominant and hence the wavelength of the final bedforms is ∼12.3 km.These results indicate that the wavelength of sfcr in the saturated state is in the range of 11.5−12.5 km, which denotes a mild lengthening compared to results obtained with the linear analysis.
A new and interesting phenomenon captured by the nonlinear model for runs including subharmonic modes is the patchiness of sfcr during their stage towards saturation.Similar patterns are also observed in the field.Stretches of the inner shelf where several ridges are present are alternated with stretches where no ridges occur.This was never obtained with the previous models.Figure 7 illustrates this patch behaviour.It shows the bottom patterns at t/T m ∼8 along a coastal stretch for runs including a different number of subharmonic modes.The patchiness is visible in all cases shown.However, the precise coastal stretches where patches of sfcr occur are definitely different and depend strongly on the choice for N.As N is a numerical parameter, the model results suggest that the detailed evolution of the sfcr is only predictable for a finite amount of time.time.Figure 8 shows the bottom patterns at different times for the case that N =10.Patch behaviour of ridges is only observed during the evolution of the system towards the saturated state.Clearly, at time t/T m =20 the entire domain is covered by sfcr.Whether sfcr occur in patches or not is strongly related to the interactions between different modes in the nonlinear analysis.Modulation behaviour will occur when (1) the spectrum of growing modes is sufficiently narrow, and (2) the dominant modes that control the evolution of the system have comparable wavelengths and amplitudes.In the experiments performed with subharmonic modes the first condition is met, whilst the second condition is only obeyed during the transient stage of the evolution.For example, Fig. 5 shows that at t/T m ∼8 the (10/10,1) mode and (9/10,1) mode have comparable amplitudes.However, at a later stage (e.g.t/T m ∼20), the amplitudes of modes with successive modenumbers are not that close.In fact, obtained patterns are a superposition of the dominant mode and its superharmonics and these imply ridge activity everywhere on the inner shelf.

Sensitivity results to wave characteristics
A series of runs were performed with the model without subharmonics (N=1) for different values of the transverse slope β (by changing the depth of the outer shelf H s ) and offshore wave parameters s , H rms,s and T .Wave parameters were varied one after another and other parameters kept their default values.All runs revealed that a critical transverse bottom slope β c has to be exceeded before sfcr start to develop, where β c depends on the model parameter values.Once β>β c the time evolution of the maximum height of the bedforms is generally characterised by growth, reduction in growth and saturation.Numerical stable solutions could be obtained up to approximately 60% of the observed value of the slope of the inner shelf in case of no subharmonics.For larger values of β solutions become singular at some point during the evolution.Instability behaviour is related to a rapid growth of the smallest length scale included in the nonlinear analysis.
Contour plots of the finite maximum bedform height and inverse of saturation time in the s −β plane were already presented and discussed in Vis-Star et al. (2008).It turns out that if waves approach the coast more obliquely, then β c increases, the final height of the ridges decreases and the saturation time increases.
Figure 9 shows the finite maximum bedform height and inverse of saturation time in the (H rms,s −β) parameter space.An increase in β c from 1.0×10 −4 to 4.3×10 −4 is obtained for an increase in the offshore root-mean-square wave height from 1.2 m to 2.0 m.Once β>β c , sfcr grow faster and become higher for smaller offshore wave heights.The dependence of the final height and the saturation time on the offshore wave height is quite strong.
Contour plots of the finite maximum bedform height and inverse of saturation time in the T −β plane (not shown) reveal that the critical transverse bed slope β c increases from 1.0×10 −4 to 2.7×10 −4 for an increase in wave period from 8 s to 14 s.For larger bed slopes sfcr evolve and it appears that the shorter the wave period, the higher bedforms become and the faster they evolve.Sensitivity of results to changes in the wave period are rather small.

Analysis of saturation behaviour
All runs discussed in the previous section show that, starting from an initial state without bedforms, sfcr develop as free instabilities of the coupled water-bottom system.Their initial growth is exponentially and can be described and understood using linear stability analysis (cf.Calvete et al., 2001a).On the long term saturation behaviour occurs and the ridges attain a finite, almost constant height.To further understand the behaviour of the bedforms their potential energy balance is analysed.Results are shown in Fig. 10.In panel b of this figure time series are shown of the gain and loss of potential energy due to bedload and suspended load transport, respectively.The time periods during which jumps occur are marked.These jumps represent significant changes in competitive behaviour of individual bottom modes included in the analysis, as is clear from Fig. 5.In Fig. 10c the jumps are also visible in the curves around |h|=0.15 m.Note that the time that the system needs to pass through the part of the curve for relative small bar amplitude is short compared to the time needed to approach the end of the curve.The alternative open and closed circles below the time axis in Fig. 10b represent five specific time values which are also plotted in Fig. 10c, d.Around saturation time the increase in |h| is rather gradual.From Fig. 10d it seems that a balance between production and dissipation is not reached.This, however, is due to numerical accuracy, as the absolute error in the computation of both the production and dissipation term is about 3×10 −7 m 2 yr −1 .
Performing the potential energy analysis for runs including different amount of subharmonic modes (results not shown) all show the same behaviour.The gross characteristics of the sfcr, e.g.instantaneous global growth rate, in-stantaneous global migration speed, alongshore spacing and finite amplitude, are a robust model result, independent of the amount of subharmonic modes included.However, changes in the competition between individual modes are quite unpredictable and also strongly dependent on the amount of subharmonic modes.As the bottom evolution is strongly dependent on parameter N the evolution of the sfcr can only be predicted for a finite amount of time.In this context it is relevant to remark that Huntley et al. (2008) demonstrated that the predictability of bottom bedforms is also strongly influenced by the presence of defects in the initial pattern.
The lengthening of patterns in the course of time is consistent with that found in other nonlinear morphodynamic models, e.g., that of Coco et al. (2004) for beach cusps, that of Murray and Thieler (2004) for sorted bedforms, thay of Ashton and Murray (2006) for shoreline sand waves and the model of Garnier et al. (2006) for nearshore sand bars.The shift in length scale is similar as in Coco et al. ( 2004), but in the other studies the lengthening is larger.

Model limitations
In reality storms only occur during a certain time fraction (about 5%) and numbers would be a factor 20 smaller.On the other hand, taking a realistic value of β would cause these numbers to be larger by approximately the same factor.
In case of using a measured value of the transverse bed slope solutions become unbounded before the saturated state is reached.The results presented in this paper were obtained with a version of the model in which interactions between waves and growing bedforms were ignored.Nonlinear experiments including wave-topography interactions revealed spurious modes in the linear analysis which have to be filtered.When incorporated in the nonlinear model sfcr continue to grow on the long term and numerical instabilities arise before the saturated state is reached.An explanation for the latter is that stirring of sediment by waves increases with increasing ridge height.This tendency will be counteracted by e.g.wave breaking, a process which is not yet implemented in the model.
The spectral method itself is also subject to limitations.The nonlinear analysis uses the solutions of the linear system as expansion modes.However, the latter are dependent on the specific parameter values used.Simulations in which strong, mild and weak storms alternate, as will be the case in reality, are feasible, but require different eigenmodes for the different realizations.This implies that several transformations between the spectral domain and the physical space have to be performed, which will cause a considerable loss of accuracy.
Note that the present model assumes a constant mean sea level, whereas in several studies (cf.Swift and Field, 1981;McBride and Moslow, 1991) it is argued that sea level changes during the Holocene play an important role in the evolution of sfcr.To assess the sensitivity of model results to different sea level stands runs were performed using different cross-shore bottom profiles.The latter were constructed by generalizing the "Bruun rule" (Bruun, 1954(Bruun, , 1962)), following concepts discussed in Masetti et al. (2008).Starting from a given profile, a sea level rise s will cause the new depth of the outer shelf to become H s + s (Fig. 11).
Second, the depth at the transition from shoreface to inner shelf remains unchanged (value H 0 ), but its location shifts landward (distance L s ).Third, the depth profile between coastline and the inner shelf maintains its equilibrium shape and it is shifted upward and landward.The distance L s is calculated by imposing that the total volume of sand per unit width between coast and outer shelf remains unchanged.This results in a new transverse bottom slope β of the inner shelf (larger for positive s ).Experiments for different values of s were performed, starting from the default situation with subharmonics (N=10).It was found that the overall behaviour of the model does not change if s is varied.The main quantitative changes are that the height of the ridges increases (from 0.65 m to 2.2 m if s =+1 m) and the saturation time scale decreases.This is mainly caused by the increase of β.
The reconstruction method we applied is highly idealised, but nevertheless yields a first idea of how sfcr might respond to changes in mean sea level.The numerical model of Masetti et al. (2008) is more advanced: it calculates the explicit time dependence of bottom profiles under the influence of sea level rise and also accounts for the underlying stratigraphy.Furthermore, the "Bruun rule" often yields results that do not comply with observations, because it ignores threedimensional and site-specific aspects (Cooper and Pilkey, 2004).Ultimately, profile reconstructions should be tested against data.Such information however is not available for Long Island shelf.and field data is therefore limited.Generally, observations show patches of 4−8 sfcr (see literature cited in the introduction).In the present study we were able to reproduce the patchiness of sfcr, which was never done before.The number of sfcr per patch are consistent with the field data.

Comparison with field observations
The spatial patterns of the finite-amplitude sfcr exhibit typical asymmetric profiles with steep slopes at the downstream sides.Extrapolating default model results to the realistic inner shelf bottom slope yields bedforms with a finite height of ∼4 m and a saturation time scale of ∼700 yr.These values have the correct order of magnitude if compared with field data (e.g.Duane et al., 1972).

Conclusions
The long-term evolution of shoreface-connected sand ridges (sfcr) has been investigated and the dependence of model results with respect to offshore wave properties and the inner shelf bed slope has been explored.For any forcing conditions (wind, waves) a critical transverse bed slope has to be exceeded before sfcr start to grow.Once this critical bed slope is exceeded, sfcr initially form as free morphodynamic instabilities and, after a stage of scale selection, their height increases exponentially in time.On the long term the sfcr attain a finite height which becomes constant and their spatial pattern becomes asymmetrical (mild stoss sides, steep lee sides).An analysis of the potential energy balance of the sfcr has been performed which shows that bed slope-induced sediment transport is crucial for the saturation process.
If subharmonic modes are included, sfcr also show saturation and the finite height increases with about 15% compared to simulations without subharmonics.The initially most preferred mode is no longer dominant in the saturated state.Due to nonlinear interactions subharmonic modes dominate and cause a 20% increase in the distance between successive sfcr in time.Furthermore, the saturation time of the amplitudes of individual modes is much longer than the time scale on which the ridge height saturates.Considering the age of the sfcr, time scales of order 10 T m are realistic, which suggests that observed sfcr have not reached their final stage of development yet.In the transient stage the sfcr occur in patches and the number of ridges per patch is in the observed range of 4−8.Results indicate that the detailed evolution of the sfcr is only predictable over a finite time interval, whereas the overall characteristics (e.g.instantaneous global growth rate, instantaneous global migration speed, finite bedform height and alongshore spacing) are also predictable for long time scales.Model results are consistent with field observations.

Appendix A Expressions for operators S, L, N
The symbolic form of the system of nonlinear partial differential equations describing the evolution of the flow, mass, sediment concentration and bottom is given in Eq. ( 15), whilst Eq. ( 13) contains its linearised version.
Here, the elements of 5×5 matrix S (containing all the temporal information of the perturbations) are all zero, except for S(5, 5)=1−p.The linear operator L involves spatial derivatives and has the matrix form The elements L 51 and L 55 of operator L are Finally, the vector N governs all nonlinear contributions.Its elements are with q =q b +q s and L(∇•q ) the linearised part of ∇•q .

Introduction
In the main paper results are shown of the linear stability analysis for the default setting of the parameters only.In this note the sensitivity of the initial growth rate, alongshore spacing and migration speed to changes in the transverse bottom slope β, the offshore wave height H rms,s and the wave period are presented.The sensitivity of these variables to variations in the offshore angle of incidence can be found in Vis-Star et al. (2008).Finally, an explanation, based on earlier literature, will be given of the physical mechanism that causes the initial growth of sand ridges.

Dependence linear results on inner shelf slope and offshore wave height
In Fig. 1 cross-shore profiles of U w are shown for different offshore root-meansquare wave heights.The larger the offshore root-mean-square wave height the larger the offshore wave orbital velocity.The onshore increase in the wave orbital velocity is stronger for higher waves due to a more intense shoaling.
In Fig. 2 contour plots of longshore spacing, growth rate and migration speed of the most preferred mode in the H rms,s −β plane are shown.Note that for offshore wave heights smaller than 1.2 m the model assumption that the wave orbital velocity is larger than the current amplitude is no longer satisfied.The longshore spacing of the most preferred mode is smaller for steeper transverse bottom slopes.The longshore spacing increases with an increase in the offshore root-mean-square wave height.The critical transverse bottom slope increases from 1.0 × 10 −4 to 4.3 × 10 −4 when the offshore root-mean-square wave height is increased from 1.2 m to 2.0 m, respectively.As long as the critical inner shelf slope is exceeded, the growth rate of the bedforms increases and the migration speed decreases with increasing transverse bottom slope of the inner shelf.The dependence of these variables on the offshore wave height is more complicated.In general, the growth rate decreases with an increase in offshore wave height.However, for offshore wave heights between 1.2 m and 1.6 m and large bed slopes an opposite trend is visible.Furthermore, it appears that for an offshore wave height of about 1.8 m the migration speed exhibits a maximum: the migration is slower for both smaller and larger offshore waves.The orientation of the bottom pattern with respect to the shoreline does not depend on the offshore wave height.

Sensitivity to wave period and inner shelf slope
As is shown in Fig. 3, the wave orbital velocity is larger for low-frequency waves.The onshore increase in the wave orbital velocity is similar for low-and highfrequency waves.The dependence of the longshore spacing, maximum growth rate and migration speed on the transverse bottom slope β and the wave period T are shown in Fig. 4. For a wave period smaller than 8 s the model assumption that a weak current limit is considered is no longer satisfied.The longshore spacing varies between about 7 km and 10 km and is larger for smaller transverse bottom slopes and low-frequency waves.In general, the initially most preferred mode evolves slower and migrates faster when the wave period is increased.Again, a critical bed slope has to be exceeded before growing bedforms are obtained.The bedforms cover the entire width of the inner shelf under all conditions and are up-current oriented.The angle between crest axis and coastline is ϕ ∼ 30 • for both low-and high-frequency waves.

Physics of the instability mechanism
The equation governing the evolution of the bed perturbations is derived from the linearised version of the bed evolution equation, the linearised continuity equation and the linearised formulations for q b and q s , respectively.The detailed expressions can be found in Eq. ( 13) and in Appendix A of the main paper.If settling lag effects are neglected, the result is where In this equation, term T0 represents the growth or decay of bedforms.Bedforms grow (decay) if ∂h /∂t > 0 (∂h /∂t < 0) above the crests.Term T1 describes the alongshore migration of the bed perturbations due to bedload transport of sediment, whereas term T2 is a consequence of the downslope sediment transport and causes diffusion of bedforms.
The different sources of instability are given by the two terms on the right-hand side of the bed evolution equation.First consider term T3 in the specific case of (3/2)ν b U 2 w → K stir = constant, which is the case discussed by Trowbridge (1995).Clearly, this term = 0 if the transverse bottom slope dH/dx = β = 0.For β > 0 sfcr grow if u and h are positively correlated, i.e., the current should exhibit an offshore deflection over the crests.As u decreases seaward, due to the increase in water depth, the current is convergent.As here sediment transport is linearly related to the current, deposition of sediment occurs on top of the crests and sfcr will grow.An obliquely oriented ridge can induce a deflection of the longshore current as a consequence of water mass continuity, which causes an increase in the cross-bank component of the flow over the crest.Only for up-current oriented sfcr the deflection in the current is directed offshore above its crest, thus only up-current rotated sfcr will grow.The offshore current deflection over sfcr is indeed reproduced by the model (results not shown here).
A more general case is considered by Calvete et al. (2001), who consider term T3 and T4 in equation ( 1), where U w now depends explicitly on H and C ∼ HU 3 w (see Eq. ( 11) of the main paper).They argue that stirring of sediment by the waves increases towards the coast, which is an additional source of sediment deposition over the crests and thus growth is enhanced.More generally, sfcr grow if the cross-shore gradient of the depth-averaged volumetric sediment concentration in the basic state is negative.The latter quantity is defined as d(C/H) dx , thus U 3 w should decrease with increasing distance from the coast.It appears that term T4 is dominant over term T3.A sketch of the Trowbridge and Calvete mechanisms is given in Fig. 5.The offshore deflection of the longshore current over an up-current oriented ridge is illustrated in the top view.The cross-section through the ridge shows (1) a convergence of the flow as it enters deeper water and (2) a decrease in the wave orbital velocity into deeper water.As a consequence sediment transport is convergent and sediment is deposited on the crests.The convergence is most effective on the downstream side of the ridge again due to flow convergence, which explains the downstream migration of the bedforms.Sensitivity of the model results (longshore spacing, growth rate and migration speed) to offshore wave characteristics can be understood by considering the magnitude of the different terms in equation (1).Whether bedforms have the tendency to grow or decay depends on the competition between advective terms (T3 and T4, of which T4 is dominant) and the diffusion term (T2) which render unstable and stable bottom perturbations, respectively.If advective terms are larger than the diffusion term then bedforms can grow, otherwise they decay.Variation in parameter values leads to changes in the growth rate if the absolute change in the magnitude of the advective terms is different from that of the diffusive term.Furthermore, changes in the longshore spacing of bedforms occur if there is a change in the magnitude of the ratio of advective and diffusive terms.Term T2 is proportional to U 5 w , whereas term T4 is proportional to d Let us now consider the case of increasing the bed slope β.Here, an increase in β corresponds with an increase in H s and thus a larger water depth at the outer shelf.
Across the major part of the inner shelf, an increase in β results in a smaller wave Table 1: Sensitivity of model results to the transverse bottom slope and offshore wave characteristics.For explanation see the text.
orbital velocity and thus T2 ∼ U 5 w decreases.On the other hand, considering Fig. 6, the cross-shore gradient in the depth-averaged sediment concentration increases for a larger slope.Together with an increase in the magnitude of V , due to reduced friction, term T4 clearly increases.Hence, the growth rate of bedforms becomes larger for larger β.The relative increase of T4/T2 with increasing β results in a preferred wavelength which is smaller on steeper inner shelf slopes.As the migration is determined by T1∼ U w /H, it explains the decrease in migration speed with increasing β.In a similar way the dependence of growth rate, migration speed and longshore spacing of bedforms on the magnitude of the offshore wave height can be understood.The wave orbital velocity is proportional to H rms,s and therefore T2 is larger for higher offshore waves.The dominant advection term T4 also increases (see Fig. 6), but not as fast as T2.The absolute increase in T2 is so large compared to the absolute increase in T4 that bedforms grow slower in case of larger offshore wave heights.As the ratio T4/T2 also increases, bedforms have longer spacings for higher offshore waves.In general, bedform migration increases with the presence of higher offshore waves, which corresponds with the increase in the magnitude of U w .Note that under conditions that term T3 and T4 become large (e.g. for high offshore waves), their imaginary parts will no longer be small compared to term T1 and will therefore also contribute to the migration speed.This might explain the decrease in migration speed when waves become very high (see Fig. 2).The sensitivity of model results for a variation in the wave period and the angle of wave incidence can be explained using similar arguments.A summary is given in Table 1.

Fig. 1 .
Fig. 1.Schematic representation of a typical time-averaged bottom topography of the continental shelf, representing the inner and outer shelf.Forcing of the water motion is due to (obliquely) incident waves and a storm-driven current.Symbols are explained in the text.

Fig. 2 .
Fig. 2. Linear growth rates for cross-shore mode n j as a function of the longshore wavenumber.The most preferred mode is indicated by (10/10, 1).The solid dots indicate modes that are included in the nonlinear analysis for N =1.The open dots represent additional modes included if N=2.

Fig. 4 .
Fig. 4. Maximum height of the bedforms versus t/T m for model runs in which (N−1) subharmonic modes are included.For the present bottom slope T m ∼1000 yr, whereas for the realistic slope T m ∼100 yr.

Fig. 5 .Fig. 6 .=Fig. 7 .
Fig. 5. Time evolution of the amplitude of the five bottom modes which have the largest amplitude at (a) t/T m ∼8, (b) t/T m ∼15, (c) t/T m ∼50 and (d) t/T m ∼100.For the present bottom slope T m ∼1000 yr, whereas for the realistic slope T m ∼100 yr.Modes are indicated by a specific colour and denoted as (j/N ,n j ), where j is the longshore modenumber, (N −1) the amount of subharmonic modes included and n j the cross-shore modenumber, respectively.Line style(solid, dotted, dashed, dot-dashed, dot-dot-dashed)  indicates the order of dominance of modes at the final time step.Here, N=10.

Fig. 8 .
Fig.8.Bottom pattern (light: crests, dark: troughs) at different times t/T m ; N =10.For the present bottom slope T m ∼1000 yr, whereas for the realistic slope T m ∼100 yr.The shoreface is at the top and left is downstream.Here, the longshore extent of the domain is ∼N p .

Fig. 9 .
Fig. 9. Contour plot of equal (a) finite maximum height (m) of the bedforms and (b) inverse of time needed for saturation (10 −3 yr −1 ) in the H rms,s − β plane.Obtained with nonlinear model (solid lines) or extrapolated (dotted lines).
Fig. 10.(a) Time evolution of the bedforms.(b) Time evolution of production and dissipation terms due to suspended load and bedload sediment transport and their sum.Shaded areas indicate time periods during which jumps occur in all the curves.The alternative closed en open circles below the time axis are positioned at t/T m ∼5, t/T m ∼10, t/T m ∼22, t/T m ∼43 and t/T m ∼65 (after t/T m ∼65 no change occurs).These times are also indicated in panels c and d.For the present bottom slope T m ∼1000 yr, whereas for the realistic slope T m ∼100 yr.(c) Square root of the production and dissipation terms as a function of |h|.(d) The sum (P + ) as a function of |h|.Here, N=10 and β=2.7×10 −4 .

Fig. 11 .
Fig. 11.Sketch of evolution of the shelf geometry under sea level rise s .In the case of Bruun (1954): A=B.In case of Masetti et al. (2008): A=B+C.Other variables are defined in the text.

{Figure 1 :Figure 2 :
Figure 1: Cross-shore profiles of the wave orbital velocity in the basic state for different offshore root-mean-square wave heights.

Figure 3 :Figure 4 :
Figure 3: Cross-shore profiles of the wave orbital velocity in the basic state for different wave periods.

Figure 5 :
Figure 5: Schematic view of the Trowbridge and Calvete mechanisms.(a) An up-current rotated ridge causes an offshore deflection of the current.(b) Flow converges when it enters deeper water, which causes sediment convergence over the ridge.This effect is enhanced by nonuniform wave stirring: stirring of sediment by the waves is stronger in shallow water compared to deep water and thus also causes convergence of sediment over the ridge.
. The latter estimate follows from water mass continuity.Note that V itself is inversely proportional to U w due to frictional effects.The magnitude of the bedform migration is determined by term T1, which is proportional to U 2 w V /H ∼ U w /H.

Figure 6 :
Figure 6: Cross-shore profiles of the cross-shore gradient of the depth-averaged volume concentration in the basic state for (a) different offshore angles of wave incidence (b) different offshore root-mean-square wave heights (c) different wave periods and (d) different values for the transverse bottom slope of the inner shelf.Other parameters have their default values.