Intra-wave-phase cross-shore profile modelling by using boundary-fitted slowly moving grid

2.1. Flow module

Fluid flow is basically three-dimensional due to three-dimensionality of turbulence. However, two-dimensional simulation might be acceptable to represent a flow field, if the flow is uniform in one redundant direction like flow within a uniform sectional flume. The flow module simulates free-surface gravity wave motion of incompressible fluid in $x^{\text{'}}$-$z$ domain.

The governing equations of the flow module are the continuity equation and two momentum equations in two directions. The continuity equation for incompressible fluid reads:

where $u$ is the fluid velocity in the $x\text{'}$ direction, $w$ is the fluid velocity in the $z$ direction, $x\text{'}$ is the longitudinal horizontal direction, and $z$ is the upward vertical direction. The Reynolds-averaged Navier-Stokes equations, and the eddy viscosity are introduced as:

where $\rho$ is the fluid density; $p$ is the pressure; ${\epsilon }_{x^{\text{'}}}$, and ${\epsilon }_z$ are the turbulent eddy viscosity in the $x^{\text{'}}$ and $z$ directions, respectively; and $g$ is the gravitational acceleration.

The position of free water surface is computed by using the volume-of-fluid (VOF) method.

The above governing equations have often been solved with a fixed, rectangular grid, see Fig. 3(a). For simulation of morphological change, fixed grid cannot follow micro-scale bed level change with a finite-sized grid height. It is useful for computation of flow around fixed bed morphology only.

Another group tried to introduce partially full grid concept at the liquid-soil interface, see Fig. 3(b)-(c). Fig. 3(b) allows partially full grid with flat or vertical border, and Fig. 3(c) allows partially full grid with a linearly sloped border [18]. These trials may produce continuous flow pattern around the partially full grid, but they have not yet been proved to produce smooth flow-related physical properties like bed shear stress, which acts an important role for a computation of sediment transport and morphologic change, an existing bed-boundary-fitted grid, sigma-coordinate fails for short period waves, especially for waves with overturning crests, which is not good for present purpose.

In order to simulate sediment transport at coastal zones a bed-boundary-fitted grid is adopted here, see Fig. 3(d). A typical coastal beach profile may be expressed as Fig. 4 when this new grid is adopted.

Fig. 3Merits and demerits of grid systems for simulation of morphologic change

a) Bed level rises by a grid size

b) Partial grid with flat level

c) Partial grid with slope

d) Grid moves as accretion or erosion height

The bottom line of the new grid follows the slowly-moving bed profile, and the other higher grid points move in parallel with the bottom grid line. The vertical grid border lines are in the $z$ direction. The bottom grid line is folded with a small angle every delta $x$. It is assumed that the grid moving speed is not large, and the folding angle at every horizontal grid point remains small, either. Coastal geography is not flat, but sloped or wavy. In most cases coastal bed slope does not exceed a certain value, unless non-erodible rocks are around. The steepest slope is often observed at swash zone along cross-shore profiles, and is smaller than the angle of repose for still water condition. The folding angle at every horizontal grid point interval is limited to 10° (slope = 0.18). To simulate wave transformation through curvy grid lines with no discontinuity in pressure distribution, gravity is modified: additional gravity, $g_x$, in the $x\text{'}$ direction is added in the momentum equation along grid columns, see Fig. 5.

The speed of morphological change due to deposition or erosion is not fast, and therefore the whole grid adjusts to slowly varying bed profile every time step. Then, the above governing equations become:

The turbulent eddy viscosity is modelled by a mixing length hypothesis, which has been known to be useful for flows near plane wall boundary [18]. Kim et al. [19] adopted the mixing length hypothesis with a substitution of the vertical gradient of horizontal fluid velocity ($\partial u/\partial z$) with fluid deformation ($\partial u/\partial z+\partial w/\partial x$) for flows over wavy bed form. The mixing length can be approximated by the vertical distance from the point of interest to the nearest bed surface, if bed slope is not steep. The eddy viscosity is assumed isotropic. Then:

where $l$ is the mixing length, $\kappa$ is the von Karman constant (= 0.4), $z^{\mathrm{*}}$ is the distance between the internal point and bed surface, and $h$ is the water depth.

Fig. 4New grid system adapted for present model system: whole grid levels change in the vertical direction depending on bed level change

a) New boundary-fitted parallelogram grid and axis $x$. Air-water interface is expressed by two kinds of grids, $F\ge$0.5 and $F<$0.5, where $F$ is volume-fill ratio

b) Distorted viewing of still water level with axis $x$

Fig. 5Artificial gravitational acceleration induced in the momentum equation

Based on operator splitting theory [20], each of the above two momentum equations are split into two equations:

where superscripts *, ** express the first and second split operators, respectively. Eqs. (6) and (8) are solved by using the first order upwind scheme, and Eqs. (7) and (9) are solved by using a time-forward, space-centered scheme.

The pressure is then computed by using the SOLA method [21] with coincident adjustment of velocity field. The present model uses a fixed coefficient $\alpha$ (= 0.8) for pressure-divergence adjustment instead of using over relaxation parameter of the SOLA:

Newton’s second law of motion is used for pressure adjustment instead of fluid compressibility for compressible fluid flows:

where $D$ is the divergence $(\partial u/\partial x+\partial w/\partial z)$, subscripts $i$, $k$ are indices for variables in grid points in a staggered grid in the $x$ and $z$ directions, respectively, see Fig. 6.

Fig. 6Staggered grid for the present flow model

2.2. Sediment transport module

Sediment transport is computed including the suspended sediment load, and the bed load. The suspended sediment transport describes sediment movement in suspended mode subject to flow conditions. On the other hand, the bed load is defined as instantaneous horizontal response to imposed shear stress on the bed surface.

The suspended sediment is often expressed as concentration, and governed by sediment mass conservation, which reads:

where $s$ is the suspended sediment concentration [kg/m³]; $w_f$ is the settling velocity of sediment particles; ${\epsilon }_{s,x}$, ${\epsilon }_{s,z}$ are the diffusion coefficients of the suspended sediment concentration in the $x$, $z$ directions, respectively.

There is an interface between water column containing the suspended sediment and the deposited bed, where deposition moves sediment from the water column onto the bed, and pickup moves sediment from the bed into the water column. The deposition rate is computed from settling of suspended sediment as:

where $D$ is the deposition rate [kg/m²s], and $C_h$ is the coefficient to take into account hindered settling due to high concentration near the bed (unity for low concentration), and $R$ is the probability function.

The bed load is related to both flow and bed sediment properties including bed shear stress. The bed shear stress acting on the bed surface exerted by flow can be approximated by fluid particle near-bed velocity under assumption of constant shear stress profile or equilibrium logarithmic velocity profile within thin boundary layer, that is:

where ${\tau }_b$ is the bed shear stress, $C_f$ is the instantaneous current friction coefficient (=${\kappa }^2\mathrm{l}\mathrm{n}(z/z_0{)}^{-2}$), $\kappa$ is the von Karman constant (= 0.4), $z_0$ is the one thirtieth of bed roughness, and $u\left(z\right)$ is the fluid horizontal particle velocity at level z from the mean bed level.

As we look in inter-wave-phase sediment movement, we need a prediction formula for instantaneous bed load.

Some existing empirical formulas for bed load are for wave-phase-average bed load at wave-current environment, while other existing formulas are for steady current. We need not distinguish empirical bed load formula for steady current or for instantaneous current, if our approach resolves the wave phase. Meyer-Peter and Muller [22], Bagnold [23], and van Rijn [24] proposed empirical formulas of bed load for current, all of which have a similar form as:

where $q_b$ is the bed load [kg/ms], $C_b$ is the scaling coefficient, ${\tau }_b$ is the bed shear stress, ${\tau }_{cri}$ is the critical shear stress needed for initiation of motion of a specific or representative grain size, and n is the empirical power. Non-dimensional form is practically used for arbitrary bed grain size, and more universal coefficient:

where $q_b^{\mathrm{\text{'}}}$ is the non-dimensional bed load; $C_{b,2}$ is the scaling coefficient; $T_b\text{,}$$T_{cri}$ are non-dimensional bed shear stress, and critical shear stress, respectively. Van Rijn’s [24] bed load formula is adopted here:

where ${\rho }_s$ is the sediment density, $\rho$ is the fluid density, $d_{50}$ is the median diameter [m], $g$ is the gravitational acceleration, $\nu$ is the fluid kinetic viscosity, ${\tau }_b$ is the bed shear stress, and ${\tau }_{cri}$ is the critical shear stress for initiation of motion.

Bed load with instantaneous response to the external force is believed to be affected by the local bed slope. Some researchers suggested modification of the critical motion-initiating bed shear stress due to the local bed slope [25-27]. Dey analyzed their laboratory data, and proposed the following formula:

where ${\theta }_x$, ${\theta }_y$ are bed slopes in the $x$, $y$ directions, respectively; and $\mathrm{Ф}$ is the angle of repose. For flows in $(x-z)$ two-dimensional vertical domain, the second part of Eq. (23) disappears.

For the present study we look in the problem from mechanical point of view instead of analysis of experimental data. Assume a situation that a cylinder is laid on a flat bed. Two moments balance i.e. overturning moment due to flow shear force, and resisting moment due to gravitational force, see Fig. 6. There is a pivot on a cylinder surface at the angle of repose on the vertical line from cylinder centre. Now we consider a case when flow is upward over a sloped bed. The resisting moment due to gravity is increased because of the expanded distance of action of the force to the vertical line from the cylinder centre. Thus, the critical overturning moment becomes larger. On the other hand for downward flow the resisting moment becomes smaller. This modification of resisting moment can be expressed as:

where $\theta$ is the upward bed slope. If bed slope is downward angle of repose (= $-\mathrm{Ф}$), the critical shear stress becomes zero as the above equation. Therefore, we modify Eq. (17) by replacing ${\tau }_{cri}$ with ${\tau }_{cri,\theta \mathrm{}}$ by using the relationship between the two in Eq. (24).

Fig. 7Effect of bed slope on critical bed shear stress for initiation of motion

The pickup rate has similarity to the bed load from the point that the initiation of motion is suppressed, if the bed shear stress is below the critical value, ${\tau }_{cri}$, and is more or less proportional to the bed load. Partheniades [28], Nielsen [29], and van Rijn [30] proposed formulas for instantaneous pickup rate in the following form as:

However, the pickup is vertical movement, and therefore empirical coefficient and unit are different from those of bed load. We use the following modified van Rijn’s formula for pickup rate by adding an additional coefficient, $C_{p,2}$ as:

where $C_r$ is the reference level concentration [volume/volume], and $a$ is the reference level above mean bed.

2.3. Bathymetric change module

Bed level change is computed from the total sediment mass conservation with sediment mass exchanges to and from water column:

where $z_b$ is the bed level [m], $t$ is time [s], $e$ is the void ratio of bed sediment, ${\rho }_s$ is the pure sediment density [kg/ m³], $D$ is the deposition rate [kg/m²s], $P$ is the pickup rate [kg/m²s], and $q_b$ is the bed load [kg/ms].

Bathymetric change is computed periodically, and extrapolated with a multiplication factor for simulation of long term change. When the bed slope at a part of the bed profile exceeds the underwater angle of repose, the bed profile is modified by using a moving boundary filter. Interim profile is further treated by using the moving-average filter multiple times as needed.

3.1. Computation results of Case 1-3

Measured new bed profile of Case 1-3 at 69.5 hours from the start is shown in Fig. 8. Computed water level and pressure distribution, and the particle velocity vector field at a wave phase after 30 waves are shown in Figs. 9-10. Computed results are presented in $x$-$z$ domain for convenience. Acceleration parameter of 1500 was applied to the bed level change speed.

Fig. 9Measured cross-shore profile change after 69.5 hours of wave action for Case 1-3 of Kajima et al. [31]

Fig. 10Computed pressure after 30 waves after 30 waves

In order to reduce computational load, the bed level change speed was accelerated by using a scaling factor (multiplication factor). The laboratory experiments lasted 69.5 hours, but only 30 waves were run for numerical modelling. After every 15 waves bed profile change was allowed with the acceleration factor.

Generated waves on flat bed propagate over sloped profile, and get broken as spilling type. Fluid particle velocity near the seabed is affected by the bed friction.

Computed phase-average horizontal velocity component (residual current) profile at $x=$80 m after 30 waves is shown in Fig. 11. The residual current near water surface is more onshore-ward, while the residual current near bed is more offshore-ward.

Asymmetry of temporal near-bed particle velocity variation with respect to mean sea level has been known as one of the important factors affecting the direction of net sediment transport on cross-shore profile. Previous research results described typical pattern as strong onshore current for a short term, and weak offshore current for a long term [29], in other words asymmetric speed with respect to the mean water level drives on-shore sediment transport. Computed near-bed velocity variation at a point shows this asymmetry, see Fig. 12. The computed near-bed velocity asymmetry at a point within a swash zone within a wave period is considered to include the effect of enhanced turbulence due to wave breaking.

Fig. 11Computed particle velocity vector field after 30 waves

Fig. 12Computed residual horizontal velocity profile at x'= 80 m after 30 waves

Fig. 13Computed near-bed particle velocity variation within a wave period after 30 waves

a)$x=$86.4 m

b)$x=$40 m

Computed wave-phase-average suspended sediment concentration distribution in the domain after 30 waves is shown in Fig. 13. The suspended sediment concentration is high near eroded zone, and it settles down quickly around the eroded zone.

Computed bed profile with measured one after 69.5 hours are shown in Fig. 14. The trend of accretion at onshore zone was well reproduced. The computed eroded position and the maximum erosion depth of 0.55 m were reasonably reproduced compared to measured erosion depth of about 0.77 m. The computed accretion position and the maximum accretion height of 0.89 m were also well reproduced, compared to measured accretion of about 1.02 m. It is thought that the model system satisfactorily reproduced the sediment transport and following bed profile change for this case.

Fig. 14Computed wave-phase-average suspended sediment concentration field after 30 waves

3.2. Computation results of Case 3-4

Case 3-4 is for a setting with higher wave height, shorter wave period, and finer grain size. Measured bed profile change of Case 3-4 at 76.1 hours is shown in Fig. 15. Waves break earlier than Case 1-3, and the asymmetry of the orbital velocity is not prominent, and therefore doesn’t much enhance the onshore sediment transport in this case.

The computed water level and pressure distribution, and the particle velocity vector field at a wave phase at the beginning stage and after 30 waves are shown in Figs. 16-17. Computed results are presented in $x\text{'}$-$z$ domain for convenience. The acceleration parameter of 4750 was applied to the bed level change speed, respectively.

Fig. 15Computed bathymetric change after 69.5 hours waves for Case 1-3

Fig. 16Measured cross-shore profile change after 76.1 hours of wave action for Case 3-4 of Kajima et al. [31]

Fig. 17Computed pressure after 30 waves

Computed phase-average velocity (residual current) vector profile at $x^{\mathrm{\text{'}}}=$80 m is shown in Fig. 18. The residual current near water surface is more onshore-ward, while that near bed is more offshore-ward.

Near-bed velocity variation within a wave period for Case 3-4 is also asymmetric with respect to the mean water level, but the asymmetry is weaker then Case 1-3, see Fig. 19.

Computed wave-phase-average suspended sediment concentration distribution in the domain after 30 waves is shown in Fig. 20. The suspended sediment concentration is high near eroded zone, and it settles down around the eroded zone.

Fig. 18Computed particle velocity vector field after 30 waves

Fig. 19Computed residual horizontal velocity profile at x= 80 m after 30 waves

Fig. 20Computed near-bed particle velocity variation within a wave period after 30 waves

a)$x=$86.4 m

b)$x=$40 m

Fig. 21Computed wave-phase-average suspended sediment concentration field after 30 waves

Computed bed profile with measured one after 76.1 hours is shown in Fig. 21. The trend of eroded onshore zone and offshore bar formation were well reproduced. The computed eroded position and the maximum erosion depth of 0.99 m were reasonably reproduced compared to measured erosion depth of about 0.81 m. The computed accretion position and the maximum accretion height of 1.82 m were also well reproduced, compared to measured accretion of about 2.04 m. It is thought that the model system satisfactorily reproduced the sediment transport and following bed profile change of Case 3-4, as well.

Fig. 22Computed bathymetric change after 76.1 hours for Case 3-4