Tide–surge interaction observed at Singapore and the east coast of Peninsular Malaysia using a semi-empirical model

Tide–surge interaction plays a substantial role in determining the characteristics of coastal water levels over shallow regions. We study the tide–surge interaction observed at seven tide gauges along Singapore and the east coast of Peninsular Malaysia, focusing on the timing of extreme non-tidal residuals relative to tidal high water. We propose a modified statistical framework using a no-tide–surge interaction (no-TSI) null distribution that accounts for asymmetry and variation in the duration of tidal cycles. We find that our modified framework can mitigate false-positive signals of tide–surge interaction in this region. We find evidence of tide–surge interaction at all seven locations, with characteristics varying smoothly along the coastline: the highest non-tidal residuals are found to occur most frequently before tidal high water in the south, both before and after tidal high water in the central region, and after tidal high water in the north. We also propose a semi-empirical model to investigate the effects of tidal-phase alteration, which is one mechanism of tide–surge interaction. Results of our semi-empirical model reveal that tidal-phase alteration caused by storm surges is substantial enough to generate significant change in the timing of extreme non-tidal residuals. To mitigate the effect of tidal-phase alteration on return level estimation, skew surge can be used. We conclude that (1) tide–surge interaction influences coastal water levels in this region, (2) our semi-empirical model provides insight into the mechanism of tidal-phase alteration, and (3) our no-TSI distribution should be used for similar studies globally.

1 Introduction

Coastal regions are vulnerable to the combined effects of tides and storm surges, which can induce significant sea-level variations and pose substantial risks to coastal communities and ecosystems (Diaz, 2016; von Storch et al., 2015; Hinkel et al., 2014). Storm surges are rises in sea level brought about by low atmospheric pressure and strong winds acting on the sea surface (Pugh and Woodworth, 2014a). Large storm tides occur when large storm surges coincide with high tides (Stephens et al., 2020; Gregory et al., 2019), which increase the risk of coastal flooding and threaten lives and livelihoods (Pugh and Woodworth, 2014c). The likelihood and impact of such destructive events are further aggravated by sea-level rise (Calafat et al., 2022; Marcos and Woodworth, 2017; Woodworth and Blackman, 2004). To identify the appropriate response to such extreme water level events, we must understand the fundamental processes of tides and storm surges and their mutual interactions.

A dependence between tides and surges has long been noticed at coastal locations (Keers, 1968) and recognised to be caused by interaction between tides and surges (Pugh and Vassie, 1978; Wolf, 1978). Understanding tidal dynamics, surge generation, and their mutual interaction is required to improve operational forecasts of sea levels and the statistical estimation of their extremes (Olbert et al., 2013; Horsburgh and Wilson, 2007; Tawn and Vassie, 1989). This tide–surge interaction is non-linear and can lead to complex coastal dynamics characterised by amplification or attenuation of water levels, which are influenced by local bathymetry, coastline geometry, and atmospheric conditions (Idier et al., 2019). Strong tide–surge interaction has been observed in shallow waters and estuaries (Wolf, 1981) and hence is expected at the Sunda Shelf where Singapore and Malaysia are located. The main mechanism behind tide–surge interaction is mutual phase alteration (Rossiter, 1961). The generation of surges over a water body is influenced by the depth of the water body. As changes in depth occur partially due to tides, the height of surges can be influenced by tides. The propagation speed of tides is also dependent on depth, which can change due to surges (Proudman, 1957, 1955). Further studies that analysed shallow-water equations along the coast found that the non-linear tide–surge interaction is caused by the advection term, non-linear bottom-friction term, and the non-linear shallow-water effects of the shallow-water equations (Zhang et al., 2021; Idier et al., 2012).

Studies of tide–surge interaction have employed a range of modelling approaches, including statistical methods (Arns et al., 2020; Haigh et al., 2010; Dixon and Tawn, 1994), numerical models (Costa et al., 2023; Horsburgh and Wilson, 2007; Prandle and Wolf, 1978), and analytical models (Horsburgh and Wilson, 2007). Dixon and Tawn (1994) proposed a statistical framework where they split the tidal range into five equiprobable bands and used a chi-squared test to determine whether non-tidal residuals above a height threshold fall uniformly into each band. Horsburgh and Wilson (2007) proposed a modified version of the framework where the tide is instead split into 13 hourly bands between 6.5 h before and 6.5 h after tidal high water (HW). Horsburgh and Wilson (2007) also provided a simple mathematical explanation for the abundance of large non-tidal residuals at timings halfway between tidal low water and HW.

Application of such frameworks revealed that extreme residuals are most often found 3–5 h before HW in the Bay of Bengal (Antony and Unnikrishnan, 2013) and in the North Sea (Horsburgh and Wilson, 2007) and typically about 2 h before HW in the English Channel (Haigh et al., 2010). In China and New Zealand, observed tide–surge interaction varies along the coastline: the frequency of extreme residuals peaks before HW at certain locations, after HW at others and is independent of tides at the remaining locations (Costa et al., 2023; Feng et al., 2019). Numerical models have shown that the inclusion of tide–surge interaction often results in better water level predictions, especially over coastal and shelf waters, whereas the omission of the interaction may lead to under- or overestimation of surges at certain locations (Fernández-Montblanc et al., 2019; Idier et al., 2012). For example, Antony et al. (2020) showed that numerically modelled peak water levels generated during Cyclone Aila at the head of the Bay of Bengal would have been overestimated if tide–surge interaction was not simulated.

Here, we focus on investigating the tide–surge interaction observed at seven tide gauge locations near Singapore and the east coast of Peninsular Malaysia using modified statistical methods and a new semi-empirical model. Our research objectives are (1) examining the tide gauge records to characterise the tide–surge interaction observed at each tide gauge location and the spatial pattern across locations and (2) explaining the observed interaction characteristics through a simple semi-empirical model.

We apply a modified version of the statistical method by Horsburgh and Wilson (2007) to determine the presence of tide–surge interaction at each tide gauge location and to characterise these interactions. The existing method groups residuals above a height threshold into 13 hourly bands between 6.5 h before and 6.5 h after HW, counts the number of residuals in each band, and compares the resulting distribution to the uniform distribution using a chi-squared test. However, due to asymmetry and variation in the duration of tidal cycles (Guo et al., 2019), the expected null distribution is not a uniform distribution over 13 hourly bands centred at HW especially at locations with mixed or diurnal tides lasting up to 25 h. Hence, our key modification to the existing methodology is to replace the uniform distribution with a “no-TSI distribution” as the null distribution. We also use an exact statistical test for hypothesis testing instead of the chi-squared test. In addition, we propose a simple approach to compare the tide–surge interaction across locations that have different tidal characteristics (diurnal, mixed, or semidiurnal).

We aim to provide an explicit explanation of the observed tide–surge interaction through our semi-empirical model by combining historical tide and surge data with winds, coastal geometry, and bathymetry. The semi-empirical model accounts for the mechanism of tidal-phase alteration: storm surges perturb the depth of the water body, which influences the propagation speed of the tide. This results in differences between observed tides and tides predicted from harmonic analysis, which are detected as non-tidal residuals. We use our model to show that this mechanism can significantly influence the timing of extreme residuals to produce signals of tide–surge interaction.

2 Singapore and the east coast of Peninsular Malaysia

The seven tide gauges are located within the Sunda Shelf, illustrated in Fig. 1a, which has a typical depth of 40–80 m. The shallowness of this shelf likely enhances the interaction between tide and surge, and its expanse allows for phase changes in the tides to compound as the tide propagates across the shelf. The eastern coastline of Peninsular Malaysia faces the South China Sea and is exposed to strong seasonal monsoon winds resulting in larger surges (Mohd Anuar et al., 2020). To the north of our study region lies the Gulf of Thailand. To the south of Singapore, the southernmost location of our study region, lie a part of the Riau Islands, Sumatra, the Karimata Strait, and the Java Sea. The Malacca Strait is west of Singapore and leads to the Indian Ocean. The diurnal and semidiurnal tides that propagate from the South China Sea and the Indian Ocean, respectively, drive the complex mixed tides in Singapore and the southern parts of the study region (van Maren and Gerritsen, 2012). At Singapore, the diurnal tides that propagate from the South China Sea are further amplified from reflection against the east coast of Sumatra (van Maren and Gerritsen, 2012).

Figure 1 (a) Bathymetry of the region of interest of this study obtained from GEBCO, the General Bathymetric Chart of the Oceans (GEBCO Compilation Group, 2023). The inset shows the location of this region within East and Southeast Asia. The seven tide gauge stations analysed are marked in red circles. The red rectangle denotes the region where 10 m winds and mean sea-level pressure are considered when calculating R_surge (Sect. 3.4). The region is roughly a rectangle of 759 km by 833 km. The unit vector $\widehat{\mathit{x}}$ used in Eqs. (5)–(7) is shown in the figure. (b–h) Close-ups of the seven tide gauge locations (OpenStreetMap contributors, 2017). © OpenStreetMap contributors 2017. Distributed under the Open Data Commons Open Database License (ODbL) v1.0.

A close-up of the bathymetry around the seven stations is shown in Fig. S1, and detailed maps of their immediate vicinity are shown in Figs. 1b–h and S2. Of our seven tide gauges, Tanjong Pagar and Johor Baharu are located within straits; Sedili at a river mouth; Tioman between Tioman Island and mainland Peninsular Malaysia, sheltered from the South China Sea; Kuantan at the mouth of a man-made bay; Cendering on the open coast; and Geting sheltered inside a bay.

Monsoonal winds are the main determining factor of hourly surges and their extremes in the Singapore Strait (Tkalich et al., 2009, 2013a). Storm surges occur along the coast throughout the year, reaching 0.35–1.0 m at Geting, the northernmost location of our study region (Abd Razak et al., 2024). The northern part of the coast is more prone to stronger storm surges due to tropical depressions (Mohd Anuar et al., 2023). On the southern part of the coast, a surge of about 1.0 m caused by Typhoon Vamei in 2001 was recorded (Mohd Anuar et al., 2018). Several storm surge events in Singapore reaching 0.3–0.8 m have also been documented since 1974 (Luu et al., 2016), particularly during the winter monsoon period.

Our analysis found that the tidal range in this region can reach 2.7–3.6 m and that the largest non-tidal residuals over the past 30 years exceeded 0.8 m. Understanding how these components interact and combine provides insight into the contributors to coastal water levels. Hydrodynamical processes, mainly those caused by regional wind forcing during the seasonal monsoons, have shown a strong influence on the water levels at the seven tide gauges analysed in this study (Tay et al., 2016; Luu et al., 2016; Kurniawan et al., 2015; Karri et al., 2014; Tkalich et al., 2013a; Chen et al., 2012).

3 Data and methods

The general methodology begins with separating hourly tide gauge data into tides and non-tidal residuals. Following the methodology of Horsburgh and Wilson (2007), we identify the largest non-tidal residuals, find their timings relative to their nearest tidal high water, and tabulate them as a frequency distribution. Next, we calculate the no-TSI distribution. We compare the frequency distribution to the no-TSI distribution using hypothesis testing to determine the presence of tide–surge interaction. We then generate model residuals using our semi-empirical model and repeat the above procedure to determine if tide–surge interactions are generated by our model. We use research-quality hourly tide gauge data obtained from the University of Hawaii Sea Level Center. Bathymetry, hourly 10 m winds, and mean sea-level pressure are obtained from ERA5.

3.1 Tide gauge data, tides, and residuals

Research-quality hourly water level from tide gauges at seven locations (Fig. 1) along Singapore and the east coast of Peninsular Malaysia are obtained from the University of Hawaii Sea Level Center (Caldwell et al., 2001). Details of the tide gauge records are tabulated in Table 1. Observations have been made over at least 29 years at each location, with a data completion rate of 95 %–99 %. The length of these records is close to the 30-year threshold typically considered the minimum length required for analysis of extreme sea levels (Rasmussen et al., 2018). Nonetheless, these stations have some of the longest records within Southeast Asia (Caldwell et al., 2001).

To compute the tidal level and non-tidal residuals, we use the equation

$$\begin{array}{}\text{(1)}& X\left(t\right)=Z\left(t\right)+T\left(t\right)+R\left(t\right).\end{array}$$

X(t) is the hourly water level recorded by a tide gauge at time t. Z(t) is the 1-year (8766 h after accounting for leap years) moving average of X(t). A 1-year window was chosen because intra-annual behaviour of residuals is well-understood to be periodic due to seasonal variations, and its influence on residuals is of interest in this study (Tkalich et al., 2013b). Detrending the hourly water level by the 1-year moving average removes trends of annual and longer timescales from the tide gauge records. The 1-year moving average is calculated only if at least 85 % of the data in the 8766 h window is available. This reduces our usable data to 87 %–99 % across the seven locations. The standard deviation and certain key quantiles of the detrended water level, X(t)−Z(t), are tabulated in Table 2. While water levels of this height do not pose an immediate threat to the region, understanding the underlying drivers is integral for applications such as early detection of extreme events, forecasting, or coastal protection design purposes, especially with sea-level rise and the increasing frequency and intensity of storms (Mohd Anuar et al., 2018). T(t) is the tidal level that we estimate using UTide, a tidal harmonic-analysis package implemented in Matlab (Codiga, 2011). X(t)−Z(t) is split into two halves based on the start and end dates in Table 1 and used as inputs to UTide, as UTide does not recommend using records longer than 18.6 years as input. UTide was used to identify the amplitudes and phases of 66 tidal constituents with periods of up to 32 d through harmonic analysis. Constituents with a signal-to-noise ratio of at least 2 are used to construct T(t). R(t), the residual, is then estimated as $R\left(t\right)=X\left(t\right)-Z\left(t\right)-T\left(t\right)$. We denote the residual obtained through this procedure as R_gauge. The time series of X(t)−Z(t), T(t) and R_gauge at each tide gauge during their largest-recorded R_gauge are illustrated in Fig. S3.

A summary of the tidal characteristics at the seven locations is tabulated in Table 3. The tides at these locations lie within the microtidal (< 2 m) and mesotidal (2–4 m) ranges, with diurnal tidal ranges of close to 2 m, except at Cendering and Geting where diurnal tidal ranges are 1.5 and 0.8 m, respectively. We found that the daily tidal range can be as large as 3.6 m at Johor Baharu and Kuantan, similar to other studies (Ismail et al., 2018; Ng and Sivasothi, 1999). On a shelf with a depth of 40–80 m, tides of 3.6 m produce a water depth deviation from its mean by up to about ± 4 %. The tides at a location can be categorised into diurnal, semidiurnal, or mixed using the tidal form factor F, which compares the relative importance of the following diurnal and semidiurnal tidal constituents (Pugh and Woodworth, 2014b):

$$\begin{array}{}\text{(2)}& F={\displaystyle \frac{K_{\mathrm{1}}+O_{\mathrm{1}}}{M_{\mathrm{2}}+S_{\mathrm{2}}}}.\end{array}$$

A common classification considers a location with a factor of < 0.25 to be semidiurnal, 0.25–1.50 to be mixed with mainly semidiurnal tides, 1.50–3.00 to be mixed with mainly diurnal tides, and > 3.00 to be diurnal (Pugh and Woodworth, 2014b). Based on this classification, the tides at all stations are identified as mixed, with Cendering and Geting mainly diurnal, while all other stations to the south are mainly semidiurnal. Out of all tidal constituents estimated across the seven locations, we find that the amplitude of most constituents does not change by more than 5 mm between the first and second half of their tide gauge records while eight constituents saw a change of 5–13 mm. The change in mean diurnal tidal range is about 0.2 cm at Tanjong Pagar to 4.9 cm at Geting. Attribution of these changes remains a difficult task and is outside the scope of this study, but some local drivers include (1) dissipation and turbulent mixing; (2) depth of channels and flats; (3) surface area, width, and convergence; (4) resonance and reflection; (5) river flow; and (6) changes in instrumentation (Haigh et al., 2020).

Table 1 Data availability of the tide gauges used in this study. The completion rate is the percentage of usable hourly observations out of the duration of records. The usable rate is the percentage of usable observations subtracting 1-year moving averages (Sect. 3.1).

Table 2 Standard deviation and certain key quantiles of detrended water level, X(t)−Z(t), at the seven tide gauges (mm).

3.2 Identifying tide–surge interaction

To determine whether tide–surge interaction is present at each of the seven tide gauge locations, we apply a modified version of the statistical method by Horsburgh and Wilson (2007) to the processed tide gauge records (Sect. 3.1). This method compares the timing of extreme residuals to the nearest HW.

To identify extreme residuals, the 99th percentile and above of R_gauge are identified. There are many clusters of R_gauge above this threshold, where a cluster refers to a collection of R_gauge consecutively measured, i.e. each measured 1 h from the last. The largest R_gauge in each cluster is identified and sorted. Starting from the largest, the sorted R_gauge are added to the set of extreme R_gauge unless the R_gauge was measured within 1 week (168 h) of another extreme R_gauge in the set. Some other thresholds used in published studies range from 12 to 60 h (Arns et al., 2020; Feng et al., 2019; Rasmussen et al., 2018; Buchanan et al., 2016; Horsburgh and Wilson, 2007). We choose a threshold of 1 week to reduce the odds of double counting long-lasting surges, as such events can last between 1 and 5 d in this region (Meteorological Service Singapore, 2024; Cannaby et al., 2016; Bricheno et al., 2015; Tkalich et al., 2009). We also found statistically significant autocorrelation in the R_gauge time series at time lags shorter than 168 h. A declustering threshold of 168 h still retains sufficient observations for our analysis.

To compare locations with predominantly semidiurnal tides to locations with predominantly diurnal tides, we split the extreme R_gauge into two groups. We define a tidal cycle as the duration from one local minimum in T(t) to the observation immediately preceding the next local minimum, as illustrated in Fig. 2b. One group of extreme R_gauge were observed during tidal cycles of 21 h or shorter, representing extreme residuals that occurred during semidiurnal cycles. The other group contains extreme R_gauge observed during tidal cycles of at least 22 h, representing extreme residuals that had occurred during diurnal tides. This separation between 21-hours-or-less and 22-hours-or-more was chosen based on the characteristics of the duration of tidal cycles at the seven locations, where the distribution of the duration of tidal cycles was found to be bimodal at each location, and the 22 h mark tends to distinguish the two modes (Fig. S4).

Table 3 Summary of the tidal characteristics at the study locations including tidal range, four tidal constituents (K₁, O₁, M₂, S₂), and the tidal form factor (F). Units, where applicable, are in metres (m). The diurnal tidal range is the difference between mean higher high water and mean lower low water and is also referred to as the great diurnal tidal range or great diurnal range (NOAA, 2000). Maximum tidal range is the greatest difference between higher high and lower low water within a single day and is much larger than the diurnal tidal range.

Figure 2 An example of how the no-TSI distribution is obtained, using three tidal cycles at Tanjong Pagar tide gauge station between 1 and 2 January 1989 (GMT). The number of hourly measurements made at h hours from HW are counted and denoted as f(h), the no-TSI distribution before normalisation and scaling by the number of events (Sect. 3.3). In this example, if an extreme residual can occur at any hour with equal probability, it is 3 times more likely to happen at HW than at 7 h before HW. This is observed in (b) where $-\mathrm{5}\le h\le \mathrm{4}$ occurs thrice, $h=-\mathrm{6}$ and h=5 occur twice, and $h=-\mathrm{7}$ and h=6 occur once, resulting from the three irregular tidal cycles.

The HW in the same tidal cycle as each extreme R_gauge is identified, and the timing difference between the extreme R_gauge and the respective HW is quantified at an hourly resolution. Across the set of extreme R_gauge, the frequency of extreme R_gauge found a certain number of hours relative to HW, h, is counted. This frequency is plotted as a frequency distribution in the form of a histogram (Fig. 3). The magnitude of each extreme R_gauge is represented by colour in the frequency distribution. Boxplots below the histograms show the median and its 95 % confidence interval, the interquartile range (IQR), the range that extends up to 1.5× IQR from the limits of the IQR, and the outliers (Fig. 3).

We use the frequency distribution to test the null hypothesis that assumes that there is no tide–surge interaction. To do so, we test whether the frequency distribution is drawn from a null distribution representing a scenario where extreme events are equally likely to occur at any stage of a tidal cycle.

Figure 3 Frequency distribution – the number of extreme events that have occurred at a certain number of hours relative to HW – of extreme residuals (R_gauge) during semidiurnal tidal cycles (Sect. 3.2). The plots are truncated at ±8 h from HW. The frequency distribution is compared to the no-TSI distribution, shown in grey, to determine the presence of tide–surge interaction (Sect. 3.3). Summary statistics of the frequency distribution are shown using the horizontal notched boxplot (Sect. 3.2). The whiskers of the boxplots at (b) Johor Baharu, (e) Kuantan, (f) Cendering, and (g) Geting extend beyond ±8 h from HW, and their full extents are shown in Fig. 4.

3.3 No-tide–surge interaction distribution

In existing studies, the null distribution is a uniform distribution with a value of $n/\mathrm{13}$ at $h=-\mathrm{6},-\mathrm{5},\mathrm{\dots },\mathrm{0},\mathrm{\dots },\mathrm{5},\mathrm{6}$, where n is the total number of extreme events (Horsburgh and Wilson, 2007; Haigh et al., 2010; Antony and Unnikrishnan, 2013; Feng et al., 2019; Costa et al., 2023). It assumes that tidal cycles are always 13 h long (i.e. semidiurnal) and that HW is always at the midpoint of this cycle (i.e. tides are always symmetric). Hence, if there are n extreme events, we expect $n/\mathrm{13}$ events to happen at each hourly band if extreme events are equally likely to occur at any stage of a tidal cycle. However, as tides are not always semidiurnal and symmetric, the uniform distribution is not the most-suitable null distribution to represent the null hypothesis. Instead, we propose a no-tide–surge interaction distribution or no-TSI distribution as the null distribution. The no-TSI distribution is the expected frequency distribution in the absence of tide–surge interactions. Figure 2 illustrates how this distribution is obtained. The no-TSI distribution is empirically derived from T(t), the tidal level obtained by applying UTide to the tide gauge records at each location. It is a distribution that depends on the local tidal characteristics and hence is location specific. For a given location, the null distribution is generally non-uniform because the length of the tidal cycle varies. For example, tidal cycles of 14 h or longer are relatively rare at Tanjong Pagar, so randomly selected times that occur at 7 h from the nearest tidal HW will also be relatively rare (Fig. 3a). This leads to non-uniform sampling of the number of hours from the nearest tidal HW. The no-TSI distribution corresponds to uniform sampling in time, which is non-uniform with respect to the number of hours from the nearest tidal HW. Thus, the no-TSI distribution is obtained by counting the number of tide-gauge observations found at a certain number of hours relative to HW, as illustrated in Fig. 2.

The no-TSI distribution allows us to account for the complex mixed tides at each location, which causes the tidal cycles at the coast to have a period of any duration from 8 to 26 h (Fig. S4). Let f(h) be the number of tide gauge measurements collected at h hours from HW. Assuming independence between tide and surge, the probability that an extreme event will be found at h is $p_h=f\left(h\right)/{\sum }_{h}f\left(h\right)$: the normalised frequency of f(h). Note that p_h is a probability mass function over the domain h: a function that tells us the probability of an extreme event happening at h is p_h. Letting n be the number of extreme events that occurred at this tide gauge over its length of records, the no-TSI distribution is simply n⋅p_h.

We can also find the 95 % confidence interval of the no-TSI distribution. Assuming that extreme events are mutually independent, the probability that k_h out of n extremes are found at h hours from HW follows the binomial distribution k_h∼Bin(np_h). We compute the 2.5th and 97.5th percentiles of the binomial distribution at each h to obtain the 95 % confidence interval of the no-TSI distribution.

For testing the null hypothesis, a bootstrapping method is used to calculate p values (Appendix A) with 1 000 000 bootstrap samples. To account for the family-wise error rate due to our multiple comparisons at seven locations, we apply the Bonferroni correction to our chosen significance level of 0.05 and require a p value below $\mathrm{0.05}/\mathrm{7}=\mathrm{0.007}$ and $\mathrm{0.05}/\mathrm{4}=\mathrm{0.0125}$ to reject the null hypothesis when testing for tide–surge interaction during semidiurnal and diurnal cycles, respectively. This method is chosen over the usual approach of the chi-squared test due to our relatively small sample size n. As the chi-squared test approximates the binomial distribution of k_h with a normal distribution, a common rule of thumb requires $n\cdot p_h\ge \mathrm{5}$ for a decent approximation. This criterion is not satisfied for most values of h at all seven locations.

3.4 Semi-empirical model

One effect of tide–surge interaction is tidal-phase alteration, where the surge-caused increase in water depth advances the timing of HW. Our semi-empirical model aims to investigate the first-order effects of tidal-phase alteration on the height and timing of residuals due to storm surges driven by wind and the inverse barometer effect. To do so, we first consider the regional context. Singapore and the east coast of Peninsular Malaysia are located within the Sunda Shelf (Fig. 1). The typical bathymetry of 40–80 m in this region of the Sunda Shelf terminates about 700 km away from the coast, declining steeply at the edge of the Sunda Shelf where it borders the South China Sea, which has depths of 4000 m and deeper. It is within the shallow region that tide–surge interaction is the strongest. Hence, for our semi-empirical model, we estimate the effects of tidal-phase alteration due to surges generated within the ocean bounded by the red rectangle in Fig. 1. A rectangular region was chosen for simplicity, with one edge passing through the Tanjong Pagar and Geting gauge locations and the remaining three edges encompassing as much of the shallow shelf as possible. The edges of the rectangle parallel to the Malaysian coast are roughly 833 km long and are separated by 6.5° of longitude, making the edges perpendicular to the coast roughly 759 km long. To compare the storm surge at the seven tide gauge locations of interest, we assume that the surges at these locations result from forcings over the same region. We model surges as

$$\begin{array}{}\text{(3)}& R_{\mathrm{sum}}=R_{\mathrm{baro}}+R_{\mathrm{wind}}+R_{\mathrm{phase}},\end{array}$$

where R_baro is the water-level response to the inverse barometer effect,  R_wind is the wind-driven surge, and R_phase is the contribution by tidal-phase alteration. R_sum corresponds to a model estimate of R_gauge.

To obtain R_baro, we assume hydrostatic balance and use the inverse barometer effect using a 25 h moving average of pressure:

$$\begin{array}{}\text{(4)}& R_{\mathrm{baro}}=-{\displaystyle \frac{\stackrel{\mathrm{‾}}{\langle \mathrm{\Delta }P{\rangle }_A}}{\mathit{\rho }g}},\end{array}$$

where ΔP is the deviation of sea-level pressure from the global mean, 〈ΔP〉_A represents a spatial average of ΔP over the ocean within the red rectangle in Fig. 1, and 〈ΔP〉_A represents a 25 h moving average of 〈ΔP〉_A. ρ is the density of seawater, g is the gravitational acceleration, and $\mathrm{1}/\mathit{\rho }g=\mathrm{9.9}\times {\mathrm{10}}^{-\mathrm{5}}$ m Pa⁻¹ (Gregory et al., 2019; Pugh and Woodworth, 2014c).

The unique topography of this region, with an extensive area that is shallow and relatively uniform in depth, has led us to adopt a simplified version of the sea-level gradient equation to estimate R_wind using wind velocity:

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\partial \mathit{\zeta }}{\partial x}}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{air}}{C}_{\mathrm{d}}\left|\mathit{W}\right|\mathit{W}\cdot \widehat{\mathit{x}}}{\mathit{\rho }gD}},\end{array}$$

where we have assumed that the coastal sea is shallow enough to keep only terms with water depth in the denominator but is deep enough to justify ignoring bottom stress (Pugh and Woodworth, 2014c). In Eq. (5), ζ is the sea level and x is the spatial displacement in a specified direction, making $\partial \mathit{\zeta }/\partial x$ the sea-level gradient along x. ρ_air is the density of air, C_d is the drag coefficient at the sea surface, ρ is the density of seawater, g is the gravitational acceleration, D is the undisturbed water depth, W is the wind velocity vector, and $\widehat{\mathit{x}}$ is a unit vector parallel to x. x is defined as perpendicular to the edge across Tanjong Pagar and Geting in Fig. 1. In this region, g=9.78 ms⁻². Bathymetry influences the height of surges through the term D in the denominator of Eq. (5), causing shallow regions to experience larger surges when subject to the same wind forcing.

We further assume that ρ_air, C_d, ρ, and g are spatially and temporally homogeneous over the bounded region (Gregory et al., 2019) and that C_d is independent of wind speed (Wróbel-Niedzwiecka et al., 2019). Surges are not only a product of instantaneous winds but also partially a result of winds over a past number of hours. To account for the time necessary for the winds over the bounded region to cause surges at the tide gauge locations, we use a 25 h running average of Eq. (5) to estimate the wind-driven surge. Based on those assumptions, we estimate R_wind by discretising Eq. (5) and numerically integrating over the 759 km along $\widehat{\mathit{x}}$:

$$\begin{array}{}\text{(6)}& R_{\mathrm{wind}}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{air}}{C}_{\mathrm{d}}}{\mathit{\rho }g}}\stackrel{\mathrm{‾}}{\int {\displaystyle \frac{\left|\mathit{W}\right|\mathit{W}\cdot \widehat{\mathit{x}}}{D}}\mathrm{d}x}=k\stackrel{\mathrm{‾}}{{〈{\displaystyle \frac{\left|\mathit{W}\right|\mathit{W}\cdot \widehat{\mathit{x}}}{D}}〉}_A},\end{array}$$

where W is the wind velocity with its zonal (u) and meridional (v) components obtained from the hourly 10 m winds of ERA5 (Hersbach et al., 2018), $\widehat{\mathit{x}}$ is a unit vector along x pointing towards the Malaysian east coast as shown in Fig. 1, $\langle \left|\mathit{W}\right|\mathit{W}\cdot \widehat{\mathit{x}}/D{\rangle }_A$ represents a spatial average of $\left|\mathit{W}\right|\mathit{W}\cdot \widehat{\mathit{x}}/D$ over the ocean within the red rectangle in Fig. 1, $\stackrel{\mathrm{‾}}{\langle \left|\mathit{W}\right|\mathit{W}\cdot \widehat{\mathit{x}}/D{\rangle }_A}$ represents a 25 h moving average of $\langle \left|\mathit{W}\right|\mathit{W}\cdot \widehat{\mathit{x}}/D{\rangle }_A$, $k={\mathit{\rho }}_{\mathrm{air}}{C}_{\mathrm{d}}{L}_{\mathrm{wind}}/\mathit{\rho }g$, and L_wind=759 km is the distance over which Eq. (5) was integrated. While most of the constants in k are known, C_d is often used as a final tuning parameter in non-tidal barotropic models (Zweers et al., 2012; Kurniawan et al., 2015). Using Eq. (3) and Eq. (6) and considering that |R_phase|≪|R_wind|, we fit $\stackrel{\mathrm{‾}}{\langle \left|\mathit{W}\right|\mathit{W}\cdot \widehat{\mathit{x}}/D{\rangle }_A}$ to R_gauge−R_baro. We use a simple linear regression to obtain the constant k and our estimate for R_wind at each location.

We obtain our model storm surge as

$$\begin{array}{}\text{(7)}& R_{\mathrm{surge}}=R_{\mathrm{wind}}+R_{\mathrm{baro}}=k\stackrel{\mathrm{‾}}{〈{{\displaystyle \frac{\left|\mathit{W}\right|\mathit{W}\cdot \widehat{\mathit{x}}}{D}}}_A〉}-{\displaystyle \frac{〈\stackrel{\mathrm{‾}}{\mathrm{\Delta }{P}_A}〉}{\mathit{\rho }g}}.\end{array}$$

To estimate R_phase, we consider the influence of depth on the propagation speed of tidal waves. With the speed of tidal waves given by $c=\sqrt{gD}$ (Pugh and Woodworth, 2014c) and treating R_surge as a perturbation of the undisturbed water depth D, the tide advancement time caused by a change in D due to R_surge is

$$\begin{array}{}\text{(8)}& \mathrm{\Delta }t={\displaystyle \frac{L_{\mathrm{tide}}}{\sqrt{gD}}}-{\displaystyle \frac{L_{\mathrm{tide}}}{\sqrt{g\left(D+R_{\mathrm{baro}}+R_{\mathrm{wind}}/\mathrm{2}\right)}}},\end{array}$$

where L_tide is the distance travelled by the tidal wave along the shelf. This considers that R_baro is the water level response to ΔP averaged over the ocean within the red rectangle in Fig. 1, while R_wind is obtained by integrating the sea-level gradient over the same area. Hence, we approximate the average depth perturbation due to R_wind over this area as $R_{\mathrm{wind}}/\mathrm{2}$. We also assume that the tides travel straight towards the coast in the direction of $\widehat{\mathit{x}}$, allowing us to equate L_tide=L_wind. We then calculate the effects on residual height due to tide advancement from R_surge as $R_{\mathrm{phase}}=T\left(t+\mathrm{\Delta }t\right)-T\left(t\right)$. R_phase can be viewed as a phase shift in the tidal levels, where extreme R_phase tend to cluster on the rising or falling tides instead of during tidal high or low water (Horsburgh and Wilson, 2007).

Finally, we use the procedure described in Sect. 3.2 to obtain the timing of extreme R_surge, R_phase, and R_sum relative to HW and plot their frequency distributions. To determine their contributions to tide–surge interaction, the bootstrapping method is used to identify the presence of tide–surge interaction in R_surge, R_phase, and R_sum. The frequency distributions of R_sum are compared to the frequency distributions of R_gauge to evaluate whether R_sum can reproduce the tide–surge interaction found in the tide-gauge records.

4 Results and discussion

4.1 Observed tide–surge interaction

The frequency distributions of R_gauge at all seven locations are shown in Fig. 3, all of which deviate significantly from their respective no-TSI distribution and provide evidence of tide–surge interaction at every location. Based on how we have defined extremes, we find that extreme R_gauge range from 362 mm at Tanjong Pagar (Fig. 3a) to 1195 mm at Geting (Fig. 3g), with the smallest extremes occurring at the eastern stations and increasing as we travel north. This is expected, as the northern part of the coast is more prone to larger storm surges due to tropical depressions (Mohd Anuar et al., 2023). The extremes, especially the largest extremes at each location, are unlikely to happen close to HW at all seven locations. This means that while the R_gauge exceeding the 99th percentile can occur close to HW, the peak of each cluster of exceedances is unlikely to be found in the time window close to HW. This also means that the largest residuals are unlikely to coincide with high tides to form large storm tides. Across the locations, this time window generally begins 2–3 h before HW and ends 3–5 h after HW. Beyond this time window, the frequency of extreme R_gauge increases, giving us frequency distributions that are mostly bimodal. At the four southernmost stations, we find the primary mode of their frequency distributions before HW: at −5 h at Tanjong Pagar (Fig. 3a) and Sedili (Fig. 3c), at −6 h at Johor Baharu (Fig. 3b), and at −7 h at Tioman (Fig. 3d). Outside Southeast Asia, similar signals have been found at Port Otago in New Zealand (Costa et al., 2023); Shijiusuo, Lianyungang, and Xiamen along the coast of China (Feng et al., 2019); Hiron Point at the Bay of Bengal, India (Antony and Unnikrishnan, 2013); and Aberdeen, North Shields, Immingham, Cromer, and Sheerness on the North Sea (Horsburgh and Wilson, 2007). At Kuantan, Cendering, and Geting (Fig. 3e–g), the primary mode is found after HW, at +6, +6, and +4 h, respectively. Outside Southeast Asia, such signals have appeared less commonly but have been found at Onehunga in New Zealand (Costa et al., 2023) and at Kaohsiung and Zhapo in China (Feng et al., 2019).

Comparing our results between the seven tide gauges along Singapore and the east coast of Peninsular Malaysia, we find a spatial pattern in the tide–surge interaction. At the southernmost stations of Tanjong Pagar, Johor Baharu, and Sedili (Fig. 3a–c), the mass of the frequency distribution is heavily concentrated around their primary modes, which are found before HW. At Tioman (Fig. 3d), the primary mode is still found before HW, but the secondary mode has a distinctly heavier weight than the previous three stations. At Kuantan (Fig. 3e), tide–surge interaction has crossed over to another regime where the primary mode occurs after HW, but the secondary mode found before HW still carries comparative weight. At Cendering and Geting (Fig. 3f–g), the northernmost stations, the primary mode after HW is much heavier than the secondary mode before HW. This spatial pattern can also be seen using the boxplots, which are compiled in Fig. 4a. The tight interquartile range at Tanjong Pagar, Johor Baharu, and Sedili shows that the mass of their frequency distributions is concentrated at −6 to −4 h relative to HW. The larger interquartile range at Tioman shows that there is a more-equal mass between the two modes, with the median at −5 h and mode at −7 h revealing that the distribution is still heavier towards the negative values. The opposite is true at Kuantan, with its similar interquartile range to Tioman but with its median at −0.5 and mode at +6 instead. At Cendering and Geting, the lower quartile is closer to HW, showing that the difference in relative mass between the two modes has increased.

Figure 4 Compilation of the boxplots in Figs. 3 and 5. Results presented in Fig. 3 on the timing of extreme residuals are compiled in panel (a), and results presented in Fig. 5 on extreme values of R_sum are compiled in panel (b). The boxplots illustrate summary statistics of the distribution k⁽⁰⁾ at each location, where orange lines indicate the medians, notches indicate the 95 % confidence interval of the medians, blue circles indicate the modes, notched rectangles indicate the interquartile range (IQR), whiskers indicate a range that extends up to 1.5× IQR from the limits of the IQR, and black circles indicate outliers outside this range.

A transition in the tide–surge interaction characteristic based on a significant (95 % confidence interval) change in the medians from negative to positive values in Fig. 4a happens between Tioman and Geting, which are located about 440 km apart. This is a relatively short distance compared to studies that have seen such transition, such as 700 km between Zhapo and Xiamen in China (Feng et al., 2019) and 1650 km between Otago and Onehunga in New Zealand (Costa et al., 2023). Feng et al. (2019) suggested that this transition may be related to the ratio between the amplitude of tidal constituents M₂ and K₁, where extremes tend towards before HW at locations with larger M₂ and after HW when M₂ is small. We observe the same transition in our study region, where $M_{\mathrm{2}}/K_{\mathrm{1}}<\mathrm{1}$ at the northernmost stations of Cendering and Geting, $M_{\mathrm{2}}/K_{\mathrm{1}}=\mathrm{1}$ at Kuantan, and $M_{\mathrm{2}}/K_{\mathrm{1}}>\mathrm{1}$ at the remaining four stations to the south (Table 3).

Using the quantitative test (hypothesis testing against the no-TSI distribution), we verify that the p value of obtaining the frequency distribution from the no-TSI distribution is below the required level at all seven locations, allowing us to reject the null hypothesis that the frequency distribution was drawn from the no-TSI distribution at a significance level of 0.05 (Fig. S5). This provides further evidence supporting the presence of tide–surge interaction during the semidiurnal tides at all seven tide gauge locations studied. The presence of tide–surge interaction at these gauges was expected given the shallow regional bathymetry and the proximity of tide gauges to narrow waterways connecting the Pacific Ocean, Indian Ocean, and Java Sea. Locations like the North Sea, English Channel, and some parts along the coast of China, where tide–surge interactions have been observed, also have similar geographical properties.

During diurnal tides, we found no extreme R_gauge at Tanjong Pagar, Johor Baharu, and Sedili. Few extreme R_gauge were found during diurnal tidal cycles at Tioman and Kuantan (Fig. S6a–b), while many were found at Cendering and Geting (Fig. S6c–d). This is expected, as Cendering and Geting experience mainly diurnal tides, in contrast to the other locations. The lack of observations at the five other stations is due to the rare occurrence of diurnal tides at these locations, which is shown in Fig. S4. The extreme R_gauge during diurnal tides shares the same spatial characteristics as the extremes during semidiurnal tides, starting from 410 mm in the south at Tioman (Fig. S6d) and increasing up north to 992 mm at Geting (Fig. S6g). No extremes are found within 2 h of HW. Two observations are available at Tioman, which are too few to confidently determine the presence of tide–surge interaction even though both observations were found at least 6 h after HW. The same can be said for Kuantan, where seven observations are available and were mostly found at least 6 h after HW. Nonetheless, we calculate their p values and compare them to $\mathrm{0.05}/\mathrm{4}=\mathrm{0.0125}$. We find that their p values are insufficient to reject the null hypothesis that the observed frequency distribution was drawn from the no-TSI distribution at a significance level of 0.05 and fail to provide evidence of tide–surge interaction during diurnal tidal cycles at these two locations. The results at these two locations are more likely due to a lack of sufficient data than to an indication of the absence of any tide–surge interaction, and we expect to see tide–surge interaction based on the properties of the geography discussed earlier. At Cendering and Geting, we continue to see the pattern where extreme R_gauge are unlikely to happen close to HW. This leads to a bimodal distribution with the primary modes at both locations found after HW, like their semidiurnal counterparts. The mode is 14 h after HW at Cendering and 13 h after HW at Geting. Respective p values provide evidence of tide–surge interaction at both locations (Fig. S7), as expected at these locations.

4.2 Semi-empirical model results

We obtain R_surge by fitting $\stackrel{\mathrm{‾}}{\langle \left|\mathit{W}\right|\mathit{W}\cdot \widehat{\mathit{x}}/D{\rangle }_A}$ to R_gauge−R_baro as described in Sect. 3.4. R_surge has a correlation of 0.7–0.8 with the tide gauge residuals (Fig. S8). This corresponds to an explained variance (coefficient of determination) of 0.5–0.6. This suggests that R_surge is an adequate proxy for storm surges.

We obtain the timing of extreme R_surge as a frequency distribution using the procedure described in Sect. 3.2 and compare it to its no-TSI distribution using hypothesis testing to determine whether there is any signal of tide–surge interaction in R_surge. The reason for doing so is to show that the observed tide–surge interaction between R_gauge and tide is not caused by any correlation between wind and tide or wind and atmospheric sea-level pressure, which are generated by independent processes. The validation of this assumption would imply that the observed dependence between R_gauge and tide is not caused by possible correlation with a common third independent variable but that tide–surge interaction is indeed present. We find that the resulting frequency distributions for R_surge do not deviate significantly from their no-TSI distribution and provide no evidence of dependence between R_surge and tide (Fig. S9–12), indicating the absence of such a confounding variable.

To estimate the influence of tidal-phase modulation, we compute Δt using Eq. (8) and then calculate $R_{\mathrm{phase}}=T\left(t+\mathrm{\Delta }t\right)-T\left(t\right)$ (Sect. 3.4). We apply the procedure of Sect. 3.2 to obtain the extremes of R_phase and find a clear dependence between R_phase and tide (Fig. S13–16). During semidiurnal tidal cycles, extreme values of R_phase are mostly found 2–4 h before HW at all seven locations (Fig. S13). During diurnal cycles, extreme values are mostly found 3–5 h before HW (Fig. S15). As with R_gauge and R_wind, we found no extremes during diurnal cycles at Tanjong Pagar, Johor Baharu, and Sedili.

The prevalence of extreme R_phase within a narrow window of time relative to HW is due to R_phase being largest at one-fourth of a tidal cycle before HW, as illustrated in Fig. S17. As the natural period of a semidiurnal tidal cycle is about 12–13 h, a sinusoidal tidal waveform has the steepest gradients about 3 h from its local maxima. This results in the tidal waveform having the greatest difference from a slightly horizontally displaced copy at close to 3 h from HW. Horsburgh and Wilson (2007) describe this mechanism in detail. However, extreme values of R_phase are not found at 6 h before HW during diurnal tides. This is because a sinusoidal wave with amplitude A and frequency ω has a gradient that is proportional to the product Aω. Since semidiurnal components of tides have about twice the frequency of their diurnal counterparts, diurnal constituents need to have at least twice the amplitude of semidiurnal constituents to have the same or stronger influence on R_phase. The tidal form factor F, which computes the ratio between the amplitudes of the main diurnal and semidiurnal constituents, can provide a measure of this influence. While the tides at Cendering and Geting are mainly diurnal, the tidal form factors at both locations are lower than 2, which indicates that the diurnal constituents have amplitudes less than twice that of the semidiurnal constituents (Table 3). Therefore, the semidiurnal signal has a stronger influence on the timing of extreme R_phase.

The frequency distributions of R_sum at the seven locations during semidiurnal tidal cycles are shown in Fig. 5. The frequency distributions and p values for R_sum suggest the presence of tide–surge interaction at Tanjong Pagar, Johor Baharu, Sedili, Tioman, Kuantan, and Cendering (Figs. 5a–f and S18a–f). No significant interaction is identified at Geting (Figs. 5g and S15g). This can also be seen in Fig. 4b where the interquartile range of the frequency distributions at Geting does not deviate from zero, while the interquartile range at the other six locations does. Figure 5 reveals that our semi-empirical model predicts that the frequency distribution of R_sum at all seven locations has a single mode as opposed to the bimodal frequency distribution of the tide gauge residuals. At the six locations where R_sum provides an indication of tide–surge interaction, the modes of their frequency distributions of R_sum lie within 2–4 h before HW, following R_phase. During diurnal tidal cycles, we find evidence of tide–surge interaction only at Cendering (Figs. S19c and 20c), while no evidence of tide–surge interactions can be seen at Tioman, Kuantan, and Geting (Figs. S19a, b, d, S20a, b, d). As with their R_gauge counterparts, the hypothesis testing results at Tioman and Kuantan are likely due to a lack of sufficient observations (Fig. S19a–b). The result at Geting, which is negative for tide–surge interactions, is likely due to its relatively weaker diurnal constituents. This can be seen from the tidal form factor F in Table 3, where the form factor of 1.56 at Geting is lower than the 1.90 at Cendering, which tested positive for tide–surge interactions.

Figure 5 The frequency distribution for extreme values of R_sum and the no-TSI distribution during semidiurnal tidal cycles, truncated at ±8 h from tidal HW. Truncated horizontal notched boxplots illustrate some summary statistics of the frequency distribution, and their full extents are shown in Fig. 4.

We find that R_phase can significantly influence the distribution of the extreme values of R_sum, indicating that the process of tidal-phase alteration – where surges perturb the depth of the water body and influence the propagation speed of the tide – produces a significant and measurable tide–surge interaction at six locations. This is despite R_phase contributing to < 1 % of the variance of R_sum at all seven tide gauge locations (Fig. S21). By an alternative metric, the ratio between the standard deviation of R_phase and the standard deviation of R_surge is only 0 %–2 %. Thus, while the magnitude of R_sum is effectively fully dependent on R_surge, the timing of its largest values is dependent on R_phase, indicating the significant contribution of tidal-phase alteration to tide–surge interaction. Our findings complement those of Chen et al. (2012), who found that the influence of tide–surge interactions on large surges is negligible, although tide–surge interaction may alter the time of tidal high and low water. Accounting for tide–surge interaction would be important in applications such as forecasting, where timing is important and has been found to improve water level predictions in this region and in other parts of the world (Antony et al., 2020; Fernández-Montblanc et al., 2019; Kurniawan et al., 2015; Idier et al., 2012).

Comparing Fig. 4b to a, our semi-empirical model (Fig. 4b and 5) is unable to accurately predict all the characteristics of tide–surge interaction found in R_gauge (Figs. 3 and 4a). At Tanjong Pagar, Johor Baharu, Sedili, and Tioman, extreme residuals typically occur even earlier than what our model suggests. Tidal-phase alteration shifts the timings of extreme residuals towards times where tides are rising or falling the quickest, and no further. Our model is also unable to produce the bimodal distribution found in the tide gauge data illustrated in Fig. 3 and fails to produce statistically significant tide–surge interaction at Geting. This suggests that the mechanism of tidal-phase alteration cannot fully explain the observed tide–surge interaction and that additional explanations are required to fully account for the observations.

One possible contributor is that shallower water during low tides can result in higher surges from near-shore winds (Pugh and Woodworth, 2014c). Changes in tides due to tidal-phase alteration can affect the height of surges through Eq. (5), while surges further contribute to tidal-phase alteration through Eq. (8), causing the mutually interacting process of mutual phase alteration between tides and surges. Kurniawan et al. (2014) found that the tidal cycle does influence non-tidal residuals in this region. This could result in a further shift in the model extremes towards tidal low water, which would be closer to the observed tide–surge interaction. The impact of non-linear effects resulting from bottom friction and momentum advection also cannot be overlooked in coastal areas and are significant in this region (Kurniawan et al., 2015). We were not able to consider the effects of water flow between the Pacific and Indian oceans through the Singapore and Johor straits on the results at Tanjong Pagar and Johor Baharu. We also could not account for the tide–surge–river interactions at Sedili, the presence of Tioman island at Tioman, and how the Gulf of Thailand affects tide–surge interactions at the northern tide gauges. Our semi-empirical model can reveal the macroscopic properties of tide–surge interactions but is unable to reveal detailed spatial and temporal properties at high resolution like alternative numerical models. To reveal such details, the use of dynamical models or partial differential equation solvers would be needed. For instance, Kurniawan et al. (2015) and Tay et al. (2016) used a multi-scale depth-integrated hydrodynamic model that simulated hydrodynamics at a higher spatial and temporal resolution and found a significant contribution by tide–surge interactions, on the order of 10 cm, to non-tidal residuals. However, such numerical models are complex, difficult to tune, and computationally intensive. While it does not replace hydrodynamical models, our semi-empirical model offers theoretical insights by providing a clear and transparent explanation of tidal-phase alteration and its impact at a low computational cost.

The semi-empirical model is less applicable elsewhere because it has been designed for the study region. However, the semi-empirical model could be modified for other coastlines situated along a shallow shelf. For instance, the physical characteristics of our study region – a semi-enclosed shallow region leading out to steep depths – can be found on a smaller spatial scale at the Yellow Sea and the North Sea, and our semi-empirical model may find similar success there.

Testing our no-TSI distribution against the uniform distribution, we found false-positive signals for tide–surge interactions in R_surge at Tanjong Pagar and Johor Baharu and in R_sum at Geting when using the uniform distribution. Using the no-TSI distribution at locations that have similar frequency distributions such as Green Island in New Zealand (Costa et al., 2023), Newlyn in the UK (Haigh et al., 2010), and several locations in China (Feng et al., 2019) may lead to a different conclusion regarding the existence of tide–surge interactions.

The existence of strong tide–surge interactions in this region discourages the use of the revised joint probability method for estimation of sea-level extremes, which works well in locations with short records but tends to overestimate return levels at locations with strong tide–surge interactions (Arns et al., 2020; Olbert et al., 2013; Haigh et al., 2010; Tawn, 1992; Tawn and Vassie, 1989). The use of skew surge for the estimation of extremes may be more suitable than surge or residuals for this region (Williams et al., 2016).

5 Conclusions

We have introduced the no-TSI distribution to be used in determining the presence of tide–surge interaction. The no-TSI distribution can account for irregular tidal cycles that can lead to non-uniform sampling with respect to the number of hours from HW. Hence, the no-TSI distribution is generally not a uniform distribution. When determining the presence of tide–surge interaction, the observed frequency distribution should be compared to the no-TSI distribution instead of to a uniform distribution.

Analysis of tide gauge records using the no-TSI distribution provides evidence of tide–surge interaction at all seven tide-gauge locations along Singapore and the east coast of Peninsular Malaysia. The observed interactions have a smooth spatial dependence along the coastline. During semidiurnal tidal cycles at the southernmost location of Tanjong Pagar, extreme residuals are mostly found around 5 h before tidal HW, with a much smaller number of extremes occurring after HW. Moving northwards, similar patterns are found at Johor Baharu and Sedili, until we reach Tioman where extremes are mostly found before HW, but many extremes can also be found after HW. Northwards from Kuantan, most extremes are found after HW. The strong tide–surge interaction in this region suggests that return levels should be estimated using skew surge, not hourly residuals.

To investigate the contribution of tidal-phase alteration, we proposed a semi-empirical model. We used 10 m winds from ERA5 to estimate the effects of tidal-phase alteration. At the six southern stations, Tanjong Pagar, Johor Baharu, Sedili, Tioman, Kuantan, and Cendering, we found that R_phase – the residual component caused by the advancement of tidal HW due to surges – can significantly alter the timing of extremes despite being responsible for less than 1 % of the variance of the total residual. This demonstrates the effects of tidal-phase alteration on the timing of extreme residuals.

Our model has explored one of the underlying mechanisms behind tide–surge interaction but is not designed to forecast water level or extreme events. A forecast model would require much more accurate modelling of storm surge. The inclusion of other underlying mechanisms of tide–surge interaction, such as the effect of tidal level on surge generation, would also be beneficial. Knowledge of the interplay between tide–surge interaction and extreme sea levels can aid in the design of effective strategies for coastal planning, risk assessment, and mitigation measures and will benefit from more comprehensive analyses. When studying extreme sea level in Southeast Asia, the relatively short length of available tide gauge records poses a challenge, providing a focus for further research.

Appendix A: Calculating the p value with bootstrapping

Here, our objective is to determine whether the frequency distribution obtained in Sect. 3.2 is drawn from the no-TSI distribution. To do so, we calculate the p value: the probability of obtaining distributions from the no-TSI distribution that are at least as extreme as the frequency distribution. We can estimate this probability by considering the following: if we draw a random sample from the no-TSI distribution, the probability of drawing an outcome that is equally as probable or less probable than the frequency distribution is, by the definition of p values, equal to p. Hence, we can obtain a good estimate of the p value by drawing many samples from the no-TSI distribution and calculating the ratio of samples that are equally as probable or less probable than the frequency distribution to be drawn from the no-TSI distribution.

To obtain one bootstrap sample following the notation in Sect. 3.3, we draw one sample of size n from the normalised no-TSI distribution p_h. To proceed, we need to calculate the probability of obtaining this specific outcome. Labelling this outcome as k, k follows a multinomial distribution with probability given by $p_{\left\{k\right\}}=\frac{n\mathrm{!}}{{\prod }_h{k}_h\mathrm{!}}{\prod }_h{p}_h^{k_h}$ where k_h is the number of times (out of n times) each possible h is drawn. We can also compute $p_{\left\{k^{\left(\mathbf{0}\right)}\right\}}$, the probability of obtaining the frequency distribution from the no-TSI distribution, in the same way using the multinomial distribution. After obtaining 1 000 000 bootstrap samples and calculating their probabilities, we can achieve a good estimate of the p value as

$$\begin{array}{}\text{(A1)}& p={\displaystyle \frac{\mathit{\#}\text{of samples where}{p}_{\left\{k\right\}}\le p_{\left\{k^{\left(\mathrm{0}\right)}\right\}}}{\mathrm{1}\mathrm{000}\mathrm{000}}}.\end{array}$$