Exploring LSTM neural networks for sea level forecasting based on two-decadal tidal and meteorological data

2.1 Analytical approach of sea level estimation

The Intergovernmental Oceanographic Commission (IOC) suggests that every single sea level recording should be considered as the summation of the Mean Sea Level (MSL) component plus the astronomical tide plus the storm surge (meteorological tide) elements, according to Equation (1) [19].

where, SL is the sea level recording, Zo is the MSL component, AT is the astronomical tide, SS is the storm surge and res refers to the sea level residual signal which includes any unmodeled influence (i.e. river runoff, wave interactions, etc.).

Astronomical tide is computed using Equation (2) [20].

where, h(t) is the sea level recording at time t, Zo is the MSL component, a is the coefficient of sea level linear trend, f _j and u _j are coefficients of amplitude and phase of nodal corrections for tidal constituent j, h _j is the amplitude of tidal constituent j, V _j is the equilibrium argument in relation to time t ₀, and g _j is the phase of tidal constituent j.

According to the 2-dimension barotropic model, storm surges depend mostly on local atmospheric conditions of air pressure and wind and is described by the equations of mass and momentum conservation [21], as in Equations (3).

where, (M, N) are the vertical integrals of the wind speed components in the horizontal level (x, y), f is the Coriolis parameter, g is the acceleration of gravity, D is the mean depth of the area, ζ is the sea surface height variation due to storm surges, ρ is the sea water density, Pa is the atmospheric pressure at sea level and (τ _x , τ _y ) are the wind stress components in (x, y) directions; x is the across and, y the along to the coast directions respectively.

Wind stress is the shear force of the 2-d barotropic model, parametrized by the wind speed at a constant height (usually 10 m) above sea surface. The across the coast components is calculated as τ x = ρ air C d W x W x 2 + W y 2 , where ρ _air is the air density and C _d is the drag coefficient between atmospheric air and the sea surface. The along the coast component is calculated by using the respective W _y component of wind speed before the square root.

Wind speed is also involved in the expressions of Equation (3); nevertheless, the component inducing most of the coastal sea surface elevation is the one perpendicular to the coastline. Let W _x be the wind speed perpendicular (across) to the coast, then it is calculated as in. Equation (4).

where, W is the wind velocity, az is the coast azimuth and ω is the wind direction.

Non-tidal residual (NTR) is the remaining value in each time instant t if the astronomical tide is subtracted by the sea level signal, as in Equation (5).

The above analytical approach of tide and storm surge analysis has been successfully implemented in the study area of Thermaikos Gulf, Greece [22], 23]. The astronomical tide computed implementing UTide software, including the constant and the linear trend components as part of the analysis and storm surges by in-house scripts developed in Matlab^®.

Recently, in addition to traditional analytical methods, the evolution in data analytics and the increase of computing power have raised the interest of many researchers studying physical phenomena described in timeseries data using various forms of Artificial Intelligence (AI) techniques, including Machine Learning (ML), Artificial Neural Networks (ANN), and more recently Deep Learning (DL) techniques.

2.2 Machine learning methods in timeseries forecasting

Contrary to standard cases of classification and regression, long-standing timeseries forecasting requires special data management in model fitting considering the additional complexity induced due to time dependance in observables. ML methods, based on the Artificial Neural Networks (ANN) principle, can be rather effective in timeseries forecasting problems as they treat multiple variables with proper manipulation, including the non-linear dependence in data and time, as well as incomplete data recordings. However, certain prerequisites necessitate to achieve best results, from which proper data pre-processing and cleaning are the most important [24]. Due to the ANN’s capability to develop non-linear models, they have evolved, throughout time, to a reliable forecasting tool for hydro-climatic sciences, water resources, physical processes in various time-scales and for a wide range of environmental applications [2].

Recurrent Neural Networks (RNN) approaches, such as the Long Short-Term Memory (LSTM) algorithm compute a projection function from input to output data. These networks benefit from the explicit manipulation of the time sequence in the data during the training stage of the projection function. Particularly, this type of ANN supports data recordings comprised from individual observation series. In effect, every additional variable time series contributes an additional dimension for the estimating function. That means that the network can learn a projection function from input to output data, over time [25]. This unique feature of LSTM has been successfully applied in problems of natural language processing, including text translation and for timeseries forecasting. Furthermore, an asset of RNNs is their capability to identify the data time dependence during training. This temporal dependence consists a very powerful feature of LSTM networks, so that the most relevant input data frame to the output data is identified, trained and can be dynamically modified. In practice, when an observation enters, the network is trained considering this observation in relevance to the past ones, and these are related to the future values as a manner of prediction. The model is trained to project the input to the output, but also, to distinguish which input data frame is useful for the projection. Because of the ability of LSTM to handle long-term correlations in a data series, the need for a predefined time window is neglected providing the capability to model complex timeseries of multiple variables [26].

At a programming level, forecasting based on LSTM networks usually implements two formulations. A static method or closed loop forecasting and a dynamic one or open loop forecasting. The static method predicts subsequent time steps of a sequence based on previous predictions as input. Thus, the model does not require true values of each time step to predict the next one. In practice, to predict the values from a time instant t to time instant t + k using values in a sequence between 1 and t − 1, for each prediction of time step i, the value of the predicted step i − 1 is used as input. This method is suitable for predicting subsequent multiple time steps for the case there are no available true data values to feed the network at pre-prediction stage. In contrast, the dynamic method predicts the next time step using only input data. Hence, in a sequence containing data values from time step 1 to t − 1, only time step t is predicted. To predict the value for time step t + 1, the network should wait until the true values of time step t is recorded for use as input. For the case of a dynamic method, the existence of observations is essential to serve as input in the RNN before each time step prediction. An illustration of closed and open loop forecasting is depicted in Figures 1 and 2.

Figure 1:

Closed loop (multistep) technique for univariate forecasting.

Figure 2:

Open loop (dynamic) technique for multivariate to univariate forecasting.

2.3 State-of-the-art in AI-based sea level estimation

A thorough investigation of the last-few-years published literature confirms that sea level forecasting studies rely extensively on ML and DL techniques. ML techniques have been used in several timescales for sea level prediction. Specifically, Extreme Level Machine (ELM) and Relevance Vector Machine (RVM) techniques based on daily mean sea level values have been used as input for a prediction of 1, 3 and 7 days ahead in two tide gauges in Taiwan [27]. Similarly, Roshni et al. [28] used sea level monthly averages from the PSMSL (Permanent Service for the Mean Sea Leve) database and climate data such as air temperature, surface air pressure, humidity and rainfall for the same period originating from NCEP/NCAR Reanalysis project (National Centers for Environmental Prediction/National Center for Atmospheric Research) for Haldia port, India. They compared ML methods, such as RVM, ELM and GPR (Gaussian Process Regression) with the latter to outperform against the others. In another study, Guillou and Chapalain [29] compare alternative ANN solutions using 3-year maxima during high tides together with coefficients of astronomical tide, atmospheric pressure, wind speed and river runoff at Elon estuary, France. In return, the objective was sea level maxima prediction suggesting that the MLP (Multilayer Perceptron) methodology results in best performance. Moreover, ANNs have been used for harbor design studies; for instance, for wave transmission coefficient prediction using experimental data [30].

Sea level modelling and forecasting based on a monthly temporal resolution is quite widespread in published literature. Sun et al. [4] combine the SARIMA (Seasonal ARIMA) method with LSTM networks to taking into account the linear and non-linear features present in the sea level signal. Their algorithm provides estimates for 24 subsequent months of Sea Surface Anomaly (SSA) using an input dataset consisting of 288 monthly values. Also, ML and DL methods were used for the purpose of integrating a set of oceanographic variables for sea level forecasting at the coastal area of the Western peninsula of Malaysia [31]. They conducted several scenarios including various methods and input variables for predicting monthly mean sea level values. Results showed that LSTM outperforms Support Vector Regression (SVR) and ARIMA when using a combination of oceanographic and atmospheric input variables. Monthly sea level modelling also applied in HIDRA 1 system [12]. It realizes a DL architecture combining LSTM and ConvNN networks for sea level prediction. It’s successor, HIDRA 2 [32], features a higher temporal resolution, using also hourly values the sea surface height (SSH) at Kopper tide gauge (Slovenia), wind speed and MSL pressure aiming at SSH predictions three days ahead.

Ishida et al. [33] proposed an LSTM network to estimate future sea level based on hourly sea level recordings considering several phenomena effecting coastal sea level such as gravitational attraction of moon and sun, seasonality and storm surges. The methodology implemented in the tide gauge of Osaka, Japan and resulted in a sea level model using historical data, while no future prediction was held. Hourly sea level values were also used by Ramadhan et al. [1] to compare the performance of an RNN and an LSTM network for sea level prediction at a span 3–14 days against tidal harmonic analyses. Three-month long datasets were used, originating from a tidal buoy near Sebesi (Indonesia). The LSTM network achieved the best accuracy. Short-term forecasting implementing DL techniques, such as the encoder-decoder LSTM (sequence to sequence) was involved by Vicens et al. [34] to predict 4 days of hourly sea level variation.

Overly, the results found in published literature suggest that high temporal resolution in tide gauge data and meteorological information are important for data processing and reliable sea level forecasting. This is of importance when sea level forecasting aims in tidal cycle characterization or its variation considering a series of cycles, which couldn’t be achieved using computers of lower performance of the near past. Also, previous studies indicate that sea level forecasting depends heavily on the number of time-steps required to achieve an adequate performance, which seems to be more important than resolution. A precise forecast of many timesteps is always welcome for cases of warning or scheduling. The combination of the above two features, that is, the number of the predicted time-steps and the length of each one indicate the overall performance of the methodology. The first feature indicates the forecasting ability, while the second one is mostly determined by the data in use.

One of the longest forecasts found in the literature accounts 288 future time-steps [4] featuring a temporal resolution of one month, which seems to be poor to identify short-term variations. Nevertheless, it is adequate to study a long scale trend for various locations. Local studies demand higher resolution data to identify the several aspects of the phenomenon that contribute to the micro scale of a location. Moreover, an important issue is the availability of high-resolution data. In this study an RNN is trained to forecast hourly values of sea level utilizing few parameters (tide gauge recordings, wind speed and air pressure), reaching to 3 weeks prediction (i.e. over 500 time-steps) implementing LSTM networks, even utilizing data with gaps. Data training took place in a full tidal cycle of 18.6 years of hourly values to precisely model the tidal behavior. Several networks were designed that could be adopted depending on the availability of data.

3.1 LSTM architecture

The general RNNs procedure is mapping an input sequence of observations (X) to corresponding projected values (Y). At each time step, input values pass through the network’s hidden state (h), which then outputs both a projected value and an updated hidden state, based on weight matrices (W) and biases (Figure 3). Backpropagation is then used to compute the weight matrices.

Figure 3:

General depiction of an RNN. The left is the recursive model and the right is the corresponding unfolded RNN based on time sequence.

The hidden state h _t depends on the current observation X _t and on the previous hidden state h _t−1 as described in equation (6).

where f is a nonlinear function. Given the recursive definition of the above equation, the hidden state of an RNN inherently accumulates information from the entire input sequence. This allows RNNs to leverage these hidden variables as a form of a memory to capture long term dependencies. Typically, the hyperbolic tangent (tanh) is used as function f (activation function), and the output is calculated as the product of a weight matrix and the hidden state plus a bias term, as shown in equations (7) and (8), respectively.

where W _hh , W _Xh , W _hY are weight matrices, b _h and b _Y are biases, which are calculated during the training of the network. The tanh function calculates as in equation (9).

To estimate model parameters, the maximum likelihood is typically employed by minimizing an objective function. When computing the derivatives of this objective function with respect to the model parameters, a summation over all time steps i is involved, which includes a multiplicative factor of ∂ h t ∂ h i Traditional RNNs need to backpropagate gradients over long sequences and thus gradient values will shrink after a few time steps. Consequently, states that are far from the current time step contribute minimally to the parameters learning. Conversely, large values of the aforementioned factor can lead to gradient exploding, which has a similar detrimental effect as gradient shrinking (or vanishing gradient) [35]. Detailed calculations of the objective function’s partial derivatives are beyond the scope of this study; interested readers can refer to the relevant literature.

This inherent drawback of the common RNNs, led to the development of LSTM networks. They are specifically designed to address the vanishing gradient problem and excel at remembering the state of the neural networks between predictions and over extended periods. Additionally, LSTM NN support input data with varying sequence lengths, which is a practical advantage. Beyond the hidden state, in LSTMs the cell state is introduced. In an LSTM NN, the hidden state at time step t contains the output of the LSTM layer of that specific time step. The cell state holds information learned from all previous time steps, acting as a long-term memory. At each time step, the layer selectively adds information from this cell state, allowing it to maintain relevant context over very long sequences. The general concept of LSTM NN is illustrated in Figure 4.

Figure 4:

LSTM procedure. At each time step t, the corresponding time step of the sequence (X _t ) and the current states of the network (c _t−1, h _t−1) are used to compute the output (Y _t ) and the updated cell state (c _t ).

LSTM operations are held in cells which are comprised of memory units. The processing is controlled using specialized components called “gates”. These include:

1. the input gate (i), which controls the level of updating the cell state with new information.

2. The forget gate (f), which is a control level that resets the cell state or determine which information is discarded.

3. The cell candidate (g) that adds information to the cell state.

4. The output gate (o), that controls the cell state added to the hidden state.

At each time step, input information passes through the gates, leading to update both the cell state and the hidden state, ultimately yielding the corresponding projected output value. The schematic diagram of an LSTM unit is presented in Figure 5.

Figure 5:

LSTM unit structure including forget gate, input gate, output gate and cell candidate.

The learnable parameters of an LSTM layer are the input weights (W), the recurrent weights (R) and the bias b. Each component or gate is a function of weights, bias, current hidden state and the input sequence. The gate activation functions are the sigmoid function (σ) and the hyperbolic tangent (tanh), as in the equations (10)–(13).

where W _i , W _f , W _g , W _o , R _i , R _f , R _g , R _o are the weight matrices and b _i , b _f , b _g , b ₀ are the biases of the input gate, forget gate, cell candidate and output gate in respect. The tanh function is as in equation (9) and the sigmoid function has as in the equation (14) below.

Considering the gate function the cell state (c _t ) and the hidden state (h _t ) at time step t calculates as in the equations (15) and (16) in respect [35].

where the symbol ⊙ stands for the Hadamard product meaning an element-wise multiplication of matrices.

Finally, the network’s output, Y _t is computed by an output function specifically based on the problem type. For classification tasks, a softmax function is typically used to produce probability distributions over classes. Conversely, for regression tasks, the output is commonly derived from a linear combination of weights and the hidden state plus a bias term. In a typical LSTM architecture, this final output stage is implemented with a fully connected layer followed by either a softmax layer (for classification) or a regression layer (for continuous predictions). The general form of the output function is presented in equation (17).

It is found in literature that input gate includes both functions in equations (10) and (12) or the name of the g _t function as input modulation gate [35], 36]. This study is implemented fully in Matlab^®, so we keep the software notation and semantics.

Let the true recorded values at time t is y _t , then the model parameters can be estimated by minimizing least square of 1 2 y t − Y t 2 . The derivatives with respect to Y _t and W _h are as in equations (18) and (19) respectively.

The gradient with respect to h _T for the last time step T has as in equation (20).

During the backpropagation process, at each time step t, the gradients for the gate functions and the cell state within the LSTM unit are computed as in equation (21).

Subsequently, by back-propagating the activation functions over the whole sequence, the gradients of weight matrices are derived according to equations (22) and (23) [35].

and the corresponding hidden state at current time step t − 1, is as in equation (24), considering the activation function.

Another source of calculating the current hidden state is after equation (24), using the objective function in equation (20).

The above equations are implemented in the back propagation through time (BPTT) technique, used in training LSTM NN on sequential data. BPTT enables the network to learn by iteratively adjusting its weights based on the discrepancy between its predictions and the actual recorded values. Crucially, unlike normal backpropagation, it explicitly accounts for the temporal dimension of sequential data.

Equations (10)–(16) can be used to a forward update of the states as a feedforward NN for time steps 1 to T. Subsequently, the error can be backpropagated from T to 1 using equations (18)–(25). Let θ = {W _i , W _f , W _g , W _o , R _i , R _f , R _g , R _o }. After these gradients are obtained using backpropagation, the model θ can be optimized using gradient-based methods. Common choices include the SGD (stochastic gradient descent) or the Adam (Adaptive moment estimation) stochastic optimizer. In this study the Adam optimizer has been selected. This method is distinguished by its ability to compute individual adaptive learning rates for different parameters, leveraging estimates of first and second moments of the gradients. The Adam optimizer updates network parameters, particularly the weight matrices, as in equation (25) [37].

where m, u are the first (mean) and second (variance) moments of the gradient, a is the learning rate and ε a denominator offset value. The estimates of the moments are calculated after equations (27) and (28).

where ∇ _θ f _t (θ) are the gradients, that is the vector of partial derivatives with respect to θ at time step t.

In practice Adam optimizer is used for minimizing the expected value of the difference between the output of the LSTM network and the recorded true values with respect to parameters θ. The algorithm updates the exponential moving averages of the gradient m _t and the squared gradient u _t . It also includes a bias correction term in the beginning of the training.

The learning rate a and the decay rates of the moments estimates β ₁ , β ₂ are hyperparameters of the LSTM NN. Default values tested for machine learning applications are a = 0.001, β ₁ = 0.9, β ₂ = 0.999 and ε = 10⁻⁸ [37].

The size of the hidden state in an LSTM cell, the hidden size (or hidden units) is a hyperparameter in an analogous concept to neurons or perceptrons in other kind of NN. Typically, it represents the amount of information that the network “remembers” between time steps. Consequently, hidden size defines the size of weight matrices and their dimensions depend only on hidden size and the number of input variables. So, if hidden size is h_s and the number of input variables is i_v then the size of weight matrices W is h_s*i_v and for R is h_s*h_s. In practice, when dealing with a large amount of data such as long time series, input data are processed in batches rather than in one sample at a time. In that way, another hyperparameter is introduced, the batch size. If batch size is b_s, then the dimensions of input is i_v*b_s, the size of hidden state is h_s*b_s and dimensions of the output of each gate is h_s*b_s.

Another hyperparameter is the depth of the network or how many LSTM layers wills be used. The number of LSTM layers implies the number of hidden layers of an LSTM NN or from how many LSTM cells the information of a sequence will pass. The depth is defined by the hidden layers which are not accessible from outside of the neural network, so input and output layers do not count for the depth [38]. The schematic diagram of a deep LSTM NN is depicted in Figure 6.

Figure 6:

Layers of a deep LSTM NN for regression problems.

All time steps get put through the first LSTM layer/cell to generate a whole set of hidden states (one per time step). These hidden states are then used as inputs for the second LSTM layer/cell to generate another set of hidden states, and so on. Afterwards, the output is computed into a “Fully connected layer” and depending on the problem (classification or regression) a “softmax” or a “regression layer” is used to compute the half-squared error.

Another popular category of RNNs is the Gated Recurrent Unit (GRU) which seems effective in predicting short-term sequences [39]. GRUs are a simpler alternative to LSTMs reducing the complexity of the model maintaining similar performance especially when short-term temporal dependencies are more prominent [40]. On the other hand, more recent models such as transformers gain more attention in Deep Learning techniques. Generative Pre-trained Transformers (GPT) initially undergoes unsupervised pre-training typically followed supervised fine-tuning. Its evolution model with enhancements like GPT-2, GPT-3, GPT 3.5, GPT 4 and GPT-5 focus on adhering to user input directives and provide growing efficiency [41]. Long time series like sea level are characterized by intense temporal interdependence and epoch changing alterations and thus, they need the strong covariance of various variables to be captured from their long-term variations. A more recent approach is the sequential-based transformers (S-Transformers) used for predictions of sequential data, based on historical events [42].

Based on the preceding analysis, the LSTM type of RNNs adopted as the most appropriate approach for timeseries forecasting sea level variations. In this regard, a list of scripts developed in Matlab ^® programming software for designing and implementing a Deep Learning architecture. Prior to data analysis, raw data recordings were pre-processed and undergone through additional testing to conclude on the most suited LSTM network formulation.

3.2 Data presentation

Sea level data originate from hourly observations of Thessaloniki (Greece) tide gauge station (Lat = 40^° 37′ 57″.15, Long = 22^° 56′ 05″.76) and 3-h observables from the meteorological station located in the neighboring airport (Lat = 40^° 31′ 38″.64, Long = 22^° 58′ 17″.40). It consists of wind speed, wind direction and atmospheric pressure. The raw observables were processed and analyzed using a tidal model approach for sea level components at Thermaikos Gulf, Greece according to Papadopoulos and Gikas [23]. Processing of these data have resulted to hourly time series of astronomical tide, storm surges, sea level and non-tidal residuals, stress of wind-sea surface for the period 1999–2019. Exported data compiled in a Matlab ^® matfile for further processing in the context of time series forecasting. Specifically, the data comprising the matfile used in this study are:

1. Astronomical tide timeseries (eq. (2))

2. Storm surges timeseries (eq. (3))

3. Sea level observation timeseries

4. Sea level residuals (eq. (1))

5. Constant Mean Sea Level from harmonic analysis for astronomical tide (eq. (2))

6. Horizontal stress components of wind – sea surface timeseries in x, y directions (eq. (3))

7. Wind speed timeseries

8. Air pressure timeseries

9. Wind direction timeseries

10. Coastline azimuth (constant value)

Figure 7 depicts a subset of the available data referred to sea level and its components of Thessaloniki tide gauge station pointed out in Eq. (1).

Figure 7:

Sea level raw observables, astronomical tide, storm surges and sea level residuals timeseries at Thessaloniki tide gauge station (1999–2019).

3.3 Research methodology

As has already been pointed out, this study relies on sea level data analyses using sea level component variables, such as astronomical tide and storm surges as well as other features/physical quantities produced using an analytical approach, such as sea level residuals and non-tidal residuals. Measured meteorological parameters (i.e., wind speed, wind direction and atmospheric pressure) and indirect ones (i.e., stress of wind-sea surface, wind components) support the analysis. In return, the goal is the sea level forecasting for several days contributing to various purposes and applications. Regarding processing methodology, three formulations of LSTM networks were built up, specifically: single variable forecasting (univariate), simultaneous multiple variables prediction (multivariate), and single variable forecasting at post-multiple input (multivariate to univariate). Both closed and open loop forecasting (described in Section 2.2) was implemented by developing several scripts from scratch in Matlab ^® .

3.4 Network architecture definition

Firstly, the original annual timeseries of hourly data (1999–2019) was split to shorter ones of 200–400 h length. It is a common strategy for better data manipulation by the network because LSTM NN handle easily and better short sequences [43]. The resultant timeseries were built up to have no gaps consisting of sequential time steps. Consequently, some information was lost. Each variable sequence has corresponding time reference with the others. The above sub-sampling is arbitrary and even overlapping sequences are allowed. In this study the sub-sampled sequences are consecutive and non-overlapping.

Following, they are separated in training and test sequences with a proportion of 90 % and 10 % respectively. This ratio balance used in many studies in past [1], 2], 4], 12] that proved it provides an adequate number of sequences for testing and evaluation. Also, more training data extended through time provide a model adequate for more predictions in the future. This might cause overfitting which is avoided by using dropout layers between LSTM ones and testing various training epochs that are evaluated.

For every training sequence, when the close loop forecasting method is used, a secondary dataset is formed up using the original data, moved by one timestep. This second sequence reflects the projection of the first one at the training stage. In this manner the network is trained to predict the next time step from a given sequence. Contrarily, for the open loop forecasting method, the input sequences are the input timeseries whilst the output sequence is the projected sea level timeseries. To obtain comparable results, both networks require sequence standardization by subtracting their respective means divided by their standard deviation, as in Equation (29). Specifically, standardization assists network performance as data are transformed to a similar scale evident due to the different physical quantities and units. Finally, variable values are restored to their actual values according to Equation (30).

where, X _train is the training data, X _in is the input data, MV is the mean value of input data, σ _Xin is the standard deviation of training data, X _out is the final output forecast and X _pre is the forecast using the standardized data.

Key statistics (i.e., mean value and standard deviation) of the standardized (Z-score) sequences are produced from the training data, whereas the same parameters are used for restoring the predicted values to their normal scale. These statistics are part of the training stage of the network, so the same parameters should be used to restore the forecasted values at the validation stage and, also, in novel datasets that the model will be used in the future. Test datasets might have different values of mean and standard deviation. Nevertheless, it is a common practice to use the values computed during the training stage to check and validate the trained network as a part of the total evaluation.

Next, the network architecture specifies the layers and parameterization options and the corresponding hyperparameters. The layers adopted consist of an input layer declaring the number of input variables, the LSTM layer(s) containing the hidden units and the gates’ operations, a dropout layer (included when network depth is 2 or higher) randomly removing the 20 % of network nodes, the output layer declaring the number of predicted variables and a regression layer. The selected algorithm for weight matrix and bias estimation is the ADAM (ADAptive Moment estimation), as the stochastic solver. Then, the network is trained, and a first check is performed by predicting the test sequences and comparing them to their original values. Following, the corresponding Root Mean Square Error (RMSE) is computed for this first network evaluation. Low RMSE indicates that architecture and training work well.

Next step refers to future forecasting. Test data sequences are used as input and prediction takes place for multiple timesteps (hours). The forecast horizon is 1, 2 and 3 weeks ahead (i.e., 168, 336 and 504 timesteps-hours) to check the reliability of every period prediction. Univariate network forecast took place for 4 weeks (672 timesteps) instead of 3. Finally, evaluation is held using metrics for the goodness of fit and the prediction ability, which will be analyzed in Section 3.6.

3.5 Optimization (tuning)

The formulation of “training – prediction” timeseries using LSTM networks presupposes the selection of certain training options and hyperparameters. Different training options combinations lead to different results. In fact, the algorithm adheres a stochastic approach so that each re-run might result in a different forecast, even of the case of the same hyperparameters [43]. Hence, the choice of hyperparameters requires extensive experience, both at a theoretical and empirical level, as it resides on the individual characteristics of the physical problem itself. The overall process asks for optimization, in a way that processing time and computational resources are well balanced.

Tuning is the optimization procedure of testing and selecting the appropriate combination of hyperparameters from a reasonable range. In the present study several scenarios of training – forecast network implementation were investigated, concerning the following hyperparameters:

1. Hidden units. Hidden units define the size of weight matrices and bias of the network indicating that carry the amount of information which the network “remembers” between subsequent epochs of training. Thereby, the tested values are 24, 168, 336, 672 and 1,344, which are considered as the information of time-steps corresponding to 1, 7, 14, 28 and 56 days of data respectively.

2. Depth. The network depth is indicated by the number of LSTM layers used during training. It was chosen to test 1, 2 and 3 layers of depth. For depths 2 and 3 after each LSTM layer a dropout layer was used reducing by 20 % the hidden units in the layer to avoid overfitting.

3. Epochs. To define the appropriate number of training epochs while avoiding overfitting the number of epochs tested are 100, 200, 500, 600, 1,000 and 2,000. This range is quite wide and helps selecting empirically the appropriate value.

The tuning strategy consists running a single loop for each network and for all selected combinations of hyperparameters. Those combinations achieved an accuracy higher of the average, pass to the next stage where the script runs again. This procedure is repeated until the best combination is identified. Finally, the best combinations have been run at least four times. The metric of accuracy accepted is the RMSE. This testing procedure reflects the main idea in Bayesian Optimization, whereas hyperparameter combinations from a user-defined range of values are tested [44].

The selected tuning procedure is computationally expensive and is faster when employing the GPU of a PC. Each training process lasted from few minutes to few hours, depending on network depth and training epochs on a standard PC (i5, 11th generation, 16 GB RAM) without engaging the GPU.

3.6 Evaluation

The key metric for the evaluation and comparison among different combinations of hyperparameter selection is the RMSE, as in Equation (31). Furthermore, is used for the overall evaluation of LSTM networks formed up by different input data by the forecasted sea level values and the corresponding tide gauge observations. RMSE computation is taken for test sequences accompanied by the corresponding forecast of each sequence after they are restored to their real-scale values using Equation (30). Other evaluation metrics that were calculated are mean absolute error, mean error and its standard deviation.

where, n is the number of timesteps, y _t is the observed values of sea level and y ̂ t is the value of sea level from forecast after restored to the scale of observations.

In addition, other fitting metrics are considered, which estimate a percentage of fit and evaluate the forecast in relation to the initially observed values of sea level. Known indices are the Willmott Skill (WS) and the Nash-Sutcliffe Efficiency (NSE), defined in Equations (32) and (33) respectively.

where, y _t is the observed value of sea level, y ̂ t is the calculated value of sea level after the forecast and O ̄ and M ̄ their mean values respectively. These indices are widely used for fitting evaluation of a model to observations, especially when hydrographic or oceanographic data are involved [16], 45], 46]. These metrics express the goodness of fit of the predictions to historical data and, usually, they are used to evaluate a model through its residuals.

In cases of forecasting, in addition to model adjustment evaluation, it is a good practice to estimate the predicted accuracy through implementing techniques, known as “out of sample”. These distinguish from others that minimize the residuals. Usually, one–step–ahead error (O–S–A) is adopted according to Equation (34) [47].

where, y ̂ t ( t − 1 ) is the predicted value of y _t which was calculated from previous timesteps.

Typical metrics of forecast accuracy are the mean error (ME), the mean absolute deviation (MAD) and the mean square error (MSE) given by Equations (35)–(37).

The mean error of the forecast is the estimate of the expected error of prediction, indicating the bias in predictions. Mean absolute deviation and mean square error show the variation of prediction error. Index MSE is the direct estimation of error standard deviation of the next timestep.

The above metrics are not always comparable among them, especially when different models or data are used. Also, models and results originated from different time periods encounter the same issue, so metrics that do not depend on the scale and units of data, are needed. To face this concern, the relative error of prediction (Equation (38)), the mean percentage error (Equation (39)) and mean percentage of absolute error (Equation (40)) are introduced as follows [47].

These metrics account for computing relative type or unitless (percentage) errors that represent the deviation from the real values.

The test sequences which are used for validation, serve as the input in the trained network and their future values are forecasted. Each test sequence followed by its forecast is evaluated comparing it to the observed data. There is a point for not evaluating just the forecast, because input and forecasted data values are restored using Equation (30) and its parameters were calculated by training data. In that way, besides the prediction itself, the whole network architecture and robustness is evaluated, too. In the present study total validation sequences cover a time period for more than two years.

Overall, several scripts developed in Matlab ^® to perform the designed model for univariate, multivariate and multivariate to univariate forecast. They differ in various details, mostly on building the datasets and on the output features. Figure 8 depicts the basic model adopted.

Figure 8:

Basic LSTM model structure for timeseries forecasting.

4.1 Univariate forecasting

The univariate forecasting module and the script developed reconstructs two formulations of LSTM networks. The first network accepts as input data the hourly sea level recordings, whilst the second one the hourly non-tidal residuals. In the later, the tidal values are restored before the evaluation of the network. Both networks implement an open loop forecasting (static) model considering that from every input sequence the trained network predicts multiple future timesteps.

For the case the algorithm accepts as input sea level recordings, extended trials among all combinations of hyperparameter values as defined in Section 3.5 revealed the best hyperparameter setup features: 168 hidden units, 2,000 training epochs and network depth of 3 layers. Initial testing of the RNN by forecasting the test sequences themselves has shown that data forecasts fit data recordings with an accuracy of about 4 cm (RMSE) when the maximum annual range is 0.833 m (i.e. less than 5 %). The future forecast lengths tested for the univariate case are 7, 14 and 28 days. Next step is the fitting check of the forecasts. This evaluation considered the input test series accompanied by their forecasted values as one data sequence (23,282, 29,575 and 41,154 total timesteps referred a prediction of 7, 14 and 28 days respectively) compared to observed values. Note, at this stage the test series are restored by the Z-score statistics of the training ones, so in that manner the early standardization is evaluated also. Satisfactory results obtained for the 7 (168 h) and 14-days (336 h) using the NSE index, which for the case of 28 days is lower than 50 % (Table 3). Value of 50 % of the NSE metric means that the forecast could also be represented by the data average value. The resulting values of the O–S–A prediction indices lie in low levels. The prediction seems unbiased since the mean error is, also, low in all cases.

The best hyperparameter combination for the case that training relies on non-tidal residuals results in 24 hidden units for data forecasting, 600 training epochs and network depth of 3 LSTM layers. The first RNN control by predicting the test data themselves have shown a good fit with an RMSE value below 3.5 cm for a time-period spanning for several years. Sea level forecasting refers to 45 sequences (using 404 training sequences), with length 200–400 timesteps (hours). Results indicate a good performance for the 7- and 14-days forecasting, contrary to the 28-day forecast that barely produces satisfactory results. These findings are supported by the high values (nearly 100 %) obtained for the WS and NSE indices for the first two cases, and lower (70 %) NSE values computed for the third case. Figure 9 depicts the RMS values obtained for the forecasting part of test sequences. As far as the O–S–A evaluation is concerned, forecasting seems to be unbiased considering bias size (mean error) is very close to zero.

Figure 9:

RMS values of the forecasting minus observed series per sequence for 7 (A), 14 (B) and 28 (C) days of prediction (univariate forecast implementing non-tidal residuals).

Considering the long forecasting period ranging from one to four weeks (i.e., 168 to 672 timesteps ahead), these results characterized as sufficient to exceptional. The best network’s performance for LSTM univariate modelling is shown in Table 3.

Furthermore, the evaluation results shown in Table 3 suggest that forecasting is more accurate when the astronomical tide is subtracted from the original sea level data. Evidently, it appears that seasonality and trend (see Eq. (2)) decrease the network’s ability to produce a model suitably adjusted to the original observations. In contrast, subtraction of seasonality and trend exports a model that demonstrates better adjustment to the observation data. Nevertheless, evaluation is retained in sea level values as astronomical tide is restored after non-tidal residuals forecasting. Figure 10 depicts a standard forecasting case of good quality.

Figure 10:

Comparison between forecasted and observed data for a test sequence of 7 (A), 14 (B) and 28 (C) days using a univariate forecast implementing non-tidal residuals.

The prediction ability of the network as evaluated by the One–Step–Ahead error index values, indicates an unbiased produced model from non-tidal residuals [47]. This ability is enhanced further from the low variance of the respective errors. Error percentage is close to zero in both cases and their mean absolute value favors the model produced by the non-tidal residuals. Figure 11 depicts the design architecture of the LSTM network employing non-tidal residuals.

Figure 11:

Univariate LSTM network architecture employing non-tidal residuals as input.

4.2 Multivariate forecasting

Considering a multivariate network setup using sea level input recordings, four network types were studied via detailed testing of several hyperparameter combinations. For these network formulations we implemented the static (open loop forecasting) setup thanks to its multiple timestep prediction capability.

The first point to note from the analysis is the satisfactory results obtained for the timespans of 7 and 14 days, albeit the 21 days forecast underperformed (NSE index less than 50 %). Interestingly, the results obtained for dataset D (Table 1) are less accurate even for the 14 days prediction length.

Considering the performance of the prediction model, the most promising forecasting results were obtained using dataset A (i.e., astronomical tide, storm surges, sea level residuals, sea level recordings) for 7 and 14 days of forecast (Table 4). This is due to the higher values of WS and NSE metrics and the low value of the respective RMSE. Training took place in 390 sequences representing 200–400 h of data length, whilst evaluation assumes the 44 test sequences supplemented by the corresponding forecasting period.

Considering the O–S–A accuracy these results are satisfactory (Table 4) for all data combinations exhibiting a low bias by their Mean Error, except for dataset D. Notably, using different performance metrics (MAD and RMSE) result in similar values for all categories, while at the same time, error percentages (MPE and MAPE) exhibit low values.

Also, LSTM network behavior concerned with non-tidal residuals was studied (see Section 3.3.2. for data description). Analysis reveals significant improvement for all input data combinations when seasonality and trend do not affect the LSTM network by removing astronomical tide before training. Dataset E1 (i.e., wind velocity across the coastline, atmospheric pressure and non-tidal residuals) provides the best fit for the 7-, 14- and 21-days prediction time (see Table 4, last 3 columns). The corresponding timeseries are shown in Figure 12. Prediction evaluation by O–S–A indices, is almost clear of biases exhibiting a very low variation. Error percentages are also very small and improved in relation to case of dataset A.

Figure 12:

Input variables from dataset E1.

Finally, as a general remark, analysis suggests that training for the data categories including non-tidal residuals require a smaller number of training epochs compared to other cases. This might be attributed to overfitting implications as the number of training epochs increases leading to decreased prediction accuracy. The overfitting can also be induced by smoothing out the residual signal as the training epochs increase.

Dataset E1, achieves the best forecast with lesser hidden units than dataset A employing only one layer depth. Figure 13 plots the forecasted values for a test sequence of average performance. The diagram of the multivariate LSTM network exhibiting best results is as in Figure 14.

Figure 13:

Comparison of forecasted values to observation data of a test sequence for 7 (A), 14 (B) and 21 (C) days (multivariate forecast, dataset E1).

Figure 14:

Multivariate LSTM network architecture. Sea level is computed by non-tidal residuals after restoring astronomical tide.

4.3 Multivariate to univariate forecasting

Networks featuring multiple input/single output variables constitute a special type of LSTM networks in this study, due to the differences underlying in the training methodology and the forecasting procedure. Seven data categories (A2, B2, C2, D2, E2, F2 and G2) were used to study the performance of these networks. Contrary to the previously studied LSTM network types, in this case, the network returns only the sea level estimate via open loop forecasting. In this setup, once all input data timeseries have entered the network, a prediction output is computed just for the next time step provided all input variables of this time step are included in the network. However, this model assumes the forecasting time steps should not bear any gaps in the input variables. Likewise, if a gap is encountered in the data, no prediction is performed for this time step, and the forecasting process ends for the remaining of that sequence. Clearly, this is a drawback of the open loop forecasting method, as the algorithm withdraws the current sequence and moves directly to the next one. During the sequence’s formation in this study, datasets have been treated so that no gaps exist in both training and test sequences; however, there is no gap control in forecasting time steps for input test variables.

Output forecasting using datasets C2 and F2 failed entirely to produce reliable predictions due to their excessive RMSE values and the negative values of fit indices (NSE). Most of the datasets lead into incomplete forecasting as data gaps in the forecasting timesteps in input series have stopped the prediction (A2, B2 and E2), whereas only datasets D2 and G2 result in complete forecasting for all future time steps (7, 14 and 21 days) requested from the algorithm.

Interestingly, despite the fact that the forecast is incomplete, the resolved part of all datasets exhibits consistently accurate results; in certain cases, even characterized by exceptionally good statistical values; specifically, up to 99 % goodness of fit using WS and NSE indices. Their accuracy figures are similar among sequences, whereas O–S–A statistics exhibit a significant dispersion considering ten times spread among the largest and lower values. Dataset A2 provides the best results for its low RMSE value and the very low error percentages in O–S–A evaluation (see Table 5). The number in parentheses next to the forecasting time refers to the mean values of days considering all the test sequences that prediction finally was achieved.

Having a closer look at the field data, the problem of incomplete forecasting sought to spotted gaps in meteorological data resulting in lack of distinct surge data.

To overcome this problem, we attempt a combined approach.

Firstly, we propose a multivariate forecasting analysis (multiple outputs), so that in addition to sea level predictions (or non-tidal residuals), all other input variables are also predicted. Investigation of the multivariate network that exhibits the best forecast for storm surges (i.e., lower RMSE) evince the network/dataset which could fill in the gaps best. A dedicated script was developed to locate the gaps in the input data (i.e. storm surges timeseries) and, subsequently, to fill them in by the corresponding predicted values of the multivariate forecast. Following, the multivariate to univariate part of the script runs again. The respective results (fill in and predict technique) are presented in last three columns of Table 5 and the training progress (Matlab ^® output) in Figure 15.

Figure 15:

Training progress of multivariate to univariate forecast utilizing dataset A2.

The results obtained from the “fill in and predict” algorithm are very promising. As shown in Table 5, the predicted values match very closely the recorded sea level values. Furthermore, this observation reflects in the very high (over 99 %) values of WS and NSE parameters and the low (of the order of 1 cm; last column of Table 6) RMSE values in the forecasted values. The fit of test sequences supplemented by the respective forecast compute around 4 cm due to the restored standardization parameters calculated from the training data. The good performance of this model resumes to the fact that the predicted section of test sequences exhibits better fit to the recorded data than those from the de-standardized test sequence itself. In all cases of univariate and multivariate forecasting, the predicted part of the test sequences exhibits worse evaluation index values than the test sequence supplemented by its prediction. Also, the mean error close to zero indicates unbiased fit. Figure 16 presents the RMS values of the forecasted values for the 7, 14 and 21 days of prediction. The forecast to observation comparison for a sequence exhibiting some of the worst results are shown in Figure 17 concerned with the 7, 14 and 21-days of prediction.

Figure 16:

RMS of the forecast minus observations per sequence for 7 (A), 14 (B) and 21 (C) days of prediction (multivariate to univariate forecast, fill in and predict, dataset A2).

Figure 17:

Comparison of forecasted values to observation data of a test sequence for 7 (A), 14 (B) and 21 (C) days (multivariate to univariate forecast, fill in and predict, dataset A2).

Besides the success of fitness, evaluation of prediction exhibits good results by the O–S–A evaluation. Mean errors are below 2 mm suggesting an unbiased prediction. Error percentages exhibit also very low values. Figure 18 depicts the multivariate to univariate LSTM network architecture leading the best performance.

Figure 18:

Multivariate to univariate LSTM network architecture.

4.4 Summary

Table 6 summarizes the results obtained from the analyses of all types of LSTM networks and datasets that achieved the best performance for the 7 days forecast. Also, it depicts the goodness of fit values for test sequences supplemented by the next time steps forecast and solely the forecasted time steps, separately. Likewise, the datasets resulted in the best univariate (univ), multivariate (multiv) and multivariate to univariate (multi2uni) forecasts are also presented. For the case of univariate forecasting the datasets consist only sea level (SL) and non-tidal residuals (NTR). Multivariate networks are represented through the dataset of best performance in full form, also assuming de-seasoning and de-trending (multi-de:se-de:tr). Note that for the case of multivariate to univariate (multi2uni) results, the days of prediction are limited to 3.7 instead of 7, due to data gaps in storm surges, so the comparison is not held in same terms. For that reason, the corresponding results implementing the technique fill in and predict (multi2uni-finp) are quoted.

The evaluation of fitness for the 7 days forecast shows similar results when the seasonality and trend (through astronomical tide) are removed. The difference in the RMSE values between the univariate and multivariate forecasts is below 1 mm, whereas the values obtained for WS and NSE metrics certify the success of fitness to the original data. The low mean error in the two cases of multivariate forecasting indicates that there is no bias matching to historical data and the univariate forecasting significantly falls behind it.

The Goodness of fit in forecasted values of test sequences, for 7 days (with similar results for longer forecasts), results in substantial success for the multivariate to univariate approach. For the case of incomplete forecast (i.e., 3.7 days) fitting is extremely good. The complete forecast of 7-days as also the longer periods, show almost perfect fit to historical data exhibiting very low RMSE value, of the order of 1 cm. The mean error is below 1 mm suggesting an unbiased fit, considering the initial data resolution of 1 cm.

Apart from the incomplete forecast, all cases present similar results in O–S–A evaluation, whilst the mean error is small with low variance. Finally, the error percentages exhibit low values with absolute error in the order of 2 % in all cases.

Additionally, the Pearson Correlation Coefficient (PCC) calculated at high rates (0.89–0.95) for the best scenario of each technique (univariate, multivariate and multivariate to univariate), while p-values were zero. When only the forecasted parts of the validation sequences were used, the PCC dropped to values of 0.74–0.80 for the univariate and multivariate cases and p-values were calculated again zero. For the multivariate to univariate case considering only the forecasted part of the validation sequences, the PCC got even higher reaching to values of 0.98–0.99.

In general, when training and forecast take place in data abridged from seasoning and trending a better fit is obtained exhibiting lower RMSE values, too. Also, dynamic forecasting, predicting each time step of sea level after inserting the corresponding variables to the network present a perfect match compared to the recorded values. These conclusions are shown clearly in Figures 19 and 20.

Figure 19:

Test sequences performance of univariate, multivariate and multivariate to univariate networks forecasting sea level values for various time spans by sea level, A and A2 datasets respectively. The dot location of each case represents the fit by Willmott Skill index and its size represents the RMSE value.

Figure 20:

Test sequences performance of univariate, multivariate and multivariate to univariate networks forecasting sea level values for various time spans by non-tidal residuals, deseason-detrend E1 and technique fill in and predict in A2 datasets respectively. The dot location of each case represents the fit by Willmott Skill index and its size represents the RMSE value.

Table 7 contains the detailed parameter information adopted for the Deep Learning NN resulted in optimal sea level forecasting associated with the network architecture depicted in Figure 18.

Conclusively, in the present study various data variables were used for LSTM NNN formation including deliverables from analytical sea level modelling and not just recordings, binding the two approaches. Also, an interchangeability between techniques (univariate, multivariate, multivariate to univariate) is provided offering a variety of choices depending on data availability. The effect of seasonality and trend removal is studied concluded to better results in many cases. The dynamic prediction implementing open loop forecasting proved to outperform significantly the static one, exhibiting stable accuracy even in longer predictions. Variable selection meditated and hyperparameter combinations, as well, concluding to a model tuning procedure. Besides sea level, other variables were predicted, also, implementing the multivariate technique assisting solving issues with data gaps which lead to incomplete forecasts. The presented procedure was cross-validated by repeating the same forecast using the described strategy for forming LSTM NN at the site for Peiraeus tide gauge [48]. The results are comparable to those of this study demonstrating slightly better accuracy.