RHYTHMIC SHAPES CAUSED BY SHORELINE INSTABILITY USING BG MODEL

The BG model (a three-dimensional model for predicting beach changes based on Bagnold’s concept) was used to simulate the shoreline evolution caused by the high-angle wave instability discussed by Ashton et al. Three calculations were carried out: the wave direction was assumed to be obliquely incident from 60 ̊ counterclockwise (Case 1) or from the directions of ±60 ̊ with probabilities of 0.5:0.5 (Case 2) and 0.65:0.35 (Case 3), while determining the incident wave direction from the probability distribution at each step. The three-dimensional development of multiple sand spits and cuspate forelands with rhythmic shapes was successfully explained using the BG model. The results of the previous study conducted by Ashton et al. were reconfirmed and reinforced.


INTRODUCTION
showed that multiple sand spits with rhythmic shapes may develop in a shallow water body such as the Azov Sea and a lagoon facing the Chukchi Sea, as shown in Fig. 1, and concluded that under oblique wave incidence with the angle between the direction normal to the shoreline and the wave direction being larger than 45˚, shoreline instability may develop, and during the development of sand spits, the wave-sheltering effect due to the sand spits themselves plays an important role.Ashton et al. (2001) adopted this mechanism in their model and successfully modeled this shoreline instability using the upwind scheme in their finite difference method to prevent the numerical instability on the basis of the conventional longshore sand transport formula.They also evaluated the wave-sheltering effect using the angular spreading method (Ashton and Murray, 2006).(Zenkovich, 1967).
In this study, we used the BG model (a three-dimensional model for predicting beach changes based on Bagnold's concept) developed by Serizawa and Uda (2011) to simulate the shoreline evolution caused by high-angle wave instability (Ashton and Murray, 2006) and showed that the threedimensional (3-D) beach changes associated with the shoreline instability can be explained using this model.In this case, although the fundamental equation of the BG model defined by Serizawa and Uda (2011) was employed as the sand transport equation, the sand transport flux was assumed to be proportional to the wave energy dissipation rate instead of the third power of the amplitude of the bottom oscillatory velocity due to waves, and the wave energy dissipation rate was given by that due to wave breaking at each point determined in the calculation of the wave field.
This sand transport equation was derived by applying the concept of the equilibrium slope introduced by Inman and Bagnold (1963) and the energetics approach of Bagnold (1963).This sand transport equation is equivalent to the CERC-type longshore sand transport formula (Komar and Inman, 1970) when the total longshore sand transport Q is calculated by integrating the longshore component of the sand transport flux in the cross-shore direction, as shown by Serizawa et al. (2006).Falqués et al. (2008) also predicted the development of sand waves caused by high-angle wave instability using equations similar to that of our model, but not the development of sand spits protruding offshore.In this study, the development of sand spits and cuspate forelands was predicted under oblique wave incidence with the angle between the direction normal to the shoreline and the wave direction being larger than 45˚, given a small perturbation in the initial topography.

NUMERICAL MODEL (BG MODEL)
We use Cartesian coordinates (x, y) and consider that the elevation at a point Z (x, y, t) is a variable to be solved, where t is the time.The beach changes are assumed to occur between the depth of closure h c and the berm height h R .The BG model proposed by Serizawa and Uda (2011) was employed to predict the beach changes.An additional term given by Ozasa and Brampton (1980) was incorporated into the fundamental equation of the BG model to evaluate the longshore sand transport due to the effect of the longshore gradient of wave height.The fundamental equation is given as follows. € Here, € q = q x ,q y ( ) is the net sand transport flux, Z (x, y, t) is the elevation, n and s are the local coordinates taken along the directions normal (shoreward) and parallel to the contour lines, respectively, € ∇Z = ∂Z ∂x, ∂Z ∂y ( ) is the slope vector, € ew is the unit vector of the wave direction, € es is the unit vector parallel to the contour lines, α is the angle between the wave direction and the direction normal to the contour lines, € tanβ = ∇Z is the seabed slope, tanβ c is the equilibrium slope, and € tanβ es = −∂Z ∂y, ∂Z ∂x ( ) .Moreover, K s and K n are the coefficients of longshore and cross-shore sand transport, respectively, K 2 is the coefficient of the term given by Ozasa and Brampton (1980), € ∂H ∂s = e s ⋅∇H is the longshore gradient of the wave height H measured parallel to the contour lines, and € tanβ is the characteristic slope of the breaker zone.In addition, C 0 is the coefficient transforming the immersed weight expression into a volumetric expression ( , where ρ is the density of seawater, ρ s is the specific gravity of sand particles, p is the porosity of sand and g is the acceleration due to gravity), h c is the depth of closure, and h R is the berm height.
The intensity of sand transport P in Eq. ( 1) is assumed to be proportional to the wave energy dissipation rate, as described by Serizawa et al. (2006), on the basis of the energetics approach of Bagnold (1963).P is given by Φ all (Eq.( 2)) in accordance with the BG model proposed by Serizawa et al. (2006), in which the intensity of sand transport is proportional to the wave energy dissipation rate due to wave breaking, instead of the assumption that it is proportional to the third power of the amplitude of the bottom oscillatory velocity u m due to waves, as described in Serizawa et al. (2011).
For the calculation of the wave field, the numerical simulation method using the energy balance equation given by Mase (2001), in which the directional spectrum of irregular waves is the variable to be solved, was employed with an additional term of energy dissipation due to wave breaking proposed by Dally et al. (1984), similarly to that reported by Serizawa and Uda (2011).Φ all in Eq. ( 2) was calculated from Eq. (3), which defines the total sum of the energy dissipation of each component wave due to breaking.
Here, f D is the energy dissipation rate, E is the wave energy, K is a coefficient expressing the intensity of wave dissipation due to breaking, h is the water depth, Γ is the ratio of the critical wave height to the water depth on a flat bottom, and γ is the ratio of wave height to the water depth H/h.In addition, a lower limit was set for the water depth h in Eq. ( 3) similarly in Serizawa and Uda (2011).
In this method, the energy dissipation rate obtained from the calculation of the plane wave field including the effect of wave dissipation due to breaking was used for the calculation of sand transport.
The same approach was employed by Katayama and Goda (2002).In the calculation of the wave field in the wave run-up zone, an imaginary depth was assumed as in the case of Serizawa and Uda (2011).Furthermore, the wave energy at locations with elevations higher than the berm height was set to 0.
In the numerical simulation of beach changes, the sand transport and continuity equations ( € ∂Z ∂t+∇•q =0 ) were solved on the x-y plane by the explicit finite-difference method using the staggered mesh scheme.In estimating the intensity of sand transport near the berm top and at the depth of closure, the intensity of sand transport was linearly reduced to 0 near the berm height or the depth of closure to prevent sand from being deposited in the zone higher than the berm height and the beach from being eroded in the zone deeper than the depth of closure, similar to in Serizawa et al. (2003).

CALCULATION CONDITIONS
As the wave conditions, considering the formation of sand spits in a shallow lagoon, we assumed H i = 1 m and T = 4 s.The wave direction was assumed to be obliquely incident from 60˚ counterclockwise (Case 1) or from the directions of ±60˚ with probabilities of 0.5:0.5 (Case 2) and 0.65:0.35(Case 3), while determining the direction from the probability distribution at each step.We considered a shallow lake with a flat solid bed, the depth of which was given by Z = -4 m, and a uniform beach with a slope of 1/20 and a berm height of h R = 1 m were considered on the landward end.At the initial stage, a small random perturbation with an amplitude of ΔZ = 0.5 m was applied to the slope.The calculation domain was a rectangle of 4 km length and 1.2 km width, and a periodic boundary condition was set at both ends.In addition, the depth of closure was assumed to be h c = 4 m.The equilibrium and repose slopes were 1/20 and 1/2, respectively.The coefficients of longshore and cross-shore sand transport were set to be K s = K n = 0.2, respectively.The calculation domain was divided with a mesh size of Δx = Δy = 20 m, and Δt was selected to be 0.5 hr.The total number of calculation steps considered was 5.5 × 10 4 (2.75 × 10 4 hr).The calculation of the wave field was carried out every 10 steps in the calculation of beach changes.Table 1 shows the calculation conditions.

CALCULATION RESULTS
Oblique wave incidence from 60° counterclockwise (Case 1) Figure 2 shows the results of the calculations at eight stages starting from the initial straight shoreline with a slope of 1/20, to which a small random perturbation with an amplitude of ΔZ = 0.5 m was applied, up to 5.5 × 10 4 steps.The small perturbation applied to the slope at the initial stage developed into eleven cuspate forelands within 5 × 10 3 steps and the shoreline projection increased with time while moving rightward owing to the wave incidence from the counterclockwise direction.Because of the periodic boundary conditions at both ends, the cuspate forelands that moved away through the right boundary reentered the calculation domain through the left boundary.After 1 × 10 4 steps, the shoreline protrusion had increased at the tip of the cuspate forelands and had developed as slender sand spits.After 2 × 10 4 steps, the small-scale sand spits located adjacent to each other had merged into large-scale sand spits and disappeared, and finally six sand spits were formed.
Two reasons for these changes are considered (Zenkovich, 1967;Ashton and Murray, 2006).( 1) Of the two sand spits of different scales, the small sand spit moves faster than the large sand spit in the absence of the wave-sheltering effect, and then the small sand spit catches up and merges with the large sand spit.(2) On the lee of sand spits with an elongated neck, a wave-shelter zone is formed and the velocity of sand spits is reduced in this zone because of wave calmness, resulting in the stoppage of the movement of the sand spits and in the merging of small sand spits with larger spits.Furthermore, the sand spits developed and protruded because (1) their tip is semicircular, meaning that the angle between the direction normal to the shoreline and the wave incident direction exceeds 45˚ at a point along the shoreline and the shoreline protrusion occurs at such a point owing to high-angle wave instability.(2) In a wave-shelter zone, sand transport is significantly reduced, whereas it is enhanced near the tip of the sand spits, and thus the derivative of the sand transport rate takes a maximum value near the boundary between the tip of the sand spits and the wave-shelter zone, inducing the protrusion of sand spits.
After 3 × 10 4 steps, the small sand spits located in the wave-shelter zone of the large-scale sand spits had stopped moving and merged into the large-scale sand spits, resulting in an increase in the interval between the sand spits and a decrease in the number of sand spits per finite length of the shoreline.After 4 × 10 4 steps, the number of sand spits had decreased to 2 and the tip of the sand spits approached closely to the original shoreline, permitting the downcoast passage of the sand of the sand spits.
After 5 × 10 4 steps, because of the movement of sand spits sweeping rightward, part of the sand deposition zone immediately offshore of the shoreline was left intact at the base of the sand spit located at y = 2350 m but almost all parts had merged with the sand spits.Although two large-scale sand spits were formed from the straight shoreline within 5.5 × 10 4 steps, the offshore contour of -4 m depth obliquely extended and a gentle seabed slope was formed upcoast of the sand spit, whereas a very steep slope was formed at the tip of the sand spits.These features are in good agreement with those measured around sand spits in lake and bay (Uda and Yamamoto, 1991).At the downcoast base of the sand spit extending from y = 2400 m, part of the sand bar formed in the previous process from 4 × 10 4 steps was left intact, implying that historical changes could be recovered at the downcoast side of the sand spit on the basis of the present topography.
Figure 3 shows the sand transport flux after 5.5 × 10 4 steps.The longshore sand transport mainly develops along the outer margin of the sand spit with a maximum value at the tip of sand spit and then rapidly decreases.Examining the sand transport flux at the neck of the sand spit in Fig. 3, cross-shore sand transport flux from the exposed side to the lee of the sand spit is also observed at a location of y = 2750 m.Thus, the neck of the sand spits is gradually eroded and moves downcoast because of this cross-shore sand transport, caused by the small difference between the given berm height h R and the actual crown height of the sandy beach comprising the neck.Sand can be directly transported from the exposed side to the lee side without the sand transport turning around the tip of the sand spits.This effect makes the movement of an entire sand spit possible.Figure 4 shows the change in the wave field around the sand spits at each stage from the initial straight shoreline to the fully developed sand spits as shown in Fig. 2. At the initial stage, waves are obliquely incident to the straight shoreline with uniform exposure to waves at all locations.With the development of the shoreline undulation over time, wave-shelter zones were formed behind them.After 2 × 10 4 steps, the formation of sand spits was clear and the wave-shelter zone expanded downcoast and the toe of the adjacent sand spit was included inside the wave-shelter zone.As a result of the inclusion of the sand spit into the wave-shelter zone formed by the upcoast sand spit, a marked reduction in wave height occurred, which in turn caused a reduction in sand transport.After 5.5 × 10 4 steps, long sand spits extended so that the toe of the slender sand spit was subject to the wave-sheltering effect of the upcoast sand spit.Owing to the development of sand spits over time, the entire original shoreline zone was included in the wave-shelter zone produced by sand spits, which is very similar to the calm wave zone protected by the extension of long port breakwaters.
The calculations of beach changes were also carried out under the conditions of oblique wave incidence with an angle of 50˚ and 40˚, as shown in Figs. 5 and 6.The sand spits developed owing to the instability mechanism under the condition of oblique wave incidence with an angle of 50˚, but no sand spits had developed at an angle of 40˚ and the shoreline undulations were smoothed out with time, resulting that the shoreline undulations did not develop unless the wave incidence to the zone shallower than h c exceeds 45˚.This result agrees with the conclusion by van den Berg et al. ( 2011) that the shoreline instability develops only if the bathymetric changes related to shoreline perturbations extend to a depth where the wave angle is greater than the critical angle of 42˚ and the potential for coastline instability is therefore limited by the wave incidence angle at the depth of closure and not the angle at deep water.A numerical simulation was carried out for the case that waves were obliquely incident from the directions of ±60˚ relative to the direction normal to the shoreline with probabilities of 0.5:0.5 on the initial straight coastline given a small random perturbation at the initial stage.Figure 7 shows the bathymetric changes between the initial stage and 4 × 10 4 steps.After 1 × 10 4 steps, approximately triangular cuspate forelands had developed and irregularly distributed.When waves were incident from one direction, asymmetric sand bars developed, as shown in Fig. 2. In contrast, when waves with the same probability were incident from the opposite direction, the cuspate forelands became symmetric, and small-scale cuspate forelands disappeared and merged with larger cuspate forelands.Because the probability of occurrence of both wave directions was the same and there was no net longshore sand transport, the unidirectional movement of the sand body did not occur, as in Case 1.
After 2 × 10 4 steps, the number of triangular cuspate forelands had been reduced to five, and the development of triangular cuspate forelands further continued and small-scale cuspate forelands merged into larger cuspate forelands.Finally, after 4 × 10 4 steps, four large-scale cuspate forelands had developed.Figures 8(a) and 8(b) show the wave height distribution around the cuspate forelands immediately after 4 × 10 4 steps under the conditions of clockwise and counterclockwise wave incidence.The wavesheltering effect due to the protruded cuspate forelands alternately extends to the bays.The importance of this effect to the development of shoreline undulations was pointed out by Ashton and Murray (2006); when the multiple cuspate forelands with a different size have developed, the effect of the highangle wave instability becomes stronger at the tip of the forelands with a large size than that at the bays, so that the positive feedback will occur.In contrast, around the cuspate forelands with a small size, the effect of the high-angle wave instability is weakened by the wave-sheltering effect by the large cuspate forelands and the cuspate forelands are gradually modified to a stable form.Furthermore, when a large cuspate foreland develops, sand composed of the small cuspate foreland is absorbed into the large forelands, resulting in the decline of the small cuspate forelands and the increase in size of the large cuspate forelands.Thus, the development of large cuspate forelands will continue while small cuspate forelands are gradually disappearing.This results in the decrease in the number of the cuspate forelands.Figure 9 shows the calculation results obtained when waves were incident from the directions of ±60˚ with probabilities of 0.65:0.35.Because the direction of net longshore sand transport was rightward, sand spits developed while moving rightward.After 2 × 10 4 steps, slender sand spits extended offshore, and a hooked shoreline had formed downcoast of the spits.After 3 × 10 4 steps, these sand spits extended further rightward with a narrow neck at the connecting point to the land.In addition, after 4 × 10 4 steps, sand spits obliquely extended rightward from the tip of the cuspate forelands with a larger angle than that in Fig. 2. In particular, the calculation results obtained after 4 × 10 4 steps and the case in a lagoon facing the Chukchi Sea shown in Fig. 1 are in good agreement.
Thus, the mechanism based on the high-angle wave instability and the evolution of threedimensional beach changes was explained well by the BG model.Using this model, not only the shoreline configuration but also the 3-D topographic changes around the sand spits and cuspate forelands can be predicted.The shoreline evolution caused by the high-angle wave instability was successfully simulated by the present model.

CONCLUSION
Regarding the development of multiple sand spits and cuspate forelands with rhythmic shapes owing to the instability mechanism along the shore of the Azov Sea (Zenkovich, 1967), Ashton et al. (2001) successfully developed a numerical model that can predict this shoreline evolution on the basis of the longshore sand transport formula of the one-line model.In this study, the BG model developed by Serizawa and Uda (2011) was employed to simulate the shoreline evolution caused by high-angle wave instability in three cases: the wave direction was assumed to be obliquely incident from 60˚ counterclockwise (Case 1) or from the directions of ±60˚ with probabilities of 0.5:0.5 (Case 2) and 0.65:0.35(Case 3), while determining the direction from the probability distribution at each step.As a result, the 3-D development of multiple sand spits and cuspate forelands with rhythmic shapes was successfully explained by the present model and the results of the previous study conducted by Ashton et al. (2001) were reconfirmed and reinforced.
Because the wave field was predicted using the energy balance equation for irregular waves in this study, the wave field including wave refraction, wave breaking and the wave-sheltering effect can be systematically predicted.In addition, because two-dimensional sand transport equations were employed in our model, in contrast to the model developed by Ashton et al. (2001) in which the longshore sand transport formula was used, this model has the advantages of the conventional 3-D model for predicting beach changes for various applications.

Figure 1 .
Figure 1.Multiple sand spits with rhythmic shapes developed in Azov-type shallow water body facing the Chukchi Sea in Russia (Zenkovich, 1967).
Wave conditionsIncident waves: HI = 1 m, T = 4 s, wave direction θI = 60° relative to direction normal to initial shoreline Berm height hR = 1 m Depth of closure hc = 4 m (still water depth) Equilibrium slope tanβc = 1/20 Angle of repose slope tanβg = 1/2 Coefficients of sand transport Coefficient of longshore sand transport Ks = 0.2 Coefficient of Ozasa and Brampton (1980) term K2 = 1.62KsCoefficient of cross-shore sand transport Kn = Ks Mesh size Δx = Δy = 20 m Time intervals Δt = 0.5 hr Duration of calculation 2.75 × 10 4 hr (5.5 × 10 4 steps) Boundary conditions Shoreward and landward ends: qx = 0, right and left boundaries: periodic boundary Calculation of wave field Energy balance equation (Mase, 2001) •Term of wave dissipation due to wave breaking: Dally et al. (1984) model •Wave spectrum of incident waves: directional wave spectrum density obtained by Goda (1985) •Total number of frequency components NF = 1 and number of directional subdivisions N θ = 8 •Directional spreading parameter Smax = 25 •Coefficient of wave breaking K = 0.17 and Γ = 0.3 •Imaginary depth between depth h0 (3 m) and berm height hR •Wave energy = 0 where Z ≥ hR •Lower limit of h in terms of wave decay due to wave breaking Φ : 0.5 m

Figure 4 .
Figure 4. Change in wave field during full development of sand spits.

Figure 5 .
Figure 5. Formation of sand spits under the condition of oblique wave incidence from +50°.

Figure 6 .
Figure 6.No shoreline instability under the condition of oblique wave incidence from +40°.