Tide-Surge Interaction in the Pearl River Estuary: A Case Study of Typhoon Hato

In this study, the characteristics and mechanisms of tide-surge interaction in the Pearl River Estuary (PRE) during Typhoon Hato in August 2017 are studied in detail using a 3D nearshore hydrodynamic model. The wind field of Typhoon Hato is firstly reconstructed by merging the Holland parametric tropical cyclone model results with the CFSR reanalysis data, which enables the model to reproduce the pure astronomical tides and storm tides well; in particular, the distinctive oscillation pattern in the measured water levels due to the passage of the typhoon has been captured. Three different types of model runs are conducted in order to separate the water level variations due to the astronomical tide, storm surge, and tide-surge interactions in the Pearl River Estuary. The results show the strong tidal modulation of the surge level, as well as alteration of the phase of surge, which also changes the peak storm tidal level, in addition to the tidal modulation effects. In order to numerically assess the contributions of three non-linear processes in the tide-surge interaction and quantify their relative significance, the widely used “subtraction” approach and a new “addition” approach are tested in this study. The widely used “subtraction” approach is found to be unsuitable for the assessment due to the “rebalance” effect, and thus the new “addition” approach is proposed along with a new indicator to represent the tide-surge interaction, from which more reasonable results are obtained. Detailed analysis using the “addition” approach indicates that the quadratic bottom friction, shallow water effect, and nonlinear advective effect play the first, second, and third most important roles in the tidal-surge interaction in the estuary, respectively.


INTRODUCTION
Storm surges are abnormal variations of sea level driven by atmospheric forcing associated with extra-tropical storms or tropical cyclones (also known as hurricanes and typhoons). When combined with the astronomical tide, storm surges often result in extreme water levels and can bring devastating damage to coastal areas, especially for those low-lying areas bordered by extensive continental shelves and exposed to the regular passage of typhoons and storms (Bertin et al., 2012;Zhang et al., 2017). To be able to predict the peak water levels, operational systems, and research studies often superpose an atmospheric-only forced storm surge onto the astronomical tide without considering the effect of tide-surge interaction (Peng et al., 2004;Bobanović et al., 2006;Graber et al., 2006). However, tide-surge interaction has long been recognized as one of the most important contributors in the storm surges and peak water levels in coastal regions (Proudman, 1955(Proudman, , 1957Doodson, 1956;Bernier and Thompson, 2007;Zhang et al., 2010). Compared with observations, errors in a simple linear superposition of the astronomical tide with a separately computed surge are found to be up to 1-2 m (Rego and Li, 2010). Therefore, quantitative insights into tide-surge interaction are very important for the prediction of storm tide level and flood risk assessment.
It has long been recognized that tide-surge interaction is a non-linear phenomenon. Previous literature broadly focused on different aspects of the interaction, e.g., the tide-induced modulation of the phase of surge and, consequently, the variations in sea level, and the different contributions from various physical processes to the surge level. Proudman (1955) was among the first few studies to develop solutions for the propagation of an externally forced tide and surge into an estuary of uniform section. It showed that due to tide-surge interaction, the peak storm surge height occurring near to high tide was less than that which occurred near to low tide for a progressive wave. Rossiter (1961) suggested that a key mechanism of tidesurge interaction was mutual phase alteration and showed how a negative surge would retard tidal propagation, whereas a positive surge would advance the high water. Horsburgh and Wilson (2007) showed that surge generation was reduced during high water and that the surge peak was less likely to occur during high water for a large amplitude tide. Rego and Li (2010) studied the effects of tide and shelf geometry under Hurricane Rita. The results indicated that for landfall at midebb or midflood, the storm tide level was less affected, but for landfall at high tide or low tide, the peak storm tide was either reduced or increased compared to a linear superposition.
It is also widely accepted that tide-surge interaction comprises three nonlinear physical processes: (a) the non-linear horizontal and vertical advection in the momentum equations, (b) the non-linear bottom friction effect associated with the quadratic parameterization, and (c) the shallow water effect arising from the non-linear terms related to the total water depth in both the continuity and momentum equations (Tang et al., 1996;Bernier and Thompson, 2007;Rego and Li, 2010;Zhang et al., 2010Zhang et al., , 2017. However, it is difficult to separate these processes and quantify their contributions to the interaction by using observation data. Therefore, numerical models have been extensively used to examine the mechanisms of tidesurge interaction. Wolf (1978) showed that the tide-surge interaction was dominated by quadratic friction, followed by the shallow water effect and advection process. Subsequently, Wolf (1981) further demonstrated that the shallow water effect became important for small tidal range and depths of <10 m. Using a two-dimensional numerical model of the shallow-water equations, Tang et al. (1996) demonstrated that, with the tides included in the storm surge model, the sea level elevation on the North Queensland coast was generally lower than that obtained by simply adding the astronomical tides to the surge due to the quadratic bottom friction law. Rego and Li (2010) suggested that the non-linear advection dominated in a realistic simulation, while the quadratic friction was the largest in an idealized simulation. Zhang et al. (2010) studied the tide-surge interaction in the Taiwan Strait and indicated that the non-linear bottom friction was a major factor in predicting the elevation, while the non-linear advective terms and the shallow water effect contributed little.
To quantify the contributions from each of the above three processes to the tide-surge interaction, a "subtraction" approach has been widely adopted in previous studies (Tang et al., 1996;Bernier and Thompson, 2007;Zhang et al., 2010Zhang et al., , 2017. Based on a standard model that includes all three processes, this approach assesses the changes to the interaction intensity by using a reduced model in which the non-linear terms associated with one of the three physical processes are linearized or eliminated in turn from the standard model. To facilitate the quantification, various indicators have been used to represent the intensity of tide-surge interaction in different studies, e.g., the maximum positive, minimum negative, or root-mean-square of the tidesurge interaction-induced residual elevation. However, such a method is found to be defective due to the so called "rebalance" effect (Zhang et al., 2017), which means that the "subtraction" approach cannot clearly separate the contributions of these three processes and quantify their relative significance to the interaction. A new approach is therefore needed to properly reveal the individual contributions to the tide-surge interaction without interference from other processes. This is fulfilled by adopting a new "addition" approach in the present study by quantifying the interaction intensity obtained from a reduced model in which only one non-linear process is included and comparing this intensity with that obtained from the standard model (see more details in section 5). Furthermore, a new indicator of the interaction intensity is also proposed in this study, which is thought to be more appropriate for quantifying the relative importance of different physical processes in studying the mechanism of tide-surge interaction.
The Pearl River Estuary (PRE), connecting with the Pearl River at its northern end, is the largest estuary in the Pearl River Delta (PRD). Its shape is that of an inverted funnel, with a narrow neck in the north and a wide mouth opening to the South China Sea. The topography of the PRE is constituted of deep channels, shallow shoals, and tidal flats, which makes the PRE extremely vulnerable to storm surges resulting from typhoons or strong tropical cyclones. According to data from the annual tropical cyclone publication of Hong Kong Observatory (HKO, 2017), fourteen typhoons inducing high storm surges over 1 m were recorded in Hong Kong (located in the south of PRE) from 1999 to 2018, two of which caused storm surge elevations of over 2 m. One of these two events, Typhoon Hato in August 2017, generated a pronounced storm surge along the coast of the PRE. The maximum storm surge reached 1.62 m at A-Ma station in Macau, a record high in Macao since records began in 1925 , and reached 2.79 and 2.42 m at Zhuhai and Tsim Bei Tsui in Hong Kong, respectively. Observations of the water levels during the passage of Hato provided a unique dataset to assess tide-surge interactions and the relative contributions from the three different processes.
The main objectives of this study are therefore to apply a three-dimensional hydrodynamic model to identify the characteristics of tide-surge interaction in the PRE during Typhoon Hato and to quantify the relative importance of the three nonlinear effects on the tide-surge interaction. In section 2, the numerical model and model configurations used in this study are briefly described. The reconstructed wind field, modelsimulated astronomical tides, and total water levels are evaluated and validated in detail by comparing with observations in section 3. The characteristics of tide-surge interaction associated with Typhoon Hato and its impact in the PRE are studied in section 4. In section 5, the relative importance of the three nonlinear effects on the tide-surge interaction is quantified by using the newly proposed "addition" approach. Finally, the results are summarized and conclusions drawn in section 6.

The Numerical Model
In this study, a prognostic, three-dimensional coastal-ocean model developed for hydrodynamic-wave coupling (Zheng P. et al., 2017) has been applied to study the tide-surge interaction in the PRE. The model is based on the Finite-Volume Community Ocean Model (FVCOM, by Chen et al., 2003). It uses nonoverlapped triangular grids in the horizontal plane (x and y) to resolve the complex shoreline and geometry and the generalized terrain-following Sigma coordinate (s) in the vertical direction to accommodate the irregular bathymetry. The modesplit approach is used for the solution of the circulation model, in which currents are divided into external and internal modes and computed using an external and internal time step, respectively (Chen et al., 2003). After the Boussinesq and hydrostatic approximations, the 3D momentum and continuity equations used in FVCOM are as follows: where u, v, and ω are the velocity components in the x, y, and s directions, respectively, the vertical s coordinate ranges from s = −1 at the bottom to s = 0 at the free surface, D = ζ + h is the total water depth, where ζ is the surface elevation and h is the resting water depth, ζ a is the sea level displacement induced by the "inverse barometer effect, " g is the gravitational acceleration, f is the Coriolis parameter, ρ 0 and ρ are the reference water density and water density, respectively, K m and ν are the vertical eddy and molecular viscosity coefficients, respectively, and (F x , F y ) represent the horizontal momentum mixing terms in the x and y directions, respectively. In the above momentum equations (i.e., Equations 1, 2), the second, third, and fourth terms on the left-hand side are the advection terms (ADV), while the second term on the right-hand side represents the baroclinic pressure gradient force (which is neglected in this study). The surface and bottom boundary conditions for u, v, ω are given as follows: in which τ sx , τ sy and τ bx , τ by are the x and y components of surface wind and bottom stresses, respectively. The quadratic law is applied in the parameterization of both the surface wind and bottom stresses, as follows: where ρ a is the air density, C ds and C db are the surface wind stress and bottom drag coefficients, respectively, and (U w , V w ) are the wind speed components at a height of 10 m above the sea surface in the x and y directions, respectively. In FVCOM, the surface drag coefficient C ds is determined with a bulk formula as follows (Large and Pond, 1981): w is the magnitude of wind velocity; the bottom drag coefficient C db is determined by matching a logarithmic bottom layer to the model at a height of z r above the bottom, i.e., Frontiers in Marine Science | www.frontiersin.org where κ = 0.4 is the von Karman constant, z 0 is the bottom roughness parameter, and z r is a reference height above the bed, normally equivalent to half of the height of the first grid cell above the bed (e.g., z r = D/[2(N − 1)]; N is the number of vertical sigma layers). It is noted that the C db calculated as above is dependent on the total water depth, which should also include a non-linear shallow water effect. This effect is eliminated by applying a constant C db of 0.0025 in order to cleanly separate the contributions of non-linear bottom friction and the shallow water effect and also due to its negligible role in affecting the tide-surge interactions (Zhang et al., 2010).

Model Configuration in the PRE
The model domain covers the whole Pearl River Delta together with part of the South China Sea shelf. The open boundary (OB) is parallel to the coast and is placed far enough away to eliminate any boundary effects on the simulation inside the PRE (Figure 1). The resolution of the horizontal grid is ∼ 50 − 200 m within the Pearl River network and ∼ 300 − 500 m inside the PRE and decreases from the coastline (∼ 500 − 1, 000m) toward offshore. The maximum grid size at the OB is approximately 15 km. The resultant horizontal mesh contains a total of 97,602 elements and 56,993 nodes ( Figure 1B). In the vertical direction, 25 sigma layers are used, with a uniform layer thickness of about 0.2 m inside the majority of the PRE.
The model was driven by tidal forcing from the OB and atmospheric forcing (i.e., wind stress and sea level pressure, SLP) at the sea surface. Eight tidal constituents (i.e., M 2 , S 2 , N 2 , K 2 , K 1 , O 1 , P 1 , Q 1 ) from the TPXO database (Egbert and Erofeeva, 2002) were used to generate tidal water level time series to drive the model at the open boundary. The atmospheric forcing consisted of hourly 10-m wind speed and SLP with a horizontal resolution of 0.2 • (latitude/longitude), which were obtained from the National Centers for Environmental Prediction (NCEP) Climate Forecast System Reanalysis (CFSR) dataset. In order to better describe the typhoon-associated wind field and SLP, blended atmospheric forcing was used in this study by inserting an idealized wind field and SLP for a tropical cyclone, which was calculated by the Holland parametric tropical cyclone model (Holland, 1980), into the original large-scale CFSR atmospheric data (see details in section 7). In addition, high temporal resolution (hourly) observed river discharge from three upstream hydrological stations (i.e., Gaoyao, Shijiao, and Boluo) were used to represent the freshwater inputs from the West River, the North River, and the East River, respectively.
Three sets of numerical experiments were conducted to assess the model performance and to analyze the mechanism of tidesurge interaction: (a) Full run (Run-Full): The model was driven by both the tidal forcing at the OB and the blended atmospheric forcing. The resultant water level from this model run is the storm tide (ζ ST ). (b) Storm-only run (Run-SO): Only the blended atmospheric forcing was used to drive the model, while the tidal forcing was turned off. The resultant water level from this model run is called pure storm surge (ζ SO ). (c) Tide-only run (Run-TO): Only the tidal forcing was included. The resultant water level is the pure astronomical tide level (ζ TO ).
All of the above experiments started from August 1, 2017, and spun up from rest (i.e., zero velocity and undisturbed water level) for the first 4 days; the simulations were then conducted continuously through the whole of August 2017. The split-mode time-stepping method is used in this model, with a 6-s internal time step and a 1-s external time step.

Wind Speed Evaluation
As shown in Figure 1, Hato formed as a tropical depression over the sea northeast of Luzon Island on August 19, 2017, and intensified to a tropical storm over the same waters on August 20. It moved westwards across the Luzon Strait, and intensified to typhoon strength over the northeastern part of the South China Sea on August 22. After that, Hato moved westnorthwest toward the coast of China, where it intensified further and became a super typhoon in the early morning of August 23 over the sea south of Hong Kong, reaching its peak intensity with an estimated sustained wind speed of 185 km/h near its center. After making landfall at Zhuhai with severe typhoon intensity, Hato gradually degenerated into a low-pressure system on August 24. Based on the above information, a reconstructed blended wind field for Typhoon Hato was created by using the Holland parametric model (see details in the Appendix).
Compared with the original CFSR wind, the blended wind field shows a much larger wind speed near the typhoon center and a more asymmetric vortex structure, which has larger wind speeds on the right-hand side of the typhoon track due to the typhoon translation motion (Figure 2). Especially at 03:00 GMT on 23rd August, when Hato intensified to a super typhoon, the blended data (Figure 2h) clearly reproduced the much stronger typhoon intensity; by contrast, no obvious vortex structure of the typhoon was found in the original CFSR data (Figure 2d). Moreover, the locations of the typhoon center in the blended data are consistent with, while those in the CFSR data deviate more or less from (e.g., Figures 2a,b), those taken from the best track data. In order to qualitatively evaluate both the CFSR and the blended winds, the observed wind speeds from six representative wind gauge stations are used in this study, including stations 59,682 and G3599, which are located near Hato's track center, and G3598, which is relatively far away, and another two locations (i.e., G1217 and 59,479), which are near the entrance of the PRE but also not far from the tropical cyclone track, and an extra one (i.e., G1211), which is located inside the PRE.
The observed and reconstructed wind speeds at the above six stations are compared in Figure 3, in which a common feature of two distinct peaks is observed in the last 10 days of August. The first peak on August 23 results from Typhoon Hato, while the other one is due to another typhoon, Pakhar. In this research, only Hato is analyzed in detail, and thus the blended wind field is only created during its passage (i.e., August 21-24), while for the rest of the time, the blended wind field is identical to the CFSR dataset. When comparing with the observations, it can be seen that the wind speeds based on CFSR are very close to the measurements when Hato's effects are minimal, e.g., between August 15 and August 21, when the typhoon is absent at all stations, and throughout the whole period at G3598, which is far away from the typhoon center. However, the CFSR data severely underestimates the wind speed during the passage of both Typhoon Hato and Typhoon Pakhar. In contrast, the blended approach reproduces both Typhoon Hato's peak wind magnitude and timing well on the whole, although some discrepancies are still observed (e.g., G1211) due to the fact that the parametric tropical cyclone model does not account for the structural changes and wind reduction caused by local land topography. These comparative results suggest that a blended approach is able to achieve a reasonably good estimation of the peak wind stresses under a typhoon condition, while the CFSR data can only be reasonably used with minimal typhoon impacts.

Water Level Validation
The root-mean-square error (RMSE), correlation coefficient (R), and model skill (Skill) metrics were used to validate the computed water level. The RMSE indicates the average deviation of the model results from the observations and is defined as where M n and C n are the measurements and model computed results, respectively, at N discrete points. CCF and Skill evaluate the coherence between the model results and observations; a CCF or Skill value of 1 indicates a perfect agreement between the model results and measurements, whereas a value of 0 indicates complete disagreement. The CCF is given by where σ C and σ M are the standard deviations of the model results and measurements, respectively; the overbar represents the mean value. Following Willmott (1981), the Skill formulation is: The computed astronomical tide was first evaluated at four hydrological stations, Guanchong, Xipaotai, Huangpu, and Nansha, over August 5-31, 2017 (Figure 4). As shown in Table 1, the model predictions follow the reconstructed astronomical tides very well: the RMSE values at all four stations are <0.17 m, and the correlation coefficient (CCF) and model skill (Skill) FIGURE 4 | Comparisons of model-predicted (lines) with the reconstructed (circles) astronomical tides over August 5-31, 2017, at the stations of Guanchong, Xipaotai, Huangpu, and Nansha. The reconstructed astronomical tides are calculated from the tidal constituents that were obtained from long-term harmonic analysis of the observed total water levels.
are generally above 0.96. The model-predicted storm tides at the above four stations were further compared with the total water level observations, as shown in Figure 5. At all of the above four stations, the observed storm tide reaches its maximum (above 2 m) on the morning of August 23, shortly after Typhoon Hato makes landfall. Among these four stations, the recorded storm tide shows a pattern of reaching a single maximum at Huangpu and Nansha, with peak water levels of 2.92 and 3.3 m, respectively. At the other two stations (i.e., Guanchong and Xipaotai), it is interesting to observe that the recorded water level shows a double-peak pattern of "abrupt decline and then rapid rise" in a short time period just before reaching the maximum value on the August 23. This is closely related to the positions of these two stations relative to the typhoon center, which determines the local wind direction, and their relative relationship with the local geometry of the coastline. When Hato moves close but has not made landfall, these two stations are located in the right front quadrant of the typhoon, with offshore winds prevailing locally; a negative storm surge is thus produced, making the local water level drop considerably. After Hato makes a landing, the local wind direction becomes onshore in a short time, with the above two stations lying at the right rear of the typhoon center. The local water level thus increases, with a positive storm surge produced. It is the strong local wind that  leads to the significant intensity of the local drop and rise of water level, whereas it is the fast translation speed of Hato that results in the sharp change in water levels from a local minimum to the maximum value. When Hato is far away (i.e., before and after the August 23) from the local stations, the model predicted storm tide from the CFSR wind field agrees well with the observations. However, the CFSR model results severely under-estimate the maximum water levels (e.g., at Nansha station) when Hato moves close; in the meantime, it totally misses the "double-peak" pattern of water level observed at Guanchong and Xipaotai. By contrast, the model-calculated water levels from the blended data agree well with the observations during the whole passage of Typhoon Hato, with both the storm tide maxima and the above "double-peak" pattern of water level well reproduced. The model discrepancies at the time when peak storm tides occur are reduced from 1.37, 1.32, 0.46, and 1.06 m when the original CFSR data is used to 0.42, 0.08, 0.18, and 0.20 m when the blended data is used at the stations of Guanchong, Xipaotai, Huangpu, and Nansha, respectively. Therefore, significant improvements in the modelpredicted water levels were obtained in this study by using the blended data shown in the Appendix. Table 1 also shows that the CCF (Skill) at all four stations is above 0.96 (0.94), indicating an overall good agreement of the model-predicted storm tide with the observations over the whole simulation period. However, when we zoom in on the validation period for August 21-24, the CCF (Skill) reduces slightly, while the RMSE increases by more than 9 cm at all four stations. This is largely due to some physical processes being missing in the present model simulations, e.g., wave-induced setup and non-hydrostatic pressure gradients (Zhang et al., 2017). Figure 6 shows the time series of the model-predicted storm tide levels (ζ ST ), astronomical tide levels (ζ TO ), and pure surge elevations (ζ SO ) at the above four tide gauges; they are water level results from the standard experiments Run-Full, Run-TO, and Run-SO, respectively. In addition, two residual water elevations, ζ PS and ζ TSI , are also included in Figure 6. The residual water elevation ζ PS is calculated by subtracting ζ TO from ζ ST and is known as the practical storm surge, as defined in most operational storm-surge monitoring systems (Horsburgh and Wilson, 2007;Zhang et al., 2010); while ζ TSI = ζ ST − ζ TO − ζ SO is the residual elevation due to the tide-surge interaction. The model results show that the magnitudes of ζ PS near the landfall of Hato are 2-3 m at the four tide gauges and are much larger than their neighboring astronomical tidal high levels. These high water elevations overtopped the coastal sea walls, bringing a large amount of flooding to the coastal areas of the PRE .

TIDE-SURGE INTERACTION AND ITS IMPACT
Without tide-surge interaction, the practical storm surge ζ PS will be equal to the pure storm surge ζ SO . However, this is generally not true, as shown in Figure 6: the ζ PS is not equal to ζ SO during most of the time period at all four stations. The comparison of ζ PS and ζ SO shows a general feature, in that the magnitudes of ζ PS are greater near low tide but smaller near high tide than ζ SO , especially in the first tidal cycle on August 23, when the storm surge maxima occur. Similar results have also been reported in many previous studies (e.g., Horsburgh and Wilson, 2007;Rego and Li, 2010;Zhang et al., 2010Zhang et al., , 2017, reflecting the effects of tidal modulation on surge generation, which can be explained by an idealized expression for the equilibrium between sea surface slope and a constant wind stress term (Pugh, 1996) as follows: Although such an equilibrium is rarely established in the real world because the wind field changes frequently, Equation (13) illustrates the fundamental point that the wind stress is more effective in producing surges in the shallower waters, e.g., during tidal low waters. In addition to the change of magnitude, the phase of the surge can also be altered by the tide-surge interaction (tidal modulation). Previous studies (e.g., Horsburgh and Wilson, 2007;Wolf, 2009;Rego and Li, 2010), have pointed out that a reduced water depth will result in reduced phase speed both directly and indirectly due to the effects of bottom friction, as it is inversely proportional to the water depth, whereas the enhanced water depth will increase the phase speed. Consistent with the above physics, the peaks of the predicted ζ PS shown in Figure 6 arrive a bit earlier than those of ζ SO . The impact of tidal modulation (tide-surge interaction) on the storm surge and total water levels in the whole PRE can be examined in detail by reference to Figure 7, in which the distribution of the differences between the maxima of ζ PS and ζ SO (i.e., ζ max PS − ζ max SO ; Figure 7A) and the differences between the maxima of ζ ST and ζ SO + ζ TO (i.e., ζ max Figure 7B) are presented. In these figures, two notable features can be observed: firstly, the spatial distributions of both differences defined above show considerable variations in the PRE, indicating that the effect of tide-surge interaction is highly localized and spatially varying; secondly, both of the tidally modulated peak water elevations (i.e., ζ PS , ζ ST ) have higher magnitudes near the east coast but smaller magnitudes close to the west coast of the PRE than those predicted without the effects of tide (i.e., ζ SO , [ζ SO + ζ TO ]), which confirms previous studies (e.g., Brown et al., 2010;Quinn et al., 2012) in showing that tidesurge interaction can either enhance or reduce the peak surge elevations. More detailed examinations of the magnitudes show that the peak water elevations at Shenzhen Bay are significantly raised by 0.1-0.5 m due to tide-surge interaction, whereas, in the coastal area of Zhuhai and Macau, the peak water elevations are reduced by 0.2-0.4 m. From a surge-protection point of view, the increase in the water level shown in Figures 7A,B is of more practical significance, as an underestimation of the peak water elevations, e.g., near the east PRE coast in this study when the effect of tide-surge interaction is not taken into account, could lead to huge economic loss and high fatalities.
The differences in the maxima of the practical storm surge ζ PS and the pure storm surge ζ SO (i.e., ζ max PS − ζ max SO ) in Figure 7A represent the tide-surge interaction-induced changes in the magnitude of the storm surge. By contrast, the differences between the maximum elevations of ζ ST and ζ SO + ζ TO (i.e., Figure 7B include the effects from the tide-surge interaction on both the magnitudes and phases of the storm surge. The fact that the tide-surge interaction not only influences the surge level but also the peak timing of the storm surge is clearly reflected in the contrast between Figures 7A,B, which is also detailed in Figure 7C. A close examination of Figure 7C suggests that the phase alteration mainly increases (i.e., positive magnitudes) the peak total water elevations (i.e., the storm tide elevation ζ ST ) in the majority of the PRE. One of the most notable areas is near the top of Shenzhen Bay, where a maximum magnitude of 0.18 m is found, which is largely caused by the phase alteration due to the nonlinear shallow water effects (see details in section 5.2). The above analysis indicates that both the tidally modulated surge generation and phase alteration contribute considerably to the peak overall water elevations; a linear superposition of the atmospheric-only forced pure storm surge (ζ SO ) with the astronomical tide (ζ TO ) can deviate  from the real conditions significantly, as shown in Figure 7, and thus the effects of non-linear tide-surge interactions are vitally important.
As noted in previous studies (Horsburgh and Wilson, 2007;Wolf, 2009;Rego and Li, 2010), the modulation of surge generation and propagation shown above represents the effect of tide on the surge, while the effect of surge on the tide is largely presented as a phase shift of the tidal signal. These mutual influences between the tide and surge contribute to the total effects of tide-surge interaction. Since the residual water elevation ζ TSI , calculated as ζ ST − ζ TO − ζ SO , is the result of tide-surge interaction, it has been taken as a direct measure of the interaction intensity in previous studies (e.g., Bernier and Thompson, 2007;Rego and Li, 2010;Zhang et al., 2010Zhang et al., , 2017. Figure 6 shows that ζ TSI is negligible before and after the passage of Typhoon Hato and that it increases greatly in magnitude as the storm surge develops at all stations. Notable oscillations with a near-tidal period are found in ζ TSI , which is very likely due to the effect of tidal modulation. To quantify the absolute intensity of tide-surge interaction, some studies (e.g. Horsburgh and Wilson, 2007;Rego and Li, 2010;Zhang et al., 2017) have used various different indicators, including the maximum positive (MAX) or minimum negative (MIN) magnitude of ζ TSI , whereas some others (e.g., Bernier and Thompson, 2007;Rego and Li, 2010;Zhang et al., 2010) have used the root-mean-square value (RMS) of ζ TSI 1 , as the representative variable. Evidently, RMS(ζ TSI ) represents the average intensity of tide-surge interaction, while MAX(ζ TSI ) and MIN(ζ TSI ) are more concerned with the maximum intensity that occurs during an entire typhoon event. For Typhoon Hato, the spatial distributions of MAX(ζ TSI ) and RMS(ζ TSI ) in the PRE are shown in Figures 8A,B, respectively. Both of these figures demonstrate the feature that the intensity of tide-surge interaction is strongest in the top of the PRE and Shenzhen Bay and gradually decreases from the estuary/bay head to the estuary/bay entrance, as the bell-shaped geometry can amplify the impact of tide-surge interaction. MAX(ζ TSI ) is about 0.18-0.6 m in the PRE, whereas the magnitude of RMS(ζ TSI ) is much smaller (0.07-0.25 m). The contrast between MAX(ζ TSI ) and RMS(ζ TSI ) indicates that the effect of tide-surge interaction varies strongly over time, which coincides with the distribution pattern of ζ TSI as shown in Figure 6, so that the majority of the energy of the tide-surge interaction is concentrated near the time when the largest storm surge occurs. Besides MAX(ζ TSI ) and RMS(ζ TSI ), a new indicator I r is also plotted in Figure 8C. It is defined as the ratio of RMS(ζ TSI ) to the square root of the product of RMS(ζ SO ) and RMS(ζ TO ), i.e., I r = RMS(ζ TSI )/ RMS(ζ SO ) * RMS(ζ TO ), and is used to reflect the total relative intensity of tide-surge interaction to pure storm surge and pure astronomical tide, similar to that in Zhang et al. (2010). A similar feature is found in Figure 8C as is shown in Figures 8A,B. As the intensity of tide-surge interaction increases in proportion to both surge height and tidal range (Horsburgh and Wilson, 2007), I r is considered to be more appropriate to quantify the relative importance of different physical processes in studying the mechanisms of tide-surge interaction (see details in section 5).
1 The root-mean-square (RMS) of ζ TSI is defined as RMS(ζ TSI ) = T ζ 2 TSI dt/ T, in which T represents the duration of the typhoon event.

The "Subtraction" Approach
To assess the contribution of each nonlinear physical process to the tide-surge interaction, previous studies (Bernier and Thompson, 2007;Zhang et al., 2010Zhang et al., , 2017 conducted numerical experiments using a reduced model approach in which the non-linear terms associated with each physical process were eliminated or linearized: (1) to quantify the nonlinear advective effect (Exp-LAdv), the advective terms were removed from Equations (1) and (2); (2) to quantify the nonlinear bottom friction effect (Exp-LBot), the quadratic form of bottom friction was linearized by using τ bx , τ by = ρ 0 C db (u, v); (3) to quantify the shallow water effect (Exp-LSW), the total water depth D = h + ζ in the governing equations was replaced by h. Therefore, this approach can be regarded as a "subtraction approach, " as it is based on a standard model that includes all three processes and assesses the changes to the interaction intensity after one of the processes is removed. Various aspects of this approach are also briefly summarized in Table 2. Following the same procedure as in the standard experiment (Exp-SD, i.e., the experiment conducted using the complete model including all three processes; section 4), three model runs (i.e., Run-Full, Run-TO, and Run-SO) were conducted in each reduced-model experiment, from which the corresponding residual elevations due to tide-surge interaction (i.e., ζ LAdv TSI , ζ LBot TSI , and ζ LSW TSI ) are calculated. The contribution from each process is then assessed by quantifying the extent to which the intensity of tide-surge interaction is reduced. For this purpose, Zhang et al. (2010) calculated a reduction ratio of RMS(ζ TSI ), i.e., I p , whereas Zhang et al. (2017) closely compared the MAX(ζ TSI ) value calculated by the reduced experiments with that obtained from the standard experiment. I p is defined as follows (Zhang et al., 2010): where RMS(ζ SD TSI ) and RMS(ζ * TSI ) are root-mean-square of ζ TSI obtained from the standard experiment and reduced experiments, respectively, and * represents either LAdv, LBot, or LSW.
Although the contribution from each process can be discerned on close comparisons of the interaction intensity between the results from a reduced model and the standard model as in Zhang et al. (2017), it is best visualized through detailed analysis of the differences obtained by subtracting the interaction intensity of a reduced model from that of the standard model. The reduction ratio R o , based on a generalized form of I p in Equation (14), is employed to quantify the reduction of tide-surge interaction intensity as follows: where P is a general indicator used to represent the intensity of tide-surge interaction, e.g., MAX(ζ TSI ), RMS(ζ TSI ), or I r ; SD   (1) and (2) To assess the non-linear advective effect Exp-LBot Linearize bottom friction by using τ bx , τ by = ρ 0 C db (u, v) To assess the non-linear bottom friction effect Exp-LSW Replace the total water depth D with h To assess the shallow water effect represents the standard experiment, and * represents either LAdv, LBot, or LSW. The calculated R o over the PRE is shown in Figures 9A-I for the reduced experiments Exp-LAdv, Exp-LBot, and Exp-LSW, respectively. All three indicators, MAX(ζ TSI ), RMS(ζ TSI ), and I r , are used to represent the intensity of tide-surge interaction and to calculate the corresponding R o . In the present approach, the contribution from each physical process is expected to lead to a nonnegative reduction ratio (R o ), with its magnitude indicating the strength of contribution. However, negative values of R o are found in all three reduced experiments based on all three intensity indicators (RMS, I r , and MAX) in Figures 9A-I. This common feature suggests that it is more likely that the "subtraction" approach is the reason for the negative reduction ratio rather than that inappropriate indicators are being used. Similar results were also observed in several previous studies (e.g., Tang et al., 1996;Zhang et al., 2017). As explained in Zhang et al. (2017), this phenomenon is due to the "rebalance" effect: in each of the three reduced experiments, when one physical process is removed, the remaining two processes will increase their strength to rebalance the governing equations; a larger intensity of tidesurge interaction induced by these two processes is thus obtained, which leads to a negative R o . Furthermore, the change in the strengths of the remaining two processes (e.g., the nonlinear bottom friction effect and shallow water effect) indicates that the tide-surge interaction intensities induced by these two processes from a reduced model (i.e., P * ) are different from those included in the standard model (i.e., those included in P SD ). Even if the value of P SD −P * is positive, it may not be correct for the intensity induced by the first process (e.g., the nonlinear advective effect). This means that in addition to negative R o values being produced, positive R o values can be influenced by the "rebalance" effect. The R o shown in Figure 9, whether positive or negative, thus cannot correctly represent the contribution from one non-linear process properly. An "addition" approach is therefore developed in the next section in order to improve the analysis.

The "Addition" Approach
Due to the defects found in the above "subtraction" numerical approach, a new method is proposed in this section in order to clearly separate the contributions of the three physical processes and quantify their relative contributions to the tide-surge interaction. As introduced in section 5.1, the "subtraction" approach quantifies the contribution of one specific process to the tide-surge interaction by removing/linearizing its corresponding momentum terms from the standard model. After this operation, each reduced model still contains two of three non-linear effects. In contrast, the present new approach takes an "addition" approach ( Table 3): (a) firstly, a base experiment (Exp-None) was conducted using a reduced model with all three non-linear effects removed; (b) three experiments (Exp-AAdv, Exp-ABot, Exp-ASW) were then carried out, each only taking one non-linear effect into account; (c) following the same method as in the standard experiment and the "subtraction" approach, the astronomical tide (ζ TO ), surge (ζ SO ), and tide-surge interaction residual  (1) and (2) To assess the non-linear advective effect Exp-ABot Exp-None + quadratic bottom friction To assess the non-linear bottom friction effect Exp-ASW Exp-None + D = ζ + h To assess the shallow water effect (ζ TSI ) values corresponding to the above four experiments were obtained.
To assess the quantification of the contribution from each physical process to the tide-surge interaction, a new ratio R n is defined as follows: where P is the general indicator used to represent the intensity of tide-surge interaction as used in Equation (15), SD represents the standard experiment, and ** represent either AAdv, ABot, or ASW. It should be noted that although the ζ TSI obtained from the base experiment (Exp-None) should theoretically be zero as all three nonlinear physical processes are removed, it, in fact, has a magnitude of O(mm) due to the existence of numerical errors. As only one process is included in a specific reduced model, the interaction intensity induced by this process will not be affected by the other two processes. Figure 10 shows the R n values calculated from the reduced experiments Exp-AAdv (Figures 10A-C), Exp-ABot (Figures 10D-F), and Exp-ASW (Figures 10G-I), respectively, by using all of the three representative intensity indicators. As expected, positive R n values were obtained in all cases. For the same reduced experiment, the spatial distribution pattern of R n calculated from RMS(ζ TSI ) is very close to that calculated from I r , indicating that these two indicators of interaction intensity, RMS(ζ TSI ) and I r , provide similar quantification of the relative contributions from the physical processes. However, the spatial distribution of R n values calculated from RMS(ζ TSI ) (or I r ) and that from MAX(ζ TSI ) are very different. This can be explained as follows: both RMS(ζ TSI ) and I r represent the average intensity, whereas MAX(ζ TSI ) represents the maximum intensity of tidesurge interaction that occurs during an entire typhoon event. The magnitudes of R n calculated from RMS(ζ TSI ) and I r also differ from each other, indicating that the pure storm surge levels (ζ SO ) and pure astronomical tide elevations (ζ TO ) in the reduced experiments are not the same as those in the standard experiment. As the tide-surge interaction increases in direct proportion to both surge height and tidal range, a larger RMS(ζ TSI ) or MAX(ζ TSI ) in the reduced experiment may be due to the larger surge height and/or the larger tidal range but not necessarily due to the corresponding nonlinear physical processes themselves. Therefore, it is not appropriate to use RMS(ζ TSI ) or MAX(ζ TSI ) to represent the contributions from the three physical processes to the tide-surge interaction. In contrast, the ratio I r , as shown in Equation (16), reflects the total relative intensity of tide-surge interaction to the pure storm surge and pure astronomical tide, thus eliminating the influences of the changes in surge height and tidal range on the interaction intensity. It is therefore more reasonable to use I r rather than RMS(ζ TSI ) or MAX(ζ TSI ) to quantify the relative contribution from the three physical processes.
The results in Figures 10D-F also show a common feature that the calculated R n in some areas of the PRE is larger than 100%, indicating that the intensity of tide-surge interaction due to one of those processes alone is larger than that obtained from the standard model, in which all three are included. This is a very interesting result, which suggests that certain interactions must have taken place between those three processes and that, for some areas in the PRE, the result of this interaction is to reduce the magnitude of the contribution from the individual processes. In addition, this phenomenon may also be one of the reasons that the "rebalance effect" described in section 5.1 occurs: when one of the three physical processes is removed from the standard model, the remaining processes in the reduced model still interact with each other in a somewhat different way; the "rebalance effect" thus occurs. This finding therefore further indicates that the "addition" approach is a better choice to avoid complication in the quantification of the tide-surge interaction.
From Figure 10, the relative contributions from the three processes to the tide-surge interaction can be directly compared based on the magnitude of R n obtained from the three reduced experiments in specific regions in the PRE. For instance, the results demonstrate that the quadratic bottom friction is most significant in the majority of the PRE, whereas, at the top of Shenzhen Bay, the shallow water effect is more significant due to the limited water depth over the tidal flat. To get a clear overview of the overall contribution from the three processes in the whole PRE, the R n values obtained from the three reduced models are firstly compared with each other and then sorted at every model grid according to their magnitudes. Subsequently, based on the RMS(ζ TSI ) indicator, the process with the largest R n value at each grid node is plotted using its specific color code in Figure 11A. Similarly, the process with the second R n is presented in Figure 11B and with the smallest R n in Figure 11C. Taking the top of Shenzhen Bay as an example, Figure 11A shows that the most important non-linear process there is the shallow water effect (represented in blue); the second most important nonlinear process, shown in Figure 11B, is the quadratic bottom friction (represented in green); and the third most important nonlinear process, shown in Figure 11C, is the non-linear advective effect (represented in red). In a similar way, the processes with the FIGURE 10 | Similar as Figure 9 but for R n calculated by the "addition" numerical approach as described in section 5.2. (A-C) Show R n values calculated by using RMS(ζ TSI ), I r , and MAX(ζ TSI ), respectively, from the reduced experiment Exp-AAdv; (D-F) show the R n values calculated by using RMS(ζ TSI ), I r , and MAX(ζ TSI ), respectively, from the reduced experiment Exp-ABot; (G-I) show the R n values calculated by using RMS(ζ TSI ), I r , and MAX(ζ TSI ), respectively, from the reduced experiment Exp-ASW. largest R n , second R n , and smallest R n based on the I r indicator are shown in Figures 11D-F, and those based on the MAX(ζ TSI ) indicator are shown in Figures 11G-I. The processes with the largest, second, and smallest contribution at any location can thus be directly identified from the corresponding color code. Meanwhile, the area covered by one specific color represents the overall relative importance in the PRE overall. Clearly, no matter which intensity index is used, the results demonstrate a common conclusion that among all the largest-contribution figures (Figures 11A,D,G), the quadratic bottom friction occupies the largest area, which means that the bottom friction contributes the most to the tide-surge interaction. In the second-contribution figures (Figures 11B,E,H), the shallow water effect is clearly the most significant, and hence it contributes to the tidesurge interaction at the second level. Non-linear advection is the third most significant contributor over the majority of the area of the PRE, as shown in Figures 11C,F,I. Similar to those shown in Figure 10, the results obtained from I r are close to those from RMS(ζ TSI ) but are different from MAX(ζ TSI ) at certain locations. For example, at the top of Shenzhen Bay, Figure 11D demonstrates that the shallow water effect dominates, whereas Figure 11G shows the quadratic bottom friction as more important. Due to its shallow water depth, this area is expected to be more significantly affected by the shallow water effect. Therefore, as demonstrated above, the I r value expressed in Equation (16) is recommended for use for the quantification of the contributions from any particular process.

SUMMARY AND CONCLUSIONS
In this study, the characteristics and mechanisms of tide-surge interaction in the Pearl River Estuary during Typhoon Hato are studied in detail by using a 3D ocean model (Zheng P. et al., 2017). Along with the use of a blended atmospheric forcing, which merged the Holland parametric model results with the CFSR reanalysis data, the model reproduces the pure astronomical tides and total sea levels reasonably well, especially at Guanchong and Xipaotai, where the distinctive "double-peak" pattern observed in the measured water levels is well-reproduced by the present model.
To study the characteristics of tide-surge interaction in the PRE, three types of model runs were conducted, from which the total water level (storm tide ζ ST ), the pure storm surge (ζ SO ), the pure astronomical tide level (ζ TO ), the practical storm surge (ζ PS ), and the residual elevation due to the tide-surge interaction (ζ TSI ) were obtained. These results show that, due to the tidesurge interaction, the storm surge is clearly modulated by the tide level, e.g., the magnitudes of ζ PS are greater near low tide but smaller near high tide than ζ SO . The timing of the surge is also altered due to the tidal modulation effect, and the peaks of the predicted ζ PS shown in Figure 6 arrive a bit earlier than those of ζ SO . The horizontal distributions of the differences between ζ max PS and ζ max SO (and the differences between ζ max ST and [ζ SO + ζ TO ] max ) show that the effect of tide-surge interaction can either enhance or reduce the peak surge elevations. In addition, the resultant phase alteration can also affect the peak total water elevation (ζ ST ). A close examination of Figure 7C indicates that the phase alteration significantly increases the peak ζ ST in the majority of the PRE. One of the most notable areas affected by such a process is near the top of Shenzhen Bay, where a maximum magnitude of 0.18 m is found. Three indicators were used to quantify the absolute intensity of tide-surge interaction, namely the previously used MAX(ζ TSI ) and RMS(ζ TSI ) and a newly proposed I r , which reflects the total intensity of tide-surge interaction relative to pure storm surge and pure astronomical tide. As I r eliminates the dependence of the interaction intensity on the magnitude of surge height and tidal range, it is considered more appropriate for use in quantifying the relative importance of different physical processes to tide-surge interaction.
A widely used "subtraction" approach and a newly proposed "addition" approach are adopted to separate the contributions of three non-linear processes to tide-surge interaction and to quantify their relative significance, respectively. In the widely used "subtraction" approach, each of the three processes is removed or linearized from a standard model that includes all processes. The contribution from each specific process to the tide-surge interaction is quantified based on the reduction ratio (R o ) of interaction intensity. However, the results show that the R o value from the "subtraction" approach is greatly affected by the "rebalance" effect ( Figure 9); Thus, it can not correctly represent the significance of its corresponding nonlinear process. An "addition" approach is therefore proposed in which one of the three processes is added onto the baseline simulation without any non-linear effects. A new general ratio R n is defined to quantify the contribution of each process, the value of which can be calculated from any one of three representative indicators of tide-surge interaction intensity. The comparison of the magnitudes of R n between those obtained from the three reduced experiments clearly shows that the quadratic bottom friction, shallow water effect, and nonlinear advective effect have the largest, second-largest, and third-largest contributions to the tide-surge interactions in the majority of the PRE, respectively. Among the three indicators that have been used to represent the intensity of tide-surge interaction, I r is suggested to be more reasonable for use to quantify the relative importance of the three non-linear effects.
Taking Typhoon Hato as a case study, the present research reveals the detailed characteristics of tide-surge interaction in the PRE. The present results are able to provide valuable information for the coastal defense management of different regions inside the PRE, although studies on more typhoon events may be needed. Furthermore, the mechanism of tidesurge interaction is examined by using a newly proposed "addition approach." This new approach is free of the problems caused by the "rebalance" effect and thus is recommended for use in future similar studies and in other regions of the world.

DATA AVAILABILITY STATEMENT
To enable use of the data from the present study by the readers and to achieve further understanding of tide-surge interaction processes, a database hosted by the University of Liverpool (UK) is available at the website: http://datacat.liverpool.ac.uk/id/ eprint/1070.

AUTHOR CONTRIBUTIONS
PZ took part in all activities in this research, including writing the first draft of the paper, model setup and experimental design, and model results analysis. He also proposed the new addition approach and the use of the I r indicator to represent the tidesurge interaction after some inspiring discussions with ML. ML modified the first draft of the paper throughout, helped PZ to design the model experiments, and also gave advice in drafting the final paper. CW designed the approach used to reconstruct the Typhoon Hato wind field by merging the CFSR winds and Holland parametric model results. She also provided the observed wind data at six stations and made a preliminary comparison between the reconstructed and observed wind at six stations. JW and XC helped PZ to analyze the model results, they also made some very useful modifications to the draft paper, which greatly improved the readability and completeness of the draft. MD and PY helped PZ to create the model mesh grid and to determine the basic model configuration. JW and ZH are the project co-PIs, they provided critical data and revision advice on the structure of the manuscript. All the authors participated in the revision of the manuscript.