Wave Energy Disbalance as Generator of Extreme Wave Occurrence in Semi-Enclosed Coastal Waters (Example of Rijeka Bay—Croatia)

Abstract

The conditions for the occurrence of high waves in front of the Rijeka port in the Rijeka Bay were analyzed. The analysis was carried out on the basis of measured data on the wave rider station located in front of the main breakwater of the port of Rijeka and the results of numerical wave generation modelings for the wider sea area on the spatial scale of the Adriatic basin. The results of the conducted analysis show that the sudden transition in wind direction from the third to the second quadrant (and vice versa), with the simultaneous rapid increase in wind speed, creates the conditions for generating the largest waves in front of the port of Rijeka. The main reason for achieving the highest wave height in these conditions is the unbalanced wind power input with non-developed surface dissipation (white-capping) and quadruplet wave interaction. Situations with a slower increase in wind speed and approximately constant wind direction resulted in the occurrence of smaller wave heights. The direct application of anemometric data for the forcing wind field in the Adriatic basin within the wave generation model results in a more accurate simulation of wave height and wave period development than application of the wind field from the prediction atmospheric model Aladin-Hr. This is due to the fact that the site is located in a semi-enclosed sea area of restricted fetch, and the spatial/temporal resolution of atmospheric data (2 km and 3 h) is not sufficient to resolve the rapid transition in the wind field. In the case of direct application of anemometric data, the white-capping parameterization should be of a non-stationary character.

1. Introduction

The east coast of the Adriatic has 1246 islands, islets and rocks with prominent channel systems, so it can be considered an archipelago area from a geographic and oceanographic point of view [1] (Figure 1). In such areas, the fetches are relatively short and represent a limiting factor in the surface wave development. In the event of simultaneous occurrence of a sudden increase of wind speed and sudden change of wind direction, the white-capping dissipative process and nonlinear wave interactions (quadruplets) fail to balance the energy introduced by wind action [2]. There is a rapid increase in wave height that can reach the extreme values for a specific location, although the terminal wind speed in this situation is significantly lower than the maximum wind speed in some other situations with a slower transition of wind speed and direction [3].

The occurrence of such a situation was recorded in the Rijeka Bay coastal waters (Figure 1), with the highest significant and maximum wave height (H_S-MAX = 2.06 m; H_MAX-MAX = 4.07 m) measured at wind speed of only 10.5 m/s, after a sudden rise in wind speed and wind direction change from 300° to 170°. Although a number of situations with higher wind speed and approximately constant wind direction (≈ 170°) were recorded during the two-years monitoring period (waves at V1 and wind at W1, Figure 1), all of them ended with a significantly lower wave height than the recorded maximum.

Port of Rijeka is located in the northern part of the Rijeka Bay, which covers 450 km² and is deep up to 67 m. The amplitude of sea surface elevation reaches 0.8 m, its average salinity is 37.5‰, and the transparency is 20 m [4]. The dominant wind is the bora, then sirocco and maestrale.

Information about wave height is of great importance for maritime structures design and planning of different operation activities related to the sea. For the purposes of design, usually, long-term prediction of wave height is provided to gain the information about rare events with return period from 5 to 100 years. Operation activities planning requires real-time and short-term prediction of wave height to gain the information about sea state in the next several hours or days. The process of transferring energy from wind to waves and the consequent generation of wave spectra has not been fully resolved yet [2]. Commonly used numerical methods with parametrically based equations or differential equations include a hint of uncertainty in converting wind energy into wave energy.

Leading contributions to the general theory of wave generation have been given by papers [5,6,7]. A great advancement in the implementation of theoretical foundations into the numerical models has been achieved by the articles [8,9,10]. Pioneer research on the development of wave heights in the field of limited fetches was carried out by Sverdrup and Munk [11] and Bretschneider [12], while Kahma and Calkoen [13] contributed significantly to the research in the more recent period.

Measurement and research of the wind-generated surface waves along the eastern part of the Adriatic Sea in coastal areas are rare [14], while systematic measurements have been carried out in the open sea [15,16].

A significant contribution to solving the problem of nonlinear wave interaction (quadruplets) is given through works [17,18]. The influence of wave energy dissipation through the surface wave breaking (white-capping) and the possibility of its implementation in numerical models were analyzed in [19,20,21].

In Chapter 2 of this paper, an overview of the measured waverider and anemometric data is presented with the selection of the most interesting situations in the monitoring period, including the comments on their specifics. Some details of the applied numerical model are also given in Chapter 2. The results of the analysis on the influence of wind field transition velocity (wind speed and direction) on the development of significant wave heights (H_S) in the case of hypothetical closed area are presented in Chapter 3, along with the results of numerical simulations for the realistic sea area in four selected real situations. The comparison of numerical results with in-situ measurement results is also a part of Chapter 3. Chapter 4 provides conclusions on the possibilities and drawbacks of numerical modeling of wave generation in a semi-enclosed area with the restricted fetch and the conditions of a sudden transition in the wind field.

2. Materials and Methods

2.1. Waverider and Anemographic Data for the Period 1.6.2009.–1.7.2011.

The wave measurements were carried out in the time interval from 1.6. 2009 to 1.7.2011, at the station V1 (φ = 45°19′37.80″ N; λ = 14°23′38.88″ E, WGS 84, Figure 1). Waverider Datawell DWR MKIII, 0.7 m in diameter, with a built-in GPS receiver and a digital data recording device was deployed on the surface, and moored at the depth of 57.5 m. Moored waverider measured wave direction, height, period, and sea surface temperature, so H_S is measured value and not derived through the relationship H_S ≈ H_m0 = 4(m₀)^1/2.

The transfer of data from the measuring position to the recorder was carried out using radio and GSM communication. Recorded data contains standard wave statistics for periods of 30 min. In the same time interval, wind speed and direction were measured with the ultrasonic anemometer (Vaisala WS425A2C3B, minute resolution of data) located at the station W1 in the vicinity of waverider station, at the roof of the Port Control Centre of the Port of Rijeka Authority (Figure 1).

The time series of H_S in half-hour resolution for the monitoring period (1.6.2009.–1.7.2011.) are shown in Figure 2. The maximum value was measured on 19.2.2010, 14:00 (H_S = 2.03 m; H_MAX = 4.07 m; T_P = 5.56 s; Dir = 157.5°). Twenty-two situations with H_S > 1 m were recognized from the complete set of data. The average value of H_MAX/H_S for those 22 situations was 1.82. Generally, value H_MAX/H_S is a function of the number of waves N₀ in a wave train or wave record, and more specifically proportional to squared logarithm of N₀. According to Goda [22], one can expect the value of H_MAX/H_S in the range of 1.6–2.0.

It should be emphasized that two years is not enough to know in detail the wave characteristics in front of the Port of Rijeka. Nevertheless, some conclusions were drawn from the two-year measurement period, primarily those reached in this paper regarding basic statistical characteristics of waves, energy imbalance processes responsible for the appearance of maximum waves, and examination of the possibility of applying numerical models that use measured and modeled values of wind direction and velocity.

Figure 3, Figure 4 and Figure 5 show H_S and wave direction time series, along with wind speed and direction measured at anemographic station W1 near wave station V1 (Figure 1). Representation is given for four situations, 30.11.–1.12.2009, 22–23.12.2009, 19.2.2009, and 6.–7.12.2010.

The first situation (30.11.2009. 12:00–1.12.2009. 0:00, Figure 3) is characterized by an increase in wind speeds from 8 to 13 m/s in the period from 12:30 to 21:30. Wind direction changes very little in the same period (160°–170°). Subsequently, there is a sudden decrease in wind speed and wind direction change (160°–220°), followed by a decrease in wave height.

At the beginning of the second situation (22.12.2009. 21:00, Figure 3), the wind direction transition from 1800 to 1650 can be observed during 2.5 h, with the wind speed increasing from 9.5 to 14.5 m/s, while significant wave heights H_S in this period grow from 1.45 to 1.95 m. In the teaching period of 23.12.2009 23:30, wind speed is reduced to 9 m/s while maintaining approximately constant direction, resulting in the fall of the wave heights.

In the third situation (February 19, 2010. Figure 4), the highest significant and maximum wave heights were measured during the measurement period (H_S = 2.03 m; H_MAX = 4.07 m). Wind direction changes from SSE (157°) to NNW (300°) for 1.5 h (9:00–10:30), while wind speed changes from starting values of 0.5 up to 1.5 m/s. After that, wind direction suddenly returns from NNW (300°) to S/SSE direction (170°) in 2 h time (10:30–12:30), while the speed in the same time period increases from 1.5 to 5.4 m/s. Afterwards, there is a sudden increase in wind speed from 5.4 to 10.3 m/s for only 0.5 h (12:30–13:00). During the period 12:30–13:30 wind direction is constant. Considering the relatively low wind speeds where the highest values of H_S and H_MAX are recorded throughout the measurement period, this situation is a kind of anomaly. It is obvious that the non-stationarity of the wind field plays an important role in the wave generation process.

In the fourth situation (6.12.2010. 05:30–7.12.2010 03:00, Figure 5), the wind changes oscillate from 224° (SW) to 147° (SE/SSE), while the speed changes in the range of 4–9 m/s. For the wind speed, one can notice a slight upward trend until 6.12.2010 17:00, and for the wind direction the trend of transition from SW (224°) to S (175°). Significant wave height also increases from 0.42 m (6.12.2010. 05:30) to the final value of 1.36 m (6.12.2010. 17:30). The increase in wind speed from 05:30 to 17:30 is mild, and it can be concluded that the increase of H_S occurs primarily due to fetch extension during wind direction change from SW to S. After 6.12.2010. 18:00 prevails the trend of reduction in wind speed and direction angle (147° on 7.12.2010 03:00) with a slight decrease in wave heights.

2.2. Applied Numerical Model of Wave Generation

For the purpose of numerical simulations, the spectral wave generation model Mike 21sw was used (www.dhigroup.com). The model relies on spectral formulation according to [23] where the spectral frequency domain is discretized by a logarithmic scale. As part of model routines, the processes of wave-generating with wind, nonlinear wave interactions (quadruplet), refraction, and surface wave breaking (white-capping) have been followed. Diffraction, depth induced wave breaking, sea currents, and surface elevation dynamics have not been taken into consideration due to the minor impact on the analyzed process in the domain of this study. The numerical integration for the source member was carried out according to the methodology presented in the work [24]. Full spectral formulation was used.

The basic equation of the analyzed process is the conservation equation for wave action N(X,t,ω,θ) = E(X,t,ω,θ)/ω, defined by the ratio of wave energy spectra density E and angular frequency ω:

$$\frac{\partial N}{\partial t}+\nabla ·\left(VN\right)=\frac{I+P}{\omega }$$

(1)

where: t time, X(x,y), Cartesian coordinates, V = (c_x, c_y, c_ω, c_θ) wave group velocity in 4-dimensional domain (X,ω,θ), I and P source and sink members, ∇ differential operator in 4-4-dimensional domain (X,ω,θ), θ wave propagation direction, ω angular frequency expressed by equation of linear wave dispersion, $\omega =\sqrt{gh\tanh (kd)}$, k = 2 $\pi$/L wave number, L wave length, d water depth.

Four characteristic velocities c_x, c_y, c_ω and c_θ are defined as follows: (c_x, c_y) = dX/dt, c_ω = dω/dt, c_θ = dθ/dt.

The source member I is defined according to the results of the research presented in the works of Janssen’s [9] in which it was shown that the wave generation intensity varies depending on the elapsed time since wave initialization:

$$I\left(f,\theta \right)=\gamma E\left(f,\theta \right)$$

(2)

where: f = ω/2π wave frequency, γ wave generation intensity.

The wave generation intensity was defined by the expression proposed by Janssen [9]:

$$\gamma =\omega \frac{{\rho }_A}{{\rho }_W}\left(\frac{1.2}{{\kappa }^2}\mu {ln}^4\mu \right)\left[\frac{u_{*}}{c}cos\left(\theta -{\theta }_W\right)\right]$$

(3)

where: ρ_A, ρ_W air and water density, κ von Karmans constant (0.4), θ_W wind direction, θ wave direction, u_* wind friction velocity, c = ω/k wave phase velocity, μ nondimensional critical wave height defined by equation μ = kz₀ exp(κ/m), z₀ wind induced sea surface roughness.

Parameter z₀ is defined with relation:

$$z_0=\frac{Z_{CHARNOCK}{u}_{*}^2}{g}{\left(1-\frac{{\tau }_W}{{\rho }_Z{u}_{*}^2}\right)}^{-1/2}$$

(4)

where: τ_W wind induced sea surface stress, z_CHARNOCK model constant (adopted value 0.01).

Sink member P includes the dissipation process caused by surface wave breaking (white-capping), which is incorporated in the model [23]. The formulation of member P is expressed by the following equation:

$$P\left(f,\theta \right)=-C_{dis}{\left(\frac{\alpha }{{\alpha }_{PM}}\right)}^m\left[\left(1-\delta \right){\left(\frac{k}{\overline{k}}\right)}^2\right]\overline{\omega }E\left(f,\theta \right)$$

(5)

where: Cdis and δ model parameters (subject to calibration), α averaged wave steepness, α_PM wave steepness in Pierson–Moskowitz wave spectra, $\overline{\omega }$ averaged angular frequency, $\overline{k}$ averaged wave number.

Numerical simulations were initially focused on the general influence of wind field transition velocity (wind speed and direction) on the development of significant wave heights H_S in the case of a hypothetical closed sea area. Thereafter, simulations were performed on the scale of the Adriatic basin for the four mentioned situations from the previous chapter. In this case, as the first step, the unstationary and homogeneous wind field with measured velocities at the anemometer station W1 was used for the model forcing. In the second step, the nonstationary and nonuniform wind field from the prediction atmospheric numerical model Aladin-Hr was used for the forcing of the wave generation model. Model Aladin-Hr was used as an operational model in the Croatian national hydrometeorological service with three-hours and two-kilometer resolution of output data [25].

3. Numerical Simulation Results

3.1. Wave Height Change in Closed Sea Area during the Transition of Wind Direction

Using the methodology according to CERC [26], the effective fetch lengths for the location of the waverider station V1 were estimated for SE (18 km), SSE (21 km), S (24 km), and SW (20 km) direction [27]. The rapid change of wind direction from SE to SW direction, or vice versa, periodically occurs in the area of the Rijeka Bay. Just during this transition, high waves were visually observed and instrumentally recorded. During the transition of wind direction from SE to SW (or vice versa), wind also passes through the SSE and S directions with longer effective margins. If the transition is sufficiently slow and without a change of wind speed, there will be an increase of wave height in time of wind passing through the SSE-S sector. On the other hand, increasing the speed of wind direction change (SE → SW or SW → SE) results in a decrease of wave height in the transitional period. Numerical analyses with a hypothetical sea area were made to confirm these statements.

The case of the hypothetical semicircular domain of a 25 km radius was first analyzed (Figure 6). The edges of the model domain represent absorbing boundaries, i.e., there are no open boundaries. The depth is homogeneous and is 50 m. The spatial domain is discretized using a calculation mesh (final volume cells) with distances between the cell centers from 400 m up to 1000 m. In this case, the effective fetches for the KT point and SE, S, SW directions are equal and amount 25 km. Furthermore, if the wind direction changes from SE to the SW direction, the wave height at the KT point will primarily depend on the speed of the wind direction change. Figure 6 shows the cases in which the change of direction SE → SW takes place during 1, 2, 4, and 8 h at a constant wind speed of 18 m/s.

From Figure 6, one can notice that the significant wave height variation at position KT is practically negligible in the case of 8 h duration for wind direction change SE → SW. As opposed to that, when wind direction changes from SE to SW during just one hour, H_S at KT drops during the transition period from 2.3 to 1.95 m.

If the spatial domain is defined by one-fourth of the circle (Figure 7), effective fetch is changed during wind transition from SE to SW direction. It should be noted that the position of KT point is moved southerly for 1.4 km (Figure 7). Effective fetch for SE and SW direction for point KT is equal, and amounts 14.6 km. The effective fetch for S direction is the longest and amounts 25 km − 1.4 km = 23.6 km.

The results of H_S dynamics during the transition periods of 1, 2, 4, and 8 h, using stationary wind speeds of 10, 14, and 18 m/s, are shown in Figure 8.

The results shown in Figure 8 express the following phenomenology. Wind direction change SE→SW through one hour (from 6th to 7th hour, full lines) generates the most intense reduction of H_S, with a more pronounced reduction at higher wind speeds. The decrease is 9% for wind speed 10 m/s, 13% for wind speed 10 m/s and 5% for wind speed 18 m/s. The extension of the transitional period allows for “recognition” of a larger fetch length during the passage of wind through the S direction. In the case of 4 and 8 h transition period, only the increase of wave height is registered in the transition period. For example, in the transition of 8 h and wind speeds of 10, 14, and 18 m/s, the increase is 8.5, 7.5, and 6.5%, respectively. In the case of a transition of 2 h, a decrease and increase of H_S occur (Figure 8). The results of these simulations suggest that the occurrence of the maximum wave height at the position V1 should be expected in situations with slower wind transition from SE or SW.

It should be emphasized that the previously presented statement is in contradiction with the registered wind speed and wave height extremes (H_S and H_MAX) in real situations during the two-year monitoring time interval. In the following chapters of this paper, attention will be focused on recognizing the limitations of numerical model simulations in the emergence of extreme transition in the wind field.

3.2. Numerical Simulations for Real Sea Area

Figure 9 shows the area covered by the spatial domain of the wave generation numerical model (the coverage of the Adriatic basin). In spatial discretization, triangular cells (final volumes) were used. The distance between the numerical nodes, located at the center of each final volume, is variable and extends from 7500 m to the deep-water area of up to 200 m in the zone of the coastline itself. The subject sea area is characterized by a large bathymetric gradient in the coastal zone and the shallow water wave deformations are not essential to the wave generation process. Therefore, the applied calculation step is rational.

The frequency domain was discretized by a logarithmic scale, from the minimum frequency of 0.08 Hz (wave period 12.5 s) to the maximum frequency of 4.53 Hz (wave period 0.22 s), through 40 discrete steps.

Within the model parameterization, homogeneous and stationary values of surface dissipation coefficients C_dis and δ were adopted. When using wind fields from the Aladin-Hr model, Cdis = 3.5 and δ = 0.5 were adopted [28]. The simulations in which the model forcing is defined by the measured wind speed and direction on anemometer W1 are carried out with the adoption of Cdis = 1 and δ = 0.5. In addition to the Otranto Strait, a zero spectrum was used (excluding the influence of wave generation in the Mediterranean Sea on the wave generation process within the Adriatic basin). The initial conditions of all simulations are defined by the wave absence in the whole modeled area.

Figure 10 shows the time series of modeled and measured significant wave height H_S and T_P (peak period) for four analyzed situations with model forcing with wind field according to anemometer data. Figure 11 shows the model fields of H_S for the term of maximum recorded H_S values at the waverider position V1 in situations 1 to 4. The shown H_S fields display only the segment of the model spatial domain in which the Rijeka Bay and the surrounding channel system are located.

The results presented in Figure 10 and Figure 11 demonstrate the fact that, in situations associated with slow wind speed and direction change (situations 1, 2, and 4), anemometric data can be directly used for modeling of wave height dynamic in the vicinity of Rijeka port, or generally in a relatively small and closed sea area. It should be emphasized that the value of the Cdis coefficient is determined by the calibration procedure for each individual aquarium. The final-adopted value of Cdis will be considerably lower than the reference one for the open sea region. The results shown in Figure 10 also indicate the possible improvement in the calculation results for situations 1, 2, and 4, that can be expected through the application of non-stationary values for the coefficient Cdis = f (wind speed).

In the event of a rapid increase in wind speed, with the simultaneous change of direction through a greater direction range (e.g., from 170° to 290°, situation 3), the numerical model forced with anemometric data fails to reproduce the real development of the wave field. Following the basic idea of the imbalance of the energy input by the effect of wind on one side, and the dissipation of energy by white-capping dissipation and energy transfer in the nonlinear quadrupole process on the other hand, the values of the Cdis coefficient and the Charnock coefficient (Equation (4)) have been additionally adjusted. By adopting Cdis = 0 one inactivates the white-capping dissipation, and by decreasing the value of the Charnock coefficient (Equation (4)), one allows an unlimited increase of the value γ in Equation (3). Namely, due to the decrease of z₀ value in Equation (4), a decrease of μ in Equation (3) takes place (μ < 1), and the member ln⁴μ becomes the dominate one in the same equation.

Figure 12 shows the comparison of H_S time series for situation 3, with the fixed value Cdis = 0 and variation of Charnock coefficient value. Model H_S field for the term of maximum H_S value recorded at waverider station V1 in situation 3 is also given in Figure 12. The comparison of Figure 11 and Figure 12 indicates a “dramatic” difference in the resultant H_S field.

The results of model simulations using the wind fields from the numerical atmospheric model Aladin-Hr are given in Figure 13 and Figure 14. Figure 13 shows the time series of modeled and measured H_S and wave direction, along with modeled (Aladin-Hr) and measured wind speed and direction at anemometer site W1. Figure 14 shows the wind fields from the Aladin-Hr model and H_S model field for the term of maximum recorded H_S values at waverider site V1 in situations 1–3.

In situations 1 and 2, model Aladin-Hr gives continuously higher wind speeds and lower angles than measured at site W1. Accordingly, model H_S values in situations 1 and 2 are continuously higher than the measured values. On the other hand, modeled H_S does not reach the measured extremes in situation 3. Apart from the foregoing reasons, the cause for this deviation is also partially due to error in the predicted wind direction obtained by Aladin-Hr model (range 125°–145° in situation 3), resulting in shortening of fetch length in comparison to the accomplished one in measured wind direction range 165°–170° (Figure 13).

4. Conclusions

The measured data from the waverider and anemometer stations, and the results of numerical modeling of wave generation have been analyzed for the purposes of identifying the conditions for the occurrence of the largest waves in front of the port of Rijeka. From the measurement and modeling results, it is recognized that a sudden wind transition from third to second quadrant, with a simultaneous sudden increase in wind speed, results in the emergence of the largest wave in front of the port of Rijeka. The reason for the occurrence of the highest wave heights in these conditions is the more intense rate of wind energy input than the rate of energy dissipation due to white-capping and energy transfer through quadruplet nonlinear wave interaction. In addition, the results of measurements and numerical simulations point to the following conclusions:

• In the case of slow transition of wind speed and direction, anemometer data can be directly used for modeling of wave height dynamic at the location of the port of Rijeka;
• The value of the Cdis coefficient should be determined by the calibration procedure for each individual sea area;
• In the case of direct application of anemometer data for the model forcing wind field, the value of Cdis is to be significantly lower than the referenced for the open sea area;
• The application of non-stationary values for the coefficient Cdis = f (wind speed) results in a more accurate calculation of wave field development;
• It is proposed to use Cdis = 0 up to wind velocity of 10 m/s;
• After the adoption of the value Cdis = 0, the Charnock coefficient becomes dominant, and its adjustment through the value range <0.01 can provide an arbitrary rise in the wave height. It should be noted that the results thus obtained have been “overforced” through the hypersensitivity of the Charnock coefficient and are not a confirmation of numerical model robustness or model resolving possibilities in the process of rapid wave field growing, induced by sudden transitions in the wind field;
• The use of the wind field from a prediction atmospheric model Aladin-Hr with a spatial resolution of 2 km and a time resolution of 3 h is not sufficient to detect extreme transitions in the wind field;
• Extreme wave heights measured in the coastal area of the eastern Adriatic coast are significantly smaller than those in the open sea area.

The processes of wave generation by wind, quadruplet wave–wave interaction, and dissipation by white-capping is not fully resolved in advanced third-generation wave-prediction models. Energy transfer is accounted for with a separate source term in the spectral energy balance equation, but phase-coupling is ignored, although it may be locally important for the short-term statistics of the waves.

In the follow-up of the research, special attention will be paid to finding an optimal strategy for model parameterization in the simulations of the entire monitoring period. Furthermore, since the monitoring period of two years is not enough to know in detail the wave characteristics in front of the Port of Rijeka, we are planning to continue collecting data in the forthcoming period (2020–2025), and validate the results of this study.