1. Introduction
According to Franco (Reference Franco2009), the maritime vertical reference for bathymetric surveys should be a chart datum (CD) corresponding to low water, enabling mariners to assess under-keel clearance directly from charted depths without accounting for tidal variations. A similar principle applies to high-water levels, which ensure safe clearance under bridges and other obstructions. However, depths reduced to these references seldom reflect actual sea level conditions during a vessel’s transit, as the CD is based on one of the worst-case tidal scenarios under average meteorological conditions.
Consequently, if only nautical chart depths are used, the safety margin for draft is often ultra-conservative (Wells et al., Reference Wells, Kleusberg and Vaníček1996). This has direct operational consequences, because even a few centimetres less clearance beneath the keel can equate to tonnes of additional cargo capacity, prompting some vessels to wait for a tidal ‘window’ – the high tide that permits entry with greater draft, occasionally using the fluid mud layer for navigation (Ramos, Reference Ramos2007). Nonetheless, reduced clearance increases the risk of grounding, potentially causing economic loss and environmental damage, in contrast to the global movement towards environmental preservation and a sustainable blue economy (IOC et al., 2011; IOC, 2020).
To address these limitations, many ports have adopted dynamic draft systems, which allow under-keel clearance (UKC) to be managed through real-time calculations based on prevailing environmental conditions during manoeuvres. Electronic Chart Display and Information Systems (ECDIS) presenting S-10X class nautical charts (IHO, 2020) can integrate real-time data on seabed type, tides, wind, currents, waves and updated depths. These systems typically rely on four inputs: (1) predicted sea level based on tide-gauge observations, buoy data and hydrodynamic model results; (2) vessel vertical positioning from high-precision Global Navigation Satellite System (GNSS) antennae; (3) seabed dynamics; and (4) nautical chart depths (FAMOS, 2020). Combined, these provide mariners with an accurate, location-specific depiction of sea level, enabling safer navigation and more efficient port operations.
To standardise the use of dynamic draft systems and ensure the navigational safety, it is essential to quantify the uncertainty associated with all factors contributing to UKC. These factors may be related to the water level, the vessel and the seabed, the latter reflecting the vertical uncertainties of hydrographic surveys and, consequently, the uncertainty in the realisation of tidal CD (PIANC, 2014; Ruggeri et al., Reference Ruggeri, Watai, Rosetti, Tannuri and Nishimoto2018). The operational benefits of this approach are illustrated by the Port of Santos, Brazil, where the implementation of a dynamic draft system in 2013 led to a 24% increase in manoeuvres involving vessels with drafts exceeding 12 m. This underscores the critical importance of precise sea level measurements and an accurate realisation of their vertical reference for optimising port efficiency and ensuring navigational safety (Argonáutica, 2020; Wells et al., Reference Wells, Kleusberg and Vaníček1996; Krueger, Reference Krueger1996; DPC, 2019; CONAPRA, 2020).
This issue is particularly relevant in navigation channels monitored by multiple tide gauges, especially where nearby secondary stations are not fully consistent with the reference station (Simon, Reference Simon2006). In such environments, spatial and temporal inconsistencies in CD, arising from short observation periods collected at different epochs, may compromise the continuity of the available water column. Comparable problems may also occur when CD is revised at an individual tide gauge station without a corresponding update of the bathymetric information portrayed on the nautical chart. Under these circumstances, the reference level adopted for sounding reduction in one section of the navigation channel may differ from that used for tide prediction, as discussed by Mononen (Reference Mononen2017). Furthermore, nautical chart depths are commonly produced using a discrete tidal zoning methodology, which leads to stepwise variations in CD between areas influenced by different tide gauge stations. This further complicates the delimitation of CD in transition zones, as noted by Iliffe et al. (Reference Iliffe, Ziebart and Turner2007) and Robin et al. (Reference Robin, Nudds, MacAulay, Godin, De Lange Boom and Bartlett2016). Taken together, these inconsistencies suggest that CD should be examined not only as a datum concept, but also as a realisation tied to a specific reference epoch. In Brazil, for example, the determination of CD is based on the difference between Mean Sea Level (MSL) and Z0, and is therefore influenced by both the observation period adopted for MSL estimation and the harmonic method used to define Z0. Mean Sea Level is computed over a specific observation period with a minimum duration of 32 days, in accordance with IHO (2024a). The vertical separation between MSL and CD, denoted by Z0, is calculated as the sum of the principal harmonic constituents, following the methods proposed by Courtier (Reference Courtier1939) and Balay (Reference Balay1952a). At Ilha Fiscal, the tide-table station for the Port of Rio de Janeiro, where the tidal regime is semidiurnal with diurnal inequality, Z0 is obtained from the sum of the amplitudes of M2, S2, O1, K1 and P1, yielding a level approximately equivalent to Mean Low Water Springs (MLWSFootnote 1 ).
Figure 1 illustrates the temporal component of CD determination by comparing the parameters used to derive CD at Ilha Fiscal for the currently adopted period (2018–2019) and the previous period (2003–2005). The results show differences of 20.7 cm in So, defined as the elevation of MSL above the tide-gauge zero, and 3.9 cm in Z0. Because water level is referenced to the tide-gauge zero, which differed by less than 1 cm between the two epochs relative to the primary benchmark (RN2-DHN), these discrepancies are mainly attributable to differences in the observation periods and epochs, resulting in a 16.8 cm difference in CD relative to the benchmark.
Comparison between MSL, Z0, and CD for the data periods used to generate the 2010 version of the tide sheet and the 2021 version.
Figure 1. Long description
Panel A: A diagram comparing maritime vertical references for the 2010 version using observations from 2003-2005. The diagram includes labels for RN2-DHN, RN-2953H, RN-93517H, Tide staff, MSL, CD, and various measurements in centimeters. The measurements include 4.3 cm, 25.2 cm, 492.4 cm, 271.1 cm, 339.9 cm, 221.3 cm, 68.8 cm, and 152.5 cm. Panel B: A diagram comparing maritime vertical references for the 2021 version using observations from 2018-2019. The diagram includes labels for RN2-DHN, RN-2953H, RN-93517H, Tide staff, MSL, CD, and various measurements in centimeters. The measurements include 4.3 cm, 24.95 cm, 491.2 cm, 249.2 cm, 321.9 cm, 242 cm, 72.7 cm, and 169.3 cm.
Since CD = MSL – Z0, the observed difference reflects variations in both Mean Sea Level (MSL) and the harmonic constituents, although the results indicate that the change in MSL is the main driver. Such variations may arise from both natural and anthropogenic forcing. Harmonic constituents may evolve over time in response to sea level fluctuations, wind, waves, river discharge, erosion-induced geomorphological changes, dredging, river regulation, reservoir construction and land reclamation (Munk and Cartwright, Reference Munk and Cartwright1966; Woodworth, Reference Woodworth2010; Müller et al., Reference Müller, Arbic and Mitrovica2011; Talke and Jay, Reference Talke and Jay2020). Mean Sea Level also varies under the influence of atmospheric pressure, wind, river discharge and other environmental drivers (Flather, Reference Flather, Steele, Turekian and Thorpe2001; Pugh and Woodworth, Reference Pugh and Woodworth2014).
As recommended by NOAA (2019), uncertainties in the determination of CD from a short time series should be treated as systematic errors in bathymetry derived from hydrographic surveys. This principle may also be relevant to under-keel clearance calculations whenever the CD used for tidal prediction differs from that adopted for bathymetric reduction. Swanson (Reference Swanson1974) proposed a method for estimating CD variability by comparing short-duration records (1–12 months) with those computed from a 19-year time series. His analysis focused on Mean Lower Low Water (MLLW), which is defined as the average of the lower low waters observed over a 19-year period. However, long-term observational records are not always available and CDs are not necessarily defined as averages of observed water levels; they may also be derived from the sum of harmonic constituent amplitudes or from tidal predictions, as in the case of Lowest Astronomical Tide (LAT) or Lower Low Water Large Tide (LLWLT).
In this context, the present study addresses two closely related questions. First, how much can CD vary when it is computed from records of different lengths and at different reference epochs? Second, how much does it vary as a function of the adopted CD concept? To investigate the temporal dimension, two complementary approaches are applied to the MLWS concept. In the first, Z0 is computed from multiple shorter observational time series at Ilha Fiscal, and its variability is quantified from the dispersion of these estimates relative to the trend defined by the short-series results. In the second, Z0 is likewise derived from harmonic analysis, but its variability is estimated from the propagated uncertainties associated with the harmonic constituents used in its computation.
In both approaches, the resulting variability of Z0 is combined with the variability of Mean Sea Level (MSL), which is also evaluated from multiple time series relative to its trend. This procedure can be applied to other regions, assuming that the reference and secondary stations are affected by similar meteoceanographic conditions.
The conceptual aspect of CD is equally important. As noted by Balay (Reference Balay1952b), the problem is not simply to choose a low reference, but to choose one that is low enough to ensure safety without becoming unnecessarily restrictive. If the datum is set too low, charted depths may systematically underestimate the water actually available, creating artificial obstacles to navigation. If it is set too high, the chart may represent conditions that are rarely reached at low tide, increasing the risk when accurate tidal predictions are not available at the moment of transit. Figure 2 illustrates this problem for the current MLWS-based CD adopted at Ilha Fiscal. From 2018 to 2021, sea level remained below the datum 4.3% of the time, corresponding to 47.2 days, and reached a maximum departure of 40.3 cm below the adopted reference. The figure also relates these occurrences to the number of vessel manoeuvres in Guanabara Bay, averaging approximately two manoeuvres per day, providing an operational perspective on the implications of the adopted datum.
Sea-level variability at Ilha Fiscal and its implications for navigation. Top: time series of sea level relative to the current Chart Datum (CD) from 2018 to 2020. Middle left: monthly percentage of occurrences in 2021 during which sea level remained below the current CD, derived from the 2018–2019 reference period. Middle right: number of vessel arrivals and departures in Guanabara Bay in 2021. Bottom: distribution of the magnitudes by which sea level fell below the current CD in 2021.
Source: Multiterminais ( 2022).
The search for an optimal percentage of occurrences in which sea level falls below the CD was recently revisited by Slobbe et al. (Reference Slobbe, Klees, Verlaan, Dorst and Gerritsen2013). However, attempts to define a universal formula for determining CD date back to the International Hydrographic Conference held in London in 1919, when one proposal was to establish it as a fixed offset below Mean Sea Level (MSL) (Franco, Reference Franco1966). This approach was later abandoned from a global perspective. Over the following decades, successive international meetings sought to establish an appropriate vertical reference for depths and tides, culminating in IHO (2025), which presents three recommendations according to the local tidal regime. In oceanic tidal areas, where astronomical forcing predominates over other forcing mechanisms, the Resolution recommends the adoption of Lowest Astronomical Tide (LAT), or of a datum as close as practicable and acceptable to Hydrographic Offices.
Alternatively, another similar datum may be used where low-water levels in a specific area frequently deviate from LAT, while a different datum may also be adopted when established by national policy. LAT is generally regarded as a more conservative datum because it represents the lowest low water predicted over a 19-year period from harmonic constants derived from at least one year of observations. By accounting for long-period tidal variations, including the nodal cycle, LAT reduces both the frequency of sea levels falling below the datum and the occurrence of negative tide predictions in tide tables.
Some countries, such as Canada, adopt a different CD: Lowest Low Water Large Tide (LLWLT), defined as the mean of the lowest predicted low water in each year over a 19-year period (Canadian Hydrographic Service, 2022). In general, one of these annual minima will correspond to the Lowest Astronomical Tide (LAT). Because LLWLT is based on an average rather than on a single extreme value, it is usually less restrictive than LAT and less sensitive to interannual variability. However, as it is generally higher than LAT, its adoption tends to produce a slightly higher number of negative tidal heights in tide tables.
This study shows that the hydrographic problem is not limited to selecting the most appropriate CD concept, as has long been debated within the hydrographic community. Equally critical are: (i) the reference epoch at which the datum is realised and (ii) the observational period used for its computation, since both may affect the resulting value. In some cases, an MLWS realisation referred to one reference epoch and derived from a given observational basis may approach, or even fall below, a LAT realisation referred to another reference epoch and/or computed from a different observation period. This practice is already reflected in some datum frameworks, such as the Baltic Sea Chart Datum 2000 (BSCD2000), in which the reference epoch forms part of the datum specification. These issues have direct implications for navigational safety, for vertical consistency in channels where different CD values have been established using different temporal references and, consequently, for S-100 products, which require the CD used to reference depths to be consistent with that used for water-level prediction (IHO, 2024b).
Against this background, the next section describes the methodology adopted in this study, including the procedures used to assess variability (Sections 2.1 and 2.2) and the methods applied in the harmonic analysis (Section 2.3). The results are then presented in two parts, addressing the temporal dimension (Section 3.1) and the conceptual dimension (Section 3.2), each with its corresponding discussion.
2. Methodology
The Ilha Fiscal tide gauge station is one of the oldest in Brazil, with annual mean sea level records available since 1963 (PSMSL, 2023). However, the longest continuous time series begins in 1989. According to BNDO (2026), from 1989 to 2009, sea level measurements were obtained using an analogue float-and-counterweight sensor with an hourly sampling interval. After 2009, this system was replaced by a digital sensor and the data began to be recorded at 5-minute intervals, with gaps filled by linear interpolation, or by a secondary sensor (radar). The Thompson (Reference Thompson1983) filter was then applied to most of the observations, with weights adjusted to remove high-frequency signals, an approach that offers the additional advantage of allowing the target frequencies to be selected a priori (Silva, Reference Silva1992; Kalil, Reference Kalil1999). Finally, the data were resampled to 60-minute intervals for harmonic analysis.
Regarding the vertical reference, the observations are referenced to the staff gauge zero, whereas the staff gauge was levelled relative to benchmark RN-2. Over the years, a maximum variation of 8 mm was observed in the elevation difference between RN-2 and RN-1, whereas the levelled difference between RN-2 and the staff gauge shows an apparent offset of approximately 10 cm between the pre-1998 and post-1998 records (Table 1).
Temporal variation in the levelled height (cm) differences between RN-2, RN-1 and the staff gauge at Ilha Fiscal.
Table 1. Long description
The table presents the levelled height differences between RN-2, RN-1, and the staff gauge at Ilha Fiscal over various years. It consists of 2 rows and 10 columns. The first row contains the headers: Reference pair, 1993, 1994, 1996, 1998, 1999, 2001, 2005, 2008, 2021, 2025. The second row contains the data for RN-2/staff gauge: 501.7, 500.95, 501.25, 487.15, 491.8, 490.40, 492.40, 492.4, 491.2, 491.4. The third row contains the data for RN-2/RN1: 137.0, -, -, 136.45, 136.5, 136.5, 136.8, 136.8, 137.3, 136.7.
This pattern suggests a possible shift in the vertical reference associated with the staff gauge over time. However, because the exact date of this shift could not be established from the available levelling records, no retrospective vertical correction was applied to the observations. The implications of this limitation for the interpretation of the results are discussed in Section 3.2.1.
Based on this dataset, the methodology adopted in this study, as outlined in Figure 3, was organised into two main components: temporal variability and conceptual variability.
Summary of the analyses performed in this study.
2.1 Temporal variability
The temporal variability analysis focused on the MLWS concept, defined here as the sum of the amplitudes of the harmonic constituents. Two temporal-variability analyses were carried out at this stage. In the first, to investigate the influence of record length, the full time series was divided into non-overlapping subsets of 6 months, 1 year, 2 years, 5 years and 10 years, and the resulting variations in MSL, Z0 and CD were evaluated. For each subset length, three approaches were applied. The first was based on the root-sum-square propagation of the uncertainties of the harmonic constituents used to compute MLWS. The second relied on the standard deviation of the residuals relative to the linear trend, represented by the fitted regression line for Z0, MSL, and CD over time. The third used the standard deviation of each estimate relative to the corresponding 19-year result for the same variables.
Although the three approaches quantify different aspects of variability, they can be propagated to CD through the same Root Sum of Squares (RSS) framework, since CD is obtained from the difference between MSL and Z0. Under the assumption of zero covariance, RSS between these terms (equation 1) provides a consistent way to compare how each uncertainty descriptor translates into CD variability.
$$\sigma {\left( {{\rm{CD}}} \right)_{{\rm{trend}}}} = \sqrt {{\rm{\sigma }}{{({\rm{MSL}})}_{{\rm{trend}}}}^2 + {\rm{\sigma }}{{\left( {{{\rm{Z}}_0}} \right)}_{{\rm{trend}}}}^2}$$
In the second temporal-variability analysis, the series was divided into multiple non-overlapping subsets ranging from 6 months to 10 years, using 1-month (30-day) increments. The uncertainty of Z0 was then compared between T-Tide and Utide, based on the harmonic analysis results for each subset and computed through the root sum square of the uncertainties of the harmonic constituents. Pacmaré was not included because, unlike the other programs, it relies on a graphical and command-line interface that does not allow multiple automated iterations (Section 2.3).
The variability statistics obtained for harmonic constituents, MSL, and Z0 from multiple time series may provide a useful methodological reference for short-term tide-gauge stations. Their applicability to other locations, however, depends on the extent to which local MSL variability, tidal amplitudes and non-tidal forcing are comparable to those at the reference site. Therefore, the main contribution of the present study lies in the estimation framework and associated statistical descriptors, rather than in the assumption that the resulting values are universally transferable. This caution is warranted even within Guanabara Bay, where Andrade et al. (Reference Andrade, Rosman and de Azevedo2019) showed that baroclinic effects may modify mean water levels by up to approximately 5 cm, depending on the sector and season. Accordingly, the variability estimates presented here should be regarded as site-specific and transferred to other stations only with due consideration of local hydrodynamic conditions.
2.2 Conceptual variability
With respect to conceptual variability, two analyses were performed. In the first, three different criteria were adopted to define the period used for MLWS computation: (i) the period yielding the highest PTVallc (Percent Total Variance Calculated), namely from December 2014 to May 2015; (ii) the current 2-year series presently used for Ilha Fiscal (2018–2019) to calculate the CD; and (iii) the complete 19-year series (2002–2021). The maximum-PTVallc criterion represents the period with the smallest meteorological residual, that is, the interval in which the astronomical signal best reproduces the observed series used in the harmonic analysis (Codiga, Reference Codiga2011). The 19-year criterion represents the ideal series, as it spans the full nodal cycle.
In the second analysis, multiple 1-year series were used to assess the variability of MSL, Z0 and CD under three different CD concepts: MLWS, LLWLT and LAT. For each annual result, metrics such as the standard deviation relative to the trend and relative to the 19-year series were computed. For LAT and LLWLT, the corresponding 1-year harmonic solutions were used to generate tidal predictions for the period 2002–2021. The 1-year period was adopted because it corresponds to the minimum duration recommended by the IHO (2025) for LAT computation in oceanic areas predominantly influenced by tides. Finally, to assess the navigational risk associated with a given CD concept, the percentage of occurrences in which the observed sea level fell below the CD was calculated. An occurrence was defined as any instant at which the observed sea level was lower than the adopted CD.
2.3 Harmonic analysis methods
The harmonic analysis methods used in this study are described below. Both the temporal and conceptual analyses rely on harmonic constants and the tidal predictions derived from them.
In this study, three different approaches to harmonic analysis were evaluated, using different software packages: T-Tide, Pacmaré and UTide. The first is based on the method proposed by Godin (Reference Godin1972), later translated into MATLAB by Pawlowicz et al. (Reference Pawlowicz, Beardsley and Lentz2002), who also introduced additional improvements. The second, initially proposed by Franco and Rock (Reference Franco and Rock1972), was programmed in Fortran and features a graphical user interface through a program called Pacmaré, developed by Franco (Reference Franco2009). According to Coelho et al. (Reference Coelho, Teixeira, Silva, Alves and Daher2008), one of the main differences between T-Tide and Pacmaré lies in the domain in which the harmonic constants are obtained using the Least Squares Method: the former operates in the time domain and the latter in the frequency domain.
According to Franco (Reference Franco2009), in the Pacmare methodology, whenever a Fourier frequency approaches an astronomical frequency, a specific harmonic constituent is identified. This also applies to the calculation of nodal corrections for time series between 2 and 19 years, which are necessary when the Moon’s position relative to the lunar node is not precisely known. In the Pacmaré methodology, these corrections are applied consistently across the entire dataset, rather than only at the midpoint of the time series, as is done in T-Tide.
Despite these practical advantages, both T-Tide and Pacmaré exhibit limitations in certain cases. As noted by Codiga (Reference Codiga2011), T-Tide tends to produce unreliable results for time series between 2 and 19 years, due to inaccurate nodal corrections, imprecise inference of constituents and overly broad confidence intervals. To address these issues, Codiga proposed a new method, UTide, which introduces several improvements: exact matrix-based computation of nodal corrections and inferred constituents; support for irregular sampling, enabling the analysis of incomplete datasets; a modified Rayleigh criterion that incorporates noise levels; and improved confidence intervals using the Iteratively Reweighted Least Squares (IRLS) algorithm instead of the traditional Ordinary Least Squares (OLS).
In addition, UTide has a practical advantage over Pacmaré because it supports iterative processing and automation. For these reasons, all analyses presented in the following section, including the estimates shown in Figures 4, 5 and 8, were carried out using UTide. A comparison between UTide, Pacmaré and T-Tide is presented in Sections 3.2 and 3.2.1. To estimate the amplitudes of the harmonic constituents and Z0, all harmonic analyses were performed after removal of the mean and linear trend, following Parker (Reference Parker2007).
Statistical analysis of the evolution of harmonic constituents over a 30-year period considering non-overlapping successive time series ranging from six months to ten years. The average standard deviation of the harmonic constituents (Avg SD HC), standard deviation relative to the trend (SD trend) and mean absolute error relative to the trend (MAE trend) are presented in centimetres.
Figure 4. Long description
Panel A: A line graph shows the evolution of harmonic constituent M2 over various time periods ranging from every 6 months to every 10 years. The x-axis represents the years from 1989 to 2021, and the y-axis represents the values in centimeters. The graph includes labels for average standard deviation of harmonic constituents (Avg SD HC), standard deviation relative to the trend (SD Trend), and mean absolute error relative to the trend (MAE Trend). Panel B: A line graph displays the evolution of harmonic constituent S2 over the same time periods. The x-axis represents the years from 1989 to 2021, and the y-axis represents the values in centimeters. The graph includes similar labels as Panel A. Panel C: A line graph illustrates the evolution of harmonic constituent O1 over various time periods. The x-axis represents the years from 1989 to 2021, and the y-axis represents the values in centimeters. The graph includes labels for Avg SD HC, SD Trend, and MAE Trend. Panel D: A line graph shows the evolution of harmonic constituent K1 over the same time periods. The x-axis represents the years from 1989 to 2021, and the y-axis represents the values in centimeters. The graph includes similar labels as the previous panels. Panel E: A line graph depicts the evolution of harmonic constituent P1 over various time periods. The x-axis represents the years from 1989 to 2021, and the y-axis represents the values in centimeters. The graph includes labels for Avg SD HC, SD Trend, and MAE Trend. Each panel shows the trends and variations in the harmonic constituents over the specified time periods, with different time intervals providing insights into the short-term and long-term changes.
Statistical analysis of the evolution of MSL, Z0 and CD over a 30-year period, considering non-overlapping successive time series ranging from six months to five years. The amplitude, the standard deviation relative to the trend (SD) and the mean absolute error (MAE) are presented in centimetres. The blue line represents the multiple series; the shaded area represents the standard deviation derived from the harmonic analysis; the red dashed line depicts the trend; and the green line corresponds to the series calculated over a 19-year period.
Analysis of the trend for each series and the standard deviation relative to the trend series, the 19-year series, and the harmonic analysis. The variability of the CD is calculated as the root-sum-square of the uncertainties of MSL and Z0, following the principle of equation (2). The variability of MSL can be expressed relative to the trend (MSL trend) or to the 19-year time series (MSL 19). The uncertainty of Z0 can be expressed relative to the trend (Z0 trend), the 19-year time series (19 Z0) or to the confidence interval of the harmonic analysis (CH Z0).
Upper panel: variation in the mean propagated uncertainty of the harmonic constituents used to derive Z0. Middle panel: variation in the mean uncertainty obtained from the harmonic analyses of each group of series used to calculate Z0. Lower panel: variation in the mean PTVallc for each group of series, grouped by month. All series were analysed using both UTide and T-Tide.
Figure 7. Long description
Panel A: A line graph shows the mean uncertainty of tidal constituents over time. The x-axis represents months, ranging from 6 to 118, and the y-axis represents the mean uncertainty in centimeters. The graph includes multiple lines representing different tidal constituents, with labels such as M2, S2, O1, K1, and P1. The uncertainty values start high and decrease over time, stabilizing after around 20 months. Panel B: A line graph displays the mean Zo uncertainty with ±1 standard deviation band over time. The x-axis represents months, ranging from 6 to 118, and the y-axis represents the mean Zo uncertainty in centimeters. Two lines represent UTide and T-Tide analyses, with UTide in blue and T-Tide in red. The uncertainty values start high and decrease over time, stabilizing after around 20 months. Panel C: A line graph shows the mean PTVallc with ±1 standard deviation band over time. The x-axis represents months, ranging from 6 to 118, and the y-axis represents PTVallc in percent. Two lines represent UTide and T-Tide analyses, with UTide in blue and T-Tide in red. The PTVallc values start high and decrease over time, stabilizing after around 20 months.
Comparison of T-Tide, UTide and Pacmaré in terms of residual energy, Z0 (MLWS), and its constituent harmonics, using three selection criteria: current two-year period (2018–2019), max PTVallc (90%) and full 19-year period.
3. Results
3.1 Temporal variability
Figure 4 illustrates the variation of the harmonic constituents M2, S2, O1, K1 and P1, which compose the Z0 of the CD at Ilha Fiscal. For each time series, the following metrics were calculated: the standard deviation of the respective constituents, the standard deviation relative to the trend and the mean absolute error relative to the trend. Additionally, the value of each constituent obtained from a 19-year time series (2002–2021) is also shown. In each panel, the shaded band represents the uncertainty from the harmonic analysis of each sub-series and the units (cm) are indicated in the y-axis. In Figure 5, a similar approach is applied to the time series of MSL, Z0 and CD (CD = MSL – Z0). In these panels, the data amplitude is also presented, defined as the difference between the maximum and minimum values over the entire set of sub-series, and the date assigned to each point corresponds to the end date of the data window used in the calculation. Subsequently, all related statistical parameters are detailed in several tables, which are grouped in Figure 6 for comparison purposes. This figure highlights the variation in the trend of the harmonic constituents used in the CD calculation, as well as the estimated CD variability, considering the MSL trend; the MSL calculated over the 19-year period and three additional methods for estimating the Z0 uncertainty according to Equation 1: relative to its trend, relative to its 19-year value and based on the uncertainties derived from harmonic analyses.
Subsequently, to estimate the uncertainty of Z0 for a given time interval, Figure 7 presents the mean uncertainty obtained from all available time series using window lengths ranging from 1 month to 10 years. For each window size, the original record was divided into consecutive non-overlapping sub-series of that duration. Figure 7 also shows the variation in the uncertainties of the harmonic constituents obtained with both T-Tide and Utide, and the mean PTVallc.
3.1.1 Discussion on temporal variability
Regarding the temporal variation in the uncertainties of the harmonic constituents, and consequently of Z0, Figures 4 and 6 demonstrate a consistent decrease in the average uncertainties of all constituents as well as in the uncertainties relative to both the trend line and to the 19-year reference series. Additionally, the trend values of all series are influenced by the lengths of the sequential multiple series. For instance, in a series progressing in six-month increments, the M2 constituent exhibits a trend of +0.049 cm per year, whereas in the ten-year series, the trend is +0.034 cm per year.
This dependence on record length also affects the temporal behaviour of the MSL, Z0 and CD series shown in Figure 5. In this study, MSL is defined as the arithmetic mean of sea level calculated over a specified time window. Under this operational definition, the estimated MSL reflects the combined influence of the processes acting within the selected period, including seasonal variability, meteorological forcing, local effects and longer-term variability. Therefore, the purpose of Figure 5 is not to isolate the secular component of mean sea level, but to assess how the estimates of MSL, Z0 and consequently CD vary as a function of the adopted update window.
As the averaging window increases from six months to ten years, non-persistent and episodic effects are progressively smoothed, resulting in more stable series with smaller amplitudes, while the long-term trends remain broadly consistent. The six-month series highlights the greater sensitivity of short records to temporary fluctuations, which may lead to MSL estimates that differ substantially from the 19-year reference. This indicates that short-duration series are more likely to produce less stable CD estimates and to over- or underestimate the level effectively experienced by navigators when compared with MSL defined from a 19-year period. Longer observation periods, in contrast, provide a more robust basis for defining a long-term CD.
Uncertainties of Z0 derived from MLWS were also evaluated for multiple time series using the T-Tide and UTide programs. As shown in Figure 7, the UTide results stabilise from three years onward, with uncertainty values converging around 0.4 cm and 0.2. In contrast, T-Tide presents oscillations throughout the entire period, which may be related to inaccuracies in the nodal corrections and contamination of the confidence interval by outliers.
As for the PTVallc, computed from observed and predicted sea levels using all harmonic constituents solved in the analysis, a sharp decline is observed up to 11 months, followed by a more gradual decrease. This behaviour may reflect the growing influence of meteorological factors as the series length increases.
Regarding the method used to calculate CD variability – whether relative to the trend, the 19-year series or the harmonic analysis, a maximum difference of less than 1.7 cm was found between the methods for the 2-year series (Figure 6). Additionally, it was found that UTide tends to underestimate the uncertainties of the harmonic constituents more than T-Tide (Figure 7), when compared with the uncertainties measured relative to the trend (Figure 4).
3.2 Conceptual uncertainties
Regarding the conceptual analysis of the CD, Figure 8 presents a comparison between T-Tide, Pacmaré and UTide for three different periods used to determine Z0: the 19-year period, the period with the max PTVallc and the two-year (2018–2019) period currently adopted for defining the CD at Ilha Fiscal. To identify the period with the least noise, six-months, one-year, two-years, five-years and ten-years series were analysed, as shown in Figure 9.
Statistical analysis of the evolution of PTVallc in successive time series ranging from six months to ten years. For each panel, the maximum and average PTVallc, as well as the standard deviation relative to the trend, are present.
Figure 9. Long description
Five line graphs depict the evolution of PTVallc over different time intervals. Panel A: The line graph shows PTVallc measured every 6 months from 1989 to 2021. The y-axis represents PTVallc in percent, ranging from 70 to 90 percent. The x-axis represents the years from 1989 to 2021. The graph shows fluctuations in PTVallc with a maximum value of 89 percent, a mean value of 83 centimeters, and a standard deviation trend of 7 centimeters. Panel B: The line graph shows PTVallc measured every 1 year from 1989 to 2021. The y-axis represents PTVallc in percent, ranging from 70 to 85 percent. The x-axis represents the years from 1989 to 2021. The graph shows fluctuations in PTVallc with a maximum value of 85 percent, a mean value of 81 centimeters, and a standard deviation trend of 4.9 centimeters. Panel C: The line graph shows PTVallc measured every 2 years from 1989 to 2021. The y-axis represents PTVallc in percent, ranging from 75 to 85 percent. The x-axis represents the years from 1989 to 2021. The graph shows fluctuations in PTVallc with a maximum value of 84 percent, a mean value of 81 centimeters, and a standard deviation trend of 4.3 centimeters. Panel D: The line graph shows PTVallc measured every 5 years from 1989 to 2021. The y-axis represents PTVallc in percent, ranging from 76 to 82 percent. The x-axis represents the years from 1989 to 2021. The graph shows an increasing trend in PTVallc with a maximum value of 83 percent, a mean value of 79 centimeters, and a standard deviation trend of 2 centimeters. Panel E: The line graph shows PTVallc measured every 10 years from 1989 to 2021. The y-axis represents PTVallc in percent, ranging from 76 to 82 percent. The x-axis represents the years from 1989 to 2021. The graph shows a steady increase in PTVallc with a maximum value of 82 percent, a mean value of 78 centimeters, and a standard deviation trend of 0.57 centimeters.
To evaluate the CD (MLWS) concept adopted in Ilha Fiscal, a comparison was made with the LAT obtained over a 19-year period (2002–2021), calculated using the UTide software. From this harmonic analysis, 68 harmonic constituents were identified with a signal-to-noise ratio (SNR) greater than two. Of these 68 constituents, those whose uncertainty was smaller than their respective amplitudes were selected, totalling 65 constituents, which were used for the LAT calculation and its associated uncertainty. The harmonic constants used for the analysis of both the annual and 19-year series are provided in the Supplementary material.
Next, the number of occurrences in which sea level fell below the CD was evaluated for the LAT and four other types of CD, namely: MLWS based on a 2-year MSL and Z0, as currently adopted at the Ilha Fiscal tide gauge station; MLWS with a 2-year MSL and the Z0 with the max PTVallc; MLWS with a 19-year MSL and the Z0 with the max PTVallc; LLWLT with a 19-year Z0 and MSL; and LAT with a 19-year Z0 and MSL (Figure 10). In addition to calculating the total percentage of time during which sea level remained below the CD, these occurrences were also distributed by month and presented as histograms.
The top graphs show the percentage of time that water level (WL) remained below the indicated CD. The bottom graphs display the distribution of these occurrences by month throughout the year.
Finally, Figure 11 presents the MLWS, LLWLT and LAT datums calculated from one-year observation windows, as indicated on the x-axis. For LLWLT and LAT, the values correspond to the minimum predicted tide over the 2002–2021 period. Variability was quantified both relative to the trend and relative to the corresponding chart datums obtained from the full 2002–2021 series, thereby also providing a measure of accuracy.
Temporal variation of mean sea level (MSL), Z0 and chart datum (CD) from 1989 to 2021. Square markers denote the annual estimates. The upper panel presents the annual MSL values and the mean for the full study period. The middle panel presents the annual Z0 values for MLWS, LLWLT and LAT, together with the corresponding values derived from the 19-year reference series. The lower panel presents the annual CD values, calculated as CD = MSL – Z0, along with the constant values obtained from the 19-year reference series.
Figure 11. Long description
Panel A: A line graph titled ‘Mean Sea Level Variation’ shows the annual mean sea level (MSL) in centimeters from 1989 to 2021. The y-axis ranges from 220 to 240 centimeters, and the x-axis represents the years. The graph includes a solid line for annual MSL and a dashed line for the period mean. The MSL shows fluctuations with notable peaks around 1994 and 2019. Panel B: A line graph titled ‘Z0 Variation’ displays annual Z0 values for MLWS, LLWLT, and LAT from 1989 to 2021. The y-axis ranges from 70 to 100 centimeters. The graph includes three lines: green for MLWS, red for LLWLT, and blue for LAT, each with corresponding 19-year reference series values. The Z0 values exhibit variations with MLWS and LLWLT showing more significant fluctuations compared to LAT. Panel C: A line graph titled ‘Chart Datum = MSL - Z0’ presents annual chart datum (CD) values in centimeters from 1989 to 2021. The y-axis ranges from 120 to 180 centimeters. The graph includes lines for MLWS 1y, MLWS 19y, LLWLT 1y, LLWLT 19y, LAT 1y, and LAT 19y. The CD values show trends and variations over the years, with different lines representing various reference periods.
3.2.1 Discussion on conceptual variability
Comparing the harmonic analyses of Pacmaré, T-Tide and UTide (Table 1) revealed variations in uncertainty. Using UTide as the reference, T-Tide showed the smallest differences in the two-year series, while Pacmaré performed best in the 19-year series. Pacmaré and UTide differed by less than 1 cm for all constituents and time series, whereas differences between T-Tide and UTide for O1 and M2 exceeded 2 cm. A possible explanation is that T-Tide’s 2003 beta version included an 18.6-year harmonic analysis, which may have been affected by outliers or data gaps (Pawlowicz et al., Reference Pawlowicz, Beardsley and Lentz2002). Additionally, the negative difference for M2 and the positive difference for O1 partially offset each other in Z0, yielding an overall difference of less than 1 cm.
In Figure 10, PTVallc reached a maximum of 90% in 2015 for the six-month series. However, Z0 from this series differed by less than 1 cm from those in the two-year (83% residual) and 19-year (81% residual) series. This suggests that residual energy is not a significant criterion for defining Z0, in the MLWS concept.
Histograms show a consistent seasonal pattern for all CD types: more occurrences in autumn and fewer in spring. This behaviour is consistent with the mechanisms discussed in the introduction, including oceanic hysteresis: the ocean tends to expand during warmer months, raising MSL, and to contract during colder months, lowering MSL and increasing the likelihood of sea levels falling below CD. Tecchio et al. (Reference Tecchio, de Souza, da Silva, de Oliveira Costa, de Camargo and Harari2024) also relate lower MSL in spring and winter to persistent northeasterly winds associated with strong high-pressure systems, whereas higher MSL often follows southerly and southwesterly winds linked to post-frontal highs. At Ilha Fiscal, MSL peaks in May–June and reaches a minimum in September–October (Figure 12). This seasonal signal is consistent with the annual constituent Sa obtained from the 19-year harmonic analysis closely matching the annual variation of approximately 5 cm shown in Figure 12, as discussed by Tsimplis and Woodworth (Reference Tsimplis and Woodworth1994).
Monthly MSL at Ilha Fiscal (blue line) for the period 1989–2021. The red line represents the overall MSL for the entire period.
Slobbe et al. (Reference Slobbe, Klees, Verlaan, Dorst and Gerritsen2013) identified a similar seasonal influence of the non-astronomical component on MSL in the Netherlands, suggesting that LAT could be defined on a seasonal basis. Although this would be impractical for paper nautical charts, it becomes feasible in the context of S-10X electronic charts, where depths can be updated dynamically. For short-term tidal stations, however, a more practical solution may be to retain a fixed CD throughout the year and apply seasonal corrections only to tidal predictions, as is already done for some ports in UKHO tide tables (UKHO, 2024). Under this approach, navigational safety depends on quantifying the possible deviation of MSL from the seasonal signal.
Even in such cases, however, short observation periods may still introduce significant variability into the determination of CD. To mitigate this effect, the IHO (2011) recommends using control stations with at least one year of continuous data. As noted in Section 3.1.1, total variability can then be estimated by combining the short-term station’s Z0 variability with the long-term station’s MSL variability relative to its trend, provided that the two locations are under the same meteorological and oceanographic influence.
Figure 11 compares the temporal variation of three chart datum (CD) concepts – LAT, LLWLT and MLWS – by explicitly showing the annual evolution of mean sea level (MSL), Z0 and the resulting CD values. One possible contributor to the abrupt MSL variation is the apparent offset of approximately 10 cm between the pre-1998 and post-1998 records between the tide-staff zero and benchmark RN-2, since the historical records were reduced relative to the tide-staff zero. Because the exact date of this shift is unknown, no reliable retrospective correction could be applied without introducing an additional assumption. Excluding the pre-1999 portion of the record would also imply the loss of approximately 10 years of observations, substantially reducing the temporal coverage available to assess the influence of reference epoch on CD realisation. This limitation would be critical if the objective were to quantify long-term sea-level rise. In the present study, however, the aim is different: to assess how the realised CD depends on the reference epoch embodied in the available records. From that perspective, a possible epoch-dependent offset in the vertical reference does not invalidate the analysis; rather, it illustrates a practical mechanism through which CD realisation may vary over time and affect vertical consistency.
Additionally, as shown in Figure 3, it was observed that over the 19-year period, sea level remained below LAT for a cumulative duration equivalent to 9 days (228 h). For LLWLT, this duration was 24 days (590 h) and for MLWS, 90 days (2168 h). Alternative methodologies for calculating CD are also found in the literature. Franco (Reference Franco2009), for example, considers the probability of extreme sea level occurrences, while Slobbe et al. (Reference Slobbe, Klees, Verlaan, Dorst and Gerritsen2013) propose defining CD based on the percentage of time that the sea level remains below a given threshold. However, the results of this study indicate that the frequency of sea level falling below CD, when LAT (with observation period from 2002 to 2021) is applied in the study area, can be considered adequate for Ilha Fiscal, since a percentile value below 1% was verified.
4. Conclusions
This study shows that chart datum (CD) should be understood not only as a conceptual surface, but also as a time-dependent realisation. Over the 30-year record at Ilha Fiscal, for a six-month series, the variability of CD (MLWS), measured by the standard deviation relative to the linear regression trend, reached 13 cm (Figure 5). Moreover, in Figure 10, when comparing CD (MLWS), computed using the Z0 value associated with the highest PTVallc and a 2-year MSL, with CD (LAT) computed using a 19-year MSL, the difference in Z0 accounts for approximately 40% of the total CD difference, whereas the difference in MSL accounts for approximately 60%. These results indicate that the final value of CD is influenced not only by its conceptual definition, but also by the temporal realisation adopted.
Seasonal variability in mean sea level is pronounced (Figure 12), with the seasonal signal reaching approximately 10 cm peak-to-peak (approximately ±5 cm relative to the annual mean), while intra-month variability, represented by the error bars, is of a comparable order, with standard deviations approaching 10 cm. Lower mean sea levels occur more frequently during winter and spring. This provides a physical basis for the marked sensitivity of short-period CD realisations to the selected epoch and helps explain why an MLWS estimated for one period may, in some cases, lie below a LAT estimated for another. It should be emphasised, however, that within a given epoch, MLWS remained consistently lower than LAT.
Among the datum concepts analysed, LAT yielded the lowest frequency of sea levels falling below the datum and therefore the highest degree of navigational safety. When CD was represented by MLWS derived from a 2-year series, sea level fell below the datum on 240 days, whereas for MLWS derived from a 19-year series, this number decreased to 90 days. In this case, the difference can be attributed almost entirely to MSL, since the corresponding difference in Z0 was only of the order of millimetres. Using the same 19-year series, the number of days below the datum decreased further to 27 for LLWLT and 10 for LAT.
However, when dynamic under-keel clearance is estimated by combining real-time water-level observations with depth information, the specific definition of CD becomes less critical, because it mainly serves as the reference separating the static bathymetric component from the time-varying water-level component. Nevertheless, it remains essential that the entire navigation channel be referred to a continuous CD, without vertical discontinuities.
The uncertainty assessment based on the root-sum-square propagation of harmonic-constituent uncertainties proved consistent with the dispersion of multiple Z0 realisations relative to the trend, but its applicability is limited to Z0 defined according to the Balay method, since in this case, Z0 is directly computed as the sum of harmonic constituent amplitudes. This approach cannot be directly extended to LAT or LLWLT, because these datums are derived from 19-year tidal predictions rather than from a direct sum of amplitudes. For these CD concepts, the standard deviation relative to the trend may be adopted as an empirical descriptor of temporal variability across multiple annual realisations.
Finally, for the MLWS concept, the comparison between T-Tide, Pacmaré and UTide showed that differences in the resulting Z0 were generally smaller than 1 cm. In the context of CD realisation at Ilha Fiscal, these inter-software differences are small when compared with the contribution of MSL variability to the final CD variability. Whether a similar level of agreement also holds for LAT and LLWLT remains to be investigated, since these concepts are based on 19-year tidal predictions and involve a larger set of harmonic constituents than the direct amplitude summation adopted for MLWS. Further research is therefore recommended to assess the sensitivity of prediction-based datums to the choice of harmonic analysis method.
Overall, the results reinforce the need to report both datum concept and realisation epoch in hydrographic products and metadata. This is particularly important in areas controlled by multiple tide gauges or updated at different epochs, where inconsistent datums may introduce vertical discontinuities with consequences for hydrographic practice and navigational safety.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0373463326101647.
Competing interests
The authors declare none.
