Sea Level Time Series Reconstruction of the California Coast Using Neural Network Residual Imputation and ERA5 Data

1. Introduction

Sea level is mainly influenced by the gravitational interaction between the Earth, Moon and Sun, which generates phenomena such as tides [1]. In addition to the astronomical component, sea level can also be modulated by atmospheric and oceanic factors, such as wind, atmospheric pressure, and the regional dynamics of currents, as well as by large scale patterns such as the El Niño Southern Oscillation (ENSO) and the Pacific Decadal Oscillation (PDO) [2]. During the 20th century, mean sea level is estimated to have risen by approximately 10 cm, and it is currently increasing at a rate of about 1–2 mm per year; this change is associated with global warming, thermal expansion, and glacier melt [3]. As a result, the risk of coastal flooding, upwelling, and shoreline erosion increases, representing a threat to the infrastructure and the safety of cities located in coastal regions [4].

Sea level time series represent an effort to collect and analyze chronologically organized information in order to understand the historical behavior of the ocean and estimate its future evolution [5]. However, missing data in these series can induce errors in the analysis, leading to inaccurate trend estimates and hindering the identification of periodic signals such as seasonal cycles or tide related components. This problem reduces the reliability of statistical models and limits the ability to generate robust projections for future scenarios [6].

The California coast has an approximate shoreline length of 5,515 km, [7] making it a region with extensive coastal exposure and pronounced sea level variability. In addition, the California coast hosts several institutions dedicated to scientific research, including the Southwest Fisheries Science Center, CalCOFI Laboratories, the California Ocean Science Trust, and the Ames Research Center. These institutions contribute to the understanding of sea level variability, marine ecosystem changes, and the development of technologies for the conservation and sustainable management of marine environments.

The present study uses sea level observations recorded by tide gauges operated by the National Oceanic and Atmospheric Administration (NOAA), considering NOAA’s international recognition as a leading authority in oceanographic observations, the high quality and consistency of its measurements, the availability of long term historical records, and the open accessibility of its data.

Different approaches can be used, for time series reconstruction, ranging from classical statistical methods to neural network based models. Among the classical methods, one of the most widely used is the autoregressive integrated moving average (ARIMA) model [8]. A more elaborate variant is the Seasonal ARIMA (SARIMA) model, which explicitly incorporates periodic components and allows a better representation of regular seasonal patterns commonly observed in time series data [5]. However, one of its main limitations is that it often assumes linear relationships between the predictor variables and the target variable. To overcome this restriction, neural network models can be implemented, as they allow more complex and nonlinear patterns present in time series dynamics to be captured [9].

Deep learning based models have become particularly relevant for time series modeling and missing data imputation, including architectures such as GRU-D [10], M-RNN [11], and BRITS [12] as well as Convolutional neural network (CNN)-based models, such as TimesNet [13] and GP-VAE [14], as well as graph neural networks (GNNs) [15]. More recently, attention-based models, including CDSA[16], Transformers [17], SAITS [18], and NRTSI [19], have demonstrated superior performance over traditional statistical approaches. Despite these advantages, such improvements often come at the cost of increased computational complexity and higher memory requirements. In contrast, multilayer perceptrons (MLPs), when combined with well designed features, provide a computationally efficient alternative while maintaining model simplicity.

An important distinction in time series modeling is that between forecasting and imputation. While forecasting estimates future values from previously available information, imputation reconstructs missing observations within an existing time series. Recurrent neural networks (RNNs) are well suited to sequential modeling but may face difficulties with long contiguous gaps due to the lack of observations within the missing interval. In contrast, the approach adopted in this study treats sea level reconstruction as an imputation problem, using an MLP to estimate the non tidal residual from meteorological and temporal predictors.

A widely used type of neural network is the multilayer perceptron (MLP), known for its application in pattern recognition. In this model, information flows only forward, from the input layer through one or more hidden layers to the output layer. During training, the model parameters (weights and biases) are adjusted using a backpropagation algorithm in order to minimize the error and produce outputs that become increasingly closer to the expected values [20]. From a theoretical perspective, MLPs are universal function approximators, meaning that they can approximate any continuous nonlinear function with arbitrary accuracy given a sufficient number of hidden units [21]. In practice, however, well designed MLP architectures often achieve competitive performance with a relatively small number of neurons, resulting in reduced computational cost during training.

Sea level can be described as the combination of two main components: the astronomical tide and the meteorological tide. The astronomical tide responds to the gravitational forcing by celestial bodies and is analyzed using the T-TIDE toolbox proposed by Pawlowicz. [22], allowing the extraction of tidal signals. On the other hand, the meteorological tide corresponds to phenomena unrelated to astronomical interactions. In this context, the main contribution of the proposed methodology consists of modeling the residual component through a multilayer perceptron (MLP), designed to predict the non tidal variability and support the reconstruction of missing data in the time series.

This paper is organized as follows: Section 3 presents the materials and methods, including sea level series data source, handling of missing and invalid data, harmonic analysis, residual determination, meteorological variables obtained from ERA5, MLP training, network architecture, training strategy, residual imputation, and model evaluation. Section 4 presents the results of the imputations for the La Jolla, Los Angeles, San Diego, and San Francisco stations, together with the statistical evaluation of the reconstructed sea level series. Section 5 discusses the results obtained. Finally, Section 6 summarizes the main conclusions and future research directions.

2. Workflow

Figure 1 presents the general workflow of the proposed methodology for sea level time series reconstruction. The process begins with the acquisition of sea level observations from NOAA tide gauge records. These data are first subjected to preprocessing, including the handling of missing or invalid values and normalization to ensure consistency for subsequent analysis. Next, harmonic analysis is applied to estimate the astronomical tidal component of the signal. Once the tidal contribution is identified, the residual component is obtained by subtracting the harmonic reconstruction from the original observations.

Afterward, meteorological variables from the ERA5 reanalysis dataset are incorporated as external predictors. These variables, together with temporal information, are used in the training stage of a multilayer perceptron (MLP) model. The trained network is then employed to estimate missing values within the residual series through a residual imputation step. Finally, the complete sea level signal is reconstructed by combining the harmonic component with the completed residual series, and the resulting time series is evaluated using statistical performance metrics.

3. Materials and Methods

3.1. Sea Level Time Series

Sea level data were obtained from tide gauge stations of the [23], while the associated atmospheric variables were extracted from the ERA5 reanalysis dataset developed by the European Centre for Medium Range Weather Forecasts (ECMWF) [24]. The experiment used sea level records from four tide gauge stations in the California Coast: San Diego, La Jolla, San Francisco, and Los Angeles. The analysis period covered January 1, 1950, to June 30, 2025 with a temporal resolution of one hour. Figure 2 shows the geographic location of the four stations considered in the study, while Table 1 summarizes the main characteristics of each dataset.

3.2. Handling of Missing and Invalidad Data

To ensure a consistent temporal structure and to properly identify missing or invalid observations, time series were preprocessed using a gap completion strategy based on an hourly time grid. First, the raw dataset was inspected to identify incomplete timestamp information. The dataset contains the fields year, month, day, hour, minute, and water level. Records with missing values in any of the timestamp fields (year, month, day, hour, or minute) were removed, as they cannot be mapped to a valid datetime index. No further preprocessing steps were required, as the dataset already met the necessary quality standards, exhibiting no redundant entries, duplicate records, or significant outliers that required correction.

After constructing the datetime index, a complete hourly sequence was generated from the earliest to the latest available timestamp. The dataset was then reindexed to this full hourly range, forcing the insertion of any missing time steps as explicit gaps. As a result, all missing observations are represented as NaN values in the water level field, while maintaining a temporally uniform sampling interval. Finally, the datetime index was converted back to separate columns (year, month, day, hour, minute) to preserve the original file structure and facilitate downstream processing.

3.3. Harmonic Analysis

In many oceanic time series, tides represent one of the most important components of sea level variability and exhibit a predominantly periodic behavior associated with astronomical forcing. Consequently, they can be modeled using classical harmonic analysis, which represents the signal as the sum of periodic waves with known frequencies, such as $M_2$ and $S_2$, among others. Under this approach, amplitudes and phases are estimated through a least squares fitting procedure, allowing tidal signals to be reconstructed and predicted even during periods with incomplete observations. Tools such as T-TIDE implement this methodology and are widely used in oceanographic studies [22]. Once the astronomical component has been removed, the residual signal can be modeled using meteorological predictors in order to capture non tidal variability.

Therefore, the proposed methodology separates the observed sea level $\eta \left(t\right)$ into an astronomical component ${\eta }_{\mathrm{tide}}\left(t\right)$ and a residual component ${\eta }_{\mathrm{res}}\left(t\right)$. The astronomical component is estimated through classical harmonic analysis, in which the tide is represented as a finite sum of cosine functions with known frequencies. In particular, the astronomical tide can be expressed as:

$${\eta }_{\mathrm{tide}}\left(t\right)=a_0+{\displaystyle \sum _{i=1}^N}{a}_{i}cos\left({\omega }_{i}t+{\varphi }_i\right),$$

(1)

where $a_0$ denotes the mean sea level, while $a_i$ and ${\varphi }_i$ represent the amplitude and phase lag of the i-th tidal constituent, respectively. Here, ${\omega }_i$ corresponds to the angular frequency associated with each tidal constituent (e.g., $M_2$, $K_1$, $S_2$, and $O_1$), and N is the total number of constituents considered [25].

Harmonic tidal analysis is a conventional approach in oceanography [22] and is widely used by operational tidal services to reconstruct and predict the astronomical tidal component [26]. In this study, this harmonic only reconstruction is therefore considered as a reference approach, while the proposed methodology extends it by incorporating the meteorologically driven non tidal residual.

3.4. Residual Determination

Once the astronomical component was reconstructed for each temporal window, the residual signal was obtained as the difference between the observed sea level and the synthesized tidal component:

$${\eta }_{\mathrm{res}}\left(t\right)=\eta \left(t\right)-{\widehat{\eta }}_{\mathrm{tide}}\left(t\right).$$

(2)

The residual represents the non astronomical contribution to sea level variability. Because the harmonic reconstruction preserves only strictly periodic components, the residual retains low frequency meteorological signals, instrumental biases, and other unknown contributions that may be either stochastic or deterministic. This separation allows the modeling framework to focus exclusively on the physically meaningful non tidal variability using atmospheric features from ERA5.

3.5. Meteorological Variables (ERA5)

Meteorological features were obtained from the ERA5 reanalysis produced by the Copernicus Climate Change Service (C3S) and accessed through the Climate Data Store (CDS) API. Data were retrieved from the dataset reanalysis era5 single levels, which provides hourly global reanalysis fields at the surface and near surface level.

For each tide gauge location, ERA5 variables were extracted over a rectangular spatial domain defined by latitude–longitude bounds in the CDS “area” format [N, W, S, E]. The bounding box was selected to represent local atmospheric and oceanic forcing conditions near the station.

The following physical variables were downloaded with hourly resolution:

• East-West wind component in meters per second [m/s] at a height of ten meters above the surface of the Earth
• North-South wind component in meters per second [m/s] at a height of ten meters above the surface of the Earth
• Mean sea level atmospheric pressure adjusted to the height of mean sea level in [Pa]
• Sea surface temperature near the surface in [K]
• Sea surface atmospheric pressure at the surface of land, sea and inland water in [Pa]

Data were obtained as monthly NetCDF files using a year month loop. Each request included all calendar days and the 24 hourly time stamps (00–23), allowing CDS to return the valid subset for each month.

3.6. MLP Training

3.6.1. Meteorological Features (ERA5)

In the temporal dimension, ERA5 variables were retrieved at hourly resolution and synchronized with the corresponding sea level observations.

$${\mathbf{x}}_{\mathbf{p}}={\left[\begin{array}{c}u\left(t\right),v\left(t\right),msl\left(t\right),sst\left(t\right),sp\left(t\right)\end{array}\right]}^{\mathsf{T}}$$

(3)

where ${\mathbf{x}}_p$ represents the meteorological forcing vector used as input to the neural network.

3.6.2. Cyclical Temporal Variables (Sine/Cosine)

To avoid artificial discontinuities (for example, between 23:00 and 00:00, or between December and January), time variables are transformed into a cyclical representation using sine/cosine pairs. In particular, the following vector is constructed:

$${\mathbf{x}}_{\mathbf{t}}={\left[\begin{array}{c}sin\left({\theta }_Y\right),cos\left({\theta }_Y\right),sin\left({\theta }_M\right),cos\left({\theta }_M\right),sin\left({\theta }_D\right),cos\left({\theta }_D\right),sin\left({\theta }_H\right),cos\left({\theta }_H\right)\end{array}\right]}^{\mathsf{T}}$$

(4)

where ${\theta }_Y$, ${\theta }_M$, ${\theta }_D$, and ${\theta }_H$ represent the angular transformations associated with the year (Y), month (M), day (D), and hour (H), respectively. Each temporal variable is encoded through its sine and cosine components, $sin\left(\theta \right)$ and $cos\left(\theta \right)$, enabling a continuous representation of cyclical temporal patterns. Finally [27]. The complete input vector is defined as:

$$\mathbf{x}={\left[\begin{array}{c}{\mathbf{x}}_{\mathbf{p}},{\mathbf{x}}_{\mathbf{t}}\end{array}\right]}^{\mathsf{T}}.$$

(5)

3.7. Network Architecture

The complete input vector $\mathbf{x}$ is provided to a fully connected feed forward artificial neural network, known as a multilayer perceptron (MLP) with three hidden layers. The architecture is defined by a hyperparameter vector $\mathbf{k}={[k_1,k_2,k_3]}^{\mathsf{T}}$, where $k_i$ denotes the number of neurons in the i-th hidden layer. In the present implementation, the network configuration was set to $\mathbf{k}={[50,20,10]}^{\mathsf{T}}$. Hidden layers employed the hyperbolic tangent sigmoid activation function (tansig), while the output layer used a linear activation function (purelin). The overall architecture of the network is illustrated in Figure 3.

3.8. Training Strategy

The network was implemented using MATLAB’s Neural Network Toolbox.Physical predictors were normalized to the range $[-1,1]$, while temporal variables were encoded using sine cosine transformations to ensure a smooth cyclical representation. 70% of available residual samples were used for training, 15% for validation (early stopping) and 15% for cross validation, implementing the dividetrain method. The loss function $\mathcal{L}$ was defined as the mean squared error (MSE).

$$\mathcal{L}={\displaystyle \frac{1}{N}}{\displaystyle \sum _{t=1}^N}{\left(Y_t-{\widehat{Y}}_t\right)}^2.$$

(6)

Model optimization was carried out using the scaled conjugate gradient (SCG) algorithm. The training process was configured with a maximum of 1000 epochs, a minimum gradient threshold of ${10}^{-7}$. To prevent overfitting, early stopping was applied: training was halted when the validation loss did not improve for 10 consecutive epochs, and the best-performing model parameters were restored.

For the reconstruction of missing values, the MLP model was applied using sliding temporal windows of five years. These windows were selected based on the presence of gaps in the time series, ensuring that each training segment contains sufficient valid data to learn the underlying patterns. Within each window, the model was trained and subsequently used to predict the missing values corresponding to that specific period.

3.9. Residual Imputation

Once the network parameters were optimized using valid residual observations, the trained model was applied to the complete input dataset in order to estimate missing residual values. Let $Y\left(t\right)$ denote the observed residual component and $\widehat{Y}\left(t\right)$ the MLP prediction used for time steps where $Y\left(t\right)$ is unavailable. The reconstructed residual series, $Y_{\mathrm{r}}\left(t\right)$, is defined as:

$$Y_{\mathrm{r}}\left(t\right)=\left\{\begin{array}{cc}Y\left(t\right),\hfill & \mathrm{if}\mathrm{available},\hfill \\ \widehat{Y}\left(t\right),\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(7)

3.10. Time Series Reconstruction

Once the residual component was completed through the imputation procedure, the full sea level time series was reconstructed by combining the harmonic prediction with the residual signal. Let $H\left(t\right)$ denote the harmonic component obtained from the classical tidal analysis and $Y_r\left(t\right)$ the reconstructed residual series. The reconstructed sea level signal, $S_r\left(t\right)$, is then obtained as

$$S_r\left(t\right)=H\left(t\right)+Y_r\left(t\right).$$

When observed residual values are available, the reconstruction uses the original residual component. Otherwise, the missing residuals are replaced by the MLP based imputation.

3.11. Model Evaluation

Model performance was evaluated in two phases, first predicting already known data and evaluating the coefficient of determination ($R^2$) and the root mean squared error (RMSE). The coefficient of determination is defined as:

$$R^2=1-{\displaystyle \frac{{\sum }_{t=1}^N{\left(Y\left(t\right)-\widehat{Y}\left(t\right)\right)}^2}{{\sum }_{t=1}^N{\left(Y\left(t\right)-\overline{Y}\right)}^2}},$$

(8)

where $Y\left(t\right)$ denotes the observed residual, $\widehat{Y}\left(t\right)$ the predicted residual, and $\overline{Y}$ the mean of the observed residual series. The $R^2$ coefficient quantifies the proportion of residual variance explained by the network model. The RMSE was computed as:

$$\mathrm{RMSE}=\sqrt{{\displaystyle \frac{1}{N}}\sum _{t=1}^N{\left(Y\left(t\right)-\widehat{Y}\left(t\right)\right)}^2},$$

(9)

The RMSE measures the average magnitude of the deviation between predicted and observed residuals.

Additionally, to assess whether the imputation preserves the fundamental statistical properties of the time series, descriptive statistics were compared between the observed and imputed data: The mean of the original residual series (${\mu }_{\mathrm{orig}}$) was compared with the mean of the imputed residual series (${\mu }_{\mathrm{imp}}$) in order to evaluate the preservation of central tendency. Similarly, the variance of the original series (${\sigma }_{\mathrm{orig}}^2$) was contrasted with that of the imputed series (${\sigma }_{\mathrm{imp}}^2$) to assess consistency in variability.

Beyond central tendency and dispersion measures, statistics including skewness and kurtosis were analyzed following the same approach. The skewness of the original residual series (${\gamma }_{\mathrm{orig}}$) was compared with that of the imputed series (${\gamma }_{\mathrm{imp}}$) to determine whether distributional asymmetry was preserved. Also, the kurtosis of the original series (${\kappa }_{\mathrm{orig}}$) was evaluated against the kurtosis of the imputed series (${\kappa }_{\mathrm{imp}}$) to verify the maintenance of tail behavior after reconstruction.

4. Results

4.1. Residual Imputation

A set of residual imputation results examples, for the selected stations, are shown in Figure 4, Figure 5, Figure 6 and Figure 7, where missing intervals are replaced using the corresponding MLP based estimates. Missing intervals are replaced using the corresponding MLP based estimates.

A summary of the statistical metrics comparing the observed residual series and the MLP based imputed values for the La Jolla, Los Angeles, San Diego and San Francisco stations are presented in Table 2, Table 3, Table 4 and Table 5.

The performance of the proposed MLP based imputed residuals was assessed using residual sea level records from the La Jolla, Los Angeles, San Diego, and San Francisco stations. A comparison between the observed and reconstructed residual series revealed that the reconstructed data consistently reproduced the overall statistical behavior of the original records, including their mean, variance, skewness, and kurtosis. The mean values of the reconstructed residuals remained close to zero and were consistent with those of the observed residuals, as expected from the centering of the residual series during preprocessing, indicating preservation of the central tendency. Likewise, the reconstructed residuals retained variance values of the same order of magnitude as the observed series, demonstrating preservation of the main dispersion characteristics of the residual data. With respect to reconstruction accuracy, all study cases yielded coefficients of determination $R^2$ greater than 0.60, with the highest values exceeding 0.85, while RMSE values remained below 0.055 m across all experiments.

4.2. Statistical Evaluation of the Reconstructed Sea Level Series

A statistical comparison between the observed sea level series and the reconstructed sea level series was performed for the four tide gauge stations. Table 6 presents the corresponding statistics for the classical harmonic analysis, while Table 7 summarizes the main statistical metrics used to evaluate the reconstruction performance, and Table 8 presents the main statistical metrics used to evaluate the reconstruction performance. For each station, the reported $R^2$ corresponds to the mean coefficient of determination obtained from the different reconstruction windows evaluated. Overall, the inclusion of the residual model reduced the RMSE by approximately 52% across the four stations.

To better highlight the improvement obtained by the MLP model with harmonics imputed residual values relative to the classical harmonic values, Table 8 presents the differences in performance metrics.

The MLP with residual reconstruction shows a significant improvement over the classical harmonic approach. As shown in Tables 6–8, the proposed method increases $R^2$ from an average of 0.9794 to 0.9955 across the four stations, while reducing the RMSE from an average of 0.0743 to 0.0348. This represents an average RMSE reduction of 53%, indicating that the MLP with residual imputation substantially improves the reconstruction accuracy of the sea level series despite the limitations imposed by the fixed training window.

5. Discussion

The combination of the astronomical tide with the imputed residual with meteorological variables produces an highly accurate reconstruction of the observed sea level series. The coefficients of determination $R^2$ for all the complete series between 1 January 1950 to 30 June 2025 are close to 0.995 at all stations (Table 7), quantifying the gain obtained by incorporating the residual modeling.

The proposed model obtains a clear improvement in predictive performance in comparison with other approaches reported in the literature, such as the feedforward model proposed by [28], which reports $R^2$ values from 0.8318 to 0.8775 and RMSE values between 0.059 and 0.079. Similar results are observed when comparing the proposed approach with several models reported by [29], including Prophet, LSTM, and GRU. These models report performance values of approximately $R^2$ 0.71 to 0.95 and RMSE values around 0.054 to 0.023.

These comparisons should consider the distinction between forecasting and imputation. Although recurrent models such as LSTM and GRU effectively capture temporal dependencies, their performance may decrease when imputing long contiguous gaps due to error propagation from sequential estimates. In contrast, the proposed MLP estimates the non tidal residual from meteorological and temporal predictors without relying on previous reconstructed values, making it suitable for contiguous missing intervals.

A critical aspect in data imputation is preserving the fundamental statistical properties of the time series. The comparative analysis of the mean, variance, skewness, and kurtosis between the observed and imputed residuals reveals important findings.

The mean of the imputed residuals is practically identical to that of the observed residuals in all cases, demonstrating that the model is essentially unbiased. However, it is consistently observed that the variance of the imputed residuals is lower than that of the observed residuals. This behavior is typical of regression based models, which tend to predict toward the conditional mean and therefore do not fully capture the complete high frequency variability present in the original series. As a result, the imputed residuals exhibit a slight reduction in variance compared to the observed residuals. This tendency has been widely documented in predictive modeling and data imputation contexts [30].

The preservation of higher order statistical moments shows more variable behavior. In general, the model correctly captures the direction of the skewness observed in the residual series; however, its magnitude is often slightly amplified in the imputed data.

Despite the promising results obtained in this study, several limitations should be acknowledged. First, the training strategy based on fixed five year windows may limit the model’s ability to capture long term climatic variability and regime shifts that occur at larger temporal scales. Second, the proposed approach was evaluated using a limited number of tide gauge stations, which may restrict the generalizability of the results to other coastal regions with different oceanographic conditions.

Although these limitations, the proposed framework offers a scalable and physically consistent solution for large scale sea level reconstruction and missing data imputation across diverse coastal environments, particularly due to its reliance on freely available data sources such as NOAA and ERA5.

Finally, the model shows a strong dependence on the combined use of meteorological and cyclical variables, since removing either component reduces the predictive accuracy by approximately 40%. Future research should explore adaptive training windows, the incorporation of additional physical predictors, and the evaluation of the method across a broader network of tide gauge stations.

6. Conclusions

The results obtained in this study strongly suggest that the proposed methodology provides an effective framework for reconstructing incomplete sea level time series. the method captures both the deterministic tidal signal and the tidal variability by combining harmonic analysis with a multilayer perceptron model applied to the residual component. Overall, the proposed approach shows clear performance advantages over traditional methods reported in the literature.

The use of freely available and highly reliable data sources such as NOAA and ERA5 is a strong point of the proposed methodology since it improves the reproducibility of the experimental portion of this work and its possible use over different oceanic regions.