Maximum power of saline and fresh water mixing in estuaries

According to Kleidon (2016), natural systems evolve towards a state of maximum power, leading to higher levels of entropy production by different mechanisms, including gravitational circulation in alluvial estuaries. Gravitational circulation is driven by the potential energy of fresh water. Due to the density difference between seawater and river water, the water level on the riverside is higher. The hydrostatic forces on both sides are equal but have different lines of action. This triggers an angular moment, providing rotational kinetic energy to the system, part of which drives mixing by gravitational circulation, lifting up heavier saline water from the bottom and pushing down relatively fresh water from the surface against gravity; the remainder is dissipated by friction while mixing. With a constant freshwater discharge over a tidal cycle, it is assumed that the gravitational circulation in the estuarine system performs work at maximum power. This rotational flow causes the spread of salinity inland, which is mathematically represented by the dispersion coefficient. In this paper, a new equation is derived for the dispersion coefficient related to density-driven mixing, also called gravitational circulation. Together with the steady-state advection–dispersion equation, this results in a new analytical model for density-driven salinity intrusion. The simulated longitudinal salinity profiles have been confronted with observations in a myriad of estuaries worldwide. It shows that the performance is promising in 18 out of 23 estuaries that have relatively large convergence length. Finally, a predictive equation is presented to estimate the dispersion coefficient at the downstream boundary. Overall, the maximum power concept has provided a new physically based alternative for existing empirical descriptions of the dispersion coefficient for gravitational circulation in alluvial estuaries.

1 Introduction

Estuaries are water bodies in which rivers with fresh water meet the open sea. The longitudinal salinity difference causes a water level gradient along the estuary. As a result, the water level at the limit of salt intrusion is set up several centimeters above sea level (about 0.012 times the estuary depth). Therefore, the hydrostatic forces from the seaside and riverside have different lines of action (a third of the setup apart). Since the hydrostatic forces at the seaside and the salinity limit are equal but opposed, this difference in the lines of action triggers an angular moment (a torque) that drives the gravitational circulation, whereby fresh water near the surface flows to the sea and saline water near the bottom moves upstream (Savenije, 2005). This density-driven gravitational circulation is one of the two most significant mixing mechanisms in alluvial estuaries; the other is the tide-driven mixing mechanism that can be dominant in the wider (downstream) part of estuaries (Fischer et al., 1979).

Kleidon (2016) described the concept of maximum power in the Earth system, implying that freely evolving systems perform work and dissipate energy at maximum power (close to or at the Carnot limit). Using this concept, gravitational circulation is assumed to take place at the maximum power limit. Earlier, the maximum power concept was used to solve saline and fresh water mixing as in a thermodynamic equilibrium system (Zhang and Savenije, 2018). It assumed that in thermodynamic terms, the freshwater flux maintains a potential energy gradient, triggering fresh and saline water mixing processes that work at depleting this gradient. Because the strength of the mixing in turn depends on this gradient, there is an optimum at which the mixing process performs at maximum power. It did not, however, account for the energy loss associated with this mixing process. The equation obtained appeared to have an analytical solution of a straight line for the longitudinal salinity distribution. Although this is not correct, it can be seen as a first order approximation, which agrees with earlier theoretical work by Hansen and Rattray (1965), who developed their theory for the central region of the salt intrusion length at which the salinity gradient is at its maximum and dominated by density-driven mixing. However, this approximate solution was not fully satisfactory for simulating the salinity distribution.

In contrast to the earlier work by Zhang and Savenije (2018), in this paper friction is taken into account. The available free energy by the angular moment is converted into work (mixing the saline and fresh water against the force of gravity) and the associated frictional dissipation. In the following sections we shall derive a new equation for density-driven mixing, which appears to compare well with observations in a range of alluvial estuaries.

Kleidon (2016) presented several examples for the application of the maximum power limit on nonthermal energy conversions. In one example, a fluid is kept in motion by an accelerating force that provides kinetic energy to the system. The velocity of the fluid is slowed down by friction and the remainder is converted into another form of energy. If the velocity is too large, the friction is large and energy dissipation dominates, then the power of the force to generate work is limited. In contrast, if the velocity is too small, the power is not enough to generate work. Hence, there is an optimum value for the product of the force and velocity to produce maximum useful energy. Estuaries are comparable to this system. In this article, we apply the maximum power concept to gravitation circulation generated by a longitudinal density gradient.

Traditionally, the empirical Van der Burgh (VDB) method has worked very well to describe the mixing in alluvial estuaries, leading to predictive equations to describe the salinity intrusion in alluvial estuaries (Savenije, 2005, 2012). The VDB method takes account of all mixing mechanisms, including density-driven (gravitational) circulation and tide-driven mixing. For application of the VDB method, there are two parameters that need to be calibrated, the empirical Van der Burgh coefficient K and the dispersion coefficient at the downstream boundary D₀. This method has performed surprisingly well around the world and has been used in this paper as the benchmark model for comparison with the maximum power approach.

2 Moment balance for an open estuary system

In an estuary, the cross-sectional average hydrostatic forces have equal values along the estuary axis. Over a segment, the forces are opposed but working on different lines of action due to the density gradient in the upstream and downstream directions. As a result, they exert an angular moment (torque) M_acc that drives the gravitational circulation, performing as accelerating torque. The velocity of the gravitational circulation kept in motion by this accelerating torque is slowed down by a friction moment M_fric, which is the product of the associated friction force and its arm. The remainder M_ex drives the circulation and executes work against gravity (Fig. 1). Hence, the balance in steady state in a segment is

$$\begin{array}{}\text{(1)}& M_{\mathrm{acc}}-M_{\mathrm{ex}}-M_{\mathrm{fric}}=\mathrm{0}.\end{array}$$

The moment due to the friction against the circulation is expressed as

$$\begin{array}{}\text{(2)}& M_{\mathrm{fric}}=F_{\mathrm{fric}}{l}_m,\end{array}$$

with F_fric being the friction force (N) and l_m the scale of the arm of the frictional forces in meters.

Figure 1 (a) Systematic salt transport in estuaries, with the seaside on the left and the riverside on the right. The water level (in blue) has a slope as a result of the salinity distributions (in red). The hydrostatic forces on both sides have different lines of action that trigger the gravitational circulation, providing an accelerating moment M_acc to the system. (b) A box model displaying the moment balance in open estuarine systems.

The friction force during the dispersive circulation is expressed as

$$\begin{array}{}\text{(3)}& F_{\mathrm{fric}}=\mathit{\tau }O,\end{array}$$

where τ is the shear stress (N m⁻²) and O is the contact area (m²). Estuarine mixing has two length scales: a vertical and a horizontal one. The horizontal length scale is the tidal excursion E, which is the distance a water particle travels on the tide; the vertical length scale is the depth h, over which saline water is moved upward to the surface and over which relatively fresh water is moved downward to the bottom. Since the process of gravitational mixing is essentially to move the saline water up and the fresher water down, the contact area for the resistance against this movement is determined by the depth (h) and the width (B). Following that reasoning, O is assumed equal to Bh. Meanwhile, the circulation cell has a dimension constrained by the depth. The circular movement hence has a diameter of the depth and l_m, the horizontal arm between the vertical frictional forces, is of the order of magnitude of the depth.

2.1 Maximum power condition in estuaries

Because the velocity of the dispersive gravitational circulation is small, the mixing flow is assumed to be laminar. The shear stress is typically a function of flow velocity (v): τ=ρqv, with ρ being the density (kg m⁻³) and q being a laminar resistance (m s⁻¹). The latter is assumed to be proportional to the tidal velocity amplitude ($q\propto E/T$), where T is the tidal period in seconds. Hence, the flow velocity representing the gravitational circulation is

$$\begin{array}{}\text{(4)}& v={\displaystyle \frac{M_{\mathrm{acc}}-M_{\mathrm{ex}}}{\mathit{\rho }qB{h}^{\mathrm{2}}}}.\end{array}$$

Power is defined by the product of a force and its velocity. The power of torque (angular moment) is defined as the product of the moment and its angular velocity. Hence, the power is defined as

$$\begin{array}{}\text{(5)}& P=M_{\mathrm{ex}}\mathit{\omega }=M_{\mathrm{ex}}{\displaystyle \frac{v}{h/\mathrm{2}}}\mathrm{2}\mathit{\pi }={\displaystyle \frac{\mathit{\pi }}{\mathit{\rho }qB{h}^{\mathrm{3}}}}\left(M_{\mathrm{acc}}-M_{\mathrm{ex}}\right)M_{\mathrm{ex}},\end{array}$$

where ω is the angular velocity or the rotational speed (s⁻¹). Figure 2 illustrates how the execution moment and the flow velocity vary. If the working moment is too large and causes fast mixing flow, the energy dissipation is large and diminishes the flow velocity. If it is too small, the mixing would also be suboptimal. In analogy with Kleidon (2016), the product of the working moment and the flow velocity has a well-defined maximum. The maximum power (MP) is then obtained by $\partial P/\partial M_{\mathrm{ex}}=\mathrm{0}$. Hence, the optimum values of the execution moment M_ex,opt and the flow velocity v_opt are

$$\begin{array}{}\text{(6)}& M_{\mathrm{ex},\mathrm{opt}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{M}_{\mathrm{acc}}\end{array}$$

and

$$\begin{array}{}\text{(7)}& v_{\mathrm{opt}}={\displaystyle \frac{M_{\mathrm{acc}}}{\mathrm{2}\mathit{\rho }qB{h}^{\mathrm{2}}}}.\end{array}$$

Figure 2 Sketch of the sensitivity of the exchange flow velocity v to the working moment M_ex.

Here, the accelerating force (F_acc) that produces the angular moment is the hydrostatic force obtained by integrating the hydraulic pressure over the depth:

$$\begin{array}{}\text{(8)}& F_{\mathrm{acc}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\rho }}_{\mathrm{0}}g{h}^{\mathrm{2}}B,\end{array}$$

where ρ₀ is the density of the seaside (kg m⁻³).

The accelerating moment has an arm Δh∕3 (Savenije, 2005). The water level gradient according to the balance of the hydrostatic pressures results in

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{\mathrm{d}h}{\mathrm{d}x}}=-{\displaystyle \frac{h}{\mathrm{2}\mathit{\rho }}}{\displaystyle \frac{\mathrm{d}\mathit{\rho }}{\mathrm{d}x}},\end{array}$$

where x is the distance in meters. Density is a function of salinity (S; psu): $\mathit{\rho }={\mathit{\rho }}_{\mathrm{f}}(\mathrm{1}+c_{\mathrm{S}}S)$, where ρ_f is the density of the fresh water (kg m⁻³) and c_S ($\approx \mathrm{7.8}\times {\mathrm{10}}^{-\mathrm{4}}$) is the saline expansivity (psu⁻¹).

Subsequently, the accelerating moment due to the density gradient driving gravitational circulation over a tidal cycle can be described as

$$\begin{array}{}\text{(10)}& M_{\mathrm{acc}}=F_{\mathrm{acc}}{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}{\displaystyle \frac{\mathrm{d}h}{\mathrm{d}x}}E=-{\displaystyle \frac{\mathrm{1}}{\mathrm{12}}}{\mathit{\rho }}_{\mathrm{0}}g{h}^{\mathrm{3}}B{c}_{\mathrm{S}}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}x}}E,\end{array}$$

where E is the horizontal length scale of the gravitational circulation in meters.

In steady state, the one-dimensional advection–dispersion equation averaged over the cross section and over a tidal cycle reads (Savenije, 2005, 2012)

$$\begin{array}{}\text{(11)}& \left|Q\right|S+AD{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}x}}=\mathrm{0},\end{array}$$

where Q is the freshwater discharge (m³ s⁻¹), A (=Bh) is the cross-sectional area (m²), and D is the dispersion coefficient (m² s⁻¹). The positive direction of flow is in the upstream direction.

Accordingly, with $q\propto E/T$, the optimum velocity is

$$\begin{array}{}\text{(12)}& v_{\mathrm{opt}}\propto {\displaystyle \frac{c_{\mathrm{S}}ghT}{\mathrm{24}}}{\displaystyle \frac{\left|Q\right|S}{AD}}.\end{array}$$

Assuming that the steady state over a tidal cycle is driven mainly by the accelerating moment, especially in the upstream part where tidal mixing is relatively small and this gravitational circulation (D_g) is proportional to the dispersive residual velocity (D_g∝v_optE),

$$\begin{array}{}\text{(13)}& D_{\mathrm{g}}\propto {\left({\displaystyle \frac{c_{\mathrm{S}}g}{\mathrm{24}}}{\displaystyle \frac{S\left|Q\right|ET}{B}}\right)}^{\mathrm{1}/\mathrm{2}}.\end{array}$$

This equation indicates that the dimensionless dispersion coefficient is proportional to the root of the estuarine Richardson number N_R:

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{D_{\mathrm{g}}}{\mathit{\upsilon }E}}\propto {N_{\mathrm{R}}}^{\mathrm{0.5}}={\left(c_{\mathrm{S}}S{\displaystyle \frac{gh}{{\mathit{\upsilon }}^{\mathrm{2}}}}{\displaystyle \frac{\left|Q\right|T}{AE}}\right)}^{\mathrm{0.5}},\end{array}$$

where υ is the tidal velocity amplitude (m s⁻¹). The Richardson number describes the balance between the potential energy of the fresh water flowing into the estuary ($\mathrm{\Delta }\mathit{\rho }gh\left|Q\right|T/\mathrm{2}$) and the kinetic energy of the tidal flood flow (ρυ²AE∕2) (Fischer et al., 1979; Savenije, 2005; Zhang and Savenije, 2017).

2.2 Analytical solution for the dispersion equation

Equations derived from the maximum power concept are obtained along the estuary axis, and hence they can be used at any segment along the estuary. Then, Eq. (13) becomes

$$\begin{array}{}\text{(15)}& D_{\mathrm{g}}\left(x\right)=C_{\mathrm{3}}{\left({\displaystyle \frac{S\left|Q\right|ET}{B}}\right)}^{\mathrm{1}/\mathrm{2}},\end{array}$$

where C₃ is a factor (psu⁻¹ m s⁻¹) and all local variables are a function of x.

The following equations are used for the tidal excursion and width in alluvial estuaries:

$$\begin{array}{}\text{(16)}& {\displaystyle }E\left(x\right)=E_{\mathrm{0}}{e}^{{\mathit{\delta }}_H(x-x_{\mathrm{0}})},\text{(17)}& {\displaystyle }B\left(x\right)=B_{\mathrm{0}}{e}^{-(x-x_{\mathrm{0}})/b},\end{array}$$

where δ_H is the tidal damping rate (m⁻¹) and b is the geometric convergence length of the width in meters. A smaller b value implies stronger convergence (a stronger funnel shape). The subscript “0” represents parameters at the geometric boundary condition (x=x₀).

At the boundary, Eq. (15) is given by

$$\begin{array}{}\text{(18)}& D_{\mathrm{g}\mathrm{0}}=C_{\mathrm{3}}{\left({\displaystyle \frac{S_{\mathrm{0}}\left|Q\right|E_{\mathrm{0}}T}{B_{\mathrm{0}}}}\right)}^{\mathrm{1}/\mathrm{2}}.\end{array}$$

Substitution of Eqs. (16)–(18) into Eq. (15) gives

$$\begin{array}{}\text{(19)}& D_{\mathrm{g}}\left(x\right)=D_{\mathrm{g}\mathrm{0}}{\left({\displaystyle \frac{S}{S_{\mathrm{0}}}}\right)}^{\mathrm{1}/\mathrm{2}}{e}^{\mathrm{\Omega }(x-x_{\mathrm{0}})},\end{array}$$

with $\mathrm{\Omega }={\mathit{\delta }}_H/\mathrm{2}+\mathrm{1}/\left(\mathrm{2}b\right)$.

Differentiating D_g with respect to x and using the steady-state Eq. (11) results in

$$\begin{array}{}\text{(20)}& {\displaystyle \frac{\mathrm{d}{D}_{\mathrm{g}}}{\mathrm{d}x}}={\displaystyle \frac{D_{\mathrm{g}}}{\mathrm{2}S}}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}x}}+\mathrm{\Omega }{D}_{\mathrm{g}}=\mathrm{\Omega }{D}_{\mathrm{g}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\left|Q\right|}{A}}.\end{array}$$

The cross-sectional area A is given by

$$\begin{array}{}\text{(21)}& A\left(x\right)=A_{\mathrm{0}}{e}^{-(x-x_{\mathrm{0}})/a},\end{array}$$

where a is the convergence length of the cross-sectional area in meters.

Substituting Eq. (21) into Eq. (20) and in analogy with Kuijper and Van Rijn (2011) and Zhang and Savenije (2017), the solution of the linear differential Eq. (20) is

$$\begin{array}{}\text{(22)}& {\displaystyle \frac{D_{\mathrm{g}}}{D_{\mathrm{g}\mathrm{0}}}}=e^{\mathrm{\Omega }(x-x_{\mathrm{0}})}-{\displaystyle \frac{\left|Q\right|\mathit{\zeta }}{\mathrm{2}{A}_{\mathrm{0}}{D}_{\mathrm{g}\mathrm{0}}}}\left[e^{(x-x_{\mathrm{0}})/a}-e^{\mathrm{\Omega }(x-x_{\mathrm{0}})}\right],\end{array}$$

with $\mathit{\zeta }=a/(\mathrm{1}-\mathrm{\Omega }a)$.

At the salinity intrusion limit (x=L), D_g=0, resulting in

$$\begin{array}{}\text{(23)}& L=\mathit{\zeta }\ln \left(\mathrm{1}+{\displaystyle \frac{\mathrm{2}{A}_{\mathrm{0}}{D}_{\mathrm{g}\mathrm{0}}}{\left|Q\right|\mathit{\zeta }}}\right)+x_{\mathrm{0}}.\end{array}$$

The solution for the longitudinal salinity distribution yields

$$\begin{array}{}\text{(24)}& {\displaystyle \frac{S}{S_{\mathrm{0}}}}={\left\{\mathrm{1}-{\displaystyle \frac{\left|Q\right|\mathit{\zeta }}{\mathrm{2}{A}_{\mathrm{0}}{D}_{\mathrm{g}\mathrm{0}}}}\left[e^{(x-x_{\mathrm{0}})/\mathit{\zeta }}-\mathrm{1}\right]\right\}}^{\mathrm{2}}.\end{array}$$

This solution is comparable to other research. It is similar to Savenije (2005) if Ω=0, although his solutions had an empirical Van der Burgh coefficient K. In addition, the solution is the same as Kuijper and Van Rijn (2011) if a equals b, which implies that the depth is constant along the estuary.

With these new analytical equations, the dispersion and salinity distribution can be obtained using the boundary conditions (D₀ and S₀).

3 Empirical validation and discussion

The boundary condition is often taken at the geometric inflection point (x=x₀) if the estuary has one. The compilation of the Muar estuary in Fig. 3 is an example. Vertical dashed lines display the inflection point. If there is no inflection point such as in the Landak estuary, the boundary condition is taken at the estuary mouth (x₀=0). Figure 3 demonstrates that the geometry of the alluvial estuaries fits well on a semilogarithmic plot, supporting the exponential functions of the cross section and the width (Eqs. 17 and 21).

Figure 3 Semilogarithmic presentation of estuary geometry, comparing simulations (lines) to the observations (symbols), including cross-sectional area (blue diamonds), width (red dots), and depth (green triangles). Vertical lines display the inflection point.

Subsequently, Eq. (24) is used by confronting the solution with observations using appropriate boundary conditions. Appendix B shows how the new equation based on the maximum power concept works in 23 estuaries around the world. The Van der Burgh (VDB) method (Savenije, 2005), which has been proven to perform well in alluvial estuaries in different parts of the world and includes all mixing mechanisms, is used for comparison. Density-driven gravitational circulation is one part of the dispersive actions in estuaries. Hence, the total dispersive process from the Van der Burgh method (D_VDB) must be larger than the gravitational dispersion from the maximum power method (D_MP). The general geometry and measurements follow the database from Savenije (2012), Gisen (2015), and Zhang and Savenije (2017). The information on the VDB and MP methods is summarized in Table 1. Often there is more than one salinity observation in a certain estuary (labeled by alphabet), and the observation chosen from each estuary with a star-marked label is represented in Appendix B.

Table 1 Summary of application results using two methods.

^∗ Note: the observation chosen from each estuary is represented in Appendix B.

It can be seen that the simulated curves by the new MP method do not perform well in the wider part of the estuary (particularly upstream from the inflection point) where tidal mixing is dominant. However, the salinity observations can be very well simulated landward from the inflection point in most estuaries. In the Bernam, the Pangani, the Rembau Linggi, and the Incomati estuaries, the central part, where D_MP closely approach D_VDB, is well represented. In these estuaries, the calibration is slightly lower than the observations near the intrusion limit. In general, the dispersion derived with the maximum power method declines upstream from the inflection point in agreement with the total dispersion of the empirical Van der Burgh method, which corresponds to the theory that gravitational circulation is the dominant mixing mechanism in the landward part of these estuaries, especially in the center regime (e.g., Hansen and Rattray, 1965).

However, in the Thames (no. 8), the Delaware (no. 20), the Scheldt (no. 21), and the Pungwe (no. 22), the new approach seems not to work for both the salinity and dispersion profiles. In these estuaries tide-driven mixing is dominant. Figure 4 shows the relation between the geometry and the Van der Burgh coefficient K values. It can be seen that estuaries with poor performances by the MP approach have lower b∕B₀ and K values. However, not all estuaries with a strongly convergent geometry perform poorly, revealing that the geometry is not the only effect. According to the expression of Ω, tidal damping can play a role. In wide estuaries with strong convergence, the role of gravitational circulation is insufficient to describe the mixing. Tidal mixing processes such as lateral circulation, tidal pumping, and tidal shear are dominant. The Scheldt, with preferential ebb and flood channels, is a case in point (Nguyen et al., 2008). In addition, the Corantijn (no. 9) is considered uncertain because it has a low b∕B₀ value and contains few observations.

Figure 4 Comparison of the geometry to the Van der Burgh coefficient. Numbers show the labels of the estuaries.

Overall, the maximum power approach in open systems is a useful tool to understand the mixing processes in most estuaries. In the upstream part where the effect of the tide is small, gravitational circulation plays the main role. There, the MP approach yields good results. At the same time, the predictions upstream are more relevant for water users. Where the salinity is high, it is less relevant since the water is already too saline for domestic or agricultural use.

This study provides an approach to define the dispersion coefficient due to gravitational circulation, which is proportional to the product of the dispersive velocity of the gravitational circulation and the tidal excursion length (which is the longitudinal mixing length over which one particle travels during a tidal cycle). The dispersive velocity actually represents the strength of the density-driven mechanism. Based on the maximum power method (Eq. 15), the dispersive velocity can be described as

$$\begin{array}{}\text{(25)}& v\propto {\left({\displaystyle \frac{S\left|Q\right|T}{BE}}\right)}^{\mathrm{1}/\mathrm{2}}.\end{array}$$

Hence, the dispersive flow due to gravitational circulation strengthens with larger freshwater discharge $\left|Q\right|$ (more stratification) and weakens with stronger tide E (less stratification).

Using the calibrated dispersion coefficient at the inflection point, C₃ can be calculated. Except in estuaries with poor performance, C₃ values range from $\mathrm{3.5}\times {\mathrm{10}}^{-\mathrm{3}}$ to $\mathrm{10.0}\times {\mathrm{10}}^{-\mathrm{3}}$ with an average of $\mathrm{6.8}\times {\mathrm{10}}^{-\mathrm{3}}$ (the relative standard deviation equals 0.26). Using the average C₃ value to predict D_g0 (Eq. 18), Fig. 5 shows how the predictive equation performs. It reveals that almost all the data fall within a factor of 2, and the maximum power method underestimates the dispersion coefficient in estuaries with low b∕B₀ values (in red) in which gravitational circulation is not enough to describe the total dispersive processes. Finally, the R² value of the regression in Fig. 5 equals 0.86. Considering all the uncertainties in the measurement, C₃ equalling $\mathrm{6.8}\times {\mathrm{10}}^{-\mathrm{3}}$ is a promising approximation to predict D_g0.

Figure 5 Comparison of calibrated and predicted D_g0 values by using $C_{\mathrm{3}}=\mathrm{6.8}\times {\mathrm{10}}^{-\mathrm{3}}$. Labels in red, indicating that the estuaries have relatively poor performance, are presented for validation.

Finally, there is uncertainty about the timescale of reaching this optimum. If this timescale is longer than the tidal period, then such an optimum is not reached. In a low-flow situation, however, which is the critical circumstance for salt intrusion, the variation of the river discharge is slow (following an exponential recession). If the timescale of flow recession is large compared to the timescale of salinity intrusion then it is reasonable to assume that the maximum power optimum is approached.

4 Combination of the MP and VDB methods

The fact that the MP method works well for density-driven mixing but not for tide-driven mixing, whereas the VDB method works well for the combination of the two, offers an excellent opportunity for the combination of the two methods. The VDB method requires two parameters: the K of Van der Burgh and the dispersion coefficient at the downstream boundary D₀, while the MP method only requires the downstream boundary condition D_g0. The dispersion of the VDB method, which deals with all mixing processes, should therefore always be larger than the dispersion determined by the MP method. Hence, the MP method can be used to impose an additional constraint on the calibration of the VDB method, which reduces the potential equifinality between K and D₀. Appendix B shows the result of this mixed calibration approach: the dispersion of the VDB method is always higher than the dispersion of the MP method, and the resulting fit by the VDB method is quite acceptable.

This combined approach also allowed for more accurate predictive equations as derived before. The correlation between K and the estuary geometry is strong, as shown in Fig. 4. This relation can be used as a predictive equation for K. Also, the predictive equation for D_g0 is powerful, as can be seen in Fig. 5, except for very wide estuaries where calibration remains necessary and where this predictive equation can be used as a 1st-order estimate for D₀.

5 Conclusions

An estuary is an open system that has a maximum power limit when the accelerating source is stable. This study has described a moment balance approach to nonthermal systems, yielding a new Eq. (15) for the dispersion coefficient due to the density-driven gravitational circulation. It shows that the dispersive action is closely related to the salinity, the freshwater discharge, the tide, and the estuarine width. This equation has been used to solve the tidal average salinity and dispersion distributions together with the steady-state Eq. (11). The maximum power model has then been validated with 50 salinity observations in 23 estuaries worldwide and compared with the Van der Burgh method. Generally, the new equation is a helpful tool to analyze the salinity distribution in alluvial estuaries, providing an alternative solution for the empirical Van der Burgh method in estuaries where gravitational circulation is the dominant mixing mechanism. A predictive equation for dispersion at the geometric boundary has also been provided.

As can be seen in Appendix B, the gravitational dispersion is always smaller than the total effective dispersion obtained by the Van der Burgh method. In all estuaries that have a wide mouth, we see substantial tide-driven dispersion, most probably as a result of interacting preferential ebb and flood channels. This tide-driven mechanism is responsible for the (sometimes pronounced) concave slope of the salinity curve near the mouth. In the middle reach where the salinity gradient is steepest, density-driven dispersion is dominant and equals the total effective dispersion. Further upstream, where the salinity gradient gradually tends to zero and the estuary becomes narrower, we see the tide-driven circulation again becoming more prominent. This is in the part of the estuary where the width-to-depth ratio becomes smaller and the bank shear results in more pronounced lateral velocity gradients and hence more pronounced lateral circulation. The tide-driven mixing mechanism is particularly strong in macro-tidal estuaries such as the Thames, the Scheldt, the Pungwe, and the Delaware.

This study is a further development of the paper by Zhang and Savenije (2018), which also considered gravitational circulation based on the maximum power concept but which did not consider the associated frictional dissipation. The approach followed in this paper maximizes the work performed by the driving gravitational torque to mix the fresh and saline water, taking account of the energy dissipation associated with this mixing. As a result, we found a solution that combines well with the empirical Van der Burgh method, providing an additional constraint for its calibration. Because the total mixing of the Van der Burgh method (D_VDB) should be larger than the gravitational mixing of the maximum power concept (D_MP), the calibration of the Van der Burgh method is more constrained. As a result, the Van der Burgh K and the dispersion at the boundary D₀ can be correlated with physically observable parameters through analytical equations, which makes the Van der Burgh method a more powerful predictive model that can be applied to alluvial estuaries worldwide. More reliable observations in other estuaries are suggested to validate these maximum power and Van der Burgh methods.

Appendix A: Notation

Table A1 Notations for symbols used in this study.

Appendix B: Application of the maximum power method

This Appendix represents the application in 23 estuaries around the world of the maximum power method for determining the dispersion coefficient and the salinity distribution using Eqs. (22) and (24), compared to salinity observations. The empirical Van der Burgh method is included as a reference.

Figure B1 Left: Application of the analytical solution from the maximum power method (solid lines) to observations (symbols) for high water slack (HWS, in red) and low water slack (LWS, in blue). The green line shows the tidal average (TA) condition. Dash dot lines reflect applications of the Van der Burgh method. Vertical dash lines display the inflection point. Right: Simulated dispersion coefficient using different methods.