Eight Types of BG Models and Discretization

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.

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, um, 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

Table 1.

Eight types of BG models and their applications.

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 $tan\overline{\beta }$. Φ is given by the wave energy flux transported per unit width of the shoreline, ${\left({\mathit{EC}}_g\right)}_{\mathrm{b}}cos{\alpha }_{\mathrm{b}}$, divided by the cross-shore distance $R=\left(h_{\mathrm{c}}+h_{\mathrm{R}}\right)/tan\overline{\beta }$ between the depth of closure, h_c, and the berm height, h_R, corresponding to the wave run-up height (Figure 1). Then, assuming $tan\overline{\beta }\approx tan{\beta }_{\mathrm{c}}cos{\alpha }_{\mathrm{b}}$ (Eq. (32) in Chapter 1), Φ becomes ${\left({\mathit{EC}}_g\right)}_{\mathrm{b}}\mathrm{co}{\mathrm{s}}^2{\alpha }_{\mathrm{b}}/\left(h_{\mathrm{c}}+h_{\mathrm{R}}\right)/tan{\beta }_{\mathrm{c}}$. Furthermore, introducing the depth distribution of the sand transport intensity ε(Z) and transforming $1/\left(h_{\mathrm{c}}+h_{\mathrm{R}}\right)$ into ε(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

$$\overrightarrow{q}=\frac{G}{tan{\beta }_{\mathrm{c}}}\left[tan{\beta }_{\mathrm{c}}\overrightarrow{e_w}-\overrightarrow{\nabla Z}\right].$$E1

Figure 1.

Zone with cross-shore sand movement and dissipation of wave energy.

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

$$q_x=\frac{G}{tan{\beta }_{\mathrm{c}}}\left[tan{\beta }_{\mathrm{c}}cos{\theta }_w-\partial Z/\partial x\right]$$E2a

$$q_y=\frac{G}{tan{\beta }_{\mathrm{c}}}\left[tan{\beta }_{\mathrm{c}}sin{\theta }_w-\partial Z/\partial y\right]$$E2b

$$G=C_0{K}_{1}P$$E3

$$P=\mathit{\Phi }=\epsilon \left(Z\right){\left({\mathit{EC}}_g\right)}_{\mathrm{b}}\mathrm{co}{\mathrm{s}}^2{\alpha }_{\mathrm{b}}tan{\beta }_{\mathrm{c}}$$E4

$${\int }_{-h_{\mathrm{c}}}^{h_{\mathrm{R}}}\epsilon \left(Z\right)\mathit{dZ}=1$$E5

$$\epsilon \left(Z\right)=\left\{\begin{array}{l}\frac{2}{h_{\mathrm{c}}^3}\left(\frac{h_{\mathrm{c}}}{2}-Z\right){\left(Z+h_{\mathrm{c}}\right)}^2\left(-h_{\mathrm{c}}\le Z\le h_{\mathrm{R}}\right)\\ 0\left(Z<-h_{\mathrm{c}},h_{\mathrm{R}}<Z\right)\end{array}\right.$$E6

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:

$$\epsilon \left(Z\right)=\left\{\begin{array}{l}1/\left(h_{\mathrm{c}}+h_{\mathrm{R}}\right)\left(-h_{\mathrm{c}}\le Z\le h_{\mathrm{R}}\right)\\ 0\left(Z<-h_{\mathrm{c}},h_{\mathrm{R}}<Z\right)\end{array}\right.$$E7

The wave energy flux at the breaking point (EC_g)_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:

$${\left({\mathit{EC}}_g\right)}_{\mathrm{b}}=C_1{\left(H_{\mathrm{b}}\right)}^{\frac{5}{2}}$$E8a

$$C_1=\frac{\mathit{\rho g}}{k_1}\sqrt{g/\gamma }\left(\gamma \approx 0.8\right)$$E8b

$$k_1=\left\{\begin{array}{l}8\text{for regular waves}\\ {\left(4.004\right)}^2\text{for irregular waves}\end{array}\right.$$E8c

Here, γ is the ratio of the breaker height to the water depth. k₁ = (4.004)² 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:

$$\begin{array}{l}{\alpha }_{\mathrm{b}}\approx \alpha ={\theta }_w-{\theta }_n\\ ={\theta }_w-\mathrm{ta}{\mathrm{n}}^{-1}\left(\frac{\mathit{\partial Z}}{\mathit{\partial y}}/\frac{\mathit{\partial Z}}{\mathit{\partial x}}\right)\end{array}$$E9

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 $\overrightarrow{q}=(q_x,q_y)$ are set in staggered meshes with a difference of 1/2 mesh. The explicit finite-difference method is used in the calculation. Figure 2 shows the flowchart for numerical simulation using the Type 1.

Figure 2.

Flowchart for numerical calculation using Typ. 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 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 $\left|cos{\alpha }_{\mathrm{b}}\right|$ in the calculation of 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²(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₂/K₁, 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):

$$\begin{array}{l}{q}_x=\frac{G_x}{tan{\beta }_{\mathrm{c}}}\left[tan{\beta }_{\mathrm{c}}cos{\theta }_w-\partial Z/\partial x\right]\\ \left(-h_{\mathrm{c}}\le Z\le h_{\mathrm{R}}\right)\end{array}$$E10a

$$\begin{array}{l}{q}_y=\frac{G_y}{tan{\beta }_{\mathrm{c}}}\left[tan{\beta }_{\mathrm{c}}\left(sin{\theta }_w-\frac{1}{tan\overline{\beta }}\frac{K_2}{K_y}\frac{\mathit{\partial }{H}_{\mathrm{b}}}{\mathit{\partial y}}\right)-\partial Z/\partial y\right]\\ \left(-h_{\mathrm{c}}\le Z\le h_{\mathrm{R}}\right)\end{array}$$E10b

$$G_x=C_0{K}_{x}P$$E11a

$$G_y=C_0{K}_{y}P$$E11b

Here, θ_w is the wave angle measured counterclockwise from the x-axis. K_x and K_y are the coefficients of cross-shore and longshore sand transport, respectively, K₂ is the coefficient of the term given by Ozasa and Brampton [1], and $tan\overline{\beta }$ is the seabed slope at the breaker point. Here, we assume $tan\overline{\beta }=tan{\beta }_{\mathrm{c}}$. ε(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

$$\begin{array}{l}\overrightarrow{q}=C_0\frac{P}{tan{\beta }_{\mathrm{c}}}\left\{\begin{array}{l}{K}_n\left(tan{\beta }_{\mathrm{c}}{\displaystyle \overrightarrow{e_w}}-\left|cos\alpha \right|{\displaystyle \overrightarrow{\nabla Z}}\right)\\ +\left\{\left(K_s-K_n\right)sin\alpha -\frac{K_2}{tan\overline{\beta }}\frac{\mathit{\partial H}}{\mathit{\partial s}}\right\}tan\beta {\displaystyle \overrightarrow{e_s}}\end{array}\right\}\\ \left(-h_{\mathrm{c}}\le Z\le h_{\mathrm{R}}\right)\end{array}$$E12

$$P=\rho u_m^3$$E13

$$u_m=\frac{H}{2}\sqrt{\frac{g}{h}}.$$E14

Here, $tan\beta \overrightarrow{e_s}=\left(-\mathit{\partial Z}/\mathit{\partial y},\mathit{\partial Z}/\mathit{\partial x}\right)$, K_s and K_n are the coefficients of longshore sand transport and cross-shore sand transport, respectively, $\mathit{\partial H}/\mathit{\partial s}=\overrightarrow{e_s}\cdot \overrightarrow{\nabla H}$ is the longshore gradient of the wave height H measured parallel to the contour lines, and ${\displaystyle \overrightarrow{\nabla H}}=\left(\mathit{\partial H}/\mathit{\partial x},\mathit{\partial H}/\mathit{\partial y}\right)$. $tan\overline{\beta }$ is the characteristic slope of the breaker zone, and h is the water depth. Furthermore, $tan\overline{\beta }=tan{\beta }_{\mathrm{c}}$ is assumed. u_m 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 = τu_t = ρC_fu_t³ for the instantaneous wave energy dissipation rate Φ_t to derive their sand transport equation, where τ is the bottom shear stress, u_t is the instantaneous velocity, and C_f 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 u_m 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, K_s and K_n in Eq. (12), respectively. u_m 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]:

$$h_{\mathrm{c}}=\mathit{KH}\left(K=2.5\right)$$E15

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):

$$\frac{\mathit{\partial }}{\mathit{\partial x}}\left({\mathit{DV}}_x\right)+\frac{\mathit{\partial }}{\mathit{\partial y}}\left({\mathit{DV}}_y\right)+\frac{\mathit{\partial }}{\mathit{\partial \theta }}\left({\mathit{DV}}_{\theta }\right)=F-\mathit{\Phi }$$E16

$$\mathit{\Phi }=\left(K/h\right){\mathit{DC}}_g\left[1-{\left(\mathit{\Gamma }/\gamma \right)}^2\right]\left(\mathit{\Phi }\ge 0\right)$$E17

Here, D is the directional spectrum, (V_x, V_y, 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, C_g is the wave group velocity ($C_g\approx \sqrt{\mathit{gh}}$ in shallow-water wave theory), Γ 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.

Figure 3.

Flowchart for numerical calculation using the Typ. 3 BG model.

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₀ and h_R is considered as follows, similarly to the 3D model of Shimizu et al. [12]:

$$h^{\prime }={\left(\frac{h_{\mathrm{R}}-Z}{h_{\mathrm{R}}+h_0}\right)}^r{h}_0\left(r=1\right)\left(-h_0\le Z\le h_{\mathrm{R}}\right).$$E18

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₂ = 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:

$$P=F\cdot \rho u_m^3$$E19

$$F=\frac{tan{\beta }_{\mathrm{w}}}{tan{\beta }_{\mathrm{c}}}$$E20

$$tan{\beta }_{\mathrm{w}}=\overrightarrow{e_w}•\overrightarrow{\nabla Z}\left(tan{\beta }_w\ge 0\right)$$E21

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₂ = 0 in Eq. (12), Eqs. (22) and (23) are derived for the cross-shore and longshore components of sand transport, q_n and q_s, respectively. Furthermore, under the condition that the seabed slope is equal to the equilibrium slope, Eq. (23) reduces to Eq. (24):

$$q_n=\overrightarrow{e_n}•\overrightarrow{q}=C_0{K}_{n}P\left|cos\alpha \right|\left(\frac{cos\alpha }{\left|cos\alpha \right|}-\frac{tan\beta }{tan{\beta }_{\mathrm{c}}}\right)$$E22

$$q_s=\overrightarrow{e_s}•\overrightarrow{q}=C_0{K}_{s}Psin\alpha \left\{\frac{tan\beta }{tan{\beta }_{\mathrm{c}}}+\frac{K_n}{K_s}\left(1-\frac{tan\beta }{tan{\beta }_{\mathrm{c}}}\right)\right\}$$E23

$$q_s\approx C_0{K}_{s}Psin\alpha \left(\because tan\beta \approx tan{\beta }_{\mathrm{c}}\right)$$E24

In Eq. (22), the cross-shore sand transport q_n becomes 0 when the local seabed slope is equal to the equilibrium slope, and the longshore sand transport q_s 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 $tan{\beta }_{\mathrm{c}}^{\prime }$, which gives q_n = 0 under the condition of oblique wave incidence, satisfies the relation $tan{\beta }_{\mathrm{c}}^{\prime }=tan{\beta }_{\mathrm{c}}cos\alpha$ in the Type 1, as in Eq. (33) in Chapter 1, and $tan{\beta }_{\mathrm{c}}^{\prime }$ is affected by α. On the other hand, in the Type 3, the relation $tan{\beta }_{\mathrm{c}}^{\prime }=tan{\beta }_{\mathrm{c}}$ stands after setting q_n = 0 in Eq. (22). Thus, $tan{\beta }_{\mathrm{c}}^{\prime }$ is independent of α. This is a result of multiplying $\overrightarrow{\nabla Z}$ in the sand transport equation of Eq. (13) in the Type 3 by the coefficient $\left|cos\alpha \right|$, assuming that the effect of gravity is proportional to the magnitude of the cross-slope velocity component that is generated by waves, similarly to the bedload equation of Dronkers [14].

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 (K_s = K_n), and the second term is the additional longshore sand transport in the case that the rates are different (K_s > K_n). 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, u_m, 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:

$$P={\mathit{\Phi }}_{\mathrm{b}}$$E25

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:

$$P={\mathit{\Phi }}_{\mathrm{b}}=f_{D}E=K\sqrt{g/h}\left[1-{\left(\mathrm{\Gamma }/\gamma \right)}^2\right]E\left(P\ge 0\right)$$E26

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 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, q_s, along with Eq. (22) for the cross-shore component of sand transport, q_n. In addition, under the condition that the seabed slope is equal to the equilibrium slope, Eq. (27) reduces to Eq. (28):

$$\begin{array}{l}{q}_s=\overrightarrow{e_s}•\overrightarrow{q}\\ =C_0{K}_{s}P\left\{\frac{tan\beta }{tan{\beta }_{\mathrm{c}}}\left(sin\alpha -\frac{K_2}{K_s}\frac{1}{tan\overline{\beta }}\frac{\mathit{\partial H}}{\mathit{\partial s}}\right)+sin\alpha \frac{K_n}{K_s}\left(1-\frac{tan\beta }{tan{\beta }_{\mathrm{c}}}\right)\right\}\end{array}$$E27

$$q_s=C_0{K}_{s}P\left(sin\alpha -\frac{K_2}{K_s}\frac{1}{tan\overline{\beta }}\frac{\mathit{\partial H}}{\mathit{\partial s}}\right)\left(\because tan\beta \approx tan{\beta }_{\mathrm{c}}\right)$$E28

The cross-shore sand transport q_n is the same as Eq. (22) in the Type 3, whereas the longshore sand transport q_s given by Eqs. (27) and (28) has the additional term of Ozasa and Brampton [1]. When the total sand transport Q_s 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 $P={\mathit{\Phi }}_{\mathrm{b}}$ (Eq. (25)) is employed, α 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:

$$\begin{array}{l}{Q}_s=\int q_s\mathit{dn}\\ =C_0{K}_s\int \left\{P\left(sin\alpha -\frac{K_2}{K_s}\frac{1}{tan\overline{\beta }}\frac{\mathit{\partial H}}{\mathit{\partial s}}\right)\right\}\mathit{dn}\\ =C_0{K}_s\int \left\{{\mathit{\Phi }}_{\mathrm{b}}\left(sin\alpha -\frac{K_2}{K_s}\frac{1}{tan\overline{\beta }}\frac{\mathit{\partial H}}{\mathit{\partial s}}\right)\right\}\mathit{dn}\left(\because P={\mathit{\Phi }}_{\mathrm{b}}\right)\\ =C_0{K}_s\left\{sin{\alpha }_{\mathrm{b}}-\frac{K_2}{K_s}\frac{1}{tan\overline{\beta }}\frac{\mathit{\partial }{H}_{\mathrm{b}}}{\mathit{\partial s}}\right\}\int {\mathit{\Phi }}_{\mathrm{b}}\mathit{dn}\left(\because \alpha \approx {\alpha }_{\mathrm{b}},\frac{\mathit{\partial H}}{\mathit{\partial s}}\approx \frac{\mathit{\partial }{H}_{\mathrm{b}}}{\mathit{\partial s}}\right)\\ =C_0{K}_s{\left({\mathit{EC}}_g\right)}_{\mathrm{b}}cos{\alpha }_{\mathrm{b}}\left\{sin{\alpha }_{\mathrm{b}}-\frac{K_2}{K_s}\frac{1}{tan\overline{\beta }}\frac{\mathit{\partial }{H}_{\mathrm{b}}}{\mathit{\partial s}}\right\}\\ \left(\because \int {\mathit{\Phi }}_{\mathrm{b}}\mathit{dn}\approx {\left({\mathit{EC}}_g\right)}_{\mathrm{b}}cos{\alpha }_{\mathrm{b}}\right)\end{array}$$E29

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:

$$\overrightarrow{q}=C_0\frac{K_{s}P}{tan{\beta }_{\mathrm{c}}}\left\{tan{\beta }_{\mathrm{c}}\overrightarrow{e_w}-\left|cos\alpha \right|\overrightarrow{\nabla Z}\right\}\left(-h_{\mathrm{c}}\le Z\le h_{\mathrm{R}}\right)$$E30

$$P=\epsilon \left(Z\right){\left({\mathit{EC}}_g\right)}_{\mathrm{b}}tan{\beta }_w\left(P\ge 0\right)$$E31

$$tan{\beta }_w=\mathit{dZ}/{\mathit{dx}}_w\left(tan{\beta }_w\ge 0\right)$$E32

$$I_{\epsilon }\left(Z\right)=\left\{\begin{array}{l}1\left(Z<-h_{\mathrm{c}}\right)\\ {\int }_Z^{h_{\mathrm{R}}}\epsilon \left(Z\right)\mathit{dZ}=\frac{h_{\mathrm{R}}-Z}{h_{\mathrm{c}}+h_{\mathrm{R}}}\left(-h_{\mathrm{c}}\le Z\le h_{\mathrm{R}}\right)\\ 0\left(Z>h_{\mathrm{R}}\right)\end{array}\right.$$E33

$$P={\left({\mathit{EC}}_g\right)}_{\mathrm{b}}\left(-\frac{{\mathit{dI}}_{\epsilon }}{{\mathit{dx}}_w}\right)\left(P\ge 0\right)$$E34

$${I_{\epsilon }}^{\left(i+1\right)}=min\left({\left.\mathrm{Eq}.\left(33\right)\right|}^{Z=Z^{\left(i+1\right)}},{I_{\epsilon }}^{\left(i\right)}\right)$$E35

$$P^{\left(i+1/2\right)}={\left({\mathit{EC}}_g\right)}_{\mathrm{b}}\left(\frac{{I_{\epsilon }}^{\left(i\right)}-{I_{\epsilon }}^{\left(i+1\right)}}{\mathrm{\Delta }{x}_w}\right)$$E36

Here, K_s is the longshore and cross-shore sand transport coefficients, and the index i in Eqs. (35) and (36) is the mesh number along the x_w-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 x_w- and y_w-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 (x_w, y_w). 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 x_w-axis is subdivided by a mesh of Δx_w, and a cumulative function of ε(Z), $I_{\epsilon }\left(Z\right)$, is introduced as in Eq. (33). Here, ε(Z) is assumed to have a uniform distribution (Eq. (7)), and $I_{\epsilon }\left(Z\right)$ takes a value of unity in the zone deeper than h_c, decreases with the water depth, and becomes 0 in the zone higher than h_R.

Figure 4.

Definition of coordinate system around an island.

Using Eq. (33), Eq. (31) is transformed into Eq. (34). Thus, the P value can be determined from the derivative of $I_{\epsilon }\left(Z\right)$ at a point along the x_w-axis using Eq. (34). $I_{\epsilon }\left(Z\right)$ is calculated along the x_w-axis from the starting point of wave incidence in the direction of wave propagation; $I_{\epsilon }\left(Z\right)$ at the (i + 1)^th point can be calculated from Eq. (35) with the given value of $I_{\epsilon }\left(Z\right)$ at the i^th point when the initial value of $I_{\epsilon }\left(Z\right)$ at the offshore end is given and the mesh location is denoted by ${x_w}^{\left(i\right)}=i\mathrm{\Delta }{x}_w$.

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_{\epsilon }\left(Z\right)$ at the (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_{\epsilon }\left(Z\right)$ at the i^th point is adopted. As a result, the value of $I_{\epsilon }\left(Z\right)$ corresponding to the minimum water depth (maximum Z) between the offshore end and a designated point along the x_w-axis can be adopted.

Using this procedure, the depth distribution along the x_w-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 x_w-axis has several uneven shapes.

Consider the case where the x_w-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 x_w-axis is always equivalent to (EC_g)_b regardless of the seabed topography. Because (EC_g)_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 x_w-axis. The P value calculated at point (x_w, y_w) 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 (x_w, y_w). The mesh intervals Δx_w and Δy_w 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.

Figure 5.

Flowchart for numerical calculation using Typ. 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:

$$\overrightarrow{q}=C_0\frac{K_{s}P}{tan{\beta }_{\mathrm{c}}}\left\{tan{\beta }_{\mathrm{c}}\overrightarrow{e_w}-\left|cos\alpha \right|\overrightarrow{\nabla Z}\right\}\left(-h_{\mathrm{c}}\le Z\le h_{\mathrm{R}}\right)$$E37

$$F^{\left(i+1\right)}=F^{\left(i\right)}+r\mathrm{\Delta }{x}_w$$E38

$$r=\left\{\begin{array}{l}1\left(Z\le 0\right)\\ 0\left(Z>0\right)\end{array}\right.$$E39

$$F^{\left(i\right)}=0\text{if}\left(Z\ge 0\text{and}\mathit{dZ}/{\mathit{dx}}_w\le 0\right)$$E40

$$\begin{array}{l}{H}_{1/3}=f\left(F,U\right)\\ =0.30\left\{1-{\left[1+0.004{\left(\mathit{gF}/U^2\right)}^{1/2}\right]}^{-2}\right\}\left(U^2/g\right)\end{array}$$E41

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 x_w-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 (x_w, y_w) 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.

Figure 6.

Definition of two coordinate systems around a lake.

Neglecting the wave refraction effect, waves are assumed to propagate in the same direction as the wind. The x_w-axis is subdivided by mesh intervals of Δx_w. The fetch F is added from upwind to downwind along the x_w-axis using Eq. (38). When a grid point is located on land and the downslope condition of dZ/dx_w ≤ 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 (x_w, y_w). The calculation of the P value is independently carried out along each x_w-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.

Figure 7.

Flowchart for numerical calculation using Typ. 6 BG model.

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 $\overrightarrow{q}$ is a linear sum of the sand transport due to waves $\overrightarrow{q_w}$ and that due to ebb-tidal currents $\overrightarrow{q_R}$,

$$\overrightarrow{q}=\overrightarrow{q_w}+\overrightarrow{q_R}.$$E42

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 $\overrightarrow{q_w}=(q_{\mathit{wx}},q_{\mathit{wy}})$ is written as

$$q_{\mathit{wx}}=\frac{G_{\mathit{wx}}}{tan{\beta }_{\mathrm{c}}}\left[tan{\beta }_{\mathrm{c}}cos{\theta }_{\mathrm{w}}-\frac{\mathit{\partial Z}}{\mathit{\partial x}}\right],$$E43a

$$q_{\mathit{wy}}=\frac{G_{\mathit{wy}}}{tan{\beta }_{\mathrm{c}}}\left[tan{\beta }_{\mathrm{c}}sin{\theta }_{\mathrm{w}}-\frac{\mathit{\partial Z}}{\mathit{\partial y}}\right].$$E43b

Here, G_wx and G_wy are given by

$$G_{\mathit{wx}}=C_0{K}_{x}P$$E44a

$$G_{\mathit{wy}}=C_0{K}_{y}P$$E44b

Here, q_x is the cross-shore component of the sediment transport flux (positive for shoreward), and q_y is the longshore component of sediment transport flux. K_x and K_y 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, K_x and K_y, 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 $\overrightarrow{q_R}=(q_{\mathit{Rx}},q_{\mathit{Ry}})$ due to ebb-tidal currents, Eq. (45) of Bailard and Inman [8] is employed, which was derived in terms of the seabed slope from the bedload transport formula of Bagnold [7] by a linear approximation:

$$q_{\mathit{Rx}}=\frac{G_R}{tan\varphi }\left[tan\varphi cos{\theta }_R-\frac{\mathit{\partial Z}}{\mathit{\partial x}}\right]$$E45a

$$q_{\mathit{Ry}}=\frac{G_R}{tan\varphi }\left[tan\varphi sin{\theta }_R-\frac{\mathit{\partial Z}}{\mathit{\partial y}}\right]$$E45b

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 G_R 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 G_R 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 G_R is given by

$$G_R=\left\{\begin{array}{l}{C}_0{K}_R{F}_w{K}_V^3\left(-h_{\mathrm{c}2}\le Z\le h_{\mathrm{R}2}\right)\\ 0\left(Z<-h_{\mathrm{c}2},h_{\mathrm{R}2}<Z\right)\end{array}\right.$$E46

$$F_w=\frac{{\left({\mathit{EC}}_g\right)}_{\mathrm{bo}}}{h_{\mathrm{c}}+h_{\mathrm{R}}}$$E47

$$K_V=\frac{V}{V_0}=K_{V1}\left(\frac{h_0}{h}\right)$$E48

$$K_{V1}=\left(\frac{V_1}{V_0}\right).$$E49

Here, K_R is the coefficient of sand transport due to ebb-tidal currents, and K_V is the ratio of the ebb-tidal current velocity to V₀, V/V₀, where V is the velocity of the ebb-tidal currents and V₀ is a reference velocity at the inlet. F_w 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, (EC_g)_b0, divided by the sum of h_c and h_R, as in Eq. (47). V₁ 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₀ 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 G_R is given as Eq. (46), G_R 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); G_R is expressed by the intensity of the ebb-tidal currents relative to the wave intensity (EC_g)_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, G_R 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 K_R 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.

Figure 8.

Flowchart for numerical calculation using Typ. 7 BG model.

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 $\overrightarrow{q}$ = (q_x, q_y) is a linear sum of the component due to waves, $\overrightarrow{q_w}$ = (q_wx, q_wy), and that due to currents, $\overrightarrow{q_c}$ = (q_cx, q_cy):

$$\overrightarrow{q}=\overrightarrow{q_w}+\overrightarrow{q_c}\left(-h_{\mathrm{c}}\le Z\le h_{\mathrm{R}}\right)$$E50

$$\overrightarrow{q_w}=C_0\frac{P_w}{tan{\beta }_{\mathrm{c}}}\left(tan{\beta }_{\mathrm{c}}\overrightarrow{e_w}-b_1{\displaystyle \overrightarrow{\nabla Z}}\right)$$E51

$$P_w=K_w\cdot \rho {u_m}^3$$E52

$$\overrightarrow{q_c}=C_0\frac{P_c}{tan\varphi }\left(tan\varphi \overrightarrow{e_c}-b_2{\displaystyle \overrightarrow{\nabla Z}}\right)$$E53

$$P_c=K_c\cdot \rho {u_m}^{2}V$$E54

Here, $\overrightarrow{q}$ is the total sand transport, $\overrightarrow{q_w}$ is the sand transport due to waves, $\overrightarrow{q_c}$ is the sand transport due to currents, $\overrightarrow{e_c}$ is the unit vector in the current direction, and u_m is the amplitude of the seabed velocity due to the orbital motion of waves defined by Eq. (14) in the Type 3 BG model. $b_1=\left|cos{\alpha }_1\right|=\left|{\displaystyle \overrightarrow{e_w}}•{\displaystyle \overrightarrow{\nabla Z}}\right|/\left|{\displaystyle \overrightarrow{\nabla Z}}\right|$, α₁ is the angle between the wave direction and the direction normal to the contour lines, V is the velocity of the currents, K_w and K_c are the coefficients of sand transport due to waves and currents, respectively, $b_2=\left|cos{\alpha }_2\right|=\left|{\displaystyle \overrightarrow{e_c}}•{\displaystyle \overrightarrow{\nabla Z}}\right|/\left|{\displaystyle \overrightarrow{\nabla Z}}\right|$, α₂ 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₂ 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²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₀ and h_R. In particular, in the calculation in Chapter 4, h₀ 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.5H, where H is the wave height at a local point, the increment due to the effect of the strong rip currents is evaluated as

$${h_{\mathrm{c}}}^{\prime }=\left\{1+a\left(V/u_m\right)\right\}{h}_{\mathrm{c}}.$$E55

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.

Figure 9.

Flowchart for numerical calculation using Typ. 8 BG model.

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.

Figure 10.

Calculation domain and meshes.

First, the calculation point of the elevation Z at $(x_i,y_j)=\left(i\mathrm{\Delta }x,j\mathrm{\Delta }y\right)$ is taken at the center of the cell, as shown in Figure 11, and $Z_{\left(i,j\right)}=Z(x_i,y_j,t)$ (1 ≤ i ≤ N_x, 1 ≤ j ≤ N_y) is evaluated at each point. Here, the subscripts i and j are the cell numbers taken in the x- and y-directions, and N_x and N_y are the numbers of cells in the x- and y-directions, respectively. Then, $q_{x\left(i-1/2,j\right)}$ (1 ≤ i ≤ N_x + 1, 1 ≤ j ≤ N_y) and $q_{y\left(i,j-1/2\right)}$ (1 ≤ i ≤ N_x, 1 ≤ j ≤ N_y + 1) are set at points $\left(\left(i-1/2\right)\mathrm{\Delta }x,j\mathrm{\Delta }y\right)$ and $\left(i\mathrm{\Delta }x,\left(j-1/2\right)\mathrm{\Delta }y\right)$, which are separated from the points Z_{(i, j)} by 1/2 mesh in the x- and y-directions, respectively.

Figure 11.

Arrangement of variables in discretization.

When the sand transport flux $\overrightarrow{q}=(q_x,q_y)$ is calculated, the seabed elevation $Z_{\left(i,j\right)}^{\prime }=Z(x_i,y_j,t+\mathrm{\Delta }t)$ after time Δt can be calculated using Eq. (56), which is the discretized form of Eq. (29) in Chapter 1:

$$\begin{array}{l}{Z}_{(i,j)}^{\prime }=Z_{(}+\left\{q_{x\left(i-1/2,j\right)}-q_{x\left(i+1/2,j\right)}\right\}\left(\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}\right)\\ +\left\{q_{y\left(i,j-1/2\right)}-q_{y\left(i,j+1/2\right)}\right\}\left(\frac{\mathrm{\Delta }t}{\mathrm{\Delta }y}\right)\end{array}$$E56

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 ${\mathit{dq}}_x/\mathit{dx}=0$ or ${\mathit{dq}}_y/\mathit{dy}=0$ is often set, where sand freely enters or passes through the boundary. In this case, the value at the boundary is evaluated to be the same value as the sand transport evaluated at an inner point.

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 (q_x, q_y) 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

$$\overrightarrow{q}=A_1\overrightarrow{e_w}+A_2\overrightarrow{\nabla Z}+A_{3}tan\beta \overrightarrow{e_s},$$E57

where $\overrightarrow{e_w}$, $\overrightarrow{\nabla Z}$, and $tan\beta \overrightarrow{e_s}$ are given by Eqs. (3), (9), and (11) in Chapter 1. When using the vector expression of Eq. (57), the coefficients of the sand transport equation of the Type 1 (Eq. (58)) are given by Eq. (59):

$$\overrightarrow{q}=C_0\frac{K_{1}P}{tan{\beta }_{\mathrm{c}}}\left[tan{\beta }_{\mathrm{c}}\overrightarrow{e_w}-\overrightarrow{\nabla Z}\right]$$E58

$$A_1=C_0{K}_{1}P$$E59a

$$A_2=-\frac{A_1}{tan{\beta }_{\mathrm{c}}}$$E59b

$$A_3=0$$E59c

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

$$\begin{array}{l}\overrightarrow{q}=C_0\frac{P}{tan{\beta }_{\mathrm{c}}}\left\{\begin{array}{l}{K}_n\left(tan{\beta }_{\mathrm{c}}\overrightarrow{e_w}-\left|cos\alpha \right|\overrightarrow{\nabla Z}\right)\\ +\left\{\left(K_s-K_n\right)sin\alpha -\frac{K_2}{tan\overline{\beta }}\frac{\mathit{\partial H}}{\mathit{\partial s}}\right\}tan\beta \overrightarrow{e_s}\end{array}\right\}\\ \left(-h_{\mathrm{c}}\le Z\le h_{\mathrm{R}}\right)\end{array}$$E60

$$A_1=C_0{K}_{n}P$$E61a

$$A_2=-\frac{\left|cos\alpha \right|}{tan{\beta }_{\mathrm{c}}}{A}_1$$E61b

$$A_3=\frac{1}{K_{n}tan{\beta }_{\mathrm{c}}}\left\{\left(K_s-K_n\right)sin\alpha -\frac{K_2}{tan\overline{\beta }}\frac{\mathit{\partial H}}{\mathit{\partial s}}\right\}{A}_1$$E61c

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:

$$q_x=A_{1}cos{\theta }_w+A_2\left(\mathit{\partial Z}/\mathit{\partial x}\right)-A_3\left(\mathit{\partial Z}/\mathit{\partial y}\right)$$E62a

$$q_y=A_{1}sin{\theta }_w+A_2\left(\mathit{\partial Z}/\mathit{\partial y}\right)+A_3\left(\mathit{\partial Z}/\mathit{\partial x}\right)$$E62b

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

$q_{x\left(i-1/2,j\right)}$ is calculated as

$$\begin{array}{l}{q}_{x\left(i-1/2,j\right)}=A_{1\left(i-1/2,j\right)}cos{\theta }_{w\left(i-1/2,j\right)}\\ +A_{2\left(i-1/2,j\right)}{\left(\mathit{\partial Z}/\mathit{\partial x}\right)}_{\left(i-1/2,j\right)}\\ -A_{3\left(i-1/2,j\right)}{\left(\mathit{\partial Z}/\mathit{\partial y}\right)}_{\left(i-1/2,j\right)},\end{array}$$E63

where

$${\left(\mathit{\partial Z}/\mathit{\partial x}\right)}_{\left(i-1/2,j\right)}=\frac{Z_{\left(i,j\right)}-Z_{\left(i-1,j\right)}}{\mathrm{\Delta }x},$$E64

$$\begin{array}{l}{\left(\mathit{\partial Z}/\mathit{\partial y}\right)}_{\left(i-1/2,j\right)}=\frac{1}{2}\left\{{\left(\mathit{\partial Z}/\mathit{\partial y}\right)}_{\left(i-1,j\right)}+{\left(\mathit{\partial Z}/\mathit{\partial y}\right)}_{\left(i,j\right)}\right\}\\ =\frac{Z_{\left(i-1,j+1\right)}-Z_{\left(i-1,j-1\right)}+Z_{\left(i,j+1\right)}-Z_{\left(i,j-1\right)}}{4\mathrm{\Delta }y}.\end{array}$$E65

$q_{x\left(i+1/2,j\right)}$ is calculated as

$$q_{x\left(i+1/2,j\right)}={\left[\mathrm{Eq}.\left(63\right)\right]}^{i=i+1}.$$E66

$q_{y\left(i,j-1/2\right)}$ is calculated as

$$\begin{array}{l}{q}_{y\left(i,j-1/2\right)}=A_{1\left(i,j-1/2\right)}sin{\theta }_{w\left(i,j-1/2\right)}\\ +A_{2\left(i,j-1/2\right)}{\left(\mathit{\partial Z}/\mathit{\partial y}\right)}_{\left(i,j-1/2\right)}\\ +A_{3\left(i,j-1/2\right)}{\left(\mathit{\partial Z}/\mathit{\partial x}\right)}_{\left(i,j-1/2\right)},\end{array}$$E67

where

$${\left(\mathit{\partial Z}/\mathit{\partial y}\right)}_{\left(i,j-1/2\right)}=\frac{Z_{\left(i,j\right)}-Z_{\left(i,j-1\right)}}{\mathrm{\Delta }y},$$E68

$$\begin{array}{l}{\left(\mathit{\partial Z}/\mathit{\partial x}\right)}_{\left(i,j-1/2\right)}=\frac{1}{2}\left\{{\left(\mathit{\partial Z}/\mathit{\partial x}\right)}_{\left(i,j-1\right)}+{\left(\mathit{\partial Z}/\mathit{\partial x}\right)}_{\left(i,j\right)}\right\}\\ =\frac{Z_{\left(i+1,j-1\right)}-Z_{\left(i-1,j-1\right)}+Z_{\left(i+1,j\right)}-Z_{\left(i-1,j\right)}}{4\mathrm{\Delta }x}.\end{array}$$E69

Finally, $q_{y\left(i,j+1/2\right)}$ is calculated as

$$q_{y\left(i,j+1/2\right)}={\left[\mathrm{Eq}.\left(67\right)\right]}^{j=j+1}.$$E70

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 $\overrightarrow{q}=(q_x,q_y)$ are reduced by multiplying by a coefficient μ as in Eq. (71), where μ is given by Eq. (72). ΔZ₀ 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.

Figure 12.

Reduction coefficient μ of sand tranport rate near berm height h_R and depth of closure h_c.

$${\displaystyle \overrightarrow{q^{\prime }}}=\mu \overrightarrow{q}$$E71

$$\mu =\mu \left(Z\right)=\frac{\mathrm{\Delta }Z}{\mathrm{\Delta }{Z}_0}\left(0\le \mu \le 1\right)$$E72

$$\mathrm{\Delta }Z=h_{\mathrm{R}}-Z\left(\text{near}\mathrm{Z}=h_{\mathrm{R}}\right)$$E73

$$\mathrm{\Delta }Z=Z-\left(-h_{\mathrm{c}}\right)\left(\text{near}Z=-h_{\mathrm{c}}\right)$$E74

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:

$$q_{x\left(i-1/2,j\right)}^{\prime }={\mu }_{\mathit{iup}}{q}_{x\left(i-1/2,j\right)}\left(\text{near}Z=-h_{\mathrm{c}}\right)$$E75a

$${\mu }_{\mathit{iup}}=\frac{\mathrm{\Delta }{Z}_{\mathit{iup}}}{\mathrm{\Delta }{Z}_0}=\frac{Z_{\left(\mathit{iup},j\right)}-\left(-h_{\mathrm{c}}\right)}{\mathrm{\Delta }{Z}_0}$$E75b

$$\mathit{iup}=\left\{\begin{array}{l}i-1\left(q_{x\left(i-1/2,j\right)}\ge 0\right)\\ i\left(q_{x\left(i-1/2,j\right)}<0\right)\end{array}\right.$$E75c

Figure 13.

Method selecting Z value near h_c depending on the direction of sand transport q_x when calculating reduction coefficient μ of sand transport rate.

Here, the subscript iup represents the location of the Z value at a point immediately updrift of q_x (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).

Figure 14.

Method selecting Z value near h_R depending on the direction of sand transport q_x when calculating reduction coefficient μ of sand transport rate.

$$q_{x\left(i-1/2,j\right)}^{\prime }={\mu }_{\mathit{\text{idown}}}{q}_{x\left(i-1/2,j\right)}\left(\text{near}Z=h_{\mathrm{R}}\right)$$E76a

$${\mu }_{\mathit{\text{idown}}}=\frac{\mathrm{\Delta }{Z}_{\mathit{\text{idown}}}}{\mathrm{\Delta }{Z}_0}=\frac{h_{\mathrm{R}}-Z_{\left(\mathit{\text{idown}},j\right)}}{\mathrm{\Delta }{Z}_0}$$E76b

$$\mathit{\text{idown}}=\left\{\begin{array}{l}i\left(q_{x\left(i-1/2,j\right)}\ge 0\right)\\ i‐1\left(q_{x\left(i-1/2,j\right)}<0\right)\end{array}\right.$$E76c

Here, the subscript idown represents the location of the Z value at a point immediately downdrift of q_x (Figure 14). A similar method can be used for the y-component of q. Here, the thickness of the sand layer ΔZ₀, 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₀ 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]:

$$\mathrm{\Delta }{Z}_0=\frac{1}{2}tan{\beta }_{\mathrm{c}}\cdot \mathrm{\Delta }L\left(\mathrm{\Delta }L=min\left(\mathrm{\Delta }x,\mathrm{\Delta }y\right)\right)$$E77

In the other calculation, ΔZ₀ is selected as a value within the following range:

$$\mathrm{\Delta }{Z}_0=\left(0.05‐0.2\right)h_{\mathrm{c}}.$$E78

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 $\overrightarrow{q_g}$ due to the effect of gravity on a steep slope is calculated by Eq. (79) (Figure 15), given the slope of the angle of repose of sand.

Figure 15.

Seaward sand transport due to gravity when local slope exceeds slope of repose angle.

$$\overrightarrow{q_g}=A_g\left(tan\beta -tan\varphi \right)\left(-{\displaystyle \overrightarrow{e_n}}\right)\left(tan\beta \ge tan\varphi \right)$$E79

Here, $\overrightarrow{e_n}$ is the unit vector normal to the contour lines (shoreward), where n is the local coordinate taken along the direction normal (shoreward) to the contour lines, $tan\beta =\left|{\displaystyle \overrightarrow{\nabla Z}}\right|$ is the local slope, and $tan\varphi$ is the slope of the angle of repose of sand particle. Here, we assume $tan\varphi =1/2$. Sand transport with flux $\overrightarrow{q_g}$ takes place when the local slope exceeds the slope of the angle of repose of sand, the direction of $\overrightarrow{q_g}$ is down the slope and normal to the contour lines, and the strength is proportional to the deviation between the local slope and the slope of the angle of repose with proportionality coefficient A_g. 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, $\overrightarrow{q_g}$ is set to 0.

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

$$\overrightarrow{\nabla Z}=tan\beta \overrightarrow{e_n}=\partial Z/\partial n$$E80

$$\overrightarrow{q_g}=A_g\left(tan\varphi /tan\beta -1\right){\displaystyle \overrightarrow{\nabla Z}}\left(tan\beta \ge tan\varphi \right)$$E81

Eq. (81) can be expressed by the basic expression for sand transport of Eq. (57), as mentioned earlier, and the coefficients A₁, A₂, and A₃ of each term of Eq. (57) are given by

$$A_1=A_3=0$$E82a

$$A_2=A_g\left(tan\varphi /tan\beta -1\right).$$E82b

Accordingly, $\overrightarrow{q_g}$ can be evaluated in the same manner as Eq. (57) using the coefficients given by Eq. (82), assuming that the coefficient 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 A_g. However, when an excessively large value is employed for A_g, it is difficult to stably carry out the numerical calculation. Here, A_g is designated as the maximum value in a range, in which the numerical simulation can stably proceed as follows.

First, $\overrightarrow{q_g}$ in Eq. (79) is reduced to Eq. (83) using Eq. (80), and the cross-shore sand transport component q_n 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 A_g. To solve Eq. (86) using the explicit finite-difference method, Eq. (87) must be satisfied as the stability condition of Eq. (86):

$$\overrightarrow{q_g}=A_g\left(tan\varphi -\partial Z/\partial n\right)\overrightarrow{e_n}$$E83

$$q_n=\overrightarrow{e_n}•\overrightarrow{q_g}=A_g\left(tan\varphi -\partial Z/\partial n\right)$$E84

$$\frac{\mathit{\partial Z}}{\mathit{\partial t}}=-\frac{\mathit{\partial }{q}_n}{\mathit{\partial n}}$$E85

$$\frac{\mathit{\partial Z}}{\mathit{\partial t}}=A_g\frac{\mathit{\partial }{Z}^2}{\mathit{\partial }{n}^2}$$E86

$$\mathrm{\Delta }t\le 0.5\frac{{\left(\mathrm{\Delta }n\right)}^2}{A_g}$$E87

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 A_g and Δn, here it is employed as a relationship to determine the upper limit of A_g given Δn and Δt. Eq. (87) can be rewritten as Eq. (88) to provide a stability condition for A_g. Thus, the right term of Eq. (88) gives the upper limit of A_g:

$$A_g\le 0.5\frac{{\left(\mathrm{\Delta }n\right)}^2}{\mathrm{\Delta }t}$$E88

The coefficient A_g that ensures the stability of Eq. (89) is given:

$$A_g=R_S\frac{{\left(\mathrm{\Delta }L\right)}^2}{\mathrm{\Delta }t}$$E89a

$$R_S=0.2$$E89b

$$\mathrm{\Delta }L=min\left(\mathrm{\Delta }x,\mathrm{\Delta }y\right).$$E89c

Finally, $\overrightarrow{q_g}$ is calculated in the same manner as Eq. (57) using Eqs. (82) and (89). In the numerical simulation, sand transport fluxes due to both waves and the effect of gravity are calculated at each time step, and the larger value at a location is employed as in Eq. (90):

$$\overrightarrow{q}=\left\{\begin{array}{l}{\displaystyle \overrightarrow{q_g}}\left(\text{if}\left|{\displaystyle \overrightarrow{q_g}}\right|>\left|{{\displaystyle \overrightarrow{q}}}_{\mathrm{due}\text{to waves}}\right|\right)\\ \\ {{\displaystyle \overrightarrow{q}}}_{\mathrm{due}\text{to waves}}\left(\text{otherwise}\right)\end{array}\right.$$E90