The Feasibility of the ERA5 Forced Numerical Wave Model in Fetch-Limited Basins

Abstract

Numerical wave models are critical in hindcasting reliable long-term time series of significant wave heights, which play a crucial role in coastal and ocean engineering activities. Although wind fields are an important input to numerical wave models, few studies have investigated the feasibility of the widely used ERA5 wind reanalysis dataset in fetch-limited basins. In this work, we investigated the feasibility of the ERA5 forced numerical wave model (SWAN) in fetch-limited basins. ERA5 wind velocities were first compared to ground-based meteorological stations, showing poorer accuracy compared to finer gridded ALADIN wind data. Subsequently, the white-capping coefficient C_ds in the Janssen white-capping formulation was calibrated separately using a surrogate model when establishing the ERA5 and ALADIN forced wave models. The calibrated ERA5 forced model showed a similar agreement to wave buoy data as the calibrated ALADIN forced wave model during the calibration period and even superior accuracy in the validation period. Overall, these results show that the wave model calibration procedure mitigates the effect of the poorer accuracy of the ERA5 wind data on the significant wave height results. Nevertheless, both ERA5 and ALADIN forced wave models showed an alarming overprediction for high simulated significant wave heights.

1. Introduction

Knowledge of long-term wave conditions (such as significant wave heights and peak wave periods) is essential for the planning, operation, and maintenance of maritime activities [1], the preparation of flood protection plans [2] the assessment of coastal vulnerability [3], and evaluating harbor oscillations [4,5]. In other words, the role of sea climate is becoming increasingly important in today’s society [6]. However, measurement devices are generally expensive, and their deployment is de facto possible only in a few locations in space and in limited time windows [7]. Since few in-situ wave measurements are available [8], numerical simulations represent an effective option to capture the long-term wind and wave climate [9,10].

Numerical wave models are a powerful tool that can reproduce wave climate in a flexible manner and are an excellent method for extending existing wave data [11,12]. In recent years, the accuracy of wave model simulations has improved considerably [13]. These improvements are due to improved mathematical formulations of the source terms in the governing equations and more accurate and detailed wind fields from atmospheric models [13]. On the other hand, numerical models forced on regional domains require local meteorological and oceanographic information, especially at the open boundaries. These data can be obtained using hindcast model results of global models; such procedure is usually referred to as “downscaling” [7,14,15,16]. It can be performed in three different ways as described by Camus, Mendez [17], with the dynamic method being the most popular of the three. In the dynamic method, the local wave conditions are computed on a regional (nested) grid by replicating all meteorological conditions, thus modeling the coastal wind-wave propagation in more detail.

Wind data is clearly a key forcing parameter for numerical modeling of the wind-generated waves [18,19], which is consistent with observations in nature where wind is the dominant driving force. As noted in previous research, the quality of the wind has a major impact on the accuracy of wave modeling results [20,21]. Therefore, accurate historical wind data is an essential input for wave climate hindcasting, whether it is a simple empirical model [22,23] or a complex numerical wave model [18,19]. Due to the large dependence of the wave model uncertainty on the wind data accuracy [24], Cavaleri [25] showed that a 10% wind error can lead to 10–20% error of the simulated wave height under fully developed conditions. Similarly, Teixeira, Abreu [26] reported that errors in the simulated wave height can be up to twice the amount of the relative error in the model wind data, and Wu, Li [27] showed that a small underestimation of wind velocity leads to large underestimations of the measured wave heights.

Fortunately, the quality of wind reanalysis products has improved significantly in recent decades due to improvements in the physical components and numerical procedures of the model, as well as greater computational resources, which in turn has had a positive impact on the results of numerical wave models [27]. At the forefront of wind data reanalysis accuracy in recent research is the European Centre for Medium-Range Weather Forecasts atmospheric reanalysis wind model ERA5 [28], which has been shown to be very powerful in coastal waves modeling [7,27,29]. Wu, Li [27] conducted a comparative analysis of seven different wind reanalysis datasets for the Yellow Sea, in China and concluded that the ERA5 model is superior to other wave reanalysis data sources tested. Similar conclusions were drawn in studies conducted in the Mediterranean Sea. For example, Amarouche, Akpınar [30], Barbariol, Davison [31] and Korres, Ravdas [32] successfully used ERA5 reanalysis wind data to force the wave numerical model of the whole Mediterranean Sea to assess the wave climate on a scale of 40 years. Further, Bellotti, Franco [5] showed the applicability of ERA5 on a regional model of a part of the Mediterranean Sea (western mid-latitude coast of Italy) and the applicability of the downscaling procedure. However, there are no islands protecting the coast and the wave measurement sites in the considered region.

However, both Barbariol, Davison [31] and Korres, Ravdas [32] cautioned against wave modeling using ERA5 wind data in semi-closed and closed basins. Korres, Ravdas [32] reported the weakest model accuracy in the Adriatic Sea compared to other regions of the Mediterranean Sea. Barbariol, Davison [31] also reported a particularly strong tendency of the wave model to underestimate the significant wave height in the Adriatic Sea. Both argue that the numerical domain of the Adriatic is shallow, closed, and constrained by complex topography with narrow and limited fetch, leading to higher error in the wave model. The spatial resolution of the wind data and the wave model itself are mentioned as the main causes [32]. In other words, dynamic downscaling is recommended in this situation [31]. Several studies [33,34] have discussed the same general principle as in [31,32], that the complex orography in coastal and semi-enclosed areas has a greater effect on the accuracy of the estimates of the local wind field, and hence the wave field, than of the offshore wind field. Son and Do [35] argued that wind data with adequate spatial resolution are especially necessary for predicting extreme significant wave heights that occur during storms. In line with the recommendations in [31,32], Christakos, Björkqvist [36] and Vannucchi, Taddei [37] downscaled the ERA5 and ERA-Interim wind product to model systems as the fjords in Norway and the Tuscan Archipelago, respectively. On the other hand, Bellotti, Franco [7] successfully downscaled only the resolution of the wave model using a nested grid while still using the ERA5 wind data.

Nevertheless, a promising, publicly available reanalysis database, such as ERA5 should not be simply dismissed, even in complex coastal basins, such as the eastern Adriatic Sea. The ERA5 reanalysis provides relatively high spatial resolution (0.25° × 0.25°), hourly temporal output, and 3-hourly uncertainty information compared to other publicly available atmospheric reanalysis models. Few of the other publicly available reanalysis models have similar or more detailed spatial and temporal resolution (e.g., JRA-55 with 1.25° × 1.25° [38], MERRA-2 with 0.5° × 0.625° [39], ERA-Interim with 0.25° × 0.25° [40]).

This study aims to evaluate the feasibility of the widely used and publicly available atmospheric reanalysis product ERA5 [28] for hindcasting significant wave heights using the numerical wave model SWAN at a regional scale by the downscaling procedure in a fetch-limited, semi-enclosed basin (study area, numerical model description, and setup are described in Section 2.1, Section 2.2, Section 2.3 and Section 2.4). The feasibility of ERA5 wind data in wave modeling in a fetch-limited area in the Adriatic Sea is evaluated by comparing wave model results during the calibration and validation periods with: 1. wave measurements, and 2. results using a more detailed ALADIN wind product (described in Section 2.5) and its results for significant wave height. The calibration procedure is performed separately for each ERA5 and ALADIN wind product by adjusting the white-capping coefficient (described in Section 2.6). This study provides important insights into the feasibility of ERA5 wind data in fetch-limited catchments, providing future users with expected margins of error when using ERA5 wind data in such areas (results and discussion in Section 3 and Section 4).

2. Materials and Methods

2.1. Study Area and Field Observations

The research area is located near Split, Croatia, on the mid-latitude Adriatic Sea, as shown in Figure 1. The topography of the region is complex, with two islands southwest of Split (Šolta and Vis), three larger and one smaller island southeast of Split (Brač, Hvar, and Korčula, and Lastovo, respectively), with several other small islands. The coastal area has developed economically, especially in terms of beaches and marinas, and has a rich ecological and archaeological value. Boats operate frequently on the navigational routes between the islands, especially during the summer season.

The wave climate in the region is mainly determined by bora winds (NE) and scirocco winds (SE). The wave buoy measurements used in this study are located off the port of Split (43.48833° N, 16.46500° E—shown with a red circle in Figure 1). The wave buoy provides hourly wave parameter data, which are retimed into daily mean and max values for clarity in Figure 2. The measured significant wave height time series in Figure 2 shows a typical year regarding the wave climate in front of Split from November 2007 to November 2008. The wave climate is harsh during the winter with calm periods in between stormy events (from November 2007 to early April 2008 and mid October 2008 to November 2008). On the other hand, during the summer there is a longer period of relatively calmer wave climate (May to October 2008). The maximum daily significant wave height for the considered measuring period reaches 1.29 m on 6 March 2008, at the most extreme for this fetch-limited location (Figure 2).

Additionally, Figure 2 shows the calibration period (25 February 2008, to 9 April 2008) for the calibration of the wave numerical model in red and validation period for testing the calibrated wave numerical model in blue (15 September 2008, to 15 November 2008). Both periods show prominent energetic events, which deem them appropriate for wave model calibration and validation. The extent of the numerical model domain is shown in Figure 1, while detailed information regarding the numerical model is presented in Section 2.2, Section 2.3, Section 2.4, Section 2.5 and Section 2.6.

The mean wave height during the measuring period for the buoy location is significantly lower at 0.21 m (Figure 2 and Figure 3—left). Most of the measured significant wave heights fall in the range of 0.1–0.2 m. The measured wave peak periods are, unfortunately, only available for the latter part of the measuring campaign (August 2008 to November 2008); subsequently Figure 3—right shows the relationship between the significant wave height and peak wave period during the validation period. For the lower range of significant wave heights (H_s < 0.25) the peak wave period is erratic, ranging from 2 to 6 s. This is a typical behavior that has also been observed in previous research. For higher significant wave heights (H_s > 0.2), the relationship shows more consistency.

DHMZ (Croatian Meteorological and Hydrological Service) maintains several automatic meteorological stations located at various points along the coast (shown by rectangles in Figure 1). Wind measurements, including hourly wind speed and wind direction, obtained from four stations located in Split, Resnik, Hvar, and Ploče, are used to verify the accuracy of the ERA5 wind reanalysis data at each location (results in Section 3.1). A complete list of stations and a description of measurement techniques can be found on the DHMZ website [41].

2.2. SWAN Model Description

The third generation of the SWAN spectral wave model (version 41.41 is used in this paper), developed in Fortran90 at Delft University of Technology, is known for its exceptional performance in simulating wave propagation and wave growth due to wind power. SWAN is based on the dynamic spectral balance equation and is represented in Cartesian coordinates as follows:

$$\frac{\partial N}{\partial t}+\frac{\partial C_{x}N}{\partial x}+\frac{\partial C_{y}N}{\partial y}+\frac{\partial C_{\sigma }N}{\partial \sigma }+\frac{\partial C_{\theta }N}{\partial \theta }=\frac{S}{\sigma },$$

(1)

$$N=\frac{E}{\sigma },$$

(2)

where the left-hand side represents the kinematic component and shows the derivatives of the active density N in the geographic and spectral space. C_x and C_y are the Cartesian coordinate components of the wave group velocity in geographic space, while C_θ and C_σ are the spectral components of the wave group velocity. The effective density N is the ratio of the wave energy spectrum E(σ, θ) and frequency σ. The term S on the right-hand side represents the total source of wave energy. The total source term is composed of the following terms:

$$S=S_{in}+S_{ds}+S_{bot}+S_{br}+S_{nl3}+S_{nl4},$$

(3)

where S_in is the energy input due to wind, S_ds is the dissipation due to white-capping, S_nl4 is the nonlinear wave energy transfer between quadruplets, S_nl3 is the triad nonlinear interaction, S_fric is dissipation due to bottom friction, and S_brk is dissipation due to depth-induced wave breaking. SWAN allows the user to turn source terms on and off, as well as choose different formulations for each source term in Equation (3). The full description of the numerical wave model can be found in the scientific and technical documents of SWAN [42].

2.3. SWAN Model Setup

We used an unstructured mesh for our downscaling procedure to accurately describe the complex shoreline geometry in the Split region (Figure 4). Similarly, studies that used grids coarser than 2 km recommend the use of unstructured grids to accurately resolve complex shoreline features, especially in the presence of islands due to their sheltering effect [43]. On the other hand, the unstructured grid allows the use of coarser elements away from the coastline to save computational costs for larger domains [44]. In particular, the Mediterranean Sea has a complicated morphology and environment, including a very complex coastline, highly variable bathymetry, and limited fetches. Therefore, an unstructured mesh allows flexible local refinement that can improve the accuracy of wave simulation for specific nearshore areas [37]. The domain of the numerical model in this study covers part of the eastern Adriatic Sea and extends over 16° E–17.6° E in longitude and 42.9° N–43.5° N in latitude.

The source of the basic bathymetric data is the latest GEBCO database [45]. In addition, local bathymetry measurements were added to the nearshore baseline data to provide finer resolution of the coastline. The resolution of the baseline bathymetric data is approximately 200 m × 200 m offshore, while the resolution of the nearshore was finer at 5 m × 5 m. To accurately represent the complex coastal geography, the spatial resolution of the nearshore data is 5 m. The unstructured grid area consists of 88,355 triangular elements, with grid elements decreasing in size to the shoreline (Figure 4). The depth of the deepest water is approximately 120 m below mean sea level (MSL) at the western open boundary (Figure 4).

SWAN makes the input parameters and the choice of empirical formulations and their constants flexible for the user. Therefore, it is possible to develop an accurate modeling system for a given region by carefully listing the formulations included in the model and calibrating them [46].

In the SWAN model setup, the formulations for the source terms from Equation (3) are adopted from the selection of formulations provided as part of the SWAN package. SWAN default formulations are chosen for the nonlinear quadruplet wave interactions that adopts the Discrete Interaction Approximation (DIA) formulation [47], the breaker index formulation [48], the JONSWAP bottom friction formulation [49], and the triad wave-wave interaction that uses the Lumped Triad Approximation [50] since our domain includes intermediate depth and shallow water areas. A summary of the empirical formulas and parameter values used is shown in Table 1.

The wind energy input term in SWAN, S_in in Equation (3), is defined as the sum of the linear and exponential growth terms. The linear term is based on the study of Cavaleri and Rizzoli [51], while there are several expressions available for the exponential term. First, Komen, Hasselmann [52] adopted the expression of Snyder, Dobson [53] to introduce a revised formulation expressed in terms of the friction velocity. Moreover, Janssen proposed a new quasi-linear wind wave formulation [54,55]. On the other hand, the wave dissipation term, S_ds in Equation (3), determines the amount of white-capping dissipation in the domain. Komen, Hasselmann [52] formulated the pressure pulse model for white-capping, which is the default formulation in SWAN. However, in this work, the white-capping formulation in conjunction with the wind energy input formulation given by Janssen [54], Janssen [55]. This choice is due to the successful application of the Janssen formulation in recent research in the Mediterranean [30,56].

Table 1. Summary of used empirical formulations and parameter values.

SWAN Parameter Setting	Formulation	Par. Value	Determination Method
Directional standard deviation	-	2	Literature [42]
Directional wave resolution	-	36 bins or 72 bins	Literature [42] and Calibration procedure (Section 2.6)
Frequency range	-	0.05–1 Hz	Literature [42]
Number of spectral frequencies	-	30 bins	Literature [42]
Triad constant trfac	Eldeberky [50]	0.05
Triad constant cutfr	Eldeberky [50]	2.5
Quadruplet constant Cnl4	Hasselmann, Hasselmann [47]	3 × 107	Literature [27,42]
Quadruplet constant λ	Hasselmann, Hasselmann [47]	0.25	Literature [27,42]
Breaking constant α	Battjes and Janssen [48]	1.00	Literature [27,42]
Breaking constant γ	Battjes and Janssen [48]	0.73	Literature [27,42]
Bottom friction Cfric	Hasselmann [49]
Number of iterations	-	6 iterations	Literature [6,29,57]
Time step size	-	20 min or 10 min	Literature [36,56,57] and Calibration procedure (Section 2.6)
White-capping δ parameter	Janssen [54], Janssen [55]	1	Literature [29,56,58]
White-capping parameter Cds	Janssen [54], Janssen [55]		Calibration procedure (Section 2.6)

Table 1. Summary of used empirical formulations and parameter values.

SWAN Parameter Setting	Formulation	Par. Value	Determination Method
Directional standard deviation	-	2	Literature [42]
Directional wave resolution	-	36 bins or 72 bins	Literature [42] and Calibration procedure (Section 2.6)
Frequency range	-	0.05–1 Hz	Literature [42]
Number of spectral frequencies	-	30 bins	Literature [42]
Triad constant trfac	Eldeberky [50]	0.05
Triad constant cutfr	Eldeberky [50]	2.5
Quadruplet constant Cnl4	Hasselmann, Hasselmann [47]	3 × 107	Literature [27,42]
Quadruplet constant λ	Hasselmann, Hasselmann [47]	0.25	Literature [27,42]
Breaking constant α	Battjes and Janssen [48]	1.00	Literature [27,42]
Breaking constant γ	Battjes and Janssen [48]	0.73	Literature [27,42]
Bottom friction Cfric	Hasselmann [49]
Number of iterations	-	6 iterations	Literature [6,29,57]
Time step size	-	20 min or 10 min	Literature [36,56,57] and Calibration procedure (Section 2.6)
White-capping δ parameter	Janssen [54], Janssen [55]	1	Literature [29,56,58]
White-capping parameter Cds	Janssen [54], Janssen [55]		Calibration procedure (Section 2.6)

The Janssen white-capping formulation requires the input of 2 parameters, δ and C_ds. The value for the δ-parameter determines the dependence of the white-capping on the wave number. Komen, Hasselmann [52] originally suggested a value of 0.5, but recent research has shown more accurate model results with a value of 1.0 [58]. The latter recommendation is used in this paper. The dissipation coefficient C_ds is determined through a calibration procedure described in Section 2.6.

2.4. SWAN Model Boundary Conditions and Run Periods

We applied wave reanalysis data to the open boundary to the west and south. Namely, the MEDSEA reanalysis data from the Copernicus database [32] (Figure 4). MEDSEA is a fine resolution (1/24°) wave model of the Mediterranean Sea that provides wave hindcasts back to 1993. Although the MEDSEA reanalysis model is the most detailed numerical wave model reanalysis for the Mediterranean Sea provided by the Copernicus project, the authors found low accuracy when verifying the reanalysis data using buoy measurements at well-protected sites. Therefore, the open boundaries of the model are positioned to have no significant protective effects from nearby islands. The shape of the spectrum (in frequency and direction) of the drive wave is defined by the JONSWAP spectrum with a peak gain parameter (gamma) of 3.3.

The model is run in non-stationary mode with a time step of 20 min and up to 6 iterations per time step, while outputs are generated and saved hourly and are written in ASCII format inside a text file. The wave simulation model for calibration purposes is run for a period of 43 days during late winter and early spring (25 February 2008, to 9 April 2008). The first 3 days are considered as model warm-up from a cold-start and are not considered in the evaluation of model performance. This is considered appropriate for the region in question because winter storms are prevalent in the Adriatic Sea, with long periods of calm between, and there were relatively strong wave storms during this period. The validation period is subsequently run with the calibrated white-capping coefficient to evaluate the wave model accuracy against the measured significant wave height and peak wave period.

The calibration period was chosen as approximately one and a half months due to the restricted computation resources and the high number of calibration model runs (75 runs), in order to make the work feasible. Due to this constraint a period, the most prominent wave energetic events are chosen from the first part of the measuring campaign (November 2007 to May 2008).

2.5. Wind Forcing Data

The accuracy of the numerical wave model is tested when it is independently forced from two different re-analysis wind data sources. These include the well-known ERA5 wind data set and the more detailed ALADIN wind data.

First, the new fifth generation atmospheric reanalysis, ERA5 [28], is tested as the forcing wind dataset. ECMWF has developed the dataset itself and distributed it through the Copernicus service, which increases the demands on the reanalysis product. Coastal wind accuracy has increased since ERA-Interim as advances have been made in model formulation and technical capabilities. The ERA5 wind dataset obtained has an hourly time resolution and a spatial resolution of 0.25° (31 km) in both latitude and longitude (Figure 5 and

Table 2. Wind data sources used in this study with corresponding spatial and temporal resolutions.

Wind Data	Type	Spatial Resolution	Temporal Resolution	Location/Region	Altitude
ERA5	Grided data	0.25° × 0.25°	1 h	Global
ALADIN	Grided data	0.07° × 0.07°	3 h	Adriatic Sea
Hvar	Wind measurements	Point data	1 h	43°10′15″ N 16°26′14″ E	20 m
Resnik	Wind measurements	Point data	1 h	43°32′22″ N 16°18′5″ E	19 m
Split	Wind measurements	Point data	1 h	43°30′30″ N 16°25′35″ E	122 m
Ploče	Wind measurements	Point data	1 h	43°2′51″ N 17°26′34″ E	2 m

).

Secondly, the ALADIN numerical weather prediction project (Aire Limitée Adaptation dynamique Développement InterNational) is a collaboration in the field of numerical weather prediction (NWP) to build limited-area models for operational weather forecasting [59]. As an ALADIN member, DHMZ has been running the ALADIN model for its area of interest, which includes the Adriatic Sea, since July 2000. The model focuses on a small area with a high spatial resolution so that fine scales can be resolved (local winds, breezes, thunderstorms, etc.). The numerical model wind dataset obtained for this study has a spatial resolution of 0.07° (8 km) in both latitude and longitude (Figure 5) and a temporal resolution of 3 h (Table 2) [60]. For the operational validation and verification procedure of the ALADIN model, the reader is referred to Tudor, Stanešić [61].

For the comparison between gridded reanalysis wind data and measured wind data at meteorological stations, the reanalysis gridded data are sampled and interpolated at the wind station locations using the numerical model SWAN methodology.

2.6. Model Calibration Procedure and Statistical Error Metrics

As a crucial parameter in the numerical model calibration procedure, the white-capping dissipation coefficient C_ds stands out. Several previous wave-hindcast studies have shown that the default C_ds in the SWAN setup often leads to results that need to be validated [8,13]. Since C_ds is the least known physical term of the SWAN wave model, it is commonly used as a tuning parameter within the calibration procedure [13,62]. This simple solution aims to minimize the discrepancy between simulated and measured significant wave height due to errors in the wind fields. Even though this method does not have a valid physical background, since the simulated white-capping dissipation may not reflect the realistic conditions, this procedure is considered to be quite efficient for solving realistic wave simulations.

A surrogate model was developed to define the automated procedure of the C_ds calibration. The algorithm of this model uses Radial Basis Functions (RBF) to construct a low-cost model for the objective function interpolation [63,64]. The procedure itself involves occasional calls to the SWAN numerical model to evaluate the objective function near the current incumbent point. The maximum number of calls to the SWAN model that the surrogate model processes is fixed at 15 after which the calibration algorithm is stopped. The surrogate model iterates the C_ds value to minimize the loss function defined by the error Equation (5) (Section 3.1). The search range for the white-capping parameter C_ds varies between 0 and 5, as the white-capping parameter is commonly limited to this range in formulation by Janssen [54], Janssen [55].

In addition, a sensitivity analysis was performed regarding wave direction resolution to test the influence of the wave direction resolution in mitigating the restricted wave diffraction formulation documented in SWAN manual. Wave diffraction could be a significant factor in a complex topography like the Split region. The number of directional bins was tested for 36 and 72 bins. Additionally, the influence of a 10 min time step was also tested. The results are presented in Section 3.2.

SWAN model uncertainty is examined by the Pearson correlation coefficient (R), the corrected indicator HH proposed by Hanna and Heinold [65], the normalized bias (NBIAS), and the normalized root mean square error (NRMSE), defined in Equations (4)–(7) respectively:

$$R=\frac{{\displaystyle \sum _{i=1}^N\left(\left(P_i-\overline{P}\right)\left(O_i-\overline{O}\right)\right)}}{{\left[\left({\displaystyle \sum _{i=1}^N{\left(P_i-\overline{P}\right)}^2}\right)\left({\displaystyle \sum _{i=1}^N{\left(O_i-\overline{O}\right)}^2}\right)\right]}^{1/2}}$$

(4)

$$HH=\sqrt{\frac{{\displaystyle \sum _{i=1}^N{\left(P_i-O_i\right)}^2}}{{\displaystyle \sum _{i=1}^N{P}_i{O}_i}}}$$

(5)

$$NBIAS=\frac{\overline{P}-\overline{O}}{\overline{O}}$$

(6)

$$NRMSE=\sqrt{\frac{{\displaystyle \sum _{i=1}^N{\left(P_i-O_i\right)}^2}}{{\displaystyle \sum _{i=1}^N{O}_i{}^2}}}$$

(7)

where P_i is the ith prediction, O_i the ith observation, the overbar denotes the mean values, as in $\overline{P}$ is the mean of all prediction values, and $\overline{O}$ is the mean of all observation values. HH is generally recommended when quantifying dispersion errors because it is unbiased toward simulations with negative biases, which is in contrast to normalized root mean square error (NRMSE) and its variants [66]. For both HH and NRMSE performance indicators, a larger value indicates a higher dispersion error, while a lower value points to a lower dispersion error when comparing measured and modeled data. Both indicators are always non-negative and show perfect alignment between measured and modeled data at a value of 0.

3. Results

The results of the paper are divided into 3 separate sections. First, in Section 3.1, a comparison is made between the wind reanalysis data from the ERA5 dataset and the ALADIN dataset and measurements from meteorological ground stations using automatic wind instruments. Then, Section 3.2 presents the results of the C_ds calibration procedure using a surrogate model (described in Section 2.6) along with the statistical parameters of the accuracy of the calibration period. Finally, using the best C_ds white-capping parameter (from Section 3.2), wave models are run separately using the ERA5 and ALADIN wind datasets for the validation period, with results presented in Section 3.3. Modeled significant wave heights and peak wave periods are compared to buoy measurements during the validation period to determine model accuracy.

3.1. Comparison of ERA5 and Aladin Wind with Meteorological Station Measurements (Accuracy of the Forcing Wind Field)

It is commonly argued that the quality of forcing wind fields has a large impact on the accuracy of the wave models, and therefore, it is appropriate to better understand the accuracy of the gridded forcing wind fields.

As shown in Figure 6, the ALADIN wind data show better agreement with the ground measurements for all performance metrics compared to the ERA5 wind data. This is not surprising since the ALADIN wind dataset has a much finer resolution, which is better suited for complex topographic regions. The finer resolution is able to model the sudden spatial changes in magnitude and direction that can occur in coastal and hill regions. The Split and Ploce wind stations are best predicted by the ALADIN reanalysis with the lowest HH and NRMSE indices, while the ERA5 reanalysis has similar difficulty modeling all tested wind stations with the HH index being similar across the board. The ERA5 wind reanalysis shows a very high overprediction in the case of the Hvar wind station at 0.52, while NBIAS is still prominent for Split and Ploce at −0.24 and 0.29, respectively. ALADIN reanalysis also shows significantly better results in terms of NBIAS with the highest overprediction at Resnik at 0.19. Overall, the ALADIN reanalysis has the most difficulty predicting wind data at Resnik, which is in the NW part of the Split region. On the other hand, the ERA5 has the highest difficulty in predicting the wind at the Hvar station, which is in the center of the domain.

In summary, the wind data from the ALADIN reanalysis agree better agreement with measurements compared to the data from the ERA5 reanalysis in the Split region, which is characterized by complex topography. Every tested performance index (NBIAS, HH, R, and NRMSE) points to the same conclusion.

3.2. Wave Model Calibration Results for ERA5 and ALADIN Wind Fields for the Calibration Period

In this section, we present the result of the white-capping parameter calibration procedure described in Section 2.6. The white-capping coefficient C_ds is calibrated separately for both wave models using ERA5 and ALADIN wind data. For the calibration procedure 75 simulations of the numerical model were run with various white-capping coefficients, number of wave direction bins (36 bins and 72 bins) and time step sizes (10 min or 20 min).

The results of the surrogate optimization algorithm reveal a relationship between the significant wave height and the white-capping parameter (Figure 7). The relationship reaches a global minimum of the HH performance metric at a value of 1.25 and 0.35 for the ERA5 and ALADIN white-capping coefficients, respectively (Figure 7a). The algorithm showed a significant difference between the value of the calibrated white-capping parameter for this particular case and the SWAN default value for the white-capping coefficient of 4.5 for the Janssen [55] formulation.

When considering the NRMSE index instead of the HH index (Figure 7b), the ERA5 wind reanalysis still shows a better performance compared to the ALADIN wind reanalysis. However, the global minimum has shifted a bit, with the NBIAS now decreasing from −0.10 to −0.18, while the optimal white-capping coefficient C_ds increased from 1.25 to 2.0 for the ERA5 wind reanalysis. The same trend is obtained for the ALADIN drive wave model with the NBIAS decreasing from −0.12 to −0.18 and the optimal white-capping coefficient C_ds increasing from 0.3 to 0.6. In summary, if relying on the NRMSE coefficient, an almost two times higher C_ds would be chosen with an even larger underprediction.

An increase in the wave directional resolution bins from 36 to 72 was considered during the calibration procedure. It was proposed that this might mitigate the error due to the limited formulation of the diffraction documented in the SWAN manual [42] (Figure 7c,d). The ERA5 driven wave model would have a slightly increased accuracy with the HH index decreasing from 0.44 to 0.43 (a decrease of 2.5%), while the NBIAS decreased from −0.11 to −0.10 (a decrease of 9%). While this increase in accuracy is noteworthy, the increase in computational resources is approximately threefold. Additionally, the ALADIN forced model did not show any increase in accuracy (Figure 7d). The effect of a smaller time step from 20 min to 10 min had a negligible effect on model accuracy (not shown here).

Figure 8 shows the scatter plots of simulated H_s with ERA5 and ALADIN wind forcing compared to wave measurements at the Split wave buoy with the best performing white-capping coefficient included (according to Figure 7a)—C_ds = 1.25 for ERA5 and C_ds = 0.35 for ALADIN) for the calibrating period. To quantify the performance of the wave model, we calculated the error metrics (R, HH, NRSME, and NBIAS, as described in Section 2.6) of the simulated H_s (Figure 8).

No significance was found between the overall accuracy of the wave numerical models forced by ERA5 and ALADIN. This is a rather surprising result. The wave model forced by ERA5 wind reanalysis data (R = 0.81, NBIAS = −0.10, NRMSE = 0.42, HH = 0.44) shows a similar overall performance to the wave model forced by ALADIN (R = 0.81, NBIAS = −0.12, NRMSE = 0.44, HH = 0.46). Both wave models show a slight underestimation of the measured significant wave heights (NBIAS = −0.08 and NBIAS = −0.11 for ERA5 and ALADIN, respectively). Scatter is slightly better for the ERA5 model, but it is not a significant difference (HH = 0.44 for ERA5 and HH = 0.46 for ALADIN).

Taken together, these performance metrics indicate that the calibration procedure was able to compensate for the higher wind errors in the ERA5 wind data (Figure 8).

Figure 9 shows the relationship between the model accuracy performance metrics (HH, NRMSE and NBIAS) and the simulated significant wave height. For both the ERA5 and ALADIN forced wave models, the accuracy of the modeled significant wave heights generally increases as the precited significant wave height itself increases (e.g., the HH metric decreases from 0.71 for low wave heights to 0.23 for high wave heights when considering the ALADIN simulation). The worst performance for the wave model forced by the ALADIN wind is found for low wave heights (HH = 0.71) and a smaller decrease in accuracy (HH = 0.39) for simulated wave heights of 0.7–0.9 m. However, the wave models forced by ERA5 show low accuracy for relatively high wave heights of 0.9–1.1 m, with the HH metric increasing to 0.48. The ERA5 performance metrics show that caution is needed because the high significant wave heights are the most important for calculating storm energy, defining significant wave heights for a given return period, etc.

Both the ERA5 and ALADIN forced wave models show an underestimation of the low simulated significant wave heights (NBIAS = −0.26 for the ALADIN forced model and −0.04 for the ERA5 forced model). An increasing trend is observed with increasing simulated significant wave height, with wave heights above 0.9 m being greatly overestimated (NBIAS = 0.25 for the ALADIN forced model and NBIAS = 0.39 for the ERA5 forced model). This overestimation can also be seen in Figure 10 for 24 March 2008. As can be seen in Figure 10, the simulated time series of significant wave height reasonably follows the trends of measured wave heights. No obvious delay in the rise and fall of the simulated wave heights could be observed when the wave models are forced with either ERA5 and ALADIN (e.g., the rise and fall of the significant wave height during the storm on 24 March 2008). The overestimation observed in the scatter plots in Figure 9 for high wave heights is also observed here on 24 March 2008 and 25 March 2008. The forced ERA5 wave model also overestimates the measured wave heights during the weaker high-energy event on 22 March 2008. The Taylor diagram on Figure 10 shows that the statistical indicators are similar for both ERA5 and ALADIN forced wave models when compared to the measured time series during the examined time period.

3.3. Accuracy of Modeled Wave Parameters Using Calibrated Whitecapping Parameters Forced by ERA5 and ALADIN Wind Fields for the Validation Period

In this section, the wave model was ran for the validation period, which was not considered in the process of calibrating the white-capping parameter. The wave model results are compared with wave buoy data in front of Split to examine the performance of wind fields in wave modeling for the coastal region with complex topography.

Figure 11 shows the scatter plots between simulated and measured significant wave height at the Split wave buoy for the validation period. The best performing white-capping coefficient from the calibration procedure is implemented into the wave numerical model (according to Figure 7a)—C_ds = 1.25 for ERA5 and C_ds = 0.35 for ALADIN). Performance metrics (R, HH, NRSME and NBIAS, as described in Section 2.6) were calculated to show agreement with measurements.

Figure 11 shows a decrease in wave model performance in the validation period when compared with the calibration period. This is expected behavior. The decrease in ERA5 wave model performance is small (NRMSE = 0.42, HH = 0.44 in the calibration period to NRMSE = 0.43, HH = 0.45 in the validation period). However, the decrease in performance for the ALADIN forced wave model is significant (NRMSE = 0.44, HH = 0.46 in the calibration period to NRMSE = 0.61, HH = 0.60 in the validation period, a 30% increase in the HH index). Both ERA5 and ALADIN forced numerical models still show an overall underestimation of the significant wave height (NBIAS = −0.19 for ERA5 and NBIAS = −0.09 for ALADIN), but this is strongly skewed by very small significant wave heights (H_s < 0.2 m), which is discussed in detail later on.

Figure 11 also shows that the ERA5 simulated wave peak periods have a higher consistency with the measured wave peak period (R = 0.58, NBIAS = 0.04, NRMSE = 0.18, HH = 0.17) compared to the ALADIN simulated wave peak period (R = 0.46, NBIAS = 0.03, NRMSE = 0.27, HH = 0.26). Nevertheless, the modeled wave peak period is still mediocre for both models. It should be noted that for this analysis, the peak wave periods with low significant wave heights (H_s < 0.25 m) are not considered in this analysis due to their inconsistency (Figure 3).

Empirical CDFs in Figure 12 show again a stronger agreement of the ERA5 forced wave numerical model wave parameters with the measured wave parameters. While CDFs for both the ERA5 and ALADIN show a similar error for significant wave height below the 80th percentile, above the 80th percentile the ERA5 model shows greater consistency. On the other hand, the ALADIN forced model CDF goes all the way to 1.8 as the most extreme predicted value. Compared to the measured CDF, the simulated CDFs wave peak periods are similar, with the ERA5 CDF showing greater agreement.

Figure 9 shows the relationship between the model accuracy performance metrics (HH, NRMSE, and NBIAS) and the simulated significant wave height.

The validation period shows the same pattern similar as for the calibration period (Figure 13), with a great underestimation at the low simulated significant wave heights (NBIAS = −0.12 for the ALADIN forced model and −0.09 for the ERA5 forced model) and great overestimations at high wave heights (NBIAS = 0.63 for the ALADIN forced model and 0.70 for the ERA5 forced model). Extremes in the validation period are more pronounced than for the calibration period. The NRMSE and HH metric are somewhat consistent across the simulated significant wave height range, with the metric showing lower model accuracy at the tails with the increasing HH and NRMSE indices (for the low, H_s < 0.3 m, and high significant have heights, H_s > 1.0 m, respectively).

In summary, the patterns shown in Figure 13 are the same as in Figure 9, with an overall deterioration of the accuracy of the wave model. The ALADIN wave model accuracy reduction is greater than the ERA5 wave model, which is visible with consistently higher NRMSE, HH and NBIAS indices for the ALADIN forced model, except for the most extreme significant wave heights (H_s = 1.2).

As shown in Figure 10, the time series of the simulated wave heights follows the trend of the measured wave heights. As was observed during the calibration period (Figure 10), there is still no obvious delay (Figure 14). The ERA5 forced simulated significant wave heights peaks during the energetic events from October 27 to 3 November are better aligned to the measured wave heights peaks than the ALADIN forced simulated wave heights. The ALADIN model tends to systematically overestimate the maximum wave heights. The peak wave period generally also follows the measured peak wave period, although there are situations where the peak wave period erratically increases to unreasonably high values, especially when the significant wave height decreases (30 October at 22:00 and 2 November at 1:00).

In addition, the simulated wave heights closely follow the wind velocity at the location of the wave buoy itself. Showing that mainly wind driven wave fields are influencing the wave climate in front of Split. The Taylor diagram for the significant wave height time series excerpt shows that the ERA5 forced model has a higher accuracy as compared to the ALADIN forced model. The ERA5 forced model has a higher correlation, lower CRMSD, and closer standard deviation to the measurements than the ALADIN forced model. On the other hand, the Taylor diagram for the peak wave period time series excerpt shows similar accuracy for the ERA5 and ALADIN forced wave models. The ERA5 forced model shows higher correlation, but less accurate standard deviation, with the CRMSD of similar values.

As can be seen in Figure 15, the significant wave heights at the boundary conditions are not able to penetrate deeply into the area in front of Split, where the wave buoy is located, due to the sheltering effect of the islands. Regardless of the direction of the wave or wind on the boundary conditions, the same pattern was observed. Therefore, the simulated wave heights originated predominantly from winds in the area just in front of Split itself.

4. Discussion

The simulated wind velocities provided by the ERA5 and ALADIN reanalysis models are not quite consistent with ground wind observations (Section 3.1). HH and NRMSE are relatively high compared to some similar studies (for example, in Wu, Li [27], the HH metric ranges from 0.17 to 0.38). However, Wu, Li [27] only considered winds above the water surface detected by satellites (scatterometer observation data). In addition, they only considered points at least 40 km from the coast, since scatterometer observations near the coast are less reliable (a limitation that would eliminate all scatterometer points in this paper). Higher HH and NRMSE indices are expected in coastal and inland areas, especially in areas with hills and mountains, such as the Split region considered in this paper [67]. Comparisons of ERA5 reanalysis data and ground measurements conducted in Apulia, Italy (approximately 250 km from the region considered here) showed NRMSE values as high as 0.65 [68]. It was shown by [27,31] that underestimation of wind velocity is common in most wind reanalysis products. However, here ERA5 showed high overprediction (0.52 in Hvar) and medium underprediction (−0.29 in Split) at different wind stations. On the other hand, ALADIN showed low bias at all but one site (0.19 in Resnik). There is no obvious geographical explanation for the biases found in the ERA5 wind data, as underestimation, overestimation, and no bias occurred, at three considered coastal sites Split, Ploce, and Resnik, respectively. This is consistent with research findings that ERA5 wind reanalysis fields should be handled with care when dealing with locations with high topographic variations, in particular [69].

After separate calibration procedures of the wave models forced by ERA5 and ALADIN wind reanalysis data, the optimal white-capping coefficient was chosen as 1.25 and 0.35 for the ERA5 and ALADIN model, respectively. The white-capping coefficient is chosen according to the lowest global HH index. The global minimum HH values for the ERA5 and ALADIN forced wave models were 0.44 and 0.46 for the ERA5 and ALADIN models, respectively. In addition, little effect on the model accuracy was observed by increasing the number of wave directional resolution bins from 36 to 72, potentially improving the diffraction through channels between sheltering islands. A higher white-capping value for the ERA5 forced model was probably due to vastly overestimated Hvar wind velocities (NBIAS = 0.52), which is in the vicinity of the wave buoy location. These calibrated white-capping coefficients are considerably lower than the default SWAN value for the Janssen white-capping formulation [54,55] of 4.5. As discussed in Wu, Li [27], the SWAN model with default settings tends to increase the underprediction of significant wave height. The same case is observed here. Considering the recommendation of Amarouche, Akpınar [56], the white-capping coefficient of 1.0 is closer to the calibrated value in this situation, especially for the calibrated value for the white-capping coefficient using ERA5 reanalysis wind data.

As discussed by Wu, Li [27], the wind error is expected to carry over to the simulated significant wave height data, but the calibration of the white-capping parameter was able to compensate for this error, through different calibrated white-capping coefficient. On the other hand, Wu tested different wind forcing datasets with the same white-capping coefficient (SWAN default at 2.36 × 10⁵ using the Komen formulation [49]), as opposed to this study where different white-capping coefficients were chosen separately for each wind forcing dataset.

There could be several reasons for the similar overall performance of the wave models forced by ERA5 and ALADIN during the calibration period. One possible explanation for this surprising result could be the choice of formulation of the wave growth and dissipation term. The Janssen white-capping formulation for these source terms is widely used in the Mediterranean region and was therefore chosen for this study. On the other hand, the newer ST6 formulation [70] showed poor and oversensitive performance in modeling wave events in narrow fetches without swell in the Norwegian fjord system [36]. In the same study, the Kommen formulation performed best in no-swell conditions. This formulation choice could support and improve the model results when considering systems with limited fetch, especially with the more detailed wind dataset ALADIN. Another possible explanation is that even the more detailed wind forcing grid (ALADIN grid) of 8 km is not yet detailed enough to make a notable difference from the coarser wind data grid (ERA5 grid). Further refinement and downscaling of the upstream global wind reanalysis data may be needed to show a significant difference in the overall accuracy of the wave models. Lastly, the documented restricted inclusion of diffraction in the wave computation native for SWAN [42] could also add to the wave errors measured in the sheltered areas, such as those presented in this work. This error would be present for both the ERA5 and ALADIN forced wave models. A more detailed formulation for the diffraction computation could considerably improve the results of the wave model.

Overall model accuracy using ERA5 and ALADIN wind data shows slightly worse performance metrics than wave models in other research, primarily because other research does not focus on fetch limited regions with complex topography, which are often cited as difficult for wind fields to portray correctly. For example, an equally large difference between simulated and measured wave data was observed by Wu, Li [27], where the HH performance index ranged from 0.21 to 0.46, depending on the wave buoy. Another example can be found in Amarouche, Akpınar [30], which showed accuracy in the range of HH = 0.18–0.4 depending on the observed buoy position location in western Mediterranean Sea, with an average overestimation of 0.12 m. Similarly, Bellotti, Franco [7] showed performance of HH = 0.24–0.27 depending on the buoy on the western Italian coast at mid-latitudes. All of these performance indices are given for wave models forced with ERA5 wind fields and for non-sheltered or non-fetch limited buoy locations.

When the calibrated wave numerical models forced by ERA5 and ALADIN wind reanalysis data were tested during the validation period, the ERA5 model showed a similar performance as during the calibration period (the HH index increased by 2%), but the accuracy of the ALADIN wave model drastically deteriorated (the HH index increased by 30%). When considering the peak wave period, the performance indices are similar to a nearby study at the western Italian coast [7] or Yellow Sea [57], although an open region was considered (R = 0.56, NRMSE = 0.26, HH = 0.25 for [7] and R = 0.58, NRMSE = 0.18 and HH = 0.17 for the ERA5 forced model). Overall, the ERA5 wave numerical wave model showed better accuracy than the ALADIN forced numerical model across the board during the validation period. Additional energy should be put into higher accuracy peak wave period predictions in future studies because of its presence in other coastal engineering applications, such as wave runup formulations [71,72] and storm energy indices [29,30].

An interesting aspect that emerged from the analysis is that both the ALADIN and ERA5 wave models tend to overestimate significant wave heights for high wave heights in both the calibration and validation period. The overestimation itself also increased in the validation period compared to the calibration period (increasing from NBIAS = 0.25 to 0.63 for the ALADIN forced model and from NBIAS = 0.39 to NBIAS = 0.70). The increase in NBIAS for significant wave heights is also observed in studies, such as [27]. Wu, Li [27] argues that this is due to poorer wind accuracy for low wind and high wind velocity situations. To mitigate this error in the high significant wave height range, the calibration procedure could be adjusted to focus more on the higher wave heights with weight functions. However, lower wave height predictions could suffer from this procedure, but in some cases, they are not of primary interest to a coastal engineer.

Nevertheless, overestimation is preferable to underestimation when significant wave heights with high return periods are calculated for the design of coastal structures. At least the designer would not underestimate the structure design parameters, which could lead to a rapid deterioration of the structure. This shows that although the calibration procedure has efficiently adjusted the entire modeled wave time series to better fit the measured significant wave heights to the data, it still has a significant error in the important area of high simulated wave heights for both ERA5 and ALADIN forced wave models.

5. Conclusions

The performance evaluations of the ALADIN and ERA5 forced wave models are not based on the same white-capping parameter of the SWAN model, but the white-capping dissipation parameter is calibrated separately for each forcing wind field. In this way, we applied a specific model parameterization to obtain optimal performance for a given forcing wind field.

The procedure for calibrating the white-capping coefficient specifically for the ERA5 forced wave model is able to bring the accuracy to a similar level as the accuracy of the wave model forced with the ALADIN wind field, although the ALADIN wind field is of considerably finer resolution (0.07° of ALADIN versus 0.25° of ERA5). The similar accuracy of the wave models forced using ERA5 and ALADIN wind fields during the calibration period could be observed with the HH = 0.44 for ERA5 and HH = 0.46 for ALADIN. For the validation period, the accuracy of the ERA5 wave model remained similar (HH = 0.45), but the accuracy of the ALADIN wave model deteriorated with a 30% higher HH index. The accuracy of the peak wave period (HH = 0.18 for the ERA5 forced wave model) is on par with other research, but it still has room for improvement.

In summary, the ALADIN wind field showed a significantly better accuracy than the ERA5 wind field when compared to ground-based wind stations, but this did not translate into significantly improved accuracy of modeled significant wave heights in the fetch-limited basin, when model white-capping coefficient is calibrated separately.

This finding opens the possibility of using ERA5, as a relatively coarse wind reanalysis dataset (0.25° of horizontal resolution) for numerical wave models if the calibration procedure for the white-capping coefficient is performed diligently. The main disadvantage of the forced ERA5 wave model in terms of accuracy is a great overprediction for high wave heights (HH = 0.70 for the ERA5 wind field in the validation period). However, the same pattern was observed for the ALADIN forced wave model (HH = 0.63), which may point to a common problem not relevant to the wind field itself. The authors suspect that the white-capping formulation itself or the diffraction formulation in SWAN are the primary reasons, but this will need to be tested in future work. To mitigate this overestimation using a cruder method, the calibration procedure could be adjusted using a weighting function to focus more on higher wave heights, which are more important for coastal applications. However, predictions for low wave heights using this procedure may be adversely affected, but in some coastal applications, only high wave heights are of interest.

Due to the limitation of the study to test only the Janssen white-capping formulation for the wave generation and wave dissipation source terms in Equation (3), we did not investigate whether other white-capping formulations might be more appropriate for fetch-limited regions. The white-capping formulation of Komen, Hasselmann [52] proved promising for narrow fetch situations [36] and could be a reasonable formulation for wave simulation in basins with complex geometry due to island sheltering, such as the eastern Adriatic coast. Additionally, an investigation into the possibility of using a cruder method with weight functions inside the calibration procedure, to force the calibration procedure to focus on high significant wave heights, could be an interesting solution for coastal engineers in fetch-limited basins.