Open access peer-reviewed chapter

Eight Types of BG Models and Discretization

By Takaaki Uda, Masumi Serizawa and Shiho Miyahara

Submitted: November 11th 2016Reviewed: September 10th 2018Published: December 19th 2018

DOI: 10.5772/intechopen.81412

Downloaded: 225

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

TypeExternal forceCharacteristicsApplicationChapter
1WavesMost fundamental model using wave parameters at breaking pointPrediction of typical beach changes in coastal engineering3
2WavesThe 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 activities3
3WavesThe 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 fieldFormation of sand spit and bay barrier5
4WavesThe intensity of sand transport P is given by the wave energy dissipation rate due to wave breaking at a local pointFormation of cuspate foreland7
5WavesThe 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 directionFormation of sand spit and bay barrier and interaction of sandy islands on flat shallow seabed owing to waves5 and 6
6Wind waves developed in closed water bodyThe height of wind waves is predicted using Wilson’s formula given the wind fetch distance and wind velocitySegmentation and merging of closed water bodies by wind waves8
7Waves and ebb-tidal currentsA 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 coefficientFormation of dynamically stable ebb-tidal delta4
8Waves and nearshore currentsThe 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 accountBeach changes around artificial reef4

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β¯. Φ is given by the wave energy flux transported per unit width of the shoreline, ECgbcosαb, divided by the cross-shore distance R=hc+hR/tanβ¯between the depth of closure, hc, and the berm height, hR, corresponding to the wave run-up height (Figure 1). Then, assuming tanβ¯tanβccosαb(Eq. (32) in Chapter 1), Φ becomes ECgbcos2αb/hc+hR/tanβc. Furthermore, introducing the depth distribution of the sand transport intensity ε(Z) and transforming 1/hc+hRinto ε(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

q=GtanβctanβcewZ.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):

qx=GtanβctanβccosθwZ/xE2a
qy=GtanβctanβcsinθwZ/yE2b
G=C0K1PE3
P=Φ=εZECgbcos2αbtanβcE4
hchRεZdZ=1E5
εZ=2hc3hc2ZZ+hc2hcZhR0Z<hchR<ZE6

Here, ε(Z) is defined so that the integral over the depth between Z = −hc and hR 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:

εZ=1/hc+hRhcZhR0Z<hchR<ZE7

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 Hb, where Hb is given by the following relation:

ECgb=C1Hb52E8a
C1=ρgk1g/γγ0.8E8b
k1=8for regular waves4.0042for irregular wavesE8c

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

αbα=θwθn=θwtan1∂Z∂y/∂Z∂xE9

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 q=(qx,qy)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 Hb, αb, hc, and hR, 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 cosαbin 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 Kd and the direction of diffracted waves θd are predicted using the angular spreading method for irregular waves [5, 6], and (ECg)b in the case of no offshore structure is reduced by multiplying by the square of the coefficient Kd so that a relation of (ECg)b′ = Kd2(ECg)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, K2/K1, and the longshore gradient of the wave height Hb multiplied by Kd. 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):

qx=GxtanβctanβccosθwZ/xhcZhRE10a
qy=Gytanβctanβcsinθw1tanβ¯K2KyHb∂yZ/yhcZhRE10b
Gx=C0KxPE11a
Gy=C0KyPE11b

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, K2 is the coefficient of the term given by Ozasa and Brampton [1], and tanβ¯is the seabed slope at the breaker point. Here, we assume tanβ¯=tanβ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 hR and hc, the Kd 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

q=C0PtanβcKntanβcewcosαZ+KsKnsinαK2tanβ¯∂H∂stanβeshcZhRE12
P=ρum3E13
um=H2gh.E14

Here, tanβes=∂Z/∂y∂Z/∂x, Ks and Kn are the coefficients of longshore sand transport and cross-shore sand transport, respectively, ∂H/∂s=esHis the longshore gradient of the wave height H measured parallel to the contour lines, and H=∂H/∂x∂H/∂y. tanβ¯is the characteristic slope of the breaker zone, and h is the water depth. Furthermore, tanβ¯=tanβcis assumed. 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. hc 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]:

hc=KHK=2.5E15

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

∂xDVx+∂yDVy+∂θDVθ=FΦE16
Φ=K/hDCg1Γ/γ2Φ0E17

Here, D is the directional spectrum, (Vx, Vy, ) 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 (Cgghin 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 h0 and hR is considered as follows, similarly to the 3D model of Shimizu et al. [12]:

h=hRZhR+h0rh0r=1h0ZhR.E18

In addition, at locations whose elevation is higher than hR, 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., K2 = 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ρum3E19
F=tanβwtanβcE20
tanβw=ewZtanβw0E21

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

qn=enq=C0KnPcosαcosαcosαtanβtanβcE22
qs=esq=C0KsPsinαtanβtanβc+KnKs1tanβtanβcE23
qsC0KsPsinαtanβtanβcE24

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 tanβc, which gives qn = 0 under the condition of oblique wave incidence, satisfies the relation tanβc=tanβccosαin the Type 1, as in Eq. (33) in Chapter 1, and tanβcis affected by α. On the other hand, in the Type 3, the relation tanβc=tanβcstands after setting qn = 0 in Eq. (22). Thus, tanβcis independent of α. This is a result of multiplying Zin the sand transport equation of Eq. (13) in the Type 3 by the coefficient cosα, 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 (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:

P=ΦbE25

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=Φb=fDE=Kg/h1Γ/γ2EP0E26

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

qs=esq=C0KsPtanβtanβcsinαK2Ks1tanβ¯∂H∂s+sinαKnKs1tanβtanβcE27
qs=C0KsPsinαK2Ks1tanβ¯∂H∂stanβtanβcE28

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 P=Φ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:

Qs=qsdn=C0KsPsinαK2Ks1tanβ¯∂H∂sdn=C0KsΦbsinαK2Ks1tanβ¯∂H∂sdnP=Φb=C0KssinαbK2Ks1tanβ¯Hb∂sΦbdnααb∂H∂sHb∂s=C0KsECgbcosαbsinαbK2Ks1tanβ¯Hb∂sΦbdnECgbcosαbE29

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:

q=C0KsPtanβctanβcewcosαZhcZhRE30
P=εZECgbtanβwP0E31
tanβw=dZ/dxwtanβw0E32
IεZ=1Z<hcZhRεZdZ=hRZhc+hRhcZhR0Z>hRE33
P=ECgbdIεdxwP0E34
Iεi+1=minEq.33Z=Zi+1IεiE35
Pi+1/2=ECgbIεiIεi+1ΔxwE36

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), IεZ, is introduced as in Eq. (33). Here, ε(Z) is assumed to have a uniform distribution (Eq. (7)), and IεZtakes a value of unity in the zone deeper than hc, decreases with the water depth, and becomes 0 in the zone higher than hR.

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εZat a point along the xw-axis using Eq. (34). IεZis calculated along the xw-axis from the starting point of wave incidence in the direction of wave propagation; IεZat the (i + 1)th point can be calculated from Eq. (35) with the given value of IεZat the ith point when the initial value of IεZat the offshore end is given and the mesh location is denoted by xwi=iΔxw.

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εZat 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εZat the ith point is adopted. As a result, the value of IεZcorresponding to the minimum water depth (maximum 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 hR, 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 hR 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.

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

q=C0KsPtanβctanβcewcosαZhcZhRE37
Fi+1=Fi+rΔxwE38
r=1Z00Z>0E39
Fi=0ifZ0anddZ/dxw0E40
H1/3=fFU=0.3011+0.004gF/U21/22U2/gE41

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.

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 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.

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 qis a linear sum of the sand transport due to waves qwand that due to ebb-tidal currents qR,

q=qw+qR.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 qw=(qwx,qwy)is written as

qwx=Gwxtanβctanβccosθw∂Z∂x,E43a
qwy=Gwytanβctanβcsinθw∂Z∂y.E43b

Here, Gwx and Gwy are given by

Gwx=C0KxPE44a
Gwy=C0KyPE44b

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 qR=(qRx,qRy)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:

qRx=GRtanϕtanϕcosθR∂Z∂xE45a
qRy=GRtanϕtanϕsinθR∂Z∂yE45b

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

GR=C0KRFwKV3hc2ZhR20Z<hc2hR2<ZE46
Fw=ECgbohc+hRE47
KV=VV0=KV1h0hE48
KV1=V1V0.E49

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 V0, V/V0, where V is the velocity of the ebb-tidal currents and V0 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 hc and hR, as in Eq. (47). V1 is the flow velocity of the ebb-tidal currents on the plane bottom. KV1 is the ratio of the velocity to the flow velocity at the inlet on the plane bottom. In addition, h0 and h in Eq. (48) are the reference depth on the plane bottom and water depth, respectively. hc2 and hR2 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 KV1, 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. KV1 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.

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 q= (qx, qy) is a linear sum of the component due to waves, qw= (qwx, qwy), and that due to currents, qc= (qcx, qcy):

q=qw+qchcZhRE50
qw=C0Pwtanβctanβcewb1ZE51
Pw=Kwρum3E52
qc=C0Pctanϕtanϕecb2ZE53
Pc=Kcρum2VE54

Here, qis the total sand transport, qwis the sand transport due to waves, qcis the sand transport due to currents, ecis the unit vector in the current direction, and 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. b1=cosα1=ewZ/Z, α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, b2=cosα2=ecZ/Z, α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, b2 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 um2V, 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 hR and the shoreline. In this case, the imaginary depth h’ can be obtained from Eq. (18) using the minimum depth h0 and hR. In particular, in the calculation in Chapter 4, h0 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 hc 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

hc=1+aV/umhc.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 (xi,yj)=iΔxjΔyis taken at the center of the cell, as shown in Figure 11, and Zij=Z(xi,yj,t)(1 ≤ 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, qxi1/2j(1 ≤ i ≤ Nx + 1, 1 ≤ j ≤ Ny) and qyij1/2(1 ≤ iNx, 1 ≤ j ≤ Ny + 1) are set at points i1/2ΔxjΔyand iΔxj1/2Δy, 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 q=(qx,qy)is calculated, the seabed elevation Zij=Z(xi,yj,t+Δt)after time Δt can be calculated using Eq. (56), which is the discretized form of Eq. (29) in Chapter 1:

Z(i,j)=Z(i,j)+qxi1/2jqxi+1/2jΔtΔx+qyij1/2qyij+1/2ΔtΔyE56

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 dqx/dx=0or dqy/dy=0is 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 Hb, αb, hc, and hR, 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

q=A1ew+A2Z+A3tanβes,E57

where ew, Z, and tanβesare 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):

q=C0K1PtanβctanβcewZE58
A1=C0K1PE59a
A2=A1tanβcE59b
A3=0E59c

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

q=C0PtanβcKntanβcewcosαZ+KsKnsinαK2tanβ¯∂H∂stanβeshcZhRE60
A1=C0KnPE61a
A2=cosαtanβcA1E61b
A3=1KntanβcKsKnsinαK2tanβ¯∂H∂sA1E61c

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:

qx=A1cosθw+A2∂Z/∂xA3∂Z/∂yE62a
qy=A1sinθw+A2∂Z/∂y+A3∂Z/∂xE62b

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

qxi1/2jis calculated as

qxi1/2j=A1i1/2jcosθwi1/2j+A2i1/2j∂Z/∂xi1/2jA3i1/2j∂Z/∂yi1/2j,E63

where

∂Z/∂xi1/2j=ZijZi1jΔx,E64
∂Z/∂yi1/2j=12∂Z/∂yi1j+∂Z/∂yij=Zi1j+1Zi1j1+Zij+1Zij14Δy.E65

qxi+1/2jis calculated as

qxi+1/2j=Eq.63i=i+1.E66

qyij1/2is calculated as

qyij1/2=A1ij1/2sinθwij1/2+A2ij1/2∂Z/∂yij1/2+A3ij1/2∂Z/∂xij1/2,E67

where

∂Z/∂yij1/2=ZijZij1Δy,E68
∂Z/∂xij1/2=12∂Z/∂xij1+∂Z/∂xij=Zi+1j1Zi1j1+Zi+1jZi1j4Δx.E69

Finally, qyij+1/2is calculated as

qyij+1/2=Eq.67j=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 hR or hc to prevent sand from being deposited in the zone higher than hR and the beach from being eroded in the zone deeper than hc [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 q=(qx,qy)are reduced by multiplying by a coefficient μ as in Eq. (71), where μ is given by Eq. (72). ΔZ0 is the thickness of the sand layer from which the reduction in sand transport begins, and ΔZ is given by Eq. (73) near hR and by Eq. (74) near hc (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 hR and depth of closure hc.

q=μqE71
μ=μZ=ΔZΔZ00μ1E72
ΔZ=hRZnearZ=hRE73
ΔZ=ZhcnearZ=hcE74

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:

qxi1/2j=μiupqxi1/2jnearZ=hcE75a
μiup=ΔZiupΔZ0=ZiupjhcΔZ0E75b
iup=i1qxi1/2j0iqxi1/2j<0E75c

Figure 13.

Method selecting Z value near hc depending on the direction of sand transport qx 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 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).

Figure 14.

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

qxi1/2j=μidownqxi1/2jnearZ=hRE76a
μidown=ΔZidownΔZ0=hRZidownjΔZ0E76b
idown=iqxi1/2j0i1qxi1/2j<0E76c

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 ΔZ0, 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, ΔZ0 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]:

ΔZ0=12tanβcΔLΔL=minΔxΔyE77

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

ΔZ0=0.050.2hc.E78

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

Regarding the seaward sand movement in regions above hR and deeper than hc, 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 qgdue 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.

qg=AgtanβtanϕentanβtanϕE79

Here, enis 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β=Zis the local slope, and tanϕis the slope of the angle of repose of sand particle. Here, we assume tanϕ=1/2. Sand transport with flux qgtakes place when the local slope exceeds the slope of the angle of repose of sand, the direction of qgis 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 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, qgis set to 0.

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

Z=tanβen=Z/nE80
qg=Agtanϕ/tanβ1ZtanβtanϕE81

Eq. (81) can be expressed by the basic expression for sand transport of Eq. (57), as mentioned earlier, and the coefficients A1, A2, and A3 of each term of Eq. (57) are given by

A1=A3=0E82a
A2=Agtanϕ/tanβ1.E82b

Accordingly, qgcan be evaluated in the same manner as Eq. (57) using the coefficients given by Eq. (82), assuming that the coefficient Ag 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, qgin Eq. (79) is reduced to Eq. (83) using Eq. (80), and the cross-shore sand transport component 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):

qg=AgtanϕZ/nenE83
qn=enqg=AgtanϕZ/nE84
∂Z∂t=qn∂nE85
∂Z∂t=AgZ2n2E86
Δt0.5Δn2AgE87

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:

Ag0.5Δn2ΔtE88

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

Ag=RSΔL2ΔtE89a
RS=0.2E89b
ΔL=minΔxΔy.E89c

Finally, qgis 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):

q=qgifqg>qdueto wavesqdueto wavesotherwiseE90

© 2018 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution-NonCommercial 4.0 License, which permits use, distribution and reproduction for non-commercial purposes, provided the original is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Takaaki Uda, Masumi Serizawa and Shiho Miyahara (December 19th 2018). Eight Types of BG Models and Discretization, Morphodynamic Model for Predicting Beach Changes Based on Bagnold's Concept and Its Applications, Takaaki Uda, Masumi Serizawa and Shiho Miyahara, IntechOpen, DOI: 10.5772/intechopen.81412. Available from:

chapter statistics

225total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Morphodynamic Model for Predicting Beach Changes Based on Bagnold's Concept and Its Applications

Authored by Takaaki Uda

Next chapter

Prediction of Typical Beach Changes Owing to Human Activities

By Takaaki Uda, Masumi Serizawa and Shiho Miyahara

Related Book

First chapter

Organisational and Individual Leadership Challenges During the Rescue Operation Following the 2004 Tsunami

By Aida Alvinius, Erna Danielsson and Gerry Larsson

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us