Hindcast of Signifi cant Wave Heights in Sheltered Basins Using Machine Learning and the Copernicus Database Kratkoročna prognoza značajnih valnih visina u zaštićenim akvatorijima koristeći se strojnim učenjem i Copernicus bazom podataka

Long-term time series of wave parameters play a critical role in coastal structure design and maritime activities. At sites with limited buoy measurements, methods are used to extend the available time series data. To date, wave hindcasting research using machine learning methods has mainly focused on fi lling in missing buoy measurements or fi nding a mapping function between two nearshore buoy locations. This work aims to implement machine learning methods for hindcasting wave parameters using only publicly available Copernicus data. Ensemble regression and artifi cial neural networks were used as machine learning methods and the optimal hyperparameters were determined by the Bayesian optimization algorithm. As inputs, data from the MEDSEA reanalysis wave model were used for the wave parameters and data from the ERA5 atmospheric reanalysis model were used for the wind parameters. The results of this study show that the normalized RMSE of the


INTRODUCTION / Uvod*
Knowledge of a long-term time series of wave climate (e.g.signifi cant wave height, peak wave period, etc.) at a location is essential for planning, operation, and maintenance of maritime activities [1], fl ood protection engineering design [2], and coastal vulnerability assessment [3].A long wave height time series is widely recognized as key to reliable long-term signifi cant wave height forecasting or hindcasting, (e.g., to defi ne the return level period of signifi cant wave height) [4].This long-term wave height time series can be * Corresponding author constructed using wave hindcast methods when data is missing or not suffi cient [5].The high return period of signifi cant wave heights (50-year, 100-year return period, etc.) is required as input for the planning phase of coastal structures, long-term morphodynamics process studies [6], or as a boundary value for numerical nearshore wave climate models [7].The capability and availability of a reliable long-term wave database are of paramount importance to the ocean and coastal engineers.This is especially true for coastal seas, for which there are much less data compared to deeper seas.
Wave buoys provide long-term time series of wave parameters at a variety of locations if national or international climate monitoring networks maintain them continuously.If this is the case, wave records could be available for 15-30 years [8].In the Mediterranean Sea, for example, the Hellenic Centre for Marine Research, and Institute of Oceanography started monitoring the wave fi eld with their buoy ATHOS in May of 2000, and the Spanish Harbor Authority with their buoy 6100197 in April of 1993 [9].In addition, there are 39 years of continuous wave measurements at the Aqua Alta oceanographic tower in the Adriatic Sea [10].More commonly, however, wave buoys are only maintained for only a few years during specifi c campaigns before they are usually recovered.This does not provide long enough time series for long-term forecasting, so methods to extend the time series should be applied [11].
One way to increase the available data set at sites where measurements have been made for only a limited duration is to use historical wind data as input to wave hindcasting to reconstruct a long-term wave time series.There is a wide range of solutions for this, ranging from simple empirical models [11; 12] to complex wave-generating numerical models [13; 14].Numerical models are limited by available computational power, detailed bathymetry data at locations sheltered by islands, their complexity, and diffi cult-to-determine coeffi cient (e.g.white-capping parameters, bed frictional dissipation, depthlimited wave breaking, etc.) [15], while the main advantage compared to the simpler empirical models is that they perform physically based calculations.Hindcasting of wave heights at a local level from global/regional reanalysis models has also been performed using locally based wave numerical models (SWAN) using the reanalysis data as input [7; 16].
Climate reanalysis products can provide data for historical wind and wave data for a specifi c location or an entire region.The website Advancing Reanalysis [17] provides a visual comparison of the various reanalysis products [18].Reanalysis is a scientifi c method of creating a complete record of changes in weather and climate over time.It usually spans decades or more and covers the entire Earth or focuses on a specifi c region.The information obtained from reanalysis is widely used for monitoring, comparison, determining the causes of climate variability, and supplementing climate predictions.The Copernicus database provides several reanalysis models that include wave data, such as ERA5 global wave climate [19], WAVERYS -Global Ocean Waves Reanalysis [20], and MEDSEA -Mediterranean Sea Waves Reanalysis [21].WAVERYS and ERA5 both have wave height grided at low resolution (0.2° and 0.25°, respectively) when compared to MEDSEA (1/24°).Furthermore, the WAVERYS model reanalysis has a temporal resolution of 3 h, compared to the 1 h temporal resolution of MEDSEA and ERA5.The disadvantage of the MEDSEA regional model compared to the other two is its relatively short time range 1993-2020 (last accessed on 27.04.2022.), while the ERA5 model has data back to 1979.A longer wave and wind history could contribute to a longer wave time series when hindcasting.In addition to the wave data provided, ERA5 also off ers wind reanalysis data at a spatial resolution of 0.25°.Although the MEDSEA reanalysis model is the most detailed Copernicus numerical model reanalysis in the Mediterranean Sea (both in spatial and temporal terms), Korres et al. [21] still observed low accuracy when validating the reanalysis data to buoy measurements in well-sheltered areas.They note that the worst model performance is observed for the Adriatic Sea.Korres et al. [21] argue that the Adriatic Sea is shallow, enclosed, and bounded by complex topography, and therefore not adequately represented by the spatial resolution of the forcing wind and possibly by the spatial resolution of the wave model.
There is an alternative to performing hindcasting using complex numerical reanalysis models, such as hindcasting using machine learning models, from simple models such as stepwise linear models to complex artifi cial neural networks (ANN), and more.Machine learning models are capable of mapping complex non-linear functions between inputs and outputs when suffi cient training data is available [22].These methods are used in many applications in the ocean and coastal engineering in many applications, such as wave forecasting with several hours of lead time [1; 15; 23; 24], wave runup [25], beach sediment transport [26; 27], beach nourishment requirements [28], etc.As for wave hindcasting, it has been mainly performed using ANN based on wind reanalysis data at a regional level to downscale it to a local level [8; 29].Machine learning models have been used when data were missing in the measured wave time series to fi ll in missing wave heights [30] or to fi nd a mapping function between wave data at a nearshore site, by using data from one or more nearby off shore sites [31].To train and test a machine learning model of this kind, data from a short-term wave buoy campaign can be used.The insights gained from machine learning techniques can additionally be used to improve the hindcast predictions given by a reanalysis wave model to better fi t the available campaign wave buoy measurements.
This study aims to present a modeling chain that uses machine learning models with state-of-the-art regional reanalysis wave data (MEDSEA) and global reanalysis wind data (ERA5) as input for long-term reconstruction of signifi cant wave heights at three diff erent locations in the Adriatic Sea, which has proven to be the most challenging region for MEDSEA.We propose a method to rapidly improve MEDSEA wave prediction at sheltered locations using machine learning tools.Consequently, a validated machine learning model can extend the wave time series beyond the duration of the buoy measurement at the site, which is initially available to the user.The paper evaluates two diff erent machine learning techniques: artifi cial neural networks (ANN) and ensemble learning with regression trees.Validity is demonstrated by comparing the newly created hindcast using a machine learning model with the MEDSEA hindcast time series as a benchmark.This research explores the possibility of using a rapid methodology to extend data from wave monitoring campaigns to a time that has not been measured in the past.A reliably extended long-term wave time series is needed for predicting the return period from 1 to 100 years and consequently for planning coastal structures.In addition, the wave time series can be extended by periodically updating the input data set as new MEDSEA and ERA5 data become available over time.
This article is organized as follows: the methodology is outlined in section 2, and the paper will then go on to compare the MEDSEA modeled data with measured data in section 3.1 and analyze the infl uence of the machine learning correction method in section 3.2.Section 4 will present the paper's conclusions.

Wave data / Podaci o valu
The Copernicus Marine Environment Monitoring Service (CMEMS) provides a 27-year wave reanalysis product, MEDSEA, covering the period from January 1993 to December 2019 for the Mediterranean Sea [21].This wave reanalysis is based on the advanced third-generation wave model WAM Cycle 4.6.2[13; 32].It explicitly solves the wave transport equations without assuming the wave spectrum shape.Included source terms are wind input, white capping dissipation, nonlinear transfer, and bottom friction.The wind and white-capping dissipation terms are based on Janssen's quasilinear theory of wind-wave generation [33; 34], while the bottom friction term is based on the empirical JONSWAP model [35].The numerical model discretizes the wave spectra using 32 frequencies covering a logarithmically scaled frequency band from 0.04177 Hz to 0.8018 Hz (with wave periods from about 1 s to 24 s) and 24 uniformly distributed directional bins (bin size of 15 deg).Winds from the ERA5 reanalysis 10 m above the sea surface (Copernicus Climate Service -ECMWF) are forcing the numerical wave model.The bathymetric map was created using the GEBCO bathymetric dataset [36].In addition, the reanalysis includes an assimilation scheme that uses the signifi cant wave heights determined from altimeters and adjusts the wave spectrum at each grid point accordingly (originally developed by [37]).
Korres et al. [21] show typical diff erences between the MEDSEA reanalysis model and the in-situ and satellite observations (RMSE) of 0.23 ± 0.012 m and 0.24 ± 0.01 m respectively, and BIAS of -0.06 ± 0.022 m (7% ± 3% relative to the observed mean) and -0.05 ± 0.011 m (4% ± 1%) for the Mediterranean Sea as a whole.BIAS is predominately negative, indicating widespread underestimation of the measured wave heights by the reanalysis.In the Adriatic Sea, the model accuracy for buoy 61217 deteriorates further (RMSE) to 0.27 m and BIAS to -0.14, as stated by the model authors (Table 1).
The MEDSEA reanalysis model provides 2D hourly instantaneous fi elds (Table 2).It generates data every hour with a horizontal resolution of 1/24º.This dataset will be used as inputs to the machine learning model to hindcast wave heights that are more accurate to measured wave data at the Croatian coast (described in section 2.2.).

Wind data / Podaci o vjetru
ERA5 is the fi fth generation of ECMWF's global climate and weather reanalysis for recent decades (starting in 1979).The reanalysis integrates model data with observations from around the world to produce a complete and consistent global dataset.This reanalysis also uses data assimilation, which optimally combines the forecast with newly available observations every 12 hours to produce an optimal estimate of atmospheric conditions.Because the reanalysis is not forced to produce timely forecasts, there is more time to collect new data and incorporate historical data to improve the quality of the reanalysis product.ERA5 provides hourly estimates for a variety of atmospheric, ocean wave, and land surface parameters.The wind data used in this work were resampled to a regular lat-lon grid with a resolution of 0.25º.Recent studies have compared the ERA5 reanalysis wind data with measured wind station data in the Adriatic and Mediterranean Seas.They showed moderate accuracy of ERA5 when compared to more detailed regional wind reanalysis, and signifi cant underprediction (2 m/s on average) for high wind velocities (larger than 10 m/s) [38; 39].Nevertheless, the ERA5 is a state-of-the-art atmosphere reanalysis model, and is publicly available, which is important from a data availability standpoint.
Lag components of signifi cant wave height from CMEMS and wind magnitude from ECMWF were added to the predictor set to account for the historical components of these variables.Thus, at a given time in the reconstruction of the measured signifi cant wave height, the model has insight into the current wave height hindcast from MEDSEA and 10 previous predictions of wave height (MEDSEA) and wind magnitude (ERA5).Similarly, some researchers have added this to help predict wave height [8; 23].In this way, wind duration is taken into account, which is expected to aff ect the wave height reconstruction accuracy.Table 1 Stations in the western Adriatic Sea (Italy) for which model validation and performance metrics were conducted in [21]; location of the wave buoys is shown in Figure 1

Field wave measurements / Terenska mjerenja valova
The measurements were made using the well-known DATAWELL Waverider DWR MKIII, which was anchored in cooperation with the Hydrographic Institute of the Republic of Croatia.The moored wave rider measures wave direction, wave height, and peak period.The measured data is stored on the internal data logger of the buoy, but also through the HF antenna connection on the buoy, the data is transmitted to the RX -C receiver on shore.The receiver is connected to a computer with a software package needed to collect and analyze the data.It is also equipped with GPS for positioning and tracking the buoy.The high-capacity batteries inside the Waverider ensure operation for up to one year without battery replacement.

Ensemble regression (ER) -Least Squares Boosting Ensemble (LSBoost) / Ansambl regresija (ER) -Boosting ansambl najmanjih kvadrata (LSBoost)
Ensembles regression aggregates a set of trained weak learners (also called individual learners) to predict an ensemble response.The weak learner in this work is a decision tree, and the Least-squares boosting (LSBoost) method was used during cross-validation training of the ensemble regression model [40][41][42].The algorithm fi rst creates an initial model of the selected weak learner.Then at each time step, LSBoost fi ts a new weak learner to the current residuals, i.e., the diff erence between the observed response and the aggregated prediction of all learners created previously.Eventually, the aggregation of the weak learners should cause the residuals to converge.The LSBoost weighting uses the least-squares function as the loss function.In addition, Bayesian optimization was used to fi nd the best possible hyperparameters values for the ensemble training: Number of learning cycles, learning rate, minimum leaf size and the maximum number of splits.The entire process was iterated 50 times to fi nd the best hyperparameter combination.The optimization fi tness objective was the mean squared error (described in Section 2.4).After training, the model was regularized using the Lasso algorithm, which fi nds an optimal set of weak learner weights and reduces overfi tting to the data: where λ is the lasso parameter, h t is a weak learner in the ensemble trained on N observations with predictors x n , responses y n , and weights w n .The λ value for regularization used in this work was 0.0008.

Artifi cial neural network / Umjetna neuralna mreža
A multi-layer feed-forward network was used to fi t the function linking the input data (MEDSEA and ERA5, as described in section 2.1.)and the response data (measured fi eld data, as described in section 2.2.) [22].This is a common type of ANNs, where the output of each node in each layer is passed only to the following layer [43].An ANN consists of input, hidden, and output nodes arranged in layers.Each input node is connected to multiple nodes, which together form the hidden layer or multiple connected hidden layers.Information is passed from the input nodes forward to nodes in the hidden layer: where x i is the input variables, h j is the responses of the hidden layer neuron, w i is the weights, aj is the biases, and f is the activation function.Finally, the hidden layer was fully connected to the output layer, which consists of 1 node corresponding to the corrected signifi cant wave height.The algorithm used the limited-memory Broyden-Fletcher-Goldfarb-Shanno quasi-Newton algorithm (LBFGS), where the mean squared error (MSE) was the optimization objective for training the weights and biases [44].
In addition, Bayesian optimization was used to fi nd the best possible hyperparameter values for ANN training: activation function (relu, tanh, or sigmoid), lambda, number of layers, and layer sizes).The whole process was iterated 50 times to fi nd the best hyperparameter combination.Lambda is a regularization strength term that counteracts the tendency of the training procedure to overfi t the network to the training data.This is evident when statistical error metrics for the training and test datasets are drastically diff erent [45].
The input and response data were preprocessed to effi ciently train the ANNs.This normalization helped to avoid very small gradients and thus long training times.A common approach was used to normalize inputs and responses to fall within the range [-1, 1].

Statistical error metrics / Statistička metrika pogrešaka
To evaluate the prediction accuracy of the machine learning model and MEDSEA model, various statistical error metrics are used such as the coeffi cient of determination (R 2 ), normalized bias (NBIAS), the normalized root mean square error (NRMSE), mean absolute percentage error (MAPE), and scatter index (SI) as defi ned in Equations ( 5)-( 9) respectively: where P i is the i th prediction, O i is the i th observation, P is the average prediction, and O is the average observation.

RESULTS AND DISCUSSION / Rezultati i rasprava
To determine if the MEDSEA modeled signifi cant wave height data need to be corrected for the Rijeka, Istra, or Split sites, comparisons of modeled using MEDSEA and measured wave data are presented in Section 3.1.Based on the NRMSE and R 2 error metrics, the comparisons are presented and discussed whether a correction is needed.No resampling of the modeled MEDSEA was done with wave data since both have the same time step of 1 h.Then, machine learning correction methods (described in Sections 2.3.1.and 2.3.2.) are applied to the modeled wave data to achieve a better fi t to the measured data.The corrected signifi cant wave height data are compared to measured wave data for Split and Rijeka in Sections 3.2.1.for ensemble regression and 3.2.2. for ANN model fi t.Section 3.3 presents the improvements by adding lag components to the predictor set as described in Section 2.1.Finally, in Section 3.4.we evaluate the importance of the predictors.This is an inherent property of ensemble regression that can be easily extracted from a trained model.

Comparison of MEDSEA reanalysis modeled data and fi eld wave measurements / Usporedba MEDSEA reanalizom modeliranih podataka i stvarnih izmjerenih valova
MEDSEA-mod eled wave data showed good agreement with measurements near the Istra buoy, with an NRMSE value of 0.32 (Figure 2).This error is lower than the statistical metrics reported in the MEDSEA report (Table 1) [21] for other locations in the Adriatic Sea.Also, the modeled data showed low overprediction bias (2%).A similar NRMSE with the MEDSEA report values is to be expected since both the Istra buoy and the MEDSEA buoys (Figure 1) are not sheltered by islands that would reduce the fetch length below 30 km.
On the other hand, a substantial diff erence was found between the MEDSEA reanalysis modeled signifi cant wave heights in comparison to the measured wave data for the locations of Split and Rijeka (Figure 2).The NRMSE increased to 0.52 and 0.74 for Split and Rijeka, respectively.This is 62% and 131% more than the NRMSE for Istra.These error metrics are also higher than the NRMSE values reported for other locations in the Adriatic (Table 1) [21].The Split buoy wave measurements are moderately underpredicted (12%), while those at the Rijeka buoy are signifi cantly underpredicted (35%).It is important to point out that the most likely reason for this discrepancy is the extremely complex orography surrounding the Rijeka and Split buoys, whose fetch length is less than 30 km.This is not the case for the Istra and MEDSEA buoys.This decrease in MEDSEA model accuracy in the case of Split and Rijeka is likely caused by unresolved topography by the wind and wave models and fetch accuracy limitations caused by the wave model resolution, as reported by Korres et al. [21].This is all due to the eastern Adriatic Sea being enclosed basins near the coast with small fetch lengths dominated by wind waves.
Furthermore, the measured wave direction is spread mainly between the SW and SE directions, which is to be expected given the proximity of the buoys to the mainland in the NW direction (less than 1 km) (Figure 3 -right).Instead, the wave rose for the MEDSEA modeled data shows a strong presence of waves from the NE direction (Figure 3 -left).Although Korres et al. [21] did not compare the results of their modeled direction results to the measured wave directions, this preliminary observation shows that modeled wave directions for nearshore locations should be used with caution.
Overall, these results suggest that correction methods are not required for the station of Istra but are required for the stations of Rijeka and Split.The performance of these correction methods is presented in section 3.

Ensemble regression (ER) / Regresija sklopa
The statistical error metrics in Figure 4 show a moderate improvement of the initial MEDSEA hindcast with the ensemble regression correction method, with a stronger correlation coeffi cient and lower NRMSE, NBIAS, MAPE, and SI.The NRMSE are 0.53 and 0.48 for Rijeka and Split, respectively.This is a 29% and 8% reduction from the initial MEDSEA hindcast for Rijeka and Split, respectively.
The NBIAS was reduced to a negligible value for both Rijeka and Split, with values of 0.01 and 0.04, respectively.Table 4 summarizes the values of the hyperparameters resulting from the Bayesian estimation of the ensemble regression training procedure.While the number of learning cycles and the learning rate is comparable for the Rijeka and Split sites, the values for the minimum leaf size and the maximum number of splits are signifi cantly larger for the Rijeka site.For both sites, the same lasso term lambda is used to prune the constructed ensemble of learning trees to avoid overfi tting (the term λ in Eq. 2).Increasing the lasso term above this value rapidly degrades the predictive power of the ensemble regression for both the training and test sets (not shown in the manuscript for brevity), hence the value 0.002.In summary, the statistical error metrics show the decent performance of the ensemble regression corrected MEDSEA hindcast of the signifi cant wave height with a slight overestimation of the signifi cant wave height.

Artifi cial neural networks (ANN) / Umjetna neuralna mreža (ANN)
In Figure 5, the statistical error metrics show a moderate improvement of the initial hindcast by MEDSEA with ANN as the correction method with higher correlation coeffi cients and lower NRMSE, NBIAS, MAPE, and SI.The NRMSE for the hindcast corrected with ANN is 0.55 and 0.46 for Rijeka and Split, respectively.This represents a 26% and 12% decrease in NRMSE, respectively, compared to the initial MEDSEA hindcast.The NBIAS was completely removed from the hindcast for the Rijeka site, while it was reduced to a negligible value of 0.02 for the Split-site.
Table 5 shows the hyperparameter values of ANN, which performed best in Bayesian optimization after 5-fold crossvalidation on the training data.The sigmoid activation function is used for both Rijeka and Split.The hidden layers are both shallow with only 1 and 2 layers for Split and Rijeka, respectively, and narrow with only a maximum of 17 knots (Rijeka).
Overall, the results of the statistical analysis show a good performance of the corrected hindcast of the signifi cant wave height with ANN with a slight overestimation of the signifi cant wave height for the Split site.

Impact of lag components and time series comparison / Utjecaj komponenti s vremenskim pomakom i usporedba vremenskih serija
Interestingly, there is little difference in the statistical error metrics due to the addition of lag components of wind magnitude from ERA5 and significant wave height from MEDSEA to the predictor set (described in section 2.1) (Table 6).Error metrics varied in the range of ±8%, while some metrics did not change with the inclusion of lag components.
Only the NBIAS changed significantly, as its values were already low before the inclusion of the lag components.Therefore, no significant difference in corrected model performance was found between models with and without lag components.However, the cross-correlation for the Split site calculated between the signals of the measured significant wave height, VHM0, and the lag components of the ERA5 wind magnitude, WIND_MAG, shows that the 1-hour lag component has the highest correlation with the measured VHM0 (0.82).The 0-hour lag component has a slightly lower cross-correlation of 0.81, while the higher lag components steadily decrease in value from the 1-hour lag component.This trend suggests a diffused cross-correlation without a dominant lag component that would have significant explanatory power.The Split location is examined here for cross-validation because in Section 3.2.4,the WIND_ MAG variable shows the highest predictor significance of all predictors used.Figure 6 shows a time series of measured, MEDSEAmodeled, and MEDSEA-corrected significant wave heights at the Rijeka site.The MEDSEA reanalysis model does not respond quickly enough to the increasingly significant wave height on October 21 and eventually underestimates by 0.6 m the largest wave height observed on October 22.The corrected models (ER and ANN) can correct the underestimation, with the ANN model overestimating the largest wave height by up to 0.2 m, but was unable to accelerate the increase in wave height at the beginning of the wave event.The corrected models with the lag components (ER-L and ANN-L) showed no significant improvement over the corrected models without lag components (ER and ANN).Table 6 Statistical error metrics of test data for the correction methods that include lag components for wind magnitude from ERA5 and signifi cant wave height from MEDSEA, as described in section 2.1; the percentage is relative to the error metrics of correction methods excluding the lag components (shown in sections 3.2.1.and 3.2.2) Tablica 6. Statistička metrika grešaka testnih podataka za metode korekcije koja uključuje komponente s vremenskim pomakom za magnitude vjetra od ERA5 i značajne visine vala od MEDESEA, kao što je opisano u sekciji 2.1; postotak je relativan u odnosu na metriku metoda korekcije isključujući komponente zaostajanja (prikazane u sekcijama 3.

Predictor importance / Važnost prediktora
Predictor importance is calculated by summing the changes in node risk due to splits on each predictor variable in a regression tree or ensemble of regression trees and then dividing the sum by the number of branch nodes.The change in node risk is the diff erence between the risk of the parent node and the combined risk for the two child nodes.Node risk is defi ned as the mean squared error of the node weighted by the node probability.In this way, the relative importance of the predictors can be extracted from trained ensemble regression models (Figure 7). Figure 7 shows that the wind magnitude from ERA5 (WIND_ MAG) is the most important predictor for Split, but signifi cant wave height from MEDSEA (VHM0) is the most important predictor for Rijeka.However, signifi cant wave height from MEDSEA (VHM0) still has a signifi cant eff ect on model performance in Split, albeit onefour times smaller than wind magnitude.Other predictors have a smaller eff ect on correction model performance.There could be several reasons for this diff erence.Both sites are located in sheltered basins where the wind is the dominant wave generator.Therefore, the wind strength modeled by ERA5 in Rijeka could be underestimated due to the complex orography, similar to the underestimation of wind magnitude underestimation in the Ligurian Sea in Italy [39].

CONCLUSION / Zaključak
This paper presents a methodology for extending hindcast wave data to sparsely measured locations based on machine learning models and reanalysis data.The advantage of a machine learning model (ANN, ensemble regression, etc.) is that no location-specifi c data are needed for hindcasting wave parameters.Only publicly available global or regional reanalysis model data could be used as input for training and eventually hindcasting wave parameters or fi lling gaps in existing wave measurements.This is particularly important in locations sheltered by complex topography, as nested wave models are typically required to properly represent wave processes from the open ocean to sheltered basins.These nested numerical models require more time to set up and compute.However, unlike Figure 7 Predictor importance in ensemble regression (name designations are described in Table 2) for Split (left) and Rijeka (right) Slika 7. Važnost pojedinog prediktora pri ansambl regresiji (imena parametara opisana su u Tablici 2) za Split i Rijeku (desno) machine learning models, numerical models are constrained by physical properties.Therefore, machine learning models should be used with caution.Because the present technique effi ciently improves local nearshore waves from MEDSEA with low computational cost to better refl ect measured data, it could be applied in locations with sparse wave observations to augment measured wave data.This work has shown that the machine learning hindcast did not improve the Split initial MEDSEA hindcast as well as the Rijeka hindcast.The NRMSE improvement for Rijeka is 29% for the test data (Rijeka-ER-TE), compared to the smaller improvement for Split of 12% (Split-ANN-TE).This smaller improvement for Split may be because the initial wave height hindcasts had higher accuracy, making it diffi cult for the machine learning models to further improve accuracy.Furthermore, the machine learning models reduced the biases in the MEDSEA hindcasts to negligible levels for both Split and Rijeka (NBIAS < 0.03).Nevertheless, the presented machine learning method could not improve the hindcasts for Split and Rijeka to the level of MEDSEA hindcasts for Istra (Figure 2) or other buoys in the open sea of the Adriatic (Table 1).
Interestingly, the results showed that wind duration (via the wind magnitude lag components from -1 h to -10 h) and wave height history (via wave height lag components from -1 h to -10 h) as input data did not signifi cantly improve the performance of the wave height hindcast.With a marginal improvement in statistical error metrics, this is considered a negligible improvement over previously established machine learning models without lag components.This however is not aligned with the observation done with the cross-correlation analysis between the measured signifi cant wave heights and the ERA5 wind magnitudes, where the 1-hour lag component showed the highest cross-correlation.These results also do not agree with those of Peres et al. [8], who found that including up to 18 hours of wind lag components improved wave height hindcast with ANN.However, Peres et al. [8] examined the performance of ANN in the western Mediterranean, i.e., not in extremely sheltered areas with small fetch lengths.MEDSEA greatly miscalculates the wave directions, as can be seen in Figure 3.
Therefore, MEDSEA wave direction should be used with caution in sheltered areas because the spatial resolution of the model is not suffi cient to accurately resolve the wave processes.
In the future, other methods should be explored, such as introducing a weight distribution into the machine learning objective function to increase the penalty for errors at extreme wave heights.This should improve the hindcast for extreme wave heights at the expense of reduced accuracy at lower wave heights.In addition, a random separation between the training and test sets was performed for this study.A carefully performed data separation and curation could lead to even further increases in accuracy.In addition, a hybrid approach of nested numerical wave modeling and machine learning is also a viable option.The physically modeled wave results would serve as an extension of the measured wave data, which would incur higher computational costs to obtain a longer time series on which to train the machine learning models.Finally, measured meteorological data could be used instead of or in addition to the reanalysis data to reconstruct measured wave heights.

Machine learning models / Modeli strojnoga učenja
Three methods were used as correctors for the MEDSEA hindcast data at the Rijeka and Split sites.The MEDSEA hindcast for the Istra site showed a favorable performance metric, so no correction methods were applied to this dataset (described in more detail in Section 3.1.).Sections 2.3.1 and 2.3.2 describe the correction methods used in this study, ensemble regression, and artifi cial neural networks, respectively.Before the training procedure, the fi rst 20% of the measured fi eld data (described in Section 2.2.) were separated and labeled as test data set for testing the trained machine learning model (from November 1, 2007, to December 29, 2007, for Split and from June 1, 2009, to October 30, 2009, for Rijeka).The remaining data (80% of the total measured data) were separated for k-fold training of ensemble regression and artifi cial neural network models in 5 folds.