## Abstract

This numerical modeling study quantifies for the first time the contribution of various processes to estuarine circulation in periodically stratified tidal flow under the impact of a constant horizontal buoyancy gradient. The one-dimensional water column equations with periodic forcing are first cast into nondimensional form, resulting in a multidimensional parameter space spanned by the modified inverse Strouhal number and the modified horizontal Richardson number, as well as relative wind speed and wind direction and the residual runoff. The along-tide momentum equation is then solved for the tidal-mean velocity profile in such a way that it is equated to the sum of the contributions of tidal straining (resulting from the temporal correlation between eddy viscosity and vertical shear), gravitational circulation (resulting from the depth-varying forcing by a constant horizontal buoyancy gradient), wind straining, and depth-mean residual flow (resulting from net freshwater runoff). This definition of tidal straining does not only account for tidal asymmetries resulting from horizontal buoyancy gradients but also from wind straining and residual runoff. For constant eddy viscosity, the well-known estuarine circulation analytical solution with polynomial residual profiles is directly obtained. For vertically parabolic and constant-in-time eddy viscosity, a new analytic solution with logarithmic residual profiles is found, showing that the intensity of the gravitational circulation scales with the horizontal Richardson number. For scenarios with realistic spatially and temporally varying eddy viscosity, a numerical water column model equipped with a state-of-the-art two-equation turbulence closure model is applied to quantify the individual contributions of the various processes to estuarine circulation. The fundamental outcome of this study is that, for irrotational flow with periodic stratification and without wind forcing and residual runoff, the tidal straining is responsible for about two-thirds and gravitational circulation is responsible for about one-third of the estuarine circulation, proportionally dependent on the horizontal Richardson number, and weakly dependent on the Strouhal number. This new and robust result confirms earlier estimates by H. Burchard and H. Baumert, who suggested that tidal straining is the major generation mechanism for estuarine turbidity maxima. However, a sensitivity analysis of the model results to details of the turbulence closure model shows some uncertainty with respect to the parameterization of sheared convection during flood. Increasing down-estuary wind straining and residual runoff reduce the quantitative contribution of tidal straining. For relatively small horizontal Richardson numbers, the tidal straining contribution to estuarine circulation may even be reversed by down-estuary wind straining.

## 1. Introduction

Along-channel transport in many estuaries is caused primarily by estuarine circulation: subtidal currents that bring freshwater seaward in the upper portion of the water column, with an associated near-bed landward flow bringing salty water. Estuarine circulation can interact with plankton or sediments to create a net landward or seaward transport that can be much larger than the mean transport associated with the river flow. Positively buoyant material tends to be transported seaward, whereas negatively buoyant materials, such as suspended sediments, are transported landward.

Early studies of estuarine circulation (Pritchard 1952, 1954, 1956; Hansen and Rattray 1965; Chatwin 1976) assumed that the estuarine circulation was directly driven by longitudinal buoyancy gradients. According to Hansen and Rattray (1965), such gravitational circulation scales with ∂_{x}bH^{3}/(48*A _{υ}*) for a flat bottom estuary, with the horizontal buoyancy gradient ∂

*, the mean water depth*

_{x}b*H*, and the mean eddy diffusivity

*A*. The latter quantity is typically used as a tuning parameter to fit observed estuarine circulation to this theory.

_{υ}The assumption of the dominant contribution of gravitational circulation to estuarine circulation has been challenged by a number of studies during the last two decades. Simpson et al. (1990) first described the process of strain-induced periodic stratification (SIPS), where flood tide straining of the density field destabilizes the water column and vice versa on ebb tide, resulting in substantially enhanced small-scale turbulence and vertical mixing during flood and suppressed vertical mixing during ebb. Jay and Musiak (1994) showed that this internal tidal asymmetry transports more momentum down to the near-bed region during flood than during ebb, thus generating residual transports in support of estuarine circulation and in addition to gravitational circulation. They suggested that this internal tidal asymmetry (which we refer to here as “tidal straining”) is a major mechanism to create near-bed residual flow convergence near the tip of the salt intrusion, such that generation of estuarine turbidity maxima (ETMs) is triggered. In a numerical model study, Burchard and Baumert (1998) underlined the importance of tidal straining for generation of estuarine circulation by demonstrating that switching off tidal straining did reduce estuarine circulation (and thus ETM generation) more substantially than elimination of gravitational circulation (by switching off the internal pressure gradient forcing). This is in accordance with observations by Stacey et al. (2001, 2008), who show that the creation of residual flow profiles is dominated by pulses because of interaction of shear, stratification, and mixing (as opposed to steady gravitational circulation forcing).

The strong vertical mixing asymmetry between ebb and flood resulting from tidal straining has been shown in a number of field studies including turbulence observations (Peters 1999; Peters and Bokhorst 2000; Rippeth et al. 2001; Simpson et al. 2005) and numerical model studies (Simpson et al. 2002; Baumert and Peters 2007; Burchard et al. 2008; Burchard 2009; Verspecht et al. 2009b). Despite this evidence, to better compare to classical theory or to allow for analytical solutions, many model studies of estuarine circulation are based on the assumption of constant eddy viscosity and eddy diffusivity, thus eliminating the effect of tidal straining (Lerczak and Geyer 2004; Hetland and Geyer 2004; Huijts et al. 2006; Valle-Levinson 2008). Recently, Stacey et al. (2009, manuscript submitted to *J. Phys. Oceanogr.*) have shown that this assumption of a tidal-mean effective eddy viscosity is an oversimplification, because the tidal-mean momentum flux and the tidal-mean shear may have the same sign in large parts of the water column, leading there to a negative effective eddy viscosity.

Various other processes have been identified to generate residual flows in estuaries, such as wind straining (Scully et al. 2005; Ralston et al. 2008; Verspecht et al. 2009b; Chen and Sanford 2009), flow curvature (Chant 2002), rotation of the earth (Lerczak and Geyer 2004; Huijts et al. 2006; Valle-Levinson 2008), and lateral advection (Lerczak and Geyer 2004; Cheng and Valle-Levinson 2009; Scully et al. 2009).

However, independently of wind forcing, bathymetry, and latitude, gravitational circulation and tidal straining are ubiquitous processes in estuaries. So far, the relative contribution of these two processes to estuarine circulation has not been quantified. This is surely due to the fact that both shape residual flow profiles in a similar manner. Direct quantification of these processes is difficult, because it would require accurate observations of time-dependent eddy viscosity profiles, which is difficult for estuarine flows (Lu and Lueck 1999; Simpson et al. 2005). A further problem in determining the relative importance of gravitational circulation and tidal straining to estuarine circulation is the large variety of estuaries and within each estuary the high variability of forcing conditions (spring–neap cycle and variations in runoff and wind forcing). Therefore, we follow here the approach of Burchard (2009), who showed that the flow characteristics of periodic estuarine and coastal circulation with constant wind forcing and horizontal buoyancy gradients depend on only a few nondimensional parameters. These parameters are mainly the horizontal Richardson number, the inverse Strouhal number, the inverse Ekman number, the relative contribution of wind to the velocity scale, and the wind direction. With this, one-dimensional water column models for estuarine circulation can be applied to explore this relevant parameter space with high resolution.

In the present study, we derive analytical expressions for all contributions to residual velocity profiles, including tidal straining, gravitational circulation, wind straining, and residual runoff. These are evaluated by means of a numerical water column model.

We have restricted our attention to periodically stratified estuaries, neglecting the earth’s rotation. This approach, then, ignores the case in which the estuarine bottom boundary layer does not penetrate through the water column, and the outflowing layer may be considered to be inviscid (Geyer et al. 2000; Stacey and Ralston 2005; Ralston et al. 2008). In this paper, we consider only flows that have reached a periodic state, which is in contradiction to the assumption of permanent stratification at constant horizontal buoyancy gradient forcing (Blaise and Deleersnijder 2008).

There are many tidally energetic estuaries and coastal areas for which periodic stratification, including full mixing during a certain phase of the tide, is typical. Several examples for such estuaries are the Conwy Estuary in North Wales (Nunes and Simpson 1985; Simpson et al. 2001), the Chesapeake Bay tributary York River during spring tides (Simpson et al. 2005), and the Skagit River in Washington State (Yang and Khangaonkar 2009). Examples for coastal areas are tidal gulleys in the Wadden Sea of the German Bight (Buijsman and Ridderinkhof 2008; Burchard et al. 2008) or in lagoon systems such as Willapa Bay in Washington State (Banas and Hickey 2005). Also in coastal embayments, periodic stratification is typical, such as in Liverpool Bay (Simpson et al. 1990; Sharples and Simpson 1995; Verspecht et al. 2009a), for which SIPS has been defined.

The paper is organized as follows: First, the theory is presented in section 2, including the presentation of the basic equations (section 2a) and their nondimensional formulation (section 2b) and the derivation of residual profiles (section 2c) and their idealized analytical solutions (section 2d). In section 3, the model simulations are described, including the model setup (section 3a) and the results for basic simulations and sensitivity studies without wind straining and residual runoff (section 3b), as well as with wind straining or residual runoff (section 3c). The results are then discussed in section 4, and major conclusions are summarized in section 5.

## 2. Theory

### a. Dynamical equations

The one-dimensional dynamic equations for rectilinear tidal motion subject to a longitudinal buoyancy gradient constant in time and space can be formulated as follows:

In these dynamic equations, *t* denotes time; *x* and −*H* ≤ *z* ≤ 0 denote along-estuary and upward Cartesian coordinates, with the constant water depth *H*; and *u* denotes the longitudinal velocity component. Buoyancy is denoted by *b*, with

where *ρ* is potential density, *ρ*_{0} = 1000 kg m^{−3} is reference density, and *g* = 9.81 m s^{−2} is the gravitational acceleration. The equations are similar to those used for the study by Burchard (2009), with the difference that the earth’s rotation is neglected here.

In Eq. (1), the first term on the right-hand side represents the effect of the internal pressure gradient. The second term on the right-hand side of (1), *p _{x}*(

*t*), is a periodic external pressure gradient function with period

*T*chosen such that the vertical mean velocity

*U*(

*t*) is calculated as

with the vertical mean velocity amplitude *U*_{max} and the residual velocity *U _{r}*. This construction guarantees that the tidal-mean vertically integrated transport equals

*U*(for details, see Burchard 1999). Any other periodic tidal velocity variation could be obtained by using this method. Sea le vel differences, which cause large parts of the external pressure gradients in tidal flows, are neglected in the present study.

_{r}The term *A _{υ}* is eddy viscosity and

*K*is eddy diffusivity, both calculated by means of a two-equation turbulence closure model (for details, see the appendix). Momentum boundary conditions at the bottom are no slip with

_{υ}*u*(−

*H*,

*t*) = 0. At the surface a constant momentum flux

*τ*= −

_{x}^{s}*A*∂

_{υ}*for*

_{z}u*z*= 0 is prescribed. The resulting bed momentum flux is denoted as

*τ*= −

_{x}^{b}*A*∂

_{υ}*for*

_{z}u*z*= −

*H*.

### b. Nondimensional form of equations

To systematically investigate a large parameter space, the dynamic Eqs. (1), (2), (A1), and (A2) are transformed into nondimensional form, following an approach by Burchard (2009). Nondimensional quantities are the time *t̃* = *t*/*T*, the vertical coordinate *z̃* = *z*/*H*, the longitudinal velocity component *ũ* = *u*/*U*^{wt}, the buoyancy *b̃* = *b*/(∂* _{x}bH*), the turbulent kinetic energy

*k̃*=

*k*/(

*U*

^{wt})

^{2}, the dissipation rate ε̃ = ɛ/[(

*U*

^{wt})

^{3}/

*H*], the eddy viscosity

*Ã*=

_{υ}*A*/(

_{υ}*U*

^{wt}

*H*), and the eddy diffusivity

*k̃*=

_{υ}*K*/(

_{υ}*U*

^{wt}

*H*), with the constant water depth

*H*; the constant horizontal buoyancy gradient ∂

*; the tidal period*

_{x}b*T*; and a velocity scale that is composed of tide- and wind-induced current speed scales,

with the 10-m wind speed *W* and the reference air density *ρ _{a}*. The relative contribution of the wind to the velocity scale is denoted as

For rectilinear tides under irrotational conditions, the dynamics only depend on the modified inverse Strouhal number,

the modified horizontal Richardson number with respect to the *x* direction,

the relative contribution of the wind to the velocity scale *r _{u}*; the wind direction (up estuary or down estuary); and the nondimensional discharge

*ũ*=

_{r}*U*/

_{r}*U*

^{wt}. A weak dependence on the nondimensional bed and surface roughness lengths,

*z̃*

_{0}

^{b}=

*z*

_{0}

^{b}/

*H*and

*z̃*

_{0}

^{s}=

*z*

_{0}

^{s}/

*H*, respectively, is also given (see Burchard 2009).

Here, the superscript wt denotes that the velocity scale is composed of contributions from wind w and tide t and that, therefore, the inverse Strouhal number and the horizontal Richardson number are modified in the sense that they are based on the wind–tide velocity scale *U*^{wt} instead of only the tidal velocity scale *U*_{max}. It should be further noted that here the current velocity scale and not the friction velocity scale is used to allow for easier comparison of this theory to field observations (see also Burchard 2009).

### c. Residual profiles

Here, the longitudinal residual velocity profile is derived for a periodic solution of (1) with *u*(*z*, *t*) = *u*(*z*, *t* + *T*) for all *z* and *t*. First, (1) is tidally averaged and then vertically integrated from the vertical coordinate *z* to the surface at *z* = 0,

All tidal averages are notated as

After decomposition of the eddy viscosity and the longitudinal velocity into a tidal mean and a fluctuating part with

such that

(see also Stacey et al. 2009, manuscript submitted to *J. Phys. Oceanogr.*), Eq. (9) is transformed to the following form:

After division by 〈*A _{υ}*〉 and integration from

*z*= −

*H*to

*z*, a diagnostic equation for the tidal residual longitudinal velocity is finally obtained under consideration of the no-slip bottom boundary condition for 〈

*u*〉,

For simplicity, we introduce the following variables:

such that (14) can be written as

the tidal-mean barotropic pressure gradient can be eliminated,

with the nondimensional residual runoff velocity profile,

For clarity, we denote 〈*u _{t}*〉 = 〈

*u*

_{1}〉 (tidal straining), 〈

*u*〉 = 〈

_{g}*u*

_{2}〉 (gravitational circulation), 〈

*u*〉 = 〈

_{w}*u*

_{3}〉 (wind straining), and 〈

*u*〉 =

_{r}*γ*(

*z*)

*U*(residual runoff profile). With (18), the residual longitudinal velocity profile is expressed as the sum of four tidal-mean velocity profiles with the following properties:

_{r}The first term on the right-hand side 〈

*u*〉 contains the correlation between eddy viscosity and vertical shear and thus represents the tidal straining part of 〈_{t}*u*〉. This term vanishes when the correlation between eddy viscosity and vertical shear is zero: for example, for the case of constant eddy viscosity.The second term on the right-hand side 〈

*u*〉 is proportional to the horizontal buoyancy gradient and therefore represents the gravitational circulation part of 〈_{g}*u*〉.The third term on the right-hand side 〈

*u*〉 represents wind-induced straining._{w}The last term on the right-hand side 〈

*u*〉 is proportional to the vertically averaged tidal-mean velocity and therefore represents the contribution to 〈_{r}*u*〉 from residual runoff.At the bottom (

*z*= −*H*), all four residual profiles individually fulfill a no-slip boundary condition with 〈*u*〉 = 〈_{t}*u*〉 = 〈_{g}*u*〉 = 〈_{w}*u*〉 = 0._{r}The vertical mean of the first three profiles is individually zero, whereas the vertical mean of the last profile equals the vertically averaged tidal-mean velocity.

With this, the individual contributions of tidal straining, gravitational circulation, wind-induced circulation, and vertically averaged mean transport to the residual longitudinal velocity profile are well formulated and can thus be analyzed from numerical models.

### d. Idealized analytical solutions

When assuming constant viscosity and no rotation, the following polynomial solution is obtained analytically from (18):

which is exactly the solution recently presented by Ralston et al. (2008) as an extension of the Wong (1994) solutions for the separate effects of along-estuary wind stress and horizontal density gradients and the MacCready (2004) solution for the effects of horizontal density gradient and mean runoff (note that our solution represents the residual velocity profiles, whereas their solution calculates the deviation of the residual velocity from the depth-averaged velocity). A more realistic parabolic and constant-in-time eddy viscosity profiles is given by

with the bottom friction velocity *u*_{*}* ^{b}*, the bottom roughness length

*z*

_{0}

*, and the van Kármán number*

^{b}*κ*= 0.4. With this, an analytical solution for the residual flow profile can be obtained for zero surface momentum flux,

with the integration constant

In (22), the contribution from the mean transport (last term on right-hand side) is, as expected for the parabolic eddy viscosity profile, the classical logarithmic law of the wall. A good estimate for the tidal-mean absolute value of the bottom friction 〈|*u*_{*}* ^{b}*|〉 can be obtained by assuming that the residual velocity is much smaller than the tidal velocity amplitude (

*U*≪

_{r}*U*

_{max}) and that the flow is always instantaneously following the law of the wall,

such that, with (4),

is obtained. With this, the analytical solution (22) can be cast into nondimensional form,

showing that the residual profiles resulting from gravitational circulation scale directly with the horizontal Richardson number Ri* _{x}^{t}* (tide only), with logarithmic dependence on the bed roughness (see Fig. 1 for a graphical representation). This is in accordance with the scaling derived by Geyer et al. (2000); see also Scully et al. (2009).

Interestingly, for

the residual velocity profile is linear with

## 3. Model experiments

### a. Model setup

Experiments are carried out with the numerical model described in section 2a and the turbulence closure model presented in the appendix. To cover a wide range of possible scenarios, the following parameter space is explored here: Ri_{x}^{wt} = 0, … ,1 × 10^{−3} and St_{i}^{wt} = 1 × 10^{−4}, … , 1.1 × 10^{−3}, with a resolution of 41 × 41 simulations. Variations of Ri_{x}^{wt} and St_{i}^{wt} are achieved by systematic combinations of *U*^{wt} and ∂* _{x}b* (for details, see Burchard 2009). The dimensionless bed roughness length and surface roughness length are held at constant values of

*z̃*

_{0}

^{b}= 5 × 10

^{−5}and

*z̃*

_{0}

^{s}= 1 × 10

^{−5}, respectively. The upper limit for Ri

_{x}

^{wt}is given here by the transition to permanently stratified regimes that cannot be treated with the assumption of constant ∂

*. It should be noted that for simulations including the earth’s rotation, the limit for permanent stratification is generally much higher (for details, see Burchard 2009). For an*

_{x}b*M*

_{2}tide with

*T*= 44 714 s, the upper limit for St

_{i}

^{wt}is equivalent to

*U*

^{wt}>

*H*/(50 s) (e.g., 0.2 m s

^{−1}for

*H*= 10 m), which means that we do not consider weak tidal flow here (which would most likely lead to permanent stratification). All simulations were carried out with a time step of Δ

*t*=

*T*/1000 and 100 equidistant vertical layers, thus providing a resolution high enough to avoid significant discretization errors [see Burchard et al. (2005) and Umlauf et al. (2005) for details of the discretization schemes].

The following measure for the intensity of the residual current is used here:

calculating the average absolute value of a residual velocity profile.

For solving the integrals in (18) numerically, 〈*A*′* _{υ}*∂

*〉 and 〈*

_{z}u*A*〉, which are defined at the cell interfaces, are first linearly interpolated between the interfaces for each grid cell and then the resulting piecewise functions are solved analytically. The values of 〈

_{υ}*A*〉 and 〈

_{υ}*A*′

*∂*

_{υ}*′〉 at the bottom are calculated according to boundary conditions derived from the law of the wall:*

_{z}u### b. Simulations without wind straining and residual runoff

Figure 2 shows various parameters during a full tidal cycle for a reference case with Ri_{x}^{wt} = 5 × 10^{−4} and S_{i}^{wt} = 5 × 10^{−4}. Although this is a case with weak stratification, the classical SIPS dynamics can be clearly identified by stable stratification after ebb flow (*t̃* = 1) and fully mixed conditions after flood (*t̃* = 0.5), by stronger eddy viscosity *Ã _{υ}* during flood than during ebb and largely positive buoyancy production

*b̃*during flood. For increasing buoyancy gradient forcing, the duration of the fully mixed phase (nonnegative

*b̃*) will be shorter (see, e.g., Fig. 2 of Burchard 2009) until the buoyancy gradient forcing exceeds a threshold value beyond which stable stratification persist permanently, thus not allowing for periodic states within the current model framework.

For the basic situation without wind forcing and residual runoff, Fig. 3 shows nondimensional residual velocity profiles 〈*ũ*〉 and their contributions from tidal straining 〈*ũ _{t}*〉 and gravitational circulation 〈

*ũ*〉 for nine different combinations of the inverse Strouhal number St

_{g}_{i}

^{wt}and the horizontal Richardson number Ri

_{x}

^{wt}. As expected, all profiles are directed up estuary in the lower half of the water column and down estuary in the upper half of the water column; that is, tidal straining and gravitational circulation enhance each other and add up to the total estuarine circulation. Estuarine circulation intensifies with increasing horizontal Richardson number (because of the increased buoyancy gradient forcing) and is weakly reduced with increasing inverse Strouhal number. The latter can be explained in the following way: An increased St

_{i}

^{wt}can be obtained with an increased tidal frequency (decreased tidal period

*T*), which in turn will not provide enough time for the tidal asymmetry to sufficiently develop, with the result that the residual circulation is weakened. All profiles in Fig. 3 clearly show that the residual circulation resulting from tidal straining is about a factor of 2 stronger than the gravitational circulation, thus contributing about two-thirds of the estuarine circulation.

This is accurately quantified in the top-right panel of Fig. 4, where the relative contribution of tidal straining to the residual profile is shown for the whole parameter space. For large parts of the space, this contribution lies between 0.6 and 0.7 (i.e., around two-thirds). Only near the transition to permanent stratification (i.e., for high Ri_{x}^{wt}), gravitational circulation is more important. The top-left panel shows that the estuarine circulation is indeed more strongly and about linearly depending on Ri_{x}^{wt} and weakly decreasing with St_{i}^{wt}. The two bottom panels in Fig. 4 show the contributions of tidal straining and gravitational circulation separately, suggesting that the gravitational circulation is almost invariant with St_{i}^{wt}. In contrast to this, the tidal straining contribution does significantly decrease with increasing St_{i}^{wt}, supporting the argument that increased tidal frequency reduces tidal asymmetry.

A quantitative scaling analysis (see Fig. 5) shows that, in the limit of weak buoyancy gradient forcing, estuarine circulation and its contributions from tidal straining circulation and gravitational circulation all scale linearly with Ri_{x}^{wt}. For larger buoyancy gradients, dependence on Ri_{x}^{wt} becomes stronger and nonlinear, as suggested by Fig. 4. Gravitational circulation is largely independent of St_{i}^{wt}, which is in good agreement with the analytical solution presented in Eq. (26) for constant-in-time parabolic eddy viscosity profiles. This analytical solution for gravitational circulation is a good approximation for weak buoyancy gradient forcing because for strong buoyancy gradient forcing nonlinearities are expected because of the dependence of the eddy viscosity on stratification. Figure 5 shows further that the dependence of estuarine circulation and tidal straining circulation on St_{i}^{wt} is weakly negative, with best-fit exponents of −0.2 and −0.28, respectively.

In the following, results from a number of model experiments are discussed. These results shed more light into the generation mechanism of estuarine circulation and also give an idea about the sensitivity of the model results to the turbulence closure model.

To estimate the importance of tidal straining circulation and gravitational circulation for estuarine circulation (Burchard and Baumert 1998) using a two-dimensional along-estuary vertically resolving model for exploring generation mechanisms of estuarine turbidity maxima switched off tidal straining in one experiment by setting *N* ^{2} to zero for all turbulence closure calculations and switched off gravitational circulation in another experiment by setting the internal pressure gradient to zero in the momentum equation. The same is done here, and the resulting estuarine circulation intensity is shown in Fig. 6. When tidal straining is technically removed (left panel), estuarine circulation is substantially decreased and varies linearly with Ri_{x}^{wt}. The remaining contribution of tidal straining falls below 0.3 but does not vanish completely because turbulence still contains some remaining tidal asymmetry resulting from asymmetric shear (and thus shear production) caused by gravitational circulation. Compared to the full model (bottom-right panel in Fig. 4), gravitational circulation is weaker here, especially for high Ri_{x}^{wt}, because tidally averaged eddy viscosity (which occurs in the denominator of 〈*ũ _{g}*〉) is not damped by stable stratification generated by tidal straining. This is also the reason why permanent stratification does not occur in the parameter space shown. These results compare well with the analytical gravitational circulation solution (26). Looking up the relative intensity of residual circulation for the same roughness length,

*z̃*

_{0}

^{b}= 5 × 10

^{−5}, results in a higher circulation intensity for the analytical solution (see Fig. 1). This can be explained by the constant-in-time eddy viscosity used in the analytical problem, which is consistent with fully developed currents but overestimates the low mixing during slack tides resulting for the eddy viscosity calculated by the turbulence closure model. Stacey et al. (2008) show that during these slack tide periods of low vertical mixing, pulses of residual circulation are triggered.

The dependence of estuarine circulation and the contribution from tidal straining on specifications in the turbulence closure model is shown in Fig. 7. A decreased steady-state Richardson number, resulting in decreased mixing efficiency (see appendix for details), slightly increases the relative importance of tidal straining and vice versa. This means that stronger turbulence damping by stable stratification (smaller steady-state Richardson number and smaller mixing efficiency) increases tidal asymmetry and thus tidal straining.

The role of convective mixing during parts of the flood (for details, see Burchard 2009) is highlighted by two experiments with modified parameterizations for convective mixing. In the first experiment, the generic choice of *c*_{3}^{+} = 1 (see appendix for details) is replaced by *c*_{3}^{+} = *c*_{3}^{−}, reducing the dissipation rate for convective mixing by converting the buoyancy forcing in the dissipation rate Eq. (A2) to a sink term and thus enhancing convective mixing. This has the effect that the tidal straining is more pronounced (see Fig. 8, top left). In the second experiment, convective mixing is reduced by setting all instable values of *N* ^{2} to zero for the turbulence closure model, thus assuming that convection does not enhance vertical mixing. With this, the tidal straining contribution to estuarine circulation is significantly reduced (see Fig. 8, top right), but it is still larger than 0.5 in most parts of the considered parameter space. This means that the intensity of convective mixing strongly influences the tidal straining component of the residual estuarine circulation. However, the strength of the estuarine circulation itself is not strongly affected by the different compositions of estuarine circulation. This can be explained by a negative feedback process: stronger convective mixing during flood increases the correlation between vertical shear and eddy viscosity and thus tidal straining, but it also decreases gravitational circulation by increasing tidal-mean viscosity [cf. 𝒜_{1} and 𝒜_{2} in Eq. (15)]. Turbulence closure models are generally not well calibrated for such sheared convective regimes, partially because of lack of data and partially because of underrepresented physical processes such as countergradient fluxes.

### c. Simulations with wind straining or residual runoff

Effects of wind straining and residual runoff on estuarine circulation can easily be assessed with the analytical method developed here. The effect of residual runoff on shaping residual velocity profiles is shown in Fig. 9 for a horizontal Richardson number of Ri_{x}^{wt} = 5.0 × 10^{−4} and an inverse Strouhal number of St_{i}^{wt} = 5.0 × 10^{−4}. With the depth-mean residual velocity *U _{r}* ranging between 1% (

*ũ*=

_{r}*U*/

_{r}*U*

^{wt}= −0.01) and 10% (

*ũ*= −0.1) of the velocity scale, significant variations of the contributions from tidal straining

_{r}*ũ*and gravitational circulation

_{t}*ũ*are seen. For a small residual runoff of 1% and 2%, near-bed residual flow is directed up estuary and 〈

_{g}*ũ*〉 is still typically a factor of 2 larger than 〈

_{t}*ũ*〉. However, with increasing runoff, the contribution from tidal straining decreases. This can be explained by the fact that the runoff itself adds extra turbulence production because of the increased bottom stress, which in turn decreases the tidal asymmetry. This increasing importance of gravitational circulation is quantified in Fig. 10, which shows that the relative contribution of tidal straining decreases with increasing runoff for large part of the considered parameter space. A near-linear residual velocity profiles is given for example in the bottom-left panel of Fig. 9 for

_{g}*ũ*= −0.05, in a similar way as predicted by the analytical solution for gravitational circulation with a balance between runoff and horizontal Richardson number [see Eq. (28)].

_{r}An interesting effect can be studied when wind straining is applied instead of residual runoff (see Fig. 11), where the impact of weak (*r _{u}* = 0.1) and moderate (

*r*= 0.2) up-estuary and down-estuary winds on residual profiles is shown for a horizontal Richardson number of Ri

_{u}_{x}

^{wt}= 2.0 × 10

^{−4}and an inverse Strouhal number of St

_{i}

^{wt}= 5.0 × 10

^{−4}. As expected, estuarine circulation is enhanced by down-estuary winds and reduced by up-estuary winds. Down-estuary wind stabilizes the water column such that, also during flood, convective overturning is suppressed, with the consequence that tidal straining is significantly reduced for scenarios with moderate down-estuary wind; thus, gravitational circulation becomes generally more important than tidal straining. Near-bed tidal straining may even be reversed by down-estuary wind, as clearly seen in the top-right panel of Fig. 11. This effect will be strongest for small Ri

_{x}

^{wt}and can be explained as follows: During ebb, down-estuary wind straining increases vertical shear and thus increases turbulence production, whereas during flood vertical shear and thus turbulence are reduced. This generates an inverse tidal asymmetry also with zero or small horizontal buoyancy gradients, slightly weakening the joint effect of wind straining and gravitational circulation.

The extent of inverse tidal straining is highlighted in Fig. 12. For relatively strong down-estuary wind straining (*r _{u}* = 0.2), for almost the entire parameter space inverse straining acts such that the near-bed residual velocity resulting from tidal straining

*ũ*is negative and thus reversed. In this situation, for most of the parameter space, gravitational circulation is stronger than tidal straining circulation. The opposite happens for up-estuary wind: for the entire parameter space under investigation, tidal straining is significantly stronger than gravitational circulation, because tidal mixing asymmetry resulting from wind straining is supporting tidal mixing asymmetry resulting from tidal straining.

_{t}## 4. Discussion

For a periodically stratified situation with negligible residual runoff and no wind forcing, the contribution of tidal straining to estuarine circulation is significantly stronger than gravitational circulation. This has been shown for a variety of parameters in the turbulence closure model and is therefore a robust result.

Also, for various wind forcing and runoff situations, tidal straining is generally more important than gravitational circulation, with a few exceptions for down-estuary winds, when convective overturning during flood tide and thus tidal asymmetry is decreased (section 3c). The remaining question is, for which suite of realistic estuarine and coastal situations is this result representative?

In some estuaries, the water column is capped by a permanently stratified layer such that the present result is not directly transferable. However, it may be hypothesized that, in the periodically stratified lower part of the water column in such estuaries, the situation is similar because also there the stable stratification is eroded during flood and restratified during ebb. This would be in accordance with the simple linear stress model developed by Stacey and Ralston (2005) and needs to be investigated in a further study in which the assumption of a constant horizontal buoyancy gradient is relaxed. Such a study could be done in a similar way as suggested by Blaise and Deleersnijder (2008); that is, expressing ∂* _{x}b* as a function of

*b*, with the consequence that the constant ∂

*would need to be replaced by the tidally averaged vertical profile of the horizontal buoyancy gradient in Eq. (15).*

_{x}bLateral advection is another process that had to be neglected in this one-dimensional parameter space study. It has already been shown by Nunes and Simpson (1985) that tidally energetic estuaries may develop a characteristic transverse circulation with surface convergence during flood and bottom convergence during ebb. Lerczak and Geyer (2004) were the first to show that this process leads to enhancement of classical estuarine circulation. Scully et al. (2009) could show that also this process scales with the horizontal Richardson number and may be of the same order as gravitational circulation. We hypothesize here that tidal straining, which has been neglected in the studies by Lerczak and Geyer (2004) and Scully et al. (2009), is important also in estuaries of limited width and must not be neglected as a process contributing to estuarine circulation.

The effect of the earth’s rotation may be relatively small for narrow estuaries. In wide estuaries and coastal regions, however, rotation has the effect of destabilizing the water column, such that periodic stratification may be seen for much larger horizontal Richardson numbers as compared to the case without rotation. Some of these dynamics are also reproduced by one-dimensional models, as discussed in detail by Burchard (2009). The present theory can easily be extended for rotational one-dimensional flows. In principle, extensions to two-dimensional (laterally homogeneous or longitudinally homogeneous) estuaries should also be possible, as well as extensions to fully three-dimensional estuaries, as long as strictly periodic solutions are analyzed. For this, the analysis would have to be extended to all terms in the longitudinal momentum equation, including the advective terms, which may play an important role (Li and O’Donnell 2005; Scully et al. 2009). Such an analysis would distinguish between the contributions of the relevant processes to estuarine circulation much more accurately than the method applied by Burchard and Baumert (1998), who switched off single processes in different numerical experiments, thus providing an invasive numerical analysis.

The question for the relative importance of tidal straining and gravitational circulation to estuarine circulation is of more interest than theoretical alone. Specifically for coastal flows with small horizontal buoyancy gradients, estimates only based on scalings for gravitational circulation may lead to a strong underestimation of estuarine circulation. Postma (1954) and Zimmerman (1976), for example, conclude that buoyancy gradients in the Marsdiep in the Dutch Wadden Sea are not sufficiently strong to drive net transport of suspended particulate matter. This is in contrast to recent model results by Burchard et al. (2008) showing that even weak buoyancy gradients can drive net SPM transport into the Wadden Sea.

The reason why the role of tidal straining for estuarine circulation is generally underestimated is that its contribution to residual profiles has a shape similar to that of gravitational circulation and that it even obeys a similar scaling. As shown in Fig. 5, estuarine circulation, tidal straining, and gravitational circulation scale linearly with the horizontal Richardson number Ri* _{x}* for the limit of weak buoyancy gradient forcing, a result that has been suggested for gravitational circulation already by Scully et al. (2009); see also Geyer et al. (2000). For estuarine circulation and tidal straining, additional weak negative dependence on the inverse Strouhal number is given. It would be a great (maybe impossible) challenge for future studies to directly observe these subtle differences in the field. A more promising strategy for supporting the present model results by field observations would be to obtain a close agreement between model results and field observations for a variety of parameters including turbulence and stress measurements and then to analyze the individual contributions to residual circulation from the model results.

Although in the past turbulence closure models have been assessed extensively for their performance in stably stratified shear flow, the level of predictability of these models for sheared convective flow is low. This may be due to the fact that strong convective mixing would approximately homogenize all tracers in the convective boundary layer, such that even crude methods such as convective adjustment lead to acceptable results. However, shown here, the rate of convective mixing of momentum is decisive for shaping residual circulation (see Fig. 8). Specifically, the relative importance of tidal straining increases with increased convective mixing during flood tide. This is because the no-slip bottom boundary condition for velocity prescribes zero-bed velocity and increased convective mixing increases near-bed shear and thus residual circulation resulting from tidal straining. It is a future challenge to turbulence closure modeling to refine parameterizations of sheared convective flow. Such an improved parameterization should fit into the general framework of two-equation turbulence closure models, and a reliable fit for the buoyancy parameter *c*_{3}^{+} in the dissipation rate Eq. (A2) or any other length scale equation would already be a great progress.

## 5. Conclusions

This study provides mathematical definitions of tidal straining, gravitational circulation, wind straining, and residual runoff contributions to estuarine circulation for periodic tidal flows. These are derived analytically from the tidally averaged one-dimensional momentum equations.

The major result of this study is that tidal straining is the dominant process in generating estuarine circulation for periodically stratified estuaries. As a rule of thumb, for situations without wind straining and residual runoff, tidal straining contributes two-thirds and gravitational circulation contributes one-third of the estuarine circulation.

It is shown that tidal straining may not only be driven by tidal asymmetries resulting from horizontal buoyancy gradients but also from other asymmetric processes such as wind straining. For weak horizontal buoyancy gradients and down-estuary wind straining, tidal straining may be reversed; that is, it opposes estuarine circulation.

The intensity of estuarine circulation scales linearly with the horizontal Richardson number and weakly with the Strouhal number, where the latter contribution results from tidal straining. Extensions of the present theory to partially mixed and strongly stratified estuaries as well as to two and three dimensions, including additional terms such as advection and the earth’s rotation, are possible.

## Acknowledgments

The work of Hans Burchard has been supported by the EU-funded project ECOOP (European COastal-shelf sea OPerational monitoring and forecasting system Contract 36355). Robert Hetland thanks the Fulbright Commission for funding a sabbatical at the Leibniz Institute for Baltic Sea Research, which made this collaborative work possible. The authors are indebted to Elisabeth Schulz (Warnemünde) for carefully correcting the manuscript.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

### APPENDIX

#### Turbulence Closure Model

Because effects of tidal straining strongly depend on the turbulence damping by stable stratification during ebb and the turbulence enhancement by instable stratification during flood (i.e., the exchange between turbulent kinetic energy and potential energy), the turbulence closure model needs to be presented here in some detail. Eddy viscosity and diffusivity are denoted by *A _{υ}* and

*K*, respectively, and are calculated by means of a turbulent kinetic energy (per unit mass)

_{υ}*k*and its dissipation ɛ (Umlauf and Burchard 2005),

with the shear frequency squared, *S*^{2} = (∂* _{z}u*)

^{2}; the buoyancy frequency squared,

*N*

^{2}= ∂

*; and the empirical parameters*

_{z}b*σ*= 1,

_{k}*c*

_{1}= 1.44,

*c*

_{2}= 1.92, and

*σ*

_{ɛ}=

*κ*

^{2}/[(

*c*

_{μ}^{0})

^{2}(

*c*

_{2}−

*c*

_{1})], where

*κ*= 0.4 is the van Kármán number and

*c*

_{μ}^{0}= 0.527 is the quasi-equilibrium stability function for momentum. Boundary conditions for

*k*and ɛ are derived from the law of the wall (Umlauf and Burchard 2005).

The eddy viscosity and diffusivity are then calculated as

with the nondimensional stability functions *c _{μ}* and

*c*′

*, which are nonlinear functions of the shear parameter*

_{μ}*α*=

_{S}*S*

^{2}

*k*

^{2}/ɛ

^{2}and the buoyancy parameter

*α*=

_{N}*N*

^{2}

*k*

^{2}/ɛ

^{2}and include the entire algebraic second-moment closure. Here, the second-moment closure by Cheng et al. (2002) has been applied, which provides realistic and robust estimates of vertical mixing (see Umlauf and Burchard 2005).

For stable stratification, the parameter *c*_{3} = *c*_{3}^{−} regulates the balance between mixing and stratification and defines the steady-state Richardson number Ri_{st} at which turbulence in homogeneous shear layers may come to a steady state. According to Shih et al. (2000), the steady-state Richardson number should be Ri_{st} = 0.25 for fully developed turbulence. For the second-moment closure by Cheng et al. (2002), this is obtained by setting *c*_{3}^{−} = −0.74 where Ri_{st} is an increasing function of *c*_{3}^{−}. A steady-state Richardson number Ri_{st} > 0.25 would lead to weaker suppression of turbulence by stable stratification and vice versa (for details, see Umlauf and Burchard 2005).

For stable stratification, the role of *c*_{3}^{−} can also be discussed in the light of the mixing efficiency Γ = *K _{υ}N*

^{2}/ɛ, which has been estimated for stable stratification by Osborn (1980) as Γ ≈ 0.2. Assuming homogeneous shear layers (constant

*S*

^{2}and

*N*

^{2}and no vertical variation in

*k*, ɛ,

*A*, or

_{υ}*K*) in steady state, the

_{υ}*k*–ɛ model (A1) and (A2) reduces to

which can be combined to

With this, the mixing efficiency in steady state is independent of the second-moment closure used (see also Umlauf 2009). Therefore, Ri_{st} = 0.25 [equivalent to *c*_{3}^{−} = −0.74 when using the Cheng et al. (2002) second-moment closure] results in Γ = 0.22, a value close to the estimate by Osborn (1980). A higher value of Ri_{st} = 0.28 (equivalent to *c*_{3}^{−} = −0.61) would result in Γ = 0.234, and a lower value of Ri_{st} = 0.22 (equivalent to *c*_{3}^{−} = −0.92) would result in Γ = 0.2. Thus, higher steady-state Richardson numbers are consistent with higher mixing efficiencies and thus with stronger mixing and vice versa. This actually is a new result, which is relevant here and supports the theories on two-equation turbulence closure modeling developed by Burchard and Baumert (1995), Baumert and Peters (2000), Burchard and Bolding (2001), and Umlauf and Burchard (2005).

For instable stratification, *c*_{3} = *c*_{3}^{+} = 1 is generally used (see, e.g., Burchard and Bolding 2001), a choice that is made to provide a source term for the dissipation rate equation in the case of free convection with *P* = 0 and *B* > 0. However, Umlauf and Burchard (2005) pointed out that *c*_{3}^{+} > 1 is not necessary and that, for a free convection scenario, the resulting buoyancy fluxes are more realistic for *c*_{3}^{+} = *c*_{3}^{−}. Therefore, this alternative option is also tested in one of the simulations discussed in section 3b.

## Footnotes

*Corresponding author address:* Hans Burchard, Leibniz Institute for Baltic Sea Research Warnemünde, Seestraße 15, D-18119 Rostock, Germany. Email: hans.burchard@io-warnemuende.de