Abstract

Tidal variations between three stations located in the Wouri estuary in Cameroon over a semi-decadal period are examined to more accurately characterize tides in the area. A unified harmonic analysis method using multiple software platforms allows for precise characterization of tidal components as well as more reliable long-term prediction. The iteratively reweighted least squares method with a Cauchy weighting function ensures resistance to outliers and increases the overall accuracy of the modeling. The models explain approximately 100% of the total variance, demonstrating exceptional data reconstruction. The high spectral resolution, obtained through 68 constituents, allows for a detailed representation of the tide. The reliability of the results is reinforced by the use of the Monte Carlo method to estimate confidence intervals, which are colored by a unitary Lambert Scargle oversampling factor. Their narrowness indicates high precision in estimating amplitudes and phases. An automatic adjustment is performed with a unitary minimal threshold without inclusion of trend and without application of pre-filtering correction. Thus, the method effectively adapts to irregularly spaced data, characteristic of tidal observations. Finally, the application of exact nodal and satellite corrections improves the accuracy of long-term predictions, especially since the phase shift for the Greenwich constituent is determined from astronomical arguments. Calculated residuals and statistical tools such as correlations allow for comparative analysis and visualization of data in temporal and spatial domains. The prediction models are finally validated by statistical analysis, and by comparison with real data acquired a posteriori and with other prediction sources such as SHOM. The absence of inference reinforces the credibility of this validation, as the models are based on direct observations rather than extrapolations. The results show better accuracy of our models and also that the tides for the three stations have similar trends and patterns with strong correlation.

1. Introduction

The Port Authority of Douala (PAD) in Cameroon is subject to tidal variations due to its proximity to the Atlantic Ocean. Sea level oscillations depend on several frequencies that need to be accurately identified and characterized. A harmonic analysis of the tide in the port of Douala can help understand fluctuations in water levels and currents, model them more effectively to make more accurate and reliable predictions. Knowledge of tides is crucial for navigation and port operations, infrastructure planning, as well as for ensuring maritime safety and managing risks. It also allows for the study of coastal ecosystems and the prediction of extreme events, thus optimizing port activities while protecting the environment and ensuring the safety of coastal areas. The results of harmonic analysis can also be used to study long-term variations in tidal levels and interactions between different tidal constituents, which can help better understand tidal dynamics in the port of Douala. These results can be scientifically interpreted in terms of sea hydrodynamics and its response to tidal forcing and provide parameters that can be mapped to describe the tidal characteristics of the region and to verify the performance of tide gauges ¹. This finds all its relevance in the fact that the main challenges related to the use of harmonic tidal analysis include the quality of collected tidal data and the complexity of the method ².

Recent studies of tides in the Wouri estuary are rare. The technical report by Koutitonsky V. ³ presents the installation of water level sensors and other oceanographic sensors in the Wouri estuary to improve the environmental monitoring system of the Port Authority of Douala. It proposes recommendations on the type and position of sensors needed for the system, as well as the development of a hydrodynamic model to predict water levels in real time. Moreover, the methodology of hydrodynamic modeling is not detailed and the results have not been confirmed by independent observations or measurements. The study by Onguene et al. ⁴ on the other hand, offers an overview of tidal characteristics along the entire Cameroonian coast, based on data collected from 2007 to 2012. It uses harmonic analysis and is mainly limited to the identification of major constituents and their amplitudes, focusing on the spatial distribution of tidal patterns across different stations along the coast. In addition, the study faces data gaps and uneven distribution of measurement stations. Fossi et al. ⁵ offer a broad historical perspective, dating back to the 19th century, cover a wider field of study including sediment dynamics, saline intrusion, impact on climate change, and develop a 3D hydro-sedimentary model. In contrast, in Fossi ⁶, the author analyzes the evolution of water heights in the estuary between 2002 and 2019 from data from three tide stations, focusing on extreme tide levels and annual averages to assess the impacts of sea level rise on coastal erosion and sedimentation. However, the study is exploratory in nature and is also based on discontinuous data. The present study focuses specifically on the Wouri estuary, targeting a more detailed and contextualized analysis of local tidal dynamics. The use of advanced harmonic analysis methods and the integration of multiple sophisticated software also help to fill gaps and obtain better temporal coverage, improving the accuracy of tide predictions. From sea level measurements acquired through radar tide gauges installed at stations SM2, SM3 and SM4, from January 1, 2018 to April 30, 2023, we attempt to numerically determine, through a unified analysis, the amplitudes and phases of the different harmonic tidal components. A unified tidal analysis being mainly motivated by the need to exploit a multi-year sequence of observations and measurements collected at irregularly spaced times ⁷. By comparing the five-year data from different stations, we try to identify tidal trends and patterns that are common to all these stations, as well as differences that may be related to local factors such as topography or currents. This comparative study uses statistical analysis, correlation analysis and regression analysis to determine relationships between different variables and achieve spatio-temporal visualization of tidal variations. Many tidal prediction methods exist including numerical models and statistical methods ⁸. Classical tidal harmonic analysis, based on least squares estimation, has been enriched by contributions from Godin ⁹ with the grouping of constituents for close frequencies, automated by Foreman and Henry ¹⁰, and improved by Leffler and Jay ¹¹ for better resistance to broadband noise. Empirical and statistical methods, such as the empirical correction of Hibbert et al. ¹², the Kalman filter of Yen et al. ¹³, and the ARMA model of Kaloop et al. ¹⁴, offer flexible approaches to handle irregular tidal data and perform short-term predictions. Frequency estimation-based methods, including simple ¹⁵ and multivariate ¹⁶ harmonic least squares estimation, and the method of inaction ¹⁷, aim to identify tidal frequencies without prior knowledge, allowing for more detailed analysis in complex marine environments. Wavelet-based methods, including continuous wavelet transform ¹⁸ and applied by Guo et al. ¹⁹, as well as wavelet decomposition for improving prediction models proposed by Seo et al. ²⁰, Meng et al. ²¹ and Shiri et al., ²², offer powerful tools for analyzing non-stationary and non-linear tidal signals. Artificial intelligence methods include artificial neural networks ^{23, 24}, genetic programming ²⁵, machine learning ²⁶, and adaptive neuro-fuzzy inference system (ANFIS) ²⁷, offering advanced approaches to model complex and non-linear relationships in tidal data. Despite all these developments, several authors agree that harmonic analysis allows for predicting tidal levels with high accuracy ^{2, 28, 29, 30, 31}.

2. Materials and Methods

Obtaining real-time tide gauge data and tide forecasts are major challenges for the Dredging and Maritime Logistics Department (DDLM) of PAD, for the Delegated Dredging Authority (RDD) as well as for port activity in general. Hence state-of-the-art equipment has been installed on existing tide stations to send the signal to a Production, Reception and Storage Center (CPRS) located at the PAD hydrographic service, and to the hydrographic survey boat M/V Hydro for real-time tide level acquisition. An extension of this network was recently carried out to increase the possibilities of oceanographic parameters in tide stations by initially integrating wind data. The four tide stations of the PAD Hydrographic Service are as follows (Figure 1): SM1, offshore from the estuary mouth, SM2 at the mouth of Moucouchou Bay, SM3 at the mouth of Modeaca Bay and SM4 in front of the Port of Douala ³. The tide gauge data transmission and reception network is a point-to-multipoint type network covering the entire estuary. As station SM1 is not yet fully rehabilitated, this network is therefore composed of the following active nodes (Figure 2):

•The three tide stations SM2, SM3 and SM4 whose core is a radar water level sensor with WIFI telemetry, allowing acquisition, recording and real-time remote data transmission. The system is rated IP67 suitable for harsh environments.

•The hydrographic survey boat with WIFI telemetry allowing real-time remote tide reception, recording and integration into hydrography software. The vessel can receive all stations.

•The CPRS reception hub which receives real-time data from all stations via UHF WIFI telemetry with the antenna on a high tower to improve reception. This allows remote monitoring, data extraction and exploitation as well as backup. The network is coordinated at this node. Station discrimination is performed by an identifier-based multiple access system.

Each station consists of a rectangular building housing the water level measurement, acquisition and data transmission modules. Solar panels installed above the shelter coupled with a powerful battery and charge controller ensure uninterrupted power supply to the equipment (Figure 3 and Figure 4).

A UHF radio wave system allows real-time data transmission to the CPRS. The stations are equipped with radar sensors for water level measurement, offering an accuracy of ±3 mm. These sensors use pulse radar technology to measure the distance between the sensor and the water surface.

For all stations, we have an average of 37,380 tidal data values. These data are irregularly sampled due to gaps, with an average time interval of 1.3 hours, ranging from 1.2 to 5.4 hours for SM2, 1.2 to 4.7 hours for SM3, and 1.2 to 5.1 hours for SM4. The covered period extends over 5.33 years, from December 31, 2017, at 00:00 to April 30, 2023, at 23:55 (GMT+1).

The reference point for the chart datum (CD) of each tide station is located at the top of the stilling well inside the station. The relevant data are presented in Table 1.

Harmonic analysis software is based on a decomposition of the tidal force potential developed by Doodson ³². The existing differences between these software packages lie primarily in the definitions of harmonic components and the number of elementary waves considered in the analysis ³³. The package developed by the University of Hawaii Sea Level Center (UHSLC) has been used in the region. Others include Tidalfit, TASK-2000 developed by the Permanent Service for Mean Sea Level (PSMSL), the T_Tide package ³⁴ developed by Professor Pawlowitz of the University of British Columbia ³⁵, VAV ³⁶, and U_Tide, which is an evolution of T_Tide ³⁷. The nodal corrections applied in T_Tide are designed such that they cannot be used accurately for measurement sets of more than one year. U_Tide uses nodal corrections that allow for multi-year time series processing, hence its preference. Tidalfit and UTide both rely on harmonic analysis and optimization by the least squares method. However, in practice, we observed that Tidalfit's nodal correction factors allow for better reconstruction of data gaps, while UTide enables more accurate tide predictions over additional multi-year periods while providing a richer statistical analysis platform. However, these software packages assume that the data have been pre-processed, formatted, and filtered, and therefore do not integrate these aspects of data exploitation. As such, all pre-processing will be performed using the MMDataProcessor software developed in our laboratories, and the rest of the processing will utilize a complementary combination of Tidalfit and UTide utilities.

Data recorded continuously since 2018 for the stations are often subject to problems related to malfunctions in acquisition systems, transmission or reception, or external interventions (maintenance, storms, floating objects, etc.). These problems can take the form of jumps, gaps, or spikes in water level recordings. Therefore, pre-filtering is performed to eliminate outliers using MMDataProcessor software, followed by data pre-processing in Matlab. Time-stamped data are sorted and classified by unique, uniformly increasing temporal labels. Small spikes and gaps will be replaced by linear interpolation. Large gaps will be filled using partial harmonic analysis.

Tidalfit and UTide software are used to calculate the amplitudes and phases of different tidal components for each station. With small data gaps and discrepancies already interpolated, TidalFit software is utilized for initial harmonic analysis and generation of a localized tidal model to interpolate large data gaps and obtain complete and clean time series. To avoid memory issues and optimize numerical processing costs, data subsampling is performed by decimation to hourly observations. The hourly tidal data from TidalFit serve as inputs for UTide tidal analysis software for final harmonic analysis and calculation of parameters such as Mean Sea Level (MSL) using the Demerliac filter, Lower Low Water (LLW), Mean Lower Low Water (MLLW), Lowest Astronomical Tide (LAT), and Highest Astronomical Tide (HAT). Specifically, UTide will successively calculate tidal coefficients, select constituents for harmonic prediction, calculate nodal correction factors, astronomical constants and arguments, calculate current ellipse parameters from cosine-sine coefficients, compare calculated residual tides using a tidal filter with residual tides calculated by subtracting measured tides from harmonically predicted tides, and compare astronomical prediction with measured tide. The process concludes with a statistical evaluation, notably the estimation of auto-spectral and cross-spectral densities based on direct and cross Lomb-Scargle periodograms, evaluation of auto-correlations, cross-correlations, and regression analysis. The choice of an Iteratively Reweighted Least Squares (IRLS) method with a Cauchy weighting function is explained by the desire to ensure resistance to outliers and increase the overall precision of the modeling. The reliability of the results is reinforced by the use of the Monte Carlo method to estimate confidence intervals, which will be colored through a unitary Lambert-Scargle oversampling factor. An automatic adjustment is then performed with an adequate minimal threshold, without trend inclusion and without application of pre-filtering correction. Finally, exact nodal and satellite corrections are applied, taking into account the latitudes of the stations to significantly improve the accuracy of long-term predictions, making this approach particularly suitable for analyzing and forecasting tides over long periods. The phase shift for the Greenwich constituent is determined from astronomical arguments, and no inference is made so that each harmonic constituent is directly determined from the observed data, thus reinforcing the reliability of the results.

The general equation for the height h of the tide at any time t can be harmonically represented by the formula ^{28, 35, 38}

(1)

h is the tide height at any time t, H₀ is the height of the mean water level above the used reference, H_i is the mean amplitude of any constituent with index i, f_i is the reduction factor of the mean amplitude H_i for the prediction year, ω_i represents the rate of phase change or angular velocity of the constituent with index i (degrees/hour), t is the time counted from an initial epoch such as the beginning of the prediction year, V₀ is the phase shift of the perturbing celestial body relative to the Greenwich meridian, the quantity (V₀+u) represents the value of the equilibrium argument of constituent A when t=0, and finally k_i is the epoch of the constituent with index i. f_i(t) and u_i(t) are also called the nodal correction factors in amplitude and phase for each component i of the spectrum ². This approach allows for the integration of the effects of nodal corrections, constituent inferences, and astronomical arguments in the one-dimensional analysis of tide gauge data. It can be further written as:

(2)

The residual time series h_r(t) is equal to the observed data series h(t) minus the predicted time series created using the n harmonic constants (H_i and k_i). It contains non-tidal water level variations due to meteorological effects such as wind, atmospheric pressure, and river discharge, as well as tidal oscillations caused by harmonic constants not included in the n resolvable constituents used for prediction. ³⁹. Each cosine term in the equation above is called a tidal constituent or component. The products A_i = f_i H_i are the amplitudes of the constituents and are derived from observed tidal data. The argument of the cosine is a uniformly varying angle, and its value at any time is called its phase. The periods and corresponding speeds of the constituents are derived from astronomical data and are independent of the tide station location. The quantities above, except for h and t, can be considered constants for a particular year and location, and when these constants are known, the value of h, or the predicted tide height, can be calculated for any time value t. By comparing successive values of h, the heights of high and low waters, as well as the times of their occurrence, can be approximately determined.

The least squares method is employed to minimize the residual difference between the data and an optimal tidal prediction ⁴⁰. Considering a point with n temporal instances, this involves minimizing the following cost function:

(3)

The quantity χ2(m), (or J) or sum of squared residuals, is a positive definite quadratic form. Its minimum is calculated by differentiation.

(4)

(5)

However, the analyzed time series must meet sampling and duration criteria. For instance, two very close frequencies will be difficult to separate. The minimum recording or simulation duration τ necessary to separate two waves with frequencies f₁ and f₂ is given by the Rayleigh criterion:

(6)

This relation expresses that two neighboring harmonic components can only be separated if one has at least one complete period more than the other over the temporal extent of the data. Concurrently, if the temporal sampling is insufficient, a loss of temporal information may occur. This is the phenomenon of spectral folding or aliasing, and the observed frequency is lower than the actual frequency contained in the observed signal. This problem arises when processing satellite observations.

Any instantaneous sea level measurement in a series can be considered as the sum of three components:

(7)

The more direct and economical method to remove tidal energy and find the mean sea level is to use filtering ³⁹. One of the oldest and most widely used tidal filters is Doodson's X0 filter ⁴¹ with a duration of 39 hours (i.e., 2ξ+1 = 39). Theoretically, the longer the vector a_ρ of filter coefficients, the more effective the tidal effect reduction filter, but the more limited its application will be in the case of tide gauge observations with gaps. Other low-pass filters exist: the W25 filter, the Munk filter, the Godin filter, and the Demerliac filter using respectively 25, 49, and 71 hours ². The Demerliac filter ⁴², recommended by the French Naval Hydrographic and Oceanographic Service, SHOM ⁴³, is a good compromise reducing less non-tidal energy at subtidal frequencies than the Doodson filter, and using a symmetric vector of 71 elements, hence its preference here. Groves ⁴⁴ used a method of combining simple filters to produce many other tidal filters. There are numerous papers on tidal filtering and many opinions on which filters work best. Emery and Thomson ⁴⁵ have a good discussion on filtering and seem to prefer the Butterworth filter. Various other filters can be experimented, allowing optimization on various high-order statistical parameters as well as on the numerical cost of filtering ⁴⁶.

This will consist of searching for local minima and maxima in the time series. The exact value of t for high and low water times will be the roots of the first derivative of formula (1) equal to zero, which can be written as:

(8)

In addition to constituent amplitudes and phases, harmonic analysis allows us to holistically examine other parameters of interest such as their confidence intervals, signal-to-noise ratios, and residual noise, as well as statistical parameters like correlations. Constituent amplitudes measure the intensity of water height variation due to each constituent. Higher amplitudes indicate a greater impact on the tide. Constituent phases, conversely, measure the advance or delay of water height variation due to each constituent relative to a common reference. Higher phases indicate a greater delay relative to the reference. The precision of these parameter estimates is provided by their confidence intervals. Narrower confidence intervals indicate more precise estimates. Furthermore, the signal-to-noise ratio (SNR) measures the quality of the tidal signal relative to noise that may be due to meteorological, oceanographic, or instrumental factors. Higher SNR indicates a clearer and more easily analyzable signal. Finally, the root mean square error (RMSE) measures the difference between observed values and those predicted by the harmonic model. Lower RMSE indicates a model better fitted to the data. As with other physical phenomena and signals that recur over time, such as fading in mobile radio communications, it is also appropriate for better statistical characterization of tides at target stations to employ second-order statistical parameters such as correlation and plot correlograms ⁴⁷. Autocorrelation is an index to evaluate the repetitive nature of the tide and the auto correlogram can hence reflect regularities and repeated patterns in the tidal signal while allowing detection of cycles or periodicities in tidal variations, even in the presence of noise or disturbances. This can also help understand how tidal heights evolve over time and predict future tides based on historical data. Examination of cross-correlograms of tide stations in pairs allows determination of time lags where data series become uncorrelated. At these specific lags, there is no longer a significant linear relationship between tidal measurements from the two stations, and tidal variations recorded by one station can no longer be predictive of variations recorded by the other. This information is particularly important in tidal data analysis as it can indicate a limit in the coherence of tidal models between two stations, and reflect the influence of local factors that disrupt their correlation. This can be used to improve tidal prediction models by accounting for site-specific characteristics. To achieve good precision in detecting tidal components, the analysis must be performed on a signal sampled over a sufficient time period, minimum 29.5 days, corresponding to a lunar month. The longer the measurement period, the more representative the results will be of tidal components. Generally, 19 years of measurements allow for optimal analysis, but as it is very rare to have quality tide gauge data over such a long period, one year of measurement is required to have good precision of harmonic components and for satisfactory prediction quality for navigation needs. However, when the tidal wave has progressed over shallow depths over large distances, as is the case in the Wouri estuary, calculation of certain harmonic components, resulting from non-linear interactions, requires more than four years of observations ². Hence our choice to base our study on semi-decadal data sampled every five minutes. Harmonic prediction is accomplished by combining elementary constituents according to astronomical relationships prevailing at the time predictions are made.

- Reconstruction Quality Parameters

These parameters include the mean error, variance, standard deviation, minimum value of R Rmin, the weighting factor K, the overall signal-to-noise ratio, the total variance explained by all constituents TVallc, the model's signal-to-noise ratio TVsnrc, the gross total variance TVraw, the percentage of total variance explained by the model PTVallc, and finally the percentage of total variance of the signal-to-noise ratio PTVsnrc. They are calculated according to the formulas proposed by Codiga et al. ³⁷. Correlations are also included.

- Correlations

The autocorrelation function depends on the lag τ and the autocorrelation coefficient of order τ, r_τ for tidal data, where the components have amplitudes A_i and phases g_i, can be expressed as follows:

(9)

This formula is used to evaluate the similarity between tidal components at different time intervals and is particularly useful for identifying periodic patterns and temporal dependencies. The most important of the r_τ is r₁. In the case of studying multiple tide stations, correlations between stations are calculated from the complex coefficients of the constituents, using the following formula:

(10)

where r_νξ is the correlation coefficient between stations ν and ξ, n is the number of constituents, A_νi and g_νi are respectively the amplitude and phase of constituent i for station ν.

It relies on a combination of harmonic regression, covariance matrix calculation, random sample generation, and quantile determination. Residuals are calculated by the formula:

(11)

and the variance of residuals by

(12)

M is the total number of observations or data points and p is the number of estimated parameters. For tide station ν, the covariance matrix of regression coefficients can be estimated from the least squares matrix:

(13)

Where is the vector of parameters estimated by harmonic regression, X_ν is the matrix containing the cosines and sines of frequencies. The diagonal elements of the covariance matrix allow calculation of the variances of amplitudes and phases of constituents. Random samples of model parameters using their estimated distribution, based on parameter covariance, are repeatedly generated to obtain an empirical distribution of parameters.

(14)

is the vector of parameter samples. From the generated samples, confidence intervals can be calculated using the quantiles of the sample distribution. For a 95% confidence interval, quantiles can be determined as follows:

(15)

where Q_(μ,ν) represents the μ^th quantile of simulated samples for station ν.

3. Results

The harmonic analysis allowed for each station to determine the main components which are waves. The most important of which have semi-diurnal and diurnal periodicities. The obtained models comprise 68 constituents (Table 2) that allow for a detailed representation of the various tidal components, thereby improving the precision of the reconstruction. The table includes the name or symbol of the tidal constituent, the amplitude of the constituent A in meters, the confidence interval of the amplitude A_ci, the phase of the constituent g in degrees, the confidence interval of the phase g_ci, the type of tide, the period T in hours, the origin of the constituent, and its description.

These columns provide essential information to characterize each harmonic tidal constituent identified at the stations, including its amplitude, phase, period, and other relevant details. Figure 6 to Figure 8 present the reconstruction of tidal levels for the stations, based on the results of the harmonic analysis. They show a comparison between the observed tidal levels (upper curves) and those reconstructed by the harmonic model (lower curves) over the same period. It is readily observable that the reconstruction has eliminated both minor and major data gaps and accurately reconstituted the tide. These figures thus illustrate the efficiency of the harmonic models, and particularly the method employed, in reproducing the observed tidal variations at the stations.

4. Analysis and Discussions

The constituents representing the tidal model for each station are listed and detailed in Table 2. The principal tidal constituent for stations SM2, SM3, and SM4 is M2, the principal lunar semi-diurnal tide, with amplitudes of 0.721m, 0.765m, 0.802m and phases of 151°, 161°, 173° respectively. The intermediate amplitudes and phases of SM3 between those of SM2 and SM4 indicate an intermediate geographical and temporal position, with a stronger M2 tide at SM4. The SSA constituent, representing the principal solar annual tide, has the smallest amplitude at SM3 (0.0386m), suggesting a weaker tide at this station. The S2 constituent, corresponding to the principal solar semi-diurnal tide, presents amplitudes of 0.231m, 0.241m, and 0.243m and phases of 188°, 203°, and 219° for SM2, SM3, and SM4 respectively, indicating a similar but slightly stronger solar influence at SM4.

Other important semi-diurnal constituents of SM2 (K2, N2, L2, etc.) have amplitudes ranging from 0.0117m to 0.0598m and phases from 93.3° to 209°, reflecting interactions between M2 and S2 as well as effects of the lunar orbit and ecliptic. The amplitudes and phases of these constituents are similar for SM3 and SM4, except for M4 at SM3 and MU2 at SM4, which show variations due to specific interactions and local topography. Important diurnal constituents of SM2 (K1, P1, O1, etc.) have amplitudes from 8.28x10^-6m to 0.136m and phases from 12.4° to 340°, influenced by Earth's rotation and both lunar and solar declination. The amplitudes and phases of diurnal constituents are similar for SM3 and SM4, except for K1, which has a higher phase at SM3 (26.9°) and lower at SM4 (34.1°), likely due to longitude differences. Long-period constituents (SA, SSA, MF, etc.) at SM2 have amplitudes from 0.000186m to 0.06720m and phases from 10.9° to 360°, related to seasonal and annual variations. The amplitudes and phases of these constituents are similar for SM3 and SM4, indicating low sensitivity to local variations. MO3 has a lower amplitude at SM3 (0.01340m) than at SM2 (0.00810m), suggesting a more attenuated lunar tri-diurnal tide at SM3. Important non-linear constituents of SM2 (M4, MS4, MN4, etc.) have amplitudes from 1.12 x10^-5m to 0.0376m and phases from 5.94° to 328°, influenced by interactions between linear constituents and effects of bathymetry and resonance. The amplitudes and phases of these constituents vary for SM3 and SM4, indicating a strong influence of local conditions. For example, MS4 has higher amplitudes at SM3 and SM4 (0.0490m and 0.0632m) than at SM2 (0.0165m), and M4 is more amplified at SM4 (0.0993m) than at SM2 (0.0376m). The semi-diurnal constituents M2 and S2, as well as the diurnal constituents K1 and O1, have the largest amplitudes in all three stations, which is typical for many coastal environments. The amplitude hierarchy is M2 > S2 > K1 > O1, with M2 having a significantly higher amplitude than other constituents. The amplitudes and phases of the main constituents are relatively consistent across the three stations, indicating similar tidal propagation in the estuary. Regarding spatial variations, slight differences in amplitudes and phases between stations are observed, which may reflect local effects of bathymetry and estuary geometry. For instance, the intermediate amplitude and phase of SM3 between those of stations SM2 and SM4 indicate that station SM3 is located between the other two stations geographically and temporally. Generally, the differences attest to the time lag between the three stations and that the principal lunar semi-diurnal tide is stronger at SM4 than at SM2 and SM3. Minor constituents including higher-order constituents and shallow water constituents generally have lower amplitudes but can play an important role in the detailed shape of the tide. The phases of the main constituents also show a coherent progression between stations, which is indicative of the propagation of the tidal wave in the estuary. This analysis shows that the tide in the Wouri estuary is primarily semi-diurnal with significant diurnal influence, and that tidal characteristics are relatively consistent between the three stations while presenting subtle spatial variations.

The analysis of confidence intervals (Table 2) shows that they (A_ci for amplitude and g_ci for phase) are generally very narrow for most major constituents in all three stations. This indicates high precision in the estimation of amplitudes and phases. The principal semi-diurnal (M2, S2) and diurnal (K1, O1) constituents present particularly narrow confidence intervals (CI), suggesting high reliability in their estimation. For example, the M2 constituent has a confidence interval of 0.000474m for amplitude and 0.0396° for phase at station SM2. A consistency in the narrowness of CIs is observed between the three stations for major constituents, indicating good overall data quality and analysis. Lower amplitude constituents tend to have relatively wider confidence intervals, which is normal given their lesser contribution to the total signal. For instance, constituents K2, MU2, L2, MSF, SSA, MM, MF, NU2, 2N2, LDA2, O1, M6, S4, M3, J1, Q1, and OO1 have wider confidence intervals at SM3 than at SM2 and SM4. This means that the parameter estimates for these constituents are less precise and less reliable at SM3 than at the other stations. Regarding overall precision, the general narrowness of confidence intervals (CIs) for most constituents across all three stations suggests that the harmonic analysis method employed (IRLS with Cauchy weighting function and Monte Carlo method) was highly effective in producing accurate estimates.

Figure 9 displays the auto correlograms of the time series associated with each tide gauge. The series exhibit a high degree of autocorrelation between adjacent tidal observations. Moreover, the shape of the correlograms shows a single recall to zero. This confirms the non-random nature of the underlying time series. Figure 10 and Table 3 present the first-order correlation coefficients obtained for the three stations. All converge to 1 at the origin. Furthermore, the curves are all positive, and the slow decrease of the correlation coefficient from the value 1 at the origin towards zero as the lag increases denotes data persistence. The correlograms also exhibit different peaks, even though they have similar periodicities, indicating a slightly different manifestation of the tide depending on the geographical position of the stations. The three stations present similar autocorrelations, confirming periodicity and concordance between stations. The decorrelation time, which indicates the duration over which tidal data are statistically linked, is approximately a quarter of a period for all three stations, or 3 hours (more precisely 2 hours and 59 minutes for SM2 and SM3, and 3 hours and 05 minutes for SM4) (See Figure 9). This may prove useful in understanding the dynamics and for predictive modeling of tides in the region.

The correlations between stations are given in Table 3 and (Figure 10.a). They are also similar, with the strongest observed between SM2 and SM3 and the second strongest between SM3 and SM4. At first glance, the observation of SM2 and SM3 suggests that SM3 is out of phase with the other two stations, but the cross-correlation between SM3 and SM4 eliminates this presumption. The same observation challenges a constant offset of one station from the others. The time lags at which the tidal data series become uncorrelated are almost similar for the SM2-SM3 and SM3-SM4 pairs (1 hour and 24 minutes and 1 hour and 26 minutes respectively) and slightly higher for the SM2-SM4 pair (1 hour and 34 minutes) (Figure 10.b). Overall, after an hour and a half, there is no longer a significant linear relationship between tidal measurements from the stations, and tidal variations recorded for one station can hardly be predictive of variations recorded by either of the other two stations. These high correlations confirm the spatial coherence between stations and the similarity of tidal characteristics in the estuary.

Examination of the diagnostic table data (Table 4, Table 5, and Table 6) reveals that the M2 component has the highest percentage of energy (PE) (>81%) and the highest signal-to-noise ratio (SNR) (>4x10⁶) across all stations. This elevated SNR suggests that the M2 estimation is highly precise. Furthermore, the S2, K1, N2, and SA components are also present in the upper quartile for all stations, although their order varies. This indicates that these components are the most significant in these three stations. Additionally, only constituents with sufficiently high SNRs (SNR≥2 dB) were retained in the final analysis. The low-frequency (loNAME) and high-frequency (hiNAME) constituents correspond to the neighboring frequencies (lower and upper) of the principal frequency. For illustration, M2 has H1 and H2 as its lower and upper components, which are components with frequencies slightly lower and higher than that of M2. The Rayleigh Ratios for the lower and upper components (loRR, hiRR) measure the sensitivity of the amplitude and phase of the principal component to variations in neighboring frequencies. High ratios indicate precise separation and analysis of each constituent's contribution to the overall tidal signal.

The M2 component has response ratios of 5.3 for H1 and H2, indicating low sensitivity to variations in these components. The noise-modified variants of the Rayleigh ratios, or noise ratios of the lower and upper components (loRNM, hiRNM), measure the variability of the principal component's amplitude and phase due to noise from neighboring frequencies. Higher RNM values indicate greater variability of the principal component. For example, M2 has noise ratios of 7.06 x10³ for H1 and H2, signifying low variability due to these components. The maximum correlations of the lower and upper components (loCorMx, hiCorMx) measure the degree of linear relationship between the amplitude and phase of the principal component and those of neighboring frequencies. High correlation between parameters may indicate multicollinearity issues that could affect the stability and reliability of model parameter estimates. Higher CorMx values indicate stronger relationships. For instance, M2 has maximum correlations of 0.04 for H1 and H2, indicating a weak relationship with these components. For all three stations, the M2 constituent is predominant, representing 81.75% of the total energy for SM2, 82.12% for SM3, and 81.75% for SM4, with respective signal-to-noise ratios (SNR) of 4.10x10⁶, 4.70x10⁶, and 4.10x10⁶.

The correlation matrices for low and high frequencies exhibit low values for M2 across all stations, indicating minimal multicollinearity and good independence of constituents. The Rayleigh ratios for M2 are similar, ensuring effective frequency separation. Constituents S2 and K1 also demonstrate notable contributions across the three stations, with comparable energy percentages and SNRs. However, SM3 slightly distinguishes itself with a higher energy contribution from S2 (8.17%) compared to SM2 (7.50%) and SM4 (7.50%). Overall, these diagnostics underscore the significance of M2 in tidal modeling and the stability of the models employed for all three stations.

The case of SM2 is considered here for illustration purpose. SM3 and SM4 exhibit similar graphs.

The first diagnostic graph of the fit (Figure 11.a: illustration for SM2) is utilized to evaluate the goodness of fit of the regression model and displays the dispersion of residuals following the adjustment (using model predictions) as well as a ranking of constituents in ascending order of frequency (Figure 11.b).

Figure 12 presents the second diagnostic graph for SM2, illustrating the constituents ranked in decreasing potential energy (PE) order (Figure 12.a) and their associated phase shifts with 95% confidence intervals in panel (Figure 12.b). This graph enables a detailed analysis of constituent amplitudes, where the error bars signify the 95% confidence interval for each amplitude measurement. A signal-to-noise ratio (SNR) exceeding 2.0 indicates that a constituent is reliably detected and distinguishable from background noise; such constituents are prioritized for inclusion in tidal reconstruction efforts. The vertical axis of the graph displays the intensity of these constituents, arranged in descending order of energy, thereby facilitating a clear understanding of their relative contributions to tidal dynamics.

Figure 13 presents a diagnostic graph of the fit showing the residual curve for SM2, which represents the difference between observed and modeled values.

The residuals are generally low, indicating a precise fit of the harmonic model to the tidal data. The residual curve is well distributed around zero, suggesting the absence of systematic bias and random errors. Some extreme values are present, but their limited number demonstrates that the IRLS method with a Cauchy weighting function has effectively minimized their impact. Finally, the low temporal variability of the residuals indicates stability and consistent performance of the model throughout the analysis period. The curves are similar for SM3 and SM4

Table 7 presents a comparative summary of tidal reconstruction metrics for stations SM2, SM3, and SM4. The parameters include mean error, variance, standard deviation, Rmin, K (weighting factor), SNRallc, TVallc, TVsnrc, TVraw, PTVallc, and PTVsnrc.

The harmonic analysis of tides for the three stations, with methodological similarities (Rmin=1.00, K=2.07, MinSNR=2.00), again demonstrates very high reconstruction performance, but notable differences in results. SM2 exhibits the best overall performance, with the highest global signal-to-noise ratio (5.99x10⁸), the lowest mean error (0.01280), the smallest variance (0.000267), and the lowest standard deviation (0.01640). These indicators confirm a more precise and less dispersed reconstruction for SM2. Conversely, SM4 displays the least favorable values, indicating a less precise and more variable reconstruction. SM3 presents intermediate values. The total variance increases from SM2 to SM4, indicating increasing variability in water levels as one approaches the port of Douala. Despite these differences, the percentage of total variance remains exceptionally high for all stations (≥99.8%), demonstrating the excellent reconstruction capability of the harmonic model, with a slightly better precision performance for SM2 (99.9%) compared to SM3 and SM4 (99.8%). These results corroborate previous observations on SNR and total variance, confirming a gradation in tidal reconstruction quality from SM2 to SM4, likely due to local differences in the hydrodynamic conditions of the Wouri estuary.

Table 8 and Figure 14 show that tidal levels, low and high waters for stations SM2, SM3, and SM4. They exhibit significant variations.

The mean sea level (MSL) above chart datum is relatively similar for the three stations, with SM2 having the highest MSL (1.722 m) and SM3 the lowest (1.588 m). SM2 presents the largest tidal range, with a difference of 2.212 m between the Lowest Astronomical Tide (LAT) and the Highest Astronomical Tide (HAT). LAT and Mean Lower Low Water (MLLW) values increase from SM3 to SM4, suggesting an influence of topography along the estuary.

The MSL-LAT offset (ML Zo) also follows this trend, increasing from SM2 (1.327 m) to SM4 (1.461 m). This is consistent with the increase in MLLW. These differences between stations have important implications for navigation, port management, and coastal infrastructure planning, highlighting the importance of site-specific analysis for accurate tidal predictions.

The mean sea level is not constant throughout the year and exhibits a cyclical evolution with an annual period (Figure 15, Figure 16 and Figure 17). It experiences similar variation patterns throughout the year for all three stations.

It remains nearly constant from January to April, then undergoes a decrease of approximately 5 cm to reach its lowest level at the beginning of June. From June to October, the variation trend is uniformly increasing, with the peak occurring in early October and an excursion of 20 cm. Beyond this, the mean sea level experiences a uniform decrease until December. Paradoxically, contrary to what geographical proximity might suggest, SM2 and SM4 exhibit better congruence, while the least congruence is observed between SM2 and SM3.

Similarly, the relative offsets between stations are not constant but periodic and dependent on tidal range amplitudes (Figure 18). The smallest and most regular offset amplitudes are observed for the SM2-SM3 pair. The largest offset amplitudes are observed for the SM2-SM4 pair, followed by SM3-SM4. These amplitudes are more strongly linked and proportional to tidal range amplitudes, with large offset amplitudes corresponding to large tidal ranges.

For the purpose of validating the models derived from harmonic analysis for the three stations, comparisons were made on one hand between the model predictions and the prediction from the Tides & Currents software, between the model predictions and the prediction from SHOM ⁴⁸, and on the other hand between the model predictions and the actual tide acquired a posteriori for various periods of 2024 depending on data availability.

The predicted tide used here is produced by the Tides & Currents software for the geographical area of Douala with the following reference coordinates:

DOUALA (WGS84): Latitude: 04° 03' N

Longitude: 009° 41' E

We observe a concordance in terms of harmonics between the stations and a concordance with the predicted tide for all three stations SM2, SM3 and SM4, thus validating the construction model used (Figure 19). The best in connection with the predicted tide is given by SM4 over the entire duration of the observations, surely due to its geographical proximity. In addition, the maximum amplitudes show offset compared to the prediction.

Another comparison was made through correlations. The cross-correlations between the model predictions and the Tides & Currents prediction are given in Figure 20 and Table 9.

They are similar with the strongest observed for SM4, surely due to its geographical location within the port of Douala followed by SM3 and finally SM2.

Subsequent acquisitions of actual tide data in 2024 allowed for another validation approach by comparison with our model predictions and the predictions of SHOM. The graphs in Figure 21 and Figure 22 show good agreement between the model predictions and the actual tide obtained retrospectively from March 21 to 22, 2024 for SM2 (Figure 21.a), June 13, 2024 for SM4 (Figure 21.b), and February 28, 2024 for both stations simultaneously (Figure 22). We also observe that the model best agreeing with SHOM predictions is that of station SM3.

These predictions are also better than those of the Tides & Currents software in both cases. SM4, which in the previous paragraph showed the best correlation with the software prediction due to its geographical position, presents an equally surprising result: our model's prediction for this station shows the closest agreement with the actual tide curve than that of both the Tides & Currents software and the prediction from SHOM. Thus, the reconstruction results from the models can be used for long-term predictions, thereby aiding navigation, planning of ship arrivals and departures, and optimization of port operations. In this case, they demonstrate the effectiveness and accuracy of harmonic analysis for the three stations, highlighting the model's ability to faithfully reproduce observed tide levels and provide more reliable long-term predictions than other existing alternatives.

5. Conclusion

This five-year multi-station study of tidal variations in the Wouri estuary in Cameroon has allowed for precise characterization of tidal fluctuations and establishment of reliable prediction models for the Douala Port Authority. The unified methodological approach, combining data from multiple stations over a long period, multiple software and analysis tools, has proven particularly effective in analyzing and modeling tides in this complex environment. Robust harmonic analysis, using the IRLS method with a Cauchy weighting function, yielded highly accurate results. The high spectral resolution, obtained through the analysis of 68 constituents over five years of data, provided a detailed representation of tides, explaining about 100% of the total variance. This approach, without resorting to inference, ensured that each harmonic constituent was directly determined from observed data, thus strengthening the reliability of the results. The harmonic analysis was therefore able to resolve most of the higher frequency and tidal components with close periods and confirmed that the tide in the Wouri estuary is a semi-diurnal tide with diurnal inequality that amplifies as it progresses into the estuary. Model validation, performed by comparison with actual data acquired retrospectively and other serious prediction sources, confirmed the accuracy and robustness of the adopted approach. The results highlighted similar trends and common patterns among the three stations studied, with a strong correlation, particularly between SM3 and SM4. However, noteworthy differences were also observed, likely due to local dynamics such as topography and currents. In this regard, it is observed that water levels in the estuary display a certain spatial gradient and consequently, it would be necessary to multiply tide gauges along the estuary to better characterize tidal level variations along the navigation channel leading to the port, in temporal and spatial domains, confirming the interest in rehabilitating station SM1. For future work, it would be relevant to extend this study by integrating other oceanographic parameters, such as currents and winds, to obtain a more comprehensive understanding of the estuary's hydrodynamics. The incorporation of satellite data and the use of machine learning and artificial intelligence techniques could also enrich the analysis and further improve prediction accuracy.

ACKNOWLEDGEMENTS

The author would like to thank the management and staff of Locatech Services Sarl as well as the personnel of the Hydrographical Office of Port Authority of Douala for their kind cooperation as the information and data used in this study were made available through their gracious diligence.

References