The physics behind Van der Burgh’s empirical equation, providing a new predictive equation for salinity intrusion in estuaries

The practical value of the surprisingly simple Van der Burgh equation in predicting saline water intrusion in alluvial estuaries is well documented, but the physical foundation of the equation is still weak. In this paper we provide a connection between the empirical equation and the theoretical literature, leading to a theoretical range of Van der Burgh’s coefficient of 1/2<K < 2/3 for density-driven mixing which falls within the feasible range of 0<K < 1. In addition, we developed a one-dimensional predictive equation for the dispersion of salinity as a function of local hydraulic parameters that can vary along the estuary axis, including mixing due to tide-driven residual circulation. This type of mixing is relevant in the wider part of alluvial estuaries where preferential ebb and flood channels appear. Subsequently, this dispersion equation is combined with the salt balance equation to obtain a new predictive analytical equation for the longitudinal salinity distribution. Finally, the new equation was tested and applied to a large database of observations in alluvial estuaries, whereby the calibratedK values appeared to correspond well to the theoretical range.


Introduction
Estuaries play an essential role in the human-earth system, affecting fresh water resources, the mixing between ocean and river water, and the health of aquatic ecosystems. This makes the functioning of estuarine systems an important field of research. A crucial element of estuarine dynamics is the interaction between saline and fresh water. The river discharges fresh water into estuaries, flushing out the salt, while saline water penetrates landward as a result of density gra-dients. The temporal and spatial distribution of salinity in an estuary is determined by the competition between fresh water flushing and penetration of saline water by gravity.
Dispersion is the mathematical reflection of the spreading of a substance (e.g., salinity s) through a fluid as a function of a gradient in the concentration of the substance (e.g., the salinity gradient ds / dx). Hence, dispersion is the mathematical description of mixing. The physical process driving dispersion differs at different scales, depending on the dominant mechanism. For instance, at the molecular scale, the dominant mechanism is the Brownian movement of water molecules. At the scale of river flow, the process is driven by the transfer of friction from the riverbed into the cross section through turbulence. At this scale, the dispersion coefficient is called hydraulic eddy viscosity (K E ) (Dyer, 1997). The most important type of mixing in estuaries is the result of salinity gradients and the non-concurrence of the velocity and salinity field (u s ) (MacCready, 2004), which is the result of gravitational and tidal mixing processes. Finally, there is mixing by residual circulation, driven by the tide, where ebb and flood flows of different densities mix (e.g., Nguyen et al., 2008).
The dispersion resulting from density gradients is closely connected to the stratification number N R , which is the balance between the potential energy resulting from the buoyancy of fresh water flowing into the estuary and the kinetic energy of the tide that provides the energy of mixing. This stratification number, also known as the estuarine Richardson number, is widely used in theoretical and practical studies (e.g., Fischer, 1972;Savenije, 1986;Kuijper and Van Rijn, 2011). If N R is large, potential energy of river discharge dominates and stratification occurs; if N R is small, the estuary is 2 Linking Van der Burgh to the traditional literature The one-dimensional mass-conservation equation averaged over the cross section and over a tidal cycle can be written as (e.g., Savenije, 2005) where A = Bh is the cross-sectional area, B is the width, h is the depth, s is the cross-sectional average salinity, t is time, Q f is the fresh water discharge, x is the distance from the estuary mouth and D is the effective longitudinal dispersion coefficient. The positive direction of flow is in the upstream direction. At steady state, where ∂s/∂t = 0, using the boundary condition at x → ∞, s = s f and ∂s/∂x = 0, integration yields where s f is the fresh water salinity, usually close to zero. Van der Burgh (1972) found an empirical equation describing the tidal average longitudinal variation of the effective dispersion: where the dimensionless coefficient K ∈ (0, 1) according to Savenije (2005). Combining Eqs.
(2) and (3) yields (Savenije, 1986(Savenije, , 1989(Savenije, , 1993a where D 1 and s 1 are the dispersion coefficient and salinity at a given point x 1 , generally taken at the inflection point in the exponential estuary geometry. This equation is special in that it links the dispersion to the salinity instead of the salinity gradient, as most other researchers do (e.g., Fischer, 1976;Prandle, 1981;Thatcher and Najarian, 1983). Interestingly, using Eqs. (2) and (4) we can derive the dispersion as a function of the salinity gradient (Savenije, 2015): which connects the dispersion coefficient to local variables (A, ds / dx), boundary conditions (D 1 , s 1 ) and K.

MacCready
where u s is the tidal average and width average exchange flow salt flux, u is the depth-varying velocity, s is the depth-varying salinity, K H is the along-channel diffusion coefficient, m 1 = 2 105 , m 2 = 19 420×48 and m 3 = 19 630×48 2 are constant values following MacCready's vertical integration, u f = |Q f |/A is the depth-averaged velocity of fresh water, K S is the effective vertical eddy diffusivity, g is the gravity acceleration, c s is the saline expansivity equal to 7.7 × 10 −4 , and K E is the effective hydraulic eddy viscosity. For the latter, we use the equation K E = 0.1 2 π u * h, with u * = √ g C υ as the shear velocity in relation to the tidal velocity amplitude υ (= πE T ; E is tidal excursion length; T is tidal period), where C = K m h 1/6 is the coefficient of Chézy, and K m is Manning's coefficient. Comparing the salt balance equation of MacCready to Eq. (2) implies that Eq. (6) is identical to −D ds dx . MacCready assumed the estuary to be narrow and rectangular, in the sense that cross-sectional shape does not basically modify the width-averaged dynamics. In the derivation, he also assumed the effective vertical eddy viscosity to be constant with depth, following Hansen and Rattray (1965), and that the salinity gradient of the depthvarying part is much smaller than the depth-averaged part, following Pritchard (1952). Additionally, other effects like salt storage, internal hydraulics and the Coriolis force were considered negligible.
After division of all terms by the salinity gradient, it becomes an equation for the dispersion coefficient D: whereby the first term is not dependent on the salinity gradient, the second is directly proportional to it, and the third term depends on the square of the salinity gradient. Based on Eq. (5) we can also derive an expression for the dispersion: where A 1 is the cross-sectional area at the inflection point (at x = x 1 ), l = L−x 1 is the distance from the inflection point to where salinity becomes the same as the fresh water salinity, and L is the total intrusion length.
Given the function F (γ ) = γ K 1−K , a Taylor series expansion near γ = 1 can be derived as where R 2 (x) is the residual term, considered to be small. To analyze the importance of the different terms of Eq. (9), Fig. 1 presents the factors g 1 = (2K−1)(3K−2) 2(1−K) 2 , g 2 = K(2−3K) (1−K) 2 and g 3 = K(2K−1) 2(1−K) 2 . g 1 is the closure term which compensates for g 2 and g 3 so as to make g i = 1 (i = 1, 2, 3). It is clear that the absolute value of the first term is much smaller than the density-driven terms. Also, the larger the value of K, the more important the third term is. This is in accordance with traditional literature. If K = 1/2, F (γ ) = A A 1 l s 1 − ds dx , and dispersion is proportional to the salinity gradient. If , then dispersion is proportional to the square of the salinity gradient, which means that the dispersion is more sensitive to the salinity gradient.
Considering only the density-dependent terms in Eqs. (7) and (9), the proportionality results in leading to an analytical expression for K: According to Eq. (10), K is not time independent, as was previously assumed by Savenije (2012); rather, it is deter-mined by the tidal excursion and the fresh water discharge. In the case of a relatively constant discharge, a larger tidal excursion implies less stratification, a larger value of w and K approaching the lower limit (1/2). On the other hand, a smaller tidal excursion implies more stratification, a smaller value of w and K approaching the higher limit (2/3), which corresponds to the situation where the dispersion is more sensitive to the salinity gradient. We have used this expression to compute K values in 18 real estuaries using the database of Savenije (2012). These K values are in a range of 0.51-0.64 (see Sect. 4.3).
Overall, there are three results for the estimation of Van der Burgh's coefficient: (1) by comparison with traditional studies (K = 1/2 or K = 2/3), (2) by comparison with Mac-Cready considering the salinity relevant terms (1/2 < K < 2/3), and (3) based on empirical calibration (see Sect. 4). These results are surprisingly close, even though the theoretical comparison is limited to density-driven mixing.

Including residual circulation in wide estuaries
In the theory about mixing in estuaries, several authors have distinguished between tide-driven and density-driven dispersion (e.g., Hansen and Rattray, 1965;Banas et al., 2004). The tide is an active hydraulic driver that creates shear stresses in the flow as momentum, resulting from friction along the boundaries, transferred to the heart of the channel by turbulence. Generally these shear stresses reduce stratification and hence reduce dispersion. However, at a large scale, the tide facilitates mixing by tidal trapping and residual circulation, which enhances dispersion. Tidal trapping results from irregularities of the channel, leading to pockets of relatively high or low salt concentrations that later reunite with the stream. The mixing length scale of tidal trapping is the tidal excursion. By using the tidal excursion as the mixing length, tidal trapping can be incorporated into a predictive equation. Residual circulation is more complicated. It can be a very powerful tide-driven mechanism in the wider parts of estuaries where the tide causes mixing by the cross-over of preferential ebb and flood channels that develop in wide estuaries, such as the Schelde, described by Nguyen et al. (2008). But how can we parameterize residual circulation? Here a different approach is followed from Nguyen et al. (2008), trying to combine this effect in the regular one-dimensional advection-dispersion equation. be described as

Model including residual circulation
where V = AE is the water volume, s i is the salinity at different locations i, and d and r are longitudinal and lateral exchange flows. The balance equation then becomes where x and y are the mixing lengths, which are taken as x = E and y = B/2.
The assumption used is that the lateral exchange is proportional to the longitudinal (Fischer, 1972): As a result, longitudinal and lateral processes can be combined into one single one-dimensional equation: Comparing Eq. (15) with the traditional salt balance equation, the effective longitudinal dispersion is Subsequently, the longitudinal exchange flow d is assumed to be proportional to the amplitude of the tidal flow (driving the circulation) ( Q t = Aυ), and to the stratification number to the power of K: with N R defined as the ratio of potential energy of the river discharge to the kinetic energy of the tide over a tidal period: Hydrol. Earth Syst. Sci., 21, 3287-3305, 2017 www.hydrol-earth-syst-sci.net/21/3287/2017/ where ρ/ρ = c s s is the relative density difference between river water and saline water. The reason why the exchange flow is a function of the stratification number to the power of K is because it is in agreement with Eq. (4), ρ/ρ being directly proportional to s.
We then obtain a simple dimensionless expression for the dispersion coefficient, simulating to the one by Gisen et al. (2015b) but incorporating lateral exchange flow: where C 1 and C 2 are constants.

Analytical solution
In almost all estuaries, the ratio of width to excursion length is quite small, particularly upstream where salinity intrusion happens. So for further analytical solutions we can focus on the first part of Eq. (19): The traditional approach by Savenije (2012) merely uses this equation as the boundary condition at x = x 1 , after which D(x) values are obtained by integration of Van der Burgh's equation along the estuary axis. But, in principle, with this equation the dispersion can be calculated at any point along the estuary, provided local hydraulic and geometric variables are known. Equation (20) can be elaborated into where all local variables are now a function of x.
The following equations are used for the tidal velocity amplitude, width and tidal excursion: where δ υ ≈ δ H are the damping/amplifying rate of tidal velocity amplitude and tidal range, and b is the width convergence length (b 1 downstream of the inflection point and b 2 upstream). At the inflection point, the predicted equation is given by where the subscript "1" means parameters are evaluated at the inflection point (x = x 1 ). Substitution of Eqs. (22)-(25) into Eq. (21) gives Differentiating D with respect to x and using Eq. (26) results in Combining the result with the time-averaged salt balance, Eq. (27) For a prismatic channel (b → ∞) with constant width and little tidal damping, = 0 and Eq. (28) becomes Van der Burgh's equation. As a result, the exponent of N R in this model represents Van der Burgh's coefficient.
The cross-sectional area A is given by where a is the cross-sectional convergence length (a 1 downstream of the inflection point and a 2 upstream). Substitution of Eq. (29) into Eq. (28) gives In analogy with Kuijper and Van Rijn (2011), the solution of this linear differential equation is with ζ = a 1− a . The maximum salinity intrusion length is obtained from Eq. (31) after substitution of D → 0 at x = L: This is the same equation as in Savenije (2005) if ζ = a. Using Eq. (26), the longitudinal salt distribution becomes This solution is similar to the solution by Kuijper and Van Rijn (2011), with the difference that Kuijper and Van Rijn used a constant value of K = 0.5 and that their value of depended on the bottom slope.
So with these new analytical equations, the local dispersion and salinity can be obtained, using the boundary condition at the inflection point. This method is limited since it only works when B/E < 1. If we want to account for residual circulation using Eq. (19), then we have to use numerical integration of Eq. (2) using Eq. (19) for D. Distance, x (m) Figure 3. Semi-logarithmic presentation of estuary geometry, comparing simulated (lines) to observations (symbols), including crosssectional area (green), width (red) and depth (blue).

Summary information
Eighteen estuaries with quite different characteristics, covering a diversity of sizes, shapes and locations, have been selected from the database of Savenije (2012). It appears that all these alluvial estuaries can be schematized in one or two segments separated by a well-defined inflection point (Savenije, 2015). As an example, Fig. 3 shows the geometry of two estuaries: Maputo with an inflection point and the Thames without an inflection point. It can be seen that the natural geometry fits well on semilogarithmic paper, indicating an exponential variation of the cross section and width. Geometric data of all 18 estuaries are presented in Appendix A.
In Table 1 the general geometry of estuaries is summarized, where B f is the bankfull stream width. It is obvious that these estuaries cover a wide range of sizes. An estuary with x 1 = 0 means there is no inflection point. In addition, the larger the convergence length a 2 (b 2 ), the slower the cross section (width) declines upstream. With a large b 2 , a relatively small value of b 1 suggests the channel with a pronounced funnel shape with fast decrease in width near the mouth. In contrast, a relatively large value of b 2 indicates estuaries with near-prismatic shape. The same values of a and b indicate that the depth is constant.
Tables 2 and 3 contain summary information of estuaries on different measurement dates, where H 1 is the tidal range at x 1 , η is tidal amplitude, α = D 1 |Q f | is the mixing coefficient and β = Ka 2 |Q f | A 1 D 1 is the dispersion reduction ratio. Tidal excursion and tidal period are more or less the same in all estuaries, except for Lalang and Chao Phraya with a diurnal tide. Most estuaries damp upstream, with negative values of δ H . In addition, most estuaries have a small tidal amplitude to depth ratio, which means relatively simple solutions of hydraulic equations are possible (Savenije, 2005). K values have been obtained by calibration of simulated salinities to observations in 18 estuaries. According to Eq. (4), the K value affects the salinity mostly in the upstream reach, where D/D 1 is small. Using an automatic solver, the best result was obtained with C 1 = 0.10, C 2 = 12 and K = 0.58. For individual es-tuaries, K values were obtained ranging between 0.45 and 0.78. The dispersion at the inflection point has a range of 50-600 m 2 s −1 in a variety of estuaries, which is consistent with Prandle (1981). The mixing coefficient demonstrates to what extent the dispersion overcomes the flushing by river flow. The larger the river discharge, the smaller the α, meaning it is difficult for the salinity to penetrate into the estuary. The dispersion reduction ratio determines the longitudinal variation of dispersion. Fischer et al. (1979) suggested that the transition from a well-mixed to a strongly stratified estuary occurs when the values of stratification number N R are in the range of 0.08-0.8. With a ratio of π between Fischer's and our expressions for the stratification number, the range becomes 0.25-2.51. It is obvious that all estuaries are partially to well mixed, with N R below 2.51.

Sensitivity to C 2
Through the use of C 2 we can use a single dispersion equation accounting for two-dimensional effects in a onedimensional model. The assumption that lateral exchange is proportional to longitudinal dispersion suggests C 2 to be independent of x. Figure 4 and Appendix B demonstrate how salinity changes with varying C 2 . Salinities were simulated by numerical solution of Eq. (2) with Eq. (19) based on the boundary condition at x = x 1 . Typically, C 2 matters mainly near the mouth, but there is almost no effect on narrow estuaries like Lalang, Limpopo, Tha Chin and Chao Phraya. Hence, the inclusion of the residual circulation improves the accuracy of salinity simulation in wide estuaries and more particularly near the mouth of the estuaries where the ratio of width to tidal excursion is relatively large.
To check the sensitivity to C 2 , values of 1, 10 and 50 have been used to calculate salinity curves. It is demonstrated that the larger the value of C 2 , the smaller the salinity gradient and the flatter the salinity curve near the estuary mouth. However, because of the interdependence of D, s and ds/dx through Eq. (2) in the upstream part, a larger value of C 2 can lead to larger salinities (e.g., Thames, Elbe, Edisto, Maputo and Corantijn). Basically, C 2 = 10 (green lines) can perform perfectly in 14 out of 18 estuaries (e.g., Maputo and the Thames). We can see that larger values than C 2 = 10 cause exaggerated salinity in the downstream part of these estuaries, which is why a general value of C 2 = 10 is recommended. The poorer results occur in estuaries that have peculiar shapes near the mouth. A larger value of C 2 applies to the Kurau. This may be because the width is underestimated in the wide estuary mouth, due to misinterpretation of the direction of the streamline (the width is determined according to a line perpendicular to the streamline). As a matter of fact, the width should be larger and dispersion should be larger with smaller salinity gradients, which would then result in a lower value of C 2 . The same applies to Endau. By contrast, a smaller value of C 2 in Perak fits better, because of overestimation of the width. Here the topographical map suggests a  ( Figure 4. Comparison between simulated and observed salinity at high water slack (thin lines) and low water slack (thick lines), scaled by the salinity s 1 at the inflection point x 1 for different C 2 values. Observations at high water slack are represented by triangles and low water slack by circles. Observe that the Thames only has low water slack observations. wider estuary mouth, whereas the tidal flow is concentrated in a much narrower main channel due to the north bank protruding into the estuary and a spur from the south projecting into the mouth. The Selangor has a similar situation. It shows that the configuration of the mouth is important for the cor-rect simulation of the salinity near the estuary mouth. But, fortunately, a relatively poor performance near the mouth of these estuaries does not affect the salinity distribution upstream as long as C 2 is not too large. In conclusion, C 2 = 10 appears to be a suitable default value as long as the trajectory of the tidal currents can be considered properly. The poor fit in the downstream parts of the Lalang and Chao Phraya, in which measured salinities are lower than simulated, can be explained by a complex downstream boundary. The Lalang estuary has a pronounced riverine character and is a tributary to the complex estuary system of the Banyuasin, sharing its outfall with the large Musi River.
So the salinity near its mouth is largely affected by the Musi. Also, pockets of fresh water can decrease the salinity near the confluence. The Chao Phraya opens to the Gulf of Thailand where the salinity is influenced by historical discharges rather than ocean salinity, remaining relatively fresh. Other measurement uncertainties may cause outliers as well.

A possible solution for K
The physical meaning of Van der Burgh's coefficient has been analyzed, linking it to traditional theoretical research. Equation (10) shows a direct relation between this coefficient and MacCready's parameters, which are measurable quantities. Hence, the coefficient is affected by tide, geometry and fresh water discharge. Shaha and Cho (2011) also found K values to depend on river discharge and considered its value to increase upstream in a range of 0-1 due to different mechanisms along the estuary. A 1 : 1 plot is presented in Fig. 5, relating the empirical K values to the predicted values using Eq. (11). The predicted K values have a smaller range (0.51-0.64) than the calibrated ones (0.45-0.78). Moreover, it can be seen that there is a steep linear relation between predictive and calibrated K values, which reveals that the predictive method overestimates the low calibrated K values and underestimates the high val- ues. Fully tide-driven dispersion would correspond to K = 0. But the predictive method does not consider tidal mixing and, as a result, the predicted K values are too high in the lower region. Smaller calibrated values imply that the tide plays a prominent role in the estuary. For the higher calibrated values, another explanation applies. The theoretical approach follows width-averaged dynamics, whereas the empirical approach relies on natural estuaries with cross-sectional variations. A K value larger than the predicted value could result from a strong lateral salinity gradient due to shearing in a complex geometry, which strengthens the sensitivity to the salinity gradient. In addition, there is quite some uncertainty in calibrating a partly empirical analytical model to data in real estuaries, as a result of a whole range of uncertain factors related with observational errors, data problems, the assumption of steady state and other factors. Some estuaries may be in non-steady state (e.g., the Thames). However, considering the K values have been obtained from different approaches, they are still quite similar. As a result, this correspondence forms, at least partly, a physical basis for the Van der Burgh coefficient. All K values are very close to 0.58, which may be a good starting value in estuaries where information on geometry and channel roughness is lacking.

Discussion and conclusion
Overall, the single one-dimensional salinity intrusion model including residual circulation appears to work well in natural estuaries with a diversity of geometric and tidal characteristics, by both analytical and numerical computation. The new equation is a simple and useful tool for analyzing local dispersion and salinity directly on the basis of local hydraulic variables. In a calibration mode, K is the only parameter to be calibrated using C 1 = 0.10 and C 2 = 10. In a predictive mode, a value of K = 0.58 can be used as a first estimate. If information on river discharge, roughness and geometry is available, K can be determined iteratively by taking K = 0.58 as the predictor and subsequently substituting s 1 and l from the first iteration by Eqs. (10) and (11) and repeating the procedure until the process converges.
The addition of the factor (1+C 2 (B/E) 2 ) in the dispersion equation proved valuable near the mouth of estuaries where residual circulation due to interacting ebb and flood channels dominates dispersion. The value C 2 = 10 was found to perform best in most estuaries, indicating that residual circulation is dominant in wide estuaries where ebb and flood currents prevail.
Van der Burgh's coefficient determines the way dispersion relates to the stratification number by a power function. Two approaches, theoretical derivation from the traditional literature and empirical validation based on observations in a large set of estuaries, provided similar estimates of Van der Burgh's coefficient. Under MacCready's assumptions, there are three ways to estimate K: 0.51 < K < 0.64 from empirical application of Eqs. (10) and (11); 1/2 < K < 2/3 as the physical boundaries of Eq. (11); and the comparison with traditional approximations (K = 1/2 or K = 2/3). After calibration of the new analytical model to the database of field observations, the values of K were in a range of 0.45-0.78 for a wide range of conditions, with an average of 0.58, close to the predicted values. MacCready's equation determines dispersion via a decomposition method, using depth-varying velocity and salinity. Although these 1-D expressions of velocity and salinity may be simplifications of reality, the good correspondence between Van der Burgh's equation and Mac-Cready's theory provides a promising theoretical basis for Van der Burgh's equation.
A previous analytical salinity intrusion model was developed by Gisen (2015a), from which the K values resulted in a range of 0.20-0.75 by calibration and 0.22-0.71 by prediction. These solutions cover a wider range than our estimates because of Gisen's assumption that K does not depend on river discharge and because of three improvements made in this paper. Firstly, we used the local hydraulic parameters to simulate the salinities, while Gisen used a constant depth and no damping ( = 0). In addition, by using an uncertainty bound of 25 % on fresh water discharge we could reduce the inaccuracy of the tail of the salinity curve and obtain a better fit (where K matters most). And finally, all geometric analyses were improved by revisiting the fit to observations. An important consequence of this research is that K depends on time and space. Where Gisen assumed K to be constant for each estuary, we find substantial variability for estuaries where a larger range of discharges is available: e.g. in the Maputo 0.57 < K < 0.70; in the Limpopo 0.61 < K < 0.72; and in the Edisto 0.48 < K < 0.58. The implication of discharge dependence needs to be tested further for predictive purposes.
Hydrol. Earth Syst. Sci., 21, 3287-3305, 2017 www.hydrol-earth-syst-sci.net/21/3287/2017/ In some particular cases, the simulated salinity with C 2 = 10 does not fit the observations near the estuary mouth. So one should be aware of peculiar configurations of streamlines and geometries near the estuary mouth when using this model. MacCready and Geyer (2010) also pointed out that the effect of irregular channel shape is important. However, a poor fit near the estuary mouth has almost no effect on the total salinity intrusion length. It is suggested that in future research, the assumption that lateral exchange is proportional to longitudinal exchange needs to be tested further. Finally, this predictive one-dimensional salinity intrusion model, having a stronger theoretical basis, may be a useful tool in ungauged estuaries. Hydrol. Earth Syst. Sci., 21, 3287-3305, 2017 www.hydrol-earth-syst-sci.net/21/3287/2017/