A novel approach for the assessment of morphological evolution based on observed water levels in tide-dominated estuaries

. Assessing the impacts of both natural (e.g., tidal forcing from the ocean) and human-induced changes (e.g., dredging for navigation, land reclamation) on estuarine morphology is particularly important for the protection and management of the estuarine environment. In this study, a novel analytical approach is proposed for the assessment of estuarine morphological evolution in terms of tidally averaged depth on the basis of the observed water levels along the estuary. The 5 key lies in deriving a relationship between wave celerity and tidal damping or ampliﬁcation. For given observed water levels at two gauging stations, it is possible to have a ﬁrst estimation of both wave celerity (distance divided by tidal travelling time) and tidal damping or ampliﬁcation rate (tidal range difference divided by distance), which can then be used to predict the morphological changes via an inverse analytical model for tidal hydrodynamics. The proposed method is applied to the 10 Lingdingyang Bay of the Pearl River Estuary, located on the southern coast of China, to analyse the historical development of the tidal hydrodynamics and morphological evolution. The analytical results show surprisingly good correspondence with observed water depth and volume in this


Introduction
Understanding the morphological changes in estuaries due to natural processes and human interventions is especially important with regard to sustainable water management and ecological impacts on the estuarine environment (Prandle, 2004;Schuttelaars et al., 2013;Winterwerp and Wang, 2013;20 Wang et al., 2015;Hoitink et al., 2017;Luan et al., 2017). In recent decades, human activities have been crucial factors affecting the morphological evolution in estuaries, and morphological responses of estuaries to human activities have attracted worldwide attention (Syvitski et al., 2009;Monge-Ganuzas et al., 2013;Du et al., 2016;van Maren et al., 2016;Jiang et al., 2017). In previous studies, aerial photographs, satellite images, topographic maps and cross-sectional profiles are the com-25 monly used field data to explore the morphological adjustments with varied resolution and accuracy by means of either data-driven studies (e.g., Zhang et al., 2015;Wu et al., 2014Wu et al., , 2016a or processbased morphodynamic models (e.g., Guo et al., 2014Guo et al., , 2016Luan et al., 2017). It is noted that these high-resolution and successive field data are rare or unavailable in most estuaries, especially because the determination of detailed morphological changes is usually time and energy consuming. Hence, 30 the development of a simple yet effective approach for the assessment of decadal morphological evolution without the need for bathymetric data collection is desirable, which is especially the case for ungauged basins.
In a frictionless, prismatic channel of a constant cross-section (or an ideal estuary with no damping or amplification), the classical formula for tidal wave propagation can be described by (e.g., Savenije,35 2012): where c 0 is the classical wave speed, g is the acceleration due to gravity, h is the tidally averaged and cross-sectionally averaged water depth, and r S is the storage width ratio accounting for the intertidal storage in tidal flats or marshes. Thus for a frictionless, prismatic channel or an ideal estu-40 ary, it is possible to derive an analytical expression for the stream depth, i.e., h = r S c 2 0 /g, which is proportional to the square of the classical wave speed (Savenije and Veling, 2005). However, Equation (1) does not apply in real alluvial estuaries due to the predominance of either morphological convergence or bottom friction. Thus, to have a first-order estimate of the stream depth, a relationship taking into account both tidal damping (or amplification) and channel convergence should be 45 adopted. Such a relationship can be easily obtained from the recently published analytical solutions for tidal hydrodynamics in convergent estuaries (e.g., Savenije and Veling, 2005;Savenije et al., 2008;Toffolon and Savenije, 2011;van Rijn, 2011;Cai et al., 2012Cai et al., , 2016Winterwerp and Wang, 2013).
Previous studies have clearly demonstrated that the behaviour of both the tidal damp-50 ing/amplification and wave speed strongly depends on the degree of channel convergence and on the intensity of friction (e.g., Jay, 1991;Savenije and Veling, 2005). In particular, analytical solutions of tidal hydrodynamics are able to predict the main trends of tidal damping/amplification and wave speed for given forcing tidal amplitude, geometry (e.g., the spatially averaged depth and width convergence length) and friction (e.g., Prandle and Rahman, 1980;Friedrichs and Aubrey, 1994; 55 Savenije et al., 2008;Toffolon and Savenije, 2011;van Rijn, 2011;Cai et al., 2012Cai et al., , 2016Winterwerp and Wang, 2013). Conversely, if tidal damping/amplification and wave speed are known from observations, analytical models can be used to estimate the unknown tidally averaged depth and other tidal dynamics (e.g., velocity amplitude, phase difference between velocity and elevation, etc.). Following this line of thought, the main objective of this study is to provide an analytical framework 60 (for an infinitely long channel) to assess the morphological evolution of estuaries based on the tidally averaged depth derived from the observed tidal hydrodynamics (i.e., tidal damping/amplification and wave speed).
It is worth noting that considerable simplifications of geometry and forcing are made in order to derive generic analytical solutions of the governing St. Venant equations. Here, we restrict consider-65 ation to tide-dominated estuaries with small tidal amplitude to depth ratios and Froude numbers. The fundamental assumption is that the geometry (cross-sectional area, width and depth) of the channel can be described by a simple exponential function. We also assume that the velocity of river discharge is small compared to the tidal velocity. In this paper, we use the analytical relation of wave celerity, first proposed by Savenije and Veling (2005), to predict the estuarine morphological dynam-70 ics in terms of the large-scale geometrical properties of the estuary, i.e., the tidally averaged depth and volume based on the observed water levels through the estuary.
The Pearl River is the second largest river in terms of water discharge in China (Liu et al., 2017) and drains into the South China Sea through Lingdingyang Bay, Modaomen Estuary and Huangmaohai Bay. Among these three outer bays, Modaomen estuary carries the largest portion of freshwater 75 discharge (Gong and Shen, 2011), while Lingdingyang Bay is the largest of these tidal systems in terms of water area and is constituted by three shoals and two troughs. Previous studies have examined the morphological evolution of Lingdingyang Bay using satellite images (e.g., Wu et al., 2014), topographic maps (e.g., Wu et al., 2016a,b) and cross-sectional profiles (Wu et al., 2016b). Due to intensive human activities (e.g., land reclamation and channel dredging), significant morphological 80 changes have occurred in Lingdingyang Bay in recent decades. As a consequence, the West and East shoals expanded towards the bay, and the two troughs became narrower and deeper (Zhang et al., 2015;Wu et al., 2016b). Furthermore, these morphological changes have resulted in changes in tidal dynamics (Zhang et al., 2010;Deng and Bao, 2011). Our study can provide a new method to explore the morphological evolution, which provides scientific guidelines for management and regulation 85 projects in estuaries.
The paper is organized as follows. Section 2 describes the analytical method for reproducing the tidal hydrodynamics and an inverse model for predicting the tidally averaged depth. An overview of the study area is provided in section 3. In section 4, the morphological evolution of Lingdingyang Bay is described, and the proposed method is applied. Subsequently, the relationship between mor-90 phology and tidal dynamics are discussed, along with the model limitations (section 5). Finally, some conclusions are drawn in section 6.
2 Methodology and material 2.1 Analytical model for tidal hydrodynamics in an infinitely long channel Neglecting the nonlinear continuity term ∂h/∂x and the advective term U ∂U/∂x, the linearized 95 depth-averaged equations for the conservation of mass and momentum in a channel with a gradually varying cross-section can be described by (e.g., van Rijn, 2011):

100
where U is the cross-sectionally averaged velocity, Z is the free surface elevation, h is the water depth, B is the tidally averaged width (hereafter, overbars denote tidal averages), t is the time, x is the longitudinal coordinate measured positive in landward direction (x=0 at the mouth), and r is the linearized friction factor defined by (Lorentz, 1926):

105
where the coefficient 8/(3π) stems from adopting Lorentz's linearization of the quadratic friction term (Lorentz, 1926) considering only one single predominant tidal constituent (e.g., M 2 ), and K is the Manning-Strickler friction coefficient.
To derive the analytical solution for the tidal hydrodynamics in convergent channels, it is assumed that the tidally averaged cross-sectional area A and width B can be described by the following 110 exponential functions: where A 0 and B 0 are the respective values at the estuary mouth and a and b are the convergence 115 lengths of the cross-sectional area and width, respectively. The other fundamental assumption is that the flow is mainly concentrated in a rectangular cross-section, with a possible influence from storage areas described by the storage width ratio r S . It directly follows from the assumption of a rectangular cross-section that the tidally averaged depth is given by h = A/B.
Considering a system that is forced by a sinusoidal tidal wave with period T and frequency ω = 120 2π/T , the wave functions of water level and flow velocity are defined as: where η and υ are the wave amplitudes of elevation and velocity, while ϕ Z and ϕ U are their corre-125 sponding phases. As the wave propagates into the estuary, it has a wave celerity c and a phase lag ε, defined as the phase difference between high water (HW) and high water slack (HWS) (or between low water (LW) and low water slack (LWS)). For a simple harmonic wave, ε = π/2 − (ϕ Z − ϕ U ).
After scaling Equations (2) and (3) where f is the dimensionless friction factor, and ζ is the dimensionless tidal amplitude defined as: 150 Making use of these dimensionless parameters, Cai et al. (2012), building on the previous work by Savenije et al. (2008), showed that the analytical solution of the tidal hydraulic equations can be obtained by solving a set of four implicit equations, i.e., the phase lag Equation (16), the scaling Equation (17), the damping Equation (18) and the celerity Equation (19): In addition, the solutions for the phases of elevation and velocity are given by (see details in Cai et al., 2016): where ℜ and ℑ are the real and imaginary parts of the corresponding term and A * and V * are 165 complex functions of amplitudes that vary along the dimensionless coordinate x * = x/(c 0 T ) as follows: where i is the imaginary unit, and Λ is a complex variable defined as: 170 It is worth considering the solution in the special case of an ideal estuary (with no damping or amplification δ = 0) since natural estuarine systems tend to reach an equilibrium state of ideal estuaries. By imposing δ = 0 in the set of Equations (16)-(19), one can easily obtain: (23) Figure 1 shows the variation of the four dimensionless parameters as a function of γ and χ ob- It is important to note that the analytical solutions of the linear model are local because they depend only on local (fixed position) quantities (i.e., the local tidal amplitude to depth ratio ζ, the local estuary shape number γ and the local friction number χ). To follow along-channel variations 180 of these local variables, a multi-reach approach (subdividing the whole estuary into short reaches) has been used, in which the damping number δ is integrated in short reaches over which the estuary shape number γ and the friction number χ are assumed to be constant. This is done by using simple explicit integration of the linear differential equation (Savenije et al., 2008):

185
where η 0 is the tidal amplitude at the origin of the axis for every short reach, while η 1 the tidal amplitude at a distance ∆x (e.g., 1 km) upstream.

Estimation of the tidally averaged depth
To address directly the correspondence between tidal dynamics and morphology, the celerity Equation (19) proposed by Savenije and Veling (2005) is considered. This analytical relationship is an 190 extension of the classical celerity equation for a progressive wave in a frictionless prismatic channel and depicts the along-channel wave celerity (λ = c 0 /c) as a function of the tidal damping (δ) and estuary shape number (γ), hence accounting implicitly for the water depth (see Equation (9)). Figure   2 shows the analytically computed λ 2 over a wide range of values of δ, from -3 to 1, and γ, from 0 to 4, directly from Equation (19). It can be clearly seen from Figure  λ < 1, since the actual wave celerity c is larger than the classical one c 0 , while it is the opposite (λ > 1) for negative values of δ (damped estuary). The underlying mechanism lies in the imbalance between channel convergence and bottom friction (i.e., stronger impact of convergence than friction in the case of amplification, and the opposite in the case of damping). Without damping nor amplifi-200 cation (λ = 1, indicating balance between convergence and friction), then the estuary corresponds to an ideal estuary, where the wave celerity equals c 0 . Interestingly, according to Equation (19), for an amplified estuary (δ > 0), the actual wave celerity c could equal to c 0 for the special case of γ = δ, which is due to the balance between channel convergence and acceleration effects (Jay, 1991).
Estimation of the tidally averaged depth h can be derived analytically by rewriting Equation (19) 205 in terms of the width convergence length b, wave celerity c and tidal amplitude damping rate δ H , which leads to (see details in Appendix B): Figure 3 illustrates the analytically computed depth with Equation (25) as a function of the wave celerity c and tidal amplitude damping rate δ H for a constant width convergence length b=200 km 210 and storage width ratio r S =1. It can be seen that the depth h tends to increase with the celerity c.
On the other hand, the depth h is decreased with the tidal damping rate δ H until a minimum value is reached at a critical δ H corresponding to the condition ∂h/∂δ H = 0. A further increase of the δ H yields an increase of h. The wave celerity c and the tidal amplitude damping rate δ H can be estimated for a reach of ∆x: where c HW and c LW are the wave celerity for high water (HW) and low water (LW), ∆t HW and ∆t HW are the corresponding travelling times over the reach, η 1 is the tidal amplitude in the seaward 220 part, and η 2 is the tidal amplitude ∆x upstream. The wave celerity c can also be estimated through harmonic analysis. In this case, Equation (26) can be rewritten as where ϕ Z1 is the phase of elevation at the seaward station, while ϕ Z2 is the corresponding values at ∆x upstream. As illustrated in the next section, these parameters can be easily obtained from tidal 225 gauge observations, providing a simple means to quickly assess the mean water depth with Equation (25).

Overview of Lingdingyang Bay
The Pearl River ( Figure  and four western outlets (Modaomen, Jitimen, Yamen and Hutiaomen). Lingdingyang Bay is the largest estuary of the Pearl River Estuary, and mainly receives water and sediment inputs from the four eastern outlets (Figure 4b). Lingdingyang Bay is a funnel-shaped subaqueous delta and has a 235 complicated geomorphology pattern with two deep channels (i.e., East and West Channels) between three shoals (i.e., East, Middle and West Shoals) (see Figure 4c). The tide in Lingdingyang Bay mainly come from the Pacific Ocean and has an irregular and semidiurnal character, with a mean tidal range between 1.0 and 1.7 m (Mao et al., 2004). Generally, tidal propagation in Lingdingyang Bay is affected by geometry (e.g., bank convergence), bottom topography and river discharge, among 240 which the impact from river discharge is minor because the freshwater discharge is small compared to the amplitude of the tidal discharge. Hence, Lingdingyang Bay is a tide-dominated estuary owing to its dramatic volume for water bodies. In response to the typical funnel shape topography, the mean tidal range is considerably amplified when a tidal wave propagates along Lingdingyang Bay, increasing from 1.1 m at Chiwan (denoted by CW hereafter) near the mouth to approximately 3.2 245 m at Sishengwei (denoted by SSW hereafter) near the head of the bay, located 58 km upstream.
The tidal wave is continuously amplified until Huangpu gauging station (24 km upstream from CW) with a mean tidal range of approximately 3.6 m, after which it is damped until vanishing .
In recent decades, intensive human activities (e.g., land reclamation, channel dredging, sand ex- (relative to the datum of lowest low water level, which is 1.7 m below mean sea level), respectively (Li, 2008). Furthermore, sand excavation is conducted in Lingdingyang Bay. According to Wu et al.
The water depths, isobaths and estuarine outlines on these maps were digitalized in order to generate 270 a digital elevation model (DEM) used to analyse the morphological changes of Lingdingyang Bay over ∼50 years (from 1965 to 2015). Using the latitude/longitude information, they were then projected to UTM-WGS84 coordinates of China and interpolated to a 50 × 50 m grid DEM with kriging interpolation in ArcGIS, which was extensively applied to analyze the evolution of morphological changes (e.g., Brunier et al., 2014;Liu et al., 2019).

275
To explore the tidal hydrodynamics in Lingdingyang Bay, tidal water level records (see Figure   5 for the records in January, the records in other months were not presented) for the periods with surveyed bathymetric maps were obtained from the Ministry of Water Resources of China (MWRC).
These observations were taken at the CW and SSW gauging stations (see their locations in Figure   4) and were used to have a first estimate of wave celerity (26) and the tidal damping/amplification 280 rate (27) over the studied channel. From Figure 5, we observe that the tide in Lingdingyang Bay is characterized by a semi-diurnal mixed tidal regime with apparent daily inequality in range and time between high and low waters. As the tide propagates from CW to SSW, the tidal range is amplified due to the strong channel convergence. The collected tidal water level records generally contained two high and two low water levels for each day. These data were then interpolated to one hour 285 intervals for harmonic analysis by using the shape preserving piecewise cubic interpolation. The tidal water levels derived from such an interpolation method can well retain the power spectra of low frequency bands and principal tides (e.g., M 2 ), while the high frequency bands may not be entirely  to 1974, an area of 36.3 km 2 was reclaimed within the bay (see Table 1), extending the coastline southward towards the open sea (see Figure 6A1, 6A2 and 6B1) and reducing the surface water area by 36.3 km 2 ( Table 1). The magnitude of bathymetric changes was mainly less than 1 m in the bay ( Figure B1), with the general dominance of deposition. The mean water depth only decreased by 0.09 m (Table 1), corresponding to a water volume decrease of 2.0×10 8 m 3 during this period (Table   295 1). From 1974 to 1989, land reclamation was the most intensive, resulting in an increase of land area (and decrease of water area) of 129.2 km 2 . During this period, significant deepening, ranging from 1 to 5 m, occurred in the West Trough due to the dredging of the navigation channel. The mean water depth only increased by 0.08 m in Lingdingyang Bay. However, the water volume continued to decrease by 0.3×10 8 m 3 in relation to the decrease in water area. Since the 1990s, dredging for the 300 maintenance of the navigation channel intensified. From 1998 to 2015, land reclamation continued to increase the land area and decrease the water area in Lingdingyang Bay (Table 1); however, the mean water depth and water volume increased by 0.75 m and 5.1×10 8 m 3 , respectively ( and pits up to ∼20 m deep occurred in its upper sector due to sand excavation ( Figure 6B5). It is evident that the morphological evolution in Lingdingyang Bay has been significantly altered by human activities since the 1990s. In this paper, the observed wave celerity is derived from the travelling time of both high and low water levels (26), while the tidal amplitude damping/amplification rate is computed according to Equation (27). To illustrate the spring-neap variabilities of tidal dynamics, Figure 7 shows the observed c and δ H as a function of the tidal amplitude at the estuary mouth η 0 (i.e., CW station) for different years. In general, we observe that both wave celerity and the damping/amplification rate 315 are decreased with increasing tidal amplitude at the estuary mouth with slightly different negative slopes (indicated by α for the wave celerity and β for tidal damping/amplification rate), suggesting a more strongly amplified yet faster wave for the neap tide than that for the spring tide. Furthermore, we observe a clear pattern of increasing average wave celerity over the period between 1965 and 2015 (Figures 7a-f), which corresponds to an increasing trend of the tidal amplification rate over the 320 period. In addition, the spring-neap variability of tidal amplification is apparently large in 1965 and decreases until 2015.
Note that both the wave celerity and the tidal damping/amplification rate reflect the imbalance between channel convergence and bottom friction (Savenije and Veling, 2005). It is evident from the definition of the friction number χ (10) that effective bottom friction experienced for the spring 325 tide is stronger than that for the neap tide due to a larger tidal amplitude to depth ratio ζ during the spring tide. On the other hand, the estuary shape number and the tidally averaged depth h are more or less the same for spring and neap tides. Hence, the tidal damping/amplification rate experienced by the estuary is larger for the neap tide than that for the spring tide.
To investigate the underlying mechanism of such a spring-neap variability of wave celerity, we 330 further rewrite the celerity Equation (19) by substituting Equation (12): It is evident from Equation (29) that the wave celerity c is increased with the tidal damping/amplification number δ for given constant values of γ and h. Hence, assuming a more or less constant tidally averaged depth over a spring-neap cycle, it is concluded that the wave celerity during 335 the neap tide is generally faster than that during the spring tide.

Performance of analytical model for reproducing the tidal hydrodynamics
Before an inverse analytical model can be used to predict the morphological changes in estuaries, it is required to well calibrate and validate the model against observations. Hence, the analytical model presented in section 2.1 is used to reproduce the historic physical properties of the tidal 340 wave (i.e., tidal damping/amplification rate and wave celerity) in Lingdingyang Bay. The estuarine system is subject to a harmonic tide at the estuary mouth (i.e., CW station). In order to calibrate and validate the analytical model, the tidal properties (including the tidal amplitude and phase) of predominant M 2 tide at a monthly averaged scale were extracted by means of a classical harmonic analysis making use of T_TIDE MATLAB Toolbox provided by Pawlowicz et al. (2002). Based on the collected topographic maps (see details in section 3.2), the geometry of Lingdingyang Bay can be described by exponential functions (5)-(6) (Figure 8), and the fitted geometric characteristics are given in Table 2. In Figure 8, we also observe a sudden decrease of cross-sectional area near the Humen outlet, which is due to the strong width convergence towards the rock-bound gorges with a relatively greater water depth.

350
In Figure 9, the analytically computed tidal amplitude η (24) and phase ϕ Z (20) of the predominant M 2 tide at SSW (x=58 km) are compared with the observed values in Lingdingyang Bay for different years. We calibrated the analytical model by adjusting the Manning-Strickler friction coefficient K and the storage width ratio r S , which are detailed in Table 3. The model performance was evaluated by the root mean squared error (RMSE), where RMSE=0 corresponds to perfect agree- Subsequently, the analytically computed tidal characteristics of Lingdingyang Bay were used to describe how the tidal hydrodynamics are affected by the morphological evolution. In Figure 10, the diagrams of the velocity number µ, damping number δ, celerity number γ and phase lag ε are shown together with the trajectory of their corresponding dimensionless parameters in Lingdingyang Bay.

365
It is worth noting that the longitudinal estuary shape number γ decreased in 1965, 1974 and 1989 due to the shallowing of the estuary in the landward direction. Conversely, we observe an increased γ in 2009 and 2015, which is due to the increased depth in the channel (see Figure 8). In these graphs,

Estimation of the tidally averaged depth
The successful reproduction of tidal hydrodynamics using an analytical model suggests a close relationship between tidal damping/amplification and wave celerity, which can be described by the 385 celerity Equation (19). Thus, it is possible to develop an inverse analytical model to estimate the tidally averaged depth h for given observed wave celerity and the tidal damping/amplification rate from observed water levels, as presented in section 2.2. In this case, we assume a constant depth, and hence, the width convergence length equals to the cross-sectional area length, i.e., a = b.
We adopted the width convergence length from topographic maps and estimated wave celerity and 390 tidal damping/amplification from the observed water levels. Combining these parameters with the calibrated storage width ratio r S (see Table 3), Equation (25) can be used for a first estimate of the tidally averaged depth. Figure 11 shows the observed wave celerity and tidal damping/amplification together with the estimated tidally averaged depth and tidal amplitude at the CW station over the studied period. The monthly variation in estimated depth is mainly due to the change in mean sea 395 level. In particular, we see a very similar variability of the derived depth with that of the wave celerity, which suggests a stronger influence of wave celerity when compared with the tidal amplification in Lingdingyang Bay. Figure 12 shows the analytically computed tidally averaged depth and water volume at an annual scale in Lingdingyang Bay in different years. Note that in Figure 12, the annual mean storage suggests that the proposed method can be a useful tool to have a first estimate of the morphological evolution in terms of depth and volume.

Influence of the morphological evolution on tidal hydrodynamics
Lingdingyang Bay is a typical funnel-shaped estuary, where the tidal dynamics are one of the main 410 factors maintaining the stability state of estuarine morphology. Generally, the tides in the Pearl River Estuary are influenced by the geometry (i.e., cross-sectional variation) and river discharge (Zhang et al., 2010). The water discharge shows insignificant change since the 1950s (Liu et al., 2017); therefore, morphological changes became the main factors influencing the tidal dynamics of Lingdingyang Bay. Due to land reclamation, the water area and water volume of Lingdingyang 415 Bay decreased from 1965 to 1974. Since the 1980s, channel dredging was conducted in the West Trough (i.e., Lingdingyang Channel) (Wu et al., 2016a), and the West Trough became deeper from 1974 to 1989 ( Figure 5B2). Therefore, although land area in Lingdingyang Bay increased by 30.5% from 1974 to 1989, water volume only decreased by 0.8%, and the mean water depth increased by 0.08 m (Table 1). In addition, since 1989, the mean water depth has increased significantly 420 because of intensive channel dredging and reduced sediment inputs (due to reservoirs operation in the Pearl River basin) in Lingdingyang Bay. Hence, it seems that the morphological evolution pattern of Lingdingyang Bay has changed since 1989. Correspondingly, temporal changes in the tidal dynamics of Lingdingyang Bay show that the tidal dynamics pattern has also changed since 1989 (see Figures 7, 9,10 and Table 4).

425
To better understand the response of the morphological evolution to tidal dynamics, we rewrite the equilibrium depth h I (subscript I indicates ideal estuary) through Equation (18) for the case of an ideal estuary, where δ = 0: which shows that the equilibrium depth h I is inversely proportional to the Manning-Strickler friction 430 coefficient K with the power of 12/11, while it is proportional to the velocity amplitude υ with the power of 6/11. If we assume a more or less constant width convergence length, then the response of the morphological evolution can be inferred from the combined impacts from the Manning-Strickler friction coefficient (representing the bottom friction) and velocity amplitude (representing the tidal dynamics). Specifically, we observed a decreasing trend of calibrated K during 1965-1989, while 435 there was an increasing trend during 1989-2015 (see Table 3). This suggests that the tidally averaged depth tends to increase before the 1990s, after which it would decrease again. However, we actually observed a consistent increase of water depth over the studied period, which has to do with the alteration of the tidal dynamics, especially due to the velocity amplitude. In Table 4, it can be seen that the dimensionless velocity number µ decreased from 1965 to 1989 (indicating a decrease trend of 440 water depth), after which it increased until 2015 (indicating an increase trend of water depth). Hence, for the period of 1965-1989, the potential impact from bottom friction on morphological evolution dominated over the alteration of tidal dynamics. In contrast, the observed channel deepening in the following decades is mainly controlled by the reinforced tidal dynamics that has dominated over the bottom friction.

Model limitations
It should be noted that several assumptions are made in order to derive the analytical solutions for tidal hydrodynamics. The fundamental assumption is that the tidal amplitude to depth ratio and the Froude number are considerably smaller than unity so that the linearized St. Venant equations can be used for the derivations. A second fundamental assumption is that both the cross-sectional 450 area and width can be described by exponential functions, following Equations (5) and (6), as is the case for most alluvial estuaries. We also assume that the cross section is most rectangular, with a possible impact from lateral storage area described by the storage width ratio r S . In this study, the r S values were obtained by calibration against the extracting harmonic constants for given geometries observed in Lingdingyang Bay. For a fully predictive model, additional satellite maps (such as 455 archived Landsat MSS/TM data, http://glovis.usgs.gov/) at different periods (i.e., the spring tide and moderate tide) are required to help define the intertidal zone (and hence the r S value). In addition, we neglect the influence of river discharge on tidal dynamics, which is not such a restriction in the downstream part of alluvial estuaries where the river flow velocity is small compared to the tidal velocity. It was shown that the tidally averaged depth varies due to the nonlinear residual friction 460 effect (e.g., Cai et al., 2014). By using a harmonic solution, the residual water level slope along the channel is assumed to be null, which means that the model cannot reproduce the spring-neap change of the tidally averaged depth along the channel, although the depth is usually larger during the spring tide (and also during the flood season) than during the neap tide (and also during the dry season) in tide-dominated estuaries. Moreover, by adopting an analytical solution for tidal hydrodynamics 465 in an infinitely long channel, the model excludes the impact of reflected waves from the closed end on tidal dynamics. Therefore, the proposed analytical approach is preferably applied to tidally dominated estuaries (or bays) free of tidal barriers or sluices at the upstream end, and where the tide dominates over the river discharge.

470
In this paper, a novel approach for estimating the tidally averaged depth was proposed to understand the morphological evolution based on the observed water levels at 2 stations (at least) in estuaries.
The linear analytical hydrodynamics model proposed by Cai et al. (2012) was used to explore the hydrodynamics changes in Lingdingyang Bay, where we observe an evident hydrodynamics pattern transition due to the substantial changes of morphology. The analytically computed tidal amplitude 475 and travelling time are compared to the observations, and the correspondence is good, which suggests that the proposed model can be a useful method for analysing the changes of physical tidal properties along the channel. Conversely, for given observed tidal water levels, it is possible to have a first estimate of wave celerity and tidal damping/amplification rate, and subsequently the tidally averaged depth and other tidal properties (e.g., velocity amplitude, phase lag), by manipulating the where there are not sufficient data (e.g., detailed navigational charts) available, while a first estimate of estuary depth and water volume are required.
Data availability. The data and source codes used to reproduce the experiments presented in this paper are available from the authors upon request (caihy7@mail.sysu.edu.cn).
Author contributions. All authors contributed to the design and development of the work. The experiments were originally carried out by HC. PZ, SH carried out the data analysis. FL and HC prepared the paper with contributions from all co-authors. EG, PM and QY reviewed the paper.
Competing interests. The authors declare that they have no conflict of interest.
Acknowledgements. We acknowledge the financial support from the National Key R&D of China (Grant No. 490 2016YFC0402601), from the National Natural Science Foundation of China (Grant No. 51979296, 51709287, 41706088, 41476073), from the Fundamental Research Funds for the Central Universities (No.18lgpy29) and from the Water Resource Science and Technology Innovation Program of Guangdong Province (Grant No. 2016-20, 2016. The work of the second author was supported by FCT research contracts Appendix A

Scaling the governing equations
We introduce a scaling on Equations (2) and (3), similar to that used by Savenije et al. (2008), to derive the dimensionless equations, with the asterisk superscript denoting dimensionless variables:

500
where η 0 and υ 0 are the tidal amplitude and velocity amplitude at the estuary mouth, L is the wavelength and T is the tidal period. Note that the scaling of tidal flow velocity and water level fluctuation are slightly different from the scaling used by Savenije et al. (2008) because they are scaled with the corresponding values at the estuary mouth. For an infinite length estuary, the velocity amplitude and the tidal amplitude are proportional: which implies that the ratio of the velocity amplitude to the tidal amplitude is constant: Making use of the assumption (A3), equations (2) and (3) may then be rewritten as: The real scales of velocity amplitude υ and wavelength L are scaled with the corresponding values for a frictionless tidal wave in a channel with zero convergence (U 0 , L 0 ) as a reference: where we introduce the unknown dimensionless velocity number µ and celerity number λ.
For the case of a frictionless estuary with zero convergence, the classical wave celerity c 0 is defined by (1), while the velocity amplitude U 0 and the wavelength L 0 are: where ζ represents the dimensionless tidal amplitude to depth ratio (see Equation (15)).
Assuming h * = 1 (i.e., h = h) in the continuity Equation (A4), then the dimensionless Equations 525 (A4) and (A5) read: where the dimensionless parameters γ and χ have been introduced as the estuary shape number 530 (9) and the friction number (10), respectively, while the corresponding definitions of the velocity number µ and the celerity number λ are defined in (11) and (12), respectively.

Appendix B
Derivation of the tidally averaged water depth from an inverse analytical model Introducing the dimensionless parameters (9), (12), (13) into the celerity equation (19) yields where a has been replaced by b in the definition of γ (9) due to the assumption of a horizontal bed.    shape number γ and friction number χ obtained with a linear model (Cai et al., 2012), in which thick red lines represent the ideal estuary, where δ = 0, λ = 1, µ = 1/ √ γ 2 + 1 and tan(ε) = 1/γ.