Effective Diahaline Diffusivities in Estuaries

Abstract

The present study aims to estimate effective diahaline turbulent salinity fluxes and diffusivities in numerical model simulations of estuarine scenarios. The underlying method is based on a quantification of salinity mixing per salinity class, which is shown to be twice the turbulent salinity transport across the respective isohaline. Using this relation, the recently derived universal law of estuarine mixing, predicting that average mixing per salinity class is twice the respective salinity times the river run‐off, can be directly derived. The turbulent salinity transport is accurately decomposed into physical (due to the turbulence closure) and numerical (due to truncation errors of the salinity advection scheme) contributions. The effective diahaline diffusivity representative for a salinity class and an estuarine region results as the ratio of the diahaline turbulent salinity transport and the respective (negative) salinity gradient, both integrated over the isohaline area in that region and averaged over a specified period. With this approach, the physical (or numerical) diffusivities are calculated as half of the product of physical (or numerical) mixing and the isohaline volume, divided by the square of the isohaline area. The method for accurately calculating physical and numerical diahaline diffusivities is tested and demonstrated for a three‐dimensional idealized exponential estuary. As a major product of this study, maps of the spatial distribution of the effective diahaline diffusivities are shown for the model estuary.

Introduction

The circulation and hydrography of estuaries is largely determined by turbulent mixing. This has been impressively demonstrated by the laboratory experiments conducted by Linden and Simpson (1986), generating lock exchange flows under the impact of air bubble‐induced turbulence. They could show that strong turbulence slowed down the exchange flow driven by gravitational forces and used this as an explanation for the tidal modulation of estuarine circulation (Linden & Simpson, 1988). Earlier, using a simple analytical model, Hansen and Rattray (1965) quantified how the intensity of estuarine exchange flow is inversely proportional to specified constant vertical eddy viscosity and diffusivity. In an idealized numerical model study, Hetland and Geyer (2004) showed how the length of an estuary, that is, the longitudinal extent of the brackish water zone, depends on a prescribed eddy viscosity and diffusivity: high eddy coefficients result in short estuaries and vice versa. Observations and model results for many estuaries reflect this process: During phases of high turbulence (spring tides), mixing reduces vertical stratification and generates more upright isohaline surfaces, resulting in shorter estuaries and vice versa during low turbulence (neap tide), see for example, Warner et al. (2005) for the Hudson River Estuary and Li et al. (2018) for the Changjiang River Estuary.

In contrast to local and instantaneous diffusivities as parameterized in ocean models by turbulence closures (e.g., Large et al., 1994; Umlauf & Burchard, 2005), effective diffusivities are often calculated to analyze diffusive transport in water bodies of different characteristics. The latter are generally calculated as the (negative) ratio of turbulent flux and tracer gradient, both individually averaged over a surface (e.g., an isopycnal) and over time, such that they result as weighted integral averages of the local and instantaneous diffusivities. Different ways in calculating mean gradients result in different definitions of effective diffusivities. The physically unique and well‐defined turbulent properties related to salinity transformation are the diahaline turbulent salinity flux and the salinity mixing (i.e., the turbulent salinity variance dissipation) defined below. However, for comparison with local and instantaneous diffusivities, the effective diffusivity is a useful quantity across a large range of salinity‐gradient scales. The substantially inhomogeneous distribution of eddy diffusivity (low in the interior, high at sloping bottoms) has the consequence that effective diffusivities in various types of basins are easily an order of magnitude larger than the local ones observed in the interior, such as found, for example, for the world ocean (Waterhouse et al., 2014), the Baltic Sea (Holtermann et al., 2012) or lakes (Goudsmit et al., 1997).

In ocean models, diffusivities are typically contaminated by numerical spurious mixing by tracer advection schemes, such that water mass transformations parameterized by carefully calibrated turbulence closure models might be overridden by numerical diffusion. Various methods have been proposed to estimate effective diffusivities in ocean models, either based on the increase of background potential energy (Griffies et al., 2000; Ilıcak et al., 2012), on numerical tracer releases (Getzlaff et al., 2010, 2012) or on water mass transformation analysis evaluating overturning stream functions (Lee et al., 2002; Megann, 2018). All these methods share the drawback that the analysis averages effective diffusivities over large regions. Methods for locally analyzing physical and numerical mixing based on the decay of tracer variance have been proposed by Burchard and Rennau (2008) and Klingbeil et al. (2014). These methods, however, do not yet allow computing effective diapycnal diffusivities.

In estuaries, effective longitudinal diffusivities, that is, the ratio of longitudinal turbulent salinity fluxes and salinity gradients individually averaged over vertical transects, are several orders of magnitude larger than vertical diffusivities (Fischer, 1976). This is mainly due to the process of shear dispersion (Taylor, 1954; Young & Jones, 1991), a combination of vertical diffusion and vertical shear. In order to reproduce salt intrusion, vertically or cross‐sectionally integrated estuarine models which do not resolve tidal shear dispersion do therefore need to parameterize this process through large values of longitudinal diffusivity (Ralston et al., 2008; Uncles & Stephens, 1990).

In early coastal ocean models, horizontal diffusivity was applied to suppress numerically induced oscillations (Blumberg & Mellor, 1987). Recent models use monotone advection schemes including implicit numerical diffusion, such that explicit horizontal diffusion is not needed (e.g., Hofmeister et al., 2011; Ralston et al., 2017). The numerical mixing caused by these advection schemes as well as the physical mixing due to the turbulence closure model can be quantified as local and instantaneous dissipation of salinity variance (Burchard & Rennau, 2008; Klingbeil et al., 2014). In some estuarine model simulations, the numerical mixing can account for a substantial part of the total mixing, such that a re‐tuning of the turbulence closure is required to obtain satisfactory salinity distributions (Ralston et al., 2017).

In estuaries, density is largely determined by salinity, although for low‐inflow estuaries, specifically in summer, also temperature might contribute to density stratification (Largier et al., 1996). Therefore, for typical strong‐inflow estuaries, temperature variations are often neglected when analyzing the dynamics, with the consequence that isohaline and isopycnal surfaces are assumed to be identical. In recent years, one focus of estuarine physical studies has been the analysis of salinity mixing, defined as the decay of salinity variance (Burchard & Rennau, 2008; Burchard et al., 2009). Based on the Total Exchange Flow (TEF) analysis framework introduced by MacCready (2011), the relation between the exchange flow and integrated estuarine salinity mixing of riverine fresh water and oceanic salt water was studied by several authors (Burchard et al., 2019; MacCready et al., 2018; Wang et al., 2017). To better understand estuarine dynamics, analysis of water mass transformations on isohaline surfaces (Hetland, 2005; MacCready & Geyer, 2001; MacCready et al., 2002; Walin, 1977) is often advantageous over the analysis in geographical coordinates, since it is not contaminated by adiabatic (advective) processes such as tidal oscillations. Using the definition of mixing per salinity class (Wang et al., 2017), that is, mixing integrated over infinitesimal isohaline volumes, Burchard (2020) were able to show that for long time integrations it converges toward twice the product of the respective salinity times the freshwater run‐off. By using the analysis of physical and numerical mixing introduced by Klingbeil et al. (2014), the total mixing per salinity class could be decomposed into physical and numerical contributions, where the latter could become negative at times when antidiffusive advection schemes are applied.

Due to the pivotal role of eddy diffusivity in estuaries, the present study aims to estimate diahaline turbulent salinity fluxes and diffusivities in an estuarine numerical model simulation. This will be achieved by exploiting the recent concept of mixing on isohaline surfaces. The method presented here resembles in parts the calculation of eddy diffusivity as the product of scalar flux averaged over isoscalar surfaces and the inverse of the diascalar gradient for the vertical coordinate z_* of the sorted scalar field, as proposed by Winters and D'Asaro (1996).

Mathematical Derivation

This section aims to review and further derive the isohaline framework within which the diahaline turbulent transport and the effective diahaline diffusivities are defined. The basic definitions of instantaneous mixing, fluxes, and diffusivities are given in Section 2.1. The isohaline geometry including isohaline surfaces and volumes is introduced in Section 2.2. Based on this, diahaline mixing, transport, and diffusivities are defined in Section 2.3. Finally, in Section 2.4, diahaline mixing relations and effective diahaline diffusivities for long‐term averaged estuarine states are discussed in the light of these new results. Note that a list of variables including their meanings and units is given in Table 1.

1 TableList of Variables Including Their Meanings and Units

Symbol	Meaning	Unit
a	TEF‐based isohaline area	m2
A	Instantaneous isohaline surface	m2
Ar	River transect area	m2
b	TEF‐based isohaline thickness	m (g/kg)−1
bloc	Local thickness per salinity class	m (g/kg)−1
Bloc	Accumulated local thickness	m
fs	Diahaline salinity flux	ms−1 (g/kg)
$f_{\text{adv}}^s$	Advective diahaline salinity flux	ms−1 (g/kg)
$f_{\text{diff}}^s$	Diffusive diahaline salinity flux	ms−1 (g/kg)
Fs	Diahaline salinity transport	m3s−1 (g/kg)
$F_{\text{adv}}^s$	Advective diahaline salinity transport	m3s−1 (g/kg)
$F_{\text{diff}}^s$	Diffusive diahaline salinity transport	m3s−1 (g/kg)
${\vec{\mathcal{F}}}_{\text{diff}}^s$	Diffusive salinity flux vector	ms−1 (g/kg)
Kh	Horizontal eddy diffusivity	m2s−1
Kn	Diahaline eddy diffusivity	m2s−1
$\overline{K_n}$	Effective total diahaline diffusivity	m2s−1
$\overline{K_n^{\text{num}}}$	Effective numerical diahaline diffusivity	m2s−1
$\overline{K_n^{\text{phy}}}$	Effective physical diahaline diffusivity	m2s−1
Kv	Vertical eddy diffusivity	m2s−1
m	Mixing per salinity class	m3s−1 (g/kg)
mnum	Numerical mixing per salinity class	m3s−1 (g/kg)
mphy	Physical mixing per salinity class	m3s−1 (g/kg)
M	Volume‐integrated mixing	m3s−1 (g/kg)2
$\vec{u}$	Velocity vector	ms−1
${\vec{u}}^s$	Velocity of point on isohaline surface	ms−1
s	Absolute salinity	gkg−1
v	Volume per salinity class	m3 (g/kg)−1
V	Accumulated estuarine volume	m3
wdia	Diahaline volume transport per unit area	ms−1
χ	Salinity mixing per unit volume	s−1 (g/kg)2
zTEF	TEF‐based isohaline position	m

Salinity Mixing and Diahaline Diffusivity

The basic physical equations from which the present theory is derived are the continuity equation (volume conservation, assuming incompressibility of sea water), [Image Omitted. See PDF]and the salinity equation (salinity conservation), [Image Omitted. See PDF]where $\vec{u}$ is the velocity vector, s is absolute salinity, and [Image Omitted. See PDF]is the turbulent (diffusive) salinity flux vector in Cartesian coordinates, where K_h and K_v are horizontal and vertical eddy diffusivities, respectively.

The salinity mixing per unit volume, χ, is defined as the local loss of salinity variance (Burchard & Rennau, 2008), [Image Omitted. See PDF]It should be noted that the single components of the turbulent salinity flux vector are downgradient, but due to the nonisotropic eddy diffusivity (K_h ≫ K_v) the flux vector is not orthogonal to the isohaline surface, see Figure 1 and the discussion below.

1 Figure. Sketch explaining the projection of the turbulent diffusive flux vector F⃗diffs to the diahaline direction. The unit normal vector n⃗ is pointing outside the domain toward higher salinity.

In ocean models, the horizontal turbulent fluxes are typically aligned with the vertical model coordinate, for example, with constant z‐levels for geopotential models, constant σ‐levels for models with terrain and surface following σ‐coordinates (with σ = (z − η)/(H + η), the sea surface height z = η and the bottom coordinate z = −H) or constant density for isopycnal models. Thus, for nongeopotential models, vertical turbulent fluxes are nonorthogonal to the lateral fluxes. In nonisopycnal models, the rotation of horizontal turbulent fluxes into isopycnal direction causes numerical problems (Beckers et al., 1998; Griffies et al., 17) and is therefore sometimes avoided.

Let ${\vec{u}}^s$ denote the velocity of a point moving with an isohaline surface with [Image Omitted. See PDF]Then, [Image Omitted. See PDF](with the normal vector $\vec{n}=\vec{\nabla }s/\left|\vec{\nabla }s\right|$, pointing toward higher salinity) is the diahaline volume transport per unit area at that point, and [Image Omitted. See PDF]is the total (advective plus diffusive) diahaline salinity flux, see Section 6.7 of Griffies (2004), Figure 3 of Groeskamp et al. (2014), and Section 7.1 of Groeskamp et al. (2019) for a full derivation. In 7, the salinity gradient normal to the isohaline surface is defined as ${\partial }_{n}s=\vec{\nabla }s\cdot \vec{n}$ and K_n is the diahaline diffusivity. Note that the diahaline turbulent salinity flux, $f_{\text{diff}}^s={\vec{\mathcal{F}}}_{\text{diff}}^s\cdot \vec{n}$, is the orthogonal projection of the turbulent salinity flux vector to the isohaline surface, see Figure 1 for an illustration. The definition 4 of the mixing per unit volume and the definition 7 of the diahaline turbulent salinity flux both imply the following expression for the diahaline diffusivity, [Image Omitted. See PDF]which means that the diahaline diffusivity is a weighted average of the horizontal and vertical eddy diffusivities depending on the salinity gradients. A consequence of this is that the diascalar diffusivities are different for each scalar tracer although the horizontal and vertical diffusivities are not.

Isohaline Volumes

We consider here a time‐averaged estuarine volume V(S) bounded in seaward direction by an isohaline of arbitrary salinity S with area A(S), that is, the volume contains estuarine water masses with a salinity s with s ≤ S. The volume is further bounded by the air‐sea interface at z = η, the river boundary area A_r located at zero salinity and zero salinity gradient, and the impermeable sea bottom at z = −H (Figure 2). Similar definitions have been used by Hetland (2005), see his conceptual Figure 2. The time‐averaging operator is introduced for any function X(t) as [Image Omitted. See PDF]where T is the averaging interval. Note that for the limit of T → 0, the nonaveraged, instantaneous value $〈X〉(t)=X(t)$ is obtained. In the following, we generally drop the argument t for time‐averaged quantities. With this, V(S) is formally defined as [Image Omitted. See PDF]with the Heaviside function [Image Omitted. See PDF]

2 Figure. Sketch illustrating the estuarine volume V(S) (blue), bounded by the isohaline of salinity S with area A(S) (red), the surface (marked by blue lines), the bottom (marked by yellow lines), and the river transect (facing to the right) with area Ar. Across the river transects, riverine water of salinity s = 0 and a rate of Qr is entering the volume. The salinity transport across the isohaline is denoted as Fs(S). As indicated, the vertical scale H is much smaller than the horizontal scale L, such that the isohaline surface is virtually flat.

The isohaline volume had been defined by Walin (1977) as the infinitesimal volume per salinity class, for which different formulations can be used: [Image Omitted. See PDF]where ${\left({\partial }_{n}s\right)}^{-1}$ is the local and instantaneous infinitesimal thickness of the isohaline. Figure 3 illustrates the definitions of volume per salinity class and thickness per salinity class by means of a finite salinity increment ΔS.

3 Figure. Sketch explaining the isohaline volume and thickness by means of a finite salinity increment ΔS. The gray area indicates the finite volume ΔV between the isohaline surfaces A(S + ΔS/2) and A(S − ΔS/2). For the limit of ΔS → 0, ΔV/ΔS → ∂SV = v, that is, the infinitesimal volume per salinity class is obtained. The local thickness of the finite volume ΔV is denoted as Δn, that for ΔS → 0 the infinitesimal local thickness per salinity class, Δn/ΔS→(∂ns)−1 is obtained.

It should be noted that the isohaline volumes defined here depend on salinity S only. To retrieve information how these properties are distributed in horizontal space, the local thickness per salinity class, b ^loc, and the accumulated local thickness with salinities ≤S, B^loc, can be defined as [Image Omitted. See PDF]We define the TEF‐based mean isohaline surface height as [Image Omitted. See PDF]which for a monotonically stable salinity stratification is the average vertical position of the isohaline. This is motivated by the TEF analysis framework proposed by MacCready (2011) which is comparable to the TRM‐II theory proposed by McDougall and McIntosh (2001). Equation 14 indicates an isohaline position inside the water column for all locations (x, y) where the salinity S has occurred during the averaging period for a finite period of time (Figure 4). The isohaline surface height z^TEF resembles the isoscalar coordinate z_* defined by Winters and D'Asaro (1996), with the difference that z^TEF is defined for each water column individually, whereas z_* is a one‐dimensional coordinate representing an entire volume, see also equivalent approaches by McDougall and McIntosh (2001) and Nurser and Lee (2004).The surface area of the TEF‐based isohaline height is denoted as a(S). With 12, [Image Omitted. See PDF]is the area‐averaged infinitesimal thicknesses per salinity class of the TEF‐based salinity distribution. The area‐averaged salinity gradients are then estimated as b⁻¹. In comparison to directly averaging the diahaline salinity gradient ∂_ns over the isohaline surface of a specific salinity class (a task that is numerically demanding in a discrete model), this method is biased toward lower values of gradients and therefore toward higher values of effective diffusivities. On the other hand, the method is easy to implement in this isohaline framework, see also Section 3. As a method to estimate mean diapycnal density gradients in ocean mixing studies, the density field may be resorted adiabatically into a state of minimum available potential energy with flat isopycnals (Griffies et al., 2000; Winters & D'Asaro, 1996), resulting in a single profile of effective diffusivity as function of density only. This method, however, is computationally demanding in situations with variable bottom topography.

4 Figure. Sketch explaining the local thickness per salinity class, bloc, the accumulated local thickness with salinities ≤S, Bloc, and the TEF‐based mean isohaline surface height, zTEF, in a discrete way with a finite salinity increment ΔS. In the limit of ΔS → 0, the discrete local thickness per salinity class converges toward the (infinitesimal) thickness per salinity class: bΔSloc(S)=(Bloc(S+ΔS/2)−Bloc(S−ΔS/2))/ΔS→bloc(S).

It should be noted that the calculations of the isohaline areas and the derived isohaline thicknesses are correct for the projections of the isohalines to geopotential surfaces. Increased areas and decreased thicknesses due to sloping isohalines are not taken into consideration. However, these errors are assumed to be small due to the small aspect ratio between vertical and horizontal scales in estuaries and the resulting small slopes of the isohalines. Additionally, isohalines are TEF‐based time‐averaged surfaces such that steeper isohaline slopes due to transient dynamics such as submesoscale or frontal dynamics are smoothed. In numerical models, it is possible to correct for these errors, but for simplicity we refrain here to do so.

Diahaline Transport

The transports of salinity through the isohaline surface are denoted as [Image Omitted. See PDF]with the total salinity transport being denoted as $F^s=F_{\text{adv}}^s+F_{\text{diff}}^s$ and the diahaline salinity fluxes defined in 7.

According to Burchard (3), the mixing per salinity class m(S) is defined as [Image Omitted. See PDF]with the volume‐integrated salinity mixing M(S). Note that 17 has already been formulated by Wang et al. (2017) for the vertical component of mixing, see their Equation 4.3. Applying a generalized form of the Leibniz theorem (see the equation above Equation 3.5 of Marshall et al., 1999), m(S) can also be formulated as [Image Omitted. See PDF]which is a type of transformation between surface integrals and volume integrals that serves as the foundation for many applications of water mass transformations (Groeskamp et al., 2019). Equation 18 is equivalent to the definition of the volume per salinity class given in 12, see Figure 3, with χ included.

Combining 7, 8, 16, and 18, the mixing per salinity class can also be formulated as [Image Omitted. See PDF]

With this, the mixing per salinity class is twice the (negative) diahaline turbulent salinity transport. This also means that the mixing integrated over the volume V(S) can be expressed by means of integrating the diahaline turbulent salinity transport through all isohalines with salinities smaller than or equal to S: [Image Omitted. See PDF]see Equation 4.4 by Wang et al. (2017). Equation 19 allows for easy calculation of diahaline turbulent transport and effective diahaline diffusivities including their decomposition into physically and numerically induced components (see below). This formulation does also allow formulating the recently proposed universal law of estuarine mixing (Burchard, 2020) without using a budget equation for the salinity variance or the integrated salinity square, see Section 2.4.

Following Klingbeil et al. (2014), the mixing per unit volume χ can be exactly decomposed into physical and numerical contributions, such that consequently also the mixing per isohaline volume can be decomposed: [Image Omitted. See PDF]

Motivated by the relations in 8, an effective total diahaline diffusivity, $\overline{K_n}$, can be calculated by dividing the negative mean diahaline turbulent flux averaged over the isohaline surface, $-F_{\text{diff}}^s/a=\frac{1}{2}m/a$, by the mean diahaline salinity gradient, b⁻¹ = a/v: [Image Omitted. See PDF]

which is the key result of the present study. When formulated as [Image Omitted. See PDF] 22 resembles Equation 12 by Winters and D'Asaro (1996) for calculating the turbulent eddy diffusivity by means of the average diascalar flux times the inverse of the scalar gradient with respect to their z_* coordinate. The formulation 22 has the advantage that the three properties m, v, and a can be calculated in an easy way in numerical models. Moreover, using 21, the effective diahaline diffusivity can be split into physical and numerical parts: [Image Omitted. See PDF]Details of the discretization of 24 are given in Section 3.

Diahaline Mixing Relations

As shown in Burchard (2020) for tidally periodic estuaries averaged over a tidal period or general estuaries averaged over a long period, the downestuarine advective salinity transport through an isohaline of salinity S is $F_{\text{adv}}^s=S{Q}_r$, where Q_r is the time‐averaged river run‐off. Under these equilibrium conditions, the downestuarine advective salinity transport must be exactly balanced by an opposing upestuarine diffusive salinity transport $F_{\text{diff}}^s=-F_{\text{adv}}^s$. Using 19, we obtain $m=-2F_{\text{diff}}^s=2F_{\text{adv}}^s=2S{Q}_r$, stating that under equilibrium conditions the mixing per salinity class is twice the product of the river run‐off times the salinity of the respective isohaline. This has been recently formulated by Burchard (2020) as the universal law of estuarine mixing, derived in a more complicated way by combining volume‐integrated budget equations for volume, salinity, and salinity‐squared. A similar relation has already been derived by Zika et al. (2015) for the fresh water cycle of the global ocean.

Discretization

The numerical model calculates for each time step n and for each grid box i, k discrete values for the time‐variable layer thickness, $h_{i,k}^n$, the salinity, $s_{i,k}^n$, the physical mixing per unit volume, ${\chi }_{i,k}^{\text{phy},n}$ (due to diffusion with eddy diffusivity calculated by the turbulence closure scheme), and the numerical mixing per unit volume, ${\chi }_{i,k}^{\text{num},n}$ (due to truncation errors of the salinity advection scheme). We calculated the latter two quantities utilizing the method proposed by Klingbeil et al. (2014) as local discrete variance decay due to diffusion and advection, respectively. Furthermore, the horizontal grid cell area is kept constant in time and is denoted as A_i. The index i is here the counter for the two‐dimensional horizontal grid which may be structured or unstructured, and k is the vertical counter of layers.

For the binning into salinity classes, the salinity range s_min ≤ s ≤ s_max is discretized into J equidistant salinity intervals Δs = (s_max − s_min)/J. The time‐averaged volume per salinity class, v_j, is then calculated for each water column i as [Image Omitted. See PDF]where s_j = s_min + jΔs and ⌈x⌉ is the ceiling function to any nonnegative real number x, j with 1 ≤ j ≤ J is the counter for the salinity class and N is the number of time steps. With 25, a vector of salinity classes j is defined for each water column i which is then subsequently filled in a binning process with subvolumes $A_i{h}_{i,k}^n$ and averaged by the number of time steps N. Division by Δs is necessary in 25, because v_i,j is volume per salinity class.

The physical and numerical mixing per salinity class, $m_{i,j}^{\text{phy}}$ and $m_{i,j}^{\text{num}}$, are calculated accordingly as [Image Omitted. See PDF]Let A = ∑_i∈IA_i be a specifically defined subarea (e.g., an estuarine segment or the entire estuary) composed of all areas A_i with the index i included in the set I. Then, following 24, with 25 and 26, the physical and numerical effective diahaline diffusivities representative for A(S = s_j) are calculated as [Image Omitted. See PDF]with the isohaline area [Image Omitted. See PDF]With 28, a water column with an empty salinity class j is associated with zero isohaline area. We speculate that it will happen, when a salinity class is outside the range of salinities occurring in a water column (such as ocean salinity in permanently brackish water). This may also happen for high numbers of salinity classes and low numbers of salinity values (due to low vertical and temporal resolution or short averaging periods). In both cases, such empty salinity classes pose no problem. If a salinity class does not occur in any water column of the chosen subarea, the effective diahaline diffusivities for that salinity class are not defined. It should be mentioned here, that the isohaline area A(S = s_j) is calculated for simplicity as the projection of the isohaline surface to geopotential surfaces, an error that is small due to the small aspect ratio of estuaries (see Figure 2 and also the results of Section 4) and it is identical for the calculation of both, the effective physical and numerical diahaline diffusivity. In regions where this approximation is not sufficiently accurate, a more exact reconstruction of the isohaline area could be made.

Idealized Model Experiments

To demonstrate the calculation of effective diahaline diffusivities and other mixing‐related properties, we simulate an idealized estuary, exponentially widening toward the open ocean (Figure 5). We placed the estuary at a latitude of 53.5°N and prescribed a length of 100 km, a central navigational channel of 15 m depth, and lateral shoals with an average depth of 3 m. At the mouth, the estuary is 81 km wide to allow for the development of a river plume, and decreases in width exponentially in landward direction. The minimum width of the estuarine channel is set to 1 km. The model is forced at the open boundary with a mono‐chromatic semidiurnal tide of 2.0 m amplitude and a constant ocean salinity of 35 g/kg. At the river end of the estuary, a constant freshwater run‐off of Q_r = 700 m³s⁻¹ is prescribed. There is no wind forcing applied.

5 Figure. Idealized convergent estuary model: salinity distribution during high water as gauged at x = 0. The upper left panel shows surface salinity, the lower left panel shows a longitudinal transect (at y = 0 km), and the right panels show a cross‐sectional transect at three different locations (note the different horizontal scales).

For the simulations, the General Estuarine Transport Model (GETM, www.getm.eu, Burchard & Bolding, 2002) has been applied, a primitive equation coastal ocean model model using general vertical coordinates and explicit mode splitting (Klingbeil et al., 2018). The discretization in GETM is based on the finite‐volume principle such that volume and salt are exactly conserved. GETM is coupled to the turbulence module of the General Ocean Turbulence Model (GOTM, www.gotm.net, Burchard & Bolding, 2001; Umlauf & Burchard, 2005), using the k‐ɛ two‐equation turbulence closure model with an algebraic second‐moment closure by Cheng et al. (2002). Explicit horizontal diffusion is not applied.

We constructed a curvilinear grid with 200 cells in longitudinal direction and 30 cells across the estuary. In the vertical, 30 terrain and surface following σ‐layers are used with some grid refinement toward the bottom. For the temporal discretization, each semidiurnal tidal cycle is resolved with 5,000 baroclinic time steps, each of them split into 10 barotropic time steps. The advection terms in the momentum, salinity, and turbulence budgets are discretized employing the TVD‐SPL‐max‐1/3 scheme (Waterson & Deconinck, 2007), known for its minimum numerical mixing (Mohammadi‐Aragh et al., 2015), combined with Strang splitting (Pietrzak, 1998). The isohaline analysis is carried out for a discrete salinity increment of ΔS = 0.1 g/kg, and it was found that the results do only weakly depend on ΔS. However, for a too small ΔS and too short averaging intervals, results may become noisy, since the number of discrete values per salinity class may be too low. A similar problem has been discussed by Lorenz et al. (2019) for TEF analysis across a fixed transect. In contrast, for too large ΔS, results will have a too low resolution in salinity space.

The simulation is started from rest (zero surface elevation, zero velocity, constant salinity s = 15 g/kg) until a quasiperiodic state is reached, still including nontidal oscillations such as internal waves of frequencies not matching the tidal frequency. Ten tidal periods of this quasiperiodic state are analyzed to calculate representative tidally averaged properties.

A snapshot of the salinity distribution at high water is shown in Figure 5. A salt wedge is reaching up to 70 km of the estuary with a strong near‐surface stratification downstream of it. A river plume veering to the north (positive y‐direction) due to Earth rotation is visible in the outer estuary. The lateral salinity structure is complex due to lateral circulation (which is not shown).

The salinity field from the TEF‐based averaging is shown in Figure 6 for the center line of the estuary (y = 0). For each water column, the TEF‐based average salinity is located at z^TEF(S) defined in 14. With this, the lowest salinity occurring in that water column is located at the time‐averaged surface, and the highest salinity is located at the bottom. Here, it can also be seen how small the error of approximating the TEF‐based time‐averaged isohaline area by its projection to geopotential coordinates is: for the example of the 25‐g/kg‐isohaline based on TEF (Figure 6), the relative error of the isohaline area is $1-L/{(H^2+L^2)}^{1/2}=2.5\cdot 10^{-7}$, with the length of the projected isohaline of L = 30 km and the centerline depth of H = 15 m, when assuming a flat bottom and a planar approximation of the isohaline.

6 Figure. Idealized convergent estuary model: averaged salinity distribution and isohaline positions at the centerline of the estuary (y = 0) resulting from TEF‐based averaging. Isohalines are shown in steps of ΔS = 1 g/kg. The red line represents the position of the isohaline with S = 25 g/kg.

Since the physical lateral diffusivity is set to zero in the model experiment, it is expected according to 8 and 22 that the effective physical diahaline diffusivities, $\overline{K_n^{\text{phy}}}$, are a weighted average of the vertical diffusivities, K_v, occurring in a specific salinity class. This should also be the case for the effective total diahaline diffusivities, $\overline{K_n}$, when physical mixing is dominating over numerical mixing. It is therefore instructive to inspect the distributions of K_v for different situations. This is shown in Figure 7 as snapshots during full flood and full ebb for the centerline of the estuary. During flood, due to a destabilization of the lower half of the water column in the salt wedge vertical diffusivity is enhanced in this region with values of around K_v = 10⁻² m²s⁻¹. During ebb, marginal shear instability is dominating in parts of the salt wedge (with salinities above 12 g/kg) such that vertical diffusivities are still elevated, but when strong stratification occurs too close to the bed, the eddy diffusivity in the water column above drops to values of below 10⁻⁵ m²s⁻¹. In the strongly stratified regions with the salt wedge, diffusivities drop down to the background value of about 10⁻⁶ m²s⁻¹ given by the turbulence closure. In the well‐mixed regions upstream (near s = 0 g/kg) and downstream (near s = 35 g/kg) of the salt wedge, diffusivities resemble parabolic profiles and reach high values of up to 10⁻¹ m²s⁻¹ in about mid water column.

7 Figure. Idealized convergent estuary model: eddy diffusivity Kv (color code) and salinity s (isolines) during (a) full flood and (b) full ebb at the centerline of the estuary (y = 0).

The analysis of mixing per salinity class shows that the universal law of estuarine mixing proposed by Burchard (2020) is closely approximated by averaging over 10 tidal periods in a quasiperiodic state (Figure 8a). For salinities below 22 g/kg (the maximum salinity that is not reaching the open boundary) physical mixing is dominating the total mixing with numerical mixing only contributing by about 30% at most. Only in the coarse resolution region of the river plume at salinities larger than 30 g/kg does numerical mixing dominate (Figure 8a). As further input to the calculation of effective diahaline diffusivity, volume per salinity class, and isohaline area (Figure 8b), and the mean salinity gradient of each salinity class, b⁻¹ (Figure 8c) are shown. Maximum values of b⁻¹ of more than 10 (g/kg) m⁻¹ are reached between salinities 16 and 23 g/kg. This means that for those salinity classes the isohalines are so much stretched out in the longitudinal direction that the thickness per salinity class decreases below 0.1 m (g/kg)⁻¹, see also Figure 6. With m, v, and a (Figures 8a and 8b), all parameters determining the effective diahaline diffusivities are present. The resulting diffusivities are shown in Figure 8d. $\overline{K_n}$ increases about linearly from 1 · 10⁻⁵ m²s⁻¹ at a salinity of 1 g/kg to a maximum of 6 · 10⁻⁵ m²s⁻¹ at 7 g/kg and then decreases to 1 · 10⁻⁵ m²s⁻¹ at 23 g/kg. A peak in effective diahaline diffusivity is generated by a combination of the mixing per salinity class and the isohaline area increasing with salinity and the volume per salinity class substantially dropping down above a salinity of 7 g/kg. For salinities larger than 23 g/kg, the effective diahaline diffusivity reaches a relatively high level of ≥2 · 10⁻⁵ m²s⁻¹, with a peak of 1 · 10⁻⁴ m²s⁻¹ at 24 g/kg, due to a high isohaline volume per salinity class and a consequently small diahaline gradient. As for the mixing per salinity class, also for the effective diahaline diffusivity, the numerical contribution is relatively high for these high salinity classes. Generally, effective diahaline diffusivities are substantially lower than the relative high instantaneous values of K_v > 10⁻² m²s⁻¹ shown in Figure 7 for the well‐mixed regions in the freshwater range, the bottom boundary layer and the well‐mixed region near the mouth. This can be explained by the fact that the mixing per unit volume, χ, is generally low in these regions, which has an impact on the mixing per salinity class, m.

8 Figure. Idealized convergent estuary model averaged over 10 tidal periods during a quasiperiodic state: properties defining effective diahaline diffusivities as function of salinity. (a) Physical, numerical, and total mixing per salinity class in comparison to the universal law m(S) = 2SQr; (b) area of isohaline S, a(S), and volume per salinity class, v(S); (c) resulting averaged salinity gradient, b−1; (d) effective physical, numerical, and total diahaline diffusivity. The vertical hatched line indicates the minimum salinity which reaches the open boundary.

To analyze the spatial distribution of total mixing per salinity class, an integration along the transverse coordinate is carried out, such that total mixing per salinity class per longitudinal distance is calculated, see Figure 9. As a result, elevated mixing occurs over a large range of distance and salinity classes. Integration of the total mixing per salinity class and meter with respect to x results in total mixing per salinity class approximating the universal law (left graph). When integrating with respect to salinity, total mixing per meter is obtained, see the upper graph in Figure 9, showing maximum values in the range −90 km ≤ x ≤ −80 km. However, this distribution does not follow a specific law. Together with the volume per salinity class per longitudinal distance (Figure 10a), the effective total diahaline diffusivity can be calculated (Figure 10b). There is a broad region in salinity and longitudinal distance where the effective total diahaline diffusivity is of the order of 10⁻⁴ m²s⁻¹. At some distinct locations (−80 km ≤ x ≤ −70 km, S ≤ 15 g/kg), elevated values of up to 5 · 10⁻⁴ m²s⁻¹ occur due to a combination of localized high mixing and high volume.

9 Figure. Idealized convergent estuary model averaged over 10 tidal periods during a quasiperiodic state: (a) mixing per salinity class and longitudinal distance, ∂xm(S), (b) mixing per salinity class, m(S) = ∂SM(S) (integrated along the estuary), and (c) mixing per longitudinal distance, ∂xM(S) (integrated over salinity).

10 Figure. Idealized convergent estuary model averaged over 10 tidal periods during a quasiperiodic state: (a) volume per salinity class and longitudinal distance, ∂xv(S), and (b) effective total diahaline diffusivity as function of salinity class and longitudinal distance, Kn¯.

Using the information on mixing per salinity class and volume per salinity class for every water column, it is possible to calculate maps of the effective diahaline diffusivity, as shown in Figures 11a–11c for the example of S = 25 g/kg. In the well‐resolved channelized region of the estuary (x ≥ −90 km), diffusivities are dominated by physical values. They are highest in the deeper parts of the channel and reach values of several 10⁻⁴ m²s⁻¹. Over the shoals, the diffusivities are about 1 order of magnitude smaller. In the region of the outer estuary (x ≤ −90 km), where the river plume veers to the right in downstream direction due to Earth rotation, the horizontal resolution becomes coarse, such that numerical mixing is comparable to physical mixing (Figures 11b and 11c). With this, the maps of Figures 11a–11c represent the spatial distribution of the effective diahaline total, physical, and numerical diffusivities that are accumulated to one value for each salinity in Figure 8d. As an example, the effective physical diahaline diffusivity is $\overline{K_n^{\text{phy}}}(S=25\text{g/kg})=2·{10}^{-5}{\mathrm{m}}^2{\mathrm{s}}^{-1}$ (Figure 8d), which can be interpreted as a weighted integral average of the spatial distribution of $\overline{K_n^{\text{phy}}}(S=25\text{g/kg},x,y)$ over all water columns x and y where the salinity of 25 g/kg is occurring at least once during the averaging interval.

11 Figure. Idealized convergent estuary model averaged over 10 tidal periods during a quasiperiodic state: (a) total, (b) physical, and (c) numerical effective diahaline diffusivity across the 25‐g/kg‐isohaline as function of x and y. Note that the x and the y axes have different scalings.

Discussion

The major goal of this study is to develop a method to calculate effective diahaline turbulent fluxes and diffusivities, including a decomposition into physical and numerical contributions. The resulting method is based on the calculation of mixing per salinity class from which diahaline transports can directly be derived, see 19, as well as on the calculation of volume per salinity class, a concept which had already been proposed by Walin (1977). The method is robust in the sense that no interpolations to isohaline surfaces in numerical models are needed and that the resulting effective turbulent diahaline transport is always negative (i.e., downgradient) and effective diahaline diffusivities are positive definite, except for unlikely situations when negative mixing by antidiffusive advection schemes dominate the effective total diahaline diffusivity. The only inaccuracy in the calculation of the effective diahaline diffusivities is the determination of the isohaline surface for which we use the projection to geopotential surfaces. However, as visible in Figure 6a, even in estuaries isohaline surfaces time‐averaged by TEF analysis are relatively flat, such that the error associated with their slope should be negligible.

The main product of the present study are maps of the effective diahaline diffusivities, as shown in Figure 11 for the isohaline of S = 25 g/kg. Such maps could be shown for all salinity classes and give a clear picture about the mixing hotspots and physical and numerical contributions. Effective diahaline diffusivities as statistical diagnostics for estuarine models are distinct from instantaneous diffusivities which in well‐mixed regions can be orders of magnitude larger. The calculation of effective diahaline diffusivities can be viewed as a complex weighted spatial and temporal averaging process, where high instantaneous diffusivities in well‐mixed regions have small weights due to their small contribution to mixing.

Although diffusivities are defined as the (negative) ratio between a turbulent flux and the associated gradient of the transported property and as such their unweighted average in time and space is not a useful property, they are key quantities for diagnosing mixing processes in the ocean. Using the Osborn and Cox (1972) and Osborn (1980) methods based on equilibrium versions of the tracer variance and turbulent kinetic energy budgets respectively, diapycnal eddy diffusivities can be calculated which in turn can be used to calculate diapycnal tracer fluxes.

In ocean models, complex and computationally demanding turbulence closure schemes are used to compute eddy diffusivities (e.g., Large et al., 1994; Umlauf & Burchard, 2005). Except for situations of double diffusion (Canuto et al., 2002) eddy diffusivities are assumed to be same for all tracers, a concept that underlines their relevance. However, considering arbitrary scalar tracers in ocean models, when explicit horizontal diffusivity is not aligned with the diascalar surface, instantaneous diascalar turbulent diffusivities are different for each scalar tracer, as seen from 8. Also, numerical mixing (and therefore numerical diffusivity) depends strongly on the tracer gradients and is thus different for all tracers (Burchard & Rennau, 2008). Moreover, effective diascalar diffusivities as calculated from 22 depend on the scalar tracer under consideration. On the other hand, a field and modeling study in the deep water of the Central Baltic Sea showed that the diagnosed profile of effective diapycnal diffusivity could explain the evolution of salinity, temperature, and an injected tracer over several years (Holtermann et al., 2012, 2014). The method developed here can only calculate the effective tracer diffusivities orthogonal to isosurfaces of the respective tracer. Still, as shown in Figure 11, it is a useful diagnostics for estuarine mixing processes and their representation in numerical models.

Whereas the salinity mixing per salinity class as well as the diahaline salinity transport averaged over a long time converge toward theoretical limits (Burchard, 2020), no such theories exist for isohaline volumes or effective diahaline diffusivities. Those depend on the individual dynamics of the estuary which in turn is mainly triggered by the bathymetry and the freshwater and tidal forcing.

This method should be applicable to all regions of the ocean. In estuarine systems, where density is dominated by the salinity contribution, the effective diahaline diffusivity is a good estimate for the effective diapycnal diffusivity. The method would require modifications in the case of salinity inversions that occur, for example, in the Baltic Sea in the surface mixed layer during summer (Burchard et al., 2017) or warm and salty inflows overlaying denser but less salty waters (Umlauf et al., 2018). The multiple isohaline surfaces that appear in such situations need to be identified and accounted for. In such situations, and also for larger scale applications, the extension of the present method toward effective diapycnal density transport and diffusivity would be beneficial, a challenging task, since contributions from mixing of temperature and salinity need to be considered. Since in most ocean models density is directly transported neither by advection nor by diffusion, physical and numerical mixing of density cannot be calculated directly. Instead, in order to diagnose diapycnal turbulent density fluxes and effective density diffusivities in ocean models, physical and numerical mixing of temperature and salinity need to be combined in a suitable way. Such a diagnostics would allow adding spatial resolution to the bulk estimates of diapycnal diffusivities proposed by various authors (Getzlaff et al., 2010, 2012; Griffies et al., 2000; Ilıcak et al., 2012; Lee et al., 2002; Megann, 2018).

Conclusions

A method has been developed to diagnose effective diahaline diffusivities in numerical models of estuaries. It is based on the calculation of volume per salinity class and mixing per salinity class, the former having been proposed already by Walin (1977) whereas the latter has only recently been developed by Wang et al. (2017) and Burchard (2020). In numerical models, effective diahaline diffusivities consist of physical (from turbulence parameterizations) and numerical (from discretization errors of advection schemes) contributions and add exactly up to the effective total diahaline diffusivity. The calculation, therefore, requires the analysis of physical and numerical mixing as introduced by Burchard and Rennau (2008) and Klingbeil et al. (2014).

The method computes generally positive values for effective physical and numerical diahaline diffusivities for each salinity class and each water column (unless antidiffusive advection schemes are used), whenever the respective salinity has occurred in the water column during the averaging period. Based on this water column information, effective diahaline diffusivities can be aggregated in various dimensions, either as horizontal maps of diffusivity for a specified salinity class, as a function of salinity and longitudinal distance, or as a function of salinity only, always showing physical and numerical values as well as the sum of both.

These diagnostics help to understand where mixing is strong in an estuary and indicate where numerical mixing is high. The latter helps to plan measures for reducing numerical artifacts by choosing a higher model resolution, better advection schemes or, as the ultimate measure, to reduce physical mixing in order to limit effective total diahaline mixing to realistic levels (Ralston et al., 2017). In the idealized simulations carried out in the present study, numerical diffusivities were on acceptable levels in the inner estuary, but due to the coarse resolution of the outer estuary, they dominated the dynamics in the region.

Acknowledgments

This paper is a contribution to the Collaborative Research Centre TRR 181 on Energy Transfers in Atmosphere and Ocean and the Research Training Group Baltic TRANSCOAST GRK 2000 funded by the German Research Foundation. The authors are grateful to the comments by Stephen Griffies and four anonymous reviewers. The revisions of a previous version of this manuscript have been made at the Bornö Research Station in Sweden.

Data Availability Statement

The tidally averaged model data for 10 tidal periods as well as the python script for the TEF data analysis can be found at https://www.io-warnemuende.de/BurchardJAMES2021.html. The model simulations have been carried out with the General Estuarine Transport Model (GETM) that is freely available at https://www.getm.eu.