Multi-Step-Ahead Forecasting of Wave Conditions Based on a Physics-Based Machine Learning (PBML) Model for Marine Operations

Short-term wave forecasts are essential for the execution of marine operations. In this paper, an efficient and reliable physics-based machine learning (PBML) model is proposed to realize the multi-step-ahead forecasting of wave conditions (e.g., significant wave height Hs and peak wave period Tp). In the model, the primary variables in physics-based wave models (i.e., the wind forcing and initial wave boundary condition) are considered as inputs. Meanwhile, a machine learning algorithm (artificial neural network, ANN) is adopted to build an implicit relation between inputs and forecasted outputs of wave conditions. The computational cost of this data-driven model is obviously much lower than that of the differential-equation based physical model. A ten-year (from 2001 to 2010) dataset of every three hours at the North Sea center was used to assess the model performance in a small domain. The result reveals high reliability for one-day-ahead Hs forecasts, while that of Tp is slightly lower due to the weaker implicit relationships between the data. Overall, the PBML model can be conceived as an efficient tool for the multi-step-ahead forecasting of wave conditions, and thus has great potential for furthering assist decision-making during the execution of marine operations.


Introduction
Marine operations such as sea transportation, offshore crane operation and mooring installation are complex operations which are highly weather sensitive. When planning these operations, the response-based operational limits of the systems and equipment are determined and used to estimate the weather window in combination with the long-term statistics of the environmental conditions. However, during the execution phase of the operation, it is important to make short-term forecasts of wave conditions that are characterized by significant wave height (H s ), peak wave period (T p ) and so on. The forecasted wave conditions are compared with the operational sea state limits to decide whether the operation should start or not as long, as the operations last no longer than 72 h [1]. In view of the specific execution position and typical execution duration of marine operations, reliable wave forecasting in a small area from a few hours to several days ahead is necessary.
Currently, wave forecasting methods are generally categorized into three main types: physics-based numerical methods, statistical methods and machine learning methods. Physics-based numerical methods are theory-driven and mimic the behavior of wave evolution based on physical processes such where D Dt is a derivative operator and t is the time. E is the wave energy density spectrum. S is the source function describing different physical processes.
Components in the source item S are the key factor for different physics-based wave models. The widely used and reliable SWAN model will be briefly introduced here to show the modeling process and primary variables of numerical wave models. SWAN is a third-generation spectral wave prediction model developed by the Delft University of Technology in 1999. Its source term includes wave generation; dissipation and quadruplet wave-wave interactions; and dissipation due to bottom friction, triad wave-wave interactions and depth-induced breaking. In this model, the action balance equation is utilized instead of the energy balance equation since the action density N (defined as N = E/σ) is conserved during wave propagation [35]. The action balance equation in the Cartesian coordinate system can be expressed as: where N is an abbreviation for N(x,y,t;σ,θ) with respect to frequency σ and propagation directions θ in space x, y and time t. c x and c y are the propagation velocities of wave energy in spatial x and y space, respectively. c σ and c θ are the propagation velocities in spectral σ and θ space, respectively. S tot is the source term that represents different physical processes, including both deep and shallow water processes. The first item on the left-hand side represents kinematic change in time, while the second and third items denote spatial change. The last two terms depict the refraction effects and shifting of the relative frequency, respectively. The source item, on the right-hand side, consists of six components describing the generation, dissipation or redistribution of wave energy, which are expressed in Equation (3).
S tot = S in + S nl3 + S nl4 + S ds,w + S ds,b + S ds,br where: • S in denotes wave growth by the wind; • S nl3 denotes the nonlinear transfer of wave energy through triad wave-wave interactions; • S nl4 denotes the nonlinear transfer of wave energy through quadruplet wave-wave interactions; • S ds,w denotes wave dissipation due to whitecapping; • S ds,b denotes wave dissipation due to bottom friction; • S ds,br denotes wave dissipation due to depth-induced wave breaking.
Regarding wave forecasting based on the SWAN model, the mean wind speed and direction at 10 m above mean sea level are the key inputs that drive the model runs. Initial and boundary conditions of the wave field are also needed in the spectrum form, such as JONSWAP spectra and Pierson-Moskowitz (PM) spectra. The finite difference method is employed together with the Gauss-Seidel technique [36] for spectral wave computation. After calculation, the wave statistics, such as significant wave height H s , peak wave period T p and mean wave period T m , can be derived by integrating the wave spectrum. In addition, from the split spectrum (i.e., the wind-generated and swell ones), wind and swell contributions can be easily revealed. It should be mentioned that the widely used peak period T p is defined from the energetic concept and is only valid in total waves. For the wind-generated wave or swell, only the mean period T m is available.
As shown above, a physics-based wave model (such as SWAN) can correctly describe the phenomena of the wave formation and propagation by the key inputs and governing equations. However, a physics-based wave model normally requires large computational and storage resources. On one hand, the space-time evolution of the wave energy needs to be solved based on a complex mathematical system, including differential and integral equations, parameterization schemes, highly nonlinear empirical expressions, etc. On the other hand, modeling usually covers a large region or global domain, which increases the storage burden of input and output information.

Machine Learning Algorithms
Recently, machine learning algorithms have drawn more attention in forecast problems as alternative methods based only on time series due to their high cost efficiency. Regarding wave forecasting, a popular type of machine learning method is based on the use of historical wave data only, not relying on the input of wind conditions. Such data can be based on measurements or simulations. It can learn the underlying behavior of a system from a large amount of historical data without having to derive the explicit relationships between them. That is, it can provide results without a governing equation by relying solely on the time history of wave data.
In this study, an artificial neural network (ANN), one of the most popular machine learning algorithms, is utilized due to its ability to capture the nonlinear relations between inputs (past and current wave data) and outputs (forecasted wave data) and to produce multiple outputs simultaneously. It mimics the information processing mechanisms and processes of neurons in the human brain [37]. An ANN consists of three types of layers: the input layer, hidden layers and output layer. A simple ANN with one hidden layer is shown in Figure 1. The first and last layers are called the input and output layers, respectively. They contain n and p neurons representing the input and output variables, respectively. The middle layer connecting the input and output layers is the hidden layer, which contains m neurons. The neuron in each layer receives input information (from the neurons in the previous layer based on the corresponding weight, bias parameter and activation function) and then further propagates it forward to the neurons in the following layer. To illustrate this process more clearly, Equations (4) and (5) show the calculations of the hidden neuron and output neuron, respectively.
where f and g are activation functions of the hidden and output layers, respectively. wij and ujk are weight values from input neuron xi to hidden neuron hj and from hidden neuron hj to output neuron yk, respectively. bj and ck are the biases of the corresponding neurons hj and yk. In this study, the back propagation neural network (BPNN) [38] is adopted, which is a model that fine-tunes the weights of a neural network based on the error rate obtained in the previous iteration. In the training phase, the initial weight and bias parameters are selected arbitrarily. Then, they are adjusted iteratively to reduce the error rate by an optimization algorithm (e.g., the gradient descent method). When the overall error is within an acceptable range, the training process stops, and the optimal neural network is determined. In the subsequent testing phase, the forecasts can be obtained by new inputs and this optimal neural network.
Generally, machine learning is an efficient tool that can handle data flexibly. However, the lack of a clear physical background presents a challenge to the performance and reliability of machine learning algorithms for wave forecast using only the historical data of wave conditions.

The Physics-Based Machine Learning Model
In view of the advantages of physics-based wave models and machine learning algorithms, a hybrid model, termed the physics-based machine learning (PBML) model, is proposed. In the model, physical knowledge from the physics-based wave model is utilized as a guide for designing inputs The first and last layers are called the input and output layers, respectively. They contain n and p neurons representing the input and output variables, respectively. The middle layer connecting the input and output layers is the hidden layer, which contains m neurons. The neuron in each layer receives input information (from the neurons in the previous layer based on the corresponding weight, bias parameter and activation function) and then further propagates it forward to the neurons in the following layer. To illustrate this process more clearly, Equations (4) and (5) show the calculations of the hidden neuron and output neuron, respectively.
where f and g are activation functions of the hidden and output layers, respectively. w ij and u jk are weight values from input neuron x i to hidden neuron h j and from hidden neuron h j to output neuron y k , respectively. b j and c k are the biases of the corresponding neurons h j and y k . In this study, the back propagation neural network (BPNN) [38] is adopted, which is a model that fine-tunes the weights of a neural network based on the error rate obtained in the previous iteration. In the training phase, the initial weight and bias parameters are selected arbitrarily. Then, they are adjusted iteratively to reduce the error rate by an optimization algorithm (e.g., the gradient descent method). When the overall error is within an acceptable range, the training process stops, and the optimal neural network is determined. In the subsequent testing phase, the forecasts can be obtained by new inputs and this optimal neural network.
Generally, machine learning is an efficient tool that can handle data flexibly. However, the lack of a clear physical background presents a challenge to the performance and reliability of machine learning algorithms for wave forecast using only the historical data of wave conditions.

The Physics-Based Machine Learning Model
In view of the advantages of physics-based wave models and machine learning algorithms, a hybrid model, termed the physics-based machine learning (PBML) model, is proposed. In the model, physical knowledge from the physics-based wave model is utilized as a guide for designing inputs and outputs, and the machine learning algorithm is adopted to learn the implicit relationships between them. Afterwards, the trained model can provide wave forecasts directly with new inputs.
The physics-based wave model reveals that the corresponding wind field is the main energy source of the wave field. Thus, wind properties such as the mean speed and direction are introduced as the input variables in the PBML model to forecast wave conditions. In addition, wave boundary conditions at the initial time step are added in the model as well. The wave-wave interaction and wave dissipation effects are implicit in the learned relationship. As a result, the physics-based machine learning (PBML) model is established as follows.

Input and Output Variables
For a typical forecast domain, as shown in the blue part in Figure 2, significant wave height H s and peak period T p at the initial time step at the highlighted grid points (by the red lines) are the initial wave boundary conditions. Wind speeds U w and directions D u at all concerned grid points are the wind forcing inputs. Similar to the physics-based wave model, the wind forcing should cover the past, the current and the whole forecast time horizon, which can be obtained from the wind forecast models such as numerical weather prediction (NWP) models. The future (forecasted) wave conditions (H s and T p ) at all grid points are the outputs. and outputs, and the machine learning algorithm is adopted to learn the implicit relationships between them. Afterwards, the trained model can provide wave forecasts directly with new inputs. The physics-based wave model reveals that the corresponding wind field is the main energy source of the wave field. Thus, wind properties such as the mean speed and direction are introduced as the input variables in the PBML model to forecast wave conditions. In addition, wave boundary conditions at the initial time step are added in the model as well. The wave-wave interaction and wave dissipation effects are implicit in the learned relationship. As a result, the physics-based machine learning (PBML) model is established as follows.

Input and Output Variables
For a typical forecast domain, as shown in the blue part in Figure 2, significant wave height Hs and peak period Tp at the initial time step at the highlighted grid points (by the red lines) are the initial wave boundary conditions. Wind speeds Uw and directions Du at all concerned grid points are the wind forcing inputs. Similar to the physics-based wave model, the wind forcing should cover the past, the current and the whole forecast time horizon, which can be obtained from the wind forecast models such as numerical weather prediction (NWP) models. The future (forecasted) wave conditions (Hs and Tp) at all grid points are the outputs.

Calculation Algorithm
Differently from the action balance equation employed in the physics-based wave model, the PBML model provides an alternative method, directly establishing the input-output relationship through machine learning algorithms. In this process, the computational cost can be significantly decreased due to no need to solve partial differential equations by an iteration process. Several machine learning algorithms are available, such as ANN, ANFIS, RNN, and SVM.

Creating the Forecast Model
After deciding the input and output variables and the calculation algorithm, the forecast model can be established. To illustrate this, a forecast model of Hs is shown as an example in Equation (6), and its architecture is shown in Figure 3.

Calculation Algorithm
Differently from the action balance equation employed in the physics-based wave model, the PBML model provides an alternative method, directly establishing the input-output relationship through machine learning algorithms. In this process, the computational cost can be significantly decreased due to no need to solve partial differential equations by an iteration process. Several machine learning algorithms are available, such as ANN, ANFIS, RNN, and SVM.

Creating the Forecast Model
After deciding the input and output variables and the calculation algorithm, the forecast model can be established. To illustrate this, a forecast model of H s is shown as an example in Equation (6), and its architecture is shown in Figure 3.
[H s1 (t + N), H s2 (t + N), . . . , H sn (t + N)] = f hN H s1 (t), T p1 (t), . . . , H sm (t), T pm (t), U w1 (t + N), D u1 (t + N), . . . , U wn (t + N), D un (t + N), . . . , U w1 (t + N − TU + 1), D u1 (t + N − TU + 1), . . . , U wn (t + N − TU + 1), D un (t + N − TU + 1)) where t is the current time and N is the forecast step. n and m are the number of grid points in the forecast domain and along the boundary lines, respectively. TU is the time range of wind speeds and directions used as wind forcing. The wind input will then include the past and the current values, and the values of the forecast time horizon. where t is the current time and N is the forecast step. n and m are the number of grid points in the forecast domain and along the boundary lines, respectively. TU is the time range of wind speeds and directions used as wind forcing. The wind input will then include the past and the current values, and the values of the forecast time horizon. As displayed in Equation (6) and Figure 3, the outputs of the model are N-step-ahead Hs in the forecast domain. The input variables include the initial wave boundary conditions (Hs and Tp) at time t and the wind forcing (Uw and Du) in the forecast domain from (t+N-TU+1) to (t+N). In addition, fhN indicates the forecast system at N-step-ahead, which depends on the selected machine learning algorithm. In this study, ANN is used.

Training the Established Model
Afterwards, a training process is performed using the training dataset to learn the relationship between inputs and outputs shown in Equation (6). The training process varies according to the selected machine learning algorithm. For instance, for ANN, the weights and bias in the model are updated in each iteration based on the optimization algorithm (e.g., gradient descent algorithm). This process repeats until the stopping criteria are reached.

Evaluating the Performance of the Trained Model
Once training is complete, the model is considered satisfactory and able to make forecasts. By using the inputs in the testing data, the forecasts can be generated from the trained model. Then, the forecast accuracy can be evaluated in terms of several error metrics by comparing forecasted and actual outputs. The metrics used for evaluating the forecast performance are described in detail in Section 2.4.

Forecast Performance Evaluation
After developing a forecast model, it is important to evaluate the model performance. In this study, conventional error measures, such as MAE (mean absolute error), RMSE (root mean square error), bias, R 2 (correlation coefficient), SI (scatter index), MAPE (mean absolute percentage error) and a new proposed forecast error factor are used to quantify the forecasting accuracy. The expressions of the error measures are shown in Equations (7)-(12). As displayed in Equation (6) and Figure 3, the outputs of the model are N-step-ahead H s in the forecast domain. The input variables include the initial wave boundary conditions (H s and T p ) at time t and the wind forcing (U w and D u ) in the forecast domain from (t + N − TU + 1) to (t + N). In addition, f hN indicates the forecast system at N-step-ahead, which depends on the selected machine learning algorithm. In this study, ANN is used.

Training the Established Model
Afterwards, a training process is performed using the training dataset to learn the relationship between inputs and outputs shown in Equation (6). The training process varies according to the selected machine learning algorithm. For instance, for ANN, the weights and bias in the model are updated in each iteration based on the optimization algorithm (e.g., gradient descent algorithm). This process repeats until the stopping criteria are reached.

Evaluating the Performance of the Trained Model
Once training is complete, the model is considered satisfactory and able to make forecasts. By using the inputs in the testing data, the forecasts can be generated from the trained model. Then, the forecast accuracy can be evaluated in terms of several error metrics by comparing forecasted and actual outputs. The metrics used for evaluating the forecast performance are described in detail in Section 2.4.

Forecast Performance Evaluation
After developing a forecast model, it is important to evaluate the model performance. In this study, conventional error measures, such as MAE (mean absolute error), RMSE (root mean square error), bias, R 2 (correlation coefficient), SI (scatter index), MAPE (mean absolute percentage error) and a new proposed forecast error factor are used to quantify the forecasting accuracy. The expressions of the error measures are shown in Equations (7)- (12).
where n is the number of data points in the testing dataset. a i and f i are ith actual and forecasted data, respectively. a and f are the mean values of the actual and forecasted data, respectively. A new proposed error factor ε M shown in Equation (13), which is defined by the difference between the forecasted and the actual value dividing the actual value, can reveal the forecast uncertainty at each time step. More information on evaluating multi-step-ahead forecasts can be found in reference [27].
where f (t) and a(t) are the forecasted and actual values, respectively, with a forecast time horizon of t.

The Case Study Area and Datasets
The North Sea is considered an area with great potential for offshore wind energy development [39]. In the present study, the North Sea center is the main focus, as highlighted in Figure 4. In this area, various marine operations have been studied and analyzed, such as the installation of offshore wind turbine monopiles and transition pieces [40], and single turbine blade [41][42][43][44]. Determining suitable weather windows based on the wave forecasts is a major challenge for decision-making during the execution of these operations. A ten-year dataset (from 2001 to 2010) of waves and wind were extracted from CERA-20C with an interval of 3 hours and a spatial resolution of 0.125° (about 13.7 km) in both longitude and latitude. The first nine years of data (from 2001 to 2009) were considered as the training data for establishing and validating the model, while the last year (2010) of data were considered the testing data for evaluating the forecast performance of the model. Regarding the typical duration of marine operations, the aim of this study was to perform one-day-ahead wave forecasts, that is, eight-stepahead forecasts for the selected data. Given that marine operations are normally carried out at a selected position, in this study, the focus is on a small area. As shown in Figure 5, the forecast domain contains nine grid points numbered from 1 to 9. .

Model Architecture Determination
It is important to determine the model architecture (including the ANN structure and the time range of wind forcing) based on the specific dataset samples before using it. Random selection may result in poor forecast performance. Therefore, the objective of this section is to determine the optimal architecture for the PBML model. For a typical ANN, three key parameters, namely, the number of hidden layers, the number of neurons in each layer and the type of activation function, need to be carefully selected to guarantee a good performance. Regarding the optimal number of hidden layers and neurons, there are no fixed selection formulas, but some rules of thumb are available to determine their scope [49,50]. In the sensitivity study, five different combinations of hidden layers and neurons In this study, the detailed site-specific wave and wind hindcast data extracted from the ECWMF CERA-20C datasets [31,45] were used. Several verification studies [31,[46][47][48] have indicated that the generated data from ECMWF are consistent with the measurement data and suitable for climate studies. However, the same method can be applied when measurement data of both wind input and wave output are available. According to the requirements of the PBML model (described in Section 2.3), the wind direction D u and wind speed U w at 10m height, significant wave height H s and peak wave period T p , are necessary.
A ten-year dataset (from 2001 to 2010) of waves and wind were extracted from CERA-20C with an interval of 3 hours and a spatial resolution of 0.125 • (about 13.7 km) in both longitude and latitude. The first nine years of data (from 2001 to 2009) were considered as the training data for establishing and validating the model, while the last year (2010) of data were considered the testing data for evaluating the forecast performance of the model. Regarding the typical duration of marine operations, the aim of this study was to perform one-day-ahead wave forecasts, that is, eight-step-ahead forecasts for the selected data. Given that marine operations are normally carried out at a selected position, in this study, the focus is on a small area. As shown in Figure 5, the forecast domain contains nine grid points numbered from 1 to 9.  A ten-year dataset (from 2001 to 2010) of waves and wind were extracted from CERA-20C with an interval of 3 hours and a spatial resolution of 0.125° (about 13.7 km) in both longitude and latitude. The first nine years of data (from 2001 to 2009) were considered as the training data for establishing and validating the model, while the last year (2010) of data were considered the testing data for evaluating the forecast performance of the model. Regarding the typical duration of marine operations, the aim of this study was to perform one-day-ahead wave forecasts, that is, eight-stepahead forecasts for the selected data. Given that marine operations are normally carried out at a selected position, in this study, the focus is on a small area. As shown in Figure 5, the forecast domain contains nine grid points numbered from 1 to 9. .

Model Architecture Determination
It is important to determine the model architecture (including the ANN structure and the time range of wind forcing) based on the specific dataset samples before using it. Random selection may result in poor forecast performance. Therefore, the objective of this section is to determine the optimal architecture for the PBML model. For a typical ANN, three key parameters, namely, the number of hidden layers, the number of neurons in each layer and the type of activation function, need to be

Model Architecture Determination
It is important to determine the model architecture (including the ANN structure and the time range of wind forcing) based on the specific dataset samples before using it. Random selection may result in poor forecast performance. Therefore, the objective of this section is to determine the optimal architecture for the PBML model. For a typical ANN, three key parameters, namely, the number of hidden layers, the number of neurons in each layer and the type of activation function, need to be carefully selected to guarantee a good performance. Regarding the optimal number of hidden layers and neurons, there are no fixed selection formulas, but some rules of thumb are available to determine their scope [49,50]. In the sensitivity study, five different combinations of hidden layers and neurons are considered. In addition, three common activation functions are used for each combination, one of which is the linear function purelin (MATLAB®), and the other two are the nonlinear functions log-sigmoid (logsig, MATLAB®) and hyperbolic tangent sigmoid (tansig, MATLAB®). They are displayed in Figure 6. For point 5 shown in Figure 5, the one-step-ahead forecast results of H s are summarized in Table 1, where error metrics RMSE and R 2 are used to measure the model's performance. It should be noted that for the sake of simplicity, all results were obtained through the PBML model with TU = 1 (see Equation (6)). are considered. In addition, three common activation functions are used for each combination, one of which is the linear function purelin (MATLAB®), and the other two are the nonlinear functions logsigmoid (logsig, MATLAB®) and hyperbolic tangent sigmoid (tansig, MATLAB®). They are displayed in Figure 6. For point 5 shown in Figure 5, the one-step-ahead forecast results of Hs are summarized in Table 1, where error metrics RMSE and R 2 are used to measure the model's performance. It should be noted that for the sake of simplicity, all results were obtained through the PBML model with TU=1 (see Equation (6)).   Table 1 that the tansig activation function is better than the others, as the RMSE is lower and R 2 is higher for all cases. However, the numbers of hidden layers and neurons have little influence on the forecast performance in all cases with the tansig activation function. By comparison, case 9 obtained the best result. Therefore, a neural network with three hidden layers, twenty neurons in each layer and the tansig type activation function was adopted for the PBML model in this study.  It is clearly seen in Table 1 that the tansig activation function is better than the others, as the RMSE is lower and R 2 is higher for all cases. However, the numbers of hidden layers and neurons have little influence on the forecast performance in all cases with the tansig activation function. By comparison, case 9 obtained the best result. Therefore, a neural network with three hidden layers, twenty neurons in each layer and the tansig type activation function was adopted for the PBML model in this study.
In addition to the ANN structure discussed above, the time range of the input wind field also deserves consideration. As shown in Equation (6), when performing the N-step-ahead forecasting, the wind forcing from (t + N − TU + 1) to (t + N) can be used as input. To investigate the impact of wind forcing at different time steps on H s forecasts, the ANN was developed in terms of the optimal structure obtained above, and the TU value was selected from one to eight. When TU is equal to 1, only the wind field at the same time as the forecasted H s is considered. With the increase in TU, an increasing number of past wind fields are taken into account for wave forecasting. The results of the sensitivity study of the time range (TU) of input wind forcing are shown in the Figure 7.  As seen in Figure 7, in general, the forecast performance will increase with increasing TU values, indicating that the wind variation should be taken into account in wave forecasting by the PBML model. However, the level of improvement in model performance gradually decreases as the TU increases from 2. Meanwhile, increasing the TU value will increase the number of inputs, which makes the forecast model more complicated and requires more computational time. Therefore, TU=2 was chosen in the PBML model in this study. That is, the dynamic characteristics of the winds are provided in the model by applying wind speeds and directions at two consecutive time steps.
In addition, it is interesting to see the importance of the future wind forcing. Figure 8 shows the forecast performance of the model using wind forcing up to (t+N-TU*) when forecasting wave conditions at (t+N), where TU* ranges from 0 to 7. From the sensitivity study, it is clear that the forecast performance will gradually decrease as TU* increases. Therefore, in order to accurately forecast the wave conditions at (t+N), the wind conditions at (t+N) should be included. As seen in Figure 7, in general, the forecast performance will increase with increasing TU values, indicating that the wind variation should be taken into account in wave forecasting by the PBML model. However, the level of improvement in model performance gradually decreases as the TU increases from 2. Meanwhile, increasing the TU value will increase the number of inputs, which makes the forecast model more complicated and requires more computational time. Therefore, TU = 2 was chosen in the PBML model in this study. That is, the dynamic characteristics of the winds are provided in the model by applying wind speeds and directions at two consecutive time steps.
In addition, it is interesting to see the importance of the future wind forcing. Figure 8 shows the forecast performance of the model using wind forcing up to (t + N − TU*) when forecasting wave conditions at (t + N), where TU* ranges from 0 to 7. From the sensitivity study, it is clear that the forecast performance will gradually decrease as TU* increases. Therefore, in order to accurately forecast the wave conditions at (t + N), the wind conditions at (t + N) should be included. provided in the model by applying wind speeds and directions at two consecutive time steps.
In addition, it is interesting to see the importance of the future wind forcing. Figure 8 shows the forecast performance of the model using wind forcing up to (t+N-TU*) when forecasting wave conditions at (t+N), where TU* ranges from 0 to 7. From the sensitivity study, it is clear that the forecast performance will gradually decrease as TU* increases. Therefore, in order to accurately forecast the wave conditions at (t+N), the wind conditions at (t+N) should be included.

Significant Wave Height
(1) Small domain In this subsection, a PBML model is developed for H s in a small domain given in Figure 5. The model architecture is set up based on the analysis in Section 3.2. Consequently, the forecast model of H s can be expressed as follows: where H s1 (t + N), H s2 (t + N), . . . , H s9 (t + N) are the outputs, which are N-step-ahead forecasts of H s for points 1 to 9. f hN is the ANN model at the N-step-ahead. The inputs of f hN include the initial wave boundary wave conditions (i.e., H s (t) and T p (t) at points 1, 2, 3, 4, 6, 7, 8 and 9) and the wind forcing at two consecutive time steps in the forecast domain (i.e., U w and D u for all nine points at N-and N−1-step-ahead). Thus, it can be seen that in this case, the model for each forecast step has a total of 52 inputs and 9 outputs. The forecast results during the testing period (one year) in the grid, point 5 (in Figure 5), are plotted in Figure 9, in which the actual and forecasted time series are represented by black and red lines, respectively. In addition, the forecast errors at each forecast lead time are summarized in Table 2. To reveal the forecast performance more clearly, two short-time windows are also shown in Figure 9. It is visible from the windows that each forecast contains the forecasting of eight points, from one-step-ahead to eight-step-ahead forecasts. The green and blue points represent one-step-ahead and eight-step-ahead forecast values, respectively, which mean the beginning and end of one forecast. The next forecast starts at the time when the previous forecast horizon ends. lines, respectively. In addition, the forecast errors at each forecast lead time are summarized in Table  2. To reveal the forecast performance more clearly, two short-time windows are also shown in Figure  9. It is visible from the windows that each forecast contains the forecasting of eight points, from onestep-ahead to eight-step-ahead forecasts. The green and blue points represent one-step-ahead and eight-step-ahead forecast values, respectively, which mean the beginning and end of one forecast. The next forecast starts at the time when the previous forecast horizon ends.   Overall, a good agreement between the forecasted and actual data is observed in Figure 9, which reveals the reliability of the PBML model. In addition, there is no significant difference in forecast performance between the two windows. This implies that the proposed PBML model is suitable for any time of year. In Table 2, one can see that the forecast accuracy decreases with the lead time, which is reflected in the decrease in R 2 and the increase in other error metrics. Nevertheless, all 8-step-ahead forecasts can be considered satisfactory because the corresponding errors are quite small. For example, the RMSE and R 2 for all steps are less than 0.3 m and higher than 96%, respectively. Furthermore, regarding the computational time, the PBML model can be quickly executed to train and make forecasts extremely fast. The training process of the PBML model on a common server typically takes few minutes.
The above results only display the forecast at the grid point 5 in the forecast domain. Figure 10 further summarizes the error measures at all nine grid points in the forecast domain done by different colors. Little difference in all error measures can be found regardless of the grid points, which indicates that the PBML model is suitable and reliable in a small forecast domain. and make forecasts extremely fast. The training process of the PBML model on a common server typically takes few minutes.
The above results only display the forecast at the grid point 5 in the forecast domain. Figure 10 further summarizes the error measures at all nine grid points in the forecast domain done by different colors. Little difference in all error measures can be found regardless of the grid points, which indicates that the PBML model is suitable and reliable in a small forecast domain.

2) Single position
As shown above, the proposed PBML model exhibits good performance for Hs forecasting in a small domain. In this subsection, the forecast performance of the PBML model for a single point is investigated. This situation is very important to explore whether the PBML model is applicable to the measurement data, since measurement data may only be available at limited positions. In practice, we cannot expect detailed measured data for a large region (which are normally available in reanalysis or hindcasting).
Since the corresponding input information is greatly reduced, the forecast model is simplified to Equation (15). In this case, the inputs are the initial wave conditions and wind forcing at the position of interest (point 5 in Figure 5 is selected as an example here), and the outputs are the future wave conditions at the same position. That is, the model is simplified to have six inputs and one output. The results of forecasted time series and forecast errors are illustrated in Figures 11 and 12, respectively.
, , , (2) Single position As shown above, the proposed PBML model exhibits good performance for H s forecasting in a small domain. In this subsection, the forecast performance of the PBML model for a single point is investigated. This situation is very important to explore whether the PBML model is applicable to the measurement data, since measurement data may only be available at limited positions. In practice, we cannot expect detailed measured data for a large region (which are normally available in reanalysis or hindcasting).
Since the corresponding input information is greatly reduced, the forecast model is simplified to Equation (15). In this case, the inputs are the initial wave conditions and wind forcing at the position of interest (point 5 in Figure 5 is selected as an example here), and the outputs are the future wave conditions at the same position. That is, the model is simplified to have six inputs and one output. The results of forecasted time series and forecast errors are illustrated in Figures 11 and 12, respectively.  Compared with the results in Figure 9, the forecast accuracy is reduced to some extent for this simple model due to the significant reduction of the input information. Nevertheless, this model is still acceptable since it exhibits good forecast performance. By observing the forecast errors in Figure  12, one can find that the R 2 and RMSE of all 8-step-ahead forecasts are higher than 94% and lower than 0.33 m, respectively. The results imply that the proposed PBML model can forecast Hs with relatively high accuracy based on the measured data of a single position.

Peak Wave Period
In addition to Hs, the forecast performance of the PBML model for Tp is investigated in this section. Similarly to Equation (14), the corresponding model of Tp is expressed in Equation (16), in which the initial wave boundary conditions and spatial wind forcing at two consecutive time steps  Compared with the results in Figure 9, the forecast accuracy is reduced to some extent for this simple model due to the significant reduction of the input information. Nevertheless, this model is still acceptable since it exhibits good forecast performance. By observing the forecast errors in Figure  12, one can find that the R 2 and RMSE of all 8-step-ahead forecasts are higher than 94% and lower than 0.33 m, respectively. The results imply that the proposed PBML model can forecast Hs with relatively high accuracy based on the measured data of a single position.

Peak Wave Period
In addition to Hs, the forecast performance of the PBML model for Tp is investigated in this section. Similarly to Equation (14), the corresponding model of Tp is expressed in Equation (16), in which the initial wave boundary conditions and spatial wind forcing at two consecutive time steps Compared with the results in Figure 9, the forecast accuracy is reduced to some extent for this simple model due to the significant reduction of the input information. Nevertheless, this model is still acceptable since it exhibits good forecast performance. By observing the forecast errors in Figure 12, one can find that the R 2 and RMSE of all 8-step-ahead forecasts are higher than 94% and lower than 0.33 m, respectively. The results imply that the proposed PBML model can forecast Hs with relatively high accuracy based on the measured data of a single position.

Peak Wave Period
In addition to H s , the forecast performance of the PBML model for T p is investigated in this section. Similarly to Equation (14), the corresponding model of T p is expressed in Equation (16), in which the initial wave boundary conditions and spatial wind forcing at two consecutive time steps are adopted as the inputs. The setting of the ANN architecture is consistent with that in Section 3.1. The forecast results are shown in Figure 13, and the errors are summarized with respect to lead times in Table 3. are adopted as the inputs. The setting of the ANN architecture is consistent with that in Section 3.1.
The forecast results are shown in Figure 13, and the errors are summarized with respect to lead times in Table 3.  In general, PBML model can obtain acceptable Tp forecasts, as the R 2 and RMSE at all lead times are higher than 71% and lower than 1.30 s respectively. However, poor performance can normally be observed at the points on steep gradients (e.g., data near 130 and 1200 in Figure 13), which is less reliable than the forecast of Hs.
Significant fluctuation in Tp might be the key reason for this phenomenon. Differently from Hs, Tp is a wave period associated with the highest energetic waves in the total wave spectrum, and it can be predominated by either wind-generated waves or swells. This uncertainty increases the instability of the series of Tp and further challenges the forecast performance.  In general, PBML model can obtain acceptable T p forecasts, as the R 2 and RMSE at all lead times are higher than 71% and lower than 1.30 s respectively. However, poor performance can normally be observed at the points on steep gradients (e.g., data near 130 and 1200 in Figure 13), which is less reliable than the forecast of H s .
Significant fluctuation in T p might be the key reason for this phenomenon. Differently from H s , T p is a wave period associated with the highest energetic waves in the total wave spectrum, and it can be predominated by either wind-generated waves or swells. This uncertainty increases the instability of the series of T p and further challenges the forecast performance.
To further explore the factors that affect the performance, in Section 3.4, the total waves are separated into wind-generated waves and swells, and the PBML model is utilized to forecast wave conditions of these two components independently.

Forecast of Separate Wind-Generated Wave Conditions and Swell Conditions
In general, a total wave comprises two components, i.e., the wind-generated wave and the swell. Approximately, the spectral component is considered to be subject to forcing by the wind when the following condition is met [51]: where u* is the friction velocity. c = c(σ, h) is the wave celerity, as derived from the linear theory of waves, in which σ is wave frequency and h is water depth. θ and φ are the wave and wind directions, respectively. Based on Equation (17), the wind-generated wave and swell can be separated, and the corresponding spectra and wave statistics can be obtained.
In this part, the data of wind-generated waves and swells were also extracted from CERA-20C with the same forecast domain, resolution and duration as those of the total waves. The extracted variables include the significant heights of wind-generated waves H sw and swells H ss ; the corresponding mean periods T mw and T ms ; and the mean wind speed U w and direction D u . The forecast performances of mean wave period (T mw and T ms ) and significant wave height (H sw and H ss ) are discussed in the following subsections, respectively.

Mean Wave Period
The PBML models of T mw and T ms are similar to Equation (14), except for the initial wave condition. In addition, the identical procedure was adopted, which utilized the data from the first nine-year and the last year to train and test the model, respectively. The corresponding results of point 5 (in Figure 5) are shown in Figures 14 and 15 as an example, given the similar forecast performance for all nine points revealed in Section 3.3.1. Meanwhile, Tables 4 and 5 give the forecast errors of T mw and T ms with respect to the lead time, respectively. To further explore the factors that affect the performance, in Section 3.4, the total waves are separated into wind-generated waves and swells, and the PBML model is utilized to forecast wave conditions of these two components independently.

Forecast of Separate Wind-Generated Wave Conditions and Swell Conditions
In general, a total wave comprises two components, i.e., the wind-generated wave and the swell. Approximately, the spectral component is considered to be subject to forcing by the wind when the following condition is met [51]: where u* is the friction velocity. c=c(σ, h) is the wave celerity, as derived from the linear theory of waves, in which σ is wave frequency and h is water depth. θ and ϕ are the wave and wind directions, respectively. Based on Equation (17), the wind-generated wave and swell can be separated, and the corresponding spectra and wave statistics can be obtained.
In this part, the data of wind-generated waves and swells were also extracted from CERA-20C with the same forecast domain, resolution and duration as those of the total waves. The extracted variables include the significant heights of wind-generated waves Hsw and swells Hss; the corresponding mean periods Tmw and Tms; and the mean wind speed Uw and direction Du. The forecast performances of mean wave period (Tmw and Tms) and significant wave height (Hsw and Hss) are discussed in the following subsections, respectively.

Mean Wave Period
The PBML models of Tmw and Tms are similar to Equation (14), except for the initial wave condition. In addition, the identical procedure was adopted, which utilized the data from the first nine-year and the last year to train and test the model, respectively. The corresponding results of point 5 (in Figure 5) are shown in Figures 14 and 15 as an example, given the similar forecast performance for all nine points revealed in Section 3.3.1. Meanwhile, Tables 4 and 5 give the forecast errors of Tmw and Tms with respect to the lead time, respectively.     Figure 14 shows good agreement between the forecasted and actual data of wind-generated waves. The small mean (0.003) and standard deviation (0.072) of the forecast error factor shown in Table 4 further verify the high reliability of the PBML model, even at the last forecast step (8-stepahead/24 h). For the swells, Figure 15 presents a result with less confidence; and a slightly higher mean (0.008) and standard deviation (0.096) of the forecast error factor at eight-steps-ahead also reveal this phenomenon in Table 5. The weak relation between the output (swells) and inputs (wind forcing) is the main reason for this forecast performance. To discuss this phenomenon in detail, Figure  16 shows scatter diagrams of the mean wind speed and different wave periods (i.e., Tp, Tmw and Tms) at the same time and position. The data from year 2001 to 2009 were used to investigate the correlations-about 26,000 data points. The color depicts the absolute frequencies of occurrence. Obviously, the correlation between the mean wind speed and mean wave period for swells ( Figure  16c) and total waves (Figure 16a) is much weaker than that for wind-generated waves (Figure 16b Figure 14 shows good agreement between the forecasted and actual data of wind-generated waves. The small mean (0.003) and standard deviation (0.072) of the forecast error factor shown in Table 4 further verify the high reliability of the PBML model, even at the last forecast step (8-step-ahead/24 h). For the swells, Figure 15 presents a result with less confidence; and a slightly higher mean (0.008) and standard deviation (0.096) of the forecast error factor at eight-steps-ahead also reveal this phenomenon in Table 5. The weak relation between the output (swells) and inputs (wind forcing) is the main reason for this forecast performance. To discuss this phenomenon in detail, Figure 16 shows scatter diagrams of the mean wind speed and different wave periods (i.e., T p , T mw and T ms ) at the same time and position.
The data from year 2001 to 2009 were used to investigate the correlations-about 26,000 data points. The color depicts the absolute frequencies of occurrence. Obviously, the correlation between the mean wind speed and mean wave period for swells ( Figure 16c) and total waves (Figure 16a) is much weaker than that for wind-generated waves (Figure 16b) since the swell stems from the propagating waves generated far away from the focused position. The time variant propagating velocity and direction increase the uncertainty of the swells and further increase the forecast complexity of T p . In contrast, wind-generated waves stemming from the local wind can be exactly described as the wind forcing used as inputs. since the swell stems from the propagating waves generated far away from the focused position. The time variant propagating velocity and direction increase the uncertainty of the swells and further increase the forecast complexity of Tp. In contrast, wind-generated waves stemming from the local wind can be exactly described as the wind forcing used as inputs.

Significant Wave Height
In line with the study of peak wave period Tp, significant wave height is also treated separately by the sources of wind-generated waves (Hsw) and swells (Hss), to further prove the findings observed in the forecast of wave period. The PBML models of Hsw and Hss were established respectively, and the corresponding forecast results of point 5 (in Figure 5) are shown in Figures 17 and 18. Meanwhile, Tables 6 and 7 summarize the forecast errors of Hsw and Hss, respectively.

Significant Wave Height
In line with the study of peak wave period T p , significant wave height is also treated separately by the sources of wind-generated waves (H sw ) and swells (H ss ), to further prove the findings observed in the forecast of wave period. The PBML models of H sw and H ss were established respectively, and the corresponding forecast results of point 5 (in Figure 5) are shown in Figures 17 and 18. Meanwhile, Tables 6 and 7 summarize the forecast errors of H sw and H ss , respectively.          Obviously, Figures 17 and 18 present similar findings to those illustrated in Figures 14 and 15, namely, that the forecast performance of the PBML model for wind-generated waves is better than that for swells. Specifically, the PBML model can generate H sw forecasts accurately and the forecast uncertainty is quite low. For example, Table 6 shows that the RMSE and R 2 for all steps are less than 0.18 m and higher than 98%, respectively. In contrast, the RMSE and R 2 for H ss forecasts are less than 0.45 m and higher than 70% for all steps respectively, as can be seen in Table 7.
In summary, limited knowledge on the source of swells persuaded us to stop at this stage, since merely increasing the control area and inputs would contribute little. In view of the swell being related to the waves propagating beyond the zones of their generation, the spatial domain of inputs is expected to be expanded, and more related input variables such as wave direction, speed and the distance from the position of the input to the output, are expected to be involved to improve the PBML model for T p forecasting in future studies.

Conclusions
In this study, a physics-based machine learning (PBML) model was proposed and developed in detail to do multi-step-ahead forecasting of wave conditions for marine operations. The model attempts to combine the advantages of both physics-based wave models and machine learning algorithms. Specifically, the PBML model takes physical background into account by applying the primary inputs (i.e., initial wave boundary conditions and the future wind forcing) in physics-based models to design the model structure. A machine learning algorithm (the ANN considered in this study is an example) is employed to train the PBML model based solely on historical data.
The PBML model attempts to provide wave forecasts for marine operations, and high computational efficiency is the primary advantage compared with widely used physics-based wave models. On one hand, the PBML model is capable of using the wind forcing parameters in a small area or even from a single position as input to generate forecasts of wave statistics, which do not need the real-time simulations for large domains using directly physics-based numerical wave models. On the other hand, the training process of machine learning is performed mainly through matrix operations rather than solving partial differential control equations. Furthermore, once the model is trained based on the site-specific data, making forecasts on new inputs is straightforward and extremely fast. Based on the forecasted wave conditions and the environmental limits of a specific marine operation, suitable weather windows can be identified, which will further assist decision-making during the execution phase of the operation.
In this study, numerical wave models were used to generate both historical input and historical output data for training the machine learning algorithms. The developed algorithms are very specific for the location in consideration. However, the same methodology can be applied if the data from field measurements or the data from other locations are used.
To assess the performance of the proposed model, three-hourly time series of wave and wind data (such as wind direction D u , mean wind speed U w , significant wave height H s and peak wave period T p ) from 2001 to 2010 at the North Sea center from CERA-20C were used. As a new hybrid model, the optimal architecture of the PBML model was determined first. Afterwards, eight-step-ahead forecasts of H s and T p were generated individually and evaluated by commonly used error measures (MAE, RMSE, SI and R 2 ) and the newly defined forecast error factor. The results indicate that the forecast performance of H s is extremely good, while that of T p is slightly worse due to the weaker implicit relationship between the selected input and output; that was further proved by forecasting wave conditions of wind-generated waves and swells separately. Overall, the PBML model can be conceived as an efficient calculation model for wave forecasting in a small area, which is particularly suitable for marine operations. The forecast uncertainty of the peak period might be further decreased by including some more related factors, such as wave direction, speed and the distance from the position of the input to the output in the PBML model. A feature selection technique such as neighborhood component analysis (NCA) will be adopted to select and weight the related factors. Moreover, how to account for the effect of weather forecast uncertainty on marine operations is a key issue in practice. A suitable method for addressing wave forecast uncertainties in the planning and execution of marine operations can also be expected in the future publications.