Evaluating the Effects of Water Circulation on the Modeling of Wave Propagation in the Southern Coast of the Iberian Peninsula

1. Introduction

The simulation of wave propagation is strongly determined by a variety of forcing conditions including wind, wave, bottom friction, bathymetry, sea surface height, and current velocity. Numerical models are powerful tools that can simulate and predict the interactions between these factors dynamically at fine spatial scales and accurately portray wave propagation. This is crucial for coastal engineering, forecasting storm events, analyzing sediment transport, and gaining a deeper understanding of shallow water processes. For forecasting purposes, besides the lateral boundary conditions, many operational wave models include only atmospheric input in their predictions of wave propagation, as modeling both hydrodynamics and wave propagation can be computationally expensive [1,2]. However, the integration of wave and hydrodynamic models has been shown to improve model performance by reducing the bias and error of simulated wave results [3,4]. Consequently, the study of wave-current interaction has been of recent interest and many wave models are commencing now to incorporate hydrodynamic forcing such as current velocity and water levels.

The interaction between waves and currents is complex as the wave spectrum consists of various components that respond differently to currents and thus there is no simple relationship between these two processes [5]. Currents are influenced by multiple variables, including wind stress, tides, and density gradients. At the same time, wind stress can be affected by the sea state, as surface waves modify air-sea momentum transfer [6,7]. The sea state, in turn, depends on wind stress and currents, creating complex feedback mechanisms between waves and hydrodynamics [7]. The interaction between hydrodynamic and wave models can allow for the exchange of information between these complex processes and can improve our understanding of the non-linear interactions between them. More explicitly, the integration between wave and hydrodynamic models allows for the computation of wave-current interaction processes, which include wave radiation stress and enhanced orbital velocities, which further influence water circulation and sediment transport [8]. Hydrodynamics directly impact the wave field by causing the Doppler effect to occur as well as refraction induced by currents [6]. This amount of refraction depends on the current velocity, the incident wave direction compared to the current, and the wave period [9]. On one hand, the Doppler effect can be seen when waves travel with a current (following current), and the observed frequency decreases. On the other hand, when a wave meets a current from the opposite direction (opposing current), the observed frequency increases and the waves appear to be moving faster [10,11]. Longuet-Higgins and Stewart [10] were some of the first researchers to discover changes in wave amplitude due to the non-linear interactions between waves and currents. More recently, several studies have found that an opposing current directly impacts the wave by increasing the wave height, while a following current increases the mean wavelength, which in turn leads to a decrease in wave height [9,12,13]. This change in the spatial distribution of the wave action has been observed by Barnes and Rautenbach [9] who used SWAN to assess the influence of currents in the strong Agulhas Current System off the coast of South Africa. These authors found increases in significant wave height of 20 – 40% with an opposing current and decreases of approximately 20% in the case of a following current. Similarly, Saruwatari et al. [14] observed changes in wave energy of up to ± 60% when implementing SWAN around the Orkney Islands in Scotland where tidal currents exceed 3 m/s. In a weaker current system off the south coast of Ireland, Calvino et al. [2] also found redistribution of the wave action when implementing WAVEWATCH III forced with surface currents. These authors stated that these significant changes in the wave field would not have been captured without hydrodynamic forcing, highlighting the need to include these parameters in wave models for an accurate representation of the wave field.

Several studies have found improvements in wave model results when adding current velocity as a forcing condition of the wave model [2,5,6,9,13]. The magnitude of improvement showed to be dependent on how accurate the modeled current fields are as well as the location of the wave/current regime. For example, Barnes and Rautenbach [9] found only minor improvements of 0.01 to 0.02 m in significant wave height accuracy when forcing their SWAN wave model with a surface current field obtained from the Operational Mercator global ocean model. In another study in the North Atlantic Ocean off the coast of Ireland, the highest change in significant wave height was approximately 10 cm when comparing results from a WAVEWATCH III simulation forced with surface currents and the same simulation without current velocity forcing [2]. Marechal and Ardhuin [5] found that forcing their WAVEWATCH III model with a well-resolved surface current field could accurately capture significant wave height gradients in the Agulhas Current System when comparing their results with satellite altimetry data. In a fully integrated wave-hydrodynamic model (NEMO and WAVEWATCH III) for the Black Sea, Causio et al. [6] found improvements in both the wave field and hydrodynamics when they considered the wave-current interaction in their model. A fully two-way coupled model between an unstructured 2D ocean circulation model Advanced CIRCulation (ADCIRC) and SWAN used depth-averaged currents for the 2D current field and was able to accurately simulate wave propagation in deep and shallow waters in the Gulf of Mexico and Louisiana during four hurricanes [15].

The majority of previous numerical studies have considered only surface currents or the full depth-averaged current velocity field, but none of these aforementioned studies have deeply investigated the depth at which current velocity impacts the wave field. This is especially true for the south coast of Portugal where very few modeled studies on the wave regime exist. Silva et al. [16] implemented the WAVEWATCH III and SWAN models to simulate waves off the Iberian Peninsula and further assess the spatial distribution of wave energy. These authors used the default SWAN settings to run a three-year hindcast and did not consider the effects of currents in their wave model [16]. Rusu et al. [17] also developed a wave forecasting system based on SWAN for western Iberia, but did not force their model with currents. In a nearshore wave-hydrodynamic study in the western margin of the south coast of Portugal, Horta et al. [18] implemented the Steady State Spectral WAVE (STWAVE) model along with field observations to understand how the orientation of the coast impacts nearshore water circulation in different wave conditions. These authors found that the angle at which the waves approach the coastline strongly determines whether the water circulation will be dominated by longshore currents or rip currents [18]. While this study highlights the significance of the effects of waves on water circulation at a specific site along the south coast, the effects of currents on wave propagation along the entire coast have yet to be investigated. To our knowledge, this is the first comprehensive study quantifying the effects of currents and water levels on wave propagation modeling in this region. The spatial variability of the effects of currents associated with the different circulation patterns in the region are also analyzed. Furthermore, the experiments implementing different hydrodynamic forcing in the wave model can be generalized to other regions. This study consists of novel results on how current velocity and water levels impact wave propagation in the south coast of Portugal and sets the stage for future work on how wave propagation influences water circulation by coupling both processes in a fully operational modeling system.

2. Materials and Methods

2.1. Study Area

The study area covers the southern coast of the Iberian Peninsula (southwestern Europe), located in the North Atlantic Ocean (Figure 1). Ocean circulation patterns in this region are complex as the Azores Current separates the oceanic gyre circulation into the Portugal Current in the north and the Canary Current in the south [19]. Furthermore, the entrance to the Mediterranean Sea allows the exchange between two different water masses, which further impacts the regional circulation [20]. The Portugal Current flows southward year-round along the west coast of Iberia [21]. In summer, predominantly northeasterly winds transport water southward along the west coast of Portugal. In autumn and winter the wind direction shifts to southwesterly, reversing the surface circulation and forming the Portugal Coastal Counter Current that flows northward for all layers in the water column from the surface to a depth of ~1500 m [19]. Mesoscale activity is the dominant factor influencing the meteo-oceanographic variability in the region, with variations in coastal alongshore wind stress occurring on a timescale of several tens of days [22].

The region is furthermore part of the Canary Current Upwelling System (CCUS), which extends from 43°N to 10°N along the west coast of Iberia to the northwest coast of Africa. Upwelling off the west coast of Portugal occurs in spring and summer when surface water along the coastline is transported offshore due to the combined effects of the northeasterly winds and the Coriolis force. [21]. Off the southern coast of Portugal, the Portugal Coastal Countercurrent transports warm surface water westward in winter [23]. Once these waters reach the furthest western point along the south coast of Portugal at Cape São Vicente, they occasionally rotate clockwise and are further transported northward along the west coast [23]. The northward transportation of surface water along the west coast produces downwelling events [22].

The wave regime in the region is characterized by dominant more energetic waves coming from W-SW and less energetic E-SE waves [24]. In the south coast of the Algarve the average significant wave height offshore is 0.9 m. Storm wave events are marked by waves reaching a significant wave height of more than 3 m and an associated peak period of 10 – 12 s [25]. Storms occur typically with waves coming from W-SW or E-SE, the former being more frequent and energetic [24]. Waves are much more energetic off the west coast, with average significant wave heights ranging between 2 and 3 m and peak periods between 9 and 11 seconds [26].

The coastline along the south coast of Portugal is highly variable as Cape Santa Maria at 36.97°N, 7.90°W splits the south coast into western and eastern margins. The western margin extends from Cape São Vicente to Cape Santa Maria and the eastern margin extends from Cape Santa Maria to the Guadiana River at the border between Portugal and Spain [27] (Figure 1). The western margin consists of a cliff-dominated, rocky coastline with pocket beaches and is fully exposed to W-SW wave conditions [18]. The eastern margin consists of mainly sandy beaches in front of the Ria Formosa barrier system and lagoon and is fully exposed to E-SE waves [27].

2.2. SWAN

The Simulating WAves Nearshore (SWAN) model, developed at Delft University of Technology, is a third-generation wave model that simulates waves in coastal and inland waters. SWAN solves the full spectral balance equation, incorporating source and sink terms. Since SWAN is a third-generation model, the wave spectrum can propagate freely without any constraints on its evolution. All processes of wave generation, dissipation, and wave-wave interaction are represented [28]. The spectral balance equation is shown below [28].

$${\displaystyle \frac{\partial N}{\partial t}}+{\nabla }_{\overrightarrow{x}}\cdot \left[\left(\overrightarrow{c_g}+\overrightarrow{u}\right)N\right]+{\displaystyle \frac{\partial c_{\mathsf{\sigma }}N}{\partial \mathsf{\sigma }}}+{\displaystyle \frac{\partial c_{\mathsf{\theta }}N}{\partial \mathsf{\theta }}}={\displaystyle \frac{S_{\mathrm{tot}}}{\mathsf{\sigma }}}$$

(1)

The terms on the left describe: 1) the action density rate of change $(N$) at a single point in space over time $t$, which is further equivalent to E/ $\mathsf{\sigma }$ where E is the energy density and $\mathsf{\sigma }$ is the relative frequency, 2) propagation of wave energy in two-dimensional geographical $\overrightarrow{x}$-space where $\overrightarrow{c_g}$ is the group velocity and $\overrightarrow{u}$is the ambient current velocity vector, 3) shifts in frequency due to currents and changes in depth, and 4) depth- and current-induced refraction. The term $S_{\mathrm{tot}}$ on the right represents the sources and sinks of wave energy which include the total wave growth by wind, the redistribution of energy from nonlinear interactions between different wave components, and wave decay caused by bottom friction, whitecapping, and depth-induced wave breaking [29]. In deep water the most influential processes impacting wind-wave growth are the atmospheric input from wind to waves, wave dissipation, and the non-linear energy transfer between wave components [30]. Processes such as bottom friction and triad interactions become more important in shallow waters with finite depth (up to a few hundred meters) [30].

2.3. Hydrodynamic Model: SOMA

The current velocity and water level used to force the wave model is produced by the Algarve Operational Modeling and Monitoring System (SOMA) [31] (www.soma.ualg.pt), which is based on the hydrodynamic model MOHID [32]. The MOHID water modeling system [33] is a three-dimensional water modeling system programmed in ANSI FORTRAN 95 that simulates physical processes in an object-oriented approach. The design of SOMA consists of two nested grids of increasing resolution. The coarser level 1 grid covers the SW Iberian region with a horizontal resolution of 2 km, with the finer level 2 grid nested with a resolution of 1 km, allowing a finer resolution of processes near the coast. The model bathymetry was created with data from the European Marine Observation and Data Network (EMODNET). The discretization of the water depth consists of 50 unevenly spaced layers of Cartesian coordinates. SOMA has been previously calibrated and validated with data from tide gauges, mooring buoys and vertical profiles [31]. This model has been fully operational since July 2019, producing four-day forecasts of current velocity, water levels, temperature, and salinity for both level 1 and level 2 grids [34]. For a complete description on the set up and validation of SOMA please refer to Janeiro et al. [31].

2.4. Standalone SWAN Model Setup

2.4.1. Model Setup

The standalone SWAN model shared the same spatial domain and resolution as the level 1 SOMA grid with 195 cells in the x-direction and 115 cells in the y-direction, spanning from 38.35°N to 36.07°N and -10.45°E to -6.55°E (Figure 2). The model uses the same bathymetric dataset from EMODNET. SWAN Cycle III version 41.31A with MPI (Message Passing Interface) parallel computing over 8 cores was used for all wave simulations. Although there is a newer version of SWAN available (version 41.51 AB), this version was selected due to its compatibility with MPI for Windows. Simulations were run in two-dimensional (2D) nonstationary mode with spherical coordinates. The Nautical convention was used for wind and wave direction. The wave spectra are computed over the entire 360° directional range, using 72 bins with a resolution of 5° each. The grid resolution in frequency-space was set to 35, ranging from a minimum frequency of 0.0345 Hz to a maximum frequency of 1 Hz. Since wave observations in the breaking zone is very limited and insufficient, authors opted to use the default values for depth-induced wave breaking, which consists of a proportionality coefficient of the rate of dissipation set to 1 and the ratio of maximum individual wave height over depth (breaker index) set to 0.73, in agreement with previous studies [17]. Bottom friction, triads, and quadruplets were also activated. All simulations were run with a time step of 10 minutes, producing output every hour.

2.4.2. Wind Input

The SWAN and MOHID models are forced with wind input from the SKIRON Atmospheric Modelling and Forecasting Group in Athens [35]. The wind data from SKIRON consists of wind velocity in the x and y directions at a height of 10 m at hourly intervals. The wind output covers well over the entire spatial domain of the present study and has a resolution of 0.05° (approximately 5 km).

2.4.3. Wave Conditions at the Boundary

Wave conditions at the boundary are provided from the Atlantic-Iberian Biscay Irish-Ocean Wave Analysis and Forecast product from the Copernicus Marine Iberia-Biscay-Ireland Monitoring and Forecasting Center (IBI-MFC). The Copernicus Marine IBI-MFC assimilates data from satellite altimeters to include near real-time wave observation data in their models [36]. The dataset consists of a high-resolution wave analysis and forecast product on a 0.0278° x 0.0278° spatial grid covering the Iberia-Biscay-Ireland area from 26°N to 56°N and -19°E to 5°E [37]. Hourly instances of significant wave height, peak period, and mean wave direction from this dataset were imposed along the boundary. Since the spatial resolution does not match perfectly the spatial resolution in this work, segments were set up along the boundary with the same spatial resolution from the IBI-MFC dataset. Each segment represents one grid point from the IBI-MFC dataset. The sensitivity of the model to directional spreading $co{s}^m\left(\theta \right)$ was tested by running the model with different coefficients of directional spreading at the boundary and assessing whether to measure as the directional standard deviation in degrees or as the power of $m$. The SWAN manual recommends using default values of 30 when measuring in degrees or 2 when measuring in power. As the southern and western boundaries are exposed to more open swell and the eastern boundary is more sheltered, the sensitivity analysis included modifying the coefficients based on boundary side. While not shown here, this was previously tested and model configurations that used directional spreading power of $m$ set to 2 around the lateral boundaries yielded overall results most similar to observations, particularly near the coast. Results of this sensitivity test are shown in the appendix.

2.4.4. Standalone SWAN Calibration and Validation

The standalone SWAN model was calibrated by testing the model with a set of formulas that represent the physical parameters influencing wave propagation. This was done for a monthlong simulation in February 2024 when waves were predominantly westerly. Several source term packages are available in SWAN, each implementing different formulas for wind growth and whitecapping. The newest source term package is called ST6, which consists of new formulas for wind and whitecapping and considers the influence of swell dissipation, wind scaling and other processes [38]. This ST6 package has been recommended by Edwards et al. [39] to use in regions exposed to swells generated in the Atlantic, Pacific, Indian, and Southern Ocean and thus it was chosen to calibrate the model with different formulas and coefficients within this package. The SWAN manual provides several ST6 physics settings to use for model calibration [40]. The various calibration settings of the ST6 physics formulas consist of modification of the dissipation coefficients $a_1$ and $a_2$ and power coefficients $L$ and $M$, normalization of the wind drag formula with the commands UP $\left(E_T\left(f\right)\right)$ or DOWN $\left(E\right(f\left)\right)$, selection of a formula for wind drag (HWANG, FAN, or ECMWF), implementation of vector (VECTAU) or scalar (SCATAU) to calculate stress, and selection of the true value for wind (TRUE10) or a scaling factor (U10PROXY) [30]. The term AGROW activates the wave growth term from Cavaleri and Malanotte [41]. The model cases used for calibration are shown in Table 1. The terms for swell were also modified with the options Ardhuin, which refers to the non-breaking dissipation from Ardhuin et al. [42], which uses a coefficient related to the laminar atmospheric boundary layer and Zieger [43], which uses a non-dimensional proportionality coefficient. Calibration was done for a monthlong simulation in February 2024, when waves coming from W-SW were observed.

Results of each case were compared with wave measurements at three mooring buoy locations with available data during the calibration period: Faro Ocean Buoy, Cadiz Buoy, and Sines Buoy (Figure 2). Wave data from these three buoys were obtained from the CMEMS In Situ Thematic Assembly Center (In Situ TAC) product Global Ocean In-Situ Near Real-Time Observations [44]. Measurements from Faro Ocean Buoy and Sines Buoy are provided by Instituto Hidrográfico (IH) whereas measurements from Cadiz Buoy are provided by Puertos del Estado. Data from two buoys closer to shore (Faro Coast and Armona Buoy) were acquired after the calibration period and only contain data for the year 2023. The 2023 dataset at Faro Coast Buoy was purchased from IH and the Armona Buoy dataset was provided by CIMA at the University of Algarve. Buoy locations and their depths are described in Table 2.

All five buoys used for model calibration and validation are equipped with high-precision sensors that accurately represent wave parameters such as significant wave height, mean wave direction, and peak period. Calibration was achieved by obtaining hourly measurements of significant wave height (Hs), mean wave direction (Dir), and peak period (Tp) from the three mooring buoys that were available during the calibration period and comparing them with results from the various simulations of altered physics formulas.

To quantify model performance several metrics were computed: root mean square error ($RMSE$), $\mathrm{BIAS}$, Pearson’s Correlation Coefficient (R), and Model Skill Score ($MSS$). $RMSE$ measures the magnitude of error between model predictions and observations and is expressed as:

$$RMSE=\sqrt{{\displaystyle \frac{{\sum }_{i=1}^n{\left(x_{\mathrm{model}}-x_{\mathrm{obs}}\right)}^2}{n}}}$$

(2)

where $x_{\mathrm{model}}$ and $x_{\mathrm{obs}}$ represent modeled data and observed data from the mooring buoys, respectively, and $n$ is the number of observations [45]. Similarly, the BIAS represents an absolute measure of how much model results differ from observations and is used to determine whether the model tends to overestimate (positive value) or underestimate (negative value) observations. It is expressed as:

$$\mathrm{BIAS}={\displaystyle \frac{1}{n}}{\sum }_{i=1}^n\left(x_{\mathrm{model}}-x_{\mathrm{obs}}\right)$$

(3)

To quantify the correlation between predicted values and observations, Pearson’s Correlation Coefficient was used. It is expressed with the following formula:

$$\mathrm{r}={\displaystyle \frac{\sum \left({\mathrm{x}}_{\mathrm{obs},\mathrm{i}}-\overline{{\mathrm{x}}_{\mathrm{obs}}}\right)\left({\mathrm{x}}_{\mathrm{model},\mathrm{i}}-\overline{{\mathrm{x}}_{\mathrm{model}}}\right)}{\sqrt{\sum {\left({\mathrm{x}}_{\mathrm{obs},\mathrm{i}}-\overline{{\mathrm{x}}_{\mathrm{obs}}}\right)}^2\sum {\left({\mathrm{x}}_{\mathrm{model},\mathrm{i}}-\overline{{\mathrm{x}}_{\mathrm{model}}}\right)}^2}}}$$

(4)

Pearson’s Correlation Coefficient measures the linearity between observed and predicted values, ranging from -1 (perfect negative correlation) to 1 (perfect positive correlation) with 0 indicating no correlation. It should be noted that this statistics measure can be influenced by outliers and may not reflect the magnitude of agreement between observations and predictions [45] and thus it should be interpreted alongside other statistics metrics.

Willmott’s Model Skill Score (MSS) provides a more quantitative measurement of model accuracy by measuring the degree to which model predictions are free from errors relative to the observed data [46]. It overcomes the sensitivity of simple correlation statistics by normalizing errors relative to the variability in the observations [45]. It is represented by the following formula:

$$MSS=1-{\displaystyle \frac{{\sum }_{i=1}^n{\left(x_{model}-x_{obs}\right)}^2}{{\sum }_{i=1}^n{\left(\left|x_{model}-{\overline{x}}_{obs}\right|+\left|x_{obs}-{\overline{x}}_{obs}\right|\right)}^2}}$$

(5)

A MSS with a value of 1 represents perfect agreement between modeled data and observed data, whereas a value of 0 indicates no agreement between the two.

The four error metric formulas described above were applied to each simulation with varying physics calibration settings and observed data from each of the three buoys. The model result that yielded the lowest RMSE, lowest magnitude of BIAS, highest correlation coefficient, and highest MSS was selected as the physics setting for future model runs of this study.

After calibration, validation was performed by running a monthlong simulation in November 2023. This time period covers a storm with Hs higher than 2.5 m and waves predominantly coming from W-SW. Validation at this time period also allows to assess model performance closer to shore, by comparing results at Faro Coast Buoy. Error metrics of RMSE, BIAS, MSS, and R were computed for model results and the four mooring buoys: Faro Ocean Buoy, Faro Coast Buoy, Cadiz Buoy, and Sines Buoy.

2.4.5. Experimental Design: Hydrodynamic Forcing on SWAN

The coupling between hydrodynamic and wave models can be implemented in two ways: offline coupling, in which the output of one model serves as input for another, and online coupling, where both models exchange information dynamically during the simulation. While online coupling provides a more complete interaction between wave and hydrodynamic processes, it requires significantly higher computational resources and an external mechanism to facilitate communication between models [6]. In contrast, offline coupling is computationally more efficient and allows for a controlled analysis of how one process influences the other. In an offline coupling approach, a one-way coupling system is achieved by running one model and using the results of that model as boundary conditions for the other model.

Various outputs from SOMA were used to test the model sensitivity to hydrodynamic forcing. The calibrated and validated standalone SWAN set-up without any current velocity or water level forcing from SOMA was used as the control run (Run 1). To test the hydrodynamic influence on wave propagation in the region several simulations were run with different forcing output from SOMA as indicated in Table 3. All outputs from SOMA that were used as initial conditions for the wave model during these experiments were time and space varying, covering the same spatial grid and timeframe of each simulation.

The surface current field consists of northward and eastward velocity components of the upper most layer of the vertical coordinate in SOMA’s grid. The base value of this first layer is 0.494 m to hydrographic zero. Therefore, the depth of this first layer varies according to the variation in the tides such that it will be 0.494 m + sea surface height (SSH). Since SWAN only takes two-dimensional input, the current velocities referring to various depths needed to be weighted-averaged. To compute the depth-averaged velocity for the first 10 m, first 20 m, and entire water column, the depth-integrated velocity of both eastward and northward velocity components was computed. A weighted average was used considering the thickness of the cells at each layer for each (x, y) point.

To evaluate the performance of the hydrodynamically forced model in various environmental conditions, three different time periods were selected to run the simulations with: 1) Westerly mean wave direction, 2) Easterly mean wave direction, and 3) a limited period with available data to validate the model in shallow water conditions (experiments 3, 4, and 5) in Table 4. The first two conditions were selected by analyzing the significant wave height and mean wave direction from each of the four buoys (Faro Ocean, Faro Coast, Cadiz, and Sines). The first timeframe corresponds to two consecutive storms in October/November 2023 when waves were predominantly coming from W-SW and the significant wave height was at least 2.5 m. The second timeframe corresponds to a storm in February 2023 when waves were predominantly coming from E-SE and the significant wave height was more than 2.5 m. The third scenario was selected based on the newly acquired buoy data at Armona Buoy, which only has data from March to July 2023. The three previously mentioned scenarios are summarized in Table 4, which further includes the simulation timeframes selected for the standalone SWAN calibration and validation.

To analyze and compare the results of these models with buoy observations, the same metrics described above of RMSE, BIAS, Model Skill Score, and R were computed. These statistics metrics were also computed for the CMEMS model that was used to force the wave conditions at the boundary to compare model performance with a well-established operational wave model.

3. Results and Discussion

3.1. Standalone SWAN Calibration

To calibrate the model for the different physics settings, various simulations were run altering the different formulas for the physics for a monthlong simulation. Time series of Hs for Faro Ocean Buoy, Cadiz Buoy, and Sines Buoy are shown in Figure 3, highlighting the observed Hs of each buoy on the left axis and the relative error between modeled results and observations (the envelope of seven runs) on the right axis. Time series showing predicted values and observations at each buoy location are shown in the appendix.

Overall, Figure 3 demonstrates that modeled Hs only slightly varied across the different ST6 physics calibration settings, demonstrating that the model was not sensitive to changes made in the formulas for wind input, quadruplet interactions, and whitecapping dissipation. The largest difference between modeled results occurred at Cadiz Buoy (Figure 3b) on 2024-02-09 directly after the highest peak of significant wave height with differences in relative error ranging from -40% to -20%. Another large difference between modeled results can be seen at Faro Ocean Buoy (Figure 3a) in the beginning of the simulation on 2024-02-03 when some modeled results underpredicted Hs right after the peak of the first storm, whereas other modeled results overpredicted Hs. At Faro Ocean Buoy, the highest peak of Hs occurred on 2024-02-08 at 23:00 with a value of 5.74 m. The highest peak of Hs for Cadiz Buoy occurred on 2023-02-09 at 4:00 with a value of 6.21 m. For both of these peaks that occurred in the beginning of the simulation, all models underpredicted Hs, but then tended to overpredict Hs for the remainder of the simulation. At Sines Buoy, the highest peak of Hs occurred on 2024-02-24 at 9:00 when observed Hs reached 7.35 m and modeled results only slightly underpredicted this value. Modeled results at Sines Buoy achieved Hs values most similar to observations compared to the results at the other buoy stations.

For a more quantitative analysis on each model’s agreement with buoy observations, error metrics computed at each buoy site are shown in Table 5. The values in bold indicate the lowest (RMSE or BIAS) or highest (MSS or R) among all cases at each site. The physics formulas implemented for each case are referenced in Table 1. The metrics for mean wave direction and peak period at each site are in the appendix.

Small differences in the error metrics at each of the buoy locations were found when comparing the results from the various physics settings. In general, values of RMSE and BIAS were quite low, with values of RMSE varying between 0.361 and 0.409 and the BIAS varying between -0.001 and 0.172. Case 6 was selected as the best setting for the model as it had the highest MSS (0.963) and Pearson’s Correlation Coefficient (0.928) and lowest RMSE (0.366) and BIAS (-0.001) at Faro Ocean Buoy. This case also resulted in the highest value in MSS and Pearson’s Correlation Coefficient at Cadiz buoy. The RMSE at this location was slightly higher compared to other simulations, but this can be influenced by the initial conditions at the boundary, since Cadiz Buoy is located very close to the boundary. The error metrics at Sines Buoy did not vary much between simulations. Case 6 was chosen as the results at Faro Ocean Buoy exceeded all other simulations and this location is a strong indicator of how the model propagates wave energy further in the domain. The BIAS value of nearly 0 and high MSS of 0.963 indicate the model’s adequacy for simulating wave propagation well into the computational domain. The calibrated set of parameters is thus the ST6 formula with dissipation coefficients 4.7E-7 and 6.6E-6 and power coefficients of 4.0 and 4.0, normalization with the command UP, wind drag formula of HWANG, vector calculation for the stress calculation with the command VECTAU, windscaling equal to 28 with the command U10PROXY 28.0, and linear growth with the command AGROW. Furthermore, swell is activated with the command ZIEGER and set to the non-dimensional proportionality coefficient 0.00025, which is default in SWAN [40]. The wind drag formula of HWANG was also found to be the best formula within the ST6 package in a recent study that used SWAN to assess wave energy off the Iberian Peninsula [47].

The results of model calibration demonstrated that the modeled wave height was not very sensitive to the selected physics parameters and maximum differences in RMSE at each of the three buoys was approximately 0.13 m. This behavior was also observed by Zhang et al. [38] and Aydoğan & Ayat [48] who found improvements to the SWAN model when calibrating with ST6 parameters, but small differences between the various ST6 calibration settings.

3.2. Standalone SWAN Validation

For model validation, a new simulation was run with the calibrated physics parameters previously found. Scatter plots showing the correlation between buoy observations and model results during the validation period of November 2023 when the dominant wave direction was W-SW are shown in Figure 4.

Figure 4 demonstrates a very strong correlation between observations and modeled results at all four buoy sites. The strongest correlation is at Faro Coast Buoy, which is located very close to shore. The results at Faro Ocean Buoy and Cadiz are slightly more spread out. The SWAN model tended to overestimate Hs at these two sites, which are located in the southern part of the domain, offshore. At Sines Buoy, the highest values of Hs were observed, with values reaching close to 10 m. The SWAN model underpredicted Hs at this high peak. For a more precise quantitative analysis on how the model results differ from buoy observations, Table 6 summarizes the error metrics of significant wave height, mean wave direction, and peak period.

Results of modeled Hs showed strong agreement with observational data across all four buoy locations for the validation period of November 2023. This agreement is reflected in the low values of RMSE and BIAS, as well as the high values of MSS and R, demonstrating the model’s ability to capture both the temporal variability and magnitude of these wave parameters. Among the four buoy locations, Faro Coast Buoy, which is located closer to shore, exhibited the best agreement with a MSS of 0.973 and RMSE of 0.156 m for Hs, demonstrating the model’s ability in accurately capturing wave propagation far into the computational domain. Faro Coast Buoy did not have data available for mean wave direction, so this parameter could not be compared with model results. Modeled mean wave direction aligned well with observations at all buoy locations, attesting the model’s adequacy in accurately simulating mean wave direction. Results of peak period yielded slightly higher errors, but still showed good agreement with observations. These results demonstrate the model’s adequacy in simulating and predicting wave parameters within the model domain and signify that it can be used to accurately represent the wave field for the south coast of Portugal. It is still crucial to constantly verify model results with observations [49] and thus all future runs of this study are validated with in-situ buoy observations.

3.3. Effects of Hydrdoynamic Forcing on Wave Propagation

3.3.1. Waves Coming from W-SW

To assess the influence of the various hydrodynamic inputs on wave propagation during a storm with predominantly W-SW waves, a time series of significant wave height is exhibited in Figure 5. Buoy observations are shown on the left axis, while the range in relative error between all modeled results (run IDs from Table 3) and buoy observations are plotted on the right. This allows for an analysis of model ability to accurately capture wave propagation throughout the simulation, and to assess the variability between the model runs that implemented various hydrodynamic forcing.

In general, it can be stated that when Hs was high (above ~3 m), all models underpredicted Hs (Figure 5). Off the south coast, the peak of the storm occurred on 2023-10-22 as shown by the results of Faro Ocean Buoy, Faro Coast Buoy, and Cadiz Buoy (Figure 5a, b, c). At Faro Ocean Buoy and Cadiz Buoy (Figure 5a, c), all models underpredicted Hs for these peaks, but then tended to overpredict Hs for the remainder of the simulation. However, at Faro Coast Buoy (Figure 5b), all models overpredicted Hs during the highest peak, with some models overpredicting Hs by nearly 1 m with a relative error of approximately 20%. At Cadiz Buoy, a relative error of 100% is observed during low wave conditions on 2023-10-29. At this time, measured Hs was 0.5 m, while model predictions overestimated it to 1 m. Although this error seems high, the absolute difference is small and within a reasonable range. Faro Coast Buoy and Sines Buoy yielded model results most similar to observations. At Faro Coast Buoy, Hs stayed relatively small throughout the simulation with values ranging between 0 and 2 m. Hs reached the highest off the west coast, at Sines Buoy with values reaching nearly 10 m on 2023-11-05. The models were not able to capture this extreme wave height. However, they matched the observations quite well during the simulation and tended to overestimate Hs when Hs was low and underestimate Hs for higher values. In terms of the variance between model results, Figure 5 shows more variability at Faro Coast (Figure 5b) and Cadiz (Figure 5c) towards the end of the simulation when Hs was low. The largest variability between models in relative error at these two locations was 37.71% and 29.84%, respectively. The smallest variability between models is seen at Sines Buoy (Figure 5d) with a maximum variance between the relative error of models of 7.51%.

To more precisely quantify the error between each specific model run and observations, Table 7 presents the metrics of RMSE, BIAS, MSS, and R. The run IDs are referenced in Table 3. The best metric value (lowest value of RMSE and BIAS, highest value of MSS and R) amongst all simulations is highlighted in bold.

Table 7 demonstrates small variability in the statistics metrics between uncoupled runs and runs that implemented the different hydrodynamic forcing from SOMA. RMSE and BIAS values stayed quite low ranging from 0.242 m to 0.594 m and 0.005 to 0.502, respectively. MSS and Pearson’s Correlation Coefficient were high ranging from 0.862 to 0.982 and 0.923 and 0.968, respectively. When comparing the SWAN runs that implemented varying hydrodynamic forcing from SOMA, Run 3 (SOMA Surface Currents and Water Level) performed better than all other simulations, by achieving the highest number of “best” metric values for all buoys. At Faro Coast Buoy, the run with SOMA water level only (Run 2) performed very well with an RMSE of 0.242 m and MSS of 0.958. This result is consistent with Wang and Elahi’s [13] finding that water levels play a more significant role in the modulation of significant wave height in shallower waters, whereas current velocity has more influence on wave propagation in the deep ocean. In comparison to the control run of the study without any hydrodynamic forcing (Run 1) all simulations performed better, implying that adding hydrodynamic forcing improves results. At Faro Ocean Buoy, the RMSE was reduced by 0.054 m or approximately 10% when surface currents and water level were implemented in the model (Run 3). The RMSE was reduced by 0.004 m at Faro Coast Buoy with the case of SOMA water level only. At this buoy location, the other runs produced less accurate results. At Cadiz Buoy the RMSE was reduced by 0.164 m (~28%) in the run with SOMA surface currents and water level. In Sines, this same run reduced the RMSE by 0.029 m. Overall, the metrics presented in Table 7 reveal improvements to the model with the addition of hydrodynamic forcing, particularly in the case of surface currents and water level. Table 7 shows better performance from the CMEMS model at all buoy locations for the simulation with westerly waves. This can be due to the fact that CMEMS assimilates wave data from satellite results, leading to a more accurate representation of the model’s wave field [36].

For a visual representation of how wave propagation in the computational domain changed due to the addition of hydrodynamic forcing, Figure 6 shows spatial maps of the difference in Hs between each simulation and the simulation without hydrodynamic forcing during the westerly storm. The 2D maps represent the average difference in Hs during the westerly storm period with overlaying velocity vectors representing the average current velocity direction as well as the mean wave direction during the storm. For a complete evaluation of the magnitude of current velocity over time, time series of current velocity at four locations within the domain are shown in Figure 7.

The magnitude of change in Hs for each simulation as shown in Figure 6 was very similar for the simulations forced with SOMA surface currents, SOMA 10 m depth-averaged velocity, and SOMA 20 m depth-averaged velocity, whereas the change in Hs for the simulation forced with depth-averaged velocity for the full water column was much smaller. The change in Hs was highest for the simulation forced with surface currents (Figure 6a) and then reduced slightly as the depth increases, as shown from Figure 6b to Figure 6c to Figure 6d. This is consistent with the time series of velocity as shown in Figure 7, as the magnitude of the depth-averaged velocity of the entire water column was always smaller compared to velocities more towards the surface. The magnitude of change in Hs was higher off the south coast compared to the change in Hs off the west coast. Figure 6a, b, and c all show decreases in Hs of approximately 0.4 m when the current velocity direction and mean wave direction were aligned. Close to the shoreline on the southern coast, increases in Hs were observed as the average current velocity direction was easterly and the mean wave direction was westerly. Along the west coast, decreases in Hs were also observed due to the southerly direction of the waves. The time series of current velocity in Figure 7 show higher velocity values at Faro Ocean, Faro Coast, and Cadiz compared to Sines, which further justify why changes in Hs are higher off the south coast compared to changes in Hs off the west coast. The magnitude of velocity was also higher in the southeast domain with values reaching close to 1 m/s as shown in the time series plot at Cadiz Buoy in Figure 7c. This further corresponded to a decrease in Hs of approximately 20 cm for all simulations forced with surface velocity as well as depth-averaged velocity up to 10 m and 20 m. The enhancement of changes in significant wave height due to stronger current velocity fields has been detected in other coupled wave-hydrodynamic models [2,9]. The increases in Hs closer to shore as shown in Figure 6 can also be attributed to the implementation of the dynamic SOMA water levels in the runs, underscoring the significance that varying water levels have on wave propagation in nearshore areas. Time series of water levels (Figure A 4) at Faro Ocean Buoy, Faro Coast Buoy, Cadiz Buoy, and Sines Buoy are shown in the appendix.

3.3.2. Waves Coming from E-SE

For the simulation with waves coming from E-SE (February 2023), Figure 8 shows a time series of significant wave height with buoy observations on the left axis and the relative error percentage between the envelope of all modeled results from Table 3 and observations on the right axis.

Similar to the simulation with westerly waves, Figure 8 shows an underestimation of Hs when the observed values were high (above 3 m) and an overprediction for lower Hs values. All models show a consistent slight underprediction at Faro Ocean Buoy and Faro Coast Buoy. The models were able to more accurately capture wave propagation at Sines Buoy, and tended to slightly overestimate Hs. It should be noted that in-situ data from Cadiz Buoy was only available for three days of this simulation (2023-02-19 – 2023-02-22) and thus was not included in this analysis. To more precisely quantify how each model differed from observations, the statistics metrics are summarized in Table 8.

The SWAN simulations performed very well in the simulation with easterly waves and achieved very low RMSE and BIAS values and very high MSS and R values at all locations. RMSE values ranged from 0.288 m and 0.352 m, while BIAS values ranged from -0.195 m to 0.170 m. The simulations that included current velocity up to a specific depth of the water column and water level forcing from SOMA (Runs 3 – 5) outperformed the run without any hydrodynamic input (Run 1). These model simulations reduced RMSE by approximately 0.2 m compared to the run without any hydrodynamic input for all buoy locations. Runs 3 – 5 achieved the highest MSS value at Faro Coast Buoy with a value of 0.974. Run 4 (SOMA depth-averaged velocity up to 10 m and water level) resulted in the lowest RMSE values at both Faro Oceanic Buoy and Faro Coast Buoy. Run 4 produced the highest number of “best” metric values at all locations, which indicates that the run with the depth-averaged velocity up to the first 10 m of the water column produced the most accurate results for the easterly wave simulation. The simulations performed better than CMEMS at all locations with the exception of Sines Buoy, despite the fact that CMEMS assimilates wave observations and our simulations do not.

To further analyze how current velocity and water levels impacted Hs over the entire spatial domain, Figure 9 shows spatial maps of the difference in Hs between each simulation and the simulation without hydrodynamic forcing during the easterly storm. Time series of current velocity at the four buoy locations are shown in Figure 7.

The simulation with waves coming from E-SE resulted in higher changes in Hs compared to the storm with waves coming from W-SW. The highest average change in Hs was observed on the southeast coast, where the average Hs during the storm increased by 0.623 m or 209% when comparing the run with no hydrodynamic input and all SWAN simulations forced with SOMA hydrodynamics. The increase in Hs along this shoreline is attributed to varying water levels from SOMA, which are shown in the appendix (Figure A 5). Furthermore, Figure 10b demonstrates that velocity magnitude closer to the coast was much higher compared to offshore. There was not much variability between the changes in Hs for the simulations forced with the varying depth-averaged velocities for the E-SE storm. All 2D spatial maps show significant increases in Hs at the coast, which highlight the impact that varying water levels have on Hs closer to the coastline. It is noted that in all spatial maps of Figure 9 where the current flowed opposite of the easterly mean wave direction, increases in Hs were observed. When the currents and waves flowed easterly, decreases in Hs were observed. In general, velocity magnitude was much smaller for the easterly simulation compared to velocity magnitude during the westerly simulation. However, the same effect still occurred where waves and currents in opposing directions resulted in increases in Hs and following currents resulted in decreases in Hs. The magnitude of these changes in Hs ranged from approximately 10-20 cm offshore to more than 0.6 m near the coast. Areas with larger magnitudes of velocity also resulted in greater changes in Hs. Opposing currents occurred more frequently off the southwest coast and following currents occurred predominantly off the south coast and west coast.

3.3.3. Wave Propagation in Shallow Waters

To evaluate how waves propagate further into the domain up to nearshore areas, Figure 11 shows the time series result at Armona Buoy for the simulation when data at this buoy was available.

Figure 11 shows a slight underprediction in Hs during the highest peaks of Hs on 2023-04-18, 2023-04-19, and 2023-04-20 and a slight overprediction of Hs for the lowest peaks of Hs on these days. There were several time instances during the simulation where the models and observations were in perfect agreement. The variability between model runs was quite small, with the variability in relative error ranging from approximately 1% to 5% throughout the simulation. This variability became larger towards the end of the simulation on 2023-04-22 when Hs started to decrease and the variability in relative error between models became more spread out, with differences between errors ranging around 30%. All models overpredicted Hs after this decrease in Hs on 2023-04-22. The statistics metrics for all buoy stations with available data during the simulation time are shown in Table 9.

Table 9 shows very good metric scores amongst all simulations, especially at Armona Buoy, which is a very good indicator of how waves propagate further into the domain up to nearshore areas. The RMSE did not vary much between simulations at Armona Buoy. At Faro Coast Buoy, the largest reduction in RMSE when comparing with the run without any hydrodynamic input was for the simulation with the SOMA depth-averaged currents up 10 m of the water column and water level with a reduction of 0.015 m. At Cadiz Buoy the simulation with surface currents and water level (Run 3) reduced the RMSE by 0.008 m in comparison to the run without any hydrodynamic input (Run 1). All simulations performed better than the CMEMS model. However, CMEMS achieved a slightly lower BIAS value at Armona Buoy. The simulation that received the highest number of the “best” error metrics was the run with SOMA surface currents and water level (Run 3).

3.3.4. Main Findings of the Hydrodynamically Forced Simulations

The results show overall improvements in wave model accuracy when forced with current velocity and water level. In the storm with waves coming from W-SW, reduction in RMSE ranged from 1.6% (Faro Coast Buoy) to 27.6% (Cadiz Buoy) and reduction in BIAS ranged from 12.5% (Faro Ocean Buoy) to 80% (Faro Coast Buoy) when comparing Hs results to the run with no hydrodynamic input (Run 1). The largest improvement was observed at Cadiz Buoy where RMSE was reduced by 16.4 cm and BIAS was reduced by 17.9 cm for the run with surface currents and water level. The magnitude of improvement at Faro Coast Buoy was smaller, with a reduction in RMSE of only 0.4 cm for the run with surface currents and water level (Run 3) and a reduction in BIAS of 3 cm for the run that implemented the depth-averaged velocity of the entire water column and water level (Run 6). At Faro Ocean Buoy the RMSE and BIAS were reduced by 5.4 cm and 5.1 cm, respectively for the run with surface currents and water level (Run 3). At Sines Buoy this same run reduced the RMSE by 2.9 cm and the run with the depth-averaged velocity for the entire water column (Run 5) reduced the BIAS by 8.6 cm.

For the storm with waves coming from E-SE, RMSE was reduced by approximately 6% at all buoy locations, corresponding to a reduction in RMSE by 2 cm. At Faro Ocean and Faro Coast buoys the run that achieved the lowest reduction in RMSE was the run with the depth-averaged velocity up to 10 m of the water column (Run 4). At Sines Buoy, the run that achieved the lowest reduction in RMSE was the run with surface currents and water level (Run 3). The BIAS had a more significant improvement where it was reduced by 17% or 3 cm at Faro Coast Buoy and Sines Buoy. At Faro Ocean Buoy it was reduced by 5% or 1.6 cm. The magnitude of improvement of just a few centimeters is a common conclusion reached in previous literature, where improvements of only a few centimeters have been observed [2,6,9].

When comparing results that incorporated current velocity from various depths of the water column, the current fields incorporating current velocity only up to a limited depth (surface, 10 m, and 20 m) produced the most accurate results. This signifies the higher influence that the upper layers of the water column have on wave propagation in this region. Including current velocity layers from the deepest parts of the water column, where the current velocity is weaker, reduces the overall magnitude of the depth-averaged 2D current field, making it less representative of the currents that actually influence wave propagation.

In this study, the simulation with waves coming from E-SE produced the most accurate results, even better than the CMEMS model, which was used to force the model with wave conditions at the boundary. The smallest differences were seen when the significant wave height was smaller, which is evident in the April 2023 simulation. This result is similar to the finding from Causio et al. [6] who found that coupling between hydrodynamic and wave models was most effective in periods with higher wave activity.

In terms of the spatial distribution of the effects of currents and water levels on wave propagation in the southwest Iberian Peninsula, the storm with waves from W-SW resulted in larger areas within the computational domain with increases in Hs. This is probably due to the westward circulation on the southern coast interacting with the westerly waves. Further offshore on the south coast, waves and currents flowed in the same direction, which resulted in decreases in Hs. In the storm with E-SE waves, decreases in Hs occurred over a broader area within the computational domain due to the water circulation and waves flowing in the same easterly direction. The highest changes in Hs for this storm were observed on the southeast coastline, due to the time varying water levels from SOMA. As the ocean circulation patterns and geomorphology of this region are both complex, accounting for all variables such as current velocity and water levels in wave propagation models is crucial. Previous studies that have developed wave models for this region have neglected these imperative forcing variables, thus overlooking crucial processes that impact the entire coastline.

Figure 12 quantifies the impact that both the magnitude of current velocity along with the wave-current angle difference had on changes in Hs for both the westerly and easterly storms. The difference between current and wave angles were computed and binned into 3-degree bins for the entire spatial domain. A wave-current angle difference of 30° or less implies a following current whereas a wave-current angle difference of 150° – 180° implies an opposing current. Current velocity magnitude was binned into 1 m/s bins. Figure 12 shows the average change in Hs during each simulation for the run with no hydrodynamic input and the run with surface current velocity and water levels.

Results of the westerly storm (Figure 12a) demonstrated larger decreases in Hs due to following currents (wave-current angle difference < 30°) and smaller increases in Hs in the case of opposing currents (wave-current angle difference between 150° and 180°). The velocity magnitude for the case of following currents was much higher compared to the velocity of opposing currents, which further justifies the larger changes in Hs. The easterly storm (Figure 12b) did not have the same effect as both the current velocity magnitude and resulting changes in Hs were smaller. The effects of the currents during this storm resulted in mostly decreases in Hs independent of the wave-current angle difference. This result is consistent with a study by Barnes et al. [9], who also found that when current velocity was less than 0.5 m/s, decreases in Hs were observed regardless of the wave-current angle difference. The increase in Hs as shown in the spatial map of Figure 9a reveal that these increases occur near the coastline, which is due to the inclusion of water levels from SOMA. Figure 9 and Figure 12 both reveal slight increases in Hs when waves and currents meet at right angles.

While the differences between simulations at all time periods were relatively small, the simulation that included surface currents and water levels consistently achieved results that were most similar to observations. The improved performance of the model when incorporating surface currents is consistent with previous studies who have used surface currents in their wave models [2,5,9]. This highlights the necessity of including the coupling when simulating waves in this study area. Future research will also run SWAN simulations over SOMA’s level 2 computational grid with a resolution of 1 km to evaluate the impacts of currents in shallower water at a higher resolution. This will further allow for a deeper investigation on the impacts of nearshore water circulation on the waves. This was not done for this study due to the high computational cost.

Limitations arise from the two-way coupling between wave and ocean circulation models that use depth-integrated velocities over an unstructured mesh. As these 2D models do not compute velocity at each vertical grid point, the computational time is much faster, but this two-way coupling ignores the influence of current velocity from the most significant layers attributing to wave propagation i.e., the surface layer. It is crucial to quantify current velocity over the vertical direction of the water column, and to properly compute the 2D current field attributed to wave propagation. The overall findings of this work demonstrate that simply computing the depth-averaged velocity over the entire water column yields less accurate results of wave propagation, and thus it is imperative to extract current velocity from a specific depth when two-way coupling with wave propagation models.

4. Conclusions

The wave model SWAN has been implemented, calibrated, and validated for the southern coast of the Iberian Peninsula. The standalone SWAN setup achieved results very similar to buoy observations, with low RMSE values ranging from 0.156 m to 0.407 m and low BIAS values ranging from -0.016 m to 0.255 m at several buoy stations during the validation period. It thus can be used to accurately represent wave characteristics in the south coast of Portugal. Improvements to this model have been investigated by assessing several model configurations that include various hydrodynamic forcing conditions from a previously developed model (SOMA). This was achieved by extracting water levels and a 2D current velocity field at different depths throughout the water column and testing the results under storm conditions coming from different directions. When the wave direction and current velocity followed the same direction, significant wave height decreased by approximately 10 cm. In cases of opposing currents, the significant wave height increased by approximately 20 cm. For the E-SE storm, Hs increased significantly along the southeast coast due to the temporally varying water levels provided by SOMA. The coupling of wave and hydrodynamic models can capture these physical processes that would otherwise be missed in standalone wave or hydrodynamic models. This is especially crucial for the south coast of Portugal, given the dynamic nature of the coastline, and the presence of rip currents in some areas which pose a hazard.

Results of this hydrodynamic forcing on the wave model showed small differences in the results that varied the 2D current velocity field, but overall improvements were made with the addition of current velocity and water level from SOMA. The configuration that produced the best overall results were SOMA surface currents and water level. This configuration will be used for all future operational runs of the wave model. The findings from this study present novel results on how water circulation impacts the wave field in the south coast of Portugal based on numerical models. Future work is planned to assess the influence of wave dynamics on water circulation by evaluating the effects of radiation stress and bottom orbital velocities. These mechanisms are crucial for gaining a deeper understanding of shallow water processes, forecasting events near the coast, and analyzing sediment transport.

Appendix A.1

Figure A1. Time series of Hs, Dir, and Tp for each model configuration during the calibration period (February 2024) at Faro Oceanic Buoy.

Figure A1. Time series of Hs, Dir, and Tp for each model configuration during the calibration period (February 2024) at Faro Oceanic Buoy.

Figure A2. Time series of Hs, Dir, and Tp for each model configuration during the calibration period (February 2024) at Cadiz Buoy.

Figure A2. Time series of Hs, Dir, and Tp for each model configuration during the calibration period (February 2024) at Cadiz Buoy.

Figure A3. Time series of Hs, Dir, and Tp for each model configuration during the calibration period (February 2024) at Sines Buoy.

Figure A3. Time series of Hs, Dir, and Tp for each model configuration during the calibration period (February 2024) at Sines Buoy.

Table A1. Error metrics of mean wave direction comparing model results of each calibration setting with buoy observations. Values in bold indicate the lowest (RMSE or BIAS) or highest (MSS or R) among all cases.

	Faro Ocean Buoy				Cadiz Buoy				Sines Buoy
	RMSE (deg)	BIAS (deg)	MSS	R	RMSE (deg)	BIAS (deg)	MSS	R	RMSE (m)	BIAS (m)	MSS	R
Case1	38.067	1.631	0.893	0.801	20.933	0.729	0.957	0.932	6.857	-1.781	0.985	0.976
Case2	41.279	3.578	0.865	0.760	21.722	0.318	0.953	0.927	6.058	-1.641	0.988	0.979
Case3	31.859	-0.904	0.930	0.866	20.648	1.005	0.959	0.932	7.607	-1.882	0.982	0.972
Case4	37.551	1.481	0.895	0.806	21.276	0.638	0.956	0.929	6.751	-1.755	0.985	0.976
Case5	30.589	-2.218	0.936	0.878	20.394	1.202	0.961	0.934	8.232	-1.953	0.979	0.968
Case6	37.340	1.282	0.897	0.809	20.718	0.684	0.958	0.934	6.700	-1.823	0.986	0.977
Case7	41.633	3.735	0.861	0.755	21.500	0.249	0.954	0.930	5.974	-1.684	0.988	0.979

Table A1. Error metrics of mean wave direction comparing model results of each calibration setting with buoy observations. Values in bold indicate the lowest (RMSE or BIAS) or highest (MSS or R) among all cases.

	Faro Ocean Buoy				Cadiz Buoy				Sines Buoy
	RMSE (deg)	BIAS (deg)	MSS	R	RMSE (deg)	BIAS (deg)	MSS	R	RMSE (m)	BIAS (m)	MSS	R
Case1	38.067	1.631	0.893	0.801	20.933	0.729	0.957	0.932	6.857	-1.781	0.985	0.976
Case2	41.279	3.578	0.865	0.760	21.722	0.318	0.953	0.927	6.058	-1.641	0.988	0.979
Case3	31.859	-0.904	0.930	0.866	20.648	1.005	0.959	0.932	7.607	-1.882	0.982	0.972
Case4	37.551	1.481	0.895	0.806	21.276	0.638	0.956	0.929	6.751	-1.755	0.985	0.976
Case5	30.589	-2.218	0.936	0.878	20.394	1.202	0.961	0.934	8.232	-1.953	0.979	0.968
Case6	37.340	1.282	0.897	0.809	20.718	0.684	0.958	0.934	6.700	-1.823	0.986	0.977
Case7	41.633	3.735	0.861	0.755	21.500	0.249	0.954	0.930	5.974	-1.684	0.988	0.979

Table A2. Error metrics of peak period comparing model results of each calibration setting with buoy observations. Values in bold indicate the lowest (RMSE or BIAS) or highest (MSS or R) among all cases.

	Faro Buoy				Cadiz Buoy				Sines Buoy
	RMSE (s)	BIAS (s)	MSS	R	RMSE (s)	BIAS (s)	MSS	R	RMSE (s)	BIAS (s)	MSS	R
Case1	2.812	1.020	0.780	0.635	3.955	1.976	0.702	0.540	1.948	0.142	0.869	0.754
Case2	2.839	1.034	0.776	0.628	3.958	1.998	0.700	0.540	1.947	0.144	0.869	0.754
Case3	2.812	1.026	0.779	0.634	3.900	1.917	0.711	0.553	1.948	0.142	0.868	0.753
Case4	2.812	1.021	0.780	0.635	3.955	1.978	0.701	0.540	1.948	0.142	0.869	0.754
Case5	2.769	0.988	0.787	0.646	3.888	1.907	0.714	0.556	1.948	0.145	0.868	0.753
Case6	2.806	1.033	0.781	0.637	3.985	2.025	0.697	0.535	1.945	0.150	0.869	0.755
Case7	2.833	1.047	0.777	0.631	4.010	2.053	0.692	0.529	1.945	0.151	0.869	0.754

Table A2. Error metrics of peak period comparing model results of each calibration setting with buoy observations. Values in bold indicate the lowest (RMSE or BIAS) or highest (MSS or R) among all cases.

	Faro Buoy				Cadiz Buoy				Sines Buoy
	RMSE (s)	BIAS (s)	MSS	R	RMSE (s)	BIAS (s)	MSS	R	RMSE (s)	BIAS (s)	MSS	R
Case1	2.812	1.020	0.780	0.635	3.955	1.976	0.702	0.540	1.948	0.142	0.869	0.754
Case2	2.839	1.034	0.776	0.628	3.958	1.998	0.700	0.540	1.947	0.144	0.869	0.754
Case3	2.812	1.026	0.779	0.634	3.900	1.917	0.711	0.553	1.948	0.142	0.868	0.753
Case4	2.812	1.021	0.780	0.635	3.955	1.978	0.701	0.540	1.948	0.142	0.869	0.754
Case5	2.769	0.988	0.787	0.646	3.888	1.907	0.714	0.556	1.948	0.145	0.868	0.753
Case6	2.806	1.033	0.781	0.637	3.985	2.025	0.697	0.535	1.945	0.150	0.869	0.755
Case7	2.833	1.047	0.777	0.631	4.010	2.053	0.692	0.529	1.945	0.151	0.869	0.754

Table A3. Sensitivity tests for directional spreading.

Sensitivity Test	Directional Spreading	Simulation Timeframe
1	Degree E30 N30 S30 W40	February 2024
2	Degree E40 N40 S40 W40	February 2024
3	Power E2 N2 S3 W4	February 2024
4	Power E1.5 N1.5 S1.5 W1.5	February 2024
5	Power E2 N2 S3 W3	February 2024
6	Power E3 N3 S3 W3	February 2024
7	Power E15 N15 S15 W15	February 2024
8	Power E2 N4 S3 W4	February 2024
9	Power E4 N4 S4 W4	February 2024
10	Power E2 N2 S2 W4	February 2024
11	Power E2 N2 S2 W3	February 2024
12	Power E2 N2 S3 W4	November 2023
13	Power E2 N2 S2 W4	November 2023
14	Power E2 N2 S2 W3	November 2023
15	Power E2 N2 S2 W4	21-10-2023 – 07-11-2023 (westerly storm)
16	Power E2 N2 S2 W4	08-02-2023 – 22-02-2023 (easterly storm)

Table A3. Sensitivity tests for directional spreading.

Sensitivity Test	Directional Spreading	Simulation Timeframe
1	Degree E30 N30 S30 W40	February 2024
2	Degree E40 N40 S40 W40	February 2024
3	Power E2 N2 S3 W4	February 2024
4	Power E1.5 N1.5 S1.5 W1.5	February 2024
5	Power E2 N2 S3 W3	February 2024
6	Power E3 N3 S3 W3	February 2024
7	Power E15 N15 S15 W15	February 2024
8	Power E2 N4 S3 W4	February 2024
9	Power E4 N4 S4 W4	February 2024
10	Power E2 N2 S2 W4	February 2024
11	Power E2 N2 S2 W3	February 2024
12	Power E2 N2 S3 W4	November 2023
13	Power E2 N2 S2 W4	November 2023
14	Power E2 N2 S2 W3	November 2023
15	Power E2 N2 S2 W4	21-10-2023 – 07-11-2023 (westerly storm)
16	Power E2 N2 S2 W4	08-02-2023 – 22-02-2023 (easterly storm)

Table A4. Error metrics of significant wave height during the calibration period (February 2024). Directional spreading refers to values for directional spreading for each segment on each side (E = east, N = north, S = south, W = west). Simulation cases with the best metrics results are highlighted in bold.

	Faro Ocean Buoy				Cadiz Buoy				Sines Buoy
Directional Spreading	RMSE (m)	BIAS (m)	MSS	R	RMSE (m)	BIAS (m)	MSS	R	RMSE (m)	BIAS (m)	MSS	R
Degree E30 N30 S30 W40	0.404	-0.070	0.955	0.916	0.370	0.082	0.958	0.923	0.421	-0.022	0.978	0.956
Degree E40 N40 S40 W40	0.415	-0.027	0.952	0.908	0.424	0.154	0.945	0.907	0.398	-0.025	0.980	0.960
Power E2 N2 S3 W4	0.348	-0.107	0.966	0.941	0.280	-0.032	0.976	0.954	0.387	0.007	0.981	0.962
Power E1.5 N1.5 S1.5 W1.5	0.390	-0.050	0.958	0.920	0.374	0.103	0.957	0.923	0.401	-0.014	0.980	0.960
Power E2 N2 S3 W3	0.359	-0.104	0.964	0.936	0.292	-0.015	0.974	0.950	0.394	0.001	0.980	0.961
Power E3 N3 S3 W3	0.359	-0.102	0.964	0.936	0.291	-0.013	0.974	0.950	0.410	-0.001	0.979	0.959
Power E15 N15 S15 W15	0.515	-0.346	0.926	0.920	0.434	-0.302	0.940	0.939	0.447	-0.004	0.975	0.951
Power E2 N4 S3 W4	0.349	-0.104	0.966	0.940	0.279	-0.031	0.976	0.954	0.414	0.003	0.979	0.958
Power E4 N4 S4 W4	0.355	-0.125	0.965	0.940	0.276	-0.066	0.977	0.957	0.414	0.001	0.979	0.958
Power E2 N2 S2 W4	0.346	-0.072	0.967	0.939	0.297	0.026	0.973	0.948	0.388	0.009	0.981	0.962
Power E2 N2 S2 W3	0.354	-0.068	0.965	0.935	0.314	0.041	0.970	0.943	0.394	0.003	0.980	0.961

Table A4. Error metrics of significant wave height during the calibration period (February 2024). Directional spreading refers to values for directional spreading for each segment on each side (E = east, N = north, S = south, W = west). Simulation cases with the best metrics results are highlighted in bold.

	Faro Ocean Buoy				Cadiz Buoy				Sines Buoy
Directional Spreading	RMSE (m)	BIAS (m)	MSS	R	RMSE (m)	BIAS (m)	MSS	R	RMSE (m)	BIAS (m)	MSS	R
Degree E30 N30 S30 W40	0.404	-0.070	0.955	0.916	0.370	0.082	0.958	0.923	0.421	-0.022	0.978	0.956
Degree E40 N40 S40 W40	0.415	-0.027	0.952	0.908	0.424	0.154	0.945	0.907	0.398	-0.025	0.980	0.960
Power E2 N2 S3 W4	0.348	-0.107	0.966	0.941	0.280	-0.032	0.976	0.954	0.387	0.007	0.981	0.962
Power E1.5 N1.5 S1.5 W1.5	0.390	-0.050	0.958	0.920	0.374	0.103	0.957	0.923	0.401	-0.014	0.980	0.960
Power E2 N2 S3 W3	0.359	-0.104	0.964	0.936	0.292	-0.015	0.974	0.950	0.394	0.001	0.980	0.961
Power E3 N3 S3 W3	0.359	-0.102	0.964	0.936	0.291	-0.013	0.974	0.950	0.410	-0.001	0.979	0.959
Power E15 N15 S15 W15	0.515	-0.346	0.926	0.920	0.434	-0.302	0.940	0.939	0.447	-0.004	0.975	0.951
Power E2 N4 S3 W4	0.349	-0.104	0.966	0.940	0.279	-0.031	0.976	0.954	0.414	0.003	0.979	0.958
Power E4 N4 S4 W4	0.355	-0.125	0.965	0.940	0.276	-0.066	0.977	0.957	0.414	0.001	0.979	0.958
Power E2 N2 S2 W4	0.346	-0.072	0.967	0.939	0.297	0.026	0.973	0.948	0.388	0.009	0.981	0.962
Power E2 N2 S2 W3	0.354	-0.068	0.965	0.935	0.314	0.041	0.970	0.943	0.394	0.003	0.980	0.961

Table A5. Error metrics of significant wave height during the validation period (November 2023). Directional spreading refers to values for directional spreading for each segment on each side (E = east, N = north, S = south, W = west). The simulation configuration with the best metrics results is highlighted in bold.

	Faro Ocean Buoy				Faro Coast Buoy				Cadiz Buoy				Sines Buoy
Directional Spreading	RMSE	BIAS	MSS	R	RMSE	BIAS	MSS	R	RMSE	BIAS	MSS	R	RMSE	BIAS	MSS	R
E2 N2 S3 W4	0.341	0.058	0.956	0.93	0.217	-0.145	0.950	0.947	0.251	0.054	0.957	0.942	0.335	-0.024	0.989	0.982
E2 N2 S2 W4	0.357	0.092	0.953	0.93	0.204	-0.121	0.956	0.948	0.300	0.118	0.942	0.933	0.322	-0.013	0.990	0.983
E2 N2 S2 W3	0.368	0.101	0.950	0.922	0.183	-0.093	0.964	0.951	0.328	0.149	0.931	0.925	0.328	-0.022	0.990	0.983

Table A5. Error metrics of significant wave height during the validation period (November 2023). Directional spreading refers to values for directional spreading for each segment on each side (E = east, N = north, S = south, W = west). The simulation configuration with the best metrics results is highlighted in bold.

	Faro Ocean Buoy				Faro Coast Buoy				Cadiz Buoy				Sines Buoy
Directional Spreading	RMSE	BIAS	MSS	R	RMSE	BIAS	MSS	R	RMSE	BIAS	MSS	R	RMSE	BIAS	MSS	R
E2 N2 S3 W4	0.341	0.058	0.956	0.93	0.217	-0.145	0.950	0.947	0.251	0.054	0.957	0.942	0.335	-0.024	0.989	0.982
E2 N2 S2 W4	0.357	0.092	0.953	0.93	0.204	-0.121	0.956	0.948	0.300	0.118	0.942	0.933	0.322	-0.013	0.990	0.983
E2 N2 S2 W3	0.368	0.101	0.950	0.922	0.183	-0.093	0.964	0.951	0.328	0.149	0.931	0.925	0.328	-0.022	0.990	0.983

Table A6. Error metrics of significant wave height during the storm with westerly waves with the directional spreading set to power E2 N2 S2 W4.

Level 1 SIMULATION 2023-10-21 12:00:00 - 2023-11-07 12:00:00 - WAVES COMING FROM WEST
	Faro Oceanic Buoy				Faro Coast Buoy				Cadiz Buoy				Sines Buoy
RUN ID	RMSE (m)	BIAS (m)	MSS	R	RMSE (m)	BIAS (m)	MSS	R	RMSE (m)	BIAS (m)	MSS	R	RMSE (m)	BIAS (m)	MSS	R
1	0.539	0.426	0.898	0.922	0.283	-0.127	0.945	0.933	0.517	0.417	0.890	0.927	0.406	-0.046	0.981	0.967
2	0.540	0.426	0.898	0.922	0.284	-0.127	0.945	0.932	0.516	0.417	0.890	0.927	0.406	-0.046	0.981	0.967
3	0.459	0.339	0.921	0.927	0.255	0.011	0.952	0.919	0.375	0.250	0.937	0.935	0.390	0.019	0.983	0.968
4	0.461	0.341	0.921	0.928	0.256	0.018	0.952	0.920	0.381	0.258	0.935	0.935	0.394	0.020	0.983	0.967
5	0.462	0.343	0.921	0.928	0.257	0.020	0.952	0.921	0.385	0.262	0.935	0.936	0.391	-0.165	0.983	0.975

Table A6. Error metrics of significant wave height during the storm with westerly waves with the directional spreading set to power E2 N2 S2 W4.

Level 1 SIMULATION 2023-10-21 12:00:00 - 2023-11-07 12:00:00 - WAVES COMING FROM WEST
	Faro Oceanic Buoy				Faro Coast Buoy				Cadiz Buoy				Sines Buoy
RUN ID	RMSE (m)	BIAS (m)	MSS	R	RMSE (m)	BIAS (m)	MSS	R	RMSE (m)	BIAS (m)	MSS	R	RMSE (m)	BIAS (m)	MSS	R
1	0.539	0.426	0.898	0.922	0.283	-0.127	0.945	0.933	0.517	0.417	0.890	0.927	0.406	-0.046	0.981	0.967
2	0.540	0.426	0.898	0.922	0.284	-0.127	0.945	0.932	0.516	0.417	0.890	0.927	0.406	-0.046	0.981	0.967
3	0.459	0.339	0.921	0.927	0.255	0.011	0.952	0.919	0.375	0.250	0.937	0.935	0.390	0.019	0.983	0.968
4	0.461	0.341	0.921	0.928	0.256	0.018	0.952	0.920	0.381	0.258	0.935	0.935	0.394	0.020	0.983	0.967
5	0.462	0.343	0.921	0.928	0.257	0.020	0.952	0.921	0.385	0.262	0.935	0.936	0.391	-0.165	0.983	0.975

Table A7. Error metrics of significant wave height during the storm with easterly waves with the directional spreading set to power E2 N2 S2 W4.

	Faro Oceanic Buoy				Faro Coast Buoy				Sines Buoy
RUN ID	RMSE (m)	BIAS (m)	MSS	R	RMSE (m)	BIAS (m)	MSS	R	RMSE (m)	BIAS (m)	MSS	R
1	0.385	-0.265	0.964	0.962	0.367	-0.266	0.954	0.958	0.356	0.162	0.901	0.878
2	0.385	-0.265	0.964	0.962	0.368	-0.266	0.954	0.958	0.356	0.162	0.901	0.878
3	0.373	-0.275	0.966	0.969	0.358	-0.254	0.955	0.956	0.333	0.118	0.912	0.880
5	0.360	-0.255	0.969	0.969	0.352	-0.241	0.957	0.955	0.341	0.135	0.908	0.880

Table A7. Error metrics of significant wave height during the storm with easterly waves with the directional spreading set to power E2 N2 S2 W4.

	Faro Oceanic Buoy				Faro Coast Buoy				Sines Buoy
RUN ID	RMSE (m)	BIAS (m)	MSS	R	RMSE (m)	BIAS (m)	MSS	R	RMSE (m)	BIAS (m)	MSS	R
1	0.385	-0.265	0.964	0.962	0.367	-0.266	0.954	0.958	0.356	0.162	0.901	0.878
2	0.385	-0.265	0.964	0.962	0.368	-0.266	0.954	0.958	0.356	0.162	0.901	0.878
3	0.373	-0.275	0.966	0.969	0.358	-0.254	0.955	0.956	0.333	0.118	0.912	0.880
5	0.360	-0.255	0.969	0.969	0.352	-0.241	0.957	0.955	0.341	0.135	0.908	0.880

Figure A4. Time series of water levels at (a) Faro Ocean Buoy, (b) Faro Coast Buoy, (c) Cadiz Buoy, and (d) Sines Buoy during the storm with westerly waves (2023-10-21 – 2023-11-07).

Figure A4. Time series of water levels at (a) Faro Ocean Buoy, (b) Faro Coast Buoy, (c) Cadiz Buoy, and (d) Sines Buoy during the storm with westerly waves (2023-10-21 – 2023-11-07).

Figure A5. Time series of water levels at (a) Faro Ocean Buoy, (b) Faro Coast Buoy, (c) Cadiz Buoy, and (d) Sines Buoy during the storm with easterly waves (2023-02-08 – 2023-02-22).

Figure A5. Time series of water levels at (a) Faro Ocean Buoy, (b) Faro Coast Buoy, (c) Cadiz Buoy, and (d) Sines Buoy during the storm with easterly waves (2023-02-08 – 2023-02-22).