Estuarine mixing

This review paper presents, explains and discusses major aspects of estuarine mixing which is defined as the destruction of salinity variance. Due to the large amounts of brackish water in estuaries produced by mixing of fresh river discharge and salty ocean water, mixing is one major characteristic of what is an estuary. In this review, mixing is quantified locally as well as on estuary-wide scales. Diagnostics of integrated mixing are given for estuarine volumes bounded by transects as well as isohalines (surfaces of constant salinity) moving with the flow. It is shown how entrainment across a moving isohaline surface depends on gradients of turbulent salt flux and mixing per salinity class. Various relations are derived that link estuarine salt mixing to other estuarine properties such as the freshwater discharge and the bulk estuarine circulation. For estuaries bounded towards the ocean by a fixed transect, the Knudsen mixing law is explained, where estuarine mixing is the product of the Knudsen salinities of inflowing and outflowing water masses and the river discharge. When the estuarine volume is bounded by a moving isohaline surface of salinity S, mixing inside the estuary is simply the product of S² and the river discharge. Major processes that drive estuarine mixing are presented on various time scales (tidal, fortnightly, weather and discharge time scales) and spatial scales (channel-shoal interaction, mixing fronts). As underlying methods for the quantification of mixing, observational concepts, as well as numerical modelling methods such as consistent turbulence closure modelling and numerical mixing analyses are presented. As an outlook, some future perspectives are sketched. Many of the concepts presented in this review are illustrated using simulation results from a numerical model setup of the Elbe River estuary.

1 Introduction

Estuaries are semi-enclosed coastal water bodies where riverine freshwater run-off from land is mixed with offshore salty ocean water to produce brackish water masses of intermediate salinities which are ejected offshore into the coastal ocean. In this sense, estuaries can be characterised as mixing machines (MacCready and Banas, 2011; Wang et al., 2017) with mixing rates far greater than in other parts of the ocean. Salt is the most characteristic constituent that is mixed in estuaries, because of (i) its significantly different concentration between rivers (typically <0.5 g kg⁻¹) and the adjacent ocean (typically >30 g kg⁻¹) and (ii) its inert character with no internal sinks and sources and no fluxes through the surface and the bottom. In addition to salt, there is a number of further properties, e.g., nutrients and pollutants, that are distinct between rivers and the ocean and which can be mixed in estuaries. This makes mixing a fundamental process in estuaries. We therefore follow here the definition of an estuary by Pritchard (1967) who stated An estuary is a semi-enclosed coastal body of water which has a free connection with the open sea and within which sea water is measurably diluted with fresh water derived from land drainage. Instead of diluted we would prefer to say mixed, which effectively has the same meaning but highlights mixing as the defining process of estuaries. Water bodies that follow this principle are classical estuaries in the sense that their functioning is based on a net freshwater water supply. The contrasting case is given by inverse estuaries that are based on a net freshwater deficit due to evaporation, leading to the export of hypersaline water masses. The present review focusses on classical estuaries (in the following just denoted as estuaries for simplicity), whereas inverse estuaries are only occasionally discussed as contrasting systems. Fundamental concepts of estuaries and estuarine circulation have already been covered by previous reviews (MacCready and Geyer, 2010; Geyer and MacCready, 2014), such that we here focus on mixing in estuaries and its physical and ecological consequences.

Much of the estuarine literature focusses on mixing, starting with Knudsen’s classic paper (Knudsen, 1900), in which he states: As the freshwater spreads out over the seawater it mixes with it so that the salinity of the surface increases seawards (translation by Burchard et al., 2018a). Fischer’s review (Fischer, 1976) is entitled Mixing and dispersion in estuaries, highlighting its fundamental importance. Notwithstanding the attention focused on mixing throughout the estuarine literature, the actual meaning of the term mixing has often been vague or ambiguous: everyone is familiar with mixing via the daily-life experience of pouring milk into a cup of tea and using a tea spoon to mix it. However, amid the complex and multi-scale processes in estuaries, that simple concept is overwhelmed by consideration of larger scale processes associated with turbulence, shear dispersion and buoyancy flux. Turbulence, diffusion, dispersion, buoyancy flux and mixing are often loosely treated as synonyms, leading to confusion as to what we mean by mixing.

In this review, we come back to mixing in a cup (or glass beaker), or more precisely to the thermodynamic definition of mixing, which is the destruction of variance of some scalar quantity (Gibbs, 1878). In the estuarine context, we focus on the destruction of salinity variance, which is defined in the oceanic turbulence literature as χ_s, due to the down-gradient diffusion of salt by molecular diffusivity

$$\begin{array}{}\text{(1)}& {\mathit{\chi }}_{\mathrm{s}}=\mathrm{2}\mathit{\kappa }\left[{\left({\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{s}}{\partial x}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{s}}{\partial y}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{s}}{\partial z}}\right)}^{\mathrm{2}}\right]\end{array}$$

(Nash and Moum, 2002; Burchard and Rennau, 2008). In Eq. (1), κ is the molecular diffusivity of salt, $\stackrel{\mathrm{̃}}{s}$ is the turbulent salinity fluctuation, and square brackets denote Reynolds averaging (see Sect. 2.1 for details). In the turbulent environments of estuaries, this molecular process occurs at sub-millimetre scales, but it is the direct result of the turbulent and shearing motions acting on the salinity gradients at a wide range of scales extending all the way up to the horizontal dimensions of estuaries. It should be noted that the process of molecular mixing is irreversible inside the water body, such that χ_s≥0. Negative mixing or un-mixing can however occur at the sea surface when evaporation takes place (Yu, 2010; Klingbeil and Lorenz, 2025) or sea ice is produced due to freezing (Notz and Worster, 2009) and in desalination plants to extract freshwater from salt water (reverse osmosis, Kim et al., 2019). In numerical models, un-mixing can also occur due to discretisation errors of advection schemes (Henell et al., 2023). Maintaining the strict thermodynamic definition of mixing turns out to be a powerful approach to examining estuarine processes, because salinity variance can be defined locally, as is often done in the turbulence literature, as well as globally, at the overall scale of the estuary. We do not question the importance of other concepts, such as vertical buoyancy flux and horizontal dispersion, but in this review we retain the strict definition of mixing to explore the processes responsible for its occurrence in estuaries as well as its quantitative relationship to estuarine exchange flow.

Throughout this review, exemplary data from a numerical simulation of the Elbe River estuary in northern Germany are used to demonstrate the different mixing theories. The Elbe River estuary is an elongated meso-tidal estuary with one major discharge source at the landward end for which several studies of estuarine mixing have been carried out (Reese et al., 2024, 2026; Burchard et al., 2025). A brief introduction into the Elbe River estuary is given in Appendix C.

This review is structured as follows: After this introduction into the topic of estuarine mixing (Sect. 1), the existing theories on estuarine mixing are defined and discussed (Sect. 2). This section is structured into micro-structure mixing (Sect. 2.1) and parameterised mixing (Sect. 2.2), where Reynolds decomposition and turbulence closure assumptions are applied. The mixing definitions from Sect. 2 as well as the Total Exchange Flow analysis framework (Sect. 3.2) will be used in Sect. 3 (Estuarine Circulation and Mixing) to quantify mixing in entire estuaries. Before discussing estuarine mixing, we give a brief introduction to estuarine hydrodynamics (Sect. 3.1). For the mixing quantification, we first present the theory using fixed transects (Knudsen theories, see Sect. 3.3). Water Mass Transformation (WMT) theories (Sect. 3.4) are explained from which mixing laws for estuarine volumes bounded by isohaline surfaces as the seaward boundary are derived (Sect. 3.5). Section 3 concludes with some remarks on the relation between estuarine mixing and estuarine circulation (Sect. 3.6) and mixing of constituents other than salt (Sect. 3.7). While Sects. 2 and 3 focus on the definition and discussion of mixing, Sect. 4 gives examples for the most important estuarine processes that drive mixing. Those processes are related to single tides (Sect. 4.1.1), the spring-neap cycle (Sect. 4.1.2), time scales of river discharge and meteorological forcing (Sect. 4.1.3), channel-shoal interaction (Sect. 4.2.1) and mixing at fronts (Sect. 4.2.2). As methods to help quantifying mixing in estuaries, techniques to observe estuarine mixing are introduced (Sect. 5.1) and numerical modelling techniques are presented (Sect. 5.2), with a focus on turbulence closure modelling (Sect. 5.2.1) and numerical mixing (Sect. 5.2.2). Finally, some future perspectives are discussed in Sect. 6. The extensive appendix contains details about an analytic illustrative example for small-scale mixing (Appendix A), some key budget equations related to salinity variance (Appendix B), information about the Elbe River estuary model used to provide estuarine mixing examples (Appendix C), a derivation of the coordinate transformation of the vertical salinity equation (Appendix D), an explanation for the calibration of two-equation turbulence closure models (Appendix E), and derivations for the numerical mixing example (Appendix F). At the end of the Appendix, a table with the most important variables, their definitions, units and defining equations is presented (Table ).

2 Quantification of mixing

While mixing occurs on the micro-scale only, its integral effects are most prominently effective on the large, estuarine scale. We therefore start our explanations with the quantification of local stirring and mixing. This will first be based on molecular diffusion on the micro-scale and Reynolds averaging on the macro-scale (Sect. 2.1) and then parameterised by means of turbulence closures as it would be calculated in numerical models of estuaries (Sect. 2.2).

2.1 Micro-structure mixing

Mixing of a tracer s (for which we use salinity here as an example) occurs at the micro-scale when tracer gradients are reduced by molecular diffusion, following the Fickian law,

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial \stackrel{\mathrm{ˇ}}{s}}{\partial t}}+{\displaystyle \frac{\partial }{\partial x_j}}\left({\stackrel{\mathrm{ˇ}}{u}}_j\stackrel{\mathrm{ˇ}}{s}-\mathit{\kappa }{\displaystyle \frac{\partial \stackrel{\mathrm{ˇ}}{s}}{\partial x_j}}\right)=\mathrm{0},\end{array}$$

where the Einstein summation convention $a_j{b}_j={\sum }_{j=\mathrm{1}}^{\mathrm{3}}{a}_j{b}_j=a_{\mathrm{1}}{b}_{\mathrm{1}}+a_{\mathrm{2}}{b}_{\mathrm{2}}+a_{\mathrm{3}}{b}_{\mathrm{3}}$ has been applied. In Eq. (2), $\stackrel{\mathrm{ˇ}}{s}$ is the instantaneous tracer concentration $\stackrel{\mathrm{ˇ}}{s}=\left[s\right]+\stackrel{\mathrm{̃}}{s}$ with the Reynolds-averaged tracer concentration [s] and the fluctuating component $\stackrel{\mathrm{̃}}{s}$, κ is the molecular diffusivity of salinity, and $\stackrel{\mathrm{ˇ}}{u}={\stackrel{\mathrm{ˇ}}{u}}_{\mathrm{1}}$ and $\stackrel{\mathrm{ˇ}}{v}={\stackrel{\mathrm{ˇ}}{u}}_{\mathrm{2}}$ are the horizontal and $\stackrel{\mathrm{ˇ}}{w}={\stackrel{\mathrm{ˇ}}{u}}_{\mathrm{3}}$ is the vertical velocity component. The terms in the brackets on the left hand side of Eq. (2) are the advective and the molecular diffusive fluxes, the divergence of which determines the change of the salinity distribution. This transport equation determines the salinity distribution on all scales ranging from the sub-millimetre scales of molecular diffusion to the global scales of meridional overturning circulation, including scales of estuarine mixing.

In turbulence theory, the Reynolds average (also called ensemble average) is defined as the average of an infinite number of macroscopically identical but microscopically different flow realisations, where the turbulent random fluctuations are averaged out (Lesieur, 2008). Consequently, $\left[\stackrel{\mathrm{̃}}{s}\right]=\mathrm{0}$. In estuarine physics, and similarly in most fields of larger-scale oceanography, the Reynolds-averaged rather than the instantaneous properties of the flow are considered. Field observations of tracer concentrations (e.g., from Conductivity-Temperature-Depth (CTD) probes) as well as numerical model results are supposed to represent Reynolds-averaged quantities. A continuity equation (incompressibility condition) is used in most ocean models:

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\partial {\stackrel{\mathrm{ˇ}}{u}}_i}{\partial x_i}}=\mathrm{0},{\displaystyle \frac{\partial \left[u_i\right]}{\partial x_i}}=\mathrm{0},{\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{u}}_i}{\partial x_i}}=\mathrm{0},\end{array}$$

applying to instantaneous and thus to Reynolds-averaged and fluctuating velocity fields. Using Eq. (2), a dynamic equation for the Reynolds-averaged salinity can be derived:

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\partial \left[s\right]}{\partial t}}+{\displaystyle \frac{\partial }{\partial x_j}}\left(\left[u_j\right]\left[s\right]+\left[{\stackrel{\mathrm{̃}}{u}}_j\stackrel{\mathrm{̃}}{s}\right]-\mathit{\kappa }{\displaystyle \frac{\partial \left[s\right]}{\partial x_j}}\right)=\mathrm{0},\end{array}$$

with the advective tracer flux [u_j][s], the turbulent tracer flux $\left[{\stackrel{\mathrm{̃}}{u}}_j\stackrel{\mathrm{̃}}{s}\right]$ and the diffusive tracer flux $-\mathit{\kappa }\partial \left[s\right]/\partial x_j$. Based on Eq. (2), it is also possible to derive an equation for the micro-structure tracer variance $\left[{\stackrel{\mathrm{̃}}{s}}^{\mathrm{2}}\right]$:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& {\displaystyle \frac{\partial \left[{\stackrel{\mathrm{̃}}{s}}^{\mathrm{2}}\right]}{\partial t}}+{\displaystyle \frac{\partial }{\partial x_j}}\left(\left[u_j\right]\left[{\stackrel{\mathrm{̃}}{s}}^{\mathrm{2}}\right]+\left[{\stackrel{\mathrm{̃}}{u}}_j{\stackrel{\mathrm{̃}}{s}}^{\mathrm{2}}\right]-\mathit{\kappa }{\displaystyle \frac{\partial \left[{\stackrel{\mathrm{̃}}{s}}^{\mathrm{2}}\right]}{\partial x_j}}\right)\\ & =\underset{P_{\mathrm{s}}}{\underbrace{-\mathrm{2}\left[{\stackrel{\mathrm{̃}}{u}}_j\stackrel{\mathrm{̃}}{s}\right]{\displaystyle \frac{\partial \left[s\right]}{\partial x_j}}}}-\underset{{\mathit{\chi }}_{\mathrm{s}}}{\underbrace{\mathrm{2}\mathit{\kappa }\left[{\left({\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{s}}{\partial x_j}}\right)}^{\mathrm{2}}\right]}}\end{array}\end{array}$$

(see, e.g., Eq. (5) by Mellor and Yamada, 1974). In Eq. (5), P_s quantifies the production of micro-structure variance due to turbulent stirring (with the Reynolds-averaged tracer gradient vector $\partial \left[s\right]/\partial x_j$) and χ_s represents destruction of microstructure variance due to molecular mixing.

Multiplying Eq. (4) by 2[s] gives a transport equation for the square of the Reynolds-averaged tracer:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}& {\displaystyle \frac{\partial [s{]}^{\mathrm{2}}}{\partial t}}+{\displaystyle \frac{\partial }{\partial x_j}}\left(\left[u_j\right][s{]}^{\mathrm{2}}+\mathrm{2}[{\stackrel{\mathrm{̃}}{u}}_j\stackrel{\mathrm{̃}}{s}\left]\right[s]-\mathit{\kappa }{\displaystyle \frac{\partial [s{]}^{\mathrm{2}}}{\partial x_j}}\right)\\ & =\underset{-P_{\mathrm{s}}}{\underbrace{\mathrm{2}\left[{\stackrel{\mathrm{̃}}{u}}_j\stackrel{\mathrm{̃}}{s}\right]{\displaystyle \frac{\partial \left[s\right]}{\partial x_j}}}}-\mathrm{2}\mathit{\kappa }{\left({\displaystyle \frac{\partial \left[s\right]}{\partial x_j}}\right)}^{\mathrm{2}},\end{array}\end{array}$$

where on the right-hand side the stirring term P_s appears as a sink term besides a destruction term due to molecular diffusivity. In contrast to Eq. (5), where the destruction of micro-structure variance occurs due to molecular diffusivity acting on micro-structure gradients $\partial \stackrel{\mathrm{̃}}{s}/\partial x_j$, in Eq. (6) the molecular diffusivity acts on the much smaller Reynolds-averaged gradient $\partial \left[s\right]/\partial x_j$ such that this term is generally negligible. This means that variance is first transferred from the Reynolds-averaged regime of Eq. (6) to the turbulent regime Eq. (5), where it is then dissipated. In turbulence closure modelling typically P_s=χ_s (see, e.g., Eq. (31) by Mellor and Yamada, 1974) is applied such that stirring equals mixing, following a local equilibrium assumption. More details on turbulence closure modelling suitable for estuaries are given in Sect. 5.2.

Figure 1 Stirring and mixing in a glass beaker: (a) Evolution of salinity $\stackrel{\mathrm{ˇ}}{s}$ for the case of no stirring, (b) little stirring, and (c) strong stirring. The initial distribution after the stirring is shown as blue lines, and the distribution after 10 min is shown as red lines. The parameters for the problem are chosen as height of fluid inside the one-dimensional glass beaker D=0.1 m, and the molecular diffusivity of salinity, $\mathit{\kappa }={\mathrm{10}}^{-\mathrm{9}}$ m² s⁻¹. The two black bars mark the area over which the local variance is estimated. More details are given in Appendix A .

In short: micro-structure tracer variance is produced by stirring P_s (increase of local micro-structure gradients due to turbulent eddies) and dissipated by mixing χ_s, while the divergence term on the left-hand side just spatially redistributes the micro-structure tracer variance. Note that the stirring term is twice the product of the turbulent flux times the Reynolds-averaged tracer gradient. The stirring term P_s is typically positive, since the Reynolds-averaged salinity gradient and the turbulent salt flux have opposite signs due to the generally down-gradient property of turbulent fluxes (classical exceptions occur in convective boundary layers, see e.g. Legay et al., 2025).

This can be explained by a laboratory experiment with salt mixing. A simple idealised model of this is given in Appendix A and results are shown in Fig. 1. After having carefully pumped saltwater of 30 g kg⁻¹ underneath freshwater (with some continuous mixing), the local salinity variance is mostly low: sufficiently small control volumes would contain water with a small salinity range (Fig. 1a). Introduction of turbulence by means of a spoon will lead to stirring (increase of P_s), such that local control volumes (marked as the area between the two black bars in Fig. 1) will contain streaks of saltwater at various ranges, with sharp gradients between them, such that the local variance $\left[{\stackrel{\mathrm{̃}}{s}}^{\mathrm{2}}\right]$ is increased. This is demonstrated as initial conditions for moderate stirring (Fig. 1b) and strong stirring (Fig. 1c). Now, mixing χ_s will be moderately or strongly enhanced due to the small (and constant) molecular diffusivity κ acting on the strong micro-scale gradients $\partial \stackrel{\mathrm{̃}}{s}/\partial z$ (which are squared in the mixing term). At the end of this process, the salt will be almost fully diluted in the water such that local variance becomes small, with $\stackrel{\mathrm{̃}}{s}\to \mathrm{0}$ and $\stackrel{\mathrm{ˇ}}{s}\to \left[s\right]=\mathrm{15}$ g kg⁻¹. Further introduction of turbulence by a spoon will not lead to further stirring (and thus not to further mixing), because the tracer gradients have vanished.

In real estuaries, stirring typically occurs due to vertical shear instabilities driven by tidal flow, generating large eddies as shown in Figs. 17 and 18 in Sect. 5.1. Via the classical turbulence downward cascade, smaller and smaller eddies are generated such that finally mixing is enhanced at the smallest scales. Since estuaries are typically narrow and friction-dominated, horizontal instabilities on the submesoscale (McWilliams, 2016) play a minor role for stirring. An exception would be large fjord-type estuaries with weak tides such as the Baltic Sea (Chrysagi et al., 2021).

Direct in-situ measurements of salinity mixing χ_s are difficult to obtain due to the small value of the molecular salinity diffusivity of $\mathit{\kappa }=\mathcal{O}({\mathrm{10}}^{-\mathrm{9}}$ m² s⁻¹) and the consequently strong gradients at small scales, but successful attempts have been reported by Nash and Moum (2002) for locations on the continental shelf. According to these authors the salinity-gradient spectrum peaks at dissipative scales ten times smaller than the temperature-gradient spectrum, such that most salinity variance decay occurs in the sub-millimetre range. Therefore, and because of estuaries having generally higher levels of turbulence than continental shelves, direct observations of χ_s in estuaries are not feasible, and indirect observations are needed (see Sect. 5.1). Instead of using turbulence observations, mixing in estuaries is mostly studied by means of well-calibrated fine-resolution numerical models equipped with accurate numerical discretisations and physically based turbulence closures (Sect. 5.2).

2.2 Mixing at resolved local scales

Whereas the irreversible process of mixing happens at very small scales, the quantification of mixing is accomplished both in observations and in models through the application of turbulence closure assumptions (Mellor and Yamada, 1974; Peters and Bokhorst, 2001; Umlauf and Burchard, 2005). On the level of numerical ocean models, the turbulent fluxes are typically parameterised by means of the eddy diffusivity assumption, resulting in down-gradient turbulent tracer fluxes:

$$\begin{array}{}\text{(7)}& \left[\stackrel{\mathrm{̃}}{u}\stackrel{\mathrm{̃}}{s}\right]=-K_{\mathrm{h}}{\displaystyle \frac{\partial s}{\partial x}};\left[\stackrel{\mathrm{̃}}{v}\stackrel{\mathrm{̃}}{s}\right]=-K_{\mathrm{h}}{\displaystyle \frac{\partial s}{\partial y}};\left[\stackrel{\mathrm{̃}}{w}\stackrel{\mathrm{̃}}{s}\right]=-K_{\mathrm{v}}{\displaystyle \frac{\partial s}{\partial z}},\end{array}$$

now using salinity s=[s] as Reynolds-averaged tracer concentration. In Eq. (7), K_h is the horizontal eddy diffusivity and K_v is the vertical eddy diffusivity. With this the Reynolds-averaged salinity budget equation Eq. (4) becomes:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& \underset{\text{change}}{\underbrace{{\displaystyle \frac{\partial s}{\partial t}}}}+\underset{\text{advection}}{\underbrace{{\displaystyle \frac{\partial \left(us\right)}{\partial x}}+{\displaystyle \frac{\partial \left(vs\right)}{\partial y}}+{\displaystyle \frac{\partial \left(ws\right)}{\partial z}}}}\\ & \underset{\text{diffusion}}{\underbrace{-{\displaystyle \frac{\partial }{\partial x}}\left(K_{\mathrm{h}}{\displaystyle \frac{\partial s}{\partial x}}\right)-{\displaystyle \frac{\partial }{\partial y}}\left(K_{\mathrm{h}}{\displaystyle \frac{\partial s}{\partial y}}\right)-{\displaystyle \frac{\partial }{\partial z}}\left(K_{\mathrm{v}}{\displaystyle \frac{\partial s}{\partial z}}\right)}}=\mathrm{0},\end{array}\end{array}$$

showing that salinity changes are exclusively determined by the divergence of advective and turbulent fluxes. Note that in Eq. (8) molecular tracer fluxes have been neglected. With Eq. (7) the production of micro-structure variance due to stirring becomes

$$\begin{array}{}\text{(9)}& P_{\mathrm{s}}=\mathrm{2}{K}_{\mathrm{h}}{\left({\displaystyle \frac{\partial s}{\partial x}}\right)}^{\mathrm{2}}+\mathrm{2}{K}_{\mathrm{h}}{\left({\displaystyle \frac{\partial s}{\partial y}}\right)}^{\mathrm{2}}+\mathrm{2}{K}_{\mathrm{v}}{\left({\displaystyle \frac{\partial s}{\partial z}}\right)}^{\mathrm{2}}\stackrel{!}{=}{\mathit{\chi }}_{\mathrm{s}},\end{array}$$

where in the last step stirring and mixing of micro-structure salinity variance are set equal, which is a typical assumption in turbulence closure modelling (see Sect. 2.1). The local variance decay χ_s is used as a local measure for the mixing of Reynolds-averaged salinity (Burchard and Rennau, 2008). This local equilibrium assumption is generally valid on the temporal and spatial scales that are resolved by numerical ocean models. In detail, χ_s appears as a sink term in both the salinity variance and salinity-square budgets. The corresponding derivation is shown in Appendix B.

In estuaries, the vertical term of the salinity variance decay Eq. (9) typically dominates over the horizontal terms due to the dominance of tidally-driven vertical shear (Li et al., 2024), such that we obtain

$$\begin{array}{}\text{(10)}& {\mathit{\chi }}_{\mathrm{s}}\approx {\mathit{\chi }}_{\mathrm{s},\mathrm{v}}=\mathrm{2}{K}_{\mathrm{v}}{\left({\displaystyle \frac{\partial s}{\partial z}}\right)}^{\mathrm{2}}.\end{array}$$

This is the case because in estuaries horizontal turbulent transports and their divergences are small compared to the vertical transports, with the consequence that in estuarine models the horizontal diffusion is often neglected in the parameterised salinity budget equation Eq. (8) as well as in the salinity variance equation Eq. (B1). Due to this dominance of vertical processes in estuaries, it is instructive to study the balance of the vertical variance,

$$\begin{array}{}\text{(11)}& {s_{\mathrm{v}}^{\prime }}^{\mathrm{2}}={\left(s-{\displaystyle \frac{\mathrm{1}}{\mathit{\eta }+H}}\underset{-H}{\overset{\mathit{\eta }}{\int }}s\text{d}z\right)}^{\mathrm{2}}={\left(s-{\overline{s}}_{\mathrm{v}}\right)}^{\mathrm{2}},\end{array}$$

for the vertical integral of which the following budget equation can be derived from Eq. (8), see Li et al. (2018) for details:

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& \underset{\text{rate of change}}{\underbrace{{\displaystyle \frac{\partial }{\partial t}}\underset{-H}{\overset{\mathit{\eta }}{\int }}{s_{\mathrm{v}}^{\prime }}^{\mathrm{2}}\text{d}z}}\underset{\text{ advection}}{\underbrace{+{\displaystyle \frac{\partial }{\partial x}}\underset{-H}{\overset{\mathit{\eta }}{\int }}u{s_{\mathrm{v}}^{\prime }}^{\mathrm{2}}\text{d}z+{\displaystyle \frac{\partial }{\partial y}}\underset{-H}{\overset{\mathit{\eta }}{\int }}v{s_{\mathrm{v}}^{\prime }}^{\mathrm{2}}\text{d}z}}\\ & =\underset{\text{horizontal straining}}{\underbrace{-\mathrm{2}\underset{-H}{\overset{\mathit{\eta }}{\int }}{u}_v^{\prime }{s}_{\mathrm{v}}^{\prime }{\displaystyle \frac{\partial {\overline{s}}_{\mathrm{v}}}{\partial x}}\text{d}z-\mathrm{2}\underset{-H}{\overset{\mathit{\eta }}{\int }}{v}_{\mathrm{v}}^{\prime }{s}_{\mathrm{v}}^{\prime }{\displaystyle \frac{\partial {\overline{s}}_{\mathrm{v}}}{\partial y}}\text{d}z}}\\ & -\underset{\text{vertical mixing}}{\underbrace{\underset{-H}{\overset{\mathit{\eta }}{\int }}{\mathit{\chi }}_{\mathrm{s},\mathrm{v}}\text{d}z}},\end{array}\end{array}$$

where η is the surface elevation. Equation (12) shows that the vertical variance balance is time-dependent and spatially variable. In contrast to the total salinity variance budget, the vertical salinity variance budget has source terms, the so-called horizontal straining terms, representing the conversion of horizontal variance (${s_{\mathrm{h}}^{\prime }}^{\mathrm{2}}$, where $s_{\mathrm{h}}^{\prime }=s_{\mathrm{tot}}^{\prime }-s_{\mathrm{v}}^{\prime }={\overline{s}}_{\mathrm{v}}-{\overline{s}}_{\mathrm{tot}}$) associated with the horizontal salinity gradient $\partial {\overline{s}}_{\mathrm{v}}/\partial x$ to vertical variance (Simpson et al., 1990). Note that horizontal straining is split into longitudinal straining and lateral straining (first and second term of horizontal straining, respectively) and can be a source (mainly during ebb straining) and a sink (mainly during flood straining) of vertical variance, whereas the effect of vertical mixing is always to reduce the vertical variance. According to Li et al. (2018), estuarine mixing is driven in a three-step process: First, horizontal variance is provided to the estuary by means of boundary variance transports from the river and the ocean, through the boundary transport term in Eq. (B2). Then the horizontal straining term in Eq. (12) converts horizontal variance into vertical variance, which is then in a third step mixed away by the vertical mixing term in Eq. (12).

3 Estuarine circulation and mixing

After a short introduction into the basics of estuarine hydrodynamics in Sect. 3.1, we here introduce mixing concepts for entire estuaries, first following the classical exchange flow theory proposed by Knudsen (1900), see Sect. 3.3, which can be quantified by using the Total Exchange Flow (TEF) analysis framework across fixed transects (Sect. 3.2). After introducing local isohaline theory (Sect. 3.4), we show how to analyse mixing in estuarine volumes bounded by an isohaline instead of a fixed transect (Sect. 3.5). Based on the local isohaline theory, the quantification of estuarine circulation is directly related to mixing (Sect. 3.6). Finally, mixing of constituents other than salt is briefly discussed (Sect. 3.7).

3.1 Basics of estuarine hydrodynamics

Although this review paper mainly focusses on estuarine salt mixing, a short introduction into the hydrodynamics determining the salt distribution and the turbulence available for salinity mixing will be given here. More detailed reviews can be found in MacCready and Geyer (2010) and Geyer and MacCready (2014).

Estuaries are characterised by a longitudinal salinity gradient extending from the region of the estuarine mouth with values typically close to ocean salinity (somewhat below 35 g kg⁻¹) to values of river salinity (typically below 0.5 g kg⁻¹) in the freshwater range. Since saline water is more dense than fresh water, this salinity gradient causes a longitudinal density gradient which results in a longitudinal pressure gradient in the momentum balance, with a barotropic and a baroclinic term. While the baroclinic term is zero at the surface and increases continuously towards the bottom (driving water in up-estuarine direction, with intensification at the bottom), the barotropic pressure-gradient term is independent of the vertical position and drives water out of the estuary. In addition oscillating tides provide small-scale turbulence resulting in vertical shear stress divergence and consequently in diffusion of the longitudinal velocity profiles. The combination of these forces results in a gravitationally-driven tidally averaged exchange flow, with near-bottom flow directed in up-estuarine direction and near-surface flow directed in down-estuarine direction, the so-called gravitational circulation. The strength of the stratifying gravitational circulation depends on the ratio of the gravitational forces due to the longitudinal density gradient and the de-stratifying vertical shear stress divergence. This ratio is expressed as the Simpson number

$$\begin{array}{}\text{(13)}& Si={\displaystyle \frac{D^{\mathrm{2}}\left|b_x\right|}{〈u_{*}^{\mathrm{2}}〉}}\end{array}$$

(Simpson et al., 1990; Monismith et al., 1996; Stacey et al., 2008) with the water depth D, the horizontal buoyancy gradient $b_x=\partial b/\partial x$, the buoyancy $b=-g/{\mathit{\rho }}_{\mathrm{0}}(\mathit{\rho }-{\mathit{\rho }}_{\mathrm{0}})$, the potential density ρ, the reference density ρ₀, the gravitational acceleration g, and the bottom friction velocity u_*. There are several other hydrodynamic processes contributing to estuarine exchange flow. One essential process is tidal straining (Simpson et al., 1990): During flood saltier ocean water is sheared over less salty estuarine water, such that the water column becomes statically unstable and thus highly turbulent such that properties are vertically homogenised. During ebb the opposite occurs, resulting in stable stratification and suppression of turbulence. This asymmetry of turbulence does not only affect the salt distribution, but also the velocity profiles: During flood up-estuarine momentum is transported downwards in a much stronger amount than down-estuarine momentum is transported downwards during ebb, which in a tidal average leads to an up-estuarine residual flow near the bottom (Jay and Musiak, 1994). This process has also been named ESCO (eddy-viscosity – shear covariance) by Dijkstra et al. (2017). In an idealised model study Burchard and Hetland (2010) showed that in tidally energetic flows the contribution of ESCO to estuarine circulation could be stronger than that of gravitational circulation. Other important hydrodynamic processes generating estuarine circulation are lateral circulation (Lerczak and Geyer, 2004), estuarine convergence (Ianniello, 1979; Burchard et al., 2014) and wind straining (Scully et al., 2005). It has been shown that gravitational circulation, ESCO, and lateral circulation are strongly scaling with Si (Burchard et al., 2011; Lange and Burchard, 2019). Since Si is larger during neap tide than during spring tide (due to smaller u_*), the estuarine exchange flow is expected to be stronger during neap tide. The estuarine circulation in concert with vertical mixing is a major process carrying salt into the estuary, against the river discharge. Often, this so-called shear dispersion process is parameterised by a horizontal diffusion term in the longitudinal salt balance equation (MacCready, 2004).

The Simpson number Si does also have a strong influence on the stratification in an estuary. For Si<0.2 (tidally energetic), the water column would be mixed throughout the tidal cycle. For Si>1 (weak tidal energy), stratification should be maintained during the entire tidal cycle. For intermediate situations with $\mathrm{0.2}\le Si\le \mathrm{1}$, however, the water column should stratify during ebb and destratify during flood, leading to strain-induced periodic stratification (SIPS, Simpson et al., 1990; Verspecht et al., 2009). Since mixing is typically proportional to the square of the vertical salinity gradient, see Eq. (10), stratification has a major impact on estuarine mixing. Specific examples of physical drivers of estuarine mixing are discussed in Sect. 4.

3.2 Total Exchange Flow

The estuarine exchange flow of water masses defined by salinity, i.e., the net inflow of high salinity ocean waters and the net outflow of low salinity estuarine waters, can be best quantified in terms of time-averaged transports in fixed salinity classes. The resulting Total Exchange Flow (TEF) provides an analysis framework based on salinity coordinates rather than geopotential (z-)coordinates which is consistently linked to the Knudsen theory as well as to estuarine mixing. As shown by several authors (MacCready, 2011; Sutherland et al., 2011; Burchard et al., 2018a), the Eulerian (z-coordinate) framework could be mapped back to time-averaged salinities, but the resulting exchange flow profiles would significantly underestimate the exchange flow. Therefore, the TEF framework has developed into a major research tool for analysing estuarine dynamics. For the Baltic Sea, approaches similar to TEF had already been developed earlier (Walin, 1977; Döös et al., 2004). Here, we briefly explain the theoretical framework for TEF and refer to the literature for the details (MacCready, 2011; Burchard et al., 2018a). Given a fixed transect T across an estuary, the time-averaged volume, salt and salt-squared transports across the transect for all salinities >S are defined as

$$\begin{array}{}\text{(14)}& \begin{array}{rl}Q\left(S\right)& =〈\underset{A\left(S\right)}{\int }u\text{d}A〉,Q_{\mathrm{s}}\left(S\right)=〈\underset{A\left(S\right)}{\int }us\text{d}A〉,\\ & Q_{{\mathrm{s}}^{\mathrm{2}}}\left(S\right)=〈\underset{A\left(S\right)}{\int }u{s}^{\mathrm{2}}\text{d}A〉,\end{array}\end{array}$$

where triangular brackets denote temporal averaging, u is the velocity normal to the transect (positive when directed into the estuary) and A(S) is the part of the transect area with instantaneous salinities >S. It should be noted that Q(S) is the streamfunction of the estuarine circulation in salinity space (MacCready, 2011; Burchard et al., 2025). When defining S_max and S_min as the maximum and minimum salinities occurring on the transect during the averaging period, respectively, then sufficiently long averaging results in $Q\left(S_{max}\right)=Q_{\mathrm{s}}\left(S_{max}\right)=Q_{S^{\mathrm{2}}}\left(S_{max}\right)=\mathrm{0}$, $Q\left(S_{min}\right)=-Q_{\mathrm{r}}$ (total volume transport equals river discharge) and $Q_{\mathrm{s}}\left(S_{min}\right)=\mathrm{0}$ (total salt transport vanishes under long-term averaging). The link to mixing is given by $Q_{{\mathrm{s}}^{\mathrm{2}}}\left(S_{min}\right)=\mathbb{M}$ (total salinity squared transport equals bulk mixing as defined in Eq. 20), see details in Burchard et al. (2019). These properties of Q(S), Q_s(S), and $Q_{{\mathrm{s}}^{\mathrm{2}}}\left(S\right)$ are demonstrated for a cross-channel transect near the mouth of the Elbe River estuary in Fig. 2a,c,d, where nearly balanced conditions are given such that the respective deviations from the expected values at $S=S_{min}$ are small.

Figure 2 TEF analysis using numerical model data from a cross-channel transect at along-channel position x_cux at Cuxhaven near the mouth of the Elbe River estuary (see Fig. C1a), averaged for the full month of April 2024. (a) Volume transport Q(S) across the transect, with the freshwater discharge Q_r and the dividing salinity s_div for reference. (b) Volume transport per salinity class, q(S), as well as the bulk inflow and outflow salinities s_in and s_out, respectively. The shaded areas and written numbers correspond to the bulk volume inflow (Q_in, red) and bulk volume outflow (Q_out, blue). (c) Salinity transport Q_s(S); (d) salinity-squared transport $Q_{{\mathrm{s}}^{\mathrm{2}}}$ across the transect, with the integrated mixing within the estuarine volume bounded by the same transect, 𝕄(x=x_cux), for reference.

Taking the S-derivative of Eq. (14) results in the volume, salinity and salinity-squared transport per salinity class (the Total Exchange Flow, TEF),

$$\begin{array}{}\text{(15)}& \begin{array}{rl}q\left(S\right)& =-{\displaystyle \frac{\partial Q\left(S\right)}{\partial S}},q_{\mathrm{s}}\left(S\right)=-{\displaystyle \frac{\partial Q_{\mathrm{s}}\left(S\right)}{\partial S}},\\ & q_{{\mathrm{s}}^{\mathrm{2}}}\left(S\right)=-{\displaystyle \frac{\partial Q_{{\mathrm{s}}^{\mathrm{2}}}\left(S\right)}{\partial S}},\end{array}\end{array}$$

where the minus sign ensures that inflow at high salinities is positive to highlight the character of the exchange flow. The connection to the Knudsen relations Eqs. (18), (21) and (22) is given by separately integrating the positive and negative contributions (denoted by superscripts ⁺ and ⁻) of the transport per salinity class,

$$\begin{array}{}\text{(16)}& \begin{array}{ll}{Q}_{\mathrm{in}}=\underset{S_{min}}{\overset{S_{max}}{\int }}{q}^{+}\left(S\right)\text{d}S,& Q_{\mathrm{out}}=\underset{S_{min}}{\overset{S_{max}}{\int }}{q}^{-}\left(S\right)\text{d}S,\\ & \\ Q_{\mathrm{s},\mathrm{in}}=\underset{S_{min}}{\overset{S_{max}}{\int }}{q}_{\mathrm{s}}^{+}\left(S\right)\text{d}S,& Q_{\mathrm{s},\mathrm{out}}=\underset{S_{min}}{\overset{S_{max}}{\int }}{q}_{\mathrm{s}}^{-}\left(S\right)\text{d}S,\\ & \\ Q_{{\mathrm{s}}^{\mathrm{2}},\mathrm{in}}=\underset{S_{min}}{\overset{S_{max}}{\int }}{q}_{{\mathrm{s}}^{\mathrm{2}}}^{+}\left(S\right)\text{d}S,& Q_{{\mathrm{s}}^{\mathrm{2}},\mathrm{out}}=\underset{S_{min}}{\overset{S_{max}}{\int }}{q}_{{\mathrm{s}}^{\mathrm{2}}}^{-}\left(S\right)\text{d}S,\end{array}\end{array}$$

and by deriving transport-weighted inflow and outflow salinities and squared salinities,

$$\begin{array}{}\text{(17)}& \begin{array}{rl}{s}_{\mathrm{in}}& ={\displaystyle \frac{Q_{\mathrm{s},\mathrm{in}}}{Q_{\mathrm{in}}}},s_{\mathrm{out}}={\displaystyle \frac{Q_{\mathrm{s},\mathrm{out}}}{Q_{\mathrm{out}}}},\\ (s^{\mathrm{2}}{)}_{\mathrm{in}}& ={\displaystyle \frac{Q_{{\mathrm{s}}^{\mathrm{2}},\mathrm{in}}}{Q_{\mathrm{in}}}},(s^{\mathrm{2}}{)}_{\mathrm{out}}={\displaystyle \frac{Q_{{\mathrm{s}}^{\mathrm{2}},\mathrm{out}}}{Q_{\mathrm{out}}}}.\end{array}\end{array}$$

An exemplary TEF profile is found in Fig. 2b for the Elbe River estuary. Zero values for q(S) occur for extreme values of Q(S), in consistency with Eq. (15). For the two-layer exchange flow shown here, there is one unique maximum of Q(S), such that the salinity S at which this occurs is the dividing salinity s_div between inflow and outflow (MacCready et al., 2018). Note that the dividing salinity s_div can also be used to calculate the Knudsen parameters Q_in, Q_out, s_in, and s_out, providing a numerically more robust method compared to the direct computation via Eq. (16) (Lorenz et al., 2019). For multi-layer flows, multiple dividing salinities may occur (Lorenz et al., 2019; Burchard et al., 2025).

The TEF analysis framework has been applied for a variety of estuarine studies, such as for tidal estuaries (MacCready, 2011; Chen et al., 2012; Wang et al., 2017; Conroy et al., 2020; Lemagie et al., 2022; Reese et al., 2024), tidal bays (Gräwe et al., 2016; Rayson et al., 2017; Xiong et al., 2021; Lemagie et al., 2022), non-tidal estuaries (Lange et al., 2020; Burchard et al., 2025), inverse estuaries (Lorenz et al., 2020), fjords (Sutherland et al., 2011; Lemagie et al., 2022; MacCready and Geyer, 2024), and regional seas (Döös et al., 2004; Burchard et al., 2018a).

3.3 Knudsen theory

To assess how entire estuaries quantitively act as mixing machines, the local relations derived in Sect. 2.2 are now integrated over estuarine volumes. For this, the continuity equation Eq. (3) and the salinity equation Eq. (4) are first integrated over the estuarine volume V which is separated from the ocean by means of a fixed vertical transect:

$$\begin{array}{}\text{(18)}& Q_{\mathrm{in}}+Q_{\mathrm{out}}+Q_{\mathrm{r}}={\displaystyle \frac{\text{d}V}{\text{d}t}},Q_{\mathrm{in}}{s}_{\mathrm{in}}+Q_{\mathrm{out}}{s}_{\mathrm{out}}={\displaystyle \frac{\text{d}\left(\overline{s}V\right)}{\text{d}t}},\end{array}$$

(see Fig. 3a) with the inflow transport and representative salinity Q_in≥0 and s_in (the ocean water inflow), outflow transport and representative salinity Q_out≤0 and s_out (the brackish water outflow), as defined in Eqs. (16) and (17). Further quantities in Eq. (18) are the river run-off Q_r≥0 (assuming zero river salinity for simplicity), the average salinity in the estuary $\overline{s}$, and the volume-integrated salinity $\overline{s}V$. Details of the derivation of Eq. (18) are given in Burchard et al. (2019). The conservation laws Eq. (18) have already been formulated by Knudsen (1900) and have become the basis for analyses of exchange flow in many estuaries (see e.g., Ji et al., 2007; MacCready, 2011; Sutherland et al., 2011; Chen et al., 2012; Burchard et al., 2018a). The positioning of the transect that separates the estuarine volume from the ocean is arbitrary, but often the geographical location of the river mouth is chosen. In Eq. (18) freshwater transports across the surface (evaporation or precipitation), the bottom (submarine groundwater discharge) as well as horizontal diffusive salt transports across the transect are neglected. Relations including freshwater transport through the sea surface can be found in Lorenz et al. (2021). Under long-term averaged conditions, the volume and salt storage terms on the right-hand side of Eq. (18) would vanish. Under such circumstances, s_out≤s_in and $Q_{\mathrm{in}}\le -Q_{\mathrm{out}}$ must hold, as illustrated for the Elbe estuary in Fig. 4a, b, where balanced conditions with a nearly vanishing volume storage $\frac{\text{d}V}{\text{d}t}$ are given. Under such circumstances, inflow from the ocean occurs at higher salinities than outflow towards the ocean due to mixing with riverine water.

Figure 3 Sketch showing the principles of volume and salt conservation as well as mixing in estuaries: (a) Estuarine volume (light blue shading) bounded by a fixed transect (bold red line), showing the classical Knudsen (1900) transports and salinities as well as the Knudsen mixing law Eq. (21) as derived by MacCready et al. (2018). (b) Estuarine volume (light blue shading) bounded by an isohaline of salinity S (bold blue line), showing the diahaline advective (Q_dia) and diffusive transports (J_S) as well as the mixing law Eq. (32) as derived by Burchard (2020). (c) Envelope of a discrete estuarine sub-volume (light blue shading) around the isohaline of salinity S (dashed blue line), bounded by the isohalines of salinities $S+\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\Delta }S$ and $S-\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\Delta }S$ (bold blue lines). The mixing per discrete salinity class is shown as well as the limit for S→0 which results in the universal law of estuarine mixing Eq. (33) for an infinitesimally thin salinity class. The bottom is marked by an orange line, the surface by a grey line and the fixed river transect by a red line. Advective volume transports are marked by blue arrows and the diffusive salt transport is marked by a yellow arrow.

The bulk mixing of an estuary is then obtained by averaging the integrated salinity variance budget Eq. (B2) in time:

$$\begin{array}{}\text{(19)}& \begin{array}{rl}{\displaystyle \frac{\text{d}}{\text{d}t}}\underset{V}{\int }{s^{\prime }}_{\mathrm{tot}}^{\mathrm{2}}\text{d}V& \approx Q_{\mathrm{r}}{\overline{s}}_{\mathrm{tot}}^{\mathrm{2}}+Q_{\mathrm{in}}({\overline{s}}_{\mathrm{tot}}-s_{\mathrm{in}}{)}^{\mathrm{2}}\\ & +Q_{\mathrm{out}}({\overline{s}}_{\mathrm{tot}}-s_{\mathrm{out}}{)}^{\mathrm{2}}-\mathbb{M},\end{array}\end{array}$$

with the estuarine bulk mixing

$$\begin{array}{}\text{(20)}& \mathbb{M}=〈\underset{V}{\int }{\mathit{\chi }}_{\mathrm{s}}\text{d}V〉,\end{array}$$

see the derivations by MacCready et al. (2018) and Burchard et al. (2019). In Eq. (19), the approximations $(s^{\mathrm{2}}{)}_{\mathrm{in}}\approx {\left(s_{\mathrm{in}}\right)}^{\mathrm{2}}$ and $(s^{\mathrm{2}}{)}_{\mathrm{out}}\approx {\left(s_{\mathrm{out}}\right)}^{\mathrm{2}}$ have been made for simplicity, where (s²)_in and (s²)_out are inflowing and outflowing salinity squares, respectively, as defined in Eq. (17). Assuming long-term averaging such that the temporal derivatives vanish and using Eq. (18), we finally obtain

$$\begin{array}{}\text{(21)}& \mathbb{M}\approx Q_{\mathrm{in}}{s}_{\mathrm{in}}^{\mathrm{2}}+Q_{\mathrm{out}}{s}_{\mathrm{out}}^{\mathrm{2}}=s_{\mathrm{in}}{s}_{\mathrm{out}}{Q}_{\mathrm{r}}=(s_{\mathrm{in}}-s_{\mathrm{out}})s_{\mathrm{in}}{Q}_{\mathrm{in}},\end{array}$$

(MacCready et al., 2018), relating the Knudsen parameters directly to estuarine mixing. While Knudsen (1900) had mentioned the role of mixing for the estuarine exchange flow qualitatively, Eq. (21) gave the first quantitative estimate of estuarine mixing as a function of the Knudsen parameters. An accurate bulk mixing estimate allowing $(s^{\mathrm{2}}{)}_{\mathrm{in}}\cancel{=}{\left(s_{\mathrm{in}}\right)}^{\mathrm{2}}$ and $(s^{\mathrm{2}}{)}_{\mathrm{out}}\cancel{=}{\left(s_{\mathrm{out}}\right)}^{\mathrm{2}}$ has been derived by Burchard et al. (2019):

$$\begin{array}{}\text{(22)}& \mathbb{M}={\displaystyle \frac{s_{\mathrm{out}}(s^{\mathrm{2}}{)}_{\mathrm{in}}-s_{\mathrm{in}}(s^{\mathrm{2}}{)}_{\mathrm{out}}}{s_{\mathrm{in}}-s_{\mathrm{out}}}}{Q}_{\mathrm{r}}.\end{array}$$

For the special case of $(s^{\mathrm{2}}{)}_{\mathrm{in}}={\left(s_{\mathrm{in}}\right)}^{\mathrm{2}}$ and $(s^{\mathrm{2}}{)}_{\mathrm{out}}={\left(s_{\mathrm{out}}\right)}^{\mathrm{2}}$, Eq. (21) is identical to Eq. (22). These equalities are only exact when the inflowing and outflowing salinities are constant in time and space during the averaging interval (Burchard et al., 2019). For estuaries with strongly fluctuating salinities at the mouth (such as for the short Merrimack estuary, see Chen et al., 2012) relation Eq. (22) has to be used instead of Eq. (21) to obtain an accurate estimate for the mixing.

Relation Eq. (21) is demonstrated in Fig. 4c for a numerical simulation of the Elbe estuary. It can be seen that the estimate is near-exact for almost the entire length of the estuary, proving its value for the study of bulk mixing in realistic estuaries, as also tested in studies by Broatch and MacCready (2022) and Reese et al. (2024).

Figure 4 Bulk parameters of the exchange flow through cross-channel transects at each along-channel position x from a numerical simulation of the Elbe River estuary, averaged for the month of April 2024. (a) Volume inflow and outflow Q_in and Q_out, respectively, compared to the freshwater discharge Q_r. (b) Inflow and outflow salinities s_in and s_out, respectively, as well as their difference $\mathrm{\Delta }s=s_{\mathrm{in}}-s_{\mathrm{out}}$. (c) Integrated mixing 𝕄(x) within the estuarine volume bounded by a cross-channel transect at along-channel position x, split into the contributions of the physical mixing 𝕄_phy due to the mixing parameterisation as well as the numerical mixing 𝕄_num due to discretization errors. The directly computed mixing (solid lines) is compared to the mixing estimate Eq. (21).

The principle of salt mixing inside an estuary bounded by a fixed transect is sketched in Fig. 3a. The first relation of Eq. (21) shows that the mixing does also balance the exchange of squared salinity with the ocean, such that mixing can also be defined as the reduction of squared salinity integrated over the estuary (which often simplifies the calculations, see Burchard et al., 2019), as it is expressed locally in Eq. (B3). The second relation of Eq. (21) shows that estuarine mixing can be estimated simply by knowing inflowing and outflowing salinities and the river run-off (MacCready et al., 2018). The third relation of Eq. (21) demonstrates the relation between estuarine circulation (quantified as strength of Q_in, see Broatch and MacCready, 2022; MacCready and Geyer, 2024) and mixing, a topic that is expanded on in Sect. 3.6.

Using the Knudsen relations Eq. (18), yet another useful reformulation of Eq. (21) has been derived by Qu et al. (2022) for estuaries in which the riverine inflow has non-zero salinity:

$$\begin{array}{}\text{(23)}& \mathbb{M}\approx Q_{\mathrm{r}}(s_{\mathrm{r}}-s_{\mathrm{out}}{)}^{\mathrm{2}}+Q_{\mathrm{in}}(s_{\mathrm{in}}-s_{\mathrm{out}}{)}^{\mathrm{2}},\end{array}$$

which is identical to Eq. (21) for a river salinity of s_r=0. In Eq. (23), the right-hand side is split into two terms representing the mixing pathways from the inflows to the outflows, with the first one leading from the river inflow to the brackish water outflow and the second one leading from the seawater inflow to the brackish water outflow (with source-water salinities put first in the brackets).

The Knudsen mixing relation Eq. (21) has been extended by Lorenz et al. (2021) for the case of non-zero freshwater fluxes through the surface, i.e., precipitation and evaporation:

$$\begin{array}{}\text{(24)}& \mathbb{M}\approx s_{\mathrm{in}}{s}_{\mathrm{out}}(Q_{\mathrm{r}}+Q_{\mathrm{surf}})-s_{\mathrm{surf}}^{\mathrm{2}}{Q}_{\mathrm{surf}},\end{array}$$

with the surface freshwater transport Q_surf (positive for net precipitation) and the representative surface salinity s_surf (square root of surface salinity variance transport divided by −Q_surf). Since $s_{\mathrm{surf}}>(s_{\mathrm{in}}{s}_{\mathrm{out}}{)}^{\mathrm{1}/\mathrm{2}}$ for evaporation and $s_{\mathrm{surf}}<(s_{\mathrm{in}}{s}_{\mathrm{out}}{)}^{\mathrm{1}/\mathrm{2}}$ for precipitation, both evaporation and precipitation are sources of mixing, in addition to the exchange flow. The example of the Persian Gulf as an inverse estuary with strong evaporation is briefly discussed in Sect. 4.1.3.

For the interpretation of the mixing relations it is instructive to consider the mixing completeness Mc (MacCready et al., 2018; Burchard et al., 2019) by non-dimensionalising Eq. (21) using Q_r and s_in:

$$\begin{array}{}\text{(25)}& \mathrm{Mc}={\displaystyle \frac{\mathbb{M}}{s_{\mathrm{in}}^{\mathrm{2}}{Q}_{\mathrm{r}}}}\approx {\displaystyle \frac{s_{\mathrm{out}}}{s_{\mathrm{in}}}},\end{array}$$

where the river run-off Q_r and inflowing salinity s_in (sometimes equated to the ocean salinity) can be considered as the external forcing of the estuary. Mixing completeness in estuaries can cover the full range of theoretically possible values of $\mathrm{0}\le \mathrm{Mc}\le \mathrm{1}$. Table 1 gives a number of examples for estuarine systems with low, medium and high mixing completeness. It should be noted that the mixing completeness is always calculated with respect to a fixed transect and that for each estuarine system the mixing completeness varies strongly with discharge and tidal intensity (e.g., during the spring-neap cycle).

Figure 5 Sketch showing estuarine conditions for low mixing causing weak estuarine circulation (a) and high mixing causing strong estuarine circulation (b). In both panels, Q_r and s_in are supposed to be identical. With prescribed low mixing 𝕄 in panel (a) and a high mixing 𝕄 in panel (b), the Knudsen relations Eqs. (18) and (21) quantify Q_in, Q_out and s_out.

In the extreme case of no mixing, the riverine freshwater would flow out at the surface with no modification and no ocean water entering the system (s_out=0 and Q_in=0), such that the mixing completeness would be zero. In estuaries with low mixing, brackish water of only low salinity is produced, with s_out≪s_in and Q_in≪Q_r, such that the mixing completeness is $s_{\mathrm{out}}/s_{\mathrm{in}}\ll \mathrm{1}$, see sketch in Fig. 5a. This would theoretically be the case for deep fjords with low tidal energy (Inall and Gillibrand, 2010). Mostly, fjords do however have large water bodies and low discharge such that the freshwater is strongly diluted by tidal mixing as, for example, for the Puget Sound where the mixing completeness is as large as about 0.97 (see Sutherland et al., 2011, and Table 1). Low mixing values of about Mc=0.18 have been observed for the tidally intense Merrimack estuary during high discharge (see Chen et al., 2012). Under these conditions, the salt intrusion length shortens considerably such that high-salinity ocean water and low-salinity river water are in close contact at the mouth of this estuary discharging directly into the coastal ocean, leading to low mixing. Other low values of mixing completeness are also observed for the Hudson river estuary during neap tide (Mc=0.36, see Wang et al., 2017) and the Elbe River estuary at high discharge (Mc=0.37, see Reese et al., 2024). In both cases, the water is relatively strongly stratified in the region of the transect such that s_out≪s_in. The large variability that estuaries have due to changes in discharge and tidal intensity is demonstrated by the fact that the Elbe for low discharge, the Merrimack for low discharge, and the Hudson for spring show values of mixing completeness as high as Mc=0.75, Mc=0.86 and Mc=0.87, respectively. Almost complete mixing is obtained for the shallow Wadden Sea of the German Bight, as sketched in Fig. 5b. In the Sylt-Rømø-Bight, where the freshwater run-off is low and tidal mixing is high, Gräwe et al. (2016) calculated values of inflowing and outflowing salinity both >30 g kg⁻¹ with only about 1 g kg⁻¹ difference, such that the mixing compleness is about $\mathrm{Mc}=s_{\mathrm{out}}/s_{\mathrm{in}}\approx \mathrm{30}/\mathrm{31}\approx \mathrm{0.97}$. These are values comparable to Puget Sound, see above. For the Baltic Sea, a non-tidal semi-enclosed sea in northern Europe, Knudsen (1900) estimated s_in=17.4 g kg⁻¹ and s_out=8.7 g kg⁻¹, such that the mixing completeness is Mc=0.5, a value that could be confirmed by a multi-decadal model simulation (Burchard et al., 2018a, 2019).

Chen et al. (2012)Wang et al. (2017)Reese et al. (2024)MacCready (2011)Knudsen (1900)Burchard et al. (2018a, 2019)Reese et al. (2024)Chen et al. (2012)Wang et al. (2017)Sutherland et al. (2011)Gräwe et al. (2016)

Table 1 List of estuarine systems with typical values of mixing completeness Mc for different tidal and runoff conditions. Note that estuaries with strong temporal variation (e.g., the Hudson River estuary during the spring-neap cycle) are not in balance, such that Eq. (25) is only a rough approximation.

3.4 Water Mass Transformation and diahaline mixing

Often, it is instructive to consider dynamics of estuaries in an isohaline framework, i.e., to evaluate transports, mixing and other properties relative to moving surfaces of constant salinity (isohalines) instead of the Eulerian framework with fixed spatial coordinates. With this, a quasi-Lagrangean perspective is added to the analysis with reference to the moving flow. In the isohaline analysis, geographical features such as a fixed transect at the mouth of the estuary do not play a central role. Since isohalines can move inside and outside the estuarine water body and extend over large areas covering parts of the estuary and the river plume, the isohaline analysis treats estuary and river plume as a dynamic continuum. This isohaline view of estuarine dynamics was first proposed by Walin (1977), with specific reference to the Baltic Sea with its isohalines extending over up to 1000 km from the Central Baltic Sea to its Western reaches (Henell et al., 2023). Later, the isohaline concept was applied to tidal estuaries (MacCready and Geyer, 2001; MacCready et al., 2002; Wang et al., 2017) and river plumes (Hetland, 2005; Muche et al., 2026). Here, we first introduce a local diahaline analysis, before we discuss the bulk analysis of estuarine dynamics across isohaline surfaces in Sect. 3.5.

Local mixing can move isohaline surfaces vertically such that a diahaline mass transport occurs relative to the moving isohaline surface. When normalised to isohaline unit surface and unit mass, this results in a so-called entrainment velocity. Starting for explanation with the one-dimensional salinity budget equation Eq. (D3), a coordinate transformation from geopotential z to salinity coordinates s (assuming a stable salinity stratification with $\partial s/\partial z<\mathrm{0}$), after time averaging in salinity coordinates, 〈⋅〉_S, a formulation for the vertical entrainment velocity u_dia,z(S) (vertical velocity relative to the vertically moving isohaline) is obtained:

$$\begin{array}{}\text{(26)}& u_{\mathrm{dia},z}\left(S\right)={〈w-{\displaystyle \frac{\partial z}{\partial t}}〉}_S={\displaystyle \frac{\partial }{\partial S}}{〈K_{\mathrm{v}}{\displaystyle \frac{\partial S}{\partial z}}〉}_S=-{\displaystyle \frac{\partial j_{\mathrm{dia},z}\left(S\right)}{\partial S}}\end{array}$$

(Wang et al., 2017; Klingbeil and Henell, 2023), where j_dia,z(S) is the time-averaged upward salinity flux through the moving isohaline. The vertical velocity of the isohaline S due to both advection and turbulent diffusion is given by $\partial z/\partial t$ (with z being the vertical position of the isohaline). Details of the derivation of Eq. (26) are given in Eqs. (D1)–(D3). The meaning of Eq. (26) is sketched in Fig.  6a: there is a maximum of vertical turbulent salinity flux j_dia,z in the entrainment layer that caps the turbulent bottom boundary layer. This maximum results from a large vertical salinity gradient $|\partial s/\partial z|$ at a still high level of turbulence originating from the boundary layer, expressed as eddy diffusivity K_v. Below this maximum, vertical salinity flux is divergent, thus lowering the local salinity which in time-average leads to an upward entrainment velocity u_dia,z. Above the entrainment layer, the opposite happens, resulting in a downward salinity flux. A similar process has been described and sketched by Ferrari et al. (2016), using density fluxes near the bottom of the ocean. The exchange flow in the bottom boundary layer itself with upwelling near the bottom and downwelling above has already been described by Garrett (1991). It should be noted that the total diahaline salt flux consists of two contributions, with one advective contribution and one diffusive contribution:

$$\begin{array}{}\text{(27)}& \begin{array}{rl}{f}_{\mathrm{dia},z}\left(S\right)& =u_{\mathrm{dia},z}\left(S\right)S+j_{\mathrm{dia},z}\left(S\right)\\ & =-{\displaystyle \frac{\partial j_{\mathrm{dia},z}\left(S\right)}{\partial S}}S+j_{\mathrm{dia},z}\left(S\right),\end{array}\end{array}$$

with the consequence that volume flux and salt flux are not proportional to each other and that the distribution of the diffusive salt flux in salinity coordinates entirely determines the total diahaline salt flux.

To relate j_dia,z to mixing χ_s, Li et al. (2022) defined the local mixing per salinity class which for a vertical water column with a monotone salinity profile reads as

$$\begin{array}{}\text{(28)}& \begin{array}{rl}m\left(S\right)& ={〈-{\displaystyle \frac{\partial z}{\partial S}}{\mathit{\chi }}_{\mathrm{s}}〉}_S={〈-{\displaystyle \frac{\partial z}{\partial S}}\mathrm{2}{K}_{\mathrm{v}}{\left({\displaystyle \frac{\partial S}{\partial z}}\right)}^{\mathrm{2}}〉}_S\\ & =\mathrm{2}{〈-K_{\mathrm{v}}{\displaystyle \frac{\partial S}{\partial z}}〉}_S=\mathrm{2}{j}_{\mathrm{dia},z}\left(S\right).\end{array}\end{array}$$

which can be seen as a thickness-weighted time-average of the local mixing χ_s, see also Klingbeil et al. (2019), Burchard et al. (2021) and Li et al. (2022). Combining Eqs. (26) and (28) results in a key relation between entrainment velocity and mixing,

$$\begin{array}{}\text{(29)}& u_{\mathrm{dia},z}\left(S\right)=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\partial m\left(S\right)}{\partial S}},\end{array}$$

which could be called the diahaline mixing-entrainment relation. Note that here upward velocities u_dia,z and fluxes j_dia,z are denoted as positive quantities. Details of the derivation of Eq. (29) for non-monotone salinity distributions in three dimensions can be found in Klingbeil and Henell (2023). The principle of Eq. (29) is sketched in Fig. 6b: as given by Eq. (28), m has local maxima in the same locations as j_dia,z, i.e., in the entrainment layers. For mixing per salinity class increasing with height, $\partial m/\partial z>\mathrm{0}\iff \partial m/\partial S<\mathrm{0}$ (for stable salinity stratification), a positive entrainment velocity $u_{\mathrm{dia},z}>\mathrm{0}$ is expected. For mixing decreasing with height, the opposite occurs. This leads to a typical pattern of diahaline exchange flow in estuaries with positive (upward, towards lower salinities into the estuary) entrainment through an isohaline near the bottom, and a negative (downward, towards higher salinities out of the estuary) entrainment through the same isohaline near the surface further seawards. For realistic estuaries this has been shown for the Hudson River estuary (Wang et al., 2017), the Pearl River estuary (Li et al., 2022, 2024), the Elbe River estuary (Reese et al., 2024), the Changjiang River estuary (Chang et al., 2024) and the Baltic Sea (Henell et al., 2023). In particular, Henell et al. (2023) and Reese et al. (2024) calculated both sides of Eq. (29) independently to demonstrate their equality (aside from small numerical differences) in real-world estuarine systems. The advantage of Eq. (29) over Eq. (26) is given by the fact that χ_s and thus m can be split into physical and numerical contributions (see Sect. 5.2.2), such that numerically generated spurious entrainment can be calculated, as shown by Henell et al. (2023) for the Baltic Sea.

The relationship between diahaline mixing m and entrainment velocity u_dia,z across the isohaline of 11 g kg⁻¹ is shown in Fig. 7 for the Elbe River estuary in northern Germany. It is clearly visible that as stated in Eq. (29), entrainment requires mixing since hotspots of the two quantities align well. In the up-estuary reach of the isohaline surface, where it is close to the bottom, upwelling (red) dominates, whereas at the down-estuary near-surface reaches of the isohaline downwelling (blue) dominates.

Figure 6 Sketch demonstrating the mechanism and distribution of diahaline exchange flow in estuaries. (a) Generation of diahaline exchange flow by means of a divergent vertical turbulent salinity flux j_dia,z, according to Eq. (26), shown for the bottom boundary layer of an estuary. High values of j_dia,z are marked by a black whirl, and low values are marked by grey whirls. The entrainment velocity is marked as red (upwards, $u_{\mathrm{dia},z}>\mathrm{0}$) and blue (downwards, $u_{\mathrm{dia},z}<\mathrm{0}$) arrows. (b) Situation of time-averaged diahaline exchange in an estuary. The bottom boundary layer and the surface boundary layer are marked by colour, both being separated from a more stratified interior via entrainment layers. Three exemplary isohalines are drawn. The entrainment velocity is again marked as red and blue arrows, where the size of the arrows corresponds to its relative magnitude. On the right side of the sketch a typical profile of the local mixing per salinity class m(S) is shown, along with consistent signs of the entrainment velocity u_dia,z, according to Eq. (29).

Figure 7 Diahaline mixing m (a) and diahaline entrainment velocity u_dia,z (b) across the isohaline surface of 11 g kg⁻¹, averaged over two spring-neap cycles during April 2024 in the lower Elbe River estuary in Germany. The line in both panels shows the 10 m isobath.

3.5 Estuarine mixing in isohaline volumes

Local diahaline mixing as introduced in Sect. 3.4 can be expanded to estuarine volumes (Walin, 1977). The local relation Eq. (27) for the total diahaline salt flux f_dia,z(S) can be integrated over the entire isohaline surface to result in

$$\begin{array}{}\text{(30)}& F_{\mathrm{dia}}\left(S\right)=Q_{\mathrm{dia}}\left(S\right)S+J_{\mathrm{dia}}\left(S\right)=-{\displaystyle \frac{\partial J_{\mathrm{dia}}\left(S\right)}{\partial S}}S+J_{\mathrm{dia}}\left(S\right),\end{array}$$

where F_dia is the total salt transport, Q_dia<0 is the diahaline volume transport, and J_dia>0 is the diffusive salt transport across the isohaline surface (see Fig. 8a by Walin, 1977). If instead of a fixed transect T a moving isohaline of salinity S is considered as the seaward boundary of the estuary (see Fig.  3b), the volume and salt budget is of this form:

$$\begin{array}{}\text{(31)}& \begin{array}{rl}& Q_{\mathrm{dia},\mathrm{in}}\left(S\right)+Q_{\mathrm{dia},\mathrm{out}}\left(S\right)+Q_{\mathrm{r}}={\displaystyle \frac{\mathrm{d}V\left(S\right)}{\mathrm{d}t}},\\ & \left(Q_{\mathrm{dia},\mathrm{in}}\left(S\right)+Q_{\mathrm{dia},\mathrm{out}}\left(S\right)\right)S+J_{\mathrm{dia}}\left(S\right)={\displaystyle \frac{\mathrm{d}\left({\overline{s}}_{\mathrm{tot}}V\right(S\left)\right)}{\mathrm{d}t}},\end{array}\end{array}$$

with the isohaline volume V(S), the average salinity inside this volume ${\overline{s}}_{\mathrm{tot}}$, the net advective inflow through the isohaline, $Q_{\mathrm{dia},\mathrm{in}}>\mathrm{0}$, and the the net advective outflow through the isohaline, $Q_{\mathrm{dia},\mathrm{out}}<\mathrm{0}$. Transformations of Eq. (31) show that long-term averaged mixing inside the estuarine volume bounded by an isohaline S is

$$\begin{array}{}\text{(32)}& \mathbb{M}\left(S\right)=\underset{V\left(S\right)}{\int }{\mathit{\chi }}_{\mathrm{s}}\text{d}V=S^{\mathrm{2}}{Q}_{\mathrm{r}},\end{array}$$

which can be seen as a special case of 𝕄≈s_ins_outQ_r from Eq. (21) with $s_{\mathrm{in}}=s_{\mathrm{out}}=S$ (Burchard, 2020). The relation Eq. (32) is exact for long-term averaging and zero freshwater transports through the surface and bottom of the estuary. A mixing relation for non-zero river salinity is shown in Eq. (F27).

With 𝕄(S) being a continuous function of S and assuming that Q_r is independent of S, we can take the derivative of M(S) with respect to S:

$$\begin{array}{}\text{(33)}& \mathbb{m}\left(S\right)={\displaystyle \frac{\partial \mathbb{M}\left(S\right)}{\partial S}}=\mathrm{2}S{Q}_{\mathrm{r}},\end{array}$$

where 𝕞(S) is the salt mixing per salinity class. It should be noted that 𝕞(S) can also be obtained by integrating the local mixing per salinity class m(S) from Eq. (29) over the projection of the isohaline surface to the horizontal (Li et al., 2022). A discrete version of Eq. (33) is sketched in Fig. 3c in order to explain the infinitesimal property of the mixing per salinity class. The linear dependence of 𝕞(S) on salinity S has been formulated and derived as the universal law of estuarine mixing (Burchard, 2020).

The relation Eq. (33) can be explained by first stating that the volume transport across the isohaline must for long-term averaged conditions equal the river runoff Q_r. Furthermore, the advective salt transport across the isohaline, $-S(Q_{\mathrm{dia},\mathrm{in}}+Q_{\mathrm{dia},\mathrm{out}})=S{Q}_{\mathrm{r}}$, must equal the diahaline diffusive salt transport, such that J_dia(S)=SQ_r, see the second equation in Eq. (31). With

$$\begin{array}{}\text{(34)}& \mathbb{m}\left(S\right)=\mathrm{2}{J}_{\mathrm{dia}}\left(S\right),\end{array}$$

which can be derived by integration of Eq. (28) over the horizontal projection of the isohaline surface, relation Eq. (33) is obtained (Burchard et al., 2021). A relation equivalent to the combination of Eqs. (33) and (34) had been derived by Hetland (2005), based on earlier work of Garvine (1999), to quantify the turbulent salt transport into river plumes due to entrainment, see the steady-state version of his Eq. (4).

To accurately reproduce the universal law by means of models of realistic estuaries such as the Pearl River estuary (Li et al., 2022), the Changjiang River estuary (Chang et al., 2024) and the Elbe River estuary (Reese et al., 2024), averaging over one spring-neap cycle is typically sufficient. In a numerical model, the mixing which a salinity field experiences is the sum of the parameterised physical mixing 𝕞_phy(S) and spurious numerical mixing 𝕞_num(S) due to the discretisation of the advection operator, see Sect. 5.2.2 for details. Both, relations Eqs. (32) and (33), are tested for the Elbe estuary in Fig. 8. There, it can be seen that the directly computed total mixing quantities 𝕄(S) and 𝕞(S), each consisting of the sum of numerical and physical mixing, closely follow their respective relation up to the point where the isohaline surfaces partly leave the model domain through the open boundary (grey-shaded areas in the Figure). Here, we first find an underestimation of the predicted mixing by relations Eqs. (32) and (33) which can likely be related to left-over stratification in the German Bight from earlier high-discharge periods entering the model domain via the open boundary and then being mixed away. For very high salinity classes, the mixing in the model is much weaker than the predicted mixing since substantial parts of the isohaline surfaces are outside of the model domain such that most of the potential mixing is not covered by the model. While the Knudsen mixing law Eq. (21) is only valid when salinity fluctuations at the open boundary are limited, the universal law of estuarine mixing, Eqs. (32) and (33), is exact for all estuaries (without relevant freshwater fluxes across the sea surface).

One interesting consequence of the universal law of estuarine mixing can be seen by reformulating Eq. (33) as

$$\begin{array}{}\text{(35)}& \mathbb{m}\left(S\right)=\mathbb{v}\left(S\right){\overline{\mathit{\chi }}}_{\mathrm{s}}\left(S\right)=\mathrm{2}S{Q}_{\mathrm{r}},\end{array}$$

where 𝕧(S) is the volume per salinity class and ${\overline{\mathit{\chi }}}_{\mathrm{s}}\left(S\right)$ is the salinity mixing averaged over the salinity class S. Since for long-term averages $\mathbb{v}\left(S\right){\overline{\mathit{\chi }}}_{\mathrm{s}}\left(S\right)$ is fixed, Eq. (35) means that at low mixing rates ${\overline{\mathit{\chi }}}_{\mathrm{s}}\left(S\right)$ the volumes per salinity class should be large, with the consequence that at low mixing rates, isohaline spacing should be wide. When comparing for example the salinity fields for the Hudson River estuary for neap tide and for spring tide (Warner et al., 2005a), the isohaline spacing at springs is wider than at neaps. Assuming that storage terms do not play an essential role for these situations, the conclusion could be drawn that average mixing is smaller at springs than at neaps. This is consistent with the findings of Warner et al. (2020) who find that maximum mixing occurs during late neap tides, see also Sect. 4.1.2. This is also in line with the study of Garvine (1999) who showed that increased (background) diffusivity would reduce the area of a river plume (see also the river plume study by Li et al., 2024, who showed that additional mixing due to islands in the plume region can reduce the plume area and volume).

Figure 8 Salt mixing from a realistic numerical simulation of the Elbe River estuary, averaged for the full month of April 2024. (a) Integrated mixing 𝕄(S) within an estuarine volume bounded by an isohaline surface of salinity S as computed directly from numerical model data (solid black line) as well as from the freshwater discharge Q_r using equation Eq. (32) (dashed line). (b) Mixing per salinity class 𝕞(S) as computed directly from numerical model data (solid black line) as well as from the universal law of estuarine mixing, Eq. (33) (dashed line). In each panel, the respective contributions of the physical mixing 𝕄_phy and 𝕞_phy due to the mixing parameterisation as well as the numerical mixing 𝕄_num and 𝕞_num due to discretisation errors to the total diagnosed mixing are shown as solid grey lines.

3.6 Relating estuarine circulation to mixing

In his famous abyssal recipes Munk (1966) fitted a vertical one-dimensional advection-diffusion equation to hydrographic observations in the central Pacific Ocean and concluded that turbulent mixing with an effective vertical diffusivity of $\mathrm{1.3}\times {\mathrm{10}}^{-\mathrm{4}}$ m² s⁻¹ would be needed to explain the global overturning circulation. In later studies, the decomposition of the underlying mixing processes into wind and tidal mixing and their regional distribution had been further specified (see e.g., Munk and Wunsch, 1998; Kuhlbrodt et al., 2007; Nikurashin and Ferrari, 2011; Cessi, 2019). On the much smaller scales of estuaries, the same must be postulated: estuarine circulation requires mixing and vice versa. Here, we discuss different concepts of this duality.

Let us first summarise what we have discussed about this issue until now. In his fundamental paper, Knudsen (1900) already stated that estuarine circulation is associated with mixing (Sect. 1). The general process is that salty ocean water entering the estuary is first mixed with fresh river water inside the estuary and then ejected as brackish water towards the ocean, making estuaries mixing machines (MacCready and Banas, 2011; Wang et al., 2017). When there is no mixing inside the estuary, then no further salty water can enter the estuary in the long term (Sect. 3.3). In that sense, the volume transport of salty water flowing into the estuary, Q_in, is a good measure for the estuarine circulation (Broatch and MacCready, 2022; MacCready and Geyer, 2024). A first quantitative estimate for the tight relationship between estuarine circulation and mixing has been established in the third relation of Eq. (21) by MacCready et al. (2018), showing a proportionality between the bulk mixing 𝕄 and Q_in, with (s_in−s_out)s_in as factor of proportionality.

The streamfunction Q(S) of the estuarine circulation as defined in Eq. (14) is the time-averaged volume transport into the estuary across a fixed transect for all salinities above S and therefore contains the information about the Total Exchange Flow, see Sect. 3.2 for details. It can be directly linked to mixing, as derived already by Walin (1977) to quantify the overturning circulation of the Baltic Sea:

$$\begin{array}{}\text{(36)}& Q\left(S\right)=Q_{\mathrm{dia},\mathrm{est}}\left(S\right)=-{\displaystyle \frac{\partial J_{\mathrm{dia},\mathrm{est}}\left(S\right)}{\partial S}}=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\partial {\mathbb{m}}_{\mathrm{est}}\left(S\right)}{\partial S}},\end{array}$$

where the subscript _est means that diahaline transport Q_dia, diahaline salt flux J_dia and diahaline mixing per salinity class 𝕞 are only considered for the part of the isohaline that is located on the estuarine side of the transect, see Fig. 9. The first equality in Eq. (36) is demonstrated in Fig. 9a: under long-term averaged conditions the volume transport Q(S) into the subvolume bounded by the transect T, the isohaline S and the bottom must be equal to the volume transport across the isohaline S, Q_dia,est(S). The second equality in Eq. (36) results from the integration of the entrainment relation Eq. (26) over the projection of the isohaline part situated upstream of the transect T. The third relation of Eq. (36) which had not been proposed by Walin (1977) is simply the S-derivative of Eq. (34) restricted to the upstream part of the estuary. Walin (1977) stated about this relation (his Eq. 7): Equation (36) represents in the most simple way how the deep water supply is related to the overall vertical (i.e. cross-isohaline) mixing properties in the basin. What Walin (1977) specifically calls the deep water supply to the Central Baltic Sea is more generally the up-estuarine part of the estuarine circulation. Note that independently of Walin (1977), Wang et al. (2017) used the first two equalities of Eq. (36) to calculate exchange flow accumulated between two estuarine segments.

When choosing the dividing salinity S=s_div for relation Eq. (36) with Q(s_div)=Q_in (see Sect. 3.2), then the quantification of the estuarine circulation is directly linked to mixing:

$$\begin{array}{}\text{(37)}& Q_{\mathrm{in}}=-{\displaystyle \frac{\partial J_{\mathrm{dia},\mathrm{est}}\left(s_{\mathrm{div}}\right)}{\partial S}}=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\partial {\mathbb{m}}_{\mathrm{est}}\left(s_{\mathrm{div}}\right)}{\partial S}},\end{array}$$

see details in Burchard et al. (2025) and Fig. 9b. The significance of Eq. (37) is that it is a direct quantification of the estuarine circulation (defined as Q_in) by means of the (S-gradient of the) diahaline mixing. This relation is directly applicable to simple estuaries with a typical two-layer exchange flow, but has also been extended to multi-layer exchange flow (Burchard et al., 2025). For the Elbe River estuary, relation Eq. (37) is demonstrated in Fig. 10.

Figure 9 Sketch showing the relation between estuarine circulation and diahaline mixing. (a) Demonstration of the relation Eq. (36) by Walin (1977). (b) Demonstration of the relation Eq. (37) by Burchard et al. (2025) using the dividing salinity isohaline to quantify the exchange flow by means of Q_in.

Figure 10 The relation Eq. (37) between estuarine circulation (quantified as the inflowing volume transport Q_in) and estuarine mixing (quantified as 𝕞_est), averaged over two spring-neap cycles during April 2024 in the lower Elbe River estuary in Germany.

3.7 Mixing of other constituents than salt

This review focuses on mixing of salt. The reason is that salinity is the defining constituent of estuaries, continuously ranging from minimum values near zero in the river water towards ocean salinity values near the mouth of the estuary or in the river plume. Therefore, salinity can be used as a coordinate in estuaries (Walin, 1977) substituting spatial coordinates. Also, salinity is conservative with basically no inner sinks and sources, and also bottom and surface fluxes of salt are negligible. However, many other constituents are mixing in estuaries, such as heat, oxygen, nutrients, and pollutants. Based on the work of Walin (1977, salinity coordinates) and Walin (1982, temperature coordinates), for larger ocean scales, theoretical Water Mass Transformation (WMT) frameworks have been developed to analyse mixing of constituents other than salt (Hieronymus et al., 2014; Groeskamp et al., 2019). To evaluate non-conservative behaviour of estuarine tracers due to sources or sinks, such tracers are often represented as function of salinity (Boyle et al., 1974). By doing so, non-conservative tracer mixing is identified by a non-linear relation between tracer concentration and salinity. However, as shown by Loder and Reichard (1981), such non-linear behaviour could also be caused by tracer variability in the freshwater source of the estuary. There is a large body of literature about effects of estuarine mixing of tracers other than salt on ecosystem dynamics. For example, Geyer (1993) proposes differential vertical mixing of suspended particulate matter (SPM) as a mechanism of creating Estuarine Turbidity Maxima. Tidal covariance between longitudinal velocity and concentration of SPM due to vertical SPM mixing can lead to up-estuary SPM transport (Scully and Friedrichs, 2007). In a similar way, Scully et al. (2022) explain the generation of local maxima of carbon dioxide partial pressure in estuaries, the so-called Estuarine Gas Exchange Maxima. Nitrogen-to-phosphate ratios in estuaries has been shown to critically depend on tidal mixing (Lui and Chen, 2011). These are just a few examples which show the essential role of mixing of tracers other than salt in estuaries. However, a general theory for such tracer mixing has not yet been proposed.

4 Major mixing processes and estuarine mixing hotspots

In the previous chapters we have presented various local and bulk theories of mixing and showed examples for the Elbe River estuary. These theories and examples prove that mixing, defined as integrated or local salinity variance decay due to turbulent processes, is an ubiquitous element in estuaries. Moreover, mixing defines what an estuary, consisting of a mixture of ocean and river water, is. While we know from the bulk mixing rules for estuaries, e.g., the Knudsen mixing law Eq. (21) or the universal law of estuarine mixing Eq. (33), how strong the overall mixing is in an estuary, we need to understand where the mixing occurs in time and space and which processes drive it. The intensity of mixing in an estuary is dictated by Eq. (10), ${\mathit{\chi }}_{\mathrm{s}}\approx \mathrm{2}{K}_{\mathrm{v}}(\partial s/\partial z{)}^{\mathrm{2}}$, indicating that mixing depends linearly on the intensity of turbulence, as expressed by the vertical diffusivity, and quadratically on the strength of the vertical salinity gradient. Typically in estuaries, the stronger the turbulence, the weaker the vertical salinity gradient, so it is not obvious a priori where and when mixing will be maximal in an estuary. In this section, the mixing in the Elbe estuary (and in one case that of the James River estuary) is used to provide an example of the various factors influencing the temporal and spatial variation of mixing in a partially mixed estuary.

4.1 Temporal variability

In estuaries, various time scales are relevant, including semi-diurnal tidal time scales and the fortnightly spring-neap cycle as well as times scales of weather and river-run-off (days to months). In the subsequent sections, the most relevant processes on these time scales are discussed.

Figure 11 Dynamics of vertical salinity variance, spatially averaged over the Elbe River estuary domain between river kilometres 85 and 160, during four tidal cycles. (a) Vertical salinity variance; (b) Terms in the variance budget Eq. (12); (c) Straining term split into longitudinal and lateral components.

4.1.1 Tidal variability

An analysis of the tidal cycle of the vertical salinity variance in the Elbe River estuary using Eq. (12) averaged over most of the estuary demonstrates the sequence of processes driving mixing (Fig. 11).

Panel (a) of Fig. 11 indicates that the vertical variance shows considerable variation through the tidal cycle, sharply rising at the end of flood, then decreasing to a minimum near the end of ebb. Panel (b) of Fig. 11 shows the individual terms in the vertical variance balance Eq. (12). Lateral straining is strongest during the flood tide (Fig. 11c). It accounts for most of the increase in stratification (as expressed by total vertical variance) toward the end of the flood tide, but it has little correspondence with mixing. This contribution of lateral straining has been observed in other estuaries (Purkiani et al., 2015; Geyer et al., 2017), with the important consequence that it produces a maximum in stratification at the beginning of the ebb tide. The longitudinal strain is actually negative during the flood (i.e., weakening stratification), but it is strongly positive during the ebb. This is the well-known signal of tidal straining, first described by Simpson et al. (1990). Particularly notable are the mid-ebb peaks in mixing, which are almost exactly in phase with the peak in longitudinal straining. This correspondence between longitudinal straining during the ebb and estuarine mixing has been found in other partially mixed estuaries including the Hudson River estuary (Wang and Geyer, 2018; Warner et al., 2020) and also the more strongly stratified Changjiang River estuary (Li et al., 2018) and Connecticut River estuary (Holleman et al., 2016).

Why is longitudinal straining so effective at increasing estuarine mixing? Going back to Eq. (10), we see that the intensity of mixing depends more sensitively on stratification than on the turbulence itself, although both are necessary to generate mixing. The positive strain that occurs during the ebb provides a continual source of stratification, which roughly balances the destruction of stratification by mixing during a significant fraction of the ebb tide (the times when the longitudinal strain and mixing have equal magnitude in Fig. 11). Paradoxically, the increased stratification during the ebb actually has an inhibitory influence on turbulence, but this inhibition of turbulence causes an enhancement of the vertical shear during the ebb. Figure 12 shows representative vertical profiles of velocity u, salinity s, eddy diffusivity K_v and mixing χ_s during flood and ebb in the Elbe estuary. The strong mixing that occurs during the ebb is found in the stratified boundary layer, in which the eddy diffusivity is actually much weaker than its value during the flood tide. The key to the mixing is the persistence of stratification, which is maintained by the strong shear that strains the along-estuary density gradient. Even though turbulence is weakened by stratification, it is not completely suppressed, due to the turbulence production originating from the bottom stress. During the flood tide, the boundary layer produces virtually no mixing, due to the absence of stratification. The only significant mixing occurs in the pycnocline when the well-mixed highly turbulent bottom boundary layer is entraining into the stratified layer above, where the turbulence is much weaker than the boundary layer but the stratification is strong. This shows that maximum mixing does not occur at the maxima of eddy diffusivity or salinity stratification, but at locations where both overlap (see also Warner et al., 2020, who report similar results for the Hudson River estuary).

Figure 12 Simulated vertical profiles of (a) the along-channel current velocity u, (b) salinity s, (c) vertical eddy viscosity K_v, and (d) local salt mixing χ_s from a near-shoal location within the inner Elbe River estuary at along-channel position x=127 km for ebb (blue) and flood (red), respectively, during a neap tidal cycle. The data was averaged over five neighbouring grid cells, corresponding to an along-channel distance of 360 m. Temporally, the data was averaged for one hour around peak ebb and peak flood, respectively.

4.1.2 Spring-neap cycle

The spring-neap cycle of tidal amplitude variation results in a large variation in the intensity of mixing, as shown in a timeseries based on the numerical model of the Elbe estuary. As often observed in partially mixed estuaries, the stratification (as represented by vertical variance, Fig. 13b) shows a large variation over the spring-neap cycle, with a sharp peak in stratification each neap tide. Again we have the paradoxical result that the peak mixing occurs during neap tides (Fig. 13c), when turbulent intensity is the weakest. Returning to Eq. (10), the mixing depends on the square of the vertical salinity gradient but only linearly on the eddy diffusivity. According to estuarine theory, stratification varies roughly as $K_{\mathrm{v}}^{-\mathrm{2}}$ (MacCready and Geyer, 2010), so the increased stratification is a much more important contributor to mixing than K_v through the spring neap cycle.

The timeseries of longitudinal strain through the spring-neap cycle (Fig. 13c) shows that it has similar spring-neap variation as mixing. The strain is a key ingredient for mixing – without it the stratification would vanish and there would be nothing to mix. The strain increases during neap tides due to the increased stratification, which augments the horizontal strain term in Eq. (12) directly by the increase in $s_{\mathrm{v}}^{\prime }$, and indirectly by the stratification-induced reduction in eddy viscosity, which increases $u_{\mathrm{v}}^{\prime }$ (MacCready and Geyer, 2010).

Other estuaries show a different phase relationship between mixing and the spring-neap cycle. For example, in the Hudson River estuary the peak mixing occurs between neaps and springs (Wang and Geyer, 2018), and in the Changjiang River estuary outflow the peak mixing occurs during spring tides (Li et al., 2018). These variations in the timing of mixing are related to the relative strength of tidal forcing to the stratifying tendency of the estuarine circulation. This balance is represented by the Simpson number as defined in Eq. (13). In the Elbe River estuary, Si remains low for most of the spring-neap cycle, with strong stratification only occurring around the time of neap tides. The Hudson River estuary has higher values of Si, and the Changjiang higher still, leading to persistent stratification through the spring-neap cycle (Li et al., 2018). A closer look into the dynamics of the Changjiang River estuary reveals that most of the mixing occurs outside the estuary in the extensive river plume, due to the high river discharge (Li et al., 2018; Chang et al., 2024). In the sense of the universal law of estuarine mixing (32), this can also be formulated as follows: Since the estuary is relatively fresh during high-flow conditions, mixing inside the estuary is small. Therefore strong mixing must occur outside the estuary, i.e. in the river plume, to amount to 𝕄(S)=S²Q_r, where here S is a salinity separating the estuary – river plume continuum from the adjacent coastal ocean. In this high Si regime, the strength of the turbulence becomes the limiting factor controlling mixing, leading to the peak mixing during spring tides.

Although these studies did not investigate the processes of mixing in detail, neap tidal turbulence might not be sufficiently strong to entrain the turbulent bottom boundary layer into the region of the surface-attached buoyant plume and cause mixing. Therefore, substantial near-surface salinity stratification remains until spring tides reduce it by mixing. In general it could be hypothesised that in tidally energetic estuaries vertical salinity variance is mixed away immediately once it is generated by straining during neap tides, as in the Elbe, Hudson and Pearl River estuaries. In stratified estuaries, this mixing process is delayed until spring-tide turbulence can efficiently mix, as in the Changjiang River estuary.

Figure 13 Dynamics of vertical salinity variance, spatially averaged over the Elbe River estuary domain between river kilometres 85 and 160, during four selected spring-neap cycles, as well as forcing conditions. (a) Runoff and wind speed; (b) vertical salinity variance and longitudinal tidal velocity amplitude; (c) terms of the vertical salinity budget according to Eq. (12).

4.1.3 Variation with variance input and direct meteorological forcing

In estuaries, there are four boundaries through which salinity variance can be introduced: the river boundary (river discharge), the seaward boundary (salinity fluctuations at the mouth), the sea surface (evaporation and precipitation) and the bottom (groundwater discharge, which we do not further consider here).

The bulk mixing laws for estuaries show clearly that under balanced conditions, mixing should be proportional to the river runoff Q_r. This is obvious from the Knudsen mixing law in an estuarine volume bounded by a fixed transect, 𝕄=s_ins_outQ_r, (see Eq. (21) as proposed by MacCready et al., 2018) and for the universal law of estuarine mixing inside a volume bounded by an isohaline surface of salinity S, 𝕄(S)=S²Q_r, (see Eq. (32) as proposed by Burchard, 2020). Since these two theories are based on long-term averaging, time lags between changes in runoff and changes in mixing are expected due to the storage of volume, salt and salinity variance (Broatch and MacCready, 2022). The dependence of mixing on runoff has most impressively been shown in the study by Broatch and MacCready (2022) for the Puget Sound: When the runoff during summer becomes four times larger than the spring runoff, mixing increases roughly by the same factor, showing a much stronger signal than the spring-neap cycle. This is specifically supported by the universal law Eq. (32): Since the salinity at the mouth of the large Puget Sound might not change much, the volume-integrated mixing should be largely proportional to the river discharge. Similarly, the Elbe River estuary simulations show higher neap-tide mixing peaks during high runoff than during low runoff (Figs. 13a, c), which is consistent with the estuarine mixing laws Eqs. (21) and (32).

In estuaries salinity fluctuations at the open ocean boundary are typically small and fluctuating with the tidal flow. In addition, salinity might vary with the dynamics of wind-driven upwelling and downwelling. An extreme example of the latter is the essentially non-tidal and highly industrialised Warnow River estuary in the Western Baltic Sea where downwelling events can decrease the offshore salinity from 20 to 8 g kg⁻¹ within a few hours (Lange et al., 2020), leading to substantial salinity variance changes in the estuary and an inversion of estuarine circulation. This variance input has also strong impacts on the mixing in the estuary (Burchard et al., 2025).

Evaporation and precipitation should have an effect on estuarine mixing by affecting the variance input available for mixing, see the extended Knudsen mixing law Eq. (24). For classical freshwater-dominated tidal estuaries we are not aware of dedicated studies of this effect, although the good agreement between the simulated mixing (including precipitation and evaporation) and the estimate Eq. (32) (which neglects precipitation and evaporation) in studies of such estuaries by Li et al. (2022) and Reese et al. (2024) indicates that its impact may be negligible. For the Persian Gulf, a large inverse estuary with some freshwater inflow, Lorenz et al. (2021) applied Eq. (24) and estimated that evaporation caused about half of the variance input for mixing, with the other half generated by the freshwater runoff. This ratio between the two mixing contributions needs to be compared to the ratio of the freshwater transports due to evaporation (−Q_surf) and river discharge (Q_r), which is 10:1 for the Persian Gulf. This implies that indeed, mixing of variance input from evaporation should generally have the tendency to be small compared to mixing caused by river discharge.

The effect of wind forcing on estuarine mixing has not yet been a focus of dedicated studies. Wind forcing can generally have two competing effects: straining of the salinity field and mixing (Scully et al., 2005; Chen and Sanford, 2009). Down-estuary wind forcing has a similar effect to ebb tidal straining by shearing less dense brackish water over dense ocean water and suppressing turbulence (Sect. 4.1.1 and Geyer, 1997). When a certain threshold is exceeded, the mixing effect of wind forcing would win over the straining effect, such that strong down-estuary winds are expected to be destratifying. In contrast to that, up-estuary winds have the same effect as flood straining and are therefore always destratifying. In a correlation analysis, Broatch and MacCready (2022) showed that wind forcing has some effect on mixing in Puget Sound, which however is dominated by the river runoff forcing. For a small weakly tidal estuary Yin et al. (2025) described the effect of wind pumping (covariance between wind stress and flow velocity) as an effective mechanism of up-estuarine salt transport which eventually will lead to increased salt mixing. However, in the tidally more energetic Elbe River estuary no obvious influence of the wind forcing is visible (compare Fig. 13a and c), although this has not yet been investigated by means of a correlation analysis.

4.2 Spatial variations

If estuaries had a flat bottom, their dynamics could be largely explained by means of one-dimensional or two-dimensional models without lateral variations. However, most estuaries are characterised by one or more deep channels (often deepened due to dredging) in longitudinal direction which carry most of the tidal flow. Shoals at the sides and between the channels lead to a typical channel-shoal structure where the channel-shoal transition leads to dynamic processes crucial to estuarine circulation and mixing (Sect. 4.2.1). But also in longitudinal direction, estuaries are not smooth. Channels often show a strong along-channel variability, e.g., due to constrictions in width and depth leading to local fronts and enhanced mixing (Sect. 4.2.2).

4.2.1 Channel-shoal interaction

Since flood and ebb currents in the tidal channels are faster than over the shoals, a significant lateral velocity gradient evolves over the channel-shoal transition (Fig. 14a). During flood, this velocity shear in conjunction with the longitudinal salinity gradient leads to a faster increase of the salinity in the channel than on the shoal (differential advection, see Huzzey and Brubaker, 1988; Geyer et al., 2020), such that a lateral salinity and thus density and internal pressure gradient is generated (Fig. 14b). The pressure gradient drives a lateral exchange flow leading to the generation of vertical salinity stratification (Fig. 14c) which is then mixed due to bottom-generated turbulence over the channel-shoal transition (Fig.  14d). A similar situation occurs during ebb, when differential advection leads to lower salinities in the channel centre as compared to the shoals. This substantially increased mixing over the channel-shoal transition has been shown by means of numerical model simulations for the Hudson River estuary by Warner et al. (2020) and for the Elbe River estuary by Reese et al. (2024, 2026). The intensified mixing over the channel-shoal transition in the Elbe River estuary is also demonstrated in Fig. 7a. Observations in San Francisco Bay (Collignon and Stacey, 2013) and the James River estuary (Huguenard et al., 2015) show enhanced mixing over the channel-shoal transition as well.

It should be noted that the lateral shear over the channel-shoal transition during flood is one leg of a lateral circulation across estuaries that leads to strong axial flow convergence near the surface (Nunes and Simpson, 1985). It was later found that lateral straining is also an important mechanism feeding back into the longitudinal estuarine circulation and up-estuarine salt transport (Lerczak and Geyer, 2004; Burchard et al., 2011; Bo et al., 2024), which indirectly leads to increased salt mixing in estuaries.

Figure 14 Sketch for explaining the effect of the channel-shoal transition on estuarine mixing during flood. The colour scaling indicates salinity, with dark blue colours representing high salinities. (a) early flood, when for simplicity no lateral salinity gradients are assumed; (b) full flood, generation of lateral salinity and thus density gradients due to differential advection; (c) full flood, lateral exchange flow driven by lateral density gradients leading to vertical stratification; (d) vertical shear generates small-scale turbulence which in concert with vertical stratification leads to vertical mixing.

4.2.2 Mixing at fronts

Bathymetrical bumps and lateral constrictions in a tidal channel in conjunction with tidal flow and longitudinal salinity gradients can lead to frontogenesis in estuaries (Geyer and Ralston, 2015) as well as increased mixing. The principle of this process is sketched in Fig. 15 for the example of an ebb current: Assuming an initially vertically well-mixed estuary with equally-spaced vertical isohaline surfaces, a sheared ebb current over a bathymetrical bump and a zone of weak recirculation downstream of the bump (Fig. 15a), the isohalines will be strained differentially. Over flat bathymetry, the classical ebb straining described by Simpson et al. (1990) will occur, leading to increased mixing (Sect. 4.1.1, and left side of Fig. 15b). On the lee side of the bump, where the ebb flow is partially blocked near the bed, the down-estuary advection of the isohalines is reduced or inverted while it is increased near the surface. This situation leads to a near-bottom retention of saline water on the lee side of the bump and thus to increased salt stratification. At the same time, the recirculation velocity near the bed and thus the vertical shear on the lee side of the bump will increase due to the backward slope of the strongly stratified isohalines. Increased shear, in turn, leads to increased turbulence production which in concert with strong stratification finally leads to strongly increased mixing. In their model simulations of the Hudson River estuary, Geyer and Ralston (2015) show that the composite Froude number (Armi and Farmer, 1986) in the frontal zone becomes supercritical in a similar way as was observed offshore of the lift-off point of a river plume at the mouth of an estuary (Horner-Devine et al., 2015). Mixing hotspots related to bathymetric bumps in tidal channels are also visible in model simulations for the Hudson River estuary (Geyer and Ralston, 2015; Warner et al., 2020).

Figure 15 Sketch explaining the effect of along-channel bathymetric bumps on increased mixing during ebb. (a) The vertical lines represent isohalines, with light grey shading indicating lower and dark grey shading indicating greater salinity. During this phase of early ebb, the salinity field is still assumed to be vertically well-mixed. The arrows represent horizontal flow velocity, showing some weak recirculation in lee of the bathymetric bump; (b) due to strong shear and blocking of the flow, the salinity field is most strongly strained in lee of the bump. The combination of strong shear-generated turbulence and strong salinity stratification then leads to a hotspot of vertical salt mixing.

A realistic example of frontal mixing comes from a numerical model representation of the partially mixed James River estuary (Fig. 16). This estuary has two pronounced constrictions, one at km 18, and the other at the mouth at km 30. The expansions on the seaward (right) side of these constrictions result in steep upward tilt of the isopycnals due to the supercritical hydraulic response to the expansions (Geyer et al., 2017). These steeply sloping isopcnals are associated with strong horizontal and vertical salinity gradients, i.e., frontal conditions. At the time of maximum ebb, mixing (as quantified by ${\mathit{\chi }}_{\mathrm{s}}\approx \mathrm{2}{K}_{\mathrm{v}}{\left(\partial s/\partial z\right)}^{\mathrm{2}}$) shows a pronounced maximum along each of the frontal zones. The energy for mixing along these frontal zones does not come from the bottom boundary layer, but rather from the baroclinic shear associated with the steeply sloping pycnocline. That strong shear also maintains the strong stratification via longitudinal straining of the horizontal salinity gradient. An observational example of the same phenomenon is discussed in Sect. 5.1 and illustrated in Fig. 17.

Figure 16 Realistic numerical simulation of conditions along the thalweg of the James River estuary during maximum ebb of spring tide, under moderate river discharge conditions, illustrating frontal zones with intensified mixing. The vertical component of the local mixing rate, ${\mathit{\chi }}_{\mathrm{s}}\approx \mathrm{2}{K}_{\mathrm{v}}{\left(\partial s/\partial z\right)}^{\mathrm{2}}$, is shown as colour scale. The along-estuary salinity distribution is indicated by black contours (contour interval 1 g kg⁻¹).

5 Methods to quantify mixing

5.1 Observational methods

Direct field measurements of the mixing of salt in estuaries is impractical, due to the microscopic scales at which the dissipation of salinity variance occurs, see Sect. 2.1. The turbulent motions that drive mixing can be resolved however, and numerous field investigations have quantified mixing based on measurements of turbulent motions at scales from metres to centimetres, then using theoretical arguments to relate the observable characteristics of the turbulence with either the mixing χ_s or the diahaline turbulent salt flux j_dia,z. A common approach is to use the relationship between the eddy diffusion coefficient for mass K_v and turbulent kinetic energy dissipation rate ε,

$$\begin{array}{}\text{(38)}& K_{\mathrm{v}}={\displaystyle \frac{\mathrm{Rf}}{\mathrm{1}-\mathrm{Rf}}}{\displaystyle \frac{N^{\mathrm{2}}}{\mathit{\epsilon }}}\end{array}$$

as proposed by Osborn (1980), where Rf is the flux Richardson number and N is the buoyancy frequency. With an estimate of K_v and a local measure of the vertical salinity gradient $\partial s/\partial z$, the vertical turbulent salt flux j_dia,z and the salinity variance dissipation χ_s are readily estimated, see Eqs. (10) and (26). Peters and Bokhorst (2000, 2001) pioneered the use of a free-falling shear probe in an estuary for the purpose of quantifying mixing; since then this method has been followed in a variety of estuarine settings (Rippeth et al., 2001; Becherer et al., 2011; Ross et al., 2019; Huguenard et al., 2019; Reese et al., 2026). This method is challenging for measuring turbulence near the bottom, due to the risk of smashing the delicate shear probe into the bottom. An alternative method is the estimation of the dissipation rate ε based on the inertial-subrange velocity spectrum measured by a turbulence-resolving current meter (Trowbridge et al., 1999; Giddings et al., 2011). A related methodology pioneered by Gargett (1994), and applied by Stacey et al. (1999), Lu and Lueck (1999) and Giddings et al. (2011) among others is to use an Acoustic Doppler Current Profiler (ADCP) to obtain a direct measurement of the Reynolds stress $\left[\stackrel{\mathrm{̃}}{u}\stackrel{\mathrm{̃}}{w}\right]$, then in combination with a measure of the vertical shear to estimate eddy viscosity A_v, and then to infer the eddy diffusivity K_v and the mixing rate χ_s.

These methods are most effective in weak stratification conditions, when dissipation rates tend to be high and the turbulent length scale is readily resolved by the sensor. As stratification gets stronger, however, the scales of turbulent motions decrease, as scaled by the Ozmidov scale

$$\begin{array}{}\text{(39)}& L_{\mathrm{O}}={\left({\displaystyle \frac{\mathit{\epsilon }}{N^{\mathrm{3}}}}\right)}^{\mathrm{1}/\mathrm{2}}\end{array}$$

making it harder to obtain a reliable estimate of ε. More problematical is that the estimation of mixing χ_s depends on the square of the local salinity gradient, see Eq. (9), which itself is a challenging quantity to measure at the small vertical scales relevant to turbulence within a stratified environment.

Micro-conductivity sensors provide a means of resolving the salinity gradient and associated overturns at scales relevant to the characterization of turbulent motions (Peters and Bokhorst, 2001). The Thorpe overturn scale L_T (Thorpe, 1977) provides an alternative means of estimating the turbulent dissipation rate:

$$\begin{array}{}\text{(40)}& \mathit{\epsilon }\approx L_{\mathrm{T}}^{\mathrm{2}}{N}^{\mathrm{3}},\end{array}$$

as shown by Peters (1997) in his turbulence measurements in the Hudson River estuary. It is particularly useful in the stratified interior, where the small vertical scales of turbulence make other methods of estimating ε more difficult (Etemad-Shahidi and Imberger, 2002; McPherson et al., 2019).

All of the above methods depend on a turbulence closure assumption as well as an estimate of mixing efficiency to link the dissipation rate ε to the actual mixing of salt, χ_s. The use of micro-conductivity sensors offers a more direct approach to estimating the mixing of salt. Following Holleman et al. (2016), high-frequency micro-conductivity time series measurements can resolve the inertial subrange and the viscous-convective subrange of the salinity variance spectrum

$$\begin{array}{}\text{(41)}& S_{\mathrm{ss}}\left(k\right)=b_{\mathrm{0}}{\mathit{\chi }}_{\mathrm{s}}{\mathit{\epsilon }}^{-\mathrm{1}/\mathrm{3}}{\mathcal{K}}^{-\mathrm{1}}min(\mathcal{K},{\mathcal{K}}_{\mathit{\eta }}{)}^{-\mathrm{2}/\mathrm{3}},\end{array}$$

where S_ss is the power spectral density of salinity variance, b₀=0.4 is the Kolmogorov constant for scalar variance, 𝒦 is the wave number and ${\mathcal{K}}_{\mathit{\eta }}=\mathrm{0.04}(\mathit{\epsilon }/{\mathit{\nu }}^{\mathrm{3}}{)}^{\mathrm{1}/\mathrm{4}}$ is the Kolmogorov wave number, with the molecular viscosity ν. Note that the last term in Eq. (41) resolves the transition from the inertial subrange to the viscous-convective subrange at length-scales of ${\mathcal{K}}_{\mathit{\eta }}^{-\mathrm{1}}$. As long as the height of the spectrum within either the inertial or viscous-convective subrange can be estimated, Eq. (41) can be used in combination with an estimate of ε to estimate mixing, i.e., the dissipation of salinity variance χ_s, without relying on turbulence closure assumptions. Moreover, the formula depends only weakly on the estimate of ε, which is often hard to estimate accurately in the stratified turbulent regime. Holleman et al. (2016) used this method to estimate mixing rates in the highly stratified Connecticut River estuary. Their analysis demonstrated peak values of χ_s in the pycnocline in association with intense shear instability during ebb tide (Fig. 17).

Figure 17 Spatial distribution of χ_s (coloured dots) in the stratified shear layer of the early ebb in the Connecticut River estuary, estimated from combined observations by a string of Acoustic Doppler Velocimetres (ADVs) and microconductivity probes. The crosses are locations where values of χ_s could not be reliably estimated. The grey scales in the background indicate the intensity of acoustic backscatter. This figure has been taken from Holleman et al. (2016).

Acoustic backscatter also provides a more direct means of estimating mixing, as demonstrated by Lavery et al. (2013). Using a broadband array of echo sounders, Lavery et al. (2013) demonstrated that the observed scattering in the pycnocline of the Connecticut River estuary has a spectral slope consistent with the viscous-convective subrange. Other scatterers, such as fish, bubbles and sediment, have distinctly different spectral slopes. If the stratification is strong enough and other scatterers do not overwhelm the signal, the acoustic backscatter intensity can be used as a nearly direct measure of χ_s. While the signal-to-noise ratio of the acoustic amplitude does not offer the same precision as microstructure measurements, it produces remarkable spatial resolution, as demonstrated by acoustic measurements within a train of Kelvin-Helmholz instabililities in the Connecticut River estuary (Fig. 18). Based on the evidence provided by the analysis of Lavery et al. (2013), echo-sounding imagery can be used to identify regions of intense mixing, even if the actual magnitude of χ_s cannot be quantified.

Figure 18 Kelvin-Helmholtz instabilities in the Connecticut River estuary. The colour shading shows the decadal logarithm of the salinity mixing χ_s, estimated from acoustic backscatter measurements. This figure has been taken from Lavery et al. (2013).

5.2 Numerical modelling techniques

Since direct observations of turbulent properties in estuaries are very tedious and noisy (see the discussion in Sect. 2.1), the analysis of mixing in estuaries largely relies on numerical models. The advantage of numerical models is certainly their coverage of the entire four-dimensional estuarine space (three spatial directions and time), whereas observations can only sparsely cover this space. However, to ensure that the numerical model results sufficiently represent real estuaries, several measures need to be taken: besides realistic input data into the model (such as bathymetry, open boundary conditions, meteorological forcing) and a thorough validation using observational data, the numerical model itself must be physically sound and numerically accurate. There is an extensive body of literature addressing these two topics (Griffies, 2004; Umlauf and Burchard, 2005; Shchepetkin and McWilliams, 2005; Zhang et al., 2016; Klingbeil et al., 2018). Here, we will present in more detail two aspects which are key to the proper assessment of mixing in estuaries: turbulence closure modelling (Sect. 5.2.1) and numerical mixing analysis (Sect. 5.2.2).

5.2.1 Turbulence closure modelling

The starting point of turbulence closures are the fundamental laws of momentum, mass and energy conservation from which transport equations based on molecular viscosities and diffusivities can be derived, see e.g. the molecular salinity equation Eq. (2). Applying the Reynolds decomposition (see Lesieur, 2008, and the discussion in Sect. 2.1) leads to transport equations for Reynolds-averaged variables, see e.g. the salinity equation Eq. (4), which include unknown second moments such as the vertical turbulent salt flux $\left[\stackrel{\mathrm{̃}}{w}\stackrel{\mathrm{̃}}{s}\right]$. In a similar manner, exact transport equations can be derived for those second moments, which however would include unknown third moments, and so forth. This infinite series of unclosed higher and higher order equations establishes the turbulence closure problem. An example for a second-moment transport equation ($\left[{\stackrel{\mathrm{̃}}{s}}^{\mathrm{2}}\right]$ in this case) is Eq. (5) which includes among others the vertical turbulent transport of the micro-structure salinity variance $\left[\stackrel{\mathrm{̃}}{w}{\stackrel{\mathrm{̃}}{s}}^{\mathrm{2}}\right]$ as a third moment. Second-moment closures use parameterisations for all unknown third moments, such that the system of second-moment transport equations is closed. To substantially simplify the solution, equilibrium assumptions are made for most second moments in such a way that the sum of the transport divergence and the time change are set to zero. For the example of the micro-structure variance equation Eq. (5), this leads to the equality of stirring and mixing (P_s=χ_s), see also the discussion in Sect. 2.1.

A central element of turbulence closures is the eddy diffusivity assumption Eq. (7), relating turbulent tracer fluxes to Reynolds-averaged tracer gradients, with the eddy diffusivity as a factor of proportionality, leading to the principle of down-gradient turbulent tracer fluxes. Similar assumptions are made for momentum and turbulent quantities. The eddy diffusivities include the entire second-moment closure. For the key quantity of the turbulent kinetic energy ($k=\frac{\mathrm{1}}{\mathrm{2}}\left[{\stackrel{\mathrm{̃}}{u}}_j^{\mathrm{2}}\right]$), a budget equation is generally solved with shear production as source term and dissipation of TKE into heat as sink term. For stable stratification, turbulence leads to mixing as well as negative buoyancy production, which leads to an increase of local potential energy and acts as a further TKE sink term. For estuarine numerical modelling, the use of so-called two-equation turbulence closure models has become a good compromise between efficiency and accuracy (Warner et al., 2005a). In such models, the first equation generally is the budget-equation for the TKE, while the second equation is related to the length scale of turbulence, such as the dissipation rate of TKE, ε (the k-ε model), or the turbulence frequency, $\mathit{\omega }=\mathit{\epsilon }/k$ (the k-ω model). If properly calibrated, these different versions of length-scale related equations perform similarly (Warner et al., 2005b). The most important aspect of the calibration is to ensure the quantitatively correct damping of vertical turbulent mixing caused by stable stratification. The principles of this calibration process are explained in Appendix E. In essence, the steady-state gradient Richardson number is set to the value of Ri_st=0.25 such that for stronger stratification (Ri>Ri_st), turbulence is suppressed and for weaker stratification it is enhanced. A properly calibrated turbulence closure model does also reproduce the canonical values of mixing efficiency Γ≈0.2 and the steady-state flux Richardson number Rf_st≈0.15 (Osborn, 1980). It should be noted that these values would only be reached in numerical models for so-called stationary homogeneous shear layers where production and destruction of TKE are balanced. In cases of strong temporal variability or locations with a substantial vertical turbulent transport of TKE (such as in active entrainment layers), significant deviations can occur (see the discussion by Holleman et al., 2016).

5.2.2 Numerical mixing analysis

Although the numerical mixing in the model of the Elbe River estuary is relatively small (see Figs. 4 and 8), it can be of considerable size in other estuarine models (Ralston et al., 2017; Henell et al., 2023). Therefore, the quantification of numerical mixing is discussed here. As demonstrated in the previous sections, the comparison of analytically derived mixing relations with diagnosed mixing from numerical simulations requires the quantification of the total variance decay experienced by a tracer in the numerical model. This does not only consist of contributions from the parameterised turbulence closure (physical mixing), but also from the discretisation of the tracer advection operator (spurious numerical mixing; Griffies et al., 2000). It is assumed that the tracer advection discretisation is conservative, monotone (in the sense that it does not generate wiggles and new tracer maxima and minima) and weakly diffusive. Many advection schemes with these properties have been developed such as the FCT (Flux-Corrected Transport) schemes (Zalesak, 1979), TVD (Total Variation Diminishing) schemes (Pietrzak, 1998) and the MPDATA (Multidimensional Positive Definite Advection Transport Algorithm) schemes (Smolarkiewicz, 2006). They all use some degree of implicit diffusion to ensure monotonicity.

In many model applications, numerical mixing has been found to explain a large portion of the total (physical plus numerical) mixing. High values of numerical mixing of typically 50 % have been found for the Baltic Sea (Hofmeister et al., 2011) and the Puget Sound (Broatch and MacCready, 2022), while Li et al. (2018) report about 33 % of numerical mixing for their simulation of the Changjiang River estuary. Low numerical mixing has been seen for simulations with high explicit horizontal diffusivity (Reese et al., 2024, for the Elbe River estuary) or for idealised estuarine models (MacCready et al., 2018). To account for the role of numerical mixing and to evaluate measures for its reduction, numerical mixing analysis methods have been developed (Burchard and Rennau, 2008; Klingbeil et al., 2014; Schlichting et al., 2023; Banerjee et al., 2024). Here, we first briefly describe general methods to quantify physical and numerical mixing in ocean models, before we give recommendations on how to reduce numerical mixing (see below).

Methods to quantify numerical mixing

To demonstrate methods to discriminate between physical and numerical mixing and to accurately quantify their sum, we give a detailed one-dimensional example of an advection-diffusion equation in Appendix F, including its discretisation (Appendix F1), the derivation of the mixing analysis (Appendix F2), a stationary solution of the one-dimensional advection-diffusion equation (Appendix F3) and the derivation of its analytical solution (Appendix F4). For this analysis, a first-order upstream (FOU) method is used as an example because it allows for an explicit discretisation of the physical and numerical mixing. Since the FOU scheme for advection shown in Eq. (F3) is inherently diffusive, it is generally not used in ocean models. Instead, mostly non-linear schemes are used for which the spurious numerical variance decay cannot be analytically derived. To separate between physical and numerical mixing, it is convenient to carry out the advection and diffusion discretisation in different steps, as an operational-split method. After the advection step, the numerical mixing is diagnosed, and after the diffusion step, the physical mixing is calculated. It is also possible to further separate the numerical and physical mixing into horizontal and vertical contributions. Different methods of numerical mixing quantification have been proposed for the pure advection step. Since the advection equation for the squared tracer, Eq. (F2) for K_h=0, is equivalent to the advection equation for the tracer itself, Eq. (F1) for K_h=0, Burchard and Rennau (2008) proposed to additionally carry out an advection step for the squared tracer (which should conserve the squared tracer) and to subtract from it the square of the advected tracer, Eq. (F1), which reduces the squared tracer. This difference should be a good estimate for the variance reduction in a particular grid point during a particular time step. Division by the time step Δt should give the local numerical variance decay (or mixing). As an alternative, Klingbeil et al. (2014) proposed to calculate the left-hand side (time derivative and advection term) of the diagnostic tracer-square advection equation as an estimate for the right-hand side (which should be the numerical mixing). Yet another method has recently been proposed by Banerjee et al. (2024). All three methods are equivalent when integrated over a larger area, and differences do only show up in the local distribution of the numerical mixing. An accurate numerical mixing quantification requires the direct implementation of the analysis into the numerical model code, since analysing numerical mixing from model output even at high output frequency has proven quite inaccurate (Schlichting et al., 2023).

There are alternative diagnostic methods to quantify numerical mixing, however not in terms of variance decay (Gibson et al., 2017; Holmes et al., 2021; Drake et al., 2025).

Measures to reduce numerical mixing

Generally, a finer resolution should lead to a reduction of numerical mixing. However, Burchard and Rennau (2008) showed that this might not be very efficient, since for an idealised example a grid-refinement by a factor of nine in the horizontal direction and by a factor of four in vertical direction (equivalent to a 144-fold increase in computational resources) led to a reduction of numerical mixing by less than a factor of two. More promising is the better alignment of grid layers with isopycnals, or adding a Lagrangean type of vertical grid motion, to reduce vertical advection with respect to moving coordinate layers. Specifically, vertically adaptive coordinates proved to significantly reduce numerical mixing (Hofmeister et al., 2011; Gräwe et al., 2015). A further reduction of numerical mixing can be achieved by aligning the horizontal grid with the major tidal flow directions. The present Elbe River estuary model (adopted from Reese et al., 2026) as well as the Weser River estuary model by Rummel et al. (2025) align the curvilinear grid with the dredged navigational channel and not, as typically done, with the lateral shoreline. This reduces numerical mixing because in estuaries lateral salinity gradients are generally much larger than longitudinal ones, such that these specially constructed horizontal coordinates align with those gradients which then do not have to be advected with the along-estuary flow across grid interfaces.

Numerical mixing is a fact in all numerical model applications and cannot be avoided. It is therefore essential to quantify its contribution to the tracer distribution, which is the result of the sum of intended physical and unintended numerical mixing. Intentionally reducing physical mixing would be one way to obtain realistic total mixing. This has been impressively demonstrated by Ralston et al. (2017) who in a model application to the Connecticut River estuary reduced the physical mixing by reducing the steady-state Richardson number Ri_st from the canonical value of 0.25 (see Sect. 5.2.1) to a very low value of 0.1. By doing so, they increased the numerical mixing from 50 % to 66 %, but reduced the total mixing such that the resulting salinity structure was more realistic. Such measures should however be handled with care since the numerical mixing generally has a different spatial and temporal distribution than the physical mixing (Klingbeil et al., 2014; Henell et al., 2023).

6 Future Perspectives

As reviewed in this paper, the work of past decades nowadays provides a consistent theoretical framework for estuarine mixing, the foundations of which have been laid in the early work by Knudsen (1900) and Walin (1977). Although their work did not explicitly define and quantify mixing, mixing theories are conveniently founded on their frameworks. In agreement with turbulence theory (Osborn and Cox, 1972; Mellor and Yamada, 1974), local mixing of a certain tracer is defined as the dissipation rate of the local variance of this tracer, χ_s, see Eqs. (9) and (B1) and Burchard and Rennau (2008). Estuarine mixing 𝕄 itself is then the integral of χ_s over estuarine volumes averaged over a certain period of time (e.g., the spring-neap cycle). Relating these definitions to bulk forcing parameters for the estuaries leads to the Knudsen mixing law Eq. (21) when the estuary is bounded by a fixed transect (MacCready et al., 2018) or to the universal law of estuarine mixing Eq. (32) when the estuary is bounded by a moving isohaline (Burchard, 2020). Therefore, to understand the estuarine mixing 𝕄, its temporal and spatial composition by means of the local variance decay χ_s has to be studied. In a few cases, this was achieved through observations (e.g., Lavery et al., 2013; Holleman et al., 2016), but is more typically done with realistic numerical models of estuaries (e.g., Warner et al., 2020; Reese et al., 2024). These models are consistent with the mixing theories only if numerical mixing due to the discretisation of the tracer advection terms is included (e.g., Li et al., 2022). Using these observations and modelling techniques, the major mixing processes in some example estuaries are now understood, see e.g. Warner et al. (2020) for the Hudson River estuary, Broatch and MacCready (2022) for the Salish Sea, Henell et al. (2023) for the Baltic Sea and Reese et al. (2026) for the Elbe River estuary.

Due to its high relevance for the understanding of estuarine dynamics and its socio-economical and ecological consequences, the research on estuarine mixing will continue. More studies for other estuaries will certainly come, most probably resulting in a different weighting of the most relevant mixing processes, maybe even describing new processes. But apart from that, the future of estuarine mixing research will likely be dominated by the technological progress that allows for ever finer-resolved, more efficient numerical modelling down to smaller and smaller scales. Here, we discuss in detail where, and to which extent, we see potential for such progress.

6.1 Increased computational resources

It is expected that the trend of ever-increasing computational resources will continue. As an important development, the ongoing replacement of the more traditional Central Processing Units (CPUs) by modern Graphics Processing Units (GPUs) will improve the overall computational efficiency. As computational codes for coastal ocean models will become available for GPUs, this increase in computational power can be used for finer resolution models, longer simulation periods or larger model domains. Finer resolution models would be able to resolve further processes of estuarine mixing with more accuracy, such as smaller-scale topographic effects. Specifically, consequences of coastal engineering such as dredging and dam-building, which are commonly under-resolved in contemporary numerical simulations, could be reproduced more realistically. Today, multi-decadal simulations of well-resolved estuaries are hardly feasible. Such longer simulation periods could, however, be used to better reproduce interannual variability and consequences of long-term trends such as sea-level rise and changes in precipitation and evaporation patterns. Larger model domains could help to include more of the tidal or non-tidal parts of rivers or larger portions of the adjacent ocean including the river plume. It has been shown that the universal law of estuarine mixing Eq. (33) is only valid for isohaline surfaces that do not leave the model domain. On the other hand, isohaline theory does not distinguish between estuarine and river plume mixing, and often the transition between the two cannot be seen from model topography and coastlines. Therefore, it is desirable to simulate the entire river-estuary-plume region within a single model setup. Similarly, the confluence of multiple estuaries and river plumes could be simulated at fine resolution when computational resources further improve.

6.2 Large Eddy Simulation modelling

More powerful computer resources could also allow for models with improved model physics, such as using higher-order turbulence parameterisations or the application of Large-Eddy Simulation (LES) models. Both improvements would have a direct effect on the computation of mixing. One example for higher-order turbulence closure models could be the use of non-local models that do not enforce the down-gradient assumption of turbulent fluxes (see the recent study by Legay et al., 2025). Furthermore, parameterisations of Langmuir Turbulence which potentially affects estuarine mixing in the presence of wind waves could be added (Harcourt, 2015). The application of LES models to estuaries would mean that the most energetic turbulent eddies could be resolved instead of being parameterised by turbulence closure models. Then, only the small-scale mixing would need to be parameterised, for which generally relatively simple closures should be sufficient. A further advantage of LES models over classical coastal ocean models is their non-hydrostatic pressure calculation which would allow for the reproduction of non-hydrostatic effects such as internal lee wave generation (Skyllingstad and Wijesekera, 2004), generation of interfacial waves at pycnoclines when the surface and the bottom boundary layers interact (Yan et al., 2022) or Langmuir Circulation effects in shallow coastal waters (Wang et al., 2022). Two specific LES model applications have been reported by Li et al. (2008) and Li et al. (2010) who calculated a tidal water-column setup with periodic boundary conditions in both horizontal directions. So far, no further estuarine LES applications have been published. Specifically efficient LES model codes such as Oceananigans (Ramadhan et al., 2020) that can be executed on GPUs open vast possibilities of coastal (see the recent study by  Huang et al., 2025) and estuarine LES applications in the future, starting with idealised setups, but also with the future potential to simulate estuarine dynamics over more realistic topographies.

6.3 Water Mass Transformation theory

The concept of Water Mass Transformation (WMT) has been strongly extended in recent years (see, e.g., Hieronymus et al., 2014; Groeskamp et al., 2019). Such multi-dimensional WMT concepts using other constituents in addition to salinity, such as temperature or biogeochemical tracer concentrations, are typically applied to large-scale ocean problems. Estuarine applications have yet to be created, but could give insight where, for example, temperature plays an important role in addition to salinity, such as shown for the Arctic Ocean (Pemberton et al., 2015), the Persian Gulf (Lorenz et al., 2020) or the Baltic Sea (Henell et al., 2023).

6.4 Machine Learning

Finally, applications of Machine Learning (ML) have entered all fields of oceanography, including coastal and estuarine research. Typical estuarine applications would comprise the calculation of river discharge from meteorological data (Börgel et al., 2025) or the estimation of the salt intrusion length inside estuaries (Rummel et al., 2025). It is not yet obvious in which way fast, efficient and well-trained ML algorithms could be exploited to support research on estuarine mixing, but it can certainly be expected that ML applications will soon extend our toolkit in addition to numerical modelling and observational methods.

Appendix A: Stirring and mixing: the case of saltwater in a glass beaker

The effects of stirring and mixing in a fluid subject to molecular diffusion are explained here quantitatively by means of a simple analytical example. We assume a fluid in a glass beaker of 0.1 m depth and a cross-sectional area of 0.01 m², thus containing 1 L of fluid. Imagine our beaker initially contains salty water underlaying freshwater with a smooth transition in between, with a horizontally homogeneous distribution such that they only depend on the vertical coordinate z and time t. The vertical diffusive spreading of the salinity $\stackrel{\mathrm{ˇ}}{s}$ in our beaker would then be described by the classical one-dimensional diffusion equation of the form

$$\begin{array}{}\text{(A1)}& {\displaystyle \frac{\partial \stackrel{\mathrm{ˇ}}{s}}{\partial t}}=\mathit{\kappa }{\displaystyle \frac{{\partial }^{\mathrm{2}}\stackrel{\mathrm{ˇ}}{s}}{\partial z^{\mathrm{2}}}}\end{array}$$

with $\mathit{\kappa }={\mathrm{10}}^{-\mathrm{9}}$ m² s⁻¹ denoting the molecular diffusivity of salt. Let us assume an initial salinity distribution of the form

$$\begin{array}{}\text{(A2)}& \stackrel{\mathrm{ˇ}}{s}\left(z\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{s}_{max}\left(\mathrm{1}+\cos \left({\displaystyle \frac{n\mathit{\pi }z}{D}}\right)\right),\end{array}$$

with the maximum salinity $s_{max}=\mathrm{30}$ g kg⁻¹ where n is a positive integer, and D the thickness of the fluid inside the beaker. For n=1, the fluid is unstirred, and we only have a single layer of salty water underneath the freshwater with mixed conditions in between (see Fig. 1a). Increasing n may be viewed as a simple model for a stirring process that creates an increasing number of interfaces between salty and fresh water. For n=30 (Fig. 1b), the stirring process has created 30 such interfaces, and for n=90 (Fig. 1c), 90 interfaces can be defined. It should be clear that in a real beaker, the stirring process induces salinity patches and streaks that are highly distorted and three-dimensional. Nevertheless, our simple one-dimensional model is sufficient to illustrate the basic effects of stirring and mixing as shown in the following.

To investigate the temporal evolution of the salinity layers as a result of mixing, we insert an ansatz of the form

$$\begin{array}{}\text{(A3)}& \stackrel{\mathrm{ˇ}}{s}(z,t)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{s}_{max}\left(\mathrm{1}+a\left(t\right)\cos \left({\displaystyle \frac{n\mathit{\pi }z}{D}}\right)\right)\end{array}$$

into Eq. (A1), where a(t) denotes the dimensionless amplitude of the salinity with $a\left(\mathrm{0}\right)=\frac{\mathrm{1}}{\mathrm{2}}$. This yields a differential equation of the form

$$\begin{array}{}\text{(A4)}& {\displaystyle \frac{\mathrm{d}a\left(t\right)}{\mathrm{d}t}}=-{\displaystyle \frac{a\left(t\right)}{\mathit{\tau }}},\end{array}$$

where $\mathit{\tau }=D^{\mathrm{2}}/\left(\mathit{\kappa }{n}^{\mathrm{2}}{\mathit{\pi }}^{\mathrm{2}}\right)$ combines all parameters of the problem. The solution of Eq. (A4) is of the form

$$\begin{array}{}\text{(A5)}& a\left(t\right)=e^{-\frac{t}{\mathit{\tau }}},\end{array}$$

which reveals that τ plays the role of a mixing time scale.

The blue lines shown in Fig. 1 (see the figure caption for the parameters chosen) illustrate the behaviour of the solution found above for different values of n. For n=1 (no stirring, Fig. 1a) the effect of mixing is seen to be negligible after a period of ten minutes. If some gentle stirring is applied (figuratively using the spoon) such that n is increased to 30, considerable mixing effects can be seen already after 10 min (Fig. 1b). Only intensive stirring (n=90) leads to the desired result of (almost) complete mixing within a reasonable amount of time (Fig. 1c). For the parameters chosen, the mixing time scale is about τ=17 000 min for the case of no stirring (n=1), about τ=19 min for little stirring (n=30), and about τ=2 min for strong stirring (n=90).

Appendix B: Salinity-variance equations

This appendix introduces some Reynolds-averaged salinity-variance or salinity-square equations to demonstrate the function of the salinity variance decay χ_s as sink term. Based on Eq. (8) and using Eq. (9), the salinity variance equation with the local variance per unit volume ${s_{\mathrm{tot}}^{\prime }}^{\mathrm{2}}=(s-{\overline{s}}_{\mathrm{tot}}{)}^{\mathrm{2}}$ and the volume-averaged salinity ${\overline{s}}_{\mathrm{tot}}$ is of the following form:

$$\begin{array}{}\text{(B1)}& \begin{array}{l}\underset{\text{ change}}{\underbrace{\frac{\partial {s_{\mathrm{tot}}^{\prime }}^{\mathrm{2}}}{\partial t}+s_{\mathrm{tot}}^{\prime }\frac{\mathrm{d}}{\mathrm{d}t}{\overline{s}}_{\mathrm{tot}}}}+\underset{\text{advection}}{\underbrace{\frac{\partial \left(u{s_{\mathrm{tot}}^{\prime }}^{\mathrm{2}}\right)}{\partial x}+\frac{\partial \left(v{s_{\mathrm{tot}}^{\prime }}^{\mathrm{2}}\right)}{\partial y}+\frac{\partial \left(w{s_{\mathrm{tot}}^{\prime }}^{\mathrm{2}}\right)}{\partial z}}}\\ \underset{\text{diffusion}}{\underbrace{-\frac{\partial }{\partial x}\left(K_{\mathrm{h}}\frac{\partial {s_{\mathrm{tot}}^{\prime }}^{\mathrm{2}}}{\partial x}\right)-\frac{\partial }{\partial y}\left(K_{\mathrm{h}}\frac{\partial {s_{\mathrm{tot}}^{\prime }}^{\mathrm{2}}}{\partial y}\right)-\frac{\partial }{\partial z}\left(K_{\mathrm{v}}\frac{\partial {s_{\mathrm{tot}}^{\prime }}^{\mathrm{2}}}{\partial z}\right)}}\\ \\ =-\underset{\text{mixing}}{\underbrace{{\mathit{\chi }}_{\mathrm{s}}}},\end{array}\end{array}$$

where the advective and diffusive flux divergences conservatively re-distribute local variance. Mixing χ_s is the sink term for the local variance. An extra term due to the non-constant volume-averaged salinity ${\overline{s}}_{\mathrm{tot}}$ is included in the change term (MacCready et al., 2018, see Burchard et al. (2018b) for hints to the derivation). The variance budget of the entire estuary results from integration of Eq. (B1) over the volume of the estuary:

$$\begin{array}{}\text{(B2)}& \begin{array}{rl}& \underset{\text{change}}{\underbrace{{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\underset{V}{\int }{s_{\mathrm{tot}}^{\prime }}^{\mathrm{2}}\text{d}V}}+\underset{\text{boundary transport}}{\underbrace{\underset{A}{\int }\left(\mathit{u}{s_{\mathrm{tot}}^{\prime }}^{\mathrm{2}}-K_{\mathrm{h}}\mathbf{\nabla }{s_{\mathrm{tot}}^{\prime }}^{\mathrm{2}}\right)\cdot \text{d}\mathit{A}}}\\ & =-\underset{\text{mixing}}{\underbrace{\underset{V}{\int }{\mathit{\chi }}_{\mathrm{s}}\text{d}V}}\end{array}\end{array}$$

(MacCready et al., 2018; Burchard et al., 2019), where A is the boundary of the estuary towards the river and the ocean, u is the velocity vector at the open boundary and dA is the normal vector orthogonal to the area element A pointing out of the estuary, see details in Li et al. (2018). While the volume integrated mixing is the only sink term in Eq. (B2), an explicit source does not exist. Instead variance may enter the estuary via the boundary transport term as freshwater from the river and as saline water from the adjacent coastal ocean.

Inserting Eq. (7) into Eq. (6), or multiplying the salinity budget equation Eq. (8) by 2s, we obtain a budget equation for the squared salinity, which has the same mixing term as the local variance equation Eq. (B1):

$$\begin{array}{}\text{(B3)}& \begin{array}{rl}& \underset{\text{change}}{\underbrace{{\displaystyle \frac{\partial s^{\mathrm{2}}}{\partial t}}}}\underset{\text{advection}}{\underbrace{+{\displaystyle \frac{\partial \left(u{s}^{\mathrm{2}}\right)}{\partial x}}+{\displaystyle \frac{\partial \left(v{s}^{\mathrm{2}}\right)}{\partial y}}+{\displaystyle \frac{\partial \left(w{s}^{\mathrm{2}}\right)}{\partial z}}}}\\ & \underset{\text{diffusion}}{\underbrace{-{\displaystyle \frac{\partial }{\partial x}}\left(K_{\mathrm{h}}{\displaystyle \frac{\partial s^{\mathrm{2}}}{\partial x}}\right)-{\displaystyle \frac{\partial }{\partial y}}\left(K_{\mathrm{h}}{\displaystyle \frac{\partial s^{\mathrm{2}}}{\partial y}}\right)-{\displaystyle \frac{\partial }{\partial z}}\left(K_{\mathrm{v}}{\displaystyle \frac{\partial s^{\mathrm{2}}}{\partial z}}\right)}}\\ & =-\underset{\text{mixing}}{\underbrace{{\mathit{\chi }}_{\mathrm{s}}}},\end{array}\end{array}$$

which is the parameterised version of Eq. (6). Often, it is more handy to diagnose the budget of the squared salinity Eq. (B3) rather than the local variance Eq. (B1), since the time-variable volume-averaged salinity ${\overline{s}}_{\mathrm{tot}}$ does not need to be considered.

Appendix C: Elbe River estuary model

C1 Study region

The Elbe River estuary is located in northern Germany and extends over 150 km from a weir south of Hamburg to the German Bight in the North Sea. It is an M₂-dominated, mesotidal estuary at a mean tidal range of 3.0 to 3.5 m (Boehlich and Strotmann, 2008). The tidal signal is further modified throughout the spring-neap cycle. Its broad salinity range (0.5 to >30 g kg⁻¹ within the model domain), medium discharge intensity of about 800 m³ s⁻¹ on a multi-decadal average (Strotmann, 2017), and relatively simple, funnel-shaped, single-channel geometry make this estuary an ideal example for the illustration of basic estuarine mixing dynamics. Further, the navigational channel is subject to intensive dredging measures, comparable to other estuaries under high anthropogenic maintenance, and surrounded by extensive tidal flats.

Figure C1 Overview of the numerical simulation of the Elbe River estuary. (a) Downstream section of the setup topography, showing the open boundary in the German Bight as a bold grey line, the distance from the upstream boundary as yellow dots in km, and the location of the cross-channel transect used for the analysis in Fig. 2 as a red line. (b) Average surface salinity. White lines indicate even salinities. (c) Average salinity along the thalweg of the navigational channel. Solid (dashed) white lines indicate even (odd) salinities. (d) Estuarine circulation in terms of the tidally averaged along-channel velocity 〈u〉 along the thalweg of the navigational channel, where red shading indicates an inflow into the estuary and blue shading indicates an outflow towards the German Bight. In (b)–(d), the properties have been averaged for the full month of April 2024.

The Elbe is a partially mixed estuary with a medium stratification intensity as shown for the analysis period of April 2024 in Fig. C1b, c, where the spring-neap averaged surface to bottom salinity difference amounts to up to 6–7 g kg⁻¹, and the salt intrusion reaches more than 50 km up-estuary. It displays a clear two-layer exchange flow pattern (see Fig. C1d).

A more detailed description of the Elbe River estuary can be found in Burchard et al. (2004), Reese et al. (2024), and Burchard et al. (2025).

C2 Numerical model and setup

The three-dimensional numerical model data used for all Elbe River examples in this paper was created with the General Estuarine Transport Model (GETM; Burchard and Bolding, 2002; Klingbeil and Burchard, 2013). For turbulence closure, GETM is coupled to the General Ocean Turbulence Model (GOTM; Burchard et al., 1999; Li et al., 2021), here using a second-order, algebraic closure for a k-ε parameterization (Umlauf and Burchard, 2005).

The specific Elbe River estuary setup used here is identical to the setup presented in Reese et al. (2026) with a high horizontal (2404×397 cells with a mean resolution of 80 m in along-channel direction and 98 m in cross-channel direction) and vertical (40 adaptive layers) grid resolution. It covers the year 2024 and uses a curvilinear grid with thalweg-following and cross-channel grid lines, resulting in horizontal grid cells of variable size, where the finest resolution is achieved within the navigational channel. The model domain is limited in upstream direction by the location of a weir that marks the end of the upstream tidal intrusion, and in downstream direction by an open boundary within the German Bight.

Overall, the Elbe River setup uses a realistic forcing. This includes initial distributions for temperature and salinity as well as tidal elevations and vertical temperature and salinity profiles along the open boundary (E.U. Copernicus Marine Service, 2024), the daily-averaged freshwater discharge at the upstream end of the model domain (Wasserstraßen- und Schifffahrtsamt Magdeburg, 2024) at a constant salinity of 0.5 g kg⁻¹, and a 1 km-resolved meteorological forcing (Norwegian Meteorological Institute, 2024).

A full setup description along with a process-oriented model performance evaluation can be found in Reese et al. (2026).

Appendix D: Coordinate transformation of vertical salinity equation

Neglecting horizontal turbulent fluxes, the salinity equation Eq. (8) can also be formulated as

$$\begin{array}{}\text{(D1)}& {\displaystyle \frac{\mathrm{D}s}{\mathrm{D}t}}{\displaystyle }={\displaystyle \frac{\partial }{\partial z}}\left(K_{\mathrm{v}}{\displaystyle \frac{\partial s}{\partial z}}\right).\end{array}$$

Assuming a stable salinity stratification with $\partial s/\partial z<\mathrm{0}$, for every isohaline the vertical position can be given by $z=z^s(x,y,s,t)$ such that Eq. (D1) can be multiplied by $\partial z^s/\partial s$ and time-averaged on a fixed isohaline, 〈〉_s to yield

$$\begin{array}{}\text{(D2)}& {\displaystyle }\underset{=u_{\mathrm{dia},z}}{\underbrace{〈{\displaystyle \frac{\partial z^s}{\partial s}}{\displaystyle \frac{\mathrm{D}s}{\mathrm{D}t}}〉}}{\displaystyle }={\displaystyle \frac{\partial }{\partial s}}\underset{=-j_{\mathrm{dia},z}}{\underbrace{〈K_{\mathrm{v}}{\displaystyle \frac{\partial s}{\partial z}}〉}}\end{array}$$

and u_dia,z (the diahaline velocity per unit horizontal area) can be obtained from

$$\begin{array}{}\text{(D3)}& {\displaystyle }〈w〉=〈{\displaystyle \frac{\mathrm{D}z}{\mathrm{D}t}}〉{\displaystyle }=〈{\displaystyle \frac{\partial z^s}{\partial t}}+u{\displaystyle \frac{\partial z^s}{\partial x}}+v{\displaystyle \frac{\partial z^s}{\partial y}}〉+\underset{=u_{\mathrm{dia},z}}{\underbrace{〈{\displaystyle \frac{\partial z^s}{\partial s}}{\displaystyle \frac{\mathrm{D}s}{\mathrm{D}t}}〉}}\end{array}$$

Note that the diahaline fluxes u_dia,z and j_dia,z are defined positive upwards here, whereas they are usually defined positive up-gradient in the literature. For almost flat isohalines $u_{\mathrm{dia},z}\approx 〈w-\frac{\partial z^s}{\partial t}〉$ can be interpreted as the vertical velocity relative to the moving isohalines. Due to the neglect of horizontal turbulent fluxes, the diahaline diffusive salt flux per unit horizontal area j_dia,z in Eq. (D2) is identical to the vertical turbulent salt flux. Henell et al. (2023) and Reese et al. (2024) demonstrated that relation Eq. (D2) approximately holds also in realistic estuarine cases where horizontal turbulent fluxes are present. Moreover, relation Eq. (D2) also holds for arbitrary salinity distributions including inversions, based on generalised definitions of u_dia,z and j_dia,z (Klingbeil and Henell, 2023).

Appendix E: Calibration of two-equation turbulence closure models

After carrying out a second-moment closure as briefly described in Sect. 5.2.1 and presented in more detail by Burchard and Bolding (2001) and Umlauf and Burchard (2005) and applying the boundary layer assumption (vertical gradients are much larger than horizontal gradients), the closure results in the down-gradient parameterisation of momentum,

$$\begin{array}{}\text{(E1)}& \left[\stackrel{\mathrm{̃}}{u}\stackrel{\mathrm{̃}}{w}\right]=-A_{\mathrm{v}}{\displaystyle \frac{\partial u}{\partial z}},\left[\stackrel{\mathrm{̃}}{v}\stackrel{\mathrm{̃}}{w}\right]=-A_{\mathrm{v}}{\displaystyle \frac{\partial v}{\partial z}},\end{array}$$

with the eddy viscosity A_v, and the down-gradient parameterisation for turbulent tracer flux Eq. (7), with the following relation for the eddy viscosity A_v and the eddy diffusivity K_v:

$$\begin{array}{}\text{(E2)}& A_{\mathrm{v}}=c_{\mathit{\mu }}\left(Ri\right){\displaystyle \frac{k^{\mathrm{2}}}{\mathit{\epsilon }}},K_{\mathrm{v}}=c_{\mathit{\mu }}^{\prime }\left(Ri\right){\displaystyle \frac{k^{\mathrm{2}}}{\mathit{\epsilon }}},\end{array}$$

with the non-dimensional stability functions c_μ and $c_{\mathit{\mu }}^{\prime }$ depending in the case of quasi-equilibrium stability functions (Galperin et al., 1988) on the non-dimensional gradient-Richardson number Ri, the turbulent kinetic energy per unit mass, k, and its dissipation rate, ε. The definition of the gradient Richardson number

$$\begin{array}{}\text{(E3)}& Ri={\displaystyle \frac{N^{\mathrm{2}}}{S_{\mathrm{v}}^{\mathrm{2}}}}\end{array}$$

is given with the squared shear frequency $N^{\mathrm{2}}=\partial b/\partial z$, the squared vertical shear $S_{\mathrm{v}}^{\mathrm{2}}=(\partial u/\partial z{)}^{\mathrm{2}}+(\partial v/\partial z{)}^{\mathrm{2}}$. The stability functions c_μ and $c_{\mathit{\mu }}^{\prime }$ include the entire second-moment closure (Burchard and Bolding, 2001).

For the case of two-equation turbulence closure models the full transport equations for k and ε are calculated (Rodi, 1987). Instead of the dissipation rate equation, other length-scale related properties such as the turbulence frequency $\mathit{\omega }=\mathit{\epsilon }/k$ could be modelled, following the generic length scale (GLS) approach by Umlauf and Burchard (2003). The k-ε model, as it is typically used in three-dimensional coastal ocean models, is of the following form:

$$\begin{array}{}\text{(E4)}& \begin{array}{l}\frac{\partial k}{\partial t}+\frac{\partial \left(uk\right)}{\partial x}+\frac{\partial \left(vk\right)}{\partial y}+\frac{\partial \left(wk\right)}{\partial z}-\frac{\partial }{\partial z}\left(\frac{A_{\mathrm{v}}}{{\mathit{\sigma }}_k}\frac{\partial k}{\partial z}\right)\\ \\ =P+B-\mathit{\epsilon },\\ \\ \frac{\partial \mathit{\epsilon }}{\partial t}+\frac{\partial \left(u\mathit{\epsilon }\right)}{\partial x}+\frac{\partial \left(v\mathit{\epsilon }\right)}{\partial y}+\frac{\partial \left(w\mathit{\epsilon }\right)}{\partial z}-\frac{\partial }{\partial z}\left(\frac{A_{\mathrm{v}}}{{\mathit{\sigma }}_{\mathit{\epsilon }}}\frac{\partial \mathit{\epsilon }}{\partial z}\right)\\ =\frac{\mathit{\epsilon }}{k}\left(c_{\mathrm{1}}P+c_{\mathrm{3}}B-c_{\mathrm{2}}\mathit{\epsilon }\right),\end{array}\end{array}$$

with the shear production $P=A_{\mathrm{v}}{S}_{\mathrm{v}}^{\mathrm{2}}$, and the vertical buoyancy flux $B=-K_{\mathrm{v}}{N}^{\mathrm{2}}$. In Eq. (E4), c₁, c₂, c₃, σ_k and σ_ε are empirical parameters. The horizontal diffusive transport of k and ε is generally ignored. To demonstrate how these two-equation turbulence closure models are calibrated to predict the correct response of mixing to stratification and shear, the transport equations for k and ε are analysed for homogenous shear-layers where the advective and diffusive transport-divergence terms vanish. Necessary conditions for equilibrium solutions of Eq. (E4) are found by setting the left-hand sides to zero, such that Eq. (E4) can be transformed to

$$\begin{array}{}\text{(E5)}& \left.\begin{array}{ll}P+B-\mathit{\epsilon }& =\mathrm{0}\\ c_{\mathrm{1}}P+c_{\mathrm{3}}B-c_{\mathrm{2}}\mathit{\epsilon }& =\mathrm{0}\end{array}\right\}\Rightarrow {\displaystyle \frac{c_{\mathrm{1}}-c_{\mathrm{2}}}{c_{\mathrm{3}}-c_{\mathrm{2}}}}={\displaystyle \frac{c_{\mathit{\mu }}^{\prime }\left(R{i}_{\mathrm{st}}\right)}{c_{\mathit{\mu }}\left(R{i}_{\mathrm{st}}\right)}}\cdot R{i}_{\mathrm{st}},\end{array}$$

with the steady-state gradient Richardson number Ri_st denoting the gradient Richardson number for the stationary solution. To obtain the result of Eq. (E5), first ε was eliminated. With the well-calibrated parameters c₁=1.44 and c₂=1.92 (Rodi, 1987) and Ri_st=0.25 (Schumann and Gerz, 1995) Eq. (E5) provides an implicit non-linear equation to calibrate c₃ which determines the damping effect of stratification on turbulence (Burchard and Bolding, 2001). For the second-moment closure by Cheng et al. (2002) this leads to $c_{\mathrm{3}}=-\mathrm{0.74}$ (Umlauf and Burchard, 2005). It can be shown that for a large gradient Richardson number with Ri>Ri_st turbulence is damped due to stable stratification and vice versa (Burchard and Bolding, 2001). This means that an increased steady-state gradient Richardson number with Ri_st>0.25 will lead to increased turbulence for a given Ri and vice versa (Umlauf and Burchard, 2005). This principle has been used by Ralston et al. (2017) to decrease physical mixing by substantially lowering Ri_st in order to account for the too high numerical mixing of an estuarine model, see Sect. 5.2.2. Interestingly, steady-state solutions of Eq. (E4) under the neglect of transport divergences are also provided by

$$\begin{array}{}\text{(E6)}& \begin{array}{rl}& \mathrm{\Gamma }={\left(-{\displaystyle \frac{B}{\mathit{\epsilon }}}\right)}_{\mathrm{st}}={\displaystyle \frac{c_{\mathrm{2}}-c_{\mathrm{1}}}{c_{\mathrm{1}}-c_{\mathrm{3}}}}\\ & \text{and}{\mathrm{Rf}}_{\mathrm{st}}={\left(-{\displaystyle \frac{B}{P}}\right)}_{\mathrm{st}}={\displaystyle \frac{c_{\mathrm{2}}-c_{\mathrm{1}}}{c_{\mathrm{2}}-c_{\mathrm{3}}}}.\end{array}\end{array}$$

For the settings of c₁, c₂ and c₃ given above, a mixing efficiency of Γ=0.21 and a steady-state flux Richardson number of Rf_st=0.17 result (both close to the generic values of Γ=0.2 and Rf_st=0.15 by Osborn, 1980), see the derivations by Umlauf (2009).

Appendix F: Derivations for numerical mixing example

F1 Discretisation of the one-dimensional advection-diffusion equation

To illustrate the quantification of physical and numerical mixing, the simple case of the one-dimensional advection-diffusion equation and its discretisation by means of a first-order upstream scheme for the advection term and a central-difference scheme for the diffusion term is examined here. The one-dimensional advection-diffusion equation is of the following form:

$$\begin{array}{}\text{(F1)}& {\displaystyle \frac{\partial s}{\partial t}}+u{\displaystyle \frac{\partial s}{\partial x}}-K_{\mathrm{h}}{\displaystyle \frac{{\partial }^{\mathrm{2}}s}{\partial x^{\mathrm{2}}}}=\mathrm{0},\end{array}$$

with the constant advection velocity u and the positive and constant physical diffusivity K_h>0. Multiplying Eq. (F1) by 2s results in the one-dimensional version of Eq. (B3):

$$\begin{array}{}\text{(F2)}& {\displaystyle \frac{\partial s^{\mathrm{2}}}{\partial t}}+u{\displaystyle \frac{\partial s^{\mathrm{2}}}{\partial x}}-K_{\mathrm{h}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{s}^{\mathrm{2}}}{\partial x^{\mathrm{2}}}}=-\underset{{\mathit{\chi }}_{\mathrm{s}}}{\underbrace{\mathrm{2}{K}_{\mathrm{h}}{\left({\displaystyle \frac{\partial s}{\partial x}}\right)}^{\mathrm{2}}}},\end{array}$$

with the physical mixing χ_s on the right-hand side as a sink term. For u>0, a straight-forward explicit-in-time discretisation for Eq. (F1) with constant time step Δt and constant spatial increment Δx is given by

$$\begin{array}{}\text{(F3)}& {\displaystyle \frac{s_i^{n+\mathrm{1}}-s_i^n}{\mathrm{\Delta }t}}+u{\displaystyle \frac{s_i^n-s_{i-\mathrm{1}}^n}{\mathrm{\Delta }x}}-K_{\mathrm{h}}{\displaystyle \frac{s_{i+\mathrm{1}}^n-\mathrm{2}{s}_i^n+s_{i-\mathrm{1}}^n}{\mathrm{\Delta }{x}^{\mathrm{2}}}}=\mathrm{0},\end{array}$$

with a first-order upstream discretisation for the advection term and a central-difference discretisation for the diffusion term. In Eq. (F3), the subscripts i indicate the spatial increment and the superscripts n indicate the number of the time step. For a negative velocity u<0, the upstream discretisation of the advection term would simply mean to exchange $s_{i-\mathrm{1}}^n$ by $s_{i+\mathrm{1}}^n$. The scheme is numerically stable for a Courant number of $\mathit{\mu }=\left|u\right|\mathrm{\Delta }t/\mathrm{\Delta }x\le \mathrm{1}$ and the diffusion number $\mathit{\nu }=K_{\mathrm{h}}\mathrm{\Delta }t/\mathrm{\Delta }{x}^{\mathrm{2}}\le \frac{\mathrm{1}}{\mathrm{2}}$. Multiplication of Eq. (F3) by $s_i^{n+\mathrm{1}}+s_i^n$ (equivalent to the multiplication of Eq. (F1) by 2s) and subsequent reorganisation results in a diagnostic equation for the advection and diffusion of s², as reproduced by the discretisation of the advection-diffusion equation Eq. (F1):

$$\begin{array}{}\text{(F4)}& \begin{array}{l}\frac{{\left(s_i^{n+\mathrm{1}}\right)}^{\mathrm{2}}-{\left(s_i^n\right)}^{\mathrm{2}}}{\mathrm{\Delta }t}+u\frac{{\left(s_i^n\right)}^{\mathrm{2}}-{\left(s_{i-\mathrm{1}}^n\right)}^{\mathrm{2}}}{\mathrm{\Delta }x}-K_{\mathrm{h}}\frac{{\left(s_{i+\mathrm{1}}^n\right)}^{\mathrm{2}}-\mathrm{2}{\left(s_i^n\right)}^{\mathrm{2}}+{\left(s_{i-\mathrm{1}}^n\right)}^{\mathrm{2}}}{\mathrm{\Delta }{x}^{\mathrm{2}}}=\\ -\underset{\text{discrete physical mixing},{\mathit{\chi }}_{\mathrm{s},\mathrm{phy},i}}{\underbrace{\begin{array}{c}\mathrm{2}{K}_{\mathrm{h}}\left\{(\mathrm{2}\mathit{\nu }+\mathrm{2}\mathit{\mu }){\left({\displaystyle \frac{s_{i+\mathrm{1}}^n-s_{i-\mathrm{1}}^n}{\mathrm{2}\mathrm{\Delta }x}}\right)}^{\mathrm{2}}+\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}-\mathit{\nu }-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\mu }\right)\right.\\ \left.{\left({\displaystyle \frac{s_{i+\mathrm{1}}^n-s_i^n}{\mathrm{\Delta }x}}\right)}^{\mathrm{2}}+\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}-\mathit{\nu }-{\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}\mathit{\mu }\right){\left({\displaystyle \frac{s_i^n-s_{i-\mathrm{1}}^n}{\mathrm{\Delta }x}}\right)}^{\mathrm{2}}\right\}\end{array}}}\\ \\ -\underset{\text{numerical mixing},{\mathit{\chi }}_{\mathrm{s},\mathrm{num},i}}{\underbrace{\mathrm{2}\frac{u\mathrm{\Delta }x}{\mathrm{2}}\left(\mathrm{1}-\mathit{\mu }\right){\left(\frac{s_i^n-s_{i-\mathrm{1}}^n}{\mathrm{\Delta }x}\right)}^{\mathrm{2}}}},\end{array}\end{array}$$

where the second and third lines are the discrete physical mixing and the last line is the numerical mixing. A step-by-step derivation of Eq. (F4) is given in Sect. F2. The left-hand side of Eq. (F4) is the discretisation of the advection and diffusion of s² and the first term on the right-hand side is a consistent discretisation of the physical mixing, where the non-dimensional numerical parameters μ and ν determine the weights for the composition of the discrete squared salinity gradient by means of the discrete values $s_{i+\mathrm{1}}^n$, $s_i^n$ and $s_{i-\mathrm{1}}^n$. Although this term partially depends on numerical parameters, we associate it with physical mixing, since it is proportional to the eddy diffusivity K_h. The second term on the right-hand side is a sink term for stable conditions with $\mathrm{0}\le \mathit{\mu }\le \mathrm{1}$, and is therefore called the numerical mixing term (Burchard and Rennau, 2008). With the numerical diffusivity $K_{\mathrm{num}}=\frac{\mathrm{1}}{\mathrm{2}}u\mathrm{\Delta }x\left(\mathrm{1}-\mathit{\mu }\right)$ it has a similar structure as the physical mixing, being twice the diffusivity times the tracer-gradient square.

F2 Derivation of physical and numerical mixing for the one-dimensional advection-diffusion equation

Here we show the step-by-step calculation of the physical and numerical mixing for the advection-diffusion equation with first-order upstream (FOU) discretisation for advection and central difference discretisation for diffusion Eq. (F3), which is the starting point here:

$$\begin{array}{}\text{(F5)}& {\displaystyle \frac{s_i^{n+\mathrm{1}}-s_i^n}{\mathrm{\Delta }t}}+u{\displaystyle \frac{s_i^n-s_{i-\mathrm{1}}^n}{\mathrm{\Delta }x}}-K_{\mathrm{h}}{\displaystyle \frac{s_{i+\mathrm{1}}^n-\mathrm{2}{s}_i^n+s_{i-\mathrm{1}}^n}{\mathrm{\Delta }{x}^{\mathrm{2}}}}=\mathrm{0}.\end{array}$$

Multiplication of Eq. (F5) by $s_i^{n+\mathrm{1}}+s_i^n$ results in

$$\begin{array}{}\text{(F6)}& \begin{array}{rl}& {\displaystyle \frac{{\left(s_i^{n+\mathrm{1}}\right)}^{\underset{\mathrm{¯}}{\mathrm{2}}}-{\left(s_i^n\right)}^{\underset{\mathrm{¯}}{\mathrm{2}}}}{\mathrm{\Delta }t}}+u\underset{\mathrm{¯}}{\left(s_i^{n+\mathrm{1}}+s_i^n\right)}{\displaystyle \frac{s_i^n-s_{i-\mathrm{1}}^n}{\mathrm{\Delta }x}}\\ & -K_{\mathrm{h}}\underset{\mathrm{¯}}{\left(s_i^{n+\mathrm{1}}+s_i^n\right)}{\displaystyle \frac{s_{i+\mathrm{1}}^n-\mathrm{2}{s}_i^n+s_{i-\mathrm{1}}^n}{\mathrm{\Delta }{x}^{\mathrm{2}}}}=\mathrm{0},\end{array}\end{array}$$

where we mark changes from step to step by an underline. Substitution of

$$\begin{array}{}\text{(F7)}& \begin{array}{ll}{s}_i^{n+\mathrm{1}}& =s_i^n-\mathit{\mu }\left(s_i^n-s_{i-\mathrm{1}}^n\right)+\mathit{\nu }\left(s_{i+\mathrm{1}}^n-\mathrm{2}{s}_i^n+s_{i-\mathrm{1}}^n\right)\\ & =\mathit{\nu }{s}_{i+\mathrm{1}}^n+(\mathrm{1}-\mathit{\mu }-\mathrm{2}\mathit{\nu })s_i^n+(\mathit{\mu }+\mathit{\nu })s_{i-\mathrm{1}}^n\end{array}\end{array}$$

with the Courant number $\mathit{\mu }=u\mathrm{\Delta }t/\mathrm{\Delta }x$ and the stability number for diffusion $\mathit{\nu }=K_{\mathrm{h}}\mathrm{\Delta }t/\mathrm{\Delta }{x}^{\mathrm{2}}$ gives

$$\begin{array}{}\text{(F8)}& \begin{array}{l}\frac{{\left(s_i^{n+\mathrm{1}}\right)}^{\mathrm{2}}-{\left(s_i^n\right)}^{\mathrm{2}}}{\mathrm{\Delta }t}\\ \\ +\frac{u}{\mathrm{\Delta }x}\left(\underset{\mathrm{¯}}{\mathit{\nu }{s}_{i+\mathrm{1}}^n+(\mathrm{2}-\mathit{\mu }-\mathrm{2}\mathit{\nu })s_i^n+(\mathit{\mu }+\mathit{\nu })s_{i-\mathrm{1}}^n}\right)\\ (s_i^n-s_{i-\mathrm{1}}^n)\\ \\ -\frac{K_{\mathrm{h}}}{\mathrm{\Delta }{x}^{\mathrm{2}}}\left(\underset{\mathrm{¯}}{\mathit{\nu }{s}_{i+\mathrm{1}}^n+(\mathrm{2}-\mathit{\mu }-\mathrm{2}\mathit{\nu })s_i^n+(\mathit{\mu }+\mathit{\nu })s_{i-\mathrm{1}}^n}\right)\\ (s_{i+\mathrm{1}}^n-\mathrm{2}{s}_i^n+s_{i-\mathrm{1}}^n)=\mathrm{0}.\end{array}\end{array}$$

$$\begin{array}{}\text{(F9)}& \begin{array}{l}\frac{{\left(s_i^{n+\mathrm{1}}\right)}^{\mathrm{2}}-{\left(s_i^n\right)}^{\mathrm{2}}}{\mathrm{\Delta }t}\\ \\ +\frac{u}{\mathrm{\Delta }x}(\underset{\mathrm{¯}}{\mathit{\nu }{s}_{i+\mathrm{1}}^n{s}_i^n+(\mathrm{2}-\mathit{\mu }-\mathrm{2}\mathit{\nu })(s_i^n{)}^{\mathrm{2}}+(\mathit{\mu }+\mathit{\nu })s_{i-\mathrm{1}}^n{s}_i^n}\\ \\ \underset{\mathrm{¯}}{-\mathit{\nu }{s}_{i+\mathrm{1}}^n{s}_{i-\mathrm{1}}^n+(-\mathrm{2}+\mathit{\mu }+\mathrm{2}\mathit{\nu })s_i^n{s}_{i-\mathrm{1}}^n+(-\mathit{\mu }-\mathit{\nu })(s_{i-\mathrm{1}}^n{)}^{\mathrm{2}}})\\ \\ -\frac{K_{\mathrm{h}}}{\mathrm{\Delta }{x}^{\mathrm{2}}}(\underset{\mathrm{¯}}{\mathit{\nu }(s_{i+\mathrm{1}}^n{)}^{\mathrm{2}}+(\mathrm{2}-\mathit{\mu }-\mathrm{2}\mathit{\nu })s_i^n{s}_{i+\mathrm{1}}^n+(\mathit{\mu }+\mathit{\nu })s_{i-\mathrm{1}}^n{s}_{i+\mathrm{1}}^n}\\ \\ \underset{\mathrm{¯}}{-\mathrm{2}\mathit{\nu }{s}_{i+\mathrm{1}}^n{s}_i^n+(-\mathrm{4}+\mathrm{2}\mathit{\mu }+\mathrm{4}\mathit{\nu })(s_i^n{)}^{\mathrm{2}}+(-\mathrm{2}\mathit{\mu }-\mathrm{2}\mathit{\nu })s_i^n{s}_{i-\mathrm{1}}^n}\\ \\ \underset{\mathrm{¯}}{+\mathit{\nu }{s}_{i+\mathrm{1}}^n{s}_{i-\mathrm{1}}^n+(\mathrm{2}-\mathit{\mu }-\mathrm{2}\mathit{\nu })s_i^n{s}_{i-\mathrm{1}}^n+(\mathit{\mu }+\mathit{\nu })(s_{i-\mathrm{1}}^n{)}^{\mathrm{2}}})=\mathrm{0}.\end{array}\end{array}$$

Reordering gives

$$\begin{array}{}\text{(F10)}& \begin{array}{l}\frac{{\left(s_i^{n+\mathrm{1}}\right)}^{\mathrm{2}}-{\left(s_i^n\right)}^{\mathrm{2}}}{\mathrm{\Delta }t}=\\ \\ -\frac{u}{\mathrm{\Delta }x}(\mathit{\nu }{s}_{i+\mathrm{1}}^n{s}_i^n+(\mathrm{2}-\mathit{\mu }-\mathrm{2}\mathit{\nu }\left)\right(s_i^n{)}^{\mathrm{2}}+(-\mathrm{2}+\underset{\mathrm{¯}}{\mathrm{2}}\mathit{\mu }+\underset{\mathrm{¯}}{\mathrm{3}}\mathit{\nu })s_{i-\mathrm{1}}^n{s}_i^n\\ \\ -\mathit{\nu }{s}_{i+\mathrm{1}}^n{s}_{i-\mathrm{1}}^n+(-\mathit{\mu }-\mathit{\nu })\left(s_{i-\mathrm{1}}^n{)}^{\mathrm{2}}\right)\\ \\ +\frac{K_{\mathrm{h}}}{\mathrm{\Delta }{x}^{\mathrm{2}}}\left(\mathit{\nu }\right(s_{i+\mathrm{1}}^n{)}^{\mathrm{2}}+(\mathrm{2}-\mathit{\mu }-\underset{\mathrm{¯}}{\mathrm{4}}\mathit{\nu })s_i^n{s}_{i+\mathrm{1}}^n+(\mathit{\mu }+\underset{\mathrm{¯}}{\mathrm{2}}\mathit{\nu })s_{i-\mathrm{1}}^n{s}_{i+\mathrm{1}}^n\\ \\ +(-\mathrm{4}+\mathrm{2}\mathit{\mu }+\mathrm{4}\mathit{\nu })(s_i^n{)}^{\mathrm{2}}+(\mathrm{2}-\underset{\mathrm{¯}}{\mathrm{3}}\mathit{\mu }-\underset{\mathrm{¯}}{\mathrm{4}}\mathit{\nu })s_i^n{s}_{i-\mathrm{1}}^n\\ +(\mathit{\mu }+\mathit{\nu })\left(s_{i-\mathrm{1}}^n{)}^{\mathrm{2}}\right).\end{array}\end{array}$$

Adding on both sides

$$\begin{array}{}\text{(F11)}& \begin{array}{rl}& +{\displaystyle \frac{u}{\mathrm{\Delta }x}}\left(\right(s_i^n{)}^{\mathrm{2}}-\left(s_{i-\mathrm{1}}^n{)}^{\mathrm{2}}\right)\\ & -{\displaystyle \frac{K_{\mathrm{h}}}{\mathrm{\Delta }{x}^{\mathrm{2}}}}\left(\right(s_{i+\mathrm{1}}^n{)}^{\mathrm{2}}-\mathrm{2}(s_i^n{)}^{\mathrm{2}}+(s_{i-\mathrm{1}}^n{)}^{\mathrm{2}})\end{array}\end{array}$$

gives

$$\begin{array}{}\text{(F12)}& \begin{array}{l}\frac{{\left(s_i^{n+\mathrm{1}}\right)}^{\mathrm{2}}-{\left(s_i^n\right)}^{\mathrm{2}}}{\mathrm{\Delta }t}\\ \underset{\mathrm{¯}}{+\frac{u}{\mathrm{\Delta }x}\left(\right(s_i^n{)}^{\mathrm{2}}-\left(s_{i-\mathrm{1}}^n{)}^{\mathrm{2}}\right)}\underset{\mathrm{¯}}{-\frac{K_{\mathrm{h}}}{\mathrm{\Delta }{x}^{\mathrm{2}}}\left(\right(s_{i+\mathrm{1}}^n{)}^{\mathrm{2}}-\mathrm{2}(s_i^n{)}^{\mathrm{2}}+(s_{i-\mathrm{1}}^n{)}^{\mathrm{2}})}=\\ \\ -\frac{u}{\mathrm{\Delta }x}(\mathit{\nu }{s}_{i+\mathrm{1}}^n{s}_i^n+(\underset{\mathrm{¯}}{\mathrm{1}}-\mathit{\mu }-\mathrm{2}\mathit{\nu }\left)\right(s_i^n{)}^{\mathrm{2}}+(-\mathrm{2}+\mathrm{2}\mathit{\mu }+\mathrm{3}\mathit{\nu })s_{i-\mathrm{1}}^n{s}_i^n\\ \\ -\mathit{\nu }{s}_{i+\mathrm{1}}^n{s}_{i-\mathrm{1}}^n+(\underset{\mathrm{¯}}{\mathrm{1}}-\mathit{\mu }-\mathit{\nu })\left(s_{i-\mathrm{1}}^n{)}^{\mathrm{2}}\right)\\ \\ +\frac{K_{\mathrm{h}}}{\mathrm{\Delta }{x}^{\mathrm{2}}}\left(\right(\underset{\mathrm{¯}}{-\mathrm{1}}+\mathit{\nu }\left)\right(s_{i+\mathrm{1}}^n{)}^{\mathrm{2}}+(\mathrm{2}-\mathit{\mu }-\mathrm{4}\mathit{\nu })s_i^n{s}_{i+\mathrm{1}}^n+(\mathit{\mu }+\mathrm{2}\mathit{\nu })s_{i-\mathrm{1}}^n{s}_{i+\mathrm{1}}^n\\ \\ +(\underset{\mathrm{¯}}{-\mathrm{2}}+\mathrm{2}\mathit{\mu }+\mathrm{4}\mathit{\nu })(s_i^n{)}^{\mathrm{2}}+(\mathrm{2}-\mathrm{3}\mathit{\mu }-\mathrm{4}\mathit{\nu })s_i^n{s}_{i-\mathrm{1}}^n\\ +(\underset{\mathrm{¯}}{-\mathrm{1}}+\mathit{\mu }+\mathit{\nu })\left(s_{i-\mathrm{1}}^n{)}^{\mathrm{2}}\right).\end{array}\end{array}$$

Equation (F12) is the diagnostic discrete equation for the advection-diffusion equation of squared salinity, where the right-hand side shows the total discrete mixing, consisting of a discretisation of the physical mixing plus contributions from numerical mixing. To differentiate numerical and physical mixing, we split the right-hand side into a purely advective contribution (only containing $u/\mathrm{\Delta }x$ and μ) and the remainder:

$$\begin{array}{}\text{(F13)}& \begin{array}{l}\frac{{\left(s_i^{n+\mathrm{1}}\right)}^{\mathrm{2}}-{\left(s_i^n\right)}^{\mathrm{2}}}{\mathrm{\Delta }t}+\frac{u}{\mathrm{\Delta }x}\left(\right(s_i^n{)}^{\mathrm{2}}-\left(s_{i-\mathrm{1}}^n{)}^{\mathrm{2}}\right)\\ -\frac{K_{\mathrm{h}}}{\mathrm{\Delta }{x}^{\mathrm{2}}}\left(\right(s_{i+\mathrm{1}}^n{)}^{\mathrm{2}}-\mathrm{2}(s_i^n{)}^{\mathrm{2}}+(s_{i-\mathrm{1}}^n{)}^{\mathrm{2}})=\\ \\ -\frac{u}{\mathrm{\Delta }x}\left(\right(\mathrm{1}-\mathit{\mu }\left)\right(s_i^n{)}^{\mathrm{2}}+(-\mathrm{2}+\mathrm{2}\mathit{\mu })s_{i-\mathrm{1}}^n{s}_i^n+(\mathrm{1}-\mathit{\mu })\left(s_{i-\mathrm{1}}^n{)}^{\mathrm{2}}\right)\\ \\ \underset{\mathrm{¯}}{-\mathit{\nu }\frac{u}{\mathrm{\Delta }x}(s_{i+\mathrm{1}}^n{s}_i^n-\mathrm{2}(s_i^n{)}^{\mathrm{2}}+\mathrm{3}{s}_{i-\mathrm{1}}^n{s}_i^n-s_{i+\mathrm{1}}^n{s}_{i-\mathrm{1}}^n-\left(s_{i-\mathrm{1}}^n{)}^{\mathrm{2}}\right)}\\ \\ +\frac{K_{\mathrm{h}}}{\mathrm{\Delta }{x}^{\mathrm{2}}}\left(\right(-\mathrm{1}+\mathit{\nu }\left)\right(s_{i+\mathrm{1}}^n{)}^{\mathrm{2}}+(\mathrm{2}-\mathit{\mu }-\mathrm{4}\mathit{\nu })s_i^n{s}_{i+\mathrm{1}}^n\\ +(\mathit{\mu }+\mathrm{2}\mathit{\nu })s_{i-\mathrm{1}}^n{s}_{i+\mathrm{1}}^n\\ \\ +(-\mathrm{2}+\mathrm{2}\mathit{\mu }+\mathrm{4}\mathit{\nu })(s_i^n{)}^{\mathrm{2}}+(\mathrm{2}-\mathrm{3}\mathit{\mu }-\mathrm{4}\mathit{\nu })s_i^n{s}_{i-\mathrm{1}}^n\\ +(-\mathrm{1}+\mathit{\mu }+\mathit{\nu })\left(s_{i-\mathrm{1}}^n{)}^{\mathrm{2}}\right).\end{array}\end{array}$$

Noting that

$$\begin{array}{}\text{(F14)}& \mathit{\nu }{\displaystyle \frac{u}{\mathrm{\Delta }x}}=\mathit{\mu }{\displaystyle \frac{K_{\mathrm{h}}}{\mathrm{\Delta }{x}^{\mathrm{2}}}}\end{array}$$

and reformulating the purely advective term, we obtain

$$\begin{array}{}\text{(F15)}& \begin{array}{l}\frac{{\left(s_i^{n+\mathrm{1}}\right)}^{\mathrm{2}}-{\left(s_i^n\right)}^{\mathrm{2}}}{\mathrm{\Delta }t}+\frac{u}{\mathrm{\Delta }x}\left(\right(s_i^n{)}^{\mathrm{2}}-\left(s_{i-\mathrm{1}}^n{)}^{\mathrm{2}}\right)\\ -\frac{K_{\mathrm{h}}}{\mathrm{\Delta }{x}^{\mathrm{2}}}\left(\right(s_{i+\mathrm{1}}^n{)}^{\mathrm{2}}-\mathrm{2}(s_i^n{)}^{\mathrm{2}}+(s_{i-\mathrm{1}}^n{)}^{\mathrm{2}})=\\ \\ -\frac{u}{\mathrm{\Delta }x}\underset{\mathrm{¯}}{(\mathrm{1}-\mathit{\mu })}\left(\right(s_i^n{)}^{\mathrm{2}}-\mathrm{2}{s}_{i-\mathrm{1}}^n{s}_i^n+\left(s_{i-\mathrm{1}}^n{)}^{\mathrm{2}}\right)\\ \\ +\frac{K_{\mathrm{h}}}{\mathrm{\Delta }{x}^{\mathrm{2}}}\left(\right(-\mathrm{1}+\mathit{\nu }\left)\right(s_{i+\mathrm{1}}^n{)}^{\mathrm{2}}+(\mathrm{2}-\underset{\mathrm{¯}}{\mathrm{2}}\mathit{\mu }-\mathrm{4}\mathit{\nu })s_i^n{s}_{i+\mathrm{1}}^n\\ +(\underset{\mathrm{¯}}{\mathrm{2}}\mathit{\mu }+\mathrm{2}\mathit{\nu })s_{i-\mathrm{1}}^n{s}_{i+\mathrm{1}}^n\\ \\ +(-\mathrm{2}+\underset{\mathrm{¯}}{\mathrm{4}}\mathit{\mu }+\mathrm{4}\mathit{\nu })(s_i^n{)}^{\mathrm{2}}+(\mathrm{2}-\underset{\mathrm{¯}}{\mathrm{6}}\mathit{\mu }-\mathrm{4}\mathit{\nu })s_i^n{s}_{i-\mathrm{1}}^n\\ +(-\mathrm{1}+\underset{\mathrm{¯}}{\mathrm{2}}\mathit{\mu }+\mathit{\nu })\left(s_{i-\mathrm{1}}^n{)}^{\mathrm{2}}\right).\end{array}\end{array}$$

Now the right-hand side will be expressed as gradient-squared terms. This is simple for the advective term. For the diffusive term, we apply the following Ansatz to express it as a linear combination of the three possible discrete gradient-square terms:

$$\begin{array}{}\text{(F16)}& \begin{array}{l}a{\left(\frac{s_{i+\mathrm{1}}^n-s_i^n}{\mathrm{\Delta }x}\right)}^{\mathrm{2}}+b{\left(\frac{s_i^n-s_{i-\mathrm{1}}^n}{\mathrm{\Delta }x}\right)}^{\mathrm{2}}+c{\left(\frac{s_{i+\mathrm{1}}^n-s_{i-\mathrm{1}}^n}{\mathrm{2}\mathrm{\Delta }x}\right)}^{\mathrm{2}}=\\ \\ \frac{\mathrm{1}}{\mathrm{\Delta }{x}^{\mathrm{2}}}\left[\left(a+\frac{c}{\mathrm{4}}\right)\right(s_{i+\mathrm{1}}{)}^{\mathrm{2}}+\left(a+b\right)(s_i^n{)}^{\mathrm{2}}+\left(b+\frac{c}{\mathrm{4}}\right)(s_{i-\mathrm{1}}{)}^{\mathrm{2}}\\ \\ -\mathrm{2}a{s}_{i+\mathrm{1}}^n{s}_i^n-\mathrm{2}b{s}_i^n{s}_{i-\mathrm{1}}^n-\frac{c}{\mathrm{2}}{s}_{i+\mathrm{1}}^n{s}_{i-\mathrm{1}}^n]\end{array}\end{array}$$

Comparison of coefficients results in

$$\begin{array}{}\text{(F17)}& \begin{array}{ll}a+\frac{c}{\mathrm{4}}=-\mathrm{1}+\mathit{\nu },& a+b=-\mathrm{2}+\mathrm{4}\mathit{\mu }+\mathrm{4}\mathit{\nu },\\ b+\frac{c}{\mathrm{4}}=-\mathrm{1}+\mathrm{2}\mathit{\mu }+\mathit{\nu },& \\ & \\ -\mathrm{2}a=\mathrm{2}-\mathrm{2}\mathit{\mu }-\mathrm{4}\mathit{\nu },& -\mathrm{2}b=\mathrm{2}-\mathrm{6}\mathit{\mu }-\mathrm{4}\mathit{\nu },\\ -\frac{c}{\mathrm{2}}=\mathrm{2}\mathit{\mu }+\mathrm{2}\mathit{\nu },& \end{array}\end{array}$$

with the solution

$$\begin{array}{}\text{(F18)}& \begin{array}{ll}a=-\mathrm{1}+\mathit{\mu }+\mathrm{2}\mathit{\nu },& b=-\mathrm{1}+\mathrm{3}\mathit{\mu }+\mathrm{2}\mathit{\nu },\\ c=-\mathrm{4}\mathit{\mu }-\mathrm{4}\mathit{\nu }.& \end{array}\end{array}$$

With this, Eq. (F15) can be expressed as

$$\begin{array}{}\text{(F19)}& \begin{array}{l}\frac{{\left(s_i^{n+\mathrm{1}}\right)}^{\mathrm{2}}-{\left(s_i^n\right)}^{\mathrm{2}}}{\mathrm{\Delta }t}+\frac{u}{\mathrm{\Delta }x}\left(\right(s_i^n{)}^{\mathrm{2}}-\left(s_{i-\mathrm{1}}^n{)}^{\mathrm{2}}\right)\\ -\frac{K_{\mathrm{h}}}{\mathrm{\Delta }{x}^{\mathrm{2}}}\left(\right(s_{i+\mathrm{1}}^n{)}^{\mathrm{2}}-\mathrm{2}(s_i^n{)}^{\mathrm{2}}+(s_{i-\mathrm{1}}^n{)}^{\mathrm{2}})=\\ \\ -u\mathrm{\Delta }x(\mathrm{1}-\mathit{\mu }){\left(\frac{s_i^n-s_{i-\mathrm{1}}^n}{\mathrm{\Delta }x}\right)}^{\mathrm{2}}\\ \\ +K_{\mathrm{h}}\left(\right(-\mathrm{1}+\mathit{\mu }+\mathrm{2}\mathit{\nu }){\left(\frac{s_{i+\mathrm{1}}^n-s_i^n}{\mathrm{\Delta }x}\right)}^{\mathrm{2}}\\ +(-\mathrm{1}+\mathrm{3}\mathit{\mu }+\mathrm{2}\mathit{\nu }){\left(\frac{s_i^n-s_{i-\mathrm{1}}^n}{\mathrm{\Delta }x}\right)}^{\mathrm{2}}\\ \\ +(-\mathrm{4}\mathit{\mu }-\mathrm{4}\mathit{\nu }){\left(\frac{s_{i+\mathrm{1}}^n-s_{i-\mathrm{1}}^n}{\mathrm{2}\mathrm{\Delta }x}\right)}^{\mathrm{2}}),\end{array}\end{array}$$

which can be reformulated as

$$\begin{array}{}\text{(F20)}& \begin{array}{l}\frac{{\left(s_i^{n+\mathrm{1}}\right)}^{\mathrm{2}}-{\left(s_i^n\right)}^{\mathrm{2}}}{\mathrm{\Delta }t}+\frac{u}{\mathrm{\Delta }x}\left(\right(s_i^n{)}^{\mathrm{2}}-\left(s_{i-\mathrm{1}}^n{)}^{\mathrm{2}}\right)\\ -\frac{K_{\mathrm{h}}}{\mathrm{\Delta }{x}^{\mathrm{2}}}\left(\right(s_{i+\mathrm{1}}^n{)}^{\mathrm{2}}-\mathrm{2}(s_i^n{)}^{\mathrm{2}}+(s_{i-\mathrm{1}}^n{)}^{\mathrm{2}})=\\ \\ -\mathrm{2}u\frac{\mathrm{\Delta }x}{\mathrm{2}}(\mathrm{1}-\mathit{\mu }){\left(\frac{s_i^n-s_{i-\mathrm{1}}^n}{\mathrm{\Delta }x}\right)}^{\mathrm{2}}\\ \\ -\mathrm{2}{K}_{\mathrm{h}}[\left(\frac{\mathrm{1}}{\mathrm{2}}-\frac{\mathit{\mu }}{\mathrm{2}}-\mathit{\nu }\right){\left(\frac{s_{i+\mathrm{1}}^n-s_i^n}{\mathrm{\Delta }x}\right)}^{\mathrm{2}}\\ +\left(\frac{\mathrm{1}}{\mathrm{2}}-\frac{\mathrm{3}}{\mathrm{2}}\mathit{\mu }-\mathit{\nu }\right){\left(\frac{s_i^n-s_{i-\mathrm{1}}^n}{\mathrm{\Delta }x}\right)}^{\mathrm{2}}\\ \\ +\left(\mathrm{2}\mathit{\mu }+\mathrm{2}\mathit{\nu }\right){\left(\frac{s_{i+\mathrm{1}}^n-s_{i-\mathrm{1}}^n}{\mathrm{2}\mathrm{\Delta }x}\right)}^{\mathrm{2}}],\end{array}\end{array}$$

and which is identical to Eq. (F4). It should be noted that for stationary problems (which must not depend on the time step Δt) such as the discrete solution of the exponential estuary in Appendic F3, the sum of all terms in Eq. (F20) containing one of the numerical parameters μ or ν is zero:

$$\begin{array}{}\text{(F21)}& \begin{array}{l}-\mathrm{2}{K}_{\mathrm{h}}\left\{(\mathrm{2}\mathit{\nu }+\mathrm{2}\mathit{\mu }){\left(\frac{s_{i+\mathrm{1}}^n-s_{i-\mathrm{1}}^n}{\mathrm{2}\mathrm{\Delta }x}\right)}^{\mathrm{2}}\right.\\ +\left(-\mathit{\nu }-\frac{\mathrm{1}}{\mathrm{2}}\mathit{\mu }\right){\left(\frac{s_{i+\mathrm{1}}^n-s_i^n}{\mathrm{\Delta }x}\right)}^{\mathrm{2}}\\ \left.+\left(-\mathit{\nu }-\frac{\mathrm{3}}{\mathrm{2}}\mathit{\mu }\right){\left(\frac{s_i^n-s_{i-\mathrm{1}}^n}{\mathrm{\Delta }x}\right)}^{\mathrm{2}}\right\}\\ \\ -\mathrm{2}\frac{u\mathrm{\Delta }x}{\mathrm{2}}\left(-\mathit{\mu }\right){\left(\frac{s_i^n-s_{i-\mathrm{1}}^n}{\mathrm{\Delta }x}\right)}^{\mathrm{2}}=\mathrm{0}.\end{array}\end{array}$$

Figure F1 Analytical and numerical solutions for the one-dimensional stationary estuary with run-off velocity u=0.05 m s⁻¹, K_h=500 m² s⁻¹, L=100 km and A=10 000 m² (such that the run-off is $Q_{\mathrm{r}}=uA=\mathrm{500}$ m³ s⁻¹ and the river salinity is $s_{\mathrm{r}}=\mathrm{1.4}\times {\mathrm{10}}^{-\mathrm{3}}$ g kg⁻¹). The numerical solution is carried out for $i_{max}=\mathrm{20}$ equidistant spatial increments of $\mathrm{\Delta }x=L/i_{max}$. (a) Analytical solution Eq. (F23) and numerical solution (bullets) for salinity s as function of non-dimensional estuarine length $x/L$. (b) Analytical solution Eq. (33) and numerical results for the salinity mixing per salinity class m(S): Total mixing (squares), physical mixing 𝕞_phy,i (bullets) and numerical mixing 𝕞_num,i (circles).

F3 One-dimensional stationary estuarine example

We demonstrate the distribution of physical and numerical mixing for a simple example of a one-dimensional stationary estuary of length L with constant cross-section A and constant discharge velocity u>0 (such that the river run-off is Q_r=uA) and constant along-estuary diffusivity K_h>0, being slightly simplified from Burchard (2020). Under those circumstances, the stationary form of Eq. (F1) describes the balance,

$$\begin{array}{}\text{(F22)}& {\displaystyle \frac{\partial }{\partial x}}\left(us-K_{\mathrm{h}}{\displaystyle \frac{\partial s}{\partial x}}\right)=\mathrm{0},\end{array}$$

of which the analytical solution for salinity is exponentially increasing from the river towards the ocean,

$$\begin{array}{}\text{(F23)}& s\left(x\right)=s_{\mathrm{o}}\exp \left({\displaystyle \frac{u}{K_{\mathrm{h}}}}(x-L)\right),\end{array}$$

with the ocean salinity s_o at x=L and the river salinity $s_{\mathrm{r}}=s_{\mathrm{o}}\exp (-uL/K_{\mathrm{h}})$ at x=0. The universal law of estuarine mixing Eq. (33) can be directly derived from Eq. (F23), see Appendix F4 for details.

The analytical solution Eq. (F23) and a numerical solution, obtained by iterating Eq. (F3) into a stationary state, are given in Fig. F1a. Using $i_{max}=\mathrm{20}$ equidistant spatial increments results in a good numerical representation of the analytical solution with some deviations at high salinities. The universal law of estuarine mixing Eq. (33) is shown as analytical solution and numerical approximation in Fig. F1b. There, the numerically approximated physical mixing 𝕞_phy,i and the numerical mixing 𝕞_num,i are calculated as

$$\begin{array}{}\text{(F24)}& \begin{array}{rl}{\mathbb{m}}_{\mathrm{phy},i}& ={\displaystyle \frac{{\mathit{\chi }}_{\mathrm{s},\mathrm{phy},i}}{\frac{\mathrm{1}}{\mathrm{2}}(s_{i+\mathrm{1}}^n-s_{i-\mathrm{1}}^n)}}A\mathrm{\Delta }x,\\ {\mathbb{m}}_{\mathrm{num},i}& ={\displaystyle \frac{{\mathit{\chi }}_{\mathrm{s},\mathrm{num},i}}{\frac{\mathrm{1}}{\mathrm{2}}(s_{i+\mathrm{1}}^n-s_{i-\mathrm{1}}^n)}}A\mathrm{\Delta }x,\end{array}\end{array}$$

with the total mixing 𝕞_tot,i being the sum of the two. In Eq. (F24), the integration over a grid cell is carried out by multiplication with the grid box volume AΔx. Note that the discrete salinity classes are not equidistant in salinity space here, but defined in terms of the stationary salinity distribution on the equidistant spatial grid such that the size of the salinity classes is $\frac{\mathrm{1}}{\mathrm{2}}(s_{i+\mathrm{1}}^n-s_{i-\mathrm{1}}^n)$. The numerical mixing amounts to about 12 % of the total mixing for this simple estuarine example with a small tidally averaged advection velocity (Fig. F1b). Clearly, the universal law of estuarine mixing Eq. (33) is only fulfilled when considering the total mixing consisting of physical and numerical mixing. Note that for this stationary solution, the results for mixing must not depend on the time step Δt. It can indeed be shown that in this case the sum of all terms in ${\mathit{\chi }}_{\mathrm{s},\mathrm{phy},i}$ and ${\mathit{\chi }}_{\mathrm{s},\mathrm{num},i}$ of Eq. (F4) containing the numerical parameters μ and ν is exactly zero, independently of Δt, see Eq. (F21) of Appendix F2. For a simulation of a realistic estuary using a high-resolution coastal ocean model, the relation of physical to numerical mixing per salinity class is shown in Fig. 8, with the numerical mixing also being of the order of 10 % of the total mixing.

F4 Analytical solution for salinity mixing per salinity class

Here, we derive the universal law of estuarine mixing directly from the analytical solution for the simple stationary one-dimensional estuary Eq. (F23). This solution differs slightly from the solution given by Burchard (2020). Inverting Eq. (F23) to resolve for x gives

$$\begin{array}{}\text{(F25)}& x\left(S\right)={\displaystyle \frac{K_{\mathrm{h}}}{u}}\ln \left({\displaystyle \frac{S}{s_{\mathrm{o}}}}\right)+L\end{array}$$

where S is now used as the salinity coordinate. The volume-integrated mixing M(S) for the monotone salinity distribution Eq. (F23) can be calculated as the integral over the local mixing (using A as cross-sectional area),

$$\begin{array}{}\text{(F26)}& \mathbb{M}\left(S\right)=A\underset{\mathrm{0}}{\overset{x\left(S\right)}{\int }}\mathrm{2}{K}_{\mathrm{h}}{\left({\displaystyle \frac{\partial s\left(x^{\prime }\right)}{\partial x^{\prime }}}\right)}^{\mathrm{2}}\text{d}{x}^{\prime }.\end{array}$$

Inserting the square of the x-derivative of Eq. (F23) into Eq. (F26) gives

$$\begin{array}{}\text{(F27)}& \begin{array}{rl}\mathbb{M}\left(S\right)& =\mathrm{2}A{s}_{\mathrm{o}}^{\mathrm{2}}{\displaystyle \frac{u^{\mathrm{2}}}{K_{\mathrm{h}}}}\underset{\mathrm{0}}{\overset{x\left(S\right)}{\int }}\exp \left(\mathrm{2}{\displaystyle \frac{u}{K_{\mathrm{h}}}}(x^{\prime }-L)\right)\text{d}{x}^{\prime }\\ & =A{s}_{\mathrm{o}}^{\mathrm{2}}u\underset{-\mathrm{2}\frac{u}{K_{\mathrm{h}}}L}{\overset{\mathrm{2}\frac{u}{K_{\mathrm{h}}}(x-L)}{\int }}\exp \left(x^{\prime }\right)\text{d}{x}^{\prime }\\ & =Q_{\mathrm{r}}{s}_{\mathrm{o}}^{\mathrm{2}}\left({\exp }^{\mathrm{2}}\left({\displaystyle \frac{u}{K_{\mathrm{h}}}}(x-L)\right)-{\exp }^{\mathrm{2}}\left(-{\displaystyle \frac{u}{K_{\mathrm{h}}}}L\right)\right)\\ & =Q_{\mathrm{r}}{s}_{\mathrm{o}}^{\mathrm{2}}\left({\displaystyle \frac{S^{\mathrm{2}}}{s_{\mathrm{o}}^{\mathrm{2}}}}-{\displaystyle \frac{s_{\mathrm{r}}^{\mathrm{2}}}{s_{\mathrm{o}}^{\mathrm{2}}}}\right)=Q_{\mathrm{r}}\left(S^{\mathrm{2}}-s_{\mathrm{r}}^{\mathrm{2}}\right),\end{array}\end{array}$$

where Eq. (F25) and Q_r=Au have been used. Finally, this yields the mixing per salinity class,

$$\begin{array}{}\text{(F28)}& \mathbb{m}\left(S\right)={\displaystyle \frac{\partial \mathbb{M}\left(S\right)}{\partial S}}=\mathrm{2}S{Q}_{\mathrm{r}},\end{array}$$

which is identical to Eq. (33). Following Lorenz et al. (2021), and in analogy to their case of freshwater input across the ocean surface due to precipitation, the additional term $-Q_{\mathrm{r}}{s}_{\mathrm{r}}^{\mathrm{2}}$ in Eq. (F27) represents the corresponding boundary flux of S² at the river end and reduces the mixing in the interior.

Table F1 List of major variables including their physical meaning and dimensions.

Table F2 List of operators demonstrated for an arbitrary variable x.