On the stability of a class of shoreline planform models

The evolution of beaches in response to the incident wave conditions has long attracted the attention of researchersandengineers.Apopularmathematicalmodeldescribingthechangeinthepositionofasingleheight contour on the coastline assumes that the beach pro ﬁ le is stable and the plan shape evolves due to wave-driven long-shore transport. Extensions of this model include more contours and allow for beach pro ﬁ le alteration through cross-shore transport of sediment. Despite this advantage, models with multiple contours remain relativelyunderused.Inthispaperweexaminethestabilityofthisclassofmodelforthecasesofonetothreecon-tours.Unstablemodesmayexistwhenthereismorethanonecontour.Theseincludeshortwaveswhosegrowth rateisstronglydependentuponwavenumber.Forthecaseofthreecontoursanadditionallongwaveinstabilityis possible.A necessary,butnot suf ﬁ cient, condition for instability isfound. Itrequiresa reversal of transportdirec- tion amongst the contours. The existence of these instabilities provides a possible explanation for the dif ﬁ culties found in implementing computational multi-line models, particularly where structures alter the natural longshore transport rates so they satisfy, locally, the condition for instability.


Introduction
The problem of beach stability, particularly with respect to shoreline response to the construction of coastal defences, has long been a familiar question in coastal research. Early work can be traced back to the 1950s, with continuing research since that time. Interest in this problem stems from the observation of instabilities in the natural environment. It is also germane to the problem of formulating numerical models that describe the time evolution of the shoreline. Whilst it was common to commission laboratory studies to test the performance of proposed coastal structures the era of cheap computing has meant that computational modelling has become an important design tool that doesn't suffer from the scaling problems inherent in modelling sediment transport in laboratory models. However, computational modelling has brought new challenges, particularly in view of the need to have reassurance that results from a computer model are reliable, stable and robust. One source of checks on the performance of a computer model is analytical solutions. As a general rule such solutions are available only for simplified situations or a restricted range of conditions. Nevertheless, analytical solutions remain one of the best forms of basic computer model testing. Checking the performance of computer models in more complex situations requires comparison against laboratory or field experiments. Current design practice typically relies on several if not all of the above approaches. In this paper we focus on analytical methods and the characteristic behaviour of beach models when some of the simplifications adopted for analytical solutions are relaxed. Whilst this might at first appear to be an exercise of solely academic interest the results have relevance for a class of computational models known as line models. The simplest of these models is the 1-line model, which describes the movement of a single contour line in response to incoming waves and sources and/or interruptions of sediment supply along the beach. This model has a fairly restrictive assumption that the beach profile remains unchanged. It was only natural that researchers wished to remove this restriction, so models describing more contour lines were developed so as to account for changes in beach profile shape. Nevertheless, these 'better' models which might have been expected to provide more realistic results, have not entered mainstream use. One reason seems to have been difficulties encountered in achieving solutions that converge. One potential reason for this is numerical instability, in which small fluctuations in the solution can grow rapidly, overwhelming the features of engineering interest, and rendering the computed solutions meaningless. In this paper we analyse the stability properties of the 1, 2 and 3 line models to investigate whether such instabilities might be an inherent property of the systems of equations, thereby rendering the computational approach rather problematic in certain situations. In Section 2 the 1-line model, which is widely known in the coastal engineering literature, is discussed briefly. In Sections 3 and 3.1 the stability of the governing equations to small perturbations in the beach position is assessed. Conclusions from the assessment are presented in Section 3.2.

The 1-line model and beyond
Pelnard-Considère (1956) undertook laboratory experiments and proposed a simplified analytical model to describe the long-shore sediment transport driven by waves. Assuming constant, uniform wave conditions and that the angle between wave crests and the shoreline contour is small, the equation for the position of the shoreline (or a single depth contour), y(x, t), from a fixed datum line takes the form of a linear diffusion equation with constant coefficient and has gained the epithet of the 'one-line equation': In Eq. (1) y is the distance of the reference contour from a datum line (usually taken to be the x-axis), t is time, and K is a diffusion coefficient that represents the factors affecting the rate at which sediment is transported along the shoreline, and may be written as 2Q 0 /D where D is the height of the active profile and Q 0 is a nominal sediment transport rate that depends upon wave height and sediment characteristics. The physical transport rate is equal to the product of Q 0 and a function of wave angle, determined empirically to be sin(2α b ) where α b is the angle between the breaking wave crests and the local shoreline contour (US Army Corps of Engineers, 2002). Specifically, where p s is the porosity of the beach material, ρ s is the density of the beach sediment (kg/m 3 ), ρ is the density of seawater (kg/m 3 ), H sb is the significant wave height at breaking (m), c gb is the wave group velocity at breaking (m/s), and κ is a dimensionless coefficient which is a function of the particle size and has a recommended value of approximately 0.4 for sandy beaches. Some of the first analytical solutions to this equation for cases of coastal engineering interest were presented by Grijm (1960), under the assumption of constant uniform waves. This assumption can be relaxed, but complicates the analysis considerably. Larson et al. (1997) presented an equation for constant but non-uniform wave conditions corresponding to Eq. (1), whilst Reeve (2006), Zacharioudaki and Reeve (2008), Walton and Dean (2011) and Valsamidis et al. (2013) have presented solutions for cases where wave conditions are uniform but time varying. The development of analytical models not only answers a pedagogical need but also broadens the range of conditions for testing computational models. Whilst some of the methods required to derive analytical solutions can appear difficult and arcane, the solutions themselves have several attractive properties. First, they can usually be evaluated in a straightforward and efficient manner. Second, no time stepping is required as the solution simply needs to be evaluated at a selected time, rather than approached iteratively through multiple time steps. Finally, analytical solutions do not suffer from numerical stability in the same way that computational models can. The need to make simplifying assumptions does restrict analytical methods and for many practical applications numerical models are preferred, particularly as they can be linked to other models that can predict, for example, nearshore wave transformation processes.
In practical applications, the 1-line model is solved using a time marching numerical solution procedure to solve the continuity, sediment transport and wave angle equations simultaneously, commonly relaxing the 'small angle approximation' made to assist in deriving analytical solutions. Wang and Le Méhauté (1980) showed that the expression for sediment transport can, under the assumption of parallel depth contours, be cast in terms of deepwater wave conditions. Instabilities can occur for certain values of offshore wave angle, although the equation governing shoreline evolution no longer takes the form in Eq. (1). This idea was used by Ashton et al. (2001) to develop a numerical model that provided an explanation for the unusual coastal spit formations found in the Azov Sea, and by Falques (2003) in an analysis of the method of estimating the diffusion coefficient.
The 1-line model has proved remarkably robust, and over the past decade or so has been used extensively, being an element in current design guidelines (eg. US Army Corps, 2002). The use of the 1-line model in more complicated situations has driven the development and expansion of a simple morphological equation (Eq. (1)) into computational prediction suites. For example, such modelling suites include elements of wave prediction, nearshore wave transformation, modifications to allow for wave diffraction, longshore variations in wave angle and height, and variations in beach slope (eg Kraus, 1989, Dabees &Kamphuis, 1999).
One criticism of the 1-line formulation is that it is not able to deal explicitly with cross-shore sediment transport and hence changes in the cross-shore beach profile. Changes in beach slope and convexity that arise from the construction of coastal defences require something more than the 1-line model. Bakker (1969) proposed a 2-line theory which predicted the simultaneous evolution of two distinct contour lines and hence changes in beach slope. Perlin and Dean (1983) subsequently presented the theory for an N-line model together with associated numerical solution methods. With a few notable exceptions (Dabees and Kamphuis, 2001;Hanson and Larson, 2001;Shibutani et al., 2009), the N-line formulation has not found widespread application. The reasons for this are not entirely clear, although there have been some reports discussing the difficulty of obtaining consistent results in practical applications. For example, Hulsbergen et al. (1976) performed laboratory experiments and found that the 2-line model worked well for a long, straight beach with a simple pattern of longshore currents but was not as accurate for more complicated situations; whilst Hanson and Larson (2001) commented upon instances of numerical difficulties in certain applications.

The 2-line model
The stability of an equation, or system of equations, can be investigated by analysing the propensity of small perturbations to grow or decay over time. Perturbations are typically taken to have a sinusoidal form. Seeking wave solutions to Eq. (1) of the form y(x,t) = Yexp{i(kx − ωt)} leads to the requirement that ω = −iKk 2 , where ω is the angular frequency of the perturbations, i is the imaginary unit and k is the wave number of the perturbations in the position of the shoreline contour. This requirement indicates that any wave solutions will be damped over time. It is a matter of convention as to over what ranges the wavenumber and frequency are allowed to vary. Here, we take k as being non-negative whilst ω may take on positive or negative values, depending upon which direction the perturbation propagates along the shoreline.
Following Bakker (1969), a 2-line model can be formulated on the argument that the beach profile can be divided into two parts, a shallow part that extends offshore to a depth D 1 and a deeper part that extends out to the depth of closure, written as D 1 + D 2 (see Fig. 1).
Again the concept of an equilibrium profile is used but in this case it is defined by the position of two contours together with the position of the depth of closure. In Fig. 1, y 1 is the distance of a chosen contour (typically the equilibrium sea level), from a fixed datum line landward of the shoreline. This datum line is usually taken as the 'x-axis'. A second contour, a distance w from the x-axis at equilibrium, is denoted by y 2 . When the profile is in equilibrium, both y 1 and y 2 are zero. When either of the contours is not in their equilibrium position a cross-shore sediment transport is inferred. The cross-shore transport of sediment between the two contours must satisfy continuity and is assumed, following Bakker (1969), to be in a sense that relaxes the beach profile towards its equilibrium form. Hence, it may be written as q y = C y (y 1 − y 2 ) where C y is a dimensional transport coefficient. It may be noted that even though y 1 is chosen to be the contour of the equilibrium water level this does not imply that longshore transport is zero, because the underlying assumption of line models is that they are applied to beaches that are predominantly drift dominated. Similarly, at the depth of closure the assumption is that there is zero net cross-shore transport, although longshore transport may occur.
The equations for the 2-line model may be written (Dean and Dalrymple, 2002): where: y 1 is the distance of the first contour from the y-axis and y 2 is the distance of the second contour from its equilibrium position; q 1 and q 2 are the alongshore sediment transport rates along the first and second contours, respectively; t is time and x is the longshore distance. Eliminating y 2 from Eqs. (2) and (3), we have: where D = D 1 + D 2 ; q i is the longshore sediment transport rates corresponding to the contour y i and q = q 1 + q 2 . Seeking wave solutions as before with the following dispersion relation may be derived: By equating the coefficients of the imaginary terms on either side of Eq. (6), it is concluded that either ω r = 0 or the sum of the imaginary terms is equal to 0.

Stationary wave condition
In this case ω r = 0, and Eq. (6) simplifies to a quadratic. We denote the solutions by ω i1 and ω i2 such that: Dean and Dalrymple (2002) noted that longwave solutions (for which k ≪ 1), have ω i b 0 and hence decaying amplitudes. However, if the full expression is retained and the sediment transport rates q 1 and q 2 have opposite signs, in other words q 1 /q 2 b 0, the dispersion equation can yield one or both roots positive for fixed values of the parameter k, and hence an unstable mode. More specifically, if q 1 , q 2 have the same signs, the dispersion equation yields both roots negative, and thus decaying perturbations. If q 1 /q 2 b 0 and (C y D + (D 2 q 1 + D 1 q 2 )k 2 ) ≥ 0 both roots of the dispersion equation are negative. On the other hand, if q 1 /q 2 b 0 and (C y D + (D 2 q 1 + D 1 q 2 )k 2 ) ≤ 0, both roots of the dispersion equation will be positive and perturbations will be unstable. If the roots have opposite signs then the criteria has to be evaluated from the full expression. Further details are provided in Appendix A. Fig. 2 presents the values of ω i1 and ω i2 plotted against the ratio q 1 /q 2 and the wave number k, in two contour graphs respectively. The values of the parameters have been chosen arbitrarily and in this case are; D 1 = D 2 = 3.20 m, C y = 4.23 × 10 −6 (s −1 ). The value of C y has been chosen to be of the same order as used by other researchers (for instance Hulsbergen et al., 1976). In this particular case, ω i1 is always negative (see Fig. 2a), whilst ω i2 may become positive for particular combinations of k and q 1 /q 2 (see Fig. 2b). For small values of wavenumber, ω i2 is negative whatever value q 1 /q 2 takes. As wavenumber increases, ω i2 becomes positive over an increasing range of values of q 1 /q 2 . A wavenumber of k = 1 corresponds to a perturbation with a wavelength of~6 m, similar to the spacing of small beach cusps. Larger values of wavenumber correspond to smaller scale features, whilst smaller values of wavenumber correspond to larger scale features. For engineering applications values of k between 0 and 1 are likely to be of greatest interest. For computational model developers the wavenumber is perhaps of less importance than the possibility of the existence of growing perturbations (of any scale) that could cause the computation to become unstable.

Propagating waves
In this case ω r ≠ 0, and the dispersion equation gives the constraint on ω i as: When the numerator of the term on the right-hand side of Eq. (8) is positive, perturbations will be unstable; more specifically if This inequality can be valid only if the sediment transport rates q 1 and q 2 have opposite signs. As C y , D, D 1 , D 2 and k 2 are all positive, if the sediment transport rates q 1 and q 2 are both positive then ω i b 0 so perturbations decay and the system is stable.

The 3-line model
In a manner analogous to the 2-line model the 3-line model considers the movement of three contours: y 1 (the shoreline one), y 2 (the distance of the second contour from its equilibrium position) , and y 3 (the distance of the seawardmost contour from its equilibrium position); and the continuity equation for each, taking into account the cross-shore sand movements q y = C y (y 1 − y 2 ) between the first and second contours and q z = C z (y 2 − y 3 ) between the second and third contours (Fig. 3).
The equations take the following form: Eliminating y 2 and y 3 from Eqs. (10)-(12) yields a 6th order equation for y 1 (additional details are provided in Appendix B). The dispersion relation becomes: This equation may be analysed as before to determine the stability criteria for perturbations to the beach.
For stationary waves (ω r = 0), the dispersion relation gives a cubic equation for ω i : where the coefficients are given in full in Appendix B. The values of ω i , corresponding to the roots of this cubic polynomial, as a function of wave number k, are presented in Fig. 4 for two different combinations of transport rates. The transport rates have been set arbitrarily to illustrate the stability properties of the equations. In practice, transport rates are rarely measured directly and are more often inferred from beach volume changes or calculated using transport formulae such as that due to US Army Corps (2002). In the 'standard' situation ( Fig. 4a) with alongshore transport rate increasing monotonically towards the shore, the solutions are all stable. The rate of 0.25 m 3 /s is approximately 8 million m 3 /year which would be at the upper end of reported values. For example, Taborda et al. (1994) estimated transport rates of almost 10 million m 3 /year on Portuguese beaches exposed to a breaking wave height of~3 m. In the second case (Fig. 4b) more modest transport rates are chosen, but with a reversal of transport direction. Such reversals are unlikely to occur on long, straight beaches with a simple lower beach structure. However, in cases of more complex bathymetry or in the presence of structures such as detached breakwaters a reversal of transport direction is certainly conceivable. From the stability analysis above it is clear that instability can occur where there is a reversal in the transport rates. (For any particular combination of parameter values the number of positive roots can be determined using Descartes' Rule of Signs.) Fig. 2. a. ω i1 as a function of wavenumber (k) and transport rate ratio. b. ω i2 as a function of wavenumber (k) and transport ratio. Positive values are shaded in red. For propagating wave solutions (ω r ≠ 0), the following expression for the frequency applies: Equating the coefficients of the imaginary terms on either side of Eq. (15) yields: which is a quadratic for ω i , with roots given by: Eq. (16) indicates that ω i has a dependency on ω r , however, for ω r ≪ 1, this dependence is weak.
The role of the factor c − aω r 2 in Eq. (17) in determining the stability of the system is analysed further below. We consider the cases where the factor is positive or negative separately, noting that the coefficient a must be positive by definition.
If c − aω r 2 b 0: Then, for ω i1 ; if b N 0 we obtain ω i1 N 0 as the term inside the square root is greater than 2b; if b b 0 we obtain ω i1 N 0 as both terms in the numerator are positive.
For ω i2 ; if b N 0 we obtain ω i2 b 0 as both terms in the numerator are negative; if b N 0 we obtain ω i2 b 0 as the term inside the square root is greater than 2b.
If c − aω r 2 N 0: Then, for ω i1 ; if b N 0 we obtain ω i1 b 0 as the term in the square root is less than 2b; if b b 0 we obtain ω i1 N 0 as the numerator is positive. Consequently, if c − aω r 2 b 0 unstable modes are possible if b b 0, which can only occur if q 1 , q 2 and q 3 are not all the same sign.
For ω i2 ; if b N 0 we obtain ω i2 b 0 as the numerator will be negative; if b b 0 we obtain ω i2 N 0 as the square root term is less than 2b.
Finally, if c − aω r 2 = 0, then ω i1 = 0 and ω i2 b 0 if b N 0, and ω i1 N 0 and ω i2 = 0 if b b 0. Fig. 5 shows the values of ω i (as determined from Eq. (13)), as a function of the wave number for various choices of transport rates.
In the 'standard' situation, where the transport rates are of the same sign and increase as the contours are closer to the shoreline, there is one stable mode and one unstable mode (see Fig. 5a). The instability has quite different characteristics to the stationary wave instability. The maximum growth rate occurs at the largest scales (|k| ≪ 1), and the growth rate has a finite maximum. A similar situation holds when the transport rate along the middle contour is halved, as is evident in Fig. 5b. For the situation with transport reversal (Fig. 5c) the small wavenumber instability remains but as a local maximum. Significantly, the high wavenumber instability is present. For large k, this has much greater growth rates than the large scale instability and would be expected to dominate the shoreline evolution.

Conclusions
In this paper the stability of a class of beach evolution models is investigated. Specifically, the line models that predict the evolving plan shape of one or more depth contours. Such models are used widely in the coastal engineering community, both for research and design. Multiple line models are, in principle, able to predict changes in beach profile slope and convexity. However, the 1-line model has proved the most robust and popular for applications, perhaps due to its stability to small perturbations in beach position.
For the case of two and three line models we establish that instabilities can exist. In the 2-line case, the system is generally stable but can exhibit instability if the sediment transport rates along each contour are in opposition. This is not very likely to occur on open straight coasts but could quite conceivably occur in the presence of coastal structures such as groynes and detached breakwaters. The instability is associated with small spatial scales of little coastal engineering interest. However its growth rate increases quadratically with wavenumber suggesting it would be difficult to control in computational models once triggered.
The 3-line system has a similar instability; a sufficient condition being that the sediment transport along one of the 3 contours is in opposition to the transport direction of the other two. Significantly, this system also has a large scale propagating instability, which does not require transport reversal. Its growth rate is finite and greatest for small spatial frequencies. This could manifest itself at scales of coastal engineering interest. However, due to its finite growth rate it is possible Fig. 4. a. ω i plotted as a function of wavenumber for longshore sediment transport rates q 1 = 0.25 m 3 /s, q 2 = 0.20 m 3 /s, q 3 = 0.15 m 3 /s. b. ω i plotted as a function of wavenumber for longshore sediment transport rates q 1 = 0.0025 m 3 /s, q 2 = −0.0025 m 3 /s, q 3 = −0.00125 m 3 /s. Fig. 5. a. ω i plotted as a function of wave number for longshore sediment transport rates q 1 = 0.25 m 3 /s, q 2 = 0.20 m 3 /s, q 3 = 0.15 m 3 /s. b. ω i plotted as a function of wave number for longshore sediment transport rates q 1 = 0.25 m 3 /s, q 2 = 0.10 m 3 /s, q 3 = 0.15 m 3 /s. c. ω i plotted as a function of wave number for longshore sediment transport rates q 1 = 1.55 m 3 /s, q 2 = −0.25 m 3 /s, q 3 = −0.0125 m 3 /s. that in computational models such modes could be controlled either with additional damping terms or a dissipative numerical time-stepping scheme.
The analysis presented in this paper provides some explanation for the difficulties encountered in implementing 2-line and 3-line computational models in all but the simplest of beach geometries. In situations where the longshore transport direction changes sign along the profile the potential for unstable oscillations to develop is likely to limit the use of such models. It should also be noted that in the treatment above we have considered longshore quasi-homogenous conditions, in which the transport rates are independent of longshore position. In reality, transport reversal may be localised, and further investigation is required to determine whether local instabilities would be similarly localised or would propagate throughout the domain of interest. The existence of these instabilities provides a possible explanation for some of the difficulties encountered with multi-line models reported in the literature. The analysis provides some guidance to model developers. Specifically, if transport reversal occurs then numerical instability can be expected. Further, the representation of cross-shore transport as a relaxation process towards equilibrium leads to the specific form of equation which supports unstable modes. This element of multi-line models may need revision in order to develop stable computational procedures for dealing with more realistic conditions.
There are three cases to consider: both roots negative; both roots positive; and one root of each sign. a. If both roots of the dispersion equation are negative: Combining Eqs. (A1) and (A2), and after some rearrangements we get condition (A3): As C y , D 1 and D 2 are positive, if q 1 , q 2 N 0 then the condition (A3) is met if q 1 , q 2 N 0. If q 1 N 0 and q 2 b 0, and then Eq. (A3) also holds. However, if D 2 q 1 + D 1 q 2 b 0 then in order for Eq. (A3) to be valid the condition that must be satisfied is a more general one: b. If both roots of the dispersion equation are positive: and ω i2 ≥ 0⇔ − C y D þ D 2 q 1 þ D 1 q 2 ð Þ k 2 − ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C y D þ D 2 q 1 þ D 1 q 2 ð Þ k 2 2 −4D 1 D 2 C y q þ q 1 q 2 k 2 k 2 r 2D 1 D 2 ≥0: Adding Eqs. (A6) to (A7), and after some rearrangements, we get condition (A8); C y D þ D 2 q 1 þ D 1 q 2 ð Þ k 2 ≤0⇔D 2 q 1 þ D 1 q 2 ≤− C y D k 2 : If q 1 ,q 2 N 0, then the inequality Eq. (A8) cannot be met. Necessarily, q 1 /q 2 b 0. c. Similarly, in the case of one positive and one negative root simplification is not possible and the condition must be determined from the sign of the numerator: − C y D þ D 2 q 1 þ D 1 q 2 ð Þ k 2 AE ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C y D þ D 2 q 1 þ D 1 q 2 ð Þ k 2 2 −4D 1 D 2 C y q þ q 1 q 2 k 2 k 2 r : If this is positive then unstable modes can exist.