Extension of the general unit hydrograph theory for the spread of salinity in estuaries

From both practical and theoretical perspectives, it is essential to be able to express observed salinity distributions in terms of simplified theoretical models, which enable qualitative assessments to be made in many problems concerning water resource utilization (such as intake of fresh water) in estuaries. In this study, we propose a general and analytical salt intrusion model inspired by Guo's general unit hydrograph theory for flood hydrograph prediction in a watershed. To derive a simple, general and analytical model of salinity distribution, we first make four hypotheses on the longitudinal salinity gradient based on empirical observations; we then derive a general unit hydrograph for the salinity distribution along a partially mixed or well-mixed estuary. The newly developed model can be well calibrated using a minimum of three salinity measurements along the estuary axis and does converge towards zero when the along-estuary distance approaches infinity asymptotically. The theory has been successfully applied to reproduce the salt intrusion in 21 estuaries worldwide, which suggests that the proposed method can be a useful tool for quickly assessing the spread of salinity under a wide range of riverine and tidal conditions and for quantifying the potential impacts of human-induced and natural changes.

1 Introduction

An estuary is the place where the fresh water meets the saline water. It is crucial to quantify the spatiotemporal salinity dynamics determined by the competition between the advective salt flux due to river flow and the dispersive salt flux caused by tidal currents, since it directly affects water quality and the related water resource management in general. It is well known that the key to quantifying the salinity distribution along an estuary is the efficient dispersion coefficient, which incorporates all mixing mechanisms that counteract the advective salt transport and regards the complex estuarine system as a whole. With the one-dimensional steady-state salt balance equation, indicating the equilibrium between the advective and dispersive transports of salt, it is possible to derive an empirical relationship for the salt intrusion in estuaries (Prandle, 1981; Savenije,1986, 1989, 1993, 2005, 2012; Lewis and Uncles, 2003; Gay and O'Donnell, 2007, 2009; Kuijper and Van Rijn, 2011; Cai et al., 2015; Zhang and Savenije, 2017, 2018). Amongst the proposed solutions, the empirical model using Van der Burgh's coefficient (e.g., Van der Burgh, 1972; Savenije, 1986) functions well in a wide range of estuaries worldwide (e.g., Savenije, 2005, 2012). In addition to practical applications, such an empirical model can be very useful from a physical perspective when its theoretical basis is well understood.

Recently, Guo (2022a, b, c) revisited the classical unit hydrograph (UH) theory, which is widely used in hydrology for predicting a flood hydrograph from a known storm in a watershed. Based on three hypotheses on instantaneous UHs derived from observations, he obtained a general and analytical expression which represents the discharge from a continuous excess rainfall occurring at a uniform rate for an indefinite period. This so-called S-hydrograph is expressed in terms of a unit volume of excess rainfall and is used to derive a UH of any storm duration. It appears that the shape of the S-hydrograph is rather similar to the salinity distribution curve along estuaries, while the instantaneous UH curve resembles the longitudinal salinity gradient. This correspondence opens the possibility that the UH method can be applied to describe the spread of salinity in estuaries.

The objective of this study is to derive a general and analytical expression of the salinity distribution and thus to derive the salinity gradient analytically following Guo's UH method (Guo, 2022a, b, c). To this end, we start with a review on Guo's general UH theory, together with the Savenije's empirical salt intrusion model, which is derived from the steady-salt balance equation and performs well against numerous salt measurements along many different estuaries (e.g., Savenije, 2005, 2012). Subsequently, we make four hypotheses based on empirical observations and follow the general UH theory, which leads to a newly developed analytical model for the spread of salinity in estuaries. The model was then applied to real estuaries with a wide range of riverine and tidal conditions. After that, the proposed model was compared with the conventional Savenije's model to discuss the physical foundation of the proposed model, which requires further study in the future.

2 Review of the general unit hydrograph and empirical salt intrusion model

2.1 General unit hydrograph theory

It was shown that a classical instantaneous UH u(t) [T⁻¹] (representing the discharge due to a unit volume input of excess rainfall) with regard to time t [T] should satisfy the following properties (Chow et al., 1988):

$$\begin{array}{}\text{(1)}& {\displaystyle }u\left(t\right)=\mathrm{0}\text{for}t\le \mathrm{0},\text{(2)}& {\displaystyle }\mathrm{0}\le u\left(t\right)\le u_{\text{p}}\text{for}t>\mathrm{0},\text{(3)}& {\displaystyle }u\left(t\right)\to \mathrm{0}\text{for}t\to \mathrm{\infty },\end{array}$$

and

$$\begin{array}{}\text{(4)}& \underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}u\left(t\right)\mathrm{d}t=\mathrm{1},\end{array}$$

where u_p [T⁻¹] in Eq. (2) represents the peak discharge of the instantaneous UH. Equation (4) is the mass conservation equation, indicating that the total outflow volume (i.e., the left side) corresponds to the unit volume input (i.e., the right side). Making use of the definition of the instantaneous UH $u\left(t\right)=\mathrm{d}U/\mathrm{d}t$ (where U represents the dimensionless S-hydrograph), Eq. (4) can be rewritten as follows:

$$\begin{array}{}\text{(5)}& U\left(\mathrm{\infty }\right)=\mathrm{1}.\end{array}$$

Nash (1957) derived an analytical expression for the instantaneous UH that satisfies Eqs. (1)–(4), which is the well-known Nash's gamma function, written as follows:

$$\begin{array}{}\text{(6)}& u\left(t\right)={\displaystyle \frac{\mathrm{1}}{{\mathit{\alpha }}_{\mathrm{1}}\mathrm{\Gamma }\left({\mathit{\alpha }}_{\mathrm{2}}\right)}}{\left({\displaystyle \frac{t}{{\mathit{\alpha }}_{\mathrm{1}}}}\right)}^{{\mathit{\alpha }}_{\mathrm{2}}-\mathrm{1}}\exp \left(-t/{\mathit{\alpha }}_{\mathrm{1}}\right),\end{array}$$

where α₁ and α₂ are model parameters, while Γ(α₂) is the gamma function. Figure 1 illustrates an arbitrary distribution (t=0–30) of the instantaneous UH for given values of ${\mathit{\alpha }}_{\mathrm{1}}={\mathit{\alpha }}_{\mathrm{2}}=\mathrm{3}$. It can be seen from Fig. 1 that the instantaneous UH consists of the following two distinct regions: the rising limb described by the power function in Eq. (6) and the recessing limb described by the exponential function in Eq. (6).

Figure 1 Illustration of the instantaneous UH analytically computed using Eq. (6) for given values of ${\mathit{\alpha }}_{\mathrm{1}}={\mathit{\alpha }}_{\mathrm{2}}=\mathrm{3}$.

Although Eq. (6) is widely used as the analytical solution of an instantaneous UH, it has the following two weaknesses: (1) it has an initial condition of zero, which is not necessarily the case for real instantaneous UHs; and (2) a general and analytical solution of the S-hydrograph does not exist. In order to remove these weaknesses, Guo (2022a) made the following three hypotheses on the instantaneous UH based on empirical observations: (1) the instantaneous UH increases exponentially along the rising limb; (2) the instantaneous UH decreases exponentially along the recessing limb; (3) the instantaneous UH tends towards 0, and the S-hydrograph tends towards 1, as t tends towards infinity. Subsequently, he derived the following general and analytical solution for the S-hydrograph:

$$\begin{array}{}\text{(7)}& U\left(t\right)=\mathrm{1}-{\left\{\mathrm{1}+{\mathit{\beta }}_{\mathrm{2}}\exp \left[{\mathit{\beta }}_{\mathrm{1}}\left(t/t_{\text{p}}\right)-\mathrm{1}\right]\right\}}^{-\mathrm{1}/{\mathit{\beta }}_{\mathrm{2}}},\end{array}$$

where β₁ is the dimensionless rising coefficient determined by the watershed characteristics, β₂ is the dimensionless recessing coefficient affected by the downstream water surface condition, while t=t_p corresponds to the inflection point with the maximum instantaneous UH as follows:

$$\begin{array}{}\text{(8)}& {\left.{\displaystyle \frac{\mathrm{d}u}{\mathrm{d}t}}\right|}_{t=t_{\text{p}}}={\left.{\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}U}{\mathrm{d}{t}^{\mathrm{2}}}}\right|}_{t=t_{\text{p}}}=\mathrm{0}.\end{array}$$

With Eq. (7), the instantaneous UH u(t) is expressed as

$$\begin{array}{}\text{(9)}& \begin{array}{rl}& u\left(t\right)=\\ & {\displaystyle \frac{{\mathit{\beta }}_{\mathrm{1}}\exp \left[{\mathit{\beta }}_{\mathrm{1}}\left(t/t_{\text{p}}\right)-\mathrm{1}\right]{\left\{\mathrm{1}+{\mathit{\beta }}_{\mathrm{2}}\exp \left[{\mathit{\beta }}_{\mathrm{1}}\left(t/t_{\text{p}}\right)-\mathrm{1}\right]\right\}}^{-\left(\mathrm{1}+\mathrm{1}/{\mathit{\beta }}_{\mathrm{2}}\right)}}{t_{\text{p}}}}\end{array}.\end{array}$$

To illustrate the general and analytical solutions of the S-hydrograph from Eq. (7) and of the instantaneous UH from Eq. (9), Fig. 2 shows the computed U(t) (solid line) and u(t) (dashed line) for given values of ${\mathit{\beta }}_{\mathrm{1}}={\mathit{\beta }}_{\mathrm{2}}=$3 and t_p=10.

Figure 2 Illustration of the S-hydrograph and instantaneous UH analytically computed using Eqs. (7) and (9), respectively, for given values of ${\mathit{\beta }}_{\mathrm{1}}={\mathit{\beta }}_{\mathrm{2}}=\mathrm{3}$ and t_p=10.

2.2 Savenije's salt intrusion model

In estuaries, the key to deriving an empirical relationship for the salinity distribution is the dispersion coefficient, which is either constant (e.g., Gay and O'Donnell, 2007) or variable (e.g., Van der Burgh, 1972; Prandle, 1981). Based on the effective tidal average dispersion under steady-state conditions, Van der Burgh (1972) proposed the following empirical relationship for the dispersion coefficient:

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{\partial D}{\partial x}}=-K{\displaystyle \frac{\mathrm{|}Q\mathrm{|}}{A}},\end{array}$$

where D [L²T⁻¹] is the longitudinal dispersion coefficient, x [L] is the longitudinal coordinate measured in the upstream direction, Q [L³T⁻¹] is the fresh-water discharge, A [L²] is the tidally averaged cross-sectional area, and K is the dimensionless Van der Burgh coefficient.

It is assumed that the longitudinal cross-sectional area follows an exponential function as follows:

$$\begin{array}{}\text{(11)}& A=A_{\mathrm{0}}\exp \left(-x/a\right),\end{array}$$

where A₀ is the cross-sectional area at the estuary mouth, and a is the convergence length. Integration of Eq. (10) and taking into account the exponential variation of the cross-sectional area using Eq. (11) yields the following analytical description of the longitudinal effective dispersion (Savenije, 2005, 2012):

$$\begin{array}{}\text{(12)}& {\displaystyle \frac{D}{D_{\mathrm{0}}}}=\mathrm{1}-{\displaystyle \frac{Ka\mathrm{|}Q\mathrm{|}}{D_{\mathrm{0}}{A}_{\mathrm{0}}}}\left[\exp \left(x/a\right)-\mathrm{1}\right],\end{array}$$

where D₀ is the dispersion coefficient at the estuary mouth.

With Eq. (10), Savenije (2005, 2012) derived a one-dimensional empirical model for salt intrusion based on the tidally averaged cross-sectional mass conservation equation:

$$\begin{array}{}\text{(13)}& F=-\mathrm{|}Q\mathrm{|}S-AD{\displaystyle \frac{\partial S}{\partial x}},\end{array}$$

where F [MT⁻¹] and S [ML⁻³] are the tidally averaged salt flux and salinity, respectively.

In a steady-state situation with no net salt flux (i.e., F=0), Eq. (13) can be rewritten as follows:

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{\mathrm{d}S}{S}}=-{\displaystyle \frac{\mathrm{|}Q\mathrm{|}}{AD}}\mathrm{d}x.\end{array}$$

We can combine Eqs. (10) and (14) into a general relationship between the dispersion coefficient and salinity through the Van der Burgh's coefficient as follows (Savenije, 2005, 2012):

$$\begin{array}{}\text{(15)}& D/D_{\mathrm{0}}={\left(S/S_{\mathrm{0}}\right)}^K,\end{array}$$

where S₀ is the salinity concentration at the estuary mouth.

Combing Eqs. (12) and (15) yields the tidally averaged salinity along an estuary as follows (Savenije, 2005, 2012):

$$\begin{array}{}\text{(16)}& {\displaystyle \frac{S}{S_{\mathrm{0}}}}={\left(\mathrm{1}-{\displaystyle \frac{Ka\mathrm{|}Q\mathrm{|}}{D_{\mathrm{0}}{A}_{\mathrm{0}}}}\left[\exp \left(x/a\right)-\mathrm{1}\right]\right)}^{\mathrm{1}/K}.\end{array}$$

To make Eq. (16) dimensionless, we introduce the following dimensionless parameters:

$$\begin{array}{}\text{(17)}& S^{\ast }={\displaystyle \frac{S}{S_{\mathrm{0}}}}\mathit{\gamma }={\displaystyle \frac{c_{\mathrm{0}}}{\mathit{\omega }a}}{D}^{\ast }={\displaystyle \frac{\mathrm{|}Q\mathrm{|}{c}_{\mathrm{0}}}{D_{\mathrm{0}}{A}_{\mathrm{0}}\mathit{\omega }}}{x}^{\ast }={\displaystyle \frac{x\mathit{\omega }}{c_{\mathrm{0}}}},\end{array}$$

where S^∗ is the dimensionless salinity that is scaled by the value at the estuary mouth, γ represents the estuary shape number describing the convergence rate of an estuary, ω is the tidal frequency, D^∗ is the dimensionless dispersion coefficient representing the downstream dispersion condition, x^∗ is the dimensionless longitudinal coordinate that is normalized by the frictionless wavelength in prismatic channels, and c₀ is the classical wave celerity of a frictionless progressive wave, which is defined as

$$\begin{array}{}\text{(18)}& c_{\mathrm{0}}=\sqrt{gh/r_{\text{S}}},\end{array}$$

in which g [LT⁻²] is the acceleration due to gravity, h [L] is the tidally averaged depth, and r_S is the storage width ratio (see Savenije et al., 2008). Here, the asterisk denotes a dimensionless variable.

Thus, Eq. (16) can be rearranged as follows (Cai et al., 2015):

$$\begin{array}{}\text{(19)}& S^{\ast }={\left(\mathrm{1}-{\displaystyle \frac{D^{\ast }K}{\mathit{\gamma }}}\left[\exp \left(x^{\ast }\mathit{\gamma }\right)-\mathrm{1}\right]\right)}^{\mathrm{1}/K}.\end{array}$$

With Eq. (19), it is possible to derive an analytical expression for the salt intrusion length (i.e., the distance from the estuary mouth to the location where the water is totally fresh), which is obtained by setting $S^{\ast }=\mathrm{0}$ in Eq. (19) as follows:

$$\begin{array}{}\text{(20)}& L^{\ast }={\displaystyle \frac{\mathrm{1}}{\mathit{\gamma }}}\ln \left({\displaystyle \frac{\mathit{\gamma }}{D^{\ast }K}}+\mathrm{1}\right).\end{array}$$

3 General unit hydrograph theory for salt intrusion

Suppose there is an ocean coupling to an estuary with the coordinate origin located at the estuary mouth. If a unit volume of excess salinity from the ocean is locally (Δx→0) released into the estuary during the time required for an equilibrium to occur, the resulting hydrograph is the instantaneous UH $\mathrm{d}{S}^{\ast }\left(x\right)/\mathrm{d}x$ that corresponds to the S-hydrograph S^∗(x) in the dimensionless form. Similar to Guo's UH method (Guo, 2022a), we make four hypotheses on the instantaneous UH (representing the instantaneous rate of change with respect to the salinity at a specific position along the estuary axis, i.e., the salinity gradient) for the salinity distribution based on depth-average observations along estuaries (i.e., data facts).

• Hypothesis 1. Along the recessing limb, the salinity gradient $\mathrm{d}{S}^{\ast }\left(x\right)/\mathrm{d}x$ decreases exponentially, which makes the salinity S^∗(x) decay exponentially in a convex shape.

• Hypothesis 2. The salinity is scaled by the almost-constant salinity in the deep ocean, i.e., approximately 36 kg m⁻³; thus, as x tends towards negative infinity, the salinity gradient tends towards zero, and the salinity S^∗(x) tends towards 1.

• Hypothesis 3. Along the rising limb, the salinity gradient $\mathrm{d}{S}^{\ast }\left(x\right)/\mathrm{d}x$ increases exponentially, which makes the salinity S^∗(x) decay exponentially in a concave shape.

• Hypothesis 4. As x tends towards infinity, the salinity gradient tends towards zero, and the salinity S^∗(x) tends towards zero.

It should be noted that the above hypotheses 1 and 3 are, in principle, valid only for well-mixed or partially mixed estuaries, where salt intrusion really matters. From a practical perspective, this is not a restrictive assumption, since the salt wedge in highly stratified conditions only occurs at the time of high river discharge, when flood protection is generally the main concern and salt intrusion is not relevant (Savenije, 2005, 2012).

According to the first hypothesis, along the recessing limb, $\mathrm{d}{S}^{\ast }\left(x\right)/\mathrm{d}x$ and S^∗(x) satisfy the following relationship:

$$\begin{array}{}\text{(21)}& {\displaystyle \frac{\mathrm{d}{S}^{\ast }}{\mathrm{d}x}}=-\mathit{\mu }{S}^{\ast },\end{array}$$

where μ represents the recessing coefficient [L⁻¹]. Meanwhile, the second hypothesis requires that $\mathrm{d}{S}^{\ast }\left(x\right)/\mathrm{d}x=\mathrm{0}$ and that $S^{\ast }\left(x\right)=$1 (representing the constant salinity in the ocean) at $x\to -\mathrm{\infty }$. To meet this requirement, we revise Eq. (21) as follows:

$$\begin{array}{}\text{(22)}& {\displaystyle \frac{\mathrm{d}{S}^{\ast }}{\mathrm{d}x}}=-\mathit{\mu }\left(\mathrm{1}-S^{\ast },\right)\end{array}$$

since $S^{\ast }(-\mathrm{\infty })=\mathrm{1}$. Integrating Eq. (22) for S^∗(x) and applying the initial condition (i.e., x=0, $S^{\ast }=S_{\mathrm{0}}^{\ast })$ results in

$$\begin{array}{}\text{(23)}& S^{\ast }=\mathrm{1}+\left(S_{\mathrm{0}}^{\ast }-\mathrm{1}\right)\exp \left(\mathit{\mu }x\right).\end{array}$$

Similarly, according to the third hypothesis, along the rising limb, we have

$$\begin{array}{}\text{(24)}& {\displaystyle \frac{\mathrm{d}{S}^{\ast }}{\mathrm{d}x}}=-{\displaystyle \frac{\mathit{\mu }}{m}}{S}^{\ast },\end{array}$$

where m represents the dimensionless rising coefficient. Integrating Eq. (24) for S^∗(x) and applying the initial condition (i.e., x=0, $S^{\ast }=S_{\mathrm{0}}^{\ast })$ results in

$$\begin{array}{}\text{(25)}& S^{\ast }=S_{\mathrm{0}}^{\ast }\exp \left(-{\displaystyle \frac{\mathit{\mu }}{m}}x\right),\end{array}$$

which satisfies the fourth hypothesis. We can combine Eqs. (22) and (24) into a generalized differential equation as follows:

$$\begin{array}{}\text{(26)}& {\displaystyle \frac{\mathrm{d}{S}^{\ast }}{\mathrm{d}x}}=-{\displaystyle \frac{\mathit{\mu }}{m}}{S}^{\ast }\left(\mathrm{1}-S^{\ast m}\right),\end{array}$$

which reduces to Eq. (22) at $x\to -\mathrm{\infty }$ where $S^{\ast }=\mathrm{1}$ and Eq. (24) at x→∞ where $S^{\ast }\to \mathrm{0}$. The inflection point x=x_p, where

$$\begin{array}{}\text{(27)}& {\left.{\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}{S}^{\ast }}{\mathrm{d}{x}^{\mathrm{2}}}}\right|}_{x=x_{\text{p}}}=\mathrm{0}\end{array}$$

corresponds to the maximum absolute value of $\mathrm{d}{S}^{\ast }/\mathrm{d}x$ (i.e., the maximum salinity gradient). Integrating Eq. (26) for S^∗(x) and applying Eq. (27) results in

$$\begin{array}{}\text{(28)}& S^{\ast }={\left\{\mathrm{1}+m\exp \left[\mathit{\mu }\left(x-x_{\text{p}}\right)\right]\right\}}^{-\mathrm{1}/m}.\end{array}$$

To make μ dimensionless, we make a transform, namely μx_p→μ; then Eq. (28) can be revised as:

$$\begin{array}{}\text{(29)}& \begin{array}{rl}{S}^{\ast }=& {\left\{\mathrm{1}+m\exp \left[\mathit{\mu }\left(x/x_{\text{p}}-\mathrm{1}\right)\right]\right\}}^{-\mathrm{1}/m}\\ =& {\left\{\mathrm{1}+m\exp \left[\mathit{\mu }\left(x^{\ast }-\mathrm{1}\right)\right]\right\}}^{-\mathrm{1}/m}\end{array},\end{array}$$

where $x^{\ast }=x$/x_p is the dimensionless distance scaled by the position of the inflection point x_p. With Eq. (29), the instantaneous UH $\mathrm{d}{S}^{\ast }/\mathrm{d}x$ is written as

$$\begin{array}{}\text{(30)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{S}^{\ast }}{\mathrm{d}{x}^{\ast }}}=-\mathit{\mu }\exp \left[\mathit{\mu }\left(x^{\ast }-\mathrm{1}\right)\right]\\ & {\left\{\mathrm{1}+m\exp \left[\mathit{\mu }\left(x^{\ast }-\mathrm{1}\right)\right]\right\}}^{-\left(\mathrm{1}+\mathrm{1}/m\right)}\end{array},\end{array}$$

which satisfies the general UH definition; i.e., ${\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\frac{\mathrm{d}{S}^{\ast }}{\mathrm{d}{x}^{\ast }}=\mathrm{1}$.

It can be seen from Eq. (28) that the theoretical salt intrusion length L^∗ is not available, since S^∗ tends towards 0 as x approaches infinity asymptotically. However, it is possible to define a specific salt intrusion length for a given salinity threshold $S_f^{\ast }$ (such as 0.01) by substituting $S_f^{\ast }$ into Eq. (29) as follows:

$$\begin{array}{}\text{(31)}& L^{\ast }={\displaystyle \frac{\mathrm{1}}{\mathit{\mu }}}\ln \left({\displaystyle \frac{S_f^{\ast -m}-\mathrm{1}}{m}}\right)+\mathrm{1}.\end{array}$$

Figure 3 illustrates the spatial distribution of the S-hydrograph (salinity distribution S^∗) and its instantaneous UH (salinity gradient $\mathrm{d}{S}^{\ast }/\mathrm{d}{x}^{\ast })$ for given values of μ=1.5 and m=1.

Figure 3 S-hydrograph S^∗ (salinity concentration) and its instantaneous UH $\mathrm{d}{S}^{\ast }/\mathrm{d}{x}^{\ast }$ (salinity gradient) as a function of the dimensionless along-estuary distance x^∗ for given values of μ=1.5, m=1.

4 Results and discussion

4.1 Sensitivity analysis of the proposed salt intrusion model

Although Eqs. (29), (30) and (31) are analytical, the sensitivity to the two controlled parameters (μ and m) is not straightforward and directly clear. Thus, it is worthwhile to have a sensitivity analysis on the two calibrated parameters. Figure 4 presents the analytical solutions of the longitudinal salinity and its gradient as a function of μ and m. It can be clearly seen from Fig. 4 that two distinct estuarine regions display a very different behavior. For $x^{\ast }<\mathrm{1}$, we see an exponential decrease of the salinity gradient until a minimum value is reached at a critical position x=x_p (or $x^{\ast }=\mathrm{1}$) defined by Eq. (27), beyond which the salinity gradient increases exponentially until zero is reached asymptotically (Fig. 4c, d). The sensitivity analysis shows that the recessing coefficient μ determines the change rate of both the rising and recessing limbs (Fig. 4a, c). With regard to the rising coefficient m, it can be seen from Fig. 4b that the coefficient m exerts more influence along the rising limb and thus affects mainly the salinity distribution after the inflection point (Fig. 4b, d). In addition, Fig. 4c and d show that m=1 gives a symmetric salinity gradient of about x=x_p ($x^{\ast }=\mathrm{1}$, $S^{\ast }=\mathrm{0.5}$), $\le m<\mathrm{1}$ gives a negatively skewed salinity gradient, and m>1 gives a positively skewed salinity gradient.

Figure 4 Sensitivity analysis of the dimensionless salinity S^∗ and its gradient $\mathrm{d}{S}^{\ast }/\mathrm{d}x$ with regard to the recessing coefficient μ and the rising coefficient m.

To understand the response of the salt intrusion length to both calibrated parameters, Eq. (31) was used to analytically compute L^∗ for a wide range of μ and m values (Fig. 5) considering a salinity threshold of $S_f^{\ast }=\mathrm{0.01}$. We can clearly see that the isolines are almost linear, converging towards the origin of the m–μ diagram. Generally, the salt intrusion length increases with m, while it decreases with μ. This suggests that the model parameter m is generally proportional to the strength of tidal dynamics that induce dispersive transport of salt in the upstream direction, while the model parameter μ is proportional to the strength of the riverine flushing seaward. With this plot, it is possible to understand the potential impacts of different hydrodynamic conditions on the salt intrusion length, which is particularly useful to decision makers for salt intrusion prevention.

Table 1 Measured salinity distribution, calibrated parameters, computed salt intrusion length and model performance in terms of RMSE.

4.2 Application to real estuaries

The proposed salt intrusion model has been applied to observations in real estuaries worldwide with a wide range of different riverine and tidal hydrodynamics. In Table 1, a selection is presented for 21 estuaries where 89 salt intrusion measurements were collected at either high water slack (HWS) or low water slack (LWS; all observations are available on the web at https://salinityandtides.com/data-sources/, last access: 3 May 2023). In this table, the three model parameters x_p, μ and m were obtained by fitting Eq. (29) to the observed longitudinal salt intrusion. This can be easily done by means of a nonlinear curve-fitting method in the least-squares sense (such as using the MATLAB “lsqcurvefit.m” function). The model performance was evaluated by using the root-mean-square error (RMSE). Figure 6 shows the comparison of the observed and computed salinity concentrations in different estuaries. It can be seen from Fig. 6 that the correspondence with the observed salt intrusion is good, with the RMSE being 1.1 kg m⁻³ on average (see also the model performance reported in Table 1 for many estuaries worldwide with distinct salt intrusion lengths).

Figure 5 Response of the salt intrusion length L^∗ to the dimensionless parameters μ and m for a given salinity threshold $S_f^{\ast }=\mathrm{0.01}$.

For illustrations, Fig. 7 shows the longitudinal computation applied to the Pungue, Incomati and Limpopo estuaries by means of both the newly proposed and Savenije's models. Generally, the results of the two models are satisfactory for the following different shapes of salt intrusion curves in well-mixed or partially mixed estuaries (Savenije, 2005, 2012): (1) a dome-shape intrusion curve (such as along the Pungue estuary – Fig. 7a), which generally occurs in strong funnel-shaped estuaries; (2) a bell-shaped intrusion curve (such as along the Incomati estuary – Fig. 7c), which generally occurs in estuaries that have a trumpet shape; (3) a recession-shape intrusion curve (such as along the Limpopo estuary – Fig. 7e), which generally occurs in narrow estuaries with a near-prismatic shape and a high river discharge. It is worth noting that these types of salt intrusion curves are very much linked to the geometry of the estuary (Savenije, 2005, 2012). However, we observe from Table 1 that the calibrated model parameters (μ and m) are rather sensitive to the varied riverine and tidal forcing for a specific estuary. Thus, further studies concerning the relationship between the forcing conditions and the model parameters (μ and m) are required in the future. It can be seen from Fig. 7 that one important difference in terms of performance between these two models lies in the rising limb when the distance approaches infinity. As x tends towards infinity, the salinity gradient of the newly proposed model asymptotically approaches zero, while it reduces to zero at a critical position corresponding to the salt intrusion length from Eq. (20) for Savenije's model. This feature allows an improved fit with observations at the toe of the salt intrusion curve (e.g., Fig. 7a). Figures S1–S8 in the Supplement show the comparison between the observed longitudinal salinity and the analytically computed salt intrusion curves along 21 estuaries worldwide (see the Supplement).

Figure 6 Comparison between the analytically computed salinity concentrations and 89 observations in 21 estuaries worldwide.

Figure 7 Observed and analytically computed salt intrusion curves using the newly proposed and Savenije's models in the Pungue estuary (a, b), in the Incomati estuary (d, e) and in the Limpopo estuary (g, h), together with the idealized shape of the estuary ((c) funnel shape, (f) trumpet shape and (i)  prismatic shape).

Figure 8 Longitudinal variations of the dimensionless salinity S^∗ (a, b) and its gradient $\mathrm{d}{S}^{\ast }\mathrm{d}{x}^{\ast }$ (c, d) along the estuary axis for different values of input parameters. The solid and dashed lines indicate the solutions obtained by Savenije's model and the newly proposed model, respectively.

4.3 Analytical difference with Savenije's salt intrusion model

In order to understand the differences between the newly proposed model and empirical solutions based on the steady salt balance equation, we have made a comparison with the widely used Savenije's salt intrusion model (Savenije, 2005, 2012) by means of a Taylor expansion. The Taylor expansion of Savenije's salt intrusion model, i.e., Eq. (19), is written as

$$\begin{array}{}\text{(32)}& S^{\ast }=\mathrm{1}-D^{\ast }{x}^{\ast }+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(D^{\ast \mathrm{2}}-D^{\ast }\mathit{\gamma }-D^{\ast \mathrm{2}}K\right)x^{\ast \mathrm{2}}+O\left(x^{\ast \mathrm{3}}.\right)\end{array}$$

To make a comparison of the proposed model (i.e., Eq. 29) with Eq. (19), we introduce $-m\to K$ (Van der Burgh's coefficient); then Eq. (29) becomes

$$\begin{array}{}\text{(33)}& S^{\ast }={\left\{\mathrm{1}-K\exp \left[\mathit{\mu }\left(x^{\ast }-\mathrm{1}\right)\right]\right\}}^{\mathrm{1}/K}.\end{array}$$

The Taylor expansion of Eq. (33) is written as

$$\begin{array}{}\text{(34)}& \begin{array}{rl}& S^{\ast }={\left(\mathrm{1}-K{e}^{-\mathit{\mu }}\right)}^{\mathrm{1}/K}-{\left(\mathrm{1}-K{e}^{-\mathit{\mu }}\right)}^{(\mathrm{1}/K-\mathrm{1})}{e}^{-\mathit{\mu }}\mathit{\mu }{x}^{\ast }\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left(\mathrm{1}-K{e}^{-\mathit{\mu }}\right)}^{(\mathrm{1}/K-\mathrm{2})}{e}^{-\mathit{\mu }}{\mathit{\mu }}^{\mathrm{2}}\left(e^{-\mathit{\mu }}-\mathrm{1}\right)x^{\ast \mathrm{2}}+O\left(x^{\ast \mathrm{3}}\right)\end{array}.\end{array}$$

It is difficult to directly compare Eqs. (32) and (34). Alternatively, the difference between Eqs. (19) and (33) can be explored by looking at the exponential function parts, which can be expanded by the Taylor series as follows:

$$\begin{array}{}\text{(35)}& {\displaystyle \frac{D^{\ast }}{\mathit{\gamma }}}\left[\exp \left(x^{\ast }\mathit{\gamma }\right)-\mathrm{1}\right]=D^{\ast }x+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{D}^{\ast }\mathit{\gamma }{x}^{\ast \mathrm{2}}+O\left(x^{\ast \mathrm{3}}\right),\text{(36)}& {\displaystyle }\exp \left[\mathit{\mu }\left(x^{\ast }-\mathrm{1}\right)\right]=e^{-\mathit{\mu }}+e^{-\mathit{\mu }}\mathit{\mu }{x}^{\ast }+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{e}^{-\mathit{\mu }}{\mathit{\mu }}^{\mathrm{2}}{x}^{\ast \mathrm{2}}+O\left(x^{\ast \mathrm{3}}\right).\end{array}$$

Interestingly, if we slightly modified Eq. (35) by removing 1 from the brackets, then the Taylor expansion of Eq. (35) can be rewritten as

$$\begin{array}{}\text{(37)}& {\displaystyle \frac{D^{\ast }}{\mathit{\gamma }}}\left[\exp \left(x^{\ast }\mathit{\gamma }\right)\right]={\displaystyle \frac{D^{\ast }}{\mathit{\gamma }}}+D^{\ast }x+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{D}^{\ast }\mathit{\gamma }{x}^{\ast \mathrm{2}}+O\left(x^{\ast \mathrm{3}}\right).\end{array}$$

In this case, Eqs. (36) and (37) are identical if they satisfy the following conditions:

$$\begin{array}{}\text{(38)}& {\displaystyle }{D}^{\ast }=e^{-\mathit{\mu }}\mathit{\mu },\text{(39)}& {\displaystyle }\mathit{\gamma }=\mathit{\mu }.\end{array}$$

Thus, for given prior conditions Eqs. (38) and (39), the main difference between the newly proposed model and Savenije's model lies in the inclusion of the term – $K{D}^{\ast }/\mathit{\gamma }=m\exp (-\mathit{\mu })$ in the braces of Eq. (33), which is closely related to the upstream river discharge, the dispersion coefficient at the estuary mouth, the tidal frequency and the geometry of the estuary according to Eq. (17). As an illustration, Fig. 8 displays the longitudinal variations of the dimensionless salinity S^∗ and its gradient $\mathrm{d}{S}^{\ast }/\mathrm{d}{x}^{\ast }$ for different values of the recessing coefficient μ and the rising coefficient m, where the solid and dashed lines represent the solutions obtained by Savenije's model and the newly proposed model, respectively. It can be seen from Fig. 8 that the additional term mexp (−μ) mainly affects the downstream part of the salt intrusion curve, while the two methods tend to be the same for larger values of x^∗ (at the toe of the salt intrusion curve).

It is worth noting that the above analysis suggests that Van der Burgh's coefficient K (being equal to −m) should be negative rather than positive, since m is generally positive (see Table 1). This indicates that the dispersion coefficient D should be increased along the estuary axis according to Eq. (10). Figure 9 shows the analytically computed longitudinal dimensionless salinity S^∗ using Eq. (19) for a wide range of the input parameters K and D^∗ when γ=1. We can see that the computed S^∗ does converge to 0 when the distance x^∗ approaches infinity when K values are negative, which is very different from the performance of the previous analytical solutions using positive K values, where the computed S^∗ generally diverges for larger values of x^∗. Consequently, from a curve-fitting perspective, the Savenije's model using negative K values can be regarded as a special case of the newly proposed salt intrusion model if we further rescale the salinity by the dimensionless salinity in the deep ocean (see Fig. S9 in the Supplement). This also suggests that an enhanced empirical relationship concerning the effective dispersion coefficient (instead of the conventional Van der Burgh's relationship) is required for deriving an accurate salt intrusion curve from a theoretic point of view.

Figure 9 Longitudinal variation of the dimensionless salinity S^∗ computed using Eq. (19) along the estuary axis for different values of K and D^∗ when γ=1.

It should be noted that, although the model fits the observations very well, the physical foundation of Eq. (29) needs further study in the future. This limitation could be relaxed by carefully comparing the proposed model with those based on the steady-state salt balance equation. Specifically, several idealized numerical models (1-D, 2-D or 3-D models) have been adopted as a first approximation to quantify the along-channel salinity dynamics (e.g., Pein et al., 2018; Dijkstra and Schuttelaars, 2021; Wei et al., 2022). With these idealized numerical models, the physics behind Eq. (29) can be understood in order to make the model fully predictive through the relation of the three model parameters (x_p, μ and m) to measurable or quantifiable variables (e.g., river discharge, tidal amplitude and cross-sectional area convergence length) by means of regression techniques or similar approaches. It should be noted that the proposed salt intrusion model is, in principle, valid only for well-mixed or partially mixed estuaries, where salt intrusion really matters. From a practical perspective, this is not a restrictive assumption, since the salt wedge in highly stratified conditions only occurs at the time of high river discharge, when flood protection is generally the main concern and salt intrusion is not relevant (Savenije, 2005, 2012). Moreover, the proposed salt intrusion model is particularly useful for quantifying the alterations in salt intrusion dynamics owing to climate change or human interventions by comparing the three calibrated model parameters for different periods with considerably different conditions.

5 Conclusions

In this paper, we revisited the empirical salt intrusion model, making use of Guo's general unit hydrograph theory (Guo, 2022a, b, c), and proposed a general and analytical model for the salinity distribution in estuaries of the partially mixed to well-mixed types. The newly developed method does not require observed or calibrated salinity at the estuary mouth and can be well calibrated using a minimum of three salt measurements along the estuary axis. This is mainly due to the fact that Eq. (29) is a monotonic function whose first derivative (i.e., Eq. 30) does not change sign and has only one minimum. In addition, the salinity converges towards zero as the along-estuary distance approaches infinity asymptotically, which might improve the model performance near the toe of the salt intrusion curve when compared to empirical solutions based on the salt balance equation. The model has been applied to numerous estuaries worldwide, and the results agree very well with the observations. This indicates that the proposed model can be a useful tool to understand the dynamics of salt intrusion in estuaries and for assessing the potential impacts of human-induced (e.g., dredging) or natural (e.g., mean sea level rise) changes. However, the underlying physical foundation of the proposed model and the physics of the model parameters need further study in the future.