## Abstract

A universal law of estuarine mixing is derived here, combining the approaches of salinity coordinates, Knudsen relations, total exchange flow, mixing definition as salinity variance loss, and the mixing–exchange flow relation. As a result, the long-term average mixing within an estuarine volume bounded by the isohaline of salinity *S* amounts to *M*(*S*) = *S*^{2}*Q*_{r}, where *Q*_{r} is the average river runoff into the estuary. Consequently, the mixing per salinity class is *m*(*S*) = ∂_{S}*M*(*S*) = 2*SQ*_{r}, which can also be expressed as the product of the isohaline volume and the mixing averaged over the isohaline. The major differences between the new mixing law and the recently developed mixing relation based on the Knudsen relations are threefold: (i) it does not depend on internal dynamics of the estuary determining inflow and outflow salinities (universality), (ii) it is exactly derived from conservation laws (accuracy), and (iii) it calculates mixing per salinity class (locality). The universal mixing law is demonstrated by means of analytical stationary and one-dimensional and two-dimensional numerical test cases. Some possible consequences for the salinity distribution in real estuaries are briefly discussed. Since the mixing per salinity class only depends on the river runoff and the chosen salinity, and not on local processes at the isohaline, low-mixing estuaries must have large isohaline volumes and vice versa.

## 1. Introduction

In classical estuaries, density gradients are dominated by salinity gradients, such that salinity coordinates offer promising perspectives to analyze estuarine physics and biogeochemistry. For this reason, Walin (1977) proposed salinity as a coordinate to analyze the overturning circulation and internal mixing in the Baltic Sea and other estuarine systems. Basic quantities defined by Walin (1977) are scalars (salinity, nutrients) integrated over all salinities smaller than the salinity of a certain isohaline *S* as well as their derivatives with respect to salinity. By doing so, he obtained tracer concentrations per salinity class. In the same way, he defined advective and diffusive fluxes of volume and salt through isohaline surfaces. The present study exploits the framework proposed by Walin (1977) to quantify estuarine mixing with respect to isohaline surfaces. Based on an isohaline perspective, MacCready and Geyer (2001) investigated the role of turbulent diahaline salt fluxes for maintaining estuarine salt intrusions. The same concept of isohaline volumes [infinitesimal volume *δV* included in an infinitesimal salinity interval *δS*, see Eq. (25) of the present manuscript] and their change due to advective and diffusive diahaline transports has been used by MacCready et al. (2002) to explain the long-term isohaline structure of tidal estuaries. Hetland (2005) applied the isohaline coordinate framework to investigate mixing in river plumes independently of the highly variable position of the plume. The new theory presented here is partially based on the isohaline theory by MacCready et al. (2002) with the difference that in the present study the salinity square budget (Burchard and Rennau 2008) is considered additionally.

Exchange flows in estuaries are typically quantified by means of the Knudsen (1900) relations, also developed for the Baltic Sea, which calculate inflow and outflow volumes and salinities for fixed transects across the estuary and relate these quantities to each other by means of conservations laws for mass and salt (see the recent review by Burchard et al. 2018). To refine the resolution of the exchange flow at a fixed cross-estuarine transect, MacCready (2011) developed the total exchange flow (TEF) analysis framework in which inflow and outflow of volume and salt across the transect are resolved in salinity coordinates, consistent with the Knudsen (1900) relations. This analysis framework proved very useful, as demonstrated by several exchange flow studies based on it (Sutherland et al. 2011; Ganju et al. 2012; Chen et al. 2012; Lemagie and Lerczak 2015; Gräwe et al. 2015b, 2016; Purkiani et al. 2016; Wang et al. 2017; Burchard et al. 2018).

To quantify salinity mixing in the ocean, the destruction of salinity microstructure variance (Nash and Moum 2002) or the decay of ensemble averaged salinity variance (Burchard and Rennau 2008) have proven to be suitable measures. In average, both measures are the same (Umlauf and Burchard 2005). The decay of ensemble averaged salinity variance has been used in several studies to quantify physical (Burchard et al. 2009; Ralston et al. 2010; Wang et al. 2017; Li et al. 2018; Wang and Geyer 2018) and numerical (Burchard and Rennau 2008; Klingbeil et al. 2014; Ralston et al. 2017) mixing in estuaries.

Recently, the relation between estuarine exchange flow based on the Knudsen (1900) relations (MacCready et al. 2018) and based on the TEF analysis framework (Burchard et al. 2019) has been formulated for a fixed transect across the estuary. For long-term averages, the total mixing within an estuary can be estimated as the product of inflow and outflow salinity and the river runoff. The physical mixing analysis proposed by Burchard et al. (2009) and the numerical mixing analysis proposed by Burchard and Rennau (2008) and improved by Klingbeil et al. (2014) give mixing values at every time step and grid box. In addition, bulk mixing estimates such as TEF and the present isohaline framework reduce complexity through spatial integration and time averaging of mixing to gain a deeper understanding of principles of estuarine physics. Results from the fully resolved numerical model are always available for a differential analysis in terms of mixing locations, isohaline positions and underlying dynamic processes.

However, the determination of the inflow and outflow salinities requires careful analysis of accurate numerical model results, and the method by MacCready et al. (2018) is not exact. Furthermore, it does not provide local mixing estimates in an easy way. Therefore, the objective of this paper is to present a universal law of estuarine mixing based on isohaline coordinates in the entire estuary rather than on isohaline coordinates on a fixed transect. An estuary in the sense of this law is defined such that the freshwater budget is fully represented by a river runoff *Q*_{r}, possibly separated into various tributaries, neglecting precipitation, evaporation and exchange with the groundwater. To derive this law, the approaches of salinity coordinates, Knudsen relations, total exchange flow, mixing definition as salinity variance loss, and the mixing–exchange flow relation are combined. The presentation focuses on the mathematical derivation of this theory and gives a few idealized examples demonstrating its validity. A few possible implications for estuarine physics are briefly discussed here as well, but a thorough analysis of the consequences for the interpretation of processes in various estuaries cannot be given at this stage.

Section 2 briefly motivates the universal law of estuarine mixing in salinity coordinates. A mathematical derivation based on volume and salt conservation is given in section 3. Examples from idealized model simulations are shown in section 4. Some possible consequences on estuarine dynamics are briefly discussed in section 5. Conclusions are drawn in section 6.

## 2. Derivation from the Knudsen mixing estimate

Here, the new universal law of estuarine mixing is briefly motivated. A mathematical derivation is provided in section 3. The volume-integrated and time-averaged salinity mixing within an estuary separated from the ocean by a transect *T* fixed in time is given as

see Burchard et al. (2009), where *χ*^{s} is salinity mixing per unit volume [see Nash and Moum (2002) and Eq. (9) below], *V*_{T} is the time-varying estuarine volume bounded by the transect *T* and angular brackets ⟨⋅⟩ denote time averaging. According to MacCready et al. (2018), the volume-integrated and time-averaged mixing in the estuary under tidally periodic or stationary conditions or for very long averaging periods can be estimated as

where *s*_{in} and *s*_{out} are the inflow and outflow salinities at *T* calculated by means of the Knudsen (1900) relations and *Q*_{r,T} is the average riverine freshwater runoff into the estuarine volume bounded by *T*. Note that the averaged volume transport *Q*_{r,T} through the transect *T* depends on the location of the transect, since runoff from more tributaries to the estuary need to be included when the transect is moved closer to the ocean.

Taking a certain isohaline *S* as the boundary of the estuary instead of the transect *T*, the mixing of the estuarine volume bounded by *S*, that is, the volume *V*(*S*) including all salinities *s* ≤ *S*, is defined here as

Under conditions of stationarity or periodicity or for long-term averaging, the averaged volume transport through the isohaline *S* equals the freshwater flux *Q*_{r} and this runoff does not depend on *S*, since all runoff is assumed to occur at zero salinity. Inflow and outflow salinities with respect to the isohaline *S* are *s*_{in} = *s*_{out} = *S*, such that the time-averaged mixing integrated over *V*(*S*) is approximated as

With Eq. (4), also the mixing per salinity class,

can be calculated. Here, *m*(*S*) is the infinitesimal mixing within an infinitesimal salinity interval *δS* centered around the isohaline *S*. Integration of *m*(*S*) over salinities ranging from 0 to *S* gives *M*(*S*). It will be shown in section 3 that the mixing estimates in Eqs. (4) and (5) are exact relations for stationary or periodic estuaries and converge to the exact integrated mixing for all estuaries when the averaging period is sufficiently long.

Equation (5) implies a fundamental law of estuarine physics: The long-term averaged mixing per salinity class in any estuary is twice the product of the salinity and the freshwater runoff.

## 3. Derivation from conservation laws

Here, a mathematical derivation of the universal mixing law [Eq. (5)] is provided. The basic physical equations from which this theory is derived are the continuity equation (volume conservation, assuming incompressibility of seawater),

and the salinity equation (salt conservation),

By multiplying Eq. (7) by 2*s* and rearranging, the salinity-squared equation is derived:

with the salt mixing per unit volume,

see Burchard et al. (2009). In Eqs. (6)–(8), (*u*, *υ*, *w*) are the velocity components in direction of the Cartesian coordinates (*x*, *y*, *z*), *s* is salinity, and *K*_{h} and *K*_{υ} are horizontal and vertical eddy diffusivities, respectively. Surface and bottom fluxes of volume and salinity such as precipitation, evaporation, ground-water discharge or influx of seawater into the sediment and consequently fluxes of squared salinity are neglected.

For any tracer *c*, we define the total (advective plus diffusive) tracer flux per square meter normal to a surface (which may for example be a vertical estuarine transect *T* or an isohaline surface *S*) as

where *u*_{n} is the outgoing velocity normal to the transect or the moving isohaline, *∂*_{n}*c* denotes the tracer gradient normal to the surface (pointing outwards), and *K*_{n} is the eddy diffusivity normal to the surface. In the case of a vertical transect *T*, *K*_{n} = *K*_{h} is obtained, and in the case of an isohaline surface *S*, *K*_{n} would be the diahaline diffusivity.

In Eq. (10), $Fac$ is the advective tracer flux per square meter and $Fdc$ is the diffusive tracer flux per square meter. Moreover, *F*^{1} = *F* = *u*_{n} is the outgoing volume flux normal to the surface. It follows from Eq. (10) that

The time-averaged boundary flux of *c* at a fixed transect *T* separating the estuary from the ocean can be expressed in salinity coordinates (MacCready 2011; Burchard et al. 2018, 2019):

with

where *A*_{T}(*S*) is the part of the open-boundary transect with salinities *s* larger than a specific salinity *S*, *s* ≥ *S*. Consequently, *A*_{T}(0) is the entire area of the transect. In Eq. (13), $QTc\u2061(S)$ is the cumulative tracer transport (negative for outgoing transport) at the transect *T* for all salinities larger than *S*, as defined by MacCready (2011). With that, $qTc\u2061(S)$ is the tracer transport per salinity class across the transect. The minus sign in the definition of $qTc\u2061(S)$ is adopted from MacCready (2011) as well to give positive values for inflow at high salinities.

Equivalently, the time-averaged mixing per salinity class *m*(*S*) can be defined:

where *V*(*S*) is the estuarine volume with salinities *s* smaller than a specific salinity *S*, *s* ≤ *S*. In Eq. (14), *M*(*S*) is the salinity mixing integrated over all salinities smaller than *S*. Properties integrated over volumes defined by isohalines and their derivatives with respect to salinity had already been proposed by Walin (1977). The only difference here is the sign convention in Eq. (14). While Walin (1977) proposed volumes defined by salinities larger than the salinity of the isohaline, here volumes for salinities smaller than the isohaline salinity are defined. This has the advantage that these volumes are located inside the estuary, such that *V*(*S*) represents a part of the estuarine volume.

Integrating Eqs. (6)–(8) over a volume *V*_{T}(*S*) [defined as the fraction of *V*(*S*) that is located inside *V*_{T}] results in

where *B*_{T}(*S*) is the part of the isohaline surface with salinity *s* = *S* that is located inside *V*_{T} (red isohaline in Fig. 1) and *A*_{r,T} is the area through which freshwater of zero salinity is discharged into the control volume. The area of the transect *T* with salinities *s* < *S* is *A*_{T}(0) − *A*_{T}(*S*) (see Fig. 1).

where

is the temporally averaged freshwater runoff into *V*_{T}(*S*). Note that for the case that the isohaline with salinity *S* is located fully outside *V*_{T} the following relations hold: *Q*_{T}(*S*) = 0 and *B*_{T}(*S*) = 0, such that Eq. (16) is identical to Eqs. (30)–(32) in Burchard et al. (2019).

Assuming long-time averaging or periodic conditions [such that the left hand sides of Eq. (16) converge to zero] and choosing an arbitrary isohaline *S* that is located fully inside *V*_{T} [such that *Q*(*S*) = *Q*(0), *Q*^{s}(*S*) = *Q*^{s}(0) and $Qs2\u2061(S)=Qs2\u2061(0)$, as indicated by the bold blue isohaline in Fig. 1], Eq. (16) simplifies to

where *B*(*S*) is the complete surface of the isohaline with salinity *S*. With Eqs. (10) and (11), the first two equations of Eq. (18) give

The meaning of Eq. (19) is that in a long-term average the diffusive transport through an isohaline is exactly balanced by an advective transport. The diffusive transport has the tendency to move the isohalines toward lower salinities (denoted as *drift* by MacCready et al. 2002) and the advective transport moves them back, such that there is an apparent advection through isohalines. The derivative of Eq. (19) with respect to *S* is identical to Eq. (4.2) by MacCready et al. (2002). With Eq. (19) and using Eqs. (10) and (11), the last equation in Eq. (18) results in

which proves the universal law in Eq. (4) and consequently also Eq. (5), that is, *m*(*S*) = 2*SQ*_{r}. Taking the divergence (in salinity coordinates) of the last equation in Eq. (18), it can also be argued that the mixing per salinity class equals the divergence of the integrated diahaline transport of squared salinity, that is,

The estuarine mixing law might be easier to interpret by expressing *m* as the product of the isohaline volume and the average mixing per unit volume and salinity class. To derive such a relation, discrete salinity classes are considered, represented by a salinity value *S* and a salinity range Δ*S*. Then, the volume-integrated mixing per discrete salinity class is

where *S*′ is the variable of integration. This mixing per discrete salinity class can also be directly calculated from Eq. (14) as

where *V*_{S,ΔS} and $V\xafS,\Delta S$ are the instantaneous and averaged volume of the discrete salinity class (denoted as *differential volume* by MacCready et al. 2002), respectively, and $\chi s\xafS,\Delta S$ is the average mixing per unit volume of the respective salinity class. Equating Eqs. (22) and (23) gives

For Δ*S* → 0, a continuous version of Eq. (24) is obtained:

where *υ*(*S*) = ∂_{S}*V*(*S*) is the isohaline volume, as defined by Walin (1977) and $\chi s\xaf\u2061(S)$ is the average mixing per salinity class. Thus, integration of isohaline volumes over the salinity range from 0 to *S* results in the estuarine volume *V*(*S*) bounded by *S*.

## 4. Idealized examples

### a. One-dimensional stationary estuary

Here an analytical solution for estuarine mixing in salinity coordinates is presented. Assume a stationary vertically mixed estuary with constant depth *H* and width *W*, constant flow velocity *u* < 0 and constant longitudinal eddy diffusivity *K*_{h}. With this, the freshwater runoff is *Q*_{r} = −*uWH*. The length of the estuary is *L*, and salinity *s*(0) = *s*_{0} at the mouth and *s*(*L*) = 0 at the river end. Under these conditions, the salinity budget in Eq. (7) reduces to

The analytical solution fulfilling Eq. (26) and the boundary conditions is

with

It should be noted that at *x* = *L* there is a salt flux into the estuary of

which vanishes for −*uL*/*K*_{h} → ∞ since then $s^0\u21920$ and ∂_{x}*s* → 0 hold at *x* = *L*. Therefore, to obtain results representative for estuaries without salinity in the freshwater inflow, scenarios with small values of −*uL*/*K*_{h} have to be considered. The volume-integrated mixing for an estuarine volume bounded by *x* is

To transform Eq. (30) into a salinity coordinate *S*,

and consequently

For the limit of a relatively long estuary with a large value of −*uL*/*K*_{h}, the mixing rates from Eqs. (4) and (5) result.

Example solutions for two different eddy diffusivities are shown in Fig. 2. The estuary is chosen to be very long (*L* = 1000 km) in order to minimize the salt inflow at *x* = *L*. For the higher eddy diffusivity (*K*_{h} = 1000 m^{2} s^{−1}), salinity diffuses much farther into the estuary than for the small diffusivity (*K*_{h} = 500 m^{2} s^{−1}), see Fig. 2a, resulting in mixing occurring farther inside the estuary (Fig. 2b). However, the mixing curves in salinity coordinates do collapse, apart from a small difference (not visible in the plots) resulting from the different salt transports into the estuary at *x* = *L* due to the second term in Eq. (33).

### b. One-dimensional periodic estuary

Using the TEF analysis framework developed by MacCready (2011), it has been shown by Burchard et al. (2019) that a simple vertically integrated tidal estuary can be characterized by an estuarine circulation in salinity coordinates at a fixed transect *T*. The numerically simulated dynamics of this one-dimensional tidally periodic estuary was in one case based on purely numerical mixing, and in another case on a combination of physical and numerical mixing. In both cases, the numerically induced mixing was intentionally dominant due to the application of a first-order upstream advection scheme (Burchard and Rennau 2008).

This one-dimensional estuary has a length of *L* = 100 km, and the bottom coordinate decreases linearly from a depth of 15 m at the mouth (*x* = 0) to 5 m at the weir (*x* = *L*), where a constant freshwater runoff of *Q*_{r} = 200 m^{3} s^{−1} is prescribed. At the mouth, the simulations are forced by a constant salinity of *s*_{0} = 30 g kg^{−1} and a semidiurnal tide with an elevation amplitude of 2 m. The numerical discretization is carried out by 100 spatial increments along the estuary and 1000 time steps per tidal period. Numerical and physical mixing are analyzed exactly using the method proposed by Klingbeil et al. (2014), based on the interfacial numerical advective and diffusive fluxes. The model simulations are first initialized with zero velocity and a constant salinity of 30 g kg^{−1} and once a periodic state is obtained, one tidal period is analyzed. For further details regarding the simulations, see Burchard et al. (2019).

In the simulation including physical mixing, the isohalines migrate about 10–20 km farther into the estuary than for the simulation based on numerical mixing only (Fig. 3a). This has the consequence that integrated mixing *M*(*x*) is significantly elevated when physical mixing is included (Fig. 3b). Nonetheless, the curves for the integrated mixing with respect to isohalines almost collapse and are very close to the theoretical value of *M*(*S*) = *S*^{2}*Q*_{r} (Fig. 3c). At salinities for which isohalines leave the model domain, that is, when these isohalines intersect with the open boundary, the integrated mixing is reduced. Since the simulations are periodic, all deviations between the curves for the numerical experiments and the theoretical result for salinities lower than these maximum salinities are due to numerical issues. Figure 3d shows the expected result that the mixing per salinity class closely follows the theoretical result of *m*(*S*) = 2*SQ*_{r}. The experimental curves are noisy, since in this one-dimensional scenario the number of mixing values per salinity class is relatively low. In contrast to this, the two-dimensional scenario discussed below leads to many more mixing values per salinity class such that the curve for *m*(*S*) is much smoother (Fig. 4d). This behavior is similar to the numerical problems in evaluating Eq. (13), see the discussion by MacCready et al. (2018) and the numerical analysis by Lorenz et al. (2019).

### c. Two-dimensional almost periodic estuary

To demonstrate the validity of the universal mixing law also for results from an ocean model, the General Estuarine Transport Model (GETM; see Burchard and Bolding 2002; Gräwe et al. 2015a) is used to simulate a two-dimensional vertically resolved tidal channel without lateral variation (Burchard et al. 2019; Klingbeil et al. 2019). The vertical turbulent fluxes are calculated by means of a two-equation turbulence closure model (Umlauf and Burchard 2005) while horizontal turbulent fluxes are ignored. As advection scheme for velocities and salinity, the monotone Superbee scheme (Pietrzak 1998) is used, which is known for its antidiffusive properties. Physical and numerical mixing are calculated by means of the analysis method proposed by (Klingbeil et al. 2014).

The estuary is *L* = 100 km long and *W* = 500 m wide and the bottom depth is linearly decreasing from 15 m at the mouth (*x* = 0) to 5 m at the end of the estuary (*x* = *L*), where a weir is assumed to limit the tidal excursion and a freshwater runoff of *Q*_{r} = 50 m^{3} s^{−1} is prescribed. At the mouth of the estuary, a constant salinity of 30 g kg^{−1} is prescribed, and the model is forced by means of a semidiurnal periodic surface elevation oscillation with an amplitude of 0.6 m. The simulations are run for 110 tidal periods of which the last 10 periods are analyzed. At that stage, the model is almost periodic apart from an internal wave of a period shorter than the tidal period. All further details of the model simulation are given by Burchard et al. (2019).

The results show a monotonically decreasing depth-mean and tidally averaged salinity, which is transported into the first 50 km of the estuary (Fig. 4a). Due to the tidal oscillations and the vertical stratification, maximum and minimum salinities are strongly deviating from the mean. The integrated mixing computed using the method by Klingbeil et al. (2014) is almost exclusively consisting of physical mixing, with numerical mixing being very small or even negative in the bulk of the estuary (Fig. 4b). The integrated mixing *M*(*S*) again follows closely the theoretical curve and only deviates from it outside the highest salinity isohaline not leaving the estuary (Fig. 4c). This is more clearly visible for the mixing per salinity class *m*(*S*), see Fig. 4d, which agrees well with the theoretical linear curve inside this critical salinity, apart from numerical inaccuracies discussed above.

## 5. Discussion

The mathematical analysis of section 3 suggests that the isohaline structure of any freshwater dominated estuary should permanently have the tendency to adjust to a state where the mixing per salinity class is twice the product of salinity and freshwater runoff. Here, a few possible consequences for estuaries are briefly discussed.

### a. Confining estuaries by isohalines

If the extent of an estuary is defined by means of a bounding isohaline (as for example given by a salinity of 20 or 30 g kg^{−1}) instead of a fixed transect (as for example near the Battery for the Hudson River Estuary or at Darss Sill and Drogden Sill for the Baltic Sea), then the long-term averaged estuarine mixing exclusively depends on the average river runoff. According to Eq. (4), the mixing equals the square of the salinity times the river runoff. In the case of a fixed transect, the exchange flow needs to be additionally considered (MacCready et al. 2018), see Eq. (2). Based on a river runoff into the Baltic Sea of *Q*_{r} = 16 098 m^{3} s^{−1} and an exchange flow defined by inflow and outflow salinities of *s*_{in} = 16.59 g kg^{−1} and *s*_{out} = 8.67 g kg^{−1}, Burchard et al. (2019) estimated a total mixing of *M*_{T} = 2.32 10^{6} m^{3} s^{−1}(g kg^{−1})^{2}, averaged over 66 years, using Eq. (2). The same value would be obtained for mixing in the estuarine volume bounded by the isohaline of *S* = (*s*_{in}*s*_{out})^{1/2} = 12 g kg^{−1}, the position of which however is highly variable in time and space. In the fjord-type Baltic Sea, the 12 g kg^{−1} isohaline is strongly stretched and might extend over the entire western and central Baltic Sea with a length along a transect through the Baltic Sea of more than 800 km, as already observed by Wüst et al. (1957). Also for tidal estuaries, mixing is low during neap tides, such that isohalines in the brackish water range are strongly distorted (Warner et al. 2005). This means that for fjord-type systems isohalines in the brackish water range are not suitable as outer boundaries. Instead, near-ocean salinity isohalines should be used to define the limit of an estuary. On the other hand, isohalines with too high salinities could extend far into the ocean and might overlap with neighboring estuaries. Although the extent of an estuary could be defined by a bounding isohaline instead of a fixed transect, but identifying the bounding isohaline optimally characterizing the estuary might be tricky.

Across an isohaline an exchange flow in the sense of the TEF analysis framework by MacCready (2011) does not exist. TEF-based exchange flows are characterized by their property that inflows and outflows at the same salinity compensate each other, and since along an isohaline surface all salinities are identical, the TEF-based exchange flow vanishes. As mentioned above, in long-term average advective and diffusive fluxes normal to any isohaline are exactly balancing each other in each point. While in average the diffusive fluxes are directed toward lower salinities (up-estuary) and the advective fluxes are directed toward higher salinities (down-estuary), there will be distinct subareas of the isohaline surface with up-estuarine and others with down-estuarine advective fluxes. This expresses the estuarine exchange flow with respect to isohaline surfaces.

### b. Variation of isohaline volumes

The formulation (25) shows that the mixing per salinity class can be expressed as the product of the isohaline volume *υ*(*S*) and the average mixing per volume and salinity class, $\chi s\xaf\u2061(S)$, which equals twice the product of salinity and freshwater runoff. This leads to the interpretation that closer to the mouth of the estuary (high salinity) the average mixing per unit volume should be much higher than in the brackish water range (low salinity) or the isohaline volumes are larger closer to the ocean. The fact that mixing per salinity class is higher for high salinities might appear counterintuitive at first sight. It is however obvious that the long-term averaged advective salinity transport through an isohaline, *SQ*_{r}, is proportional to salinity and directed down the estuary. On long time scales, this advective transport is exactly compensated by an up-estuarine diffusive salt flux, see Eq. (19) of this manuscript and Eq. (2.4) by MacCready et al. (2002). In turn, increased diffusive fluxes result in increased mixing. As already stated by Hetland (2005) for river plumes, *strong mixing requires only a small isohaline area, whereas weaker mixing requires a larger isohaline area to maintain the same total freshwater flux across the isohaline.*

An example of this is given in Fig. 5 for the two-dimensional almost periodic numerical estuary presented in section 4c. Interestingly, in this case the isohaline volume *υ*(*S*) does not vary a lot except for very low salinities (Fig. 5a), but $\chi s\xaf\u2061(S)$ varies almost linearly with salinity (Fig. 5b). These dependencies might be entirely different in real estuaries, but the product of *υ*(*S*) and $\chi s\xaf\u2061(S)$ must be equal to *m*(*S*) = 2*SQ*_{r} for long-term averages. Specifically for estuaries with strongly variable cross-sectional area, a stronger spatial variation in average mixing per isohaline volume and thus higher variation in isohaline volumes would be expected.

The high diversity of estuarine systems had been ordered by Geyer and MacCready (2014) by placing a variety of estuaries into a two-dimensional parameter space, spanned by nondimensional parameters for mixing and for freshwater runoff. Estuaries range from fjords with low mixing and low runoff to partially mixed estuaries and to time-dependent salt wedge estuaries at high mixing and high runoff. Although the spread of the estuaries in their diagram is wide, some positive correlation between mixing and runoff can be seen. When selecting a fixed bounding isohaline for all estuarine systems (e.g., 30 g kg^{−1}), then the universal mixing law (4) establishes a positive linear relationship between runoff and mixing valid for all estuaries. This mixing law does not distinguish between the various types of estuaries classified by Geyer and MacCready (2014). However, when considering mixing per salinity class as product of isohaline volume and average mixing per salinity class, some interesting insights can be gained.

When the mixing per unit volume in an estuary is reduced, such as during the transition from spring to neap tide, the isohaline volumes should increase, that is, isohalines should become more widely spaced or their surfaces should increase, and vice versa during the transition from neap to spring tides. This can be seen for the Hudson River Estuary in the numerical model study by Warner et al. (2005) and the observations by Lerczak et al. (2006). During spring tide, the isohalines are far more upright than during neap tide, where they intersect with the surface tens of kilometers downstream of their bottom location. This means that the isohaline volumes are substantially increased during neap tides, to compensate for the lower mixing. In low-mixing estuaries such as fjords, isohalines cover substantial areas in order to provide sufficient isohaline volumes, as shown in the case of the 12 g kg^{−1} isohaline of the Baltic Sea (Wüst et al. 1957), see the discussion above.

In systems with side channels and harbor basins such as San Francisco Bay (Fischer 1976) or Rotterdam Waterway (de Nijs et al. 2011) isohaline volumes might get entrapped sideways at certain phases of the tide, a process strongly contributing to longitudinal dispersion. In case that this dispersion is related to enhanced local mixing, it will contribute to the decrease of the respective isohaline volumes, no matter, if they are divided into subvolumes or not.

Increases in freshwater runoff should lead to increased mixing at all salinity classes, by increasing the isohaline volumes or by increasing average mixing per isohaline volume. There are various mechanisms for adjusting the mixing to increased river runoff. Increased isohaline volumes could be obtained by moving isohalines farther down the estuary into locations with wider cross section or by pushing even low-salinity isohalines outside the estuary into the river plume (Ralston et al. 2010) or by increasing the isohaline spacing. Another adjustment mechanism would be to increase mixing per volume and salinity class, $\chi s\xaf\u2061(S)$, for example, by displacing isohaline volume to locations with smaller cross-sectional area or by increased bottom friction due to higher residual runoff velocity.

Estuarine adjustment time scales depend on a variety of processes (MacCready 1999, 2007). With the isohaline structure of all estuaries permanently trying to converge toward a state fulfilling the universal mixing law, they may be distinct in the time scales at which they are lagging behind this law, because of the variability in external forcing. Instantaneous states of estuaries with long time scales and high variability in forcing might always be far away from their mean state prescribed by the universal law.

The isohaline structure reproduced by numerical models is a result of the total mixing, that is, the sum of physical and numerical mixing. Since numerical mixing may be high in estuarine studies (Hofmeister et al. 2011; Ralston et al. 2017), the numerical mixing certainly needs to be included when testing the validity of the universal mixing law by means of numerical model simulations.

### c. The estuary–river plume continuum

The present isohaline theory implicitly treats estuaries, river plumes and the connected regions of freshwater influence (ROFIs; see Simpson 1997) as a dynamic continuum and not separately as in MacCready et al. (2002) (estuaries) and Hetland (2005) (river plumes). The higher the salinity of a specific isohaline under consideration, the more volume is included, such that the study region may be extended to outside the estuary. The present theory is developed for closed isohalines. This implies that for high salinity isohalines several neighboring estuaries of a ROFI system might have to be included into the consideration, including their combined river runoff. Still, in the long-term mean, these systems should follow the universal mixing law.

## 6. Conclusions

The universal law of estuarine mixing derived here differs from recently developed formulations in various aspects. The Knudsen-based hierarchy of volume-integrated mixing estimates for estuaries developed by MacCready et al. (2018) and Burchard et al. (2019) considers estuarine volumes separated from the ocean by transects fixed in space. The selection of the transect locations is typically based on geographical considerations with preference to river mouths or sills, that is, locations of specifically high turbulent mixing and high variability. All these estimates require the calculation of inflowing and outflowing salinity values, and some need additional calculation of the representative inflow and outflow values of the squared salinity. These values can only be calculated by means of a careful analysis of numerical model results, as described by Burchard et al. (2019). The simple estimate by MacCready et al. (2018), see Eq. (2), is useful but not exact, with typical errors of about 10% for idealized examples (Burchard et al. 2019). The estimates considering inflow and outflow of squared salinity are exact, but more complex and not intuitive.

In contrast, the new mixing law (4) for long-term averaged estuaries based on isohaline volumes is exact (when neglecting freshwater fluxes due to evaporation, precipitation or exchange with the groundwater) and easy to calculate. Moreover, it is universal in the sense that the resulting mixing does not depend on the internal dynamics of the estuary itself, but only on the selected isohaline confining the estuary. With a confining salinity of *S* = 30 g kg^{−1}, and an average river runoff of *Q*_{r} = 1000 m^{3} s^{−1}, the long-term estuarine mixing is exactly *M* = 900.000 m^{3} s^{−1} (g kg^{−1})^{2}, a calculation not requiring a supercomputer.

Other than for the Knudsen-based mixing estimates, the universal law is formulated as an analytical function, such that the derivative (with respect to *S*) can be taken, allowing to calculate the mixing per isohaline volume. The result is that the long-term mixing per salinity class increases linearly with salinity in all estuaries. Expressing the mixing per salinity class as the product of isohaline volume and average mixing per salinity class opens perspectives for interpreting and predicting the behavior of any estuary, as discussed in section 5. One example is that a decrease (increase) in local mixing leads to an increase (decrease) in isohaline volume, for example, by flattening (erecting) the isohaline surfaces or by moving the isohaline volume to a wider (narrower) cross section.

Viewing an estuary in isohaline space rather than in Eulerian space allows for a hydrodynamic interpretation rather than a geographical interpretation, considering the estuary and its river plume as a continuum not separated in its most dynamic zone at the mouth of the estuary. In their isohaline theory, MacCready et al. (2002) differentiate between low-salinity classes remaining inside the estuary all the time, midsalinity classes being partially ejected out of the estuary and high-salinity classes being injected from the ocean. For these three classes creation and destruction mechanisms by river runoff and isohaline drift due to mixing are different. The present theory is simpler in the sense that no difference is made in terms of the value of the salinity itself, and no ejection or injection of isohalines is considered since the estuary has no limit and continues into the river plume.

This theoretical paper intends to present the derivation of the new universal law of estuarine mixing. In future studies, the new theory needs to be empirically challenged for a range of diverse estuaries covering the relevant parameter space as for example proposed by Geyer and MacCready (2014). Such empirical studies will require well-calibrated numerical models validated by means of observational data. Studies entirely based on observations are not promising due to the inherent sparseness of data especially for measurements of turbulent mixing. It is expected that for sufficiently long averaging times, using models analyzing physical as well as numerical mixing, the universal law will be confirmed. One key question to be answered is how long the averaging period needs to be. For energetic tidal estuaries averaging over a spring–neap cycle might be sufficient, if the wind and freshwater forcing do not vary too much. Weakly tidal estuaries such as fjords might require much longer averaging periods. Such comparative studies should also investigate adjustment time scales of actual mixing toward the universal law. It is expected that the new theory in addition to other estuarine theories such as TEF helps to gain a deeper understanding of estuarine processes and their time scales.

## Acknowledgments

This paper is a contribution to the Collaborative Research Centre TRR 181 on Energy Transfers in Atmosphere and Ocean and the Research Training Group Baltic TRANSCOAST GRK 2000 funded by the German Research Foundation. The author is grateful to Rainer Feistel, Knut Klingbeil, Xaver Lange, Marvin Lorenz, and Lars Umlauf (all from IOW) and Carsten Eden (University of Hamburg) for critically reading an earlier version of this manuscript and to Parker MacCready (University of Washington) for constructive comments on the submitted version. I am also grateful to Sarah Giddings (Scripps Institution of Oceanography) and two anonymous reviewers for useful comments during the review process.

## REFERENCES

## Footnotes

© 2019 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).