Eight types of BG models and their applications.

## Abstract

Eight types of the BG models are introduced in this chapter. The Type 1 is a model using wave parameters at the breaking point. In the Type 2, the effect of longshore sand transport due to the effect of the longshore gradient of breaker height is included with an additional term given by Ozasa and Brampton. In the Type 3, the intensity of sand transport P is assumed to be proportional to the third power of the amplitude of the bottom oscillatory velocity um due to waves, and in the Type 4, P is given by the wave energy dissipation rate due to wave breaking at a local point. In the Type 5, wave power is calculated using the coordinate system different from that for the calculation of beach changes to predict the topographic changes of an island or a cuspate foreland in a shallow water body under the action of waves randomly incident from every direction. In the Type 6, the height of wind waves is predicted using Wilson’s formula using the wind fetch distance and wind velocity, and then sand transport fluxes are calculated. The Type 7 is a model for predicting the formation of the ebb-tidal delta under the combined effect of waves and ebb-tidal currents with an analogy of the velocity distribution of ebb-tidal currents to the wave diffraction coefficient, which can be calculated by the angular spreading method for irregular waves. In the Type 8, the effect of the nearshore currents induced by forced wave breaking is incorporated into the model by calculating the nearshore currents, taking both the wave field and the current velocity at a local point into account.

### Keywords

- eight types of BG models
- discretization method

## 1. Introduction

Eight types of the BG models to be used in the following chapters are introduced (Table 1). The Type 1 is a model using wave parameters at the breaking point. In the Type 2, the effect of longshore sand transport due to the effect of the longshore gradient of breaker height is included with an additional term given by Ozasa and Brampton [1], whereas the sand transport equations with the coefficients of longshore and cross-shore sand transport are employed. In the Type 3, the intensity of sand transport, *P*, is assumed to be proportional to the third power of the amplitude of the bottom oscillatory velocity, *u*_{m}, due to waves predicted in the calculation of the plane wave field. In the Type 4, *P* is given by the wave energy dissipation rate due to wave breaking at a local point. To calculate the wave field, a numerical simulation method using the energy balance equation is employed. In the Type 5, wave power is calculated using the coordinate system different from that for the calculation of beach changes in order to predict the topographic changes of an island or a cuspate foreland in a shallow water body under the action of waves randomly incident from every direction. In the Type 6, the height of wind waves is predicted using Wilson’s formula given the wind fetch distance and wind velocity, and the segmentation of a closed water body is predicted. The Type 7 is a model for predicting the formation of the ebb-tidal delta under the combined effect of waves and ebb-tidal currents with an analogy of the velocity distribution of ebb-tidal currents to the wave diffraction coefficient, which is calculated by the angular spreading method for irregular waves. In the Type 8, the effect of the nearshore currents induced by forced wave breaking is incorporated into the model by calculating the nearshore currents, taking both the wave field and the current velocity at a local point into account.

Type | External force | Characteristics | Application | Chapter |
---|---|---|---|---|

1 | Waves | Most fundamental model using wave parameters at breaking point | Prediction of typical beach changes in coastal engineering | 3 |

2 | Waves | The effect of longshore sand transport due to the effect of the longshore gradient of breaker height was included with an additional term given by Ozasa and Brampton [1] | Prediction of beach changes by human activities | 3 |

3 | Waves | The intensity of sand transport P is assumed to be proportional to the third power of the amplitude of the bottom oscillatory velocity, u_{m}, due to waves predicted in the calculation of the plane wave field | Formation of sand spit and bay barrier | 5 |

4 | Waves | The intensity of sand transport P is given by the wave energy dissipation rate due to wave breaking at a local point | Formation of cuspate foreland | 7 |

5 | Waves | The wave power is calculated using the coordinate system different from that for the calculation of beach changes under the action of waves randomly incident from every direction | Formation of sand spit and bay barrier and interaction of sandy islands on flat shallow seabed owing to waves | 5 and 6 |

6 | Wind waves developed in closed water body | The height of wind waves is predicted using Wilson’s formula given the wind fetch distance and wind velocity | Segmentation and merging of closed water bodies by wind waves | 8 |

7 | Waves and ebb-tidal currents | A model for predicting the formation of the ebb-tidal delta under the combined effect of waves and ebb-tidal currents with an analogy of the velocity distribution of ebb-tidal currents to the wave diffraction coefficient | Formation of dynamically stable ebb-tidal delta | 4 |

8 | Waves and nearshore currents | The effect of the nearshore currents induced by forced wave breaking is incorporated into the model, taking both the wave field and the current velocity at a local point into account | Beach changes around artificial reef | 4 |

## 2. Eight types of BG models

### 2.1 Type 1 BG model

A simple, practical Type 1 BG model is derived here using wave parameters at the breaking point on the basis of Eq. (26) in Chapter 1. The energy dissipation ratio *Φ* is evaluated by the energy flux of waves at the breaking point as in the contour-line change model [2]. Assume that waves are obliquely incident to a coast with a uniform slope of *Φ* is given by the wave energy flux transported per unit width of the shoreline, *h*_{c}, and the berm height, *h*_{R}, corresponding to the wave run-up height (Figure 1). Then, assuming *Φ* becomes *ε*(*Z*) and transforming *ε*(*Z*), the coefficient *A* of Eq. (23) in Chapter 1 is expressed in terms of the wave energy flux at the breaking point, and the sand transport flux is reduced to

Eq. (1) can be written in components as Eqs. (2a) and (2b), and *G* in Eq. (1) is expressed by Eqs. (3) and (4):

Here, *ε*(*Z*) is defined so that the integral over the depth between *Z* = −*h*_{c} and *h*_{R} is equal to 1, as given by Eq. (5), and the depth distribution of the longshore sand transport, Eq. (6), given by Uda and Kawano [3] is used. This distribution takes a peak value on the shoreline (*Z* = 0), and it monotonically decreases landward and seaward. In addition, when it is given by a uniform distribution, Eq. (7) can be used:

The wave energy flux at the breaking point (*ECg*)_{b} can be calculated by small-amplitude wave theory in shallow water using the breaker height *H*_{b}, where *H*_{b} is given by the following relation:

Here, *γ* is the ratio of the breaker height to the water depth. *k*_{1} = (4.004)^{2} in Eq. (8c) is a constant in the relationship between the wave energy *E* and the significant wave height when the probability of the wave height of irregular waves is assumed to be given by the Rayleigh distribution [4]. Furthermore, regarding the breaker angle *α*_{b} in Eq. (4), the following approximation is assumed:

Here, *α* is the angle between the wave direction at each point *θw* and the direction (shoreward positive) normal to the contours *θn* as expressed by Eq. (12) in Chapter 1. In a numerical simulation, the mean beach slope of the coast before the construction of the structures is given as the equilibrium slope to consider the long-term prediction of beach changes, similarly to the contour-line change model [2]. Here, the measured slope is employed as the equilibrium slope to make the prediction of beach changes on real coasts possible. The measured slope is given a priori because the real seabed topography includes every effect of past events, and it has a stable form, except for seasonal short-period variations, in the long term.

In the calculation, the coastal domain is discretized using 2D elements with widths Δ*x* and Δ*y*. The calculation points of the seabed elevation *Z* and sand transport rates

Given the initial seabed topography, the distribution of *H*_{b}, *α*_{b}, *h*_{c}, and *h*_{R}, and the equilibrium slope, sand transport fluxes and the change in seabed elevation after time Δ*t* are calculated using Eq. (2) and Eq. (28) in Chapter 1, respectively. These procedures are repeated recurrently. As the boundary conditions, sand transport is set to 0 along the outer boundaries of the calculation domain and the boundary along the structures. Here, a lower minimum of 0.5 was set for *P* value in Eq. (4) to avoid the occurrence of a local discontinuity in the topography.

In the calculation of beach changes in the wave-shelter zone behind an offshore breakwater, the wave diffraction coefficient *K*_{d} and the direction of diffracted waves *θ*_{d} are predicted using the angular spreading method for irregular waves [5, 6], and (*EC*_{g})_{b} in the case of no offshore structure is reduced by multiplying by the square of the coefficient *K*_{d} so that a relation of (*EC*_{g})_{b}′ = *K*_{d}^{2}(*EC*_{g})_{b} is satisfied. On the other hand, the direction of the diffracted waves *θ*_{d} is used as the wave direction *θw* to evaluate the wave energy flux at any point. In addition, the effect of the longshore sand transport induced by a longshore change in the breaker height as given by Ozasa and Brampton [1] can be included using the method in [5], although the term expressing this effect is not included in the sand transport equation of the Type 1 BG model. In the method, the wave direction, *θw*, in the wave-shelter zone is modified in response to the ratio of the coefficients of longshore sand transport, *K*_{2}/*K*_{1}, and the longshore gradient of the wave height *H*_{b} multiplied by *K*_{d}. When modified wave directions are employed in the sand transport equation, the additional effect given by Ozasa and Brampton [1] in the CERC-type longshore sand transport formula can be implicitly included [5].

### 2.2 Type 2 BG model

In the Type 1 BG model, the fundamental equation for sand transport is expressed in terms of the wave energy at the breaking point. In the Type 2 BG model, the equation in the Type 1 is improved to an equation including the coefficients of longshore and cross-shore sand transports and the additional term given by Ozasa and Brampton [1] to evaluate the effect of the longshore gradient of the breaker height.

We use the Cartesian coordinates (*x*, *y*), in which the *x*- and *y*-axes are taken in the cross-shore (shoreward positive) and longshore directions, respectively. For the sand transport equation, Eq. (10), expressed in terms of the wave energy at the breaking point, is used with the variables given by Eq. (11) together with Eq. (4):

Here, *θw* is the wave angle measured counterclockwise from the *x*-axis. *Kx* and *Ky* are the coefficients of cross-shore and longshore sand transport, respectively, *K*_{2} is the coefficient of the term given by Ozasa and Brampton [1], and *ε*(*Z*) is the depth distribution of the intensity of longshore sand transport, as defined by Eq. (6) or (7) in the Type 1. Beach changes are obtained by solving the continuity equation. The flowchart for numerical simulation using the Type 2 is the same as that in the Type 1, as shown in Figure 2.

Note that the *x*- and *y*-axes must be taken shoreward from an offshore point and in longshore direction, respectively, although the sand transport equations of the Type 1 do not depend on the directions of the orthogonal coordinates (*x*, *y*). Regarding *h*_{R} and *h*_{c}, the *K*_{d} value is also multiplied to cover the reduction in wave height if necessary.

### 2.3 Type 3 BG model

The Type 3 BG model is employed in the calculation of the formation of a sand spit and a bay barrier to be described in Chapter 5. Although the sand transport equations in the Types 1 and 2 BG model are expressed by using the wave energy flux at the breaking point, the sand transport equation in the Type 3 is expressed using local wave parameters which can be evaluated from the calculation of the plane wave field. Because a sand spit and cuspate foreland may change their configuration significantly, their effects on the wave refraction, wave breaking, and the wave-sheltering effect must be calculated by the recurrent calculations to evaluate the time changes in the wave field. The additional term given by Ozasa and Brampton [1] is also incorporated into the fundamental equation of the BG model to evaluate the longshore sand transport due to the effect of the longshore gradient of the wave height with the inclusion of two coefficients to evaluate the cross-shore and longshore sand transports. The fundamental equation is given by

Here, *Ks* and *Kn* are the coefficients of longshore sand transport and cross-shore sand transport, respectively, *H* measured parallel to the contour lines, and *h* is the water depth. Furthermore, *um* is the amplitude of the seabed velocity due to the orbital motion of waves given by Eq. (14).

The intensity of sand transport *P* in Eq. (12) is assumed to be proportional to the wave energy dissipation rate *Φ* based on the energetics approach of Bagnold [7]. In the Type 1, *P* is formulated using the wave energy at the breaking point, but here it is combined with the wave characteristics at a local point. Bailard and Inman [8] used the relationship *Φt = τut = ρCfut3* for the instantaneous wave energy dissipation rate *Φt* to derive their sand transport equation, where *τ* is the bottom shear stress, *ut* is the instantaneous velocity, and *Cf* is the drag coefficient. We basically follow their study but assume that *Φ* is proportional to the third power of the amplitude of the bottom oscillatory velocity *um* due to waves instead of the third power of the instantaneous velocity. The intensity of sand transport *P* is then given by Eq. (13), and its coefficient is assumed to be included in the coefficients of longshore and cross-shore sand transports, *Ks* and *Kn* in Eq. (12), respectively. *um* can be calculated by small-amplitude wave theory in shallow water using the wave height *H* at a local point in Eq. (14), which can be obtained by the numerical calculation of the plane wave field. *h*_{c} is assumed to be proportional to the wave height *H* at a local point and is given by Eq. (15), referring to the relationship given by Uda and Kawano [3]:

In the numerical simulation of beach changes, the sand transport equation and the continuity equation are solved on the *x*-*y* plane by the explicit finite-difference method employing a staggered mesh scheme.

The wave field is calculated using the energy balance equation given by Mase [9], in which the directional spectrum *D* (*f*, *θ*) of the irregular waves is the variable to be solved, with the energy dissipation term due to wave breaking in [10]. Here, *f* and *θ* are the frequency and wave direction, respectively. In this method, wave refraction, wave breaking, and wave diffraction in the wave-shelter zone can be calculated with a small calculation load. The energy dissipation term due to wave breaking *Φ*, which is incorporated into the energy balance equation of Eq. (16), is given by Eq. (17):

Here, *D* is the directional spectrum, (*Vx*, *Vy*, *Vθ*) is the energy transport velocity in the (*x*, *y*, *θ*) space, *F* is the wave diffraction term in [9], *K* is the coefficient of the wave-breaking intensity, *h* is the water depth, *Cg* is the wave group velocity (*Γ* is the ratio of the critical breaker height to the water depth on the horizontal bed, and *γ* is the ratio of the wave height to the water depth. As the spectrum of the incident waves, a combination of the Bretschneider-Mitsuyasu-type frequency spectrum and the Mitsuyasu-type directional function is used as in Goda [11]. Figure 3 shows the flowchart for numerical simulation using the Type 3.

To prevent the location where the berm develops from being excessively seaward compared with that observed in the experiment or the field, a lower limit is considered for *h* in Eq. (17). As a result of this procedure, wave decay near the berm top is reduced, resulting in a higher landward sand transport rate. In the calculation of the wave field on land, the imaginary depth *h’* between the minimum depth *h*_{0} and *h*_{R} is considered as follows, similarly to the 3D model of Shimizu et al. [12]:

In addition, at locations whose elevation is higher than *h*_{R}, the wave energy is set to 0. Similarly to the Type 1, beach changes are obtained by solving the continuity equation. In the calculation of the formation of a bay barrier in [13], the sand transport equation without Ozasa and Brampton’s term in Eq. (12) is employed, i.e., *K*_{2} = 0, and in the calculation of the *P* value in Eq. (12), Eq. (13) multiplied by the coefficient *F* expressed by Eq. (20) is employed:

Here, tan*βw* is the seabed slope measured along the direction of wave propagation. In Eq. (19), *F* is introduced so that the wave dissipation rate can respond to the seabed slope when the wave dissipation in the surf zone is considered.

Eq. (12) shows that the sand transport flux can be expressed as the sum of the component along the wave direction and the components due to the effect of gravity normal to the contours and the effect of longshore currents parallel to the contours. To investigate the physical meaning of Eq. (12), neglecting the additional term given by Ozasa and Brampton [1], that is, assuming *K*_{2} = 0 in Eq. (12), Eqs. (22) and (23) are derived for the cross-shore and longshore components of sand transport, *qn* and *qs*, respectively. Furthermore, under the condition that the seabed slope is equal to the equilibrium slope, Eq. (23) reduces to Eq. (24):

In Eq. (22), the cross-shore sand transport *qn* becomes 0 when the local seabed slope is equal to the equilibrium slope, and the longshore sand transport *qs* becomes 0 when the wave direction coincides with the normal to the contour lines, as shown in Eqs. (23) and (24). When a discrepancy from these conditions arises, sand transport is generated by the same stabilization mechanism as that given by the sand transport equation in the Type 1.

Comparing Eqs. (22) and (24) in the Type 3 with Eqs. (31) and (34) in Chapter 1, which are the fundamental equations for the cross-shore and longshore components of sand transport in the Type 1, the sand transport equations in the Type 1 have a single coefficient, whereas the equations in the Type 3 are expressed using two coefficients, i.e., the intensities of cross-shore and longshore sand transports can be independently determined. The equilibrium slope *qn* = 0 under the condition of oblique wave incidence, satisfies the relation *α*. On the other hand, in the Type 3, the relation *qn* = 0 in Eq. (22). Thus, *α*. This is a result of multiplying

Taking the above into account, the first term in the parentheses in Eq. (12) gives the sand transport in the case that the rates of longshore and cross-shore sand transports are equal (*Ks* = *Kn*), and the second term is the additional longshore sand transport in the case that the rates are different (*Ks* > *Kn*). The physical meaning of the second term is that longshore sand transport is generated by the small angular shift that occurs when the wave direction is incompletely reversed in the oscillatory movement due to waves, and the second term also models the additional longshore sand transport due to the effect of longshore currents, which is only partially included in the first term.

### 2.4 Type 4 BG model

This model is used for the calculation of the development of sand spits and cuspate forelands with rhythmic shapes and the formation of a cuspate foreland in Chapter 7. Although the intensity of sand transport *P* in the Type 3 was assumed to be proportional to the third power of the amplitude of the bottom oscillatory velocity, *um*, due to waves, as shown in Eq. (13), *P* was given by the wave energy dissipation rate *Φ*_{b} due to wave breaking at a local point in the Type 4, which can be determined by the calculation of the plane wave field:

To calculate the wave field, the numerical simulation method using the energy balance equation given by Mase [9], in which the directional spectrum of irregular waves is the variable to be solved, is employed, similarly to the Type 3, with an additional term of energy dissipation due to wave breaking [10], similarly to Eq. (16). *Φ*_{b} in Eq. (25) is calculated using Eq. (26), which defines the sum of the energy dissipation of each component wave due to breaking:

Here, *fD* 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 the wave height to the water depth *H/h*. In addition, similarly to the calculation in the Type 3, a lower limit is set for the water depth *h* in Eq. (18) in the calculation of the plane wave field. To calculate the wave field in the wave run-up zone, an imaginary depth is assumed, similarly to the Type 3. Furthermore, the wave energy at locations with elevations higher than *h*_{R} is set to 0. The flowchart for numerical simulation using the Type 4 is the same as that in the Type 3, as shown in Figure 3.

When investigating the physical meaning of Eq. (12) in the Type 3, the effect of the additional term given by Ozasa and Brampton [1] was neglected. Here, the effect of the additional term is included, and Eq. (27) is derived for the longshore component of sand transport, *qs*, along with Eq. (22) for the cross-shore component of sand transport, *qn*. In addition, under the condition that the seabed slope is equal to the equilibrium slope, Eq. (27) reduces to Eq. (28):

The cross-shore sand transport *qn* is the same as Eq. (22) in the Type 3, whereas the longshore sand transport *qs* given by Eqs. (27) and (28) has the additional term of Ozasa and Brampton [1]. When the total sand transport *Qs* is calculated by integrating Eq. (28) in the cross-shore direction, it coincides with the CERC-type formula [15] with the additional term [1] as in Eq. (29) [4]. Here, the relation *α* is approximately given by the breaker angle *αb*, and the longshore gradient of the breaker height is used in place of the longshore gradient of the wave height in the Ozasa and Brampton’s term together with the assumption that the integral of *Φ* in the cross-shore direction is equal to the energy flux per unit length of the coastline at the breaking point:

### 2.5 Type 5 BG model

To predict the topographic changes of an island or a cuspate foreland in a shallow water body under the action of waves randomly incident from every direction, the wave power *P* is calculated using a coordinate system different from that for the calculation of beach changes, assuming that waves propagate in a straight line by neglecting the effect of the wave refraction. Since the wave field itself significantly changes with time in response to the deformation of the topography, the calculation of the wave power *P* and beach changes is carried out recurrently every time step. This model is employed in the numerical simulation of the elongation and merging of bay mouth sand spits in Chapter 5 and the interaction of islands in Chapter 6.

The sand transport equation employed is Eq. (30), which uses the expression of the wave energy evaluated at the breaking point. The variables in Eq. (30) are given by Eqs. (31)–(36) along with Eqs. (5), (7), (8), and (9) in the Type 1 BG model:

Here, *Ks* is the longshore and cross-shore sand transport coefficients, and the index *i* in Eqs. (35) and (36) is the mesh number along the *xw*-axis, and min denotes the selection of the smaller of either value in the parentheses. In the calculation, the local beach slope measured along the wave ray is used as the beach slope in Eq. (31), as shown in Eq. (32).

For the calculation of the *P* value, another coordinate system different from that for the calculation of beach changes is used as in [16], in which the *xw*- and *yw*-axes are taken along the wave direction and the direction normal to the wave direction, respectively (Figure 4). The fixed coordinate system (*x*, *y*) is used for the calculation of beach changes with a rectangular calculation domain ABCD, whereas the *P* value is calculated in the rectangular domain A′B′C′D′ including the domain ABCD with the coordinate system (*xw*, *yw*). In the calculation of the wave field, wave refraction is neglected, and waves are assumed to propagate in a straight line while maintaining the wave incident angle. The distance along the *xw*-axis is subdivided by a mesh of Δ*xw*, and a cumulative function of *ε*(*Z*), *ε*(*Z*) is assumed to have a uniform distribution (Eq. (7)), and *h*_{c}, decreases with the water depth, and becomes 0 in the zone higher than *h*_{R}.

Using Eq. (33), Eq. (31) is transformed into Eq. (34). Thus, the *P* value can be determined from the derivative of *xw*-axis using Eq. (34). *xw*-axis from the starting point of wave incidence in the direction of wave propagation; *i +* 1)^{th} point can be calculated from Eq. (35) with the given value of *i*^{th} point when the initial value of

Furthermore, the *P* value at the (*i + 1/2*)^{th} point is calculated from Eq. (36), which is derived from Eq. (34). Note that in calculating *i*+1)^{th} point, the smaller value calculated from Eq. (35) given the elevation *Z*(*i*+1) at the (*i*+1)^{th} point and *i*^{th} point is adopted. As a result, the value of *Z*) between the offshore end and a designated point along the *xw*-axis can be adopted.

Using this procedure, the depth distribution along the *xw*-axis between the offshore end and a certain point is automatically taken into account in the calculation of the *P* value. For example, when there is a location with an elevation higher than *h*_{R}, the *P* value in the shoreward zone is automatically reset as 0 regardless of the elevation and water depth. This procedure becomes important when the beach profile along the *xw*-axis has several uneven shapes.

Consider the case where the *xw*-axis passes through their tips of two sand spits. The wave energy is reduced when waves pass through near the tip of the first sand spit, resulting in the reduction of the wave energy reaching the second sand spit. Using this procedure, the wave-sheltering effect due to the existence of multiple sand spits can be automatically evaluated, in contrast to the method where the *P* value is calculated by substituting the local elevation *Z* into *ε*(*Z*) in Eq. (31). Moreover, the *P* value on an impermeable breakwater is set to 0 to take the wave-sheltering effect of the breakwater into account. The *P* value integrated from a location on land whose elevation exceeds *h*_{R} to an offshore end along the *xw*-axis is always equivalent to (*ECg*)_{b} regardless of the seabed topography. Because (*ECg*)_{b} corresponds to the entire power of the incident waves, the exact satisfaction of this condition is reasonable.

The calculation of the *P* value is independently carried out along each the *xw*-axis. The *P* value calculated at point (*xw*, *yw*) is memorized, and the *P* value at point (*x*, *y*) necessary for the calculation of beach changes is interpolated from the memorized value at point (*xw*, *yw*). The mesh intervals Δ*xw* and Δ*yw* are the same as Δ*x* and Δ*y*. Here, Δ*x* and Δ*y* are the mesh intervals of the coordinate system for the calculation of beach changes, and we assume the equivalent condition of Δ*x =* Δ*y*. The beach changes are calculated by explicitly solving the continuity equations on the staggered meshes using the sand transport fluxes obtained from Eq. (30). Figure 5 shows the flowchart for numerical simulation using the Type 5 BG model.

The incident wave direction at each time step in the calculation of beach changes is selected to be a value determined randomly so as to satisfy the probability distribution function of the incident wave direction, although the incident wave height is assumed to be constant. For example, in Chapter 5, the energy distribution for multidirectional irregular waves with a directional spreading parameter of *S*_{max} = 10 is used, and in Chapter 6, the wave incidence with the direction ranging between 0 and 360° with the same probability is employed as the probability distribution function of the wave direction. The *P* value, which is subject to change with the propagation of waves, is recalculated at every time step in the calculation of the beach changes.

### 2.6 Type 6 BG model

The Type 6 is employed in the numerical simulation of the segmentation of a closed water body, as mentioned in Chapter 8. The height of wind waves in a closed body is predicted using Wilson’s formula on the basis of the fetch distance and wind velocity, and then sand transport fluxes are calculated. To calculate the topographic changes in a lake under the action of wind waves randomly incident from all directions, the wave power *P* is calculated using a coordinate system different from that for the calculation of beach changes. Since the wave field itself significantly changes with time in response to the deformation of the topography, the calculation of the height of the wind waves and beach changes was carried out recurrently every time step.

As the sand transport equation, Eq. (37) is used, which is expressed using the wave energy at the breaking point, similarly to the Type 5. The variables in Eq. (37) are given by Eqs. (38)–(41) along with the relations of Eqs. (6)–(8) in the Type 1 and Eqs. (31) and (32) in the Type 5. Although the *P* value in Eq. (31) was calculated from Eqs. (33)–(36) in the Type 5, the *P* value was evaluated directly from Eq. (31) in the Type 6:

Here, *F* is the local fetch distance and *U* is the wind velocity. The index *i* in Eq. (38) is the mesh number along the *xw*-axis. Prior to the calculation of beach changes, the significant wave height at a point is calculated using Wilson’s equation (Eq. (41)), described in Wilson [17] and Goda [18], given the local fetch *F* at a point and wind velocity *U*.

In this calculation, a fixed coordinate system (*x*, *y*) is adopted for the calculation of beach changes with the rectangular calculation domain ABCD, as shown in Figure 6, whereas another coordinate system (*xw*, *yw*) is set corresponding to the wave direction, and the wave height is calculated in the rectangular domain A′B′C′D′ including the domain ABCD.

Neglecting the wave refraction effect, waves are assumed to propagate in the same direction as the wind. The *xw*-axis is subdivided by mesh intervals of Δ*xw*. The fetch *F* is added from upwind to downwind along the *xw*-axis using Eq. (38). When a grid point is located on land and the downslope condition of *dZ/dxw* ≤ 0 is satisfied, the local fetch is reset to *F* = 0 (Eq. (40)). Also we reset *F* to 0 on structures. When the grid point is located in the lake, *F* is recalculated. By this procedure, the wave-sheltering effect by a sand spit or a structure on wind-wave development is taken into account. Then, the significant wave height is calculated with Wilson’s [17] equation (Eq. (41)) using the wind fetch *F* and wind velocity *U* [18], and this wave height is assumed to be equal to the breaker height. Simultaneously, the wave power *P* was calculated and assigned to each grid point in the coordinate system (*xw*, *yw*). The calculation of the *P* value is independently carried out along each *xw*-axis, similarly to the Type 5. The beach changes are calculated by explicitly solving the continuity equations on the staggered meshes using the sand transport fluxes obtained from Eq. (37). The flowchart for numerical simulation using the Type 6 is shown in Figure 7.

The wind direction at each time step in the calculation of beach changes is selected to be a value determined randomly so as to satisfy the probability distribution function of the wind direction, although the wind velocity is assumed to be constant. For example, in Chapter 8, wind is assumed to blow uniformly from all directions between 0 and 360°, that is, a symmetric circular distribution, together with the calculation with asymmetric probability distribution of occurrence of wind direction. In every time step of the calculation of beach changes, the wind direction is reset randomly, and the distribution of the *P* value is recalculated.

### 2.7 Type 7 BG model

Tung et al. [19] predicted the evolution of ebb-tidal deltas using the Delft 3D model. In their model, the full equations of waves and nearshore and tidal currents were solved to predict three-dimensional bathymetric changes. Beck and Kraus [20] have also carried out the numerical simulation of the development of an ebb-tidal delta by solving the equations of waves and nearshore and ebb-tidal currents. For practical applications, however, the development of a model by which bathymetric changes can be more easily predicted and used to investigate measures against beach erosion is also desirable. Therefore, the model described below was developed [21, 22].

A model for predicting the formation of ebb-tidal deltas under the combined effect of waves and ebb-tidal currents was developed on the basis of the BG model [23] with an analogy of the velocity distribution of ebb-tidal currents to the wave diffraction coefficient, which was calculated using the angular spreading method for irregular waves [5, 6].

Assume that the total sand transport

For the equation of sediment transport due to waves, the equation in the Type 1 is improved to take the differences in the intensity of cross-shore and longshore sand transports into account. Consider the Cartesian coordinates (*x*, *y*), where the *x*- and *y*-coordinates are taken to be the cross-shore distance (positive for shoreward) and longshore distance parallel to the shoreline, respectively. Assume that waves are obliquely incident to a coast with a slope of tan*β*. Then, the net sand transport flux due to waves

Here, *Gwx* and *Gwy* are given by

Here, *qx* is the cross-shore component of the sediment transport flux (positive for shoreward), and *qy* is the longshore component of sediment transport flux. *Kx* and *Ky* are the coefficients of cross-shore and longshore sand transports, respectively. tan*β*_{c} is the equilibrium slope for which zero net cross-shore transport occurs when waves are incident from the direction normal to the slope.

Regarding Eq. (43), Serizawa et al. [21] used the sediment transport equation in the Type 1. This implicitly assumes that the coefficients of cross-shore and longshore sand transports are equivalent. Here, to take the difference in the intensities of cross-shore and longshore sand transport into account, the coefficients of cross-shore and longshore sand transports, *Kx* and *Ky*, are determined independently. This equation is the same as the sand transport equation in the Type 2 BG model without the additional term given by Ozasa and Brampton [1].

Regarding the sand transport flux

Here, the subscript *R* denotes the ebb-tidal currents, *θR* is the angle between the direction of the ebb-tidal currents and the *x*-axis, and tan*ϕ* is the angle of repose of the sand. In the original equation of Bailard and Inman [8], the coefficient *GR* in Eq. (45) is expressed in terms of the instantaneous velocity, the angle of repose of the sand, and a friction factor, but here sand transport due to currents is assumed to satisfy an equation of the same form as Eq. (45), and all these effects are included in the coefficient *GR* in Eq. (45). The first and second terms in the parentheses in Eq. (45) correspond to the action of currents and the downslope action due to gravity, respectively.

In fact, the sand transport due to ebb-tidal currents is comprised of bedload and suspended load, and Eq. (45), therefore, should be expressed by the ratio of ebb-tidal currents to the falling velocity as in [24, 25]. However, in this study, it was assumed that Eq. (45) could be employed as the total load formula, and *GR* is given by

Here, *KR* is the coefficient of sand transport due to ebb-tidal currents, and *KV* is the ratio of the ebb-tidal current velocity to *V*_{0}, *V/V*_{0}, where *V* is the velocity of the ebb-tidal currents and *V*_{0} is a reference velocity at the inlet. *Fw* is a characteristic value of the intensity of wave action and is given by the wave energy flux at the breaking point of the reference point, (*ECg*)_{b0}, divided by the sum of *h*_{c} and *h*_{R}, as in Eq. (47). *V*_{1} is the flow velocity of the ebb-tidal currents on the plane bottom. *K*_{V1} is the ratio of the velocity to the flow velocity at the inlet on the plane bottom. In addition, *h*_{0} and *h* in Eq. (48) are the reference depth on the plane bottom and water depth, respectively. *h*_{c2} and *h*_{R2} are the lower and upper limit depths of bathymetric changes due to ebb-tidal currents, respectively.

Regarding the reason why *GR* is given as Eq. (46), *GR* in Eq. (45) represents the work done by the ebb-tidal currents in sand movement and is proportional to the energy dissipation rate of the ebb-tidal currents [7, 24]. Here, the action of the ebb-tidal currents is evaluated in macroscopically, and Eq. (45) is assumed to have the same form as Eq. (43); *GR* is expressed by the intensity of the ebb-tidal currents relative to the wave intensity (*ECg*)_{b0} and is assumed to be zero outside the zone between the upper and lower limit depths of bathymetric changes due to the ebb-tidal currents. Furthermore, *GR* is assumed to be proportional to the third power of the velocity ratio so that the energy dissipation rate is proportional to the third power of the velocity [7, 24].

When the distribution of *K*_{V1}, as expressed by Eq. (49), on a flat bottom is given, the effect of a change in current velocity associated with a depth change is evaluated using Eq. (48) so that the mass conservation of the fluid is satisfied. *K*_{V1} is calculated by applying the angular spreading method for irregular waves [5, 6], taking into consideration the fact that the jet-like velocity distribution of ebb-tidal currents at an inlet is very similar to the distribution of the wave diffraction coefficient at the opening of offshore breakwaters. Finally, the coefficient of sand transport due to ebb-tidal currents *KR* in Eq. (46) becomes a coefficient including the ratio of the action due to the currents to the wave action. In the calculation, the mean beach slope before the construction of artificial structures is assumed to be the equilibrium slope, considering the long-term beach changes. Although the actions of waves and ebb-tidal currents both simultaneously occur in a tidal inlet, it is assumed that there is no sand supply during the ebb-tidal currents. Figure 8 shows the flowchart for numerical simulation using the Type 7.

### 2.8 Type 8 BG model

When an artificial reef (submerged breakwater) is constructed on a coast with detached breakwaters, strong shoreward currents are generated by the forced wave breaking, producing rip currents at the opening between the artificial reef and the existing detached breakwaters. The effects of these currents are not considered in the Type 3, which employs only the wave field without calculating the nearshore currents, although the effect of longshore currents is implicitly considered in the sand transport equation described using the local wave characteristics. Here, the effect of the nearshore currents induced by forced wave breaking is incorporated into the model by calculating the nearshore currents. Zanuttigh [26] and Kuroiwa et al. [27] predicted the topographic changes around the low-crested structures by numerical simulation using the so-called 3D model. In this study, an improved BG model is proposed to predict the topographic changes around the artificial reefs. In this study, the sand transport equation based on the concept of the equilibrium slope is employed.

The Type 3 BG model is improved by taking both the wave field and the current velocity at a local point into account. The fundamental equations of the model are given by Eqs. (50)–(54) in the Cartesian coordinates (*x*, *y*), assuming that the sand transport flux *qx*, *qy*) is a linear sum of the component due to waves, *qwx*, *qwy*), and that due to currents, *qcx*, *qcy*):

Here, *um* is the amplitude of the seabed velocity due to the orbital motion of waves defined by Eq. (14) in the Type 3 BG model. *α*_{1} is the angle between the wave direction and the direction normal to the contour lines, *V* is the velocity of the currents, *Kw* and *Kc* are the coefficients of sand transport due to waves and currents, respectively, *α*_{2} is the angle between the direction of the currents and the direction normal to the contour lines, and tan*ϕ* is the slope of the angle of repose of sand. Here, *b*_{2} is assumed to be 0 in the calculation.

The sand transport equations of Eqs. (51) and (53) were derived on the basis of Bagnold’s concept as the sand transport equation in the Type 7 BG model. Eq. (51), expressing the effect due to waves, is the same as the equation when assuming that the coefficient of cross-shore sand transport is equivalent to that of longshore sand transport, and the Ozasa and Brampton’s term [1] is excluded from the sand transport equation in the Type 3.

In Eq. (53), the wave power by which sand transport is caused by nearshore current is assumed to have the form of *u*_{m}^{2}*V*, referring the previous studies on 3D beach change models in the nearshore zone [4, 12, 28, 29, 30, 31]. In the calculation of beach changes, the wave field and nearshore currents are calculated first, and then beach changes are predicted.

Similarly to the calculation in the Type 3 BG model, the wave field is calculated using the energy balance equation [9] with the energy dissipation term due to wave breaking [10]. As the directional spectrum of the incident waves, the combination of the Bretschneider-Mitsuyasu-type frequency spectrum and the Mitsuyasu-type directional function is used [11]. To calculate the nearshore currents, the upwind finite-difference scheme is used for the 2-D momentum equation [4], along with the lateral diffusion term of Larson and Kraus [32]. Because recurrent feedback calculations are necessary in response to the topographic changes in this calculation, the energy balance equation method is used to reduce the calculation load.

In the calculation of the wave field on land, the imaginary depth is assumed between *h*_{R} and the shoreline. In this case, the imaginary depth *h’* can be obtained from Eq. (18) using the minimum depth *h*_{0} and *h*_{R}. In particular, in the calculation in Chapter 4, *h*_{0} is assumed to be 1 m. The wave energy is set to be 0 in the area with an elevation higher than the berm height. Although *h*_{c} is assumed to be 2.5*H*, where *H* is the wave height at a local point, the increment due to the effect of the strong rip currents is evaluated as

Here, *V* is the velocity of the nearshore currents, and in the calculation in Chapter 4, a coefficient of 0.5 is employed. Similar to the calculation in the Type 1 BG model, beach changes are obtained by solving the continuity equation. Figure 9 shows the flowchart for the Type 8.

## 3. Discretization method

### 3.1 Discretization of mass conservation equation

Beach changes can be calculated using the mass conservation equation (Eq. (28) in Chapter 1). In the calculation, the coastal domain is discretized using 2D elements with widths Δ*x* and Δ*y*, as schematically shown in Figure 10. The calculation points of the seabed elevation *Z* and sand transport rates are set in staggered meshes with a difference of 1/2 mesh, and the equations are solved by the explicit finite-difference method.

First, the calculation point of the elevation *Z* at *i* ≤ *Nx*, 1 ≤ *j* ≤ *Ny*) is evaluated at each point. Here, the subscripts *i* and *j* are the cell numbers taken in the *x*- and *y*-directions, and *Nx* and *Ny* are the numbers of cells in the *x*- and *y*-directions, respectively. Then, *i* ≤ *Nx* + 1, 1 ≤ *j* ≤ *Ny*) and *i* ≤ *Nx*, 1 ≤ *j* ≤ *Ny* + 1) are set at points *Z*_{(i, j)} by 1/2 mesh in the *x*- and *y*-directions, respectively.

When the sand transport flux *t* can be calculated using Eq. (56), which is the discretized form of Eq. (29) in Chapter 1:

In the calculation of Eq. (56) in a closed system, sand transport is set to 0 along the outer boundaries of the domain as the boundary conditions. A free boundary condition of

As a boundary condition of the structure, sand transport is set to 0 along the boundary. In the calculation of the Type 1 BG model, given the initial seabed topography, the distribution of *H*_{b}, *α*_{b}, *h*_{c}, and *h*_{R}, and the equilibrium slope, sand transport fluxes (*qx*, *qy*) are calculated using the equations given later, and the change in seabed elevation after time Δ*t* is calculated using Eq. (56). These procedures are repeated recurrently.

### 3.2 Discretization of sand transport equation

All the BG models except the Types 2 and 7 are expressed as

where

Similarly, the coefficients of the sand transport equation of the Type 3 in Eq. (60) are given by Eq. (61):

Since the sand transport equations of the various types can be written as the form in Eq. (57), the sand transport flux of Eq. (57) can be written as *x*- and *y*-component expressions for use in the calculation procedure as follows:

Using Eq. (62), the four components of the sand transport equation in Eq. (56) are expressed as follows:

where

where

Finally,

### 3.3 Procedure to estimate sand transport near berm top and depth of closure

In the estimation of the intensity of sand transport near the berm top and at the depth of closure, the intensity of sand transport is linearly reduced to 0 by multiplying by the reduction ratio near *h*_{R} or *h*_{c} to prevent sand from being deposited in the zone higher than *h*_{R} and the beach from being eroded in the zone deeper than *h*_{c} [33]. In this method, the evaluation method for sand transport on an exposed rock bed in the 3-D beach change model proposed by Ikeno et al. [34] was employed as in [35]. The sand transport fluxes *μ* as in Eq. (71), where *μ* is given by Eq. (72). Δ*Z*_{0} is the thickness of the sand layer from which the reduction in sand transport begins, and Δ*Z* is given by Eq. (73) near *h*_{R} and by Eq. (74) near *h*_{c} (Figure 12). On a fixed bed comprising a coral reef, Eq. (74) is employed near the depth of the solid bed.

The calculation points of the sand transport and seabed elevation have a 1/2 mesh difference because of the staggered meshes with a 1/2 mesh interval. Here, on the basis of the direction of sand transport, the *Z* value at a point immediately downcoast and upcoast of the calculation point of sand transport is used in Eqs. (73) and (74), respectively. For example, in the calculation of the *x*-component in Eq. (74), the *Z* value at a point immediately upcoast of the calculation point of sand transport is used and is calculated using Eq. (75) as in Figure 13:

Here, the subscript *iup* represents the location of the *Z* value at a point immediately updrift of *qx* (Figure 13). A similar method can be used for the *y*-component of *q*. In the calculation of the *x*-component in Eq. (73), the *Z* value at a point immediately downcoast of the calculation point of sand transport is used and calculated using Eq. (76) (Figure 14).

Here, the subscript *idown* represents the location of the *Z* value at a point immediately downdrift of *qx* (Figure 14). A similar method can be used for the *y*-component of *q*. Here, the thickness of the sand layer Δ*Z*_{0}, from which the reduction in sand transport begins, in Eq. (72) was empirically determined. In the calculation in Chapter 3 using the Type 1 BG model, Δ*Z*_{0} is determined using Eq. (77), considering that the reduction of sand transport begins when the slope angle coincides with the angle of the slope in a right-angled triangle whose base has a length of one mesh and whose height is equal to the depth difference, similarly to [2, 6]:

In the other calculation, Δ*Z*_{0} is selected as a value within the following range:

### 3.4 Procedure for calculating cross-shore sand movement due to the effect of gravity

Regarding the seaward sand movement in regions above *h*_{R} and deeper than *h*_{c}, the falling of sand to the deep sea bottom when the local slope exceeds the slope of the angle of repose of sand was evaluated by calculating the sand transport flux due to the gravity [35]. In this calculation, the sand transport flux

Here, *n* is the local coordinate taken along the direction normal (shoreward) to the contour lines, *Ag*. Thus, Eq. (79) expresses the effect of the seabed slope, that is, the steeper the seabed slope, the larger the sand transport. Moreover, when the local slope is gentler than the slope of the angle of repose,

Using Eq. (80), Eq. (79) can be written as Eq. (81):

Eq. (81) can be expressed by the basic expression for sand transport of Eq. (57), as mentioned earlier, and the coefficients *A*_{1}, *A*_{2}, and *A*_{3} of each term of Eq. (57) are given by

Accordingly, *A*_{g} is known. Since sand transport due to the effect of gravity is much faster than that due to waves, a large value can be selected for *Ag*. However, when an excessively large value is employed for *Ag*, it is difficult to stably carry out the numerical calculation. Here, *Ag* is designated as the maximum value in a range, in which the numerical simulation can stably proceed as follows.

First, *qn* becomes Eq. (84). Assuming the one-dimensional problem in the *n*-direction, Eq. (84) is substituted into the continuity equation of Eq. (85), and we obtain the one-dimensional diffusion equation of Eq. (86) in the *n*-direction with a diffusion coefficient of *Ag*. To solve Eq. (86) using the explicit finite-difference method, Eq. (87) must be satisfied as the stability condition of Eq. (86):

Here, Δ*n* is the mesh interval along the *n-*axis and Δ*t* is the time interval. Although Eq. (87) is normally used to determine the upper limit of Δ*t* given *Ag* and Δ*n*, here it is employed as a relationship to determine the upper limit of *Ag* given Δ*n* and Δ*t*. Eq. (87) can be rewritten as Eq. (88) to provide a stability condition for *Ag*. Thus, the right term of Eq. (88) gives the upper limit of *Ag*:

The coefficient *A*_{g} that ensures the stability of Eq. (89) is given:

Finally,