Analytical model captures intratidal variation in salinity in a convergent, well-mixed estuary

Knowledge of the processes governing salt intrusion in estuaries is important, since it influences the eco-environment of estuaries as well as its water resource potential in many ways. Analytical models of salinity variation offer a simple and efficient method for studying salt intrusion in estuaries. In this paper, an unsteady analytical solution is presented to predict the spatio-temporal variation in salinity in convergent estuaries. It is derived from a one-dimensional advection–diffusion equation for salinity, adopting a constant mixing coefficient and a single-frequency tidal wave, which can directly reflect the influence of the tidal motion and the interaction between the tide and runoff. The deduced analytical solution is illustrated with an application to the Humen estuary of the Pearl River Delta (PRD) and proves to be an efficient and accurate approach for predicting the salt intrusion in convergent estuaries. The unsteady analytical solution is tested against observations from six study sites to validate its capability to predict intratidal variation in salt intrusion. The results show that the proposed unsteady analytical solution can be successfully used to reproduce the spatial distribution and temporal processes governing salinity dynamics in convergent, well-mixed estuaries. The proposed method provides a quick and convenient approach for deciding on water-fetching methods to make good use of water resources.

1 Introduction

Salt intrusion in a river connecting to the sea is largely controlled by the river flow (Keulegan, 1966). The salinity of estuary waters is the result of the balance between river and tidal fluxes and mixing between them. The natural variability in river and tidal inputs to estuaries has been greatly disrupted as a result of the impact of global climate change and sea level rise as well as of local human activities, such as dam construction and channel dredging. These changes cause salt intrusion to become a serious problem in estuaries. It influences water quality, agricultural development in lowland areas, water utilization in upstream catchments, and the aquatic environment in estuaries (Han et al., 2010; Mo et al., 2007; Savenije, 1992). To address this issue worldwide, research efforts devoted to salt intrusion have been conducted in laboratory tanks, with numerical models, and using analytical approaches.

Nowadays, numerical models have become the most popular tool for studying salinity distribution in estuaries because they can provide visible results presenting the spatio-temporal variation in detail (e.g. Gong et al., 2012; Lerczak et al., 2006; MacCready, 2004; Wu and Zhu, 2010). However, the application of a numerical model is not an easy task, since it requires detailed data of the bathymetry and of hydrological boundary conditions, which are not available for all estuaries in the world. Here, a comparatively simple and convenient analytical model is developed as an efficient method for studying the salt intrusion in well-mixed estuaries. Analytical models are widely used because they are simple yet retain the basic physical characteristics involved. In the early 1960s, when systematic methods were developed to explore the factors controlling the instantaneous longitudinal salinity distribution, an expression was developed to compute the salt intrusion length as a function of the estuary length, mean depth, tidal amplitude, tidal period, fresh-water discharge, ocean salinity, and estuary roughness (Ippen and Harleman, 1961). In a subsequent period, analytical models of increasing complexity were developed based on the one-dimensional advection diffusion equation (Cameron and Pritchard, 1963), and on two-dimensional equations (Hansen and Rattray, 1965), capturing the dynamics of buoyancy-driven exchange flow and tidal mixing, satisfying salt conservation. Since the 1970s, numerous empirical and semi-empirical one-dimensional analytical models were put forward that correlated the salt intrusion length to the estuarine dynamical conditions and geomorphology based on the flume experiments and field measurements (e.g. Brockway et al., 2006; Fischer, 1974; Gay and O'Donnell, 2007, 2009; Kuijper and Van Rijn, 2011; Lewis and Uncles, 2003; Prandle, 1981, 1985; Rigter, 1973; Savenije, 1986). Although the literature on salt intrusion is vast, most studies concentrate on salt-water intrusion in a prismatic flume for reasons of convenience. However, the majority of estuaries in the world converge in width. The topography of the estuary is crucial to salt intrusion because the two main drivers (i.e. river flow and the tidal motion) both depend on the topography. The cross-sectional area determines the amount of the salt water entering the estuary and the efficiency of fresh water carrying salt out of the estuary. Savenije (1986) developed a fully analytical and predictive model to predict salt intrusion that applies to the natural topography of alluvial estuaries. It has been validated well in numerous estuaries where the width converges exponentially (e.g. Savenije, 1989; Savenije and Pagès, 1992; Nguyen and Savenije, 2006; Eaton, 2007; Ervine et al., 2007; Nguyen et al., 2008). In the years 2000–2010, another analytical approach (Brockway et al., 2006) was put forward which can be considered to be a modified and simplified version of the method presented in earlier studies (Prandle, 1981; Savenije, 1986). The dispersion coefficient in Brockway's model is assumed to be constant along the estuary, while it is assumed to be proportional to the spatial integral of the subtidal axial velocity in Savenije's model. In the theoretical models described above, the salt intrusion is predicted as a steady-state solution during slack water. Few studies have focussed on analysing the intratidal variation in salinity analytically. Song et al. (2008) proposed an unsteady-state model applicable to laboratory flumes and artificial channels where the cross section is assumed to be constant along the channel. Elaborating on the work of Song et al. (2008), here, an unsteady-state model is developed to predict the intratidal salinity intrusion dynamics in alluvial estuaries where the cross-sectional area typically converges. The aim of this study is to introduce a simple, unsteady analytical solution to the problem of predicting the intratidal variation in salt intrusion in convergent, well-mixed estuaries.

2 Methods

The cross-sectional area in this paper is described as an exponential function:

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

where A₀ is the cross-sectional area at the mouth (x=0), x is the distance along the estuary, and a is the convergence length of the cross-sectional area. The x axis has its origin at the mouth of the estuary, and the upstream direction is taken as positive. The one-dimensional advection–diffusion equation for salinity can be written as follows:

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial As}{\partial t}}+{\displaystyle \frac{\partial Aus}{\partial x}}={\displaystyle \frac{\partial }{\partial x}}\left(AD{\displaystyle \frac{\partial s}{\partial x}}\right),\end{array}$$

where s is salinity averaged over the cross section, t is time, u is the velocity, and D is the longitudinal dispersion coefficient. Although the assumption of a variable coefficient seems to be more reasonable, models with a constant dispersion coefficient have also proved to be capable of satisfactorily reproducing the salinity distribution (Lewis and Uncles, 2003; Brockway et al., 2006; Gay and O'Donnell, 2007, 2009). Under the assumption that D is independent of time, the salinity can be expanded in a Fourier series and be expressed as (Song et al., 2008)

$$\begin{array}{}\text{(3)}& s=\stackrel{\mathrm{‾}}{s}+s_{\mathrm{1}}\cos (\mathit{\omega }t+\mathit{\phi })+s_{\mathrm{2}}\sin (\mathit{\omega }t+\mathit{\phi }),\end{array}$$

where $\stackrel{\mathrm{‾}}{s}$ is the tide-averaged salinity, and s₁ and s₂ are coefficients. For the case of a simple harmonic wave with river discharge, the instantaneous flow velocity u is considered to consist of a time-dependent component $u_t=\mathit{\upsilon }\cos \left(\mathit{\omega }t+\mathit{\phi }\right)$ created by the tide and a steady component $u_{\mathrm{f}}=Q_{\mathrm{f}}/A$ contributed by the river flow:

$$\begin{array}{}\text{(4)}& u=u_{\mathrm{f}}+\mathit{\upsilon }\cos (\mathit{\omega }t+\mathit{\phi }),\end{array}$$

where the value of the runoff velocity u_f is negative. Introducing Eqs. (3) and (4) into Eq. (2), and using $\mathit{\theta }=\mathit{\omega }t+\mathit{\phi }$, yields

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& \left(u_{\mathrm{f}}{\displaystyle \frac{\partial s_{\mathrm{2}}}{\partial x}}-\mathit{\omega }{s}_{\mathrm{1}}\right)\sin \mathit{\theta }+\left(\mathit{\omega }{s}_{\mathrm{2}}+u_{\mathrm{f}}{\displaystyle \frac{\partial s_{\mathrm{1}}}{\partial x}}+\mathit{\upsilon }{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{s}}{\partial x}}\right)\cos \mathit{\theta }\\ & +u_{\mathrm{f}}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{s}}{\partial x}}+{\displaystyle \frac{\mathit{\upsilon }}{\mathrm{2}}}{\displaystyle \frac{\partial s_{\mathrm{1}}}{\partial x}}=\left(D{\displaystyle \frac{{\partial }^{\mathrm{2}}{s}_{\mathrm{2}}}{\partial x^{\mathrm{2}}}}-{\displaystyle \frac{D}{a}}{\displaystyle \frac{\partial s_{\mathrm{2}}}{\partial x}}\right)\sin \mathit{\theta }\\ & +\left(D{\displaystyle \frac{{\partial }^{\mathrm{2}}{s}_{\mathrm{1}}}{\partial x^{\mathrm{2}}}}-{\displaystyle \frac{D}{a}}{\displaystyle \frac{\partial s_{\mathrm{1}}}{\partial x}}\right)\cos \mathit{\theta }+D{\displaystyle \frac{{\partial }^{\mathrm{2}}\stackrel{\mathrm{‾}}{s}}{\partial x^{\mathrm{2}}}}-{\displaystyle \frac{D}{a}}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{s}}{\partial x}}.\end{array}\end{array}$$

As the equation should hold for all values of θ, Eq. (5) yields the following set of equations:

$$\begin{array}{}\text{(6)}& \left\{\begin{array}{l}-\mathit{\omega }{s}_{\mathrm{1}}+u_{\mathrm{f}}\frac{\partial s_{\mathrm{2}}}{\partial x}=D\frac{{\partial }^{\mathrm{2}}{s}_{\mathrm{2}}}{\partial x^{\mathrm{2}}}-\frac{D}{a}\frac{\partial s_{\mathrm{2}}}{\partial x},\\ \mathit{\omega }{s}_{\mathrm{2}}+u_{\mathrm{f}}\frac{\partial s_{\mathrm{1}}}{\partial x}+\mathit{\upsilon }\frac{\partial \stackrel{\mathrm{‾}}{s}}{\partial x}=D\frac{{\partial }^{\mathrm{2}}{s}_{\mathrm{1}}}{\partial x^{\mathrm{2}}}-\frac{D}{a}\frac{\partial s_{\mathrm{1}}}{\partial x},\\ u_{\mathrm{f}}\frac{\partial \stackrel{\mathrm{‾}}{s}}{\partial x}+\frac{\mathit{\upsilon }}{\mathrm{2}}\frac{\partial s_{\mathrm{1}}}{\partial x}=D\frac{{\partial }^{\mathrm{2}}\stackrel{\mathrm{‾}}{s}}{\partial x^{\mathrm{2}}}-\frac{D}{a}\frac{\partial \stackrel{\mathrm{‾}}{s}}{\partial x},\end{array}\right.\end{array}$$

where $\stackrel{\mathrm{‾}}{s}$, s₁, and s₂ can be further assumed to be

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& \stackrel{\mathrm{‾}}{s}=c_{\mathrm{0}}\exp \left(m\left(\exp \left({\displaystyle \frac{x}{a}}\right)-\mathrm{1}\right)\right),\\ & s_{\mathrm{1}}=c_{\mathrm{1}}\exp \left(m\left(\exp \left({\displaystyle \frac{x}{a}}\right)-\mathrm{1}\right)\right),\\ & s_{\mathrm{2}}=c_{\mathrm{2}}\exp \left(m\left(\exp \left({\displaystyle \frac{x}{a}}\right)-\mathrm{1}\right)\right),\end{array}\end{array}$$

with $c_{\mathrm{0}}={\stackrel{\mathrm{‾}}{s}}_{\mathrm{0}}$ being the tide-averaged salinity at the mouth of the estuary. Substitution of the Eq. (7) into the Eq. (6) yields

$$\begin{array}{}\text{(8)}& \left\{\begin{array}{l}-\mathit{\omega }{s}_{\mathrm{1}}+\left(u_{\mathrm{f}}+\frac{D}{a}\right)\frac{m}{a}\exp \left(\frac{x}{a}\right)s_{\mathrm{2}}\\ =D\frac{m}{a}\exp \left(\frac{x}{a}\right)\left(\frac{\mathrm{1}}{a}+\frac{m}{a}\exp \left(\frac{x}{a}\right)\right)s_{\mathrm{2}},\\ \mathit{\omega }{s}_{\mathrm{2}}+\mathit{\upsilon }\frac{m}{a}\exp \left(\frac{x}{a}\right)\stackrel{\mathrm{‾}}{s}+\left(u_{\mathrm{f}}+\frac{D}{a}\right)\frac{m}{a}\exp \left(\frac{x}{a}\right)s_{\mathrm{1}}\\ =D\frac{m}{a}\exp \left(\frac{x}{a}\right)\left(\frac{\mathrm{1}}{a}+\frac{m}{a}\exp \left(\frac{x}{a}\right)\right)s_{\mathrm{1}},\\ \left(u_{\mathrm{f}}+\frac{D}{a}\right)\frac{m}{a}\exp \left(\frac{x}{a}\right)\stackrel{\mathrm{‾}}{s}+\frac{\mathit{\upsilon }}{\mathrm{2}}\frac{m}{a}\exp \left(\frac{x}{a}\right)s_{\mathrm{1}}\\ =D\frac{m}{a}\exp \left(\frac{x}{a}\right)\left(\frac{\mathrm{1}}{a}+\frac{m}{a}\exp \left(\frac{x}{a}\right)\right)\stackrel{\mathrm{‾}}{s}.\end{array}\right.\end{array}$$

Then, further elaboration yields

$$\begin{array}{}\text{(9)}& \left\{\begin{array}{l}-\mathit{\omega }{c}_{\mathrm{1}}+\left(u_{\mathrm{f}}-\frac{Dm}{a}\exp \left(\frac{x}{a}\right)\right)\frac{m}{a}\exp \left(\frac{x}{a}\right)c_{\mathrm{2}}=\mathrm{0},\\ \mathit{\omega }{c}_{\mathrm{2}}+\mathit{\upsilon }\frac{m}{a}\exp \left(\frac{x}{a}\right)c_{\mathrm{0}}\\ +\left(u_{\mathrm{f}}-\frac{Dm}{a}\exp \left(\frac{x}{a}\right)\right)\frac{m}{a}\exp \left(\frac{x}{a}\right)c_{\mathrm{1}}=\mathrm{0},\\ \left(u_{\mathrm{f}}-\frac{Dm}{a}\exp \left(\frac{x}{a}\right)\right)\frac{m}{a}\exp \left(\frac{x}{a}\right)c_{\mathrm{0}}\\ +\frac{\mathit{\upsilon }}{\mathrm{2}}\frac{m}{a}\exp \left(\frac{x}{a}\right)c_{\mathrm{1}}=\mathrm{0}.\end{array}\right.\end{array}$$

Hence, the following solutions can be obtained:

$$\begin{array}{}\text{(10)}& \left\{\begin{array}{l}{c}_{\mathrm{1}}=\frac{\left(u_{\mathrm{f}}-\frac{Dm}{a}\exp \left(\frac{x}{a}\right)\right)\frac{m}{a}\exp \left(\frac{x}{a}\right)}{\mathit{\omega }}{c}_{\mathrm{2}},\\ c_{\mathrm{2}}=\frac{-\mathit{\upsilon }\frac{m}{a}\exp \left(\frac{x}{a}\right)\mathit{\omega }}{{\mathit{\omega }}^{\mathrm{2}}+{\left(u_{\mathrm{f}}-\frac{Dm}{a}\exp \left(\frac{x}{a}\right)\right)}^{\mathrm{2}}\frac{m^{\mathrm{2}}}{a^{\mathrm{2}}}\exp \left(\frac{\mathrm{2}x}{a}\right)}{c}_{\mathrm{0}},\\ m=\frac{u_{\mathrm{f}}a\exp \left(-\frac{x}{a}\right)}{D}=\frac{Q_{\mathrm{f}}a}{D{A}_{\mathrm{0}}}.\end{array}\right.\end{array}$$

The analytic solution of the unsteady-state salinity distribution is therefore represented as

$$\begin{array}{}\text{(11)}& \begin{array}{rl}& s={\stackrel{\mathrm{‾}}{s}}_{\mathrm{0}}\exp \left({\displaystyle \frac{Q_{\mathrm{f}}a}{D{A}_{\mathrm{0}}}}\left(\exp \left({\displaystyle \frac{x}{a}}\right)-\mathrm{1}\right)\right)\\ & \cdot \left(\mathrm{1}-{\displaystyle \frac{\mathit{\upsilon }{u}_{\mathrm{f}}}{\mathit{\omega }D}}\sin \left(\mathit{\omega }t+\mathit{\phi }\right)\right).\end{array}\end{array}$$

By integrating this unsteady salinity expression over the tidal period T, the salt intrusion under tidal average (TA) conditions, as defined by Brockway et al. (2006), can be obtained by

$$\begin{array}{}\text{(12)}& {\stackrel{\mathrm{‾}}{s}}_x={\stackrel{\mathrm{‾}}{s}}_{\mathrm{0}}\exp \left({\displaystyle \frac{Q_{\mathrm{f}}a}{D{A}_{\mathrm{0}}}}\left(\exp \left({\displaystyle \frac{x}{a}}\right)-\mathrm{1}\right)\right).\end{array}$$

A graph of the logarithm of salinity $\ln \left(\stackrel{\mathrm{‾}}{s}/{\stackrel{\mathrm{‾}}{s}}_{\mathrm{0}}\right)$ against exp (x∕a) should be a straight line, with the slope inversely proportional to the longitudinal dispersion coefficient D (Brockway et al., 2006). The coefficient D can then be calculated from

$$\begin{array}{}\text{(13)}& D={\displaystyle \frac{Q_{\mathrm{f}}a}{k{A}_{\mathrm{0}}}},\end{array}$$

where k is the slope of the fitted line. This approach makes it possible to estimate the longitudinal dispersion coefficient D based on the measurements of salinity made during a survey.

The tidal velocity amplitude υ can be estimated as $\mathit{\upsilon }=E\mathit{\pi }/T$ (Savenije, 1993), where E is the tidal excursion and the harmonic constant ω is given as $\mathit{\omega }=\mathrm{2}\mathit{\pi }/T$. Introducing $\mathit{\upsilon }=E\mathit{\pi }/T$ and $\mathit{\omega }=\mathrm{2}\mathit{\pi }/T$ into Eq. (11), and using $u_{\mathrm{f}}=Q_{\mathrm{f}}/A$, yields

$$\begin{array}{}\text{(14)}& \begin{array}{rl}& s={\stackrel{\mathrm{‾}}{s}}_{\mathrm{0}}\exp \left({\displaystyle \frac{Q_{\mathrm{f}}a}{D{A}_{\mathrm{0}}}}\left(\exp \left({\displaystyle \frac{x}{a}}\right)-\mathrm{1}\right)\right)\\ & \cdot \left(\mathrm{1}+{\displaystyle \frac{E}{\mathrm{2}a}}\left(-{\displaystyle \frac{Q_{\mathrm{f}}a}{D{A}_{\mathrm{0}}}}\exp \left({\displaystyle \frac{x}{a}}\right)\right)\sin (\mathit{\omega }t+\mathit{\phi })\right).\end{array}\end{array}$$

This expression can be used to describe the temporal and spatial variation in salinity, including high water slack (HWS) and low water slack (LWS), when the tidal discharge is zero by definition. Since the maximum salinity is reached at HWS and the minimum salinity is reached at LWS (Savenije, 2005), Eq. (14) can be simplified for HWS into

$$\begin{array}{}\text{(15)}& \begin{array}{rl}& s_{max}={\stackrel{\mathrm{‾}}{s}}_{\mathrm{0}}\exp \left({\displaystyle \frac{Q_{\mathrm{f}}a}{D{A}_{\mathrm{0}}}}\left(\exp \left({\displaystyle \frac{x}{a}}\right)-\mathrm{1}\right)\right)\\ & \cdot \left(\mathrm{1}+{\displaystyle \frac{E}{\mathrm{2}a}}\left(-{\displaystyle \frac{Q_{\mathrm{f}}a}{D{A}_{\mathrm{0}}}}\exp \left({\displaystyle \frac{x}{a}}\right)\right)\right),\end{array}\end{array}$$

and can be simplified for LWS into

$$\begin{array}{}\text{(16)}& \begin{array}{rl}& s_{min}={\stackrel{\mathrm{‾}}{s}}_{\mathrm{0}}\exp \left({\displaystyle \frac{Q_{\mathrm{f}}a}{D{A}_{\mathrm{0}}}}\left(\exp \left({\displaystyle \frac{x}{a}}\right)-\mathrm{1}\right)\right)\\ & \cdot \left(\mathrm{1}-{\displaystyle \frac{E}{\mathrm{2}a}}\left(-{\displaystyle \frac{Q_{\mathrm{f}}a}{D{A}_{\mathrm{0}}}}\exp \left({\displaystyle \frac{x}{a}}\right)\right)\right).\end{array}\end{array}$$

The tidal excursion E, the distance over which a water particle travels up and down the estuary with the flooding and ebbing tide, is assumed to decrease exponentially along the channel:

$$\begin{array}{}\text{(17)}& E=E_{\mathrm{0}}\exp \left(-x/e\right),\end{array}$$

where E₀ is the tidal excursion at the mouth (x=0), and e is the damping length of the tidal excursion. Thus, the combination of Eqs. (15), (16), and (17) yields

$$\begin{array}{}\text{(18)}& E={\displaystyle \frac{a\left(s_{max\mathrm{0}}-s_{min\mathrm{0}}\right)}{{\stackrel{\mathrm{‾}}{s}}_{\mathrm{0}}\left(-\frac{Q_{\mathrm{f}}a}{D{A}_{\mathrm{0}}}\right)}}\exp \left(-x/e\right),\end{array}$$

where s_max0 is the maximum salinity at the estuary mouth and s_min0 is the minimum salinity at the estuary mouth.

Since the tidal flow is assumed to vary as a simple harmonic wave, the unsteady salinity model is here presented in its simplest form, with a single frequency. As the tidal propagation celerity in the estuary is assumed to be constant, the tidal phase at each site can be made up of an initial phase φ₀ at the mouth of the estuary and a phase difference that is the travel time of the tide from the mouth to the study site. Therefore, Eq. (14) can be modified as

$$\begin{array}{}\text{(19)}& \begin{array}{rl}& s={\stackrel{\mathrm{‾}}{s}}_{\mathrm{0}}\exp \left({\displaystyle \frac{Q_{\mathrm{f}}a}{D{A}_{\mathrm{0}}}}\left(\exp \left({\displaystyle \frac{x}{a}}\right)-\mathrm{1}\right)\right)\\ & \cdot (\mathrm{1}+\left(-{\displaystyle \frac{E_{\mathrm{0}}}{\mathrm{2}a}}{\displaystyle \frac{Q_{\mathrm{f}}a}{D{A}_{\mathrm{0}}}}\right)\exp \left({\displaystyle \frac{x}{a}}-{\displaystyle \frac{x}{e}}\right)\\ & \cdot \sin \left(\mathit{\omega }\left(t-{\displaystyle \frac{x}{c}}\right)+{\mathit{\phi }}_{\mathrm{0}}\right)),\end{array}\end{array}$$

where c is the tidal propagation celerity.

We note that in the approach presented above, the tidal excursion at the mouth is inferred from salinity data, whereas an alternative theoretical approach may be applicable that is less dependent on in situ data. Tidal-wave propagation can be described analytically by a set of four implicit equations (Cai et al., 2012), namely the phase lag equation $\tan \left(\mathit{\epsilon }\right)=\mathit{\lambda }/\left(\mathit{\gamma }-\mathit{\delta }\right)$, the scaling equation $\mathit{\mu }=\sin \left(\mathit{\epsilon }\right)/\mathit{\lambda }$, the damping equation $\mathit{\delta }=\mathit{\gamma }/\mathrm{2}-\mathrm{4}\mathit{\chi }\mathit{\mu }/\left(\mathrm{9}\mathit{\pi }\mathit{\lambda }\right)-\mathit{\chi }{\mathit{\mu }}^{\mathrm{2}}/\mathrm{3}$, and the celerity equation ${\mathit{\lambda }}^{\mathrm{2}}=\mathrm{1}-\mathit{\delta }\left(\mathit{\lambda }-\mathit{\delta }\right)$, where λ is the celerity number $\mathit{\lambda }=c_{\mathrm{0}}/c$, μ is the velocity number $\mathit{\mu }=\mathit{\upsilon }\stackrel{\mathrm{‾}}{h}/\left(r_{\mathrm{s}}\mathit{\eta }{c}_{\mathrm{0}}\right)$, δ is the damping number $\mathit{\delta }=c_{\mathrm{0}}\mathrm{d}\mathit{\eta }/\left(\mathit{\eta }\mathrm{d}x\mathit{\omega }\right)$, and ε is the phase lag between high water (HW) and HWS $\mathit{\epsilon }=\mathit{\pi }/\mathrm{2}-\left({\mathit{\varphi }}_Z-{\mathit{\varphi }}_U\right)$. Here three dimensionless parameters control the tidal hydrodynamics (Savenije et al., 2008), i.e. the dimensionless tidal amplitude $\mathit{\zeta }=\mathit{\eta }/\stackrel{\mathrm{‾}}{h}$, the estuary shape number $\mathit{\gamma }=c_{\mathrm{0}}/\left(\mathit{\omega }a\right)$, and the friction number $\mathit{\chi }=r_{\mathrm{s}}g{c}_{\mathrm{0}}\mathit{\zeta }/\left(K_{\mathrm{s}}^{\mathrm{2}}\mathit{\omega }{\stackrel{\mathrm{‾}}{h}}^{\mathrm{4}/\mathrm{3}}\right){\left[\mathrm{1}-{\left(\mathrm{4}\mathit{\zeta }/\mathrm{3}\right)}^{\mathrm{2}}\right]}^{-\mathrm{1}}$, where η is the tidal amplitude, K_s is the Manning–Strickler friction coefficient, r_s the storage width ratio, $\stackrel{\mathrm{‾}}{h}$ is the tide-averaged depth, and c₀ is the classical wave celerity $c_{\mathrm{0}}=\sqrt{g\stackrel{\mathrm{‾}}{h}/r_{\mathrm{s}}}$. Then with available geometry and friction data at the estuary mouth, the tidal propagation celerity and the tidal amplitude (or the tidal excursion) can be obtained by solving the set of four equations. Rather than proceeding with this analytical model for tidal hydrodynamics, hereafter we employ Eq. (18) to close the set of equations.

3 Study area and data

3.1 Overview of study area

The Pearl River estuary (PRE) is located midway along the northern boundary of the South China Sea. It receives a large amount of fresh water from the Pearl River, which has three major branches (i.e. the West River, the North River, and the East River) in the upper drainage basin. The annual river discharge, with 80 % occurring in the wet season, empties into the South China Sea via eight outlets (Zhao, 1990). The Lingding bay is created by the inflows of fresh water from the Pearl River through four major discharge outlets, namely Humen, Jiaomen, Hongqimen, and Hengmen. Historically, about 50 %–55 % of the river flow enters the Lingding bay, while the remaining fresh water directly flows into the South China Sea through the four southwestern outlets (i.e. Modaomen, Jitimen, Hutiaomen, and Yamen).

The Humen is the largest river outlet in the Lingding bay and contributes 34.6 % of the water discharge, i.e. about 603×10⁸ m³ in terms of annual water discharge (Ren et al., 2006). The fresh-water input into Lingding bay through the Humen outlet comes from three sources: the East River, the Liuxihe River, and the North River. The annual river discharge with a peak of 1870 m³ s⁻¹, measured at the Niuxinling station in Liuxihe River, is about 10 times less than that with a flood peak in excess of 12 000 m³ s⁻¹, measured in the other two rivers (Luo et al., 2002).

The tide in the Pearl River estuary has a mixed semidiurnal–diurnal character (Zhang et al., 2012). Among the eight outlets of the Pearl River estuary, Humen is most strongly dominated by the tide, with an annual average and maximum tidal range of 1.63 and 2.59 m, respectively, at the mouth of the estuary (Li and Lei, 1998).

As a major tributary of the Pearl River, the Humen estuary can be divided into two waterways: the Guangzhou channel (the upper reach), with an average width of 431 m, and the Shiziyang channel (the lower reach), which is about 2200 m wide (Mai et al., 2001). It is a NW–SE branch of the Pearl River estuary, with a width of about 4 km at the mouth, resembling an inverted funnel with a narrow neck in the north and a wide mouth opening to the south. The Humen outlet has the highest tidal prism in the Pearl River estuary due to the large width of the mouth, resulting in strong tidal motion. Especially during spring tide in the dry season, when the river discharge is lowest, the downstream area becomes saline.

3.2 Data

The information available for the model application in this study includes data on topography, salinity, river discharge, and the tidal flow. A field survey for salt intrusion was conducted during the dry season in 2005. It was a project carried out by Guangdong Province Hydrology Bureau and the Pearl Hydrology Bureau from the River Conservancy Commission. In this paper, the field data from 29 January to 3 February were used, which were measured at six gauge stations along the channel (Fig. 1). Considering the impact of shipping, the measuring positions were near the banks, with certain distances ranging from 605 to 70 m. A Global Positioning System was used to confirm the exact measuring locations (Table 1). The Humen estuary is well-mixed under normal flow conditions during the dry season (Ou, 2009; Luo et al., 2010). Due to 3 years of drought, the river discharge decreased by 50 % during the study period in 2005 compared to a normal year (Liao et al., 2008). Thus, there is no doubt that well-mixed conditions prevailed during the calibration and validation. The average salinity of vertical profiles was calculated based on the hourly water samples. At each location, the saline water was sampled at two different elevations: at one-fifth and four-fifths of the depth of channel from the bed, and salinity was obtained using a salimeter. The water discharge at stations was provided by the hydrology bureaus during the field survey. The cross section was measured at mean sea level, with the help of an ultrasonic echo sounder.

Table 1 General information of hydrological stations in the Humen waterway.

^∗ The coordinate system's origin is set at 22^∘05^′12.9894^′′ N, 113^∘27^′34.9899^′′ E.

Figure 1 Map of the Humen estuary, showing the gauging stations where salinity concentration was measured during the field survey from 29 January to 3 February 2005.

Because of the complex river network upstream of the Humen area in the Pearl River estuary, the river discharge is difficult to determine. The total flux through the Humen outlet is composed of three parts which come from three main sources: the East River, the North River, and the Liuxihe River. The river discharge used in this paper was measured at upstream stations (Sanshui for the North River, Boluo for the East River, and Laoyagang for the Liuxihe River) from 29 January to 3 February. These data were collected from the official databases of the Hydrology Bureaus mentioned above. In the lower reach of the East River and the Liuxihe River, respectively, the Boluo station and Laoyagang station are located about 80 km upstream from the Humen outlet. The daily discharge measured at the Boluo station varied from 260 to 400 m³ s⁻¹ during the survey period, while it was about 20 m³ s⁻¹ at the Laoyagang station. The discharge of the East River and the Liuxihe River entirely flows toward the South China Sea through the Humen estuary (Ren et al., 2011). The North River is an important source for the river discharge to the Humen outlet. River discharge from the North River reaches the Humen channel through a network of channels which connects to the western part of the Pearl River delta. The Sanshui station is the primary hydrological station in the North River. About 11 % of the measured discharge flowed into the Humen estuary during the survey in 2005. The Sanshui station is located further upstream than the other two stations (Boluo and Laoyagang). The response lag of salinity variation at the estuary mouth to discharge variation at the Sanshui station is about 2 d, while the river flow spends about 1 d to travel from the Boluo and Laoyagang station to the estuary mouth.

4 Results

4.1 Model calibration

To demonstrate the practical application of the proposed analytical solution, the model has been used to simulate and analyse the spatio-temporal variation in salt intrusion in the Humen estuary. In the following, the parameters of the analytical solution are obtained from calibration.

The spatial decay of the cross-sectional area of the Humen estuary can be described by the exponential function expressed in Eq. (1). The field data (triangles) and the best-fit line are shown in Fig. 2. The cross-sectional area at the mouth at mean tide, A₀, is calculated as 37 822 m², and the convergence length of cross section a is obtained by curve fitting as 16.7 km.

Figure 2 Shape of the Humen estuary, showing the correlation between the cross-sectional area A (m²) and the distance from the estuary mouth x (km). The coefficient of determination R² is 0.92. The triangles represent observations, and the line represents the fit to Eq. (1), where the area at the estuary mouth is A₀ = 37 822 m² and the area convergence length (a) is 16.7 km.

The relative salinity is plotted as ln (s∕s₀) against exp (x∕a) in Fig. 3. There is a straight fit between these two variables, confirming the constancy of the dispersion coefficient. The dispersion coefficient D can be computed from the slope of the fitted lines according to Eq. (12), where k is the slope. This approach has been shown previously to be efficient (e.g. Brockway et al., 2006; Fang et al., 2006; Zhang et al., 2010). Table 2 shows the results of the fit for all these surveys carried out between 29 January to 3 February (the slope in column 4 and the dispersion coefficient in column 6). The dispersion coefficient estimates obtained from this fitting procedure can be interpreted as a spatial average, representing the entire reach. The coefficient of determination (R²) lies in the range between 0.85 and 0.92, with a mean value of 0.89. The assumption of the dispersion coefficient independent of distance is demonstrated to be reasonable and acceptable in the present case. The dispersion coefficient from data on 29 January is therefore used as the calibrated value.

Figure 3 Relative salinity concentration along the Humen estuary. The circlers represent observations, and the lines represent the fit to Eq. (12). s is the salinity at distance x from the estuary mouth, s₀ is the salinity at the mouth, and a is the convergence length of the cross-sectional area.

Table 2 Dispersion coefficient of salt intrusion in Humen estuary.

The tidal excursion at the mouth of the estuary is obtained through Eq. (18). For each tidal excursion at each day, the period-averaged value, maximum value, and minimum value of salinity at the mouth are obtained by statistical analysis, and the longitudinal dispersion coefficient D is computed by linear fitting, as shown previously. Moreover, the damping length of the tidal excursion e is calibrated using the observed salinity along the estuary. Similar to the tidal excursion, the value of the propagation celerity c is also obtained by calibration based on observations. The initial tidal phase φ₀ is calculated via a reverse procedure by calibration on the salinity at the mouth of estuary. Data from 31 January are used to verify the change of salinity over a tidal cycle.

The six calibration parameters (i.e. the convergence length of cross section a, the dispersion coefficient D, the tidal excursion E₀, the damping length of the tidal excursion e, the initial phase φ₀, and the tidal celerity c) are obtained based on the measurements at the mouth of the estuary, as shown in Table 3. Based on the observed data on 29 January, the results of the model calibration can be seen in Fig. 4.

Table 3 Calibrated values of parameters.

Figure 4 Comparison between calibration results and measured salinity concentration along the river on 29 January 2005, showing values of measured salinity at high water slack (circles) and low water slack (inverted triangles) and the calibrated salinity curves at high water slack (red curve) and low water slack (blue curve).

4.2 Model validation

A validation of the unsteady model is offered in two separate parts, i.e. the longitudinal distribution of salinity along the channel and the temporal variation in salinity during the tidal period. In the first part, observations during two characteristic conditions (i.e. HWS and LWS) are chosen for comparison against the calculated results of the salinity distribution. In the second part, a model for expressing the change process of salinity during tidal periods is established, according to the measurement on 31 January.

4.2.1 Longitudinal distribution of salinity

Based on the field measurements from 30 January to 3 February, Eqs. (15), (16), and (18) are used to calculate the longitudinal variation in salinity. Conditions of neap tide are considered to last from 31 January to 2 February. The calibration results are presented in Fig. 5. The good agreement between the computation and the measured data indicates that the performance of the unsteady analytical model is to a certain extent satisfactory in the Humen estuary. The analytical model is found to better reproduce the distribution of salinity at high water (HW) than at low water (LW). This can be attributed to different degrees of mixing, which is stronger at HW. As the estuary is assumed to be well-mixed, the analytical model undoubtedly will perform better when mixing is higher. Fluctuations around the theoretical curve may partly be caused by the unequal distribution of salinity over the cross section or by the indirect derivation of the salinity at HWS and LWS, which is replaced with the daily maximum and minimum values, respectively.

Figure 5 Comparison between validation result and measured salinity concentration along the river from 30 January to 3 February 2005.

It can be seen that the analytical model substantially overestimates the salinity in the downstream part of the estuary partly because of the special locations of the stations (some are located at the confluence of rivers). The expression for the distribution analysis of salinity, Eq. (11), is multiplied by the tidal average salinity with an extra component that reflects the effect of the tide and the interaction of the tide and river flow. This time-dependent component is a sine function, namely $\left(-\mathit{\upsilon }{u}_{\mathrm{f}}/\mathit{\omega }D\right)\sin (\mathit{\omega }t+\mathit{\phi })$; thus the calculated salinity at HWS and LWS is always symmetrical about the average values. The symmetry property of salinity has been demonstrated by Savenije (1989) under the assumption that the tidal excursion is independent of distance. After a transformation of variables, the sine function mentioned above is expressed as $\left(E/\mathrm{2}a\right)(-Q_{\mathrm{f}}a/DA)\sin (\mathit{\omega }t+\mathit{\phi })$, where the dispersion coefficient plays an important role. To simplify and clarify the interaction between the tidal motion and the river flow, the parameters D and Q_f are combined into one single calibration variable, the mixing coefficient α:

$$\begin{array}{}\text{(20)}& \mathit{\alpha }=-D/Q_{\mathrm{f}},\end{array}$$

which can be obtained in the same way as the dispersion coefficient. In this paper, the mixing coefficient is assumed to be constant along the channel to develop a comparatively simple analytical solution within acceptable levels. It is calibrated by the measurements from the Dahu station to the Huangpuyou station, located at the junction of two reaches in the Humen estuary.

4.2.2 Periodic variation in salinity

The observations of salinity at hourly intervals along the Humen estuary are used to calibrate the dispersion coefficient in the model and to analyse the change of salinity with time. The results indicate that the calibrated unsteady analytical model fits the observations well. Figure 6, where the analytical solution is compared with observation, demonstrates that the proposed unsteady analytical solution is able to reflect the change process of salinity over a tidal cycle. Additionally, the simplification and assumption of the tidal celerity (c) and the initial phase at the mouth of estuary (φ₀) in Eq. (19) prove to be realistic.

Figure 6 Comparison between the predicted and measured salinity concentration over time on 31 January (neap tide) at each study site, showing that the analytical model captures the temporal variation in salinity. The hourly salinity measurements are represented by rectangles, while the simulated salinity varying with time is represented by the red solid line. In the figure, x is the distance from the estuary mouth.

The theoretical results of the periodic variation in salinity are not always consistent with the observations. As can be seen in Fig. 6, the analytical model for simulating the temporal process of salinity has a relatively poor performance at the sites near the mouth of estuary, such as the Dahu station. By comparing the variation in salinity at different sites (Fig. 6), it shows that salinity variation is more symmetrical further away from the study site. The discrepancies near the mouth may have three reasons. Firstly, lateral residual circulation usually exists at the mouth of an estuary, where the cross section is widest. Secondly, the mouth of estuary is close to Lingding bay, where the salt dynamics are influenced by coastal and ocean currents. Thirdly, near the outlet, comprehensive salinity measurements are much more difficult to take due to the impact of tidal flats and complex hydrodynamics, influenced by the Coriolis force and wind effects. All the influences above are related to the width of the channel, which gradually decreases in the landward direction.

The observations at the Machong station show some non-periodic variation, which may relate to the proximity of the confluence of the East River and the Shiziyang channel. At the Dasheng station, about 2.6 km upstream from the Machong station and near another confluence, the simulated temporal process of salinity shows fairly good agreement with the observations. To understand the irregular changes of salinity at the Machong station, the daily averaged discharges at the Machong and Dasheng stations are analysed by integrating over the tidal period. The results are presented in Fig. 7, where the positive values represent the mean discharge transporting in the seaward direction. At the Machong station, the mean discharge is directed inland, which can be attributed to Stokes transport (Buschman et al., 2010; Hoitink and Jay, 2016). At the Dasheng station, only a few kilometres upstream, the mean discharge is seaward, as expected. The tide-averaged discharge thus converges in the estuary during low river flow, which will increase the total water volume in the estuary, and create a mean water level rise. We expect that this process has an impact on the mean salt balance, which explains part of the observed discrepancies.

Figure 7 Subtidal discharge measured at Machong station and Dasheng station from 29 January to 3 February. Positive values indicate discharge in the seaward direction.

For comparison, the result obtained by Song's model is also presented here. The unsteady analytical model developed by Song et al. (2008) can reproduce the salinity process in an idealized estuary with constant depth and constant width, which is expressed as

$$\begin{array}{}\text{(21)}& s={\stackrel{\mathrm{‾}}{s}}_{\mathrm{0}}\exp \left({\displaystyle \frac{u_{\mathrm{f}}}{D_S}}x\right)\left(\mathrm{1}-{\displaystyle \frac{u_{\mathrm{f}}\mathit{\upsilon }}{D_S\mathit{\omega }}}\sin \left(\mathit{\omega }t+\mathit{\phi }\right)\right).\end{array}$$

Therefore, in fact, it is more suitable for use in prismatic channels. The dispersion coefficient of Song's model is assumed to be independent of the distance and can be estimated by

$$\begin{array}{}\text{(22)}& D_S=-{\displaystyle \frac{\stackrel{\mathrm{‾}}{s}{u}_{\mathrm{f}}\mathit{\upsilon }}{\mathrm{0.5}\mathit{\omega }\left(s_{max}-s_{min}\right)}}.\end{array}$$

When an estimation for tidal velocity is made according to the relation $\mathit{\upsilon }=E\mathit{\pi }/T$, and the value of the runoff velocity is obtained using $u_{\mathrm{f}}=Q_{\mathrm{f}}/A$, then the dispersion coefficient can be calculated based on the measured salinity at the mouth of the estuary. The data which have been used for modelling in the Humen estuary can be found in Table 4. As shown in Fig. 8, the performance of Song's model for the Humen estuary is satisfactory at the study sites close to the estuary mouth, e.g. the Dahu station and the Sishengwei station. However, the salt intrusion is underestimated by the model at the Zhangpeng station (Fig. 8c) and the Dasheng station (Fig. 8d) in the upstream part of the estuary. A likely reason for the underestimation can be the fundamental assumption that the channel has a constant cross section. The river width convergence at the Humen estuary can actually be described by an exponential function. Simplifying this estuary geometry can result in the underestimation of the mixing coefficient. It indicates that topography is a key driver of the salt intrusion along the Humen estuary.

Figure 8 Comparison between observed and computed salinity concentration over time on 31 January (neap tide) at study sites along the Humen estuary.

Table 4 Calibration results of Song's model.

4.3 Sensitivity analysis

The amplitude of salinity can be described by

$$\begin{array}{}\text{(23)}& \widehat{s}={\stackrel{\mathrm{‾}}{s}}_x\cdot I_{\mathrm{s}},\end{array}$$

where ${\stackrel{\mathrm{‾}}{s}}_x$ is the tide-averaged salinity along the estuary and is a function of the river discharge, i.e. Eq. (12). The parameter I_s is the salinity amplitude coefficient that is defined as

$$\begin{array}{}\text{(24)}& I_{\mathrm{s}}=-{\displaystyle \frac{E{Q}_{\mathrm{f}}}{\mathrm{2}DA}},\end{array}$$

representing the interaction between tides and the river discharge. To investigate the longitudinal salinity distribution and intratidal salinity variation for different discharge and tidal dynamic conditions in the Humen estuary, Eqs. (12) and (23) are used to plot the longitudinal salinity curve and intratidal variation in salinity, respectively. The implemented parameters are the same as shown in Table 3; only the river discharge and the tidal excursion are variable.

Three constant discharge values of 200, 600, and 1800 m³ s⁻¹ are used to evaluate the impact of the river discharge on the salinity variation. The discharge values are chosen because the minimum discharge in the dry season is around 600 m³ s⁻¹ in the Humen estuary, and low salinity can be measured at the Huangpuyou station when the discharge is larger than 1800 m³ s⁻¹. In addition, the discharge in the extreme dry season is set to be 200 m³ s⁻¹. The longitudinal salinity curve can be seen in Fig. 9a. At tidal average conditions, the salt intrusion length becomes smaller when the discharge increases. The steepest salinity gradient can be found at the highest discharge (Q_f=1800 m³ s⁻¹). It is clear from Fig. 9b that the salinity amplitude increases firstly and then decreases as the river discharge increases. This is because during periods of low river discharge (Q_f=200 m³ s⁻¹), the tide-averaged salinity is larger but the salinity amplitude coefficient I_s is smaller, which indicates the weaker interaction between the river flow and the tides. However, the tide-averaged salinity decreases rapidly with the increasing river discharge, as we can see from Fig. 9a, resulting in a smaller amplitude of salinity during periods of high river discharge (Q_f=1800 m³ s⁻¹).

Figure 9 (a) Longitudinal salt intrusion curve at tidal average, considering different river discharge values; (b) intratidal variation in salinity at Huangpuyou station on 31 January 2005, considering different river discharge values.

The tidal effect is studied using three different tidal excursions. The tidal excursion values result in the plots that are shown in Fig. 10. The longitudinal salinity distribution at tidal average conditions is independent of the tidal excursion, as can be seen in Fig. 10a. From Eq. (23), since the salinity amplitude coefficient I_s is in direct proportion to the tidal excursion, the amplitude of the salinity shows a linearly increasing trend with the increased tidal excursion (Fig. 10b).

Figure 10 (a) Longitudinal salt intrusion curve at tidal average, considering different tidal excursion values; (b) intratidal variation in salinity at Huangpuyou station on 31 January 2005, considering different tidal excursion values.

5 Discussion

5.1 Time lag between salinity extremes and slack water

In estuaries, it is noticed that the maximum salinity appears after HWS and the minimum salinity appears before LWS. However, often, the salinity at HWS and LWS corresponds approximately to the maximum and minimum salinity, respectively. The accuracy of this approximation cannot be inferred from existing steady-state models for salt intrusion, as time variation is neglected. As shown in Fig. 11, the unsteady analytical solution proposed in this paper demonstrates that the phase lag between tidal velocity and salinity transportation is π∕2, which means that the extreme values of salinity appear when the tidal velocity is zero. Our unsteady equation for salinity (i.e. Eq. 11) demonstrates the influence of the river discharge on the occurrence of maximum and minimum salinity relative to HWS and LWS, respectively.

Figure 11 Salinity and tidal flow velocity over a tidal cycle at Huangpuyou station. The measured salinity is represented by triangles, and the measured flow velocity is indicated by circles (on 31 January 2005). The dashed line is the calculated tidal velocity, while the dashed–dotted line is the total velocity of tidal flow and river flow. The red solid curve represents salinity simulated by the unsteady analytical solution, which reproduces the time lag HWS and maximum salinity.

More generally, Eq. (14) offers a simple expression yielding qualitative insight into the role of the river discharge in the spatio-temporal variation in salinity in a well-mixed estuary. The time lag between salinity extremes and slack water is determined by the strength of the river flow in a way that is consistent with the previous observations in which the maximum salinity appears after HWS and the minimum salinity appears before LWS. The estimated river flow velocity at the Huangpuyou station is about one-sixth of the tidal flow amplitude, resulting in a time lag between HW (at maximum salinity) and HWS (when total velocity is zero) of less than 30 min. At this station, it is acceptable to assume that the salinity reaches the maximum value at HWS and the minimum value at LWS.

5.2 Optimizing water intake

Estuaries are crucial feeding and breeding grounds for many life forms and are a source of drinking water. Intrusion of salt water can temporarily halt the production of drinking water and put stress on plant and animal species that have adapted to the typical salt concentrations along the estuary. In China, a value of 0.5 ‰ salinity is considered to be the upper limit of drinking water (SWEQ PRC, 2002), while turbot farmed in man-made ponds needs to live in water with no less than 12 ‰ salinity. The unsteady solution proposed in this paper shows reproducing the intratidal variation in salt intrusion, which allows estimating the window of opportunity for drinking-water intake and has the potential of application in aquaculture and water-fetching methods in estuaries.

Due to the serious increase in salt intrusion in recent years, the water intake from Humen estuary is more suitable for saline-water aquaculture than residential use. However, the salinity along the estuarine channel is changing all the time according to the variations in the tides as well as the fresh-water discharge. This makes it important to capture the temporal variation in salinity for optimizing the water intake of the man-made ponds around the estuary. The analytical model proposed in this study provides a simple and efficient approach for predicting the variation in salinity, which is economical and practical, with the limited amount of data available.

Close to the Sishengwei station, there is an aquaculture area with many man-made ponds of different sizes. Optimizing water intake is a key issue here. The applicability of the analytical model is illustrated by focussing on turbot farming, which requires salinity of no less than 12 ‰. The observed salinity data on 29 January are used to calibrate the model, where the determination of three parameters is needed, i.e. tide-averaged salinity at mouth ${\stackrel{\mathrm{‾}}{s}}_{\mathrm{0}}$, the slope $k=\left(Q_{\mathrm{f}}a/D{A}_{\mathrm{0}}\right)$, and the tidal excursion E. The decreasing trend of subtidal salinity is close to a linear relation from 29 January to 3 February (Fig. 12b). Thus the tide-averaged salinity value of the predicted model is set as 90 % of that on 29 January, considering the slight change of the subtidal salinity in the 5 d after 29 January. Moreover, the slope and the tidal excursion are assumed to be constant during the whole period from 29 January to 3 February. As shown in Fig. 12c, the prediction by the model is in good agreement with the observation in this case. Furthermore, if more observed data are available to calibrate the tide-averaged salinity covering the period from 29 January to 3 February, Eq. (14) performs better, as shown in Fig. 12d. The available time for water intake can be obtained from the prediction, when the salinity concentration reaches a value higher than 12 ‰.

Figure 12 Time for water intake of given salinity that is higher than 12 ‰. (a) Slight changes of the subtidal discharge. (b) Decreasing trend of subtidal salinity. (c) Predicted salinity on the basis of observed data on 29 January. (d) Calibrated salinity on the basis of observed data from 29 January to 3 February 2005.

Since the fresh-water discharge influences the slope (Brockway et al., 2012), it is reasonable to assume that the slope remains constant in a short timescale, since the fresh-water discharge variation has a timescale of days to months (Fig. 12a). The tidal excursion is the integral over time of the tidal velocity between the low water slack and high water slack. It varies from day to day as the tidal wave changes from neap tide to spring tide (Savenije, 2005). Therefore, the tidal excursion is assumed to be independent of time in the neap cycle from 29 January to 3 February. Besides, Eq. (18) is demonstrated to be a useful equation for the calculation of the tidal excursion, which offers an approach for estimating the tidal excursion with salinity data. The predicted salinity fits well with observed values, indicating that the estimation of the tide-averaged salinity during the neap tide is acceptable. However, the prediction accuracy of the model can be higher if more observed tide-averaged salinity data are available.

6 Conclusions

An unsteady-state analytical solution of salt intrusion is proposed based on the one-dimensional advection–diffusion equation for salinity, assuming a harmonic tidal wave with a single-frequency and a constant mixing coefficient. The predictive skill of the model has been illustrated from an application to the Humen estuary, which shows that it can offer an efficient approach for calculating the variation in salinity in a well-mixed estuary where the channel area is convergent. The results show that the analytical model is able to reproduce the intratidal variation in salt intrusion and can be a useful tool for computing the time windows in which salinity remains below a critical threshold in an estuary.