The aim of this one-dimensional water column study is to combine modified versions of three characteristic parameters for periodic tidal flow under the influence of a longitudinal buoyancy gradient—the horizontal Richardson number, the inverse Strouhal number, and the inverse Ekman number—into a parameter space study, including constant wind forcing from various directions. It is shown how the underlying dynamical equations can be cast into nondimensional form, depending mainly on these three nondimensional parameters plus the relative wind speed and the wind direction. Idealized model simulations are carried out for the whole realistic range of horizontal Richardson and inverse Strouhal numbers for various latitudes, showing the amplitude of the tidally induced stratification for a wide range of scenarios. It is found that classical threshold values for the horizontal Richardson number, indicating the switch from periodic stratification to permanent stratification, are valid only for special cases, and that this switch also strongly depends on the inverse Strouhal number, the inverse Ekman number, and the wind vector. The transverse residual flow is close to the thermal wind balance for a variety of parameters, and nondimensionalized longitudinal residual flow shows the classical estuarine exchange flow pattern with little variation in the near-bed onshore component. Wind straining is confirmed as an important estuarine and coastal process, enhancing estuarine circulation for offshore (down estuary) winds and vice versa. Agreement with field data from Liverpool Bay is good, including the explanation of a flood tide local maximum of the dissipation rate in the upper half of the water column. An equation for the second time derivative of the potential energy anomaly is derived for quantifying the dynamical processes leading to stratification due to the straining of the horizontal density gradient.
Tidal straining is the interaction between periodically fluctuating tidal flow and slowly varying horizontal density gradients (Simpson et al. 1990), that is, a superposition of an oscillating external pressure gradient with a unidirectional internal pressure gradient. The underlying horizontal density gradients are caused by nonzero net surface and lateral buoyancy fluxes into coastal waters. Differential effects due to bottom slopes convert spatially homogeneous surface buoyancy fluxes into horizontal buoyancy gradients, with lower densities toward the coast for net freshwater or heat fluxes into the water and vice versa. Classical estuaries and coastal areas with less dense coastal waters compared to offshore waters are located in humid midlatitude areas, but inverse estuaries with evaporation larger than precipitation are characteristic features in arid regions. For the ease of notation and without loss of generality, it is assumed for this study that the coast is to the east (positive x axis) and the ocean to the west, and that density is decreasing toward the coast (or up the estuary), that is, with buoyancy increasing toward the coast. Furthermore, the rotational effects are related to the Northern Hemisphere.
The major dynamic effect of tidal straining is that stratification is increased, and mixing reduced, when the currents are directed toward higher densities by shearing less dense water over more dense water. After a current reversal, when dense waters are sheared over less dense waters, stratification is eroded and mixing increased. Depending on the balance of the forces, the waters are permanently mixed, permanently stratified, or periodically stratified. The latter regime has been denoted as strain-induced periodic stratification (SIPS) by Simpson et al. (1990). An important effect of the periodic tidal straining process is its potential to enhance onshore (up estuary) near-bottom residual currents by mixing more momentum down to the near-bed region during flow into the direction of positive buoyancy gradients (typically flood) than in the opposing direction (Jay and Musiak 1994; Scully and Friedrichs 2007). The tidal straining contribution to the residual flow profiles adds to the contribution of gravitational circulation, and has been identified as a major cause for onshore suspended matter transport (Burchard et al. 2008) and the subsequent generation of estuarine turbidity maxima in estuaries (Jay and Musiak 1994; Burchard and Baumert 1998).
By using analytical flow profiles, Simpson et al. (1990) derived threshold values for the transition between the permanently mixed and the periodically mixed states, (with the horizontal buoyancy gradient, dxb; the water depth, H; and the tidal flow amplitude, umax), and between the periodically mixed state and the permanently stratified state, . Simpson et al. (1990) and Rippeth et al. (2001) could confirm these theoretically derived threshold values by field observations in Liverpool Bay. This critical parameter used by Simpson et al. (1990) differs from the horizontal Richardson number, Rx = (∂xbH2)/u*2, with the amplitude of the bed friction velocity, , and the drag coefficient, cd, defined by Monismith et al. (1996) only by using a tidal velocity amplitude instead of a friction velocity scale. According to Monismith et al. (1996), values of Rx larger than order of one would lead to mixing that is too weak for overcoming stratification by horizontal density gradients, a result that is in agreement with the threshold estimate given by Simpson et al. (1990), given a drag coefficient of cd ≈ 2 × 10−3. Stacey et al. (2001) succeeded in relating strong observed pulses of onshore near-bottom residual flow to high values of Rx. Recently, Simpson et al. (2005) interpreted Rx as a measure for the ratio of the tidal straining to the tidal mixing terms in the potential energy anomaly, ϕ, which is a suitable measure for the water column stability (see also Simpson et al. 1990; Burchard and Hofmeister 2008, and references therein).
Baumert and Radach (1992) identified the inverse Strouhal number, Si = H/(Tu*), that is, the ratio of the vertical mixing time scale, H/u*, and the tidal period, T, as another characteristic parameter for mixing associated with tidal flow. By demonstrating that apart from bottom and surface roughness lengths, Si is the only parameter defining the dynamics of the well-mixed irrotational pressure gradient driven tidal flow, they could show how the relative time lag of turbulent parameters with respect to the bed stress increases with Si. It should be noted that there are many definitions for the Strouhal number, which originally was defined by Strouhal (1878) as the ratio of the velocity at which a wire is moved through the air and the product of its thickness and the frequency of the generated sound. In a theoretical study of estuarine circulation, Ianniello (1977) showed that the inverse Strouhal number [his parameter d02, see his Eq. (5)] determines the shape of velocity profiles and residual currents in well-mixed tidal flow in closed estuaries. This concept has later been extended to include a salinity equation by Jay and Smith (1990). In shelf sea oceanography, the Strouhal number has been defined as the ratio of the clockwise or anticlockwise tidal velocity components to the product of the water depth with the sum or difference between the tidal and inertial frequencies [multiplied by π, see Prandle (1982)]. It should be noted that a horizontal Strouhal number (ratio of tidal velocity amplitude to tidal frequency and length of the embayment) has been used by Schuttelaars and de Swart (1999) for scaling the dynamics of homogeneous tidal flow in coastal embayments.
In one-dimensional studies of estuarine circulation, the earth’s rotation is generally neglected (Prandle 2004; Scully and Friedrichs 2007; Stacey et al. 2008), although it is evident from one- (Heaps 1972; Sharples and Simpson 1995) and two-dimensional longitudinally averaged studies that the planet’s rotation plays a major role in generating transverse circulation (Lerczak and Geyer 2004; Huijts et al. 2006; Valle-Levinson 2008). The inverse Ekman number, E = fH/u*, is the characteristic parameter quantifying the influence of rotation, with the Coriolis parameter, f. Here, Ei may be interpreted as the ratio of the water depth H to the Ekman depth u*/f or, similarly, to the inverse Strouhal number, as proportional to the ratio of the mixing time scale to the inertial time scale, 2π/f.
Parameterized forms of the potential energy anomaly consider energy inputs into the water column due to wind stress as an agent for mixing, without taking wind-straining effects into account (Wiles et al. 2006). However, observations of subtidal velocity profiles and stratification carried out by Scully et al. (2005) in the James River estuary revealed a correlation between wind stress and residual shear, which under conditions of strong stratification and wind stress easily overcompensated for the tidal straining effects. Offshore (down estuary) winds were shown to have a strong tendency to intensify estuarine circulation and stratification, which is a clear sign of the wind straining effects being dominant over the wind mixing effects. Onshore (up estuary) winds tended to reduce and partially even reverse the estuarine circulation and, consequently, had the potential to reduce stratification and intensify vertical mixing.
It is the aim of the present investigation to combine modified versions of these three characteristic parameters, Rx, Si, and Ei, including constant wind forcing from various directions for an extensive parameter space study of strain-induced periodic stratification dynamics, using a water column model including a well-tested second-moment turbulence closure scheme (Cheng et al. 2002; Umlauf and Burchard 2005). This type of turbulence closure model includes a number of physical processes that are neglected in many theoretical studies of estuarine circulation (e.g., Lerczak and Geyer 2004; MacCready 2004; Huijts et al. 2006), such as tidal straining and mixing asymmetries due to stable or instable stratification, as well as mixing time lag effects. The water column model simulates periodic tidal flows with sinusoidal tidal transport under the influence of a constant longitudinal buoyancy gradient and surface stress. With this, transverse circulation effects due to differential advection or the earth’s rotation (Lerczak and Geyer 2004), coastal upwelling or downwelling (Valle-Levinson 2008), salt intrusions from branching channels (Stacey et al. 2008), or flow curvature (Chant 2002) are excluded, restricting the results of this study to relatively wide straight estuaries (large Kelvin numbers, i.e., a large ratio between the width and internal Rossby radius) and offshore coastal areas outside several internal Rossby radii from the coast.
The interest in these physical processes in estuaries and coastal seas is motivated by the fact that all solute and particulate matter that is exchanged between the coast and the open sea has to pass through these regions being dominated by horizontal density gradients. Whereas solutes are drifting passively with the tidal residual flow, particulate matter characterized by the additional process of settling is subject to the tidal dynamics in a more complex way, with substantial impacts on coastal morphodynamics (see, e.g., Elias 2006). In addition, the migration of fish larvae and juvenile fish is strongly determined by tidal dynamics and residual transports (North and Houde 2003). Another important issue in coastal zones with substantial oxygen consumption due to the mineralization of organic matter and due to the growth and metabolism of benthic fauna is stratification persisting for extended periods, a process that may lead to oxygen depletion and subsequent extinction of benthic fauna.
In the following section, the one-dimensional dynamical equations are first introduced (section 2a) and then transformed into nondimensional notation (section 2b). A nondimensional form of the dynamical equation for the potential energy anomaly ϕ is derived in section 2c. Afterward, the results for the amplitude of tidal stratification are presented as functions of a modified horizontal Richardson number and a modified inverse Strouhal number for various latitudes, relative surface wind stresses, and wind directions, together with a few examples for residual flow profiles and resolved tidal dynamics for a standard scenario (section 3). In section 4, field data obtained in Liverpool Bay in 1999 by Rippeth et al. (2001) are cast into nondimensional form and compared to the idealized results for a similar standard scenario. In section 5, the results are discussed in the light of the second time derivative of the water column stability parameter, ϕ, including the straining term only. Finally, some conclusions from the results are drawn in section 6.
a. Dynamical equations
The one-dimensional dynamic equations for harmonic rectilinear tidal motion, subject to constant wind stress, zero surface buoyancy flux, and a constant in time and space longitudinal density gradient aligned with the tide, can be formulated as follows:
where t denotes time and x, y, and −H ≤ z ≤ 0 denote eastward, northward, and upward Cartesian coordinates, with the constant water depth H. Here, u and υ denote the eastward and northward velocity components, respectively. Buoyancy is denoted by b, with
where ρ is the potential density, ρ0 = 1000 kg m−3 is the reference density, and g = 9.81 m s−2 is the gravitational acceleration.
The earth’s rotation is taken into account by the Coriolis parameter f = 2ω sin(ϕ) with the angular frequency of the earth, ω = 7.29 × 10−5 s−1, and the latitude ϕ.
In Eq. (1), the first term on the right-hand side represents the effect of the internal pressure gradient on account of the horizontal buoyancy gradient. The second term on the right-hand side of (1), pgx(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 umax. This construction guarantees that the tidal mean vertically integrated transport is zero [see Burchard (1999) for details]. To ensure zero transverse volume transport for all times, a pressure gradient function pgy(t) is constructed accordingly. Any nonharmonic tidal longitudinal velocity variation could be obtained by using this method, including nonzero residual transports caused by river runoff. This is more efficient than prescribing surface slope functions and later calibrating a constant pressure gradient correction term as recently suggested by Blaise and Deleersnijder (2008).
Eddy viscosity and diffusivity are denoted by νt and ν′t, respectively, and are calculated by means of the turbulent kinetic energy (per unit mass), k, and its dissipation, ɛ:
with the dimensionless stability functions cμ and c′μ depending on the nondimensional shear and buoyancy parameters and including all details of the second-moment closure applied. It should be noted that the ratio of these two numbers, cμ/c′μ, equals the turbulent Prandtl number, which is thus an integral part of the turbulence closure. The turbulent kinetic energy is calculated by means of Eq. (4), with the constant turbulent Schmidt number, σk. The three terms on the right-hand side of the TKE equation are the shear production, P; the buoyancy production, B; and the dissipation rate, ɛ. For the dissipation rate ɛ, various closure equations have been suggested. Here, a dynamic equation for ɛ is used directly. It has been shown in several studies that the type of the length scale equation is of less importance than is the second-moment closure applied (Burchard et al. 1998; Umlauf and Burchard 2003), such that the dissipation rate equation is not shown here. All details about the turbulence closure models used here have been reported in detail by Umlauf and Burchard (2005). In the configuration described therein, the model has been successfully applied to reproduce a number of current, stratification, and turbulence observations in coastal seas, shelf seas, and lakes [see, e.g., Simpson et al. (2002) for tidal straining in Liverpool Bay; Burchard et al. (2002) for autumn stratification in the northern North Sea; Arneborg et al. (2007) for near-bottom saline inflow into Baltic Sea; Stips et al. (2002) for convection in Lago Maggiore; and Stips et al. (2005) for near-surface turbulence in Lake Neuchatel].
At the bottom boundary, a no-slip condition is imposed for the horizontal velocity components, and a no-flux condition is employed for the buoyancy and the turbulent kinetic energy. The surface momentum flux is given by the wind stress vector normalized by ρ0:
with the surface friction velocity,
and with the direction of the wind, α, with α = 0° corresponding to wind from north and α = 90° corresponding to wind from east. In (9), W is the constant wind speed, cdw is the surface drag coefficient with
[see Large and Pond (1981)], and ρa = 1.25 kg m−3 is the constant density of air.
The surface buoyancy flux is assumed to be negligible in this study. For the turbulent kinetic energy, k, no-flux boundary conditions are used for the surface and the bed, thus neglecting energy transfers from the wave field to the turbulence. Bottom and surface boundary conditions for the dissipation rate are given as
with the surface and bottom roughness lengths, z0s and z0b, respectively; the van Kármán parameter, κ = 0.4; and the absolute value of the bottom friction velocity, u*b. It should be noted that u*b results from the flow dynamics including a no-slip bottom boundary condition and is, thus, not an externally imposed parameter.
The surface roughness length is parameterized as follows:
with ac = 1400, a formula motivated by the work of Charnock (1955) [see Craig and Banner (1994)], connecting the surface roughness to the wind-driven wave state. For the bottom roughness length, the following formulation is used (see, e.g., Kagan 1995):
with the molecular kinematic viscosity ν and the height of the bottom roughness elements, h0b.
For the physical and numerical details of the bottom friction velocity calculation and the turbulence boundary conditions, see Burchard (2002).
b. Nondimensional form of dynamic equations
The underlying dynamical equations are now transformed to nondimensional equations by defining nondimensional state variables. The scales used for this transformation are the tidal period, T; the water depth, H; the horizontal buoyancy gradient, ∂xb; and a velocity scale, composed of tide- and wind-induced current speed scales:
The form of the wind contribution to uwt has been chosen such that it is consistent with the tidally induced velocity: as the surface friction velocity divided by the square root of the drag coefficient, compare to (9). In principle, any constant velocity could be chosen as the velocity scale. In many applications, friction velocities are preferred over current velocities as velocity scales (Baumert and Radach 1992; Stacey et al. 2001). From the practical point of view, however, the tidal velocity amplitude is easier to observe (Simpson et al. 1990) and also easier to impose in a one-dimensional model (Burchard 1999), such that for the present study this option, extended by a wind velocity scale, is used.
With these scalings, the following nondimensional state variables are obtained:
With the modified inverse Strouhal number (cf. with Baumert and Radach 1992),
the modified inverse Ekman number
where the second term on the right-hand side is the nondimensional buoyancy production, B̃. The nondimensional boundary conditions read as
for the momentum flux at the surface, where
denotes the relative contribution of the wind-induced currents to the velocity scale. The boundary conditions for the dissipation rate are of the following nondimensional form:
the horizontal modified Richardson number, ;
the modified inverse Ekman number, ;
the modified inverse Strouhal number, ;
the relative contribution of the wind-induced currents to the velocity scale, ru;
the direction of the wind relative to the buoyancy gradient, α;
the relative bottom roughness length, z̃0b; and
the relative surface roughness length, z̃0s;
Not all of these degrees of freedom are relevant in midlatitude coastal sea applications. For a fixed location, the ratio between the inverse Strouhal number and the Ekman number, , is constant, assuming dominance of one tidal constituent. It can furthermore be assumed that both the surface and the bottom roughness are of relatively minor importance for the dynamics of wind-, tide-, and buoyancy-driven coastal systems, such that the number of the relevant degrees of freedom is reduced to four. It should be noted that since both the surface and the bottom roughness lengths are parameterized instead of being directly imposed [see Eqs. (12) and (13)], the relevant nondimensional parameters for these roughness lengths would be (inverse Froude number squared), ν/(uwtH) (inverse Reynolds number), and h̃0b = h0b/H (relative height of the bottom roughness elements).
c. Potential energy anomaly
As a convenient measure for the stability of the water column, the potential energy anomaly ϕ has been defined by Simpson (1981) as the amount of mechanical energy (per m3) required to instantaneously homogenize the water column with a given density stratification. Within the notation of the present study, ϕ is of the following form:
with the depth mean buoyancy
Recently, a complete dynamical equation for ϕ has been derived by Burchard and Hofmeister (2008), which for the present idealized problem simplifies to
with the depth mean velocity 〈u〉.
In (28) the first term on the right-hand side represents straining and the second term represents vertical mixing. This equation will be used in section 5 for deriving the second derivative for ϕ in order to help us understand the processes leading to stratification.
The nondimensionalized potential energy anomaly is of the following form:
The stratification parameter ϕ̃ will be used in this study to quantify the degree of stratification for various flow situations.
3. Idealized model simulations
A large number of idealized nondimensional simulations are carried out to cover a wide range of possible coastal flow situations with interaction between tides, horizontal density gradients, and wind forcing.
For all simulations, the tidal period is fixed to reproduce the semidiurnal M2 tide with T = 44714 s, the water depth is kept to a constant value of H = 20 m, and the latitude is set to 30°, 40°, 50°, and 60°N, respectively. Various stationary wind-forcing situations are considered, with no wind, and the wind varying such that the relative contributions of the wind to the velocity scale, ru are 0.1, 0.2, and 0.4. The wind direction is varied between wind from the north (α = 0°), from the east (α = 90°, offshore, down estuary), from the south (α = 180°), and from the west (α = 270°, onshore, up estuary). Furthermore, the sensitivity of the model results to the nondimensional bottom roughness length, h̃0b = h0b/H, is studied for four different values of h̃0b. With this, only the buoyancy gradient, ∂xb, and the velocity scale, uwt, are varied continuously. All simulations presented here are run for 20 tidal periods and results are shown for the 20th tidal period only. The numerical resolution is high enough for excluding significant numerical effects, with a time step of Δt = 10−3 T and an equidistant vertical grid size of Δz = 10−2 H.
To test the dependence of the model results on the nondimensional numbers only, test simulations with identical ratios, H/uwt, but different water depths have been carried out for various forcing configurations. The results of these test simulations confirmed the soundness of the present nondimensional approach. With prescribed values for the horizontal Richardson number, , and the inverse Strouhal number, , these two values are varied as follows:
Figure 1 shows for midlatitudes (50°N) and no wind forcing the results of the maximum amplitude of the nondimensional potential energy anomaly, Δϕ̃, for a range of and , with a resolution of 41 × 41 simulations. Assuming a fixed water depth of H = 20 m and the M2 tidal period, typical values for the inverse Strouhal number would be between (resulting in umax ≈ 0.9 m s−1) and (resulting in umax ≈ 0.3 m s−1). For a fixed horizontal Richardson number , the dependence of tidal straining on is strong. However, for a fixed horizontal buoyancy gradient ∂xb, the variation of straining is small, a result of that is clearly seen from Fig. 1 when following the parabolic line for constant ∂xb = 10−6 s−2 [equivalent to a salinity gradient of about 1.46 × 10−4 (g kg−1) m−1]. Increasing for fixed is obtained by increasing ∂xb, and it is obvious from Fig. 1 that the nondimensional stratification amplitude increases with . Increasing for fixed is obtained by decreasing the tidal period T. Assuming a completely mixed water column after flooding, the decreased tidal period allows less time for destratifying the water column, such that the nondimensional stratification amplitude is also decreased. For weak mixing and strong stratification (relatively large and relatively small ), a periodic steady state is not obtained, since the stratification is increasing from tidal period to tidal period. This phenomenon is called runaway stratification and is due to the physically wrong assumption of a constant horizontal density gradient in the presence of strong stratification, see e.g., Blaise and Deleersnijder (2008). The classical threshold value for separating periodic stratification from complete stratification, [Simpson et al. (1990); see the red line in Fig. 1], is only valid for . This is equivalent to ∂xb ≈ 1.6 × 10−6 s−2 [corresponding to a salinity gradient of about 2.4 × 10−4 (g kg−1) m−1] and a depth to tidal velocity ratio of H/u ≈ 36 s or, for a water depth of H = 20 m, a tidal velocity of about umax ≈ 0.56 m s−1.
In Fig. 2, the water column dynamics for one specific situation (, , Ewt = 5.0 × 10−3 ≡ 50°N) is shown in detail. Being at the margin of permanent stratification, this situation is subject to a strong ebb–flood asymmetry, with a flood current for 0 < t̃ < 0.5 and an ebb current for 0.5 < t̃ < 1. By the end of the ebb, a stable stratification has evolved, which is eroded during flooding such that at the end of the flood phase slightly unstable stratification can be seen. This is clearly visible in the positive buoyancy production B̃ covering the entire water column. During that period, the eddy diffusivity ν̃′t has a distinct maximum. The dissipation rate ε̃ shows maxima during full flood and during full ebb, with a time lag relative to the bottom stress increasing with distance from the bed. Remarkable here is a local maximum in the dissipation rate in the upper half of the water column occurring toward the end of the flood. This is a feature that has already been observed in the field [see Rippeth et al. 2001) for results from Liverpool Bay (see also section 4) and Fisher et al. (2002) for results from the Rhine outflow region]. Both of these datasets have been successfully reproduced by numerical models [by Simpson et al. (2002) for Liverpool Bay and by Souza et al. (2008) for the Rhine outflow region]. This secondary dissipation maximum (which has not yet been discussed for situations with longitudinal density gradients as in Liverpool Bay) occurs during significantly stable stratification and negative buoyancy production, mainly due to a high longitudinal shear production, . The transverse velocity component υ̃ that is forced to a zero vertical mean has a significant positive vertical gradient, with maximum values of about ∂z̃υ̃ ≈ 1 and a mean value of about ∂z̃υ̃ ≈ 0.5 (see also Fig. 3). The latter value is close to the thermal wind solution of (19), that is, the vertical derivative of the balance between the rotation and the pressure gradient, which can be written as . This ratio has a value of for the present scenario. Although the transverse shear production, , is relatively small, it does at times (e.g., at full flood with t̃ = 0.25) indirectly contribute to the shear production, when in the upper half of the water column the transverse shear is exceeding the thermal wind balance and is thus increasing the longitudinal shear via the time derivative term in the ũ term in Eq. (19); see also the discussion in section 5.
The sensitivity of the tidal cycle of stratification, shown as the nondimensional potential energy anomaly ϕ̃, and the residual profiles of the longitudinal and transverse velocity, ∫ ũ dt̃ and ∫υ̃ dt̃, are shown in Fig. 3, for a fixed horizontal buoyancy gradient, ∂xb = 10−6 s−2, at 50°N and a number of modified horizontal Richardson numbers, . Clearly, the amplitude of ϕ̃ does directly increase with . For small , unstable stratification persists during the entire second half of the flood, with the onset of unstable stratification delayed for increasing . The switch to stable stratification for all simulations occurs around the beginning of the ebb (t̃ ≈ 0.5). For , no unstable stratification is obtained throughout the tidal cycle, although simulation results are still periodic. Residual longitudinal velocity profiles show a pronounced positive (onshore, up estuary) maximum near the bottom with values of around 0.05 for . This produces an upstream-directed residual current of about 5% of the tidal velocity amplitude. Residual transverse velocity profiles all show a strong positive shear in the upper half of the water column, which for is similar to the respective thermal wind shear (not shown). For smaller values of , the mean shear is smaller than the thermal wind shear, due to the increased role of vertical mixing.
Figures 4 –6 show the strong dependence of the tidal dynamics under the impact of horizontal density gradients on the latitude with no wind (Fig. 4), with weak westerly (onshore, up estuary) wind (Fig. 5) and with weak easterly (offshore, down estuary) wind (Fig. 6). For lower latitudes, mixing is weaker and periodic simulations are obtained only for smaller . The reason for this may be the reduced transverse shear, which is a consequence of the reduced rotation in the thermal wind balance. Under weak westerly (onshore, up estuary) winds, the destabilization of the water column for relatively high is also strong enough to overcome the stratifying effects of tidal straining. This effect is not due to extra mixing (the velocity scale uwt is kept constant, such that a weak wind results in reduced tidal velocity amplitude) but rather is due to shear induced by the surface wind stress, shearing less buoyant ocean water over more buoyant estuarine water. For weak easterly wind (offshore, down estuary; see Fig. 6), the situation is opposite: estuarine water is sheared over oceanic water due to the wind stress, such that the water column is stabilized and permanent stratification is supported. For easterly winds, there is a wide range of situations with periodic stratification without the occurrence of unstable stratification. These results are in good qualitative agreement with the observations of Scully et al. (2005).
This can be seen in more detail in Fig. 7, where the tidal evolution of the stratification and residual longitudinal velocity profiles are shown for relatively weak winds from east and west. For westerly winds, all values of the modified horizontal Richardson number up to lead to periodic stratification with the occurrence of unstable stratification, whereas for easterly winds this is obtained only for . Residual longitudinal profiles for westerly wind are S shaped, with slightly reduced (compared to no wind) landward residual flow, and residual profiles for easterly wind show significantly increased mean shear between the surface and the bottom, such that stratification increased through interaction with the horizontal density gradient. The transverse residual shear in the upper half of the water column is reduced by westerly winds and increased by easterly winds.
In Fig. 8, the effects of changed bed roughness length are studied for a situation with relatively weak westerly (onshore, up estuary) winds. Except for a very high roughness length of h̃0b = 0.1, the periodic stratification is only weakly influenced. However, for these moderate roughness lengths, the residual velocity profiles are also significantly different, indicating that the roughness length plays an important role in transporting the residual suspended matter. Also, for extremely rough beds with h̃0b = 0.1, the stratification is strongly reduced due to stronger mixing.
For relatively moderate winds (ru = 0.2), the dependence of the dynamics on the wind direction is shown in Fig. 9. This demonstrates that the wind direction has a substantial impact on the periodic stratification. For westerly and easterly winds, the effects seen in Figs. 5 and 6 for ru = 0.1 are strongly enhanced. It is interesting to note the effects for northerly and southerly winds, which are strongly reducing or increasing the transverse shear, and thus reducing or enhancing the vertical mixing, such that northerly winds strongly enhance stratification whereas southerly winds also support periodic stratification for high . Although the impacts of northerly and easterly winds are similar, just as the westerly and southerly winds are similar, the mechanisms determining how these winds change the periodic stratification are different. Longitudinal winds (from east or west) directly change the stratification since they are along the density gradient, whereas transverse winds (from north or south) change the mixing through the transverse shear. These effects, which have been discussed here for relatively moderate winds, are strongly enhanced for relatively strong winds (see Fig. 10).
4. Comparison to field data
High-resolution field data have been presented and analyzed for Liverpool Bay (53°28.4′N, 3°39.2′W) over 25 h on 5–6 July 1999, as reported by Rippeth et al. (2001) for a tidal straining situation with strong SIPS. Observed water column parameters are longitudinal and transverse velocity components, temperature, and salinity, as well as the turbulent dissipation rate. The tide was almost rectilinear and parallel to the direction of the horizontal density gradient, and the wind forcing was weak, such that these observations provide a suitable test for the present theory. For numerical simulations of these data, see Simpson et al. (2002) and Burchard (2002).
Here, the dynamics of one tidal cycle starting at 2200 UTC 5 July are reanalyzed in terms of nondimensional parameters. The tidal velocity amplitude is umax = 0.54 m s−1, and the wind forcing is weak with ru < 0.05, such that uwt ≈ umax. The mean water depth is H = 32 m with a tidal amplitude of 2.6 m. The horizontal buoyancy gradient is estimated from a simplified advection equation for the buoyancy:
where the brackets denote the vertical and temporal averaging. The resulting horizontal buoyancy gradient is ∂xb = 5.36 · 10−7 s−2, such that the relevant nondimensional parameters for this flow situation are , and .
For this situation, the idealized simulations analyzed in Fig. 4 predict a strong SIPS with an amplitude of Δϕ̃ ≈ 15 (see the two bottom panels in Fig. 4 for the latitudes ϕ = 50° and ϕ = 60°, between which the present scenario is situated).
It should be noted that the estimation of the horizontal buoyancy gradient is of limited accuracy: Rippeth et al. (2001) estimated a value of ∂xb ≈ 6.8 × 10−7 s−2 based on CTD data along a transect through Liverpool Bay; Simpson et al. (2002) corrected this value for tidal displacement to ∂xb ≈ 8.2 × 10−7 s−2, resulting in elevated values of and , respectively. Following Fig. 4, the latter value would be near the margin of permanent stratification. There is also some ambiguity in estimating the tidal velocity amplitude, since the flood vertical mean of the longitudinal velocity is 0.59 m s−1 and the ebb vertical mean is 0.49 m s−1, as well as in the water depth, varying between 29.4 and 34.6 m, with both adding uncertainty to the estimates of , , and Ewt.
It is instructive to study the nondimensionalized observations of ũ, υ̃, b̃, and ε̃ shown in Fig. 11, which are directly comparable to the idealized model simulations (top four panels) in Fig. 2. These results have been calculated relative to the actual water depth and not to the mean water depth in order to better account for the dynamical effects of the water depth. When compared to the idealized solution, the velocity observations show longer, but less intense, ebb current. The complex transverse flow pattern is similar to the idealized model simulations, with a slightly more pronounced near-bed southward component. The buoyancy observations show an extended period of low stratification between full flood and full ebb, and a weaker maximum stratification after ebb. The dissipation rate observations show qualitative and quantitative results that are strikingly similarity to the idealized model simulations. Near the bed, dissipation rates are comparable during flood and during ebb, but in the upper half of the water column, dissipation rates are much stronger during flood than during ebb. During flood, a pronounced local dissipation rate maximum is visible in the observations, a pattern that is clearly present in the model simulations (Fig. 2), which has already been discussed in section 3. The observed tidal asymmetry is stronger in the observations, probably due to the stronger flood current. In addition, the observed phase lag of the dissipation rate with increasing distance from the bed is shorter than in the idealized model simulations, probably due to the model assumption of a constant buoyancy gradient.
The differences in stratification are also evident from the time series of nondimensional stratification, Δϕ̃ (see Fig. 12). The observed stratification is much stronger than in comparable idealized simulations (see Fig. 3), if based on the horizontal buoyancy gradient of ∂xb = 5.36 × 10−7 s−2. It would look more reasonable if the nondimensional stratification were based on the higher estimates discussed above. One reason for the discrepancy may be the increased stratification resulting from the net heat flux into the water column, which peaked at 800 W m−2 during the noon before the start of the analyzed tidal period. The observed longitudinal tidal mean velocity profile shows the expected onshore flow near the bed, and an offshore tendency increasing with distance from the bed, but in the upper half of the water column, it is directed onshore again. This effect could be due to the weak westerly wind component, the tidal asymmetry, or a number of other processes not reproduced in the idealized model simulations. Some similarities in this observed longitudinal residual velocity profile to the idealized simulations with weak westerly wind (see Fig. 7) can be seen. The residual transverse velocity profile is similar to the idealized simulations with and without weak westerly wind, showing a near-bed southward velocity and a near-surface northward velocity of up to 0.15 m s−1.
The potential energy anomaly equation (28) does not include the direct effects of the earth’s rotation and surface stresses. Instead, the straining term includes the velocity profile as a diagnostic quantity, with onshore shear decreasing and offshore shear increasing stratification. Processes changing the shape of the velocity profiles and therefore changing the rate of growth or destruction of the stratification can be quantified by taking the time derivative of the potential energy anomaly equation, which after neglecting the mixing term and inserting the momentum Eq. (1) is of the following form:
with the vertical mean Reynolds stress
and the vertically averaged northward velocity component, 〈υ〉.
Processes changing the shape of the velocity profiles do not directly change the stability of the water column, but change the trend in stratification or destratification. This is expressed by ∂ttϕstrain. A positive value of ∂ttϕstrain implies that an increase of ϕ is accelerated by one of the three forcing terms or a decrease is slowed down. A negative value of ∂ttϕstrain implies that an increase of ϕ is slowed down or a decrease is enhanced. In this sense, the three terms in (32) have a clear meaning, which is discussed here for ∂xb > 0 (estuary to the right) and f > 0 (Northern Hemisphere):
Nonlinear stress forcing (term A)—Under stably stratified conditions (ϕ > 0), an onshore wind (τsx < 0) causes a decrease of the water column stress τx in the near-surface region, which is larger in the upper half than in the lower half of the water column, such that term A becomes negative and the water column stability tendency is decreased. For offshore winds (τsx > 0), the opposite happens such that the water column stability tendency is increased. This is in accordance with the fact that offshore winds stabilize and onshore winds destabilize the water column [see Figs. 5 –10 and Scully et al. (2005)]. With the same arguments, decreases in onshore flow and increases in offshore flow have the tendency to increase the rate of water column stabilization, and vice versa. All these processes are strongly enhanced for stably stratified water columns, since the transfer of bottom or surface stress changes through the water column is strongly reduced, such that the stress profile is of nonlinear shape, resulting in a nonzero term A.
Thermal wind forcing (term B)—As discussed in section 3, the thermal wind balance forces a tidal mean positive vertical gradient of the transverse velocity, such that term B is generally negative and the stability tendency of the water column is decreased. This is in accordance with the fact that for fixed external flow parameters, the water column stability is higher for lower latitudes; that is, the threshold horizontal Richardson number for permanent stability is decreased (see Fig. 4).
Indirect density gradient forcing (term C)—This term is always positive and thus increases the tendency for water column stability. This can be easily understood for a situation with irrotational flow at rest and no wind forcing: in this situation, a horizontal density gradient will increase the water column stability indirectly for all density gradient directions and could also be labeled as the estuarine gravitational circulation term. For irrotational ( f = 0) stationary conditions this term C is exactly in balance with term A.
The major result of the present study is that tidal straining dynamics do not only strongly depend on the ratio of tidal straining to mixing (as expressed here by the modified horizontal Richardson number), but also on the ratio of the vertical mixing time scale to the tidal period (proportional to the modified inverse Strouhal number) and the ratio of the vertical mixing time scale to the inertial period (proportional to the modified inverse Ekman number). Furthermore, the induced mean shear by the surface wind stress seems to play a far more important role than the additional mixing provided by the wind stress.
It has been shown that the classical threshold value for identifying the shift from periodic to permanent stratification expressed as a critical horizontal Richardson number only holds for a small area in the parameter space. The role of the latitude for the tidal dynamics in coastal and estuarine regions has been investigated in many previous studies (e.g., Prandle 1982; Sharples and Simpson 1995; Lerczak and Geyer 2004; Valle-Levinson 2008). However, the present study shows for the first time that the Coriolis acceleration tends to decrease the water column stability in tidal flow under the effect of horizontal density gradients. It could be shown by analyzing the second time derivative of the potential energy anomaly, ϕ [see Eq. (32)], that the transverse velocity shear generated by the thermal wind balance decelerates the tendency for increasing water column stability. From this equation, it is also clear that the straining term includes three basic processes: nonlinear straining (including effects of surface stress and tidal mixing asymmetry), Coriolis rotation, and gravitational circulation. The relative contribution of these processes can however not be analyzed with the present theory.
The dominance of the wind stress over wind mixing already suggested by the field observations carried out by Scully et al. (2005) is an important outcome of the present study. Since significant wind straining effects are already present at relatively low wind speeds (cf. Figs. 4 –6), how an estuary is oriented toward the prevailing wind direction clearly matters for the stratification and residual flows. Estuaries and coastal seas with prevailing onshore (up estuary) wind should have a stronger tendency toward decreased estuarine circulation and subsequently reduced stratification and increased mixing, while estuaries with prevailing offshore (down estuary) winds should have a tendency toward increased estuarine circulation and reduced mixing.
The near-surface flood maximum of the dissipation rate, which has been observed in two field studies so far (Rippeth et al. 2001; Fisher et al. 2002), has also been revealed here as a characteristic feature for tidal flow interacting with longitudinal density gradients, and not only for flow with transverse density gradients, as discussed by Souza et al. (2008). Therefore, a generalization of the present study to flows where the density gradient is not aligned with the tidal flow direction and also to tidal flows with nonrectilinear tides would be necessary.
Transient regimes switching from spring tide to neap tide or vice versa could not be included in this study since it would have extended the parameter space by at least one or two dimensions. It has been recognized in the literature that such transient regimes play an important role for mixing and stratification in coastal seas (Sharples and Simpson 1995). The manner in which the present results can be interpreted for such transients is that two periodic situations with different values of umax can be interpreted as spring and neap tide, respectively.
The effects of surface waves on the tidal flow due to wave–current interaction in the bottom boundary layer have not been included in this study. The sensitivity of tidal residual flow profiles on the relative bed roughness gives an idea about how strong the dynamical effects of surface waves could be. Since in coastal areas the wave height and length strongly depend on the wind direction, the significant differences in residual flows and stratification due to mean shear generated by wind from different directions could be further enhanced.
The present study only considers estuaries and coastal sea regions where the density gradient is generally aligned with the direction of the tide. It should however be noted that there are coastal regions where this is typically not the case, such as in the Rhine region of freshwater influence [see, e.g., Simpson and Souza (1995)], where the dynamics are substantially different than in the cases considered here.
It has often been stated that field observations are challenging numerical modeling activities (Rippeth et al. 2001; Stacey et al. 2001), but the results of the present study should motivate a number of field studies in tidally dominated coastal seas under nonstandard conditions, such as low or high latitudes, dominant diurnal tides, or winds from various directions.
This work has been supported by a Kirby Laing Fellowship at the School of Ocean Sciences of the University of Wales in Bangor in 2008 and the EU-funded European Coastal-Shelf Sea Operational Monitoring and Forecasting System project (ECOOP; Contract 36355). The author is grateful to Florence Verspecht (Menai Bridge, Wales), and Tom Rippeth (Northop, Wales), Alejandro Souza (Llanfairpwllgwyngyllgogerychwyrndrobwllllantysiliogogogoch, Wales), Lars Umlauf (Rostock, Germany) and Robert Hetland (College Station, Texas), as well as to two anonymous referees for critical comments on the manuscript and to Karsten Bolding (Asperup, Denmark) for providing scripts for the execution of multiple simulations with GOTM (information online at www.gotm.net). A first version of this theory was presented in April 2006 at the Coastal and Shelf Seas—Present Understanding and Future Challenges conference in Bangor (Wales) marking the retirement of Prof. John Simpson.
Corresponding author address: Hans Burchard, Leibniz Institute for Baltic Sea Research, Warnemünde, Seestraße 15, D-18119 Rostock, Germany. Email: email@example.com