Sediment Transport Equivalent Waves for Estimating Annually Averaged Sedimentation and Erosion Trends in Sandy Coastal Areas

Abstract

In this paper, a simple approach to determine representative offshore wave characteristics for estimating the annually averaged sedimentation and erosion trends in sandy coastal areas is presented. Given the offshore wave climate, the proposed approach breaks down the climate into fixed 22.5-degree bins and based on the sediment transport potential it determines the equivalent wave characteristics for each bin, i.e., a significant wave height, a peak period, a mean wave direction, and a corresponding frequency of occurrence. The approach is validated in idealized cases of uniformly sloping beaches with the presence of a breakwater, for various sediment diameters, sea bottom slopes, and different offshore wave characteristics. The performance of the proposed approach is evaluated against the full climate, returning good results. Furthermore, the proposed approach is applied in a real-life challenge, in the coastal area of Therma in the Island of Samothraki in Greece, where the presence of a fishing shelter has led to sedimentation and erosion problems. The performance of the proposed approach is very satisfactory, given the complexity of the problem. The generic nature of the proposed methodological approach allows it to be applied in numerous sandy coastal regions to estimate the sedimentation and erosion trends, reducing the amount of input parameters and thus requiring significantly less computational efforts.

1. Introduction

With climate change unfolding and human interventions in coastal zones increasing, erosion is one of the greater risks to coastal communities. Therefore, over the past decades, coastal engineering scientists have significantly contributed to the development of powerful numerical models for predicting the evolution of coastal morphodynamics and empowered engineers to test and optimize schemes of coastal protection works. These coastal area models are applicable at a range of spatial scales–from small to macro scale coastal engineering problems. Despite the significant progress that has been made on 3D- and quasi 3D models [1,2,3], prohibitive computational efforts are required when simulating real coastal areas of several kilometers, and this is why 2DH Models [4,5,6] are widely used in real-life coastal engineering challenges. However, even these models usually require small grid cell sizes and small time-steps to resolve the dominant processes and thus they remain computationally intensive, especially when the prediction of morphological changes is required over a longer time than a storm duration, e.g., on an annual basis. Therefore, to investigate the coastal morphodynamics at a time horizon of one year, many engineers and researchers resort to the concept of Initial Sedimentation and Erosion (ISE) modelling, as named [7]. This concept does not rely on full morphodynamic simulations (i.e., the bed level is not updated after each morphological time-step) but goes only once through the implemented numerical models, i.e., the wave, hydrodynamic, sediment transport, and morphological models. All simulations are based on the same initial bathymetry and the rates of bed level change are calculated for each considered incident wave scenario. Subsequently, each resulting rate of bed level change field can be multiplied by the annual frequency of occurrence of the corresponding wave scenario, and the sum of these fields yields the annually averaged rates of bed level changes, providing a reliable estimate of sedimentation and erosion trends in a coastal study area. ISE modelling is useful in comparative investigations, e.g., for predicting the behavior of alternative coastal protection works or for testing the contributions of various individual sea state conditions [7].

Nevertheless, despite the computation burden relief that ISE modelling can provide, most of the coastal engineering projects require numerous simulations to thoroughly investigate the processes which take place in a coastal zone, followed by additional simulations to test and conclude with the most sustainable and efficient alternative of coastal protection works (e.g., system of groins, breakwaters, artificial nourishment, etc.). Moreover, additional simulations are required to optimize the key dimensions of the proposed works (e.g., the length and orientation of breakwaters/groins, the gap between them, the distance from the shoreline, etc.) and on top of that, climate change scenarios (with various sea water levels and varying incident wave characteristics) should also be considered to test the stability and operability of the proposed structures during their design lifetime.

To this end, offshore wave climate schematization (alternatively called wave input reduction) is required to reduce the input wave data to force the process-based models. The most commonly implemented wave schematization methods [8,9,10] divide the multivariate wave climate in non-equidistant directional and wave height bins. The boundaries of each bin are defined based on the principle that each bin contains an equal portion of a proxy quantity considered vital in influencing morphological bed evolution in the medium term, such as the wave energy flux or bulk longshore sediment transport rates. Then, a representative sea-state is calculated for each bin either by taking the mean value of the wave characteristics in each class [11,12] or considering a non-linear relation to select the representative significant wave height [8,10]. The performance of the abovementioned “binning” wave schematization methods depends on a multitude of parameters such as the number of representatives to be selected (usually a number between 10–12 sea-states at minimum should be selected to represent annual morphological bed evolution [8,9,10]), the duration of the wave climate to be reduced, and the subdivisions in wave height and directional bins. It should be noted that simplistic computer codes must be developed and applied to automatically define the directional bins when in possession of a large array of offshore sea-state wave characteristics. A disadvantage of the binning wave schematization methods is that they often overly estimate the contribution of lowly-energetic sea-states that are present in the wave climate dataset, hence they are usually combined with arbitrary elimination of sea-states below a certain threshold [9]. However, this approach depends mainly on the coastal modeler’s discretion, and adopting a high threshold can inadvertently lead to the deterioration of results [8]. The authors in [12] introduced a binning wave schematization method which incorporates a systematic filtration of lowly energetic sea-states based on the well-known criterion of incipient sediment motion [13] and utilizing the sediment pick-up rate [14,15] as a proxy quantity to define the wave representatives. Although this method further accomplished the acceleration of morphological modelling simulations, it is rather complex since it requires additional wave transformation simulations to estimate nearshore wave characteristics. As a general principle, wave schematization methods should both be based on solid theoretical principles on the main driving factors influencing medium-term morphological bed evolution and be relatively straightforward and easy to apply.

The present work adopts the ISE modelling concept and aims to contribute to the offshore wave climate input reduction methods by proposing an approach to determine sediment transport equivalent waves for estimating the sedimentation and erosion trends (i.e., the rates of bed level changes) in a coastal zone due to the combined action of waves and currents, with significantly less computational efforts. The approach to be presented herein considers directional bins of fixed size on which the calculation of the sediment transport equivalent waves will subsequently be carried out. Consequently, it is even less computationally intensive and simple than the classical binning input reduction methods, and calculations can readily be carried out in a single spreadsheet while the obtained results are deemed reliable and satisfactory.

The paper structure has as follows: Section 2 describes the “classical” approach to schematize a dataset of offshore sea-state characteristics whereas Section 3 describes the proposed methodology to estimate annual equivalent wave representative conditions; Section 4 gives a brief background of the implemented numerical models and performance evaluation metrics; Section 5 presents the application and results of the proposed approach in both an idealized and real-life case study; and ultimately Section 6 discusses results and draws conclusions.

2. Classical Approach for Determining the Offshore Wave Characteristics

The Classical Approach (CA hereafter) exploits hindcast wave climate data from open databases (e.g., [16,17]), which provide hourly wave records of the offshore wave characteristics, i.e., the significant wave height, $H_{soi}$ (m), the peak period, $T_{pi}$ (s), and the mean wave direction, $MW{D}_i$ (°, clockwise from North).

In order to quantify the offshore wave climate, the wave characteristics are usually, in the context of a CA, grouped by wave height and mean wave direction, in particular said offshore wave climate is divided into equally spaced groups of 0.5 m or 1.0 m, and the mean wave directions is distributed into 45.0°, 30.0°, or 22.5° bins (i.e., 8, 12, or 16 sectors, respectively). The frequency of occurrence, $f_g$ (%), per each wave height and direction group can be then calculated as the fraction of the number of records, $N$, with characteristics that lie within the boundaries of each group, to the total number of records, $N$, of the obtained time series.

For instance, hindcast simulations wave data, which will be used for the purposes of the present paper, were obtained from the data package entitled MEDSEA_MULTIYEAR_WAV [18] of the Copernicus Database, for a period spanning from 1993 to 2019. The wave data extraction point is located at 40.469° N, 19.3940° E. The obtained parameters include the significant wave height, the peak wave period, and the mean wave direction, on an hourly basis. The wave heights were classified into equally spaced groups of wave heights with a step of 0.50 m, and the $MWD$ was broken down into 22.5° bins, i.e., 16 sectors (N, NNE, NE, ENE, E, ESE, SE, SSE, S, SSW, SW, WSW, W, WNW, NW, and NNW). The frequencies of occurrence, $f_g$ (%), per each wave height group and sector are given in

Table 1. Mean annual frequencies of occurrence per wave height group and sector (the deeper blue colors denote the classes with the highest frequency of occurrence).

$\mathrm{Group}{H}_s\left(m\right)$	S	SSW	SW	WSW	W	WNW	NW	NNW	N	Other	Per Group
[0–0.5]	0.576%	7.784%	14.705%	4.461%	4.495%	16.002%	19.015%	0.377%	0.036%	0.121%	67.572%
(0.5–1]	0.461%	4.722%	4.396%	1.096%	1.193%	9.247%	5.178%	0.115%	0.000%	0.006%	26.415%
(1-1–5]	0.040%	0.752%	1.122%	0.200%	0.199%	2.002%	0.687%	0.003%	0.000%	0.000%	5.005%
(1.5–2]	0.000%	0.045%	0.211%	0.023%	0.037%	0.438%	0.053%	0.000%	0.000%	0.000%	0.807%
(2–2.5]	0.000%	0.000%	0.0169%	0.003%	0.006%	0.130%	0.002%	0.000%	0.000%	0.000%	0.158%
(2.5–3]	0.000%	0.000%	0.005%	0.000%	0.000%	0.031%	0.000%	0.000%	0.000%	0.000%	0.036%
(3–3.5]	0.000%	0.000%	0.000%	0.000%	0.000%	0.006%	0.000%	0.000%	0.000%	0.000%	0.006%
(3.5–4]	0.000%	0.000%	0.000%	0.000%	0.000%	0.001%	0.000%	0.000%	0.000%	0.000%	0.001%
Per sector	1.077%	13.303%	20.457%	5.783%	5.930%	27.857%	24.935%	0.495%	0.036%	0.127%	100.000%

and a rose diagram is illustrated in Figure 1. Sectors NNE, NE, ENE, E, ESE, SE, and SSE are omitted from Table 1 (since they will not be taken into account for the purposes of this work) and are banded together under the “Other” label.

Adopting the CA, numerous wave, hydrodynamic, and sediment transport simulations must be carried out, for each of the incident wave scenarios based on the offshore wave climate, in order to assess the annually-averaged rates of bed level changes in a coastal field. Subsequently, the results of rate of bed level changes are integrated by assigning the frequencies of occurrence, $f_g$, as weights to the different incident wave scenarios for which simulations have been performed. Given the specific offshore wave climate (Table 1) and by assuming a coastal area, with a hypothetical shoreline orientation from North to South, implying that the area is exposed to the seven sectors–SSW, SW, WSW, W, WNW, NW, and NNW–36 incident wave scenarios should be carried out in this example to conclude with the annually-averaged rates of bed level changes. An illustrative flow chart of the CA is given in Figure 2 for this case. It is noted that the significant wave height and the peak period, assigned for each scenario to be simulated, have been taken as the respective average of each wave height group and sector resulting from the obtained dataset.

It becomes apparent that a substantial number of simulations are required. Τhe computational burden increases more when several coastal protection alternatives are required to be investigated, as is usually the case in coastal engineering projects in order to find the most efficient one, and consequently the aforementioned approach should be repeated and equal amount of times for each of the considered alternatives. Furthermore, additional scenarios accounting for climate change conditions and extreme storm events should be considered as well, further increasing the required times.

3. Proposed Approach for Determining the Sediment Transport Equivalent Waves

The Proposed Approach (PA hereafter) aims to determine a single representative wave for each sector of interest based on the sediment transport potential. The quantities that need to be determined for a representative wave is the equivalent wave height $H_e$ (m), the equivalent period $T_{pe}$ (s), the equivalent frequency $f_e$ (%), and the equivalent mean wave direction $MW{D}_e$ (°), per each sector.

Many researchers have endeavored to provide sediment transport equations (e.g., [19,20,21,22,23]) based on dimensional analysis, force-balance, or energetic methods. These formulas are widely used in engineering challenges and in research efforts. The formula of [19] as modified and improved by [24] is used herein and reads:

$$Q=\frac{0.149}{\left({\rho }_s-\rho \right)\left(1-p\right)}{H}_{sb}^{2.75}{T}_p^{0.89}{m}_b^{0.86}{d}_{50}^{-0.69}{sin}^{0.5}({2a}_b)$$

(1)

where $Q$ (m³/s) is the longshore sediment transport rate, ${\rho }_s$ is the sediment density, $\rho$ is the ambient water density, $p$ is the porosity, $H_{sb}$ (m) is the significant wave height at the breaker line, $T_p$ (s) is the peak wave period, $m_b$ is the breaking zone beach slope, $d_{50}$ (m) is the median grain size, and $a_b$ (degrees) is the angle between wave propagation direction at the breaker line and shore normal direction. The breaking wave characteristics can be calculated following the methodology presented in [25].

The equivalent wave height, $H_e$, for each sector can be determined as a weighted average with respect to the frequency of occurrence and the sediment transport of each wave record, as follows

$$H_e=\frac{{{\displaystyle \sum }}_{i=1}^{N_{Dir}}\left(w_i{H}_{soi}\right)}{{{\displaystyle \sum }}_{i=1}^{N_{Dir}}\left(w_i\right)}$$

(2)

where $N_{Dir}$ is the number of records appearing in the dataset and that lie within the sector. The weights $w_i$ represent the sediment transport potential of each offshore wave and are given by

$$w_i{=f}_i{Q}_i{=f}_i\frac{0.149}{\left({\rho }_s-\rho \right)\left(1-p\right)}{H}_{sbi}^{2.75}{T}_{pi}^{0.89}{m}_{bi}^{0.86}{d}_{50i}^{-0.69}{sin}^{0.5}({2a}_{bi})$$

(3)

Assuming that the porosity and sediment characteristics, ${\rho }_s$ and $d_{50}$, remain constant along the coastal area and that the sea bottom slope is approximately uniform within the breaker zone, Equation (2) becomes

$$H_e=\frac{{{\displaystyle \sum }}_{i=1}^{N_{Dir}}\left(f_i{H}_{sbi}^{2.75}{T}_{pi}^{0.89}{sin}^{0.5}{(2a}_{bi}{)H}_{soi}\right)}{{{\displaystyle \sum }}_{i=1}^{N_{Dir}}\left(f_i{H}_{sbi}^{2.75}{T}_{pi}^{0.89}{sin}^{0.5}({2a}_{bi})\right)}$$

(4)

Given that the frequency of occurrence, $f_i$, is the same for each record appearing in the dataset (equal to $1/N$), it can be omitted. Otherwise, if the entire dataset is not available and the only known offshore wave climate data are given in a tabular form (akin to Table 1), the frequency of occurrence of each group, $f_g$, should be used in Equation (4) along with its respective wave characteristics, $H_{sbg}$, $H_{sog}$, $T_{pg}$ and $a_{bg}$. The same holds for the following equations.

Similarly, the equivalent peak period, $T_{pe}$, for each sector can be determined as a weighted average with respect to the frequency of occurrence and the sediment transport of each wave record, as follows

$$T_{pe}=\frac{{{\displaystyle \sum }}_{i=1}^{N_{Dir}}\left(w_i{T}_{pi}\right)}{{{\displaystyle \sum }}_{i=1}^{N_{Dir}}\left(w_i\right)}$$

(5)

The equivalent mean wave direction per sector is calculated as the weighted average of the directions which appear in the obtained dataset and lie within each sector. Thus, it is obtained for each sector as follows:

$${MWD}_e=\frac{{{\displaystyle \sum }}_{i=1}^{N_{Dir}}\left(w_i{MWD}_i\right)}{{{\displaystyle \sum }}_{i=1}^{N_{Dir}}\left(w_i\right)}$$

(6)

Finally, the equivalent frequency $f_e$ per each sector can be determined by assuming that the sediment transport caused cumulatively by all the incident waves, within the sector, is equal to the one caused by the equivalent wave, and hence

$$Q_e{f}_e={{\displaystyle \sum }}_{i=1}^{N_{Dir}}\left(f_i{Q}_i\right)$$

(7)

where $Q_e$ is the sediment transport rate which can be found through Equation (1) with the equivalent characteristics (Equations (4)–(6)). Each equivalent frequency of occurrence can be regarded as a quantitative weighting factor in order to conserve the sediment transport potential in each directional sector through Equation (7). It is worth mentioning that the resulting equivalent frequencies of occurrence ($f_e)$of the PA, per each directional sector, are expected to remain smaller than the corresponding sum of the CA for the same directional sector. This behavior is expected to occur in the vast majority of real offshore wave climates, due to the fact that low-energy waves, i.e., waves with small heights and periods, produce low sediment transport rates but occur more frequently, while high-energy waves generate high sediment transport rates but occur less frequently. Given that the equivalent wave is a weighted average with respect to the frequency of occurrence and the sediment transport, it is expected to lie between the low- and high-energy waves. Equations (4)–(7) can be used to determine the sediment transport equivalent waves for estimating the rates of bed level change (i.e., sedimentation and erosion trends) in sandy coastal areas. In Figure 3, the flow chart of the proposed approach is illustrated along with the equivalent characteristics per each sector of interest, representing the considered offshore wave climate of Table 1. It is noted that CA, in the considered example, requires 36 scenarios to be simulated (Figure 2), while the PA requires only 7 (Figure 3), significantly minimizing the required computational efforts.

4. Implemented Numerical Models and Performance Metrics

4.1. Numerical Models

CA and PA rely on the ISE modelling, as mentioned in the introduction. In particular, through this analysis, for each considered incident wave scenario, the spatial distribution of nearshore wave heights and radiation stresses are first simulated with the nonlinear mild-slope wave model of parabolic approximation, given the initial bathymetry. Subsequently, the radiation stresses are given as input to the hydrodynamic model to simulate the wave-driven coastal currents, providing the surface elevation and the current velocity components in the two horizontal dimensions. Finally, the wave height and hydrodynamic results are given as input to the sediment transport and morphological model to ultimately calculate the rates of bed level changes. This procedure is applied for each incident wave scenario, and each resulting field of rates of bed level changes is multiplied by the corresponding annual frequency of occurrence of each wave scenario. These fields are summed to produce the desired annually averaged rate of bed level change field (as illustrated in the flow charts of Figure 2 and Figure 3). Hence, it is noted that the chronology in the selection of representative wave conditions does not play a role. The implemented numerical models responsible for simulating the coastal processes are presented below.

4.1.1. Nearshore Wave Model

The nearshore wave model implemented herein has been developed by [26] based on the work of [27,28], and it is a nonlinear mild-slope wave model of parabolic approximation. The governing parabolic equation is given by:

$$C_g{A}_x+i\left(\overline{k}-a_{0}k\right)C_{g}A+\frac{1}{2}{\left(C_g\right)}_{x}A+\frac{i}{\omega }\left({\alpha }_1-b_1\frac{\overline{k}}{k}\right){\left(C{C}_g{A}_y\right)}_y-\frac{b_1}{\omega k}{\left(C{C}_g{A}_y\right)}_{yx}+\frac{b_1}{\omega }\left(\frac{k_x}{k^2}+\frac{{\left(C_g\right)}_x}{2k{C}_g}\right){\left(C{C}_g{A}_y\right)}_y+\frac{i\omega k^2}{2}D{\left|A\right|}^{2}A+\frac{w}{2}A=0$$

(8)

where $C_g$ [m/s] is the group celerity, $A$ [m] is the complex amplitude related to the water surface displacement given by ${\eta =Ae}^{-i\left(kx-\omega t\right)}$, $k$ [rad/m] is the local wave number related to the angular frequency of the waves, $\omega$ [rad/s], $\overline{k}$ [rad/m] is a reference wave number taken as the average wave number along the y-axis, coefficients $a_0$, ${\alpha }_1$ and $b_1$ depend on the aperture width [29], $C$ [m/s] is the phase celerity, the parameter $D$ is given by $D=\frac{\left({cosh4kh+8-2tanh}^{2}kh\right)}{{8sinh}^{4}kh}$, with $h$ [m] being the water depth and $w$ [s⁻¹] is a dissipation factor due to depth-induced wave breaking [30] and bottom friction [31]. The subscripts $x$ and$y$ denote spatial derivatives in x-and y-direction, respectively.

The model is capable of simulating the propagation of uni-directional irregular waves by discretizing a JONSWAP wave energy spectrum into separate wave components and performing corresponding simulations for each one. The linear superposition principle is then used to obtain the wave characteristics at each cell of the computational domain.

Moreover, an approximate non-linear amplitude dispersion relationship, as in [32], has been introduced in order to improve the model’s results in the nearshore area. The computations are performed in a regular grid and the governing equation is solved through the finite difference method, employing the Crank-Nicholson scheme for parabolic differential equations.

4.1.2. Nearshore Hydrodynamic Model

The hydrodynamic model is based on the incompressible depth-averaged Reynolds Averaged Navier Stokes equations [33], consisting of the continuity and the momentum equations as given below:

$$\frac{\partial \overline{\eta }}{\partial t}+\frac{\partial \left(Uh\right)}{\partial x}+\frac{\partial \left(Vh\right)}{\partial y}=0$$

(9)

$$\frac{\partial U}{\partial t}+U\frac{\partial U}{\partial x}V\frac{\partial U}{\partial y}+g\frac{\partial \overline{\eta }}{\partial x}=-\frac{1}{\rho h}\left(\frac{\partial S_{xx}}{\partial x}+\frac{\partial S_{xy}}{\partial y}\right)+\frac{1}{h}\frac{\partial }{\partial x}\left({\nu }_{h}h\frac{\partial U}{\partial x}\right)+\frac{1}{h}\frac{\partial }{\partial y}\left({\nu }_{h}h\frac{\partial U}{\partial y}\right)$$

(10)

$$\frac{\partial V}{\partial t}+U\frac{\partial V}{\partial x}V\frac{\partial V}{\partial y}+g\frac{\partial \overline{\eta }}{\partial y}=-\frac{1}{\rho h}\left(\frac{\partial S_{yy}}{\partial y}+\frac{\partial S_{xy}}{\partial x}\right)+\frac{1}{h}\frac{\partial }{\partial x}\left({\nu }_{h}h\frac{\partial V}{\partial x}\right)+\frac{1}{h}\frac{\partial }{\partial y}\left({\nu }_{h}h\frac{\partial V}{\partial y}\right)$$

(11)

where: $\overline{\eta }$ [m] is the mean sea surface elevation, $U$ [m/s] and $V$ [m/s] are the depth-averaged current velocities in the x and y axis, respectively, $\rho$ [kg/m³] is the seawater density, $h$ [m] is the total water depth, $g$ [m/s²] is the acceleration of gravity, $v_h$ is the horizontal turbulent eddy viscosity coefficient and $S_{xx}$, $S_{xy}$, $S_{yy}$ are the radiation stress components. The numerical calculations are performed on a staggered numerical grid in which the velocity components are calculated between adjacent sea-level grid points. The solution of the governing equations in the temporal domain is carried out through an explicit Euler scheme.

4.1.3. Sediment Transport and Morphological Model

A non-cohesive sediment transport model is implemented, capable of simulating the sediment transport field and the subsequent bed evolution in coastal areas due to the combined effect of waves and currents. The total load sediment transport is calculated based on [34,35], consisting of the individual contributions of bedload and suspended load transport. The bedload sediment transport rate under the combined action of waves and currents is given by [34]:

$$q_b=0.015\overline{U}h{\left(\frac{d_{50}}{h}\right)}^{1.2}{M}_e^{1.5}$$

(12)

where: ${{\rm M}}_e=\frac{U_e{-U}_{cr}}{\sqrt{{gd}_{50}(s-1)}}$ is the sediment mobility parameter, $U_e=\overline{U}+\gamma ·U_w$ is the effective flow velocity, with $\gamma =0.8$ for monochromatic and $\gamma =0.4$ for irregular waves, and $s=\frac{{\rho }_s}{{\rho }_w}$ the relative grain density.

The critical flow velocity is a function of the individual values due to waves ($U_{cr,w}$) and currents ($U_{cr,c}$), respectively, and is calculated as follows:

$$U_{cr}{=\beta U}_{cr,c}{+(1-\beta )U}_{cr,w}$$

(13)

where: $\beta =\frac{\overline{U}}{\overline{U}{+U}_w}$.

The individual critical velocities are calculated through the relationships given by [36] for waves and [37] for currents.

The suspended sediment transport rate is calculated as shown in [35]:

$$q_s=0.012\overline{U}{d}_{50}{\left(D_{*}\right)}^{-0.6}{M}_e^{2.4}$$

(14)

Finally, the governing equation solved by the model is the sediment mass balance or “Exner” equation, which calculates the rate of bed level change at each time-step:

$$\left(1-p\right)\frac{\partial z_b}{\partial t}+\nabla ·\overrightarrow{q_t}=0$$

(15)

The calculations are carried out in a regular grid and the rates of bed level change are obtained by discretizing the Equation (15) through an Upwind scheme.

4.2. Performance Metrics

In order to evaluate the performance of the PA, the following four measures are considered herein for the variable of annually averaged rate of bed level change $\frac{\partial z_b}{\partial t}$. First, the linear correlation coefficient is used to check the relation between the CA and the predictions of the PA. The correlation coefficient is given by:

$$R=\frac{s_{XY}}{{\sigma }_X{\sigma }_Y}$$

(16)

where $s_{XY}$ is the covariance between CA and PA and ${\sigma }_X{,\sigma }_Y$ are the respective standard deviations.

The Mean Squared Error [(m/day)²] is taken into account for measuring the accuracy, defined as:

$$MSE=\stackrel{╌╌╌╌╌╌╌╌╌╌╌╌}{{\left(Y-X\right)}^2}$$

(17)

where $Y$ is the PA’s results and $X$ is the CA’s results; bar denotes the mean.

Moreover, the Brier Skill Score [38] is calculated through the following relationship:

$$BSS=1-\frac{MSE(Y,X)}{MSE(B,X)}$$

(18)

where $Y$ denotes the PA’s results, $X$ denotes CA’s results, and $B$ is a baseline prediction where zero bed level change rates are considered. Following [38], the classification scores for the BSS to evaluate the performance of a specific morphological evolution model are: 1.0–0.5 Excellent; 0.5–0.2 Good, 0.2–0.1 Reasonable/Fair, 0.1–0.0 Poor, and <0.0 Bad.

Finally, in order to reveal any tendency toward under- or over- estimation, the following equation is used for calculating the average bias [m/day]:

$$BIAS=\overline{Y}-\overline{X}$$

(19)

where $\overline{Y}$ is the mean of the PA’s results and $\overline{X}$ is the mean of the CA’s results.

5. Application of the Proposed Approach in Idealized and Real-Life Cases

In order to assess the performance of the PA, a number of simulations have been carried out, initially in idealized cases, for various sea bottom slopes, median grain sizes, and different incident wave characteristics, and subsequently in a real-life case study.

5.1. Idealized Cases

Two bathymetric regular grids have been constructed with a uniform bottom slope, $m$, of 1:50 and 1:20, respectively. Both grids consist of 1000 cells in the cross-shore and 800 cells in the alongshore dimension with an equal spatial step of dx = dy = 2.5 m, as shown in Figure 4. In both grids the presence of a detached breakwater has been deliberately considered in order to affect the developing longshore wave-induced current and to produce the anticipated patterns of deposition and erosion trends on either side of the breakwater shadow. The breakwaters have been placed at a distance equal to their length seaward from the shore–approximately at a depth of 6.0 m. The median grain sizes considered herein are 0.1 mm, 1 mm, and 2.0 mm, representing fine, medium, and coarse sand. The performance metrics are calculated in a rectangular area incorporating the shadow area of the breakwater, which extends $67.5°$ from each breakwater tip (as shown in Figure 4).

The offshore wave climate presented in Table 1 is considered as the incident to the idealized beaches, with which to carry out the present analysis. Moreover, in order to test the performance of the PA in more than one case of incident wave climate, it is assumed herein that each of the seven considered sectors of the aforementioned wave climate (i.e., from SSW to NNW) could hypothetically represent an independent wave climate. Therefore, the PA is firstly applied considering all the seven sectors in which the area is exposed to [scenarios S01–S06] and subsequently for each individual sector [scenarios S07–S013]. Collectively, the scenarios under investigation are shown in

Table 2. Mean annual frequencies of occurrence per wave height group and sector.

Scenario ID	Incident Wave Climate	$d_{50} \left(\mathrm{mm}\right)$	$m$
S01	Seven sectors	0.1	1:50
S02	Seven sectors	1.0	1:50
S03	Seven sectors	2.0	1:50
S04	Seven sectors	0.1	1:20
S05	Seven sectors	1.0	1:20
S06	Seven sectors	2.0	1:20
S07	SSW	1.0	1:20
S08	SW	1.0	1:20
S09	WSW	1.0	1:20
S10	W	1.0	1:20
S11	WNW	1.0	1:20
S12	NW	1.0	1:20
S13	NNW	1.0	1:20

.

In order to assess the performance of each Scenario based on the PA, the integrated results of the rate of bed level change obtained by implementing the CA was considered as a benchmark. As has been previously mentioned, the CA was composed of 36 sea-states propagating from the seven dominant wave directions as presented in Section 2.

The results obtained by implementing the CA and the PA are showcased indicatively for case S02 in Figure 5.

Inspecting the results visually, it can be deduced that the PA very satisfactorily reproduces the corresponding spatial patterns of erosion/accretion induced by the breakwater’s presence with respect to the CA results considered as a benchmark. The PA slightly overpredicts the magnitude of the erosion rates at the shallower depths behind the tips of the breakwater, however the accretion pattern at the lee of the structure is very adequately reproduced both in extend and magnitude.

The obtained metrics, i.e., R, MSE, BSS, and BIAS, for all the examined scenarios are compiled and shown in

Table 3. Performance metrics of the morphological model adopting the PA.

Scenario ID	$R$	MSE [(m/day)2]	$BSS$	$BIAS$ [m/day]
S01	0.94	4.76 × 10−3	0.82	−3.29 × 10−4
S02	0.95	7.32 × 10−6	0.84	−2.52 × 10−5
S03	0.94	7.79 × 10−7	0.81	−1.09 × 10−5
S04	0.97	9.62 × 10−2	0.71	2.37 × 10−3
S05	0.96	2.72 × 10−4	0.79	−1.21 × 10−4
S06	0.97	8.80 × 10−6	0.85	−3.26 × 10−5
S07	0.93	9.98 × 10−6	0.76	−5.61 × 10−6
S08	0.96	6.70 × 10−5	0.76	−6.95 × 10−5
S09	0.93	8.14 × 10−7	0.73	−3.20 × 10−6
S10	0.90	6.33 × 10−7	0.60	−2.43 × 10−6
S11	0.94	1.42 × 10−4	0.61	−1.94 × 10−5
S12	0.97	1.05 × 10−5	0.82	−2.06 × 10−5
S13	0.96	7.51 × 10−10	0.90	4.46 × 10−8

pertaining to the calculation of the annually averaged rate of bed level change. As has been previously stated, the results of the CA are considered as a benchmark.

As can be observed in Table 3, all the examined Scenarios exhibit relatively high correlation values, with the maximum value reaching 0.97 and the minimum being 0.90. No particular trend can be distinguished for the correlation values with respect to the beach slope ($m$) and the median sediment diameter ($d_{50}$). Most crucially, the obtained BSS values, the most widely used metric to assess a morphological model’s performance (e.g., [12,39,40]), are classified as “Excellent” throughout all Scenarios. This further suggests that the PA is generally not so sensitive to the particular values of $m$ and $d_{50}$ while simultaneously it is well-equipped to assess erosion/accretion patterns from all the dominant wave directions separately. Regarding the model bias, a general underprediction of the bed level change rates is consistently present throughout all scenarios (except S04 and S13), which can be attributed more to the sheer reduction of the simulated sea-states between the PA and the CA. It should be noted that in general the obtained BIAS and MSE values are rather inconsequential compared to the maximum absolute rates of bed level change of each distinct scenario. Hence, larger MSE and BIAS values are observed on S01 and S04, which have the smallest values of $d_{50}$ and are therefore associated with the largest magnitudes of bed level change rates. In summary, all the scenarios carried out in the idealized test cases validate the use of the PA as a means of accelerating the morphological model simulations while simultaneously providing reliable results.

5.2. Real-Life Case Study

The PA will be thereafter implemented and evaluated in a “real-life” case study to assess performance in a more complex situation. The study area named “Therma” is located in the Island of Samothraki, Greece (a general overview of the study area is illustrated in the upper left panel of Figure 6). The area is exposed to waves arriving from WNW to ENE, while the most dominant ones are the ENE and NE sectors. In this area, a shelter has been constructed for accommodating fishing vessels. The fishing shelter is located within the basic littoral transport zone and therefore significant amounts of sediment have been accreted at the east side of the windward breakwater; sedimentation problems occur at its entrance and erosion problems have been observed at the adjacent shore, west of the shelter. Hence, the main focus of the particular case study is the assessment of the erosion/accretion patterns in the vicinity of the said shelter. The existing port infrastructure is shown in Figure 6.

Offshore sea-state wave characteristics were, similarly to the idealized case presented in Section 5.1, obtained by the Copernicus database [16] and specifically the regional product MEDSEA_MULTIYEAR_WAV [18], for a period spanning from 1993 to 2019. The wave characteristics were extracted at an offshore point with coordinates of 40.517942° N, 25.604150° E. As previously mentioned, the study area is predominantly affected by waves propagating from the ENE and NE directions, taking into account the shoreline orientation of W to E, as shown from the wave rose plot in Figure 7.

With respect to the shoreline orientation, seven dominant directions are considered to influence the sediment transport regime and the morphological bed evolution in the study area, namely concerning waves generated from WNW, NW, NNW, N, NNE, NE, ENE. By dividing the wave heights in fixed groups of 0.5 m intervals for each dominant directional sector, the CA is composed of 36 sea-states and will be once again considered a benchmark for the implementation of the PA. The annually averaged initial rates of bed level change obtained by sequentially implementing the nearshore wave, hydrodynamic, and morphological model by the PA will be compared to the corresponding ones obtained by the CA in order to assess the performance of the proposed methodology.

The bathymetry was resolved with a constant spatial step of 2.5 m in both directions, totaling 600 cells in the x direction and 300 cells in the y direction. The bed is comprised mostly of coarse sand material, hence a constant value of $d_{50}$ was considered throughout the numerical domain with a value of 1.8 mm. Model performance evaluation was undertaken in two distinct areas situated to the west and east of the existing fishing shelter. The bathymetry of the study area, as well as the two areas where model results are evaluated, are showcased in Figure 8.

The 36 sea-states of the CA are showcased in

Table 4. Selected 36 sea-states by implemented the CA.

Scenario	Hs (m)	Tp (s)	MWD (°)		$f_g [-]$
1	0.20	2.36	WNW	291.66	4.278640%
2	0.62	4.13		289.33	0.120424%
3	1.23	5.67		286.36	0.011831%
4	1.55	5.21		284.94	0.001268%
5	0.20	2.30	NW	313.42	1.620018%
6	0.66	3.79		314.29	0.049860%
7	1.14	4.45		317.78	0.005070%
8	0.21	2.26	NNW	337.27	0.893250%
9	0.68	3.51		339.39	0.068029%
10	1.21	4.27		337.60	0.016057%
11	1.51	4.74		344.58	0.000423%
12	0.25	2.32	N	1.98	0.943109%
13	0.73	3.45		2.18	0.276764%
14	1.21	4.29		0.81	0.118311%
15	1.68	4.88		0.42	0.040564%
16	2.21	5.55		0.38	0.013521%
17	2.61	5.89		1.66	0.002958%
18	0.28	2.41	NNE	293.90	1.660582%
19	0.75	3.48		295.41	1.322128%
20	1.24	4.28		296.30	1.009871%
21	1.72	4.94		295.67	0.678599%
22	2.23	5.54		295.38	0.289440%
23	2.63	5.86		296.68	0.081550%
24	3.08	6.30		302.36	0.000423%
25	0.31	2.60	NE	316.98	3.925819%
26	0.75	3.54		317.18	5.997955%
27	1.22	4.26		315.81	3.834973%
28	1.70	4.89		314.58	1.401988%
29	2.18	5.38		314.62	0.324088%
30	2.61	5.74		313.34	0.033803%
31	0.30	2.72	ENE	339.34	12.032671%
32	0.72	3.59		336.65	8.999679%
33	1.19	4.31		335.36	2.742284%
34	1.69	4.93		333.14	0.568739%
35	2.18	5.44		332.50	0.144086%
36	2.60	6.11		335.88	0.003803%

, while the 7 representatives obtained through the PA are as shown in

Table 5. Sediment transport equivalent waves for estimating the annually averaged sedimentation and erosion trends in the vicinity of the fishing shelter.

Equivalent Wave	$H_e \left(m\right)$	$T_{pe} \left(s\right)$	Sector	$f_e (\%)$
1	0.70	4.04	WNW	0.27
2	0.60	3.41	NW	0.17
3	0.89	3.72	NNW	0.09
4	1.56	4.67	N	0.20
5	1.87	5.06	NNE	1.66
6	1.52	4.61	NE	5.92
7	1.27	4.36	ENE	6.32

.

To provide a more detailed and concise overview of the simulated wave and hydrodynamic field which constitute the main forcing factors driving the sediment transport and morphological bed evolution of the study area, maps containing the spatial distribution of the significant wave height and current velocity are indicatively shown in Figure 9. The scenario shown concerns the equivalent wave coming from NE with characteristics of $H_e$ = 1.52 m, $T_{pe}$ = 4.61 s and ${MWD}_e$ = 45.29°.

As can be observed from Figure 9, waves propagating from the NE direction generated longshore current gradients directed towards the west. A decrease in current velocities can be observed at the east and west side of the shelter and a relative enhancement of the current velocities is also observed in the front of the shelter’s breakwater.

A comparison between the annually averaged rates of bed level change obtained by assessing the impact of the 36 sea-states of the CA and the 7 representatives obtained through the PA (Table 4) is showcased in Figure 10. In general, alternating patters of erosion/accretion are observed throughout the domain. The model is capable of predicting the accretion area at the east side of the windward fishing shelter’s breakwater and a distinct accretion area can be distinguished at the tip of the breakwater signaling the accumulation of sediment at the port’s entrance. In addition, the erosion area on the western adjacent coast is evident. The PA can very satisfactorily reproduce the accretion and erosion patterns obtained in the benchmark CA. Small differences at the magnitude of the calculated rates are observed in the tip of the breakwater, however the visual inspection of the results between the two methods signals that the PA leads to similar results to the more computationally intense CA.

To verify the conclusions drawn from the visual inspection of the results, and to further evaluate the performance of the morphological model predictions, the calculated statistical metrics for the two evaluation areas are shown in

Table 6. Performance metrics of the morphological model adopting the PA for the real-life case study of the fishing shelter.

Evaluation Area	$R$	$MSE \left[\right(m/\mathrm{day}{)}^2]$	$BSS$	$BIAS [m/\mathrm{day}]$
1	0.87	1.81 × 10−6	0.72	1.09 × 10−5
2	0.95	1.01 × 10−6	0.89	−1.25 × 10−4

.

Overall, the performance of the model implementing the PA is deemed very satisfactory, especially with regards to the BSS classification of [38], rendering the results in both evaluation areas as “Excellent”. The BSS values are slightly lower for Area 1, which may be attributed to the effect of wave diffraction in the calculated patterns of erosion/accretion in the shadow area of the windward breakwater. Overall, the performance of the PA given the complexity of the case study in question is considered very satisfactory.

6. Discussion

In the present paper, a simple yet effective approach to select offshore wave representative conditions to estimate the sedimentation and erosion trends in coastal areas under the combined action of waves and currents was presented. The proposed approach divides the multivariate wave climate in directional bins of fixed size and calculates a sole “annual equivalent” for each directional bin based on the sediment transport potential of each wave record.

The proposed approach (PA) was firstly implemented in an idealized case study of an alongshore uniform beach protected by a detached breakwater. Several simulations were carried out to assess the sensitivity of the methodology on a variety of parameters such as the median sediment diameter, bed slope, and nature of the wave climate. The rate of bed level change results of the seven “equivalent” wave representatives were compared to a classical approach containing a detailed wave climate composed of 36 sea-states considered as a benchmark. Overall, the results of the PA were deemed very satisfactory, and all simulations of the alternative scenarios were classified as “Excellent” with respect to the BSS classification, commonly applied to assess the morphological model performance. Furthermore, the PA exhibited a consistent performance throughout the examined scenarios validating its utilization for a variety of cases.

The PA was thereafter applied in a real-life case study area in Therma, Island of Samothraki, Greece to assess the model performance in a more complex case, with the presence of port infrastructure. Once again, predictions of the rates of bed level change were deemed very satisfactory with respect to the benchmark Classical Approach, rendering the PA suitable for assessment of erosion/accretion patterns in complex real-life cases.

Summarily, the proposed approach is relatively easy to apply and also requires the prescription of a small number of offshore wave representatives, i.e., seven equivalent waves were chosen for our case studies. Therefore, it may be considered as a valuable tool for coastal engineers desiring to accelerate the simulations of the sedimentation/erosion patterns in coastal areas under the combined effect of waves and currents while also maintaining the reliability of the results. As the performance of the PA in both an idealized and a real-life case study was found to be satisfactory and promising, future research will focus on evaluating the capability of the PA in predicting the coastal bed evolution, by performing full morphodynamic simulations.