The development of a VSEP model involves data integration from offshore and coastal areas. In this study, we used satellite altimetry to develop the VSEP model over the offshore area. Meanwhile, the VSEP model near the coastal area is usually relied on tide gauge measurements. The research framework of VSEP development is illustrated in Figure 1.

The data from multi-mission satellite altimeter are extracted and processed to obtain SSH using RADS. Subsequently, the SSH data are further processed to develop MSS and MDT and to derive LAT and HAT from satellite altimetry tidal modeling. The derivation of surface-level parameters from satellite altimetry covering offshore areas is performed before assimilation with coastal tidal data. Generally, MSS has a similar definition to MSL. For example, to calculate a tidal datum, both are denoted as the arithmetic mean of hourly tidal observations over 19 years. However, MSS is a term that describes the average of satellite-derived SSH over a period of time (Andersen and Scharroo, 2011; Yuan et al., 2020). The determination of MSS is critical in various scientific research studies, particularly in the domains of oceanography, geoscience, and environmental science.

In this study, the MSS model is determined by using a temporal average method, which depends on the subsequent processes, namely, data selection and pre-processing, ERM mean track derivation from collinear analysis, removal of GM sea-level variability, crossover minimization, and data gridding. After applying the geophysical corrections and reducing the bias in the pre-processing section, the next step is to remove the sea-level variability from ERM and GM data (Hamden et al., 2021). To minimize time variation, it is necessary to handle SSH observations in the GM of satellite altimetry in a different way. Time-averaging is inappropriate for GM data due to its lack of repeating period characteristics. The approach utilized to address the seasonal signal in GM data involves the interpolation of sea-level anomalies (SLAs), as proposed by Schaeffer et al. (2012) and Jin et al. (2016). The spatial and temporal positions of the GM data are referenced using the delayed-time Developing Use of Altimetry for Climate Studies (DUACS) Level 4 gridded SLA maps (Dibarboure et al., 2012), obtained from the European Copernicus Marine Environment Monitoring Service (CMEMS) via https://resources.marine.copernicus.eu/. Hourly gridded SLA time series corresponding to the GM data are computed to make necessary adjustments. The general processing flow in establishing the UTM20 MSS model is illustrated in Figure 2.

MDT can be defined as the separation value spanning the geoid and MSS. Different MDT models employed different geoid models in the computation. The practical approach renders the process of quantifying MDT from MSS, and a geoid model is conceptually straightforward. A precise local gravimetric geoid, the so-called MyGeoid_2017, is selected to compute the UTM20 MDT model and represent the marine geoid model in the VSEP database. In this study, a regional UTM20 MDT is determined from the difference between the MyGeoid_2017 and regional UTM20 MSS models, as derived in Equation 1 (Hamden et al., 2021): U T M 20   M D T = U T M 20   M S S − M y G e o i d _ 2017 . (1)

The UTM20 MSS model is stated in the mean tide system, while the permanent tide system of MyGeoid_2017 is not stated. According to Keysers et al. (2013), geoid undulations can be used in any system, but EGMs are generally issued as both tide-free and zero-tide systems. Generally, regional geoids acquire their tidal system from the EGM, but the system should be specified. Therefore, it is assumed that MyGeoid_2017 is provided in a tide-free system. To precisely compute MDT, both models must be in the same permanent tide. Here, UTM20 MSS is converted to a tide-free system using the conversion formula (in cm) as expressed in Equation 2 (Ekman, 1989; Keysers et al., 2013): N n = N m + 1 + k 9.9 − 29.6   sin 2 ⁡   ϕ , (2) where N n is the tide-free system, N m is the mean tide system, k is a variable called the Love number, which depends on the mass distribution within the planet (usually taken as 0.3), and ϕ is the latitude of the point. In addition, comprehensive filtering of the differences is necessary to eliminate short-scale geoid signals and to obtain a good MDT estimation. This filtering process will also eradicate the presence of residual noise in the MSS field due to unmodelled tide and the ground track striation (Farrell et al., 2012). Spatial averaging filtering methods are likely to be more decisive and precise for regional MDT applications compared with spectral filtering methods (Losch et al., 2007; Knudsen and Andersen, 2013). Spatial filters with a Gaussian-like roll-off provide more accurate results than those with sharp space cut-offs. The 2D isotropic Gaussian function is expressed as given in Equation 3 (Fisher et al., 2003; Hamden et al., 2021): G x , y = 1 2 π σ 2 e − x 2 + y 2 2 σ 2 , (3) where x is the distance from the origin in the horizontal axis, y is the distance from the origin in the vertical axis, and σ is the standard deviation of the distribution, which is also defined as the filter radius. Hence, the unfiltered regional MDT is computed by differentiating between UTM20 MSS and MyGeoid_2017, and a spatial filter is applied to smooth the MDT surfaces. The spatial filter is an average filter at which the kernel is an n x -by- n y matrix. Variables n x and n y are the number of kernel points in the east–west and north–south directions, respectively. An optimal sigma (σ) should be determined to preserve true physical signal data. A higher sigma will lead to the overfiltered signal data. Nevertheless, a lower sigma will not eliminate the error from the signal data. Here, isotropic Gaussian smoothing kernel with a standard deviation (σ) of 6 is adopted to filter the noise in the 1.5-min gridded MDT. The smoothing is performed using the imgaussfilt MATLAB function. Finally, the final regional MDT (filtered) is developed by re-interpolating it into a similar grid size, as the UTM20 MSS model. The size of the study area for the UTM20 MDT model is reduced, equal to the size of MyGeoid_2017 obtained from the DSMM, which is 0^oN ≤ latitude ≤ 9^oN and 98^oE ≤ longitude ≤121^oE.

SSH derived from satellite altimetry is utilized in developing a regional tidal datum. Computation of LAT and HAT from satellite altimetry is carried out from only four satellite missions, namely, TOPEX, Jason-1, Jason-2, and GFO. The extracted SSH data span over 24 years for the combined TOPEX class and 9 years for the GFO mission. In theory, additional missions such as ERS-1, ERS-2, Envisat, SARAL/AltiKa, and Sentinel-3A hold the potential to contribute substantial datasets to enhance ocean tide models. However, due to their sun-synchronous sampling, the observation of solar tides, particularly the principal semidiurnal S2 constituent, is severely restricted when employing the harmonic analysis method. SSH computations are pre-processed based on the best range and geophysical correction in the Malaysian region, without applying ocean tide correction. This is to prevent the elimination of the tidal signal in SSH data. In addition, SSH induced by an ocean tide signal is the signal of interest in altimetry tidal datum modeling.

Temporal offset requires considerable attention to merge TOPEX, Jason-1, and Jason-2 time series. TOPEX was launched in September 1992 and ended in August 2002. Jason-1 then continued its mission from January 2002 until January 2009. Jason-2 was then launched in July 2008 to replace Jason-1 and ended its mission in October 2016. TOPEX SSH time series are arranged ahead of Jason-1 SSH time series and followed by Jason-2 SSH time series. The end date of TOPEX SSH data is identified, and the initial data of the Jason-1 SSH time series overlapped with the near-end TOPEX SSH data are truncated. This method is also applied between Jason-1 and Jason-2 SSH time series. It is noted that only the ERM track of TOPEX class missions are used in this study for tidal datum modeling. Meanwhile, the GFO mission is a single mission with no continuity of time series from other missions. Thus, merging the SSH time series from the GFO mission is not necessary.

In theory, a satellite altimeter that travels on the repeated orbit will coincide at the same point after completing its cycle. Nevertheless, the ground tracks practically do not accurately coincide with each other for every cycle. The orbit tracks may slide within 1–2 km for each cycle. The formation of the SSH time series is performed by extracting the SSH of each cycle at discrete points. The along-track point from the first cycle of each altimetry pass is treated as the reference for the subsequent cycles. A diameter of 7 km from each point of the reference track is plotted to identify the nearest altimetry points from the reference. SSH points from subsequent cycles within the circle of 7 km radius are selected to form time series. This is because the spatial interval of 1 Hertz (Hz) satellite altimetry footprint is approximately 7 km. Thus, the altimetry points located beyond the circle are removed as outliers.

Aliasing is an effect that distinguishes or distorts other signals when a signal reconstructed from samples is different from the original continuous signal. The TOPEX class and GFO missions suffer from tidal aliasing effects due to the long temporal sampling of the satellite altimeter (9.9156 days and 17.0505 days, respectively). The aliasing period from both missions are computed as shown in Equation 4 (Lindsley, 2013; Pirooznia et al., 2016): T a = 2 π ∆ s 2 π f ∆ s − f ∆ s + 0.5 , (4) where f is the frequency of the tidal component, ∆ s is the period sample, and the bracket [.] in the formula of [ f ∆ s + 0.5 ] is the fix function that returns the greatest integer less than argument. The aliasing effect on the tidal signals as measured by satellite altimetry shows that each tidal constituent’s period seems to be longer than its actual period. Following that, the calculated aliasing period of each tidal constituent is utilized in harmonic analysis to derive amplitude and phase. The harmonic expression of the tidal height, H at time t and location φ , λ can be expressed as shown in Equation 5 (Wijaya, 2012; Ainee, 2016; Alihan et al., 2019): H φ , λ , t = Z 0   φ , λ + ∑ k A k φ , λ   cos ω k t − θ k   φ , λ , (5) where Z 0   φ , λ is the MSL of analyzed data and k represents the entire tidal constituents. A k φ , λ and θ k   φ , λ are, respectively, the tidal amplitudes and Greenwich phase lags of the tidal constituent k that is estimated via the least-square method. ω k represents the tidal frequency of k th tidal constituent. Twelve harmonic constants, namely, M2, S2, K1, O1, N2, K2, P1, Q1, MF, MM, SSA, and SA, are derived to predict the along-track tidal datums, UTM20 LAT and UTM20 HAT. To generate LAT and HAT models, the tides from both missions are predicted for a minimum duration of 19 years as proposed by Byun and Hart (2019), who acknowledge the importance of a sufficiently long tidal prediction period. Meanwhile, International Hydrographic Organization (2018) recommends using a minimum timeframe for reliable results in computing LAT and HAT.

Various approaches are used to derive the connection between a vertical datum and reference ellipsoid. Typically, most Malaysian tide gauges have no information regarding the relationship between the vertical datum and ellipsoid. Conventionally, the vertical datum is referenced either to a local benchmark (BM) or the national datum. To achieve the heights of the vertical datum with respect to the ellipsoid, UTM, in collaboration with Universiti Teknologi MARA (UiTM) Arau Campus, Perlis, has implemented Tide Gauge GNSS Campaign 2019 around Peninsular Malaysia to “tie-in” the DSMM tide gauge stations with high-precision GNSS observations. This includes the GNSS observation fieldwork carried out at Pulau Pinang, Lumut, Tg. Keling, Kukup, Tg. Sedili, Pulau Tioman, Tg. Gelang, Cendering, Geting, Pulau Langkawi, and Port Klang. Supplementary Figure S1 illustrates a total of 11 DSMM tide gauge stations involved in Tide Gauge GNSS Campaign 2019. Vertical land motion (VLM) along coastal areas has emerged as a significant concern in examining a sea-level rise across multiple decades to centuries. During this extended period, a diverse array of natural- and human-induced factors can induce VLM (Wöppelmann and Marcos, 2016). The available evidence suggests that subsidence in Peninsular Malaysia was triggered by the 2004 Mw 9.2 Sumatra–Andaman earthquake. According to Simons et al. (2023), this seismic event initiated a process that resulted in gradual sinking of the land, with a recorded decrease of 3–5 cm over the past 17 years. In East Malaysia, the rates of subsidence are comparatively lower, ranging from −0.5 to −1.0 mm per year, leading to a total landfall of 2–3 cm during the same period. However, the VLM around the Malaysian coast is not considered in this study.

The tide gauge GNSS campaign was conducted from 7 to 23 February 2019, involving 11 DSMM tide gauge stations. The GNSS observations have been performed between 5 and 6 h at each station using Trimble 5700 geodetic receiver and Zephyr Geodetic antenna, as shown in Supplementary Figure S2. Trimble Business Centre Version 5.0 (TBC v5.0) has been used to process the GNSS observations data with the cutoff angle and data interval set at 13° and 10 s, respectively. The relative positioning technique is used to resolve the GNSS data where at least 3 to 4 MyRTKnet stations have been adopted as reference stations. Due to poor sky view and multipath issues, six out of 11 tide gauges are unsuitable for direct GNSS observations on the tide gauge benchmark (TGBM). Thus, a new GNSS point, denoted as temporary benchmark (TBM), is established near the TGBM with a good sky view clearance, as shown in Supplementary Figure S2A. After GNSS observation, a precise levelling survey has been conducted using Trimble DiNi to transfer the ellipsoidal height from the new GNSS point to the TGBM. Here, the height difference obtained from precise levelling is applied to the ellipsoidal height of new GNSS point to obtain TGBM point relative to ellipsoid. Supplementary Table S4 lists the height difference of precise levelling from the new GNSS point to the six TGBMs involved.

The ellipsoidal height of MSL is determined from GNSS observations at TGBM. These data are processed in the International Terrestrial Reference Frame 2014 (ITRF2014). To achieve this, all MyRTKnet stations involved have been realized into ITRF2014 using Bernese GNSS Software Version 5.2. Then, the final coordinates of MyRTKnet stations in ITRF2014 are adopted for TBC processing to compute the GNSS solutions at TGBM. The MSL at tide gauges is derived from an observation period of 23 years from 1993 to 2015 using a simple averaging method. Here, the local MSL is relative to the zero-tide gauge. Thus, the local MSL must be transformed into MSL relative to the reference ellipsoid. Supplementary Figure S3 illustrates the relationship of various vertical surfaces to achieve tidal measurement with respect to the ellipsoidal surface (Hamden et al., 2021). The formula to obtain MSL with respect to the Geodetic Reference System 1980 (GRS80) ellipsoid is shown in Equation 6: h M S L = h G N S S − ∆ H L E V + ∆ H M S L , (6) where h M S L is the MSL height above the reference ellipsoid, h G N S S is the ellipsoidal height from GNSS observation, ∆ H L E V indicates the TGBM height above zero-tide gauge, and ∆ H M S L is the MSL height above zero-tide gauge. In determining the MDT or sea surface topography (SST), the geoid height value must be in the same reference ellipsoid and same permanent tide system with the MSL ellipsoid at tide gauges. The ellipsoid of MSL is given in a tide-free system. This has been stated by Keysers et al. (2013) where the ITRF and its realizations including the GNSS ellipsoidal height system employ the tide-free system. Apart from that, MyGeoid_2017 is given in a tide-free system and interpolated to the tide gauge locations using the bilinear interpolation method. Therefore, the ellipsoid of MSL height can be directly subtracted with MyGeoid_2017 to obtain geodetic MDT at tide gauges.

A tide gauge is a fixed point of the reference surface to observe the water level and is set up with a tide staff or other reference to continuously estimate the time series of the sea level. It has a long and high-frequency water level record (hourly data) which is beneficial for ocean investigation, including tidal analysis, prediction, and climate change. The hourly tidal data obtained via the UHSLC website (https://uhslc.soest.hawaii.edu/) are filtered using the statistical control quality method. This method implements the multiple sigma control limit, in which the multiple gross errors can be identified (Elgazooli and Ibrahim, 2012; Alihan, 2018). The sea-level time series from tide gauges are filtered using three-sigma (σ) control limits. When all measurement points are within 3σ control limits, the time series is statistically controlled, indicating that there is no gross error in the data. Consequently, the tidal data beyond the limit are excluded using the moving average technique. Then, the filtered tidal data are used to perform the tidal harmonic analysis where a tidal time series is treated as the composition of several cosine waves to estimate the presence of the tidal signal. The harmonic analysis coefficient parameters are used to estimate the amplitude and phase of principal tidal constituents.

A total of 68 tidal constituents have been estimated at each tide gauge, and these tidal constituents are subsequently used to predict the tides. The tides are predicted from 1993 to 2019, spanning more than 19 years. This is to complete the precession cycle of the moon nodes, where it takes 18.61 years to go completely around and back to their original position. Hence, LAT and HAT at tide gauges are generated from the lowest and highest predicted tides, respectively. Supplementary Figure S4 illustrates the tidal analysis and prediction, including the residuals (Top) and the 19-year tidal prediction (bottom) to determine the LAT and HAT.

A continuous vertical datum is created by combining the coastal and offshore datasets. The coastal datasets referred to the vertical reference point computed from the DSMM coastal tide gauges, while the offshore datasets referred to the vertical reference surfaces derived from satellite altimetry. MSS, MDT, LAT, and HAT are the vertical datums involved in developing the continuous VSEP model in the Malaysian region. It is noted that the integration of the vertical datum has only been conducted around the Peninsular Malaysia region, and Sabah and Sarawak are excluded due to the limitation of coastal datasets.

In this study, tidal datums derived from satellite altimetry focus only on LAT and HAT. Satellite-derived SSH time series from 24 years of TOPEX ERM mission and 9 years of GFO mission are analyzed using the harmonic approach to estimate 12 selected tidal constituents (M2, S2, N2, K2, K1, O1, P1, Q1, MF, MM, SSA, and SA). Supplementary Figure S6 shows the along-track amplitudes of the 12 major tidal constituents. Apparently, four tidal constituents (M2, S2, K1, and O1) show the most dominant tidal characteristic within the study area and have recorded higher amplitude values than the other eight constituents. These four tidal constituents are the most primary components that explain a simple model of the Earth–Sun–Moon system. The fundamental forces compelling the ocean tide signal are due to the gravitational effect induced by the moon and the sun (Lindsley, 2013).

Tidal data from seven selected tide gauges are utilized to validate the accuracy of the derived tidal constituents from satellite altimetry. The selected tide gauges include Tg. Sedili, P. Tioman, Tg. Gelang, Geting, Bintulu, Miri, and Kota Kinabalu. These tide gauges are selected based on the nearest approximation between satellite altimetry tracks and the tide gauges. In addition, the comparison between the nearest grids of DTU10 global tide models and the seven tide gauges is also conducted to determine whether the best accuracy of tidal constituents is obtained from the SDTC model or DTU10 global ocean tide model or otherwise. The development of the DTU ocean tide model is based on the finite element solution (FES 2004) and the adoption of the response method technique established by Munk and Cartwright (1966). This model uses 17 years of sea-level data from multi-mission satellite altimetry from September 1992 to December 2009. Although this study uses a similar platform as DTU10, which is satellite altimeter, it differs in terms of tidal analysis method and temporal altimetry data. This study also involves the GFO mission to derive tidal constituents. The DTU10 ocean tide model is obtained from the FTP DTU website (https://ftp.space.dtu.dk/pub/DTU10).

The nearest along-track points of SDTC to the tide gauge stations are identified and extracted. The locations and distances between tide gauge stations and the nearest point of the altimetry track are tabulated in Supplementary Table S6. In this study, two tracks of GFO altimetry and five tracks of TOPEX class are involved in the validation process. The nearest distance is at the Miri tide gauge, which is 9.14 km, and the longest distance between both points is located at the Bintulu tide gauge, which is 66.61 km. It has perceived that the maximum range between the two points is closer than that in the previous studies by Daher et al. (2015) and Fu et al. (2020), where the maximum ranges from both studies are 124 km and 77 km, respectively. The RSS values of eight tidal constituents between two models (SDTC and DTU10) and each tide gauge station are computed and presented in Figure 5. The figure depicts that SDTC provide promising RSS results compared to DTU10 in all tide gauges involved. There are significant differences in RSS values between SDTC and DTU10 located at Tg. Sedili, Geting, and Bintulu tide gauges. The differences might be due to the interpolation error of the DTU10 grid model to the tide gauge station and the insufficient time span of altimetry time series in estimating tidal constituents.

Based on the previous study by Fu et al. (2020), the RSS values for along-track satellite altimetry points were set less than 8 cm in the deep ocean (depth >200 m). Thus, this study follows the criterion of 8 cm to determine whether the SDTC are well estimated or not. Two out of seven tide gauge stations, which are P. Tioman and Geting, have passed the criterion and have high-precision outcomes with RSS values of 5.3 cm and 5.0 cm, respectively. The Bintulu station has the highest RSS value, which is 35.5 cm. This is expected due to the longest distance between the altimetry point and the tide gauge stations. Other stations have recorded the RSS values within 10 cm. The outcomes are considered good even though the RSS values do not achieve the 8 cm criterion. For instance, the nearest SDTC to the Tg. Sedili and Tg. Gelang are from the GFO mission, where the time series data for tidal constituent estimation are limited. In addition, the depth in the study area is mainly less than 200 m. Therefore, it is expected that the RSS values could exceed 8 cm. It can be inferred that high precision of satellite-derived tidal constituents has been achieved with the TOPEX class satellite mission. This is because the satellite altimetry time series from TOPEX class spans over 24 years where all the constituents have met the criterion for the 19-year time series, and all tidal constituents can be appropriately estimated. Therefore, the precision of tidal constituents estimated from satellite altimetry data is relatively acceptable in most shallow water locations.

Tides in this study are predicted based on the estimated tidal constituents. Supplementary Figure S7 illustrates the predicted time series at each along-track altimetry point. The study area in this figure is reduced to follow the size of the study area from the UTM20 MDT model. This is because each VSEP model produced must be in the same size as the study area before it is archived in the database. The black triangles in Supplementary Figure S7 indicate the coastal DSMM tide gauge stations used to assess the estimated offshore tidal datum. Each point from the along-track predicted SSH time series has been randomly selected (labelled as red dots in Supplementary Figure S7) at each of the Malaysian seas. The selection is made to visualize the modelled time series. It is noted that the point must be selected in the offshore area.

Supplementary Figures S8A–D show the time series of selected modelled SSH and its residuals from the TOPEX class on the left column and the GFO mission on the right column. The observed and predicted SSH time series are plotted in blue and green lines, respectively. These residuals are derived by computing the differences between the predicted and observed SSH time series. Subsequently, the quality of the predicted SSH can be determined based on these residuals using RMSE computation. Apparently, the SSH time series data from the TOPEX class are denser than the GFO mission. This is because the repeat period of the TOPEX class is shorter, which is 9.9156 days, compared to GFO, which is 17.0505 days. The highest RMSE value between observed and predicted SSH from TOPEX is recorded at Malacca Strait, followed by the Celebes Sea, the Sulu Sea, and the South China Sea, which is 10.1 cm, 9.8 cm, 7.6 cm, and 7.3 cm, respectively. This might be because Malacca Strait is a closed sea with tidal characteristics that is likely to have a large gradient. Nevertheless, the highest RMSE value from the GFO mission is recorded at the Celebes Sea, followed by Malacca Strait, the Sulu Sea, and the South China Sea, with values 10.9 cm, 7.9 cm, 7.0 cm, and 6.5 cm, respectively. Therefore, it can be inferred that the precision of tidal prediction from both missions are reasonable at the offshore area since the RMSE values are within the range 6.5–10.9 cm.

LAT is defined as the lowest water level, which can be predicted to occur under average meteorological conditions. Meanwhile, HAT is defined as the highest water level, which normally occurs when any astronomical conditions are combined. Both tidal datums can be derived by analyzing several years of tidal data or tidal predictions, usually 18.6 years to account for a complete nodal cycle. In this study, the SSH time series from the TOPEX class and GFO missions are predicted for at least 19 years or more. The predicted time series from the TOPEX class is from 1993 until 2019, while from GFO, it is from 2000 until 2009. The lowest and highest time series of predicted tides indicate the LAT and HAT, respectively. Subsequently, the derived LAT and HAT from the along-track satellite altimetry are interpolated and validated against the selected coastal DSMM tide gauges as distributed in Supplementary Figure S7. The derived tidal datums are converted relative to MSL prior to the results to allow comparison with the tide gauges.

The combination of TOPEX class and GFO along-track missions (LAT_MSL and HAT_MSL) is interpolated into a regular grid of 0.125° using ordinary kriging and minimum curvature spline methods, as shown in Supplementary Figure S9. This figure clearly indicates that the middle part of the Malacca Strait has the greatest tidal range of LAT_MSL and HAT_MSL compared to other regions, with values up to −2.6 m and 3.0 m, respectively. The tidal range of LAT_MSL and HAT_MSL in the middle of the South China Sea, Sulu Sea, and Celebes Sea depicts the values within −0.9 m to −1.1 m for LAT_MSL and within 0.9–1.2 m for HAT_MSL. The greatest tidal range of LAT_MSL and HAT_MSL can also be observed at the southwest of East Malaysia in the South China Sea, with values up to −2.1 m and 2.2 m, proportionately. Meanwhile, near the Gulf of Thailand, the northwestern part of the South China Sea depicts the lowest tidal range of LAT_MSL and HAT_MSL, with values near zero meters. Apparently, the tidal ranges mainly are larger near the coastal areas than the offshore areas. It is visible at the coastal part of the Celebes Sea, as shown in Supplementary Figure S9. The reason is that the satellite altimetry data acquired near the coastlines are contaminated by the inclusion of land in the footprint signal or the fact that the tide is on the ebb (Fok, 2012). Moreover, the tidal datum models interpolated using the ordinary kriging method generate smoother contour lines than the models interpolated using the minimum curvature spline method. These gridded offshore tidal datum models are then validated with the selected coastal tide gauges to determine the best interpolation method.

Generally, the optimal method to validate satellite altimetry-derived LAT and HAT is demonstrated by comparing them with the offshore tide gauges. However, the lack of establishment of offshore tide gauges and the difficulty in obtaining offshore tidal data from offshore authorities hinder this validation method. Thus, this paper adopts the statistical assessment by validating the UTM20 LAT and HAT models with 10 selected DSMM tide gauges, as shown in Supplementary Figure S7. The assessment results of LAT_MSL and HAT_MSL are described in Supplementary Table S7 and Supplementary Table S8, respectively. For UTM20 LAT, the RMSE values obtained between the minimum curvature spline model and in situ data are lower than those of the ordinary kriging model, which recorded 25.5 cm and 31.8 cm, respectively. Meanwhile, the RMSE values of UTM20 HAT between minimum curvature spline and in situ data are also lower than those of the ordinary kriging model, which yield 17.4 cm and 33.8 cm, respectively. Thus, it can be inferred that the offshore tidal datum models generated using the minimum curvature spline method have better agreement with coastal tide gauges compared to ordinary kriging models, despite the ordinary kriging producing smooth contour surfaces.

Four types of vertical datums computed at coastal tide gauges are MSL, MDT, LAT, and HAT. These datums are then integrated with the vertical reference surfaces derived from satellite altimetry data. Only 11 coastal tide gauges along Peninsular Malaysia are used to compute the vertical datum in this research, since Tide Gauge GNSS Campaign 2019 could only cover the areas around Peninsular Malaysia. MSL values are obtained at each tide gauge by conventionally averaging the hourly tidal data. It is noted that the averaged values are relative to the zero-tide gauge. Thus, these values must be shifted to the reference ellipsoid by applying Equation 6. Subsequently, MDT values are computed by subtracting the derived MSL values with the interpolated MyGeoid_2017 at tide gauge stations. The bilinear interpolation method is adopted to interpolate the MyGeoid_2017 to the tide gauge stations.

Supplementary Table S9 tabulates the GNSS processing results of Tide Gauge GNSS Campaign 2019 using TBC v5.0. The stations marked as TBM1 until TBM6 indicate that the GNSS observation are not directly conducted on the TGBM due to poor sky-view conditions and might be induced by multipath signals. Instead, the GNSS observations are conducted on the temporarily established points to obtain ellipsoidal height. Then, the final ellipsoidal height is transferred to the TGBM through precise levelling, as described in Section 2.4. The findings show that the ellipsoidal height at P. Langkawi, P. Pinang, Lumut, and Geting are located below the ellipsoidal surface and thus the negative values. Good accuracies are displayed on the latitude and longitude of GNSS results, where the standard deviation is recorded less than 1 cm. Furthermore, the standard deviation of ellipsoidal height ranges from 0.7 cm to 4.8 cm.

Table 3 shows the computation of MSL and MDT at 11 DSMM coastal tide gauges. The MSL values relative to zero-tide gauge (hereinafter ∆H_MSL) are computed for 23 years from 1993 to 2015. It is found that the MSL values vary at each location, which range from 2.233 m to 4.030 m along Peninsular Malaysia. Here, the separation offset between the computed MSL and MyGeoid_2017 is computed to obtain MDT values, as shown in Supplementary Figure S10. The results show that P. Tioman displays the largest separation offset between MSL and MyGeoid_2017, which recorded 1.053 m, while P. Langkawi records the lowest separation offset with 0.651 m. It also can be inferred that the slope of MDT values increases from north to the south of Peninsular Malaysia at both west and east coast areas. In general, along the west coast, the separation offset between UTM20 MSS and the MyGeoid_2017 model shows a progressive latitude-dependent trend consistent with the findings from Mohamed (2003) and Pa’suya (2020). The computation of MDT at coastal tide gauges is very significant in this research as it can facilitate the integration between the vertical reference surface at coastal and offshore areas to obtain a smooth continuous model of MSS. This MSS model is established by recombining the integrated MDT with the local gravimetric geoid model, MyGeoid_2017.

For the computation of LAT and HAT, the filtered hourly tidal data are analyzed and predicted using the harmonic analysis approach. The computed LAT and HAT values at each tide gauge, including the RMSE of residuals, are tabulated in Supplementary Table S10. The RMSE of residuals is the accuracy of the predicted sea level, which is calculated between the predicted and observed sea levels using Equation 11. The results show that the accuracy of the predicted sea level at each tide gauge is between 9.5 cm and 13.8 cm, where Tg. Keling records the highest accuracy and Port Klang displays the lowest accuracy. The computation of LAT and HAT relative to the reference ellipsoid is based on the modification of Equation 3.9, in which the parameter of ∆H_MSL is replaced by ∆H_LAT and ∆H_HAT. Here, the findings show that the LAT values from tide gauges are within the range −17.778 m–7.675 m. Meanwhile, the HAT values are within the range −14.246 m–10.788 m.

The establishment of UTM20 MSS, MDT, LAT, and HAT is successfully conducted to represent the offshore datasets of the vertical datum. In addition, these four vertical datums are also well developed at each tide gauge station involved in this research to represent the coastal datasets. A continuous VSEP model is developed by combining offshore and coastal surfaces. The ellipsoid-based transformation approach utilizes a set of gridded surfaces, in which each surface defines the separation of one vertical datum from the WGS84 ellipsoid. The four gridded vertical separation surfaces are MSL-WGS84, MyGeoid_2017-WGS84, LAT-WGS84, and HAT-WGS84.

The tide gauge-derived MDT is interpolated to a surface extending from the coastline to 10 km offshore. Three interpolation techniques, namely, IDW, ordinary kriging, and thin-plate spline, have been adopted to determine the optimal interpolation for coastal datasets. Table 4 shows the statistical differences between the actual and predicted values of tide gauge MDT using three interpolation techniques. The spline technique displays the optimal accuracy of MDT surfaces, followed by IDW and Ordinary Kriging, which recorded 0.2 cm, 0.48 cm, and 0.52 cm, respectively. Evidently, all interpolation techniques yield better than 1 cm accuracy when compared to the input tidal data.

For further investigations on each interpolation method, one tide gauge at a time is removed, and the MDT surface is then re-interpolated. The actual value for the removed tide gauge is compared to its predicted values. Supplementary Table S11 shows the statistical results of the differences between the actual and predicted values using the removal test for tide gauge-derived MDT from three interpolation techniques. The findings show that the IDW and ordinary kriging techniques are rejected based on their higher mean than the spline method. The standard deviation of the spline method is substantially better than the standard deviation of IDW and ordinary kriging. Apart from that, spline also generates the smoothest surface and fits well with satellite altimetry data, as shown in Figure 6. Therefore, the spline technique is well suited to create gradually changing surfaces, such as elevation and sea-level heights. The UKHO also uses it to represent LAT relative to MSL (Turner et al., 2010).

Here, the interpolation techniques disregard islands or any form of shorelines; instead, the data are interpolated assuming all surfaces are ocean. Incorporating a technique that determines the correlation of tidal data by distance over sea and utilizing a coastline polygon to account for islands’ impacts and bending shorelines may lead to future advancements in this research. Enhancement by adopting the suggested approach would be ideal with the expansion of tide gauge distribution to describe MSL variation along the coast better. Furthermore, establishing a coastal thread or hydrodynamic stream, as recommended by Keysers et al. (2013), could be adopted for future improvements. After creating the MDT surfaces extending 10 km from the coastline, these coastal datasets are integrated with 20 km offshore of UTM20 MDT by adopting ordinary kriging and minimum curvature spline methods. The statistical results of the difference between in situ and two interpolation MDT models are shown in Table 5. The findings present that the minimum curvature technique records slightly better accuracy than ordinary kriging with the standard deviation values of 1.8 cm and 1.9 cm, respectively. Both interpolation models yield the mean difference below 1 cm.

These two interpolation MDT models are then combined with the MyGeoid_2017 model to develop a smooth continuous MSS model. To determine the best interpolation techniques in establishing the final integrated MSS, validations of these surfaces are carried out by comparing them with the input tide gauge values and the original altimetric UTM20 MSS. This comparison analysis is solely demonstrated to determine the accuracy of the interpolation procedures since no other known combined model spans within the study area for comparison. Table 6 displays the results of the comparison analysis. The integrated MSS models from ordinary kriging and minimum curvature methods have lower mean difference with in situ data of −0.003 m and −0.001 m, respectively, than the UTM20 MSS, which recorded −0.079 m. In addition, the accuracy of integrated MSS models also displays better agreement than pure altimetric MSS. In addition, ordinary kriging yields a slightly higher accuracy of 0.013 m than the minimum curvature spline technique with an accuracy of 0.014 m. Thus, the surface resulting from the ordinary kriging interpolation is selected as the final integrated MSS known as iUTM20. This demonstrates that the final integrated MSS closely matches the tide gauge values compared to UTM20 MSS, implying more accurate MSS in theory.

It is challenging to distinguish the difference between the integrated MSS and altimetric MSS patterns, as shown in Supplementary Figure S11. However, the integrated MSS has a better resolution and a distinct height pattern due to the assimilated tidal gauge data. The notable difference can be further examined through the statistical analysis shown in Table 7. Significant improvement is displayed from the integrated MSS model, where most of the differences at each tide gauge are below 10 cm compared to altimetric UTM20 MSS. In areas with no tide gauges, the integrated MSS cannot depend on its current form. Thus, increasing the density of tide gauge distribution is recommended to enhance the accuracy of the integrated MSS in this region. The final integrated MSS model, iUTM20 MSS, is illustrated in Figure 7. The integration can be enhanced when better and more tide gauge data became available, as it would be possible to analyze the spatial behavior of MSL between tide gauge data and satellite altimetry. Hence, it could improve the quality of the integration. Techniques adopted by other projects can be used or serve as a guide for Malaysia’s future techniques. For instance, VORF has adopted a combination of least-square collocation and specialized algorithms based on coastal topography to interpolate the reference surfaces between tide gauge and altimetry data (Iliffe et al., 2007). Other than that, the GNSS buoy surveys can be performed in the gap between tide gauge and altimetry data to estimate MSL relative to the ellipsoid and interpolate the datasets using the least-squares technique. This method is adopted by French BATHYELLI (Pineau-Guillou and Dorst, 2013).

Integrated UTM20 LAT and HAT surfaces are the combination of the tide gauge data and altimetric-derived tidal modelling. Since the integration between tide gauges and altimetric-derived tidal datums is conducted along Peninsular Malaysia, the separation values in the other regions only depend on altimetric data. The results of LAT and HAT surfaces are shown in Figure 8. The figure clearly shows that the middle part of Malacca Strait has the lowest and the highest values of LAT and HAT, respectively. This indicates that this area has the optimal tidal range of sea level compared to other regions. In addition, there are closer contour lines at Malacca Strait, which indicate possible drastic changes in LAT and HAT values. This is probably influenced by the shape of Malacca Straits, which is narrow closed sea and shallow water.

A gridded point of 1.5 min, with approximately 2.7 km spacing, is created over the study area to establish the ellipsoidal tidal datum separation surfaces. The space of the gridded point should be denser near the coastal area to capture the coastal variations. Unfortunately, this study did not adopt the hydrodynamic model along the coastal regions due to the intensive computational nature of estimating the model, and it requires a high-performance computing system. Each gridded tidal datum surface is added independently to the final integrated MSS, iUTM20 MSS, to generate gridded ellipsoidal tidal datum separation surfaces. The final statistics of tidal datum separations are tabulated in Table 7. The statistics indicate that the RMS agreements between the established LAT and HAT relative to MSL with coastal tide gauges are 1.8 cm and 2.0 cm, respectively. For the ellipsoidal LAT and HAT, the RMS agreements yield 1.8 cm and 2.5 cm, respectively. The findings also show that the iUTM20 tidal models have significant improvement compared to the difference of LAT_MSL and HAT_MSL, which only derived using altimetric data, as shown in Supplementary Table S8 and Supplementary Table S9, respectively. Since all the available data have been utilized to compute the surface, it has noted that the results are merely a measure of precision.