**2. Materials and Methods**

During the field experiments at the Shkorpilovtsy study site (Black Sea, Bulgaria) [11], we noticed the festoon-shaped features of the shoreline. Satellite images show that the onshore festoon-pattern is accomplished by crescent underwater bars (Figure 1). According to long-term observations [11,12], the bars migrate slightly, depending on variations in the wave conditions. The shoreline shape has inter-annual variations that are possibly associated with the features of the underwater bottom relief. This fact induced us to prove the idea that the location of the underwater bar defines the shape of shoreline to some extent.

**Figure 1.** Satellite image of the crescent underwater bars at the Shkorpilovtsy study site (Image © 2017 Digital Globe). Dotted line: the crest of the bar; solid: the water's edge.

Numerical computation by using XBeach [13] has been chosen as a main tool, because this is a well-developed and popular open-source hydrodynamic and morphology modelling package [14,15]. The non-hydrostatic mode of XBeach, which we used for wave and bottom changes modelling, resolves short waves and provides an accurate reproduction of wave propagation in shallow water [16]. In contrast to stationary or surf-beat mode, the non-hydrostatic mode resolves the wave profile, and it does not require additional wave asymmetry correction.

Although the main goal of the research involves a purely numerical experiment, it is still important to set reasonable realistic boundary conditions. For this purpose, we used field observations (bathymetry and sediment properties) from the Shkorpilovtsy study site [17], as well as information about the wave climate in the north-western part of the Black Sea [18,19].

A numerical grid (1D) was built, based on a set of 12 cross-shore profiles, measured in a frame of the international field experiment "Shkorpilovtsy-2007". All of the observed profiles that were made along the beach had a bar that was located at different distances from shore. From the observed data, we obtained a characteristic shape of the bottom profile and the underwater bar.

The average bottom profile on the Shkorpilovtsy coast has no bar. It has a slope of 0.022, a slight increase of the slope in the upper part, and a small terrace at 2–3 m depth. This profile was considered to be used for modelling, and it was a basis for the creation of a set of barred profiles.

The characteristic shape of the underwater bar was superimposed with a mean profile at a different distance from the shoreline. Thus, five synthetic profiles were created with different bar positionings (Table 1, Figure 2).


**Table 1.** Parameters of the initial profiles.

**Figure 2.** The cross-shore profiles used as the model bathymetry input.

The spatial resolution of the grid was set to 2 m, as calculated using Matlab Toolbox, which has been created and recommended specifically for this purpose, according to the XBeach developers. The Toolbox helps XBeach users to choose appropriate grid settings, taking into account the wave parameters and the relief characteristics. The Black Sea is a non-tidal sea, and so the initial water level was set at 0 m for all of the simulations.

The sediment on the Shkorpilovtsy beach are anisogamous quartz sands. More than 95% of the bottom sediments in the upper part of profile (till 2.5 m) are coarse-grained or medium. For the modelling, we used medium grain, with a diameter D50 (d50) of 0.2 mm.

The validation of XBeach, which has been made by developers and users, shows that this package can be successfully used in non-hydrostatic mode with the default settings [20]; however, some studies have shown that XBeach overestimates coastal erosion [21], which we also noticed from our analysis of the field data and the modelling results. For study site conditions, we tested XBeach with stationary and non-stationary (Joint North Sea Wave Project Spectrum - JONSWAP) wave inputs, and compare these with the synchronous observations carried out during field experiment <Skorpilovtsy-2007> (Figure 3). The results could be considered reasonable.

**Figure 3.** Results of the XBeach verification for the Shkorpilovtsy study site, for cases of 7 h wave action (significant wave height Hs = 0.9 m, and spectral peak period Tp = 7 s).

The goal of the research was to investigate the wave-bottom adjustment on the time scale of one storm. For this reason, we defined the wave input in accordance with typical storm conditions, presented in the wave climate descriptions of Black Sea [18,19]. The wave input was set at the sea boundary of the numerical grid ( ≈860 m from shore) in the form of JONSWAP wave spectra, with the following parameters: peak enhancement factor γ = 3.3, significant wave height Hs = 1.5 m, spectral peak period Tp = 10.5 s, wave energy-spreading angle δ = 2.5◦. A storm with wave heights of 1.5 m was considered to be a dangerous hydrometeorological phenomenon in the Azov–Black Sea region [22]. A wave period of 10.5 s corresponded to extreme storms with return periods of 25 years [18]. The duration of a single storm event was set to 20 h, in accordance with the criteria developed for the Black region by Belberov [23], and for Atlantic coasts by Lozano [24], also taking into account the World Meteorological Organization recommendations for meteorological observations [25].

The modelling process was organized according to the following algorithm. In first hour of wave action, there was no morphology changes, but there was a wave output with very fine time resolution (5 Hz). Modelling was continued for another 20 h, with relief changes enabled. The output of the bathymetry was set every hour.

From the first step of the computation, we obtained free surface elevation data along the profile between the coordinates 620–820 m (see Figure 4). This part corresponded with depths of 2–6 m. The sampling frequency of the time series was 5 Hz. Chronograms were used for the calculation of the wave spectra, and for the following wave parameters [26]:

1. Significant wave height (in m), calculated as:

$$H\_{\mathbb{S}} = 4 \cdot \sqrt{m\_0} \tag{1}$$

where:

$$m\_0 = \int\_0^\infty S(\omega)d\omega \tag{2}$$

and *<sup>S</sup>*(*ω*) are the spectra, *ω* is the angular frequency, with linear frequency filters: 0–0.05 Hz chosen to account for a significant wave height of low frequency range, including infragravity waves ( *HIGW*).

2. The mean wave period(s) is/are as follows:

$$T\_{mean} = \frac{\int^{\infty} S(\omega) d\omega}{\int^{\infty} \omega \cdot S(\omega) d\omega} \tag{3}$$

3. The wave asymmetry coefficient (with asymmetry relative to the vertical axis):

$$As = \frac{\left< \zeta\_H^3 \right>}{\sigma^3} \tag{4}$$

where is the averaging operator, *ζ* are the free surface elevations, *σ* is the standard deviation of the free surface elevations, and *ζH* is the Hilbert transform.

4. The wave-skewness coefficient (relative to the horizontal axis):

$$Sk = \frac{\langle \zeta^3 \rangle}{\sigma^3} \tag{5}$$

From the second step of computation, we obtained the hourly data of the computed morphology changes for over 20 h of wave action. Based on a set of calculated profiles, the coastline retreat and the change in the underwater bottom relief were evaluated.

The parameters of the underwater bar for the initial profiles (input conditions) are shown in Table 1. The wave length was determined by the dispersion relation of the linear waves' theory for the spectral peak period.
