The effect of density gradient on boundary flow

Estuarine, Coastal and Shelf Science 183 (2016) 163e178

Contents lists available at ScienceDirect

Estuarine, Coastal and Shelf Science journal homepage: www.elsevier.com/locate/ecss

The effect of density gradient on boundary flow Zhuo Zhang a, b, *, Zhiyao Song a, b, c, Cheng Chen d, Fei Guo a, b, Dong Zhang a, d, Di Hu a, b a

Key Laboratory of Virtual Geographic Environment (Ministry of Education), Nanjing Normal University, Nanjing, 210023, PR China Jiangsu Center for Collaborative Innovation in Geographical Information Resource Development and Application, Nanjing, 210023, PR China c Jiangsu Provincial Key Laboratory for Numerical Simulation of Large Scale Complex Systems, School of Mathematical Sciences, Nanjing Normal University, Nanjing, 210023, PR China d State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering, Nanjing Hydraulic Research Institute, Nanjing, 210029, PR China b

a r t i c l e i n f o

a b s t r a c t

Article history: Received 29 February 2016 Received in revised form 12 October 2016 Accepted 18 October 2016 Available online 19 October 2016

Based on Prandtl's mixing-length theory, this study proposes new theoretical formulae of velocity profiles resulting from the forcing by a constant horizontal buoyancy gradient in unstratified and stably stratified flows. Based on the one-dimensional water column momentum equation, the vertical turbulent shearing stress profile is found to deviate from a linear distribution and follow a parabolic distribution, differing from that in neutral flow. The shearing stress curves upward with the current following the density gradient, and curves downward with the currents opposite to the density gradient. For a constant eddy viscosity, the well-known estuarine circulation is obtained through the parabolic shearing distribution. For a vertically parabolic eddy viscosity, the new-proposed velocity profile by Burchard & Hetland (2010) is obtained. In this paper, we estimate the viscosity profile based on Prandtl's mixing-length theory and then derive the new formulae of the velocity profiles. Through comparison with the numerical turbulence model, the velocity profiles presented in the paper agree well with the numerical results. In addition, the new velocity profiles can be applied to determine the valid range and evaluate errors of the log-fit in baroclinic boundary flows. © 2016 Elsevier Ltd. All rights reserved.

Keywords: Velocity profiles Turbulent shearing stress Eddy viscosity Density gradient Boundary layer

1. Introduction In estuarine and coastal areas, numerous velocity profiles through the bottom boundary layer (BBL) indicate a deviation from the classic logarithmic depth dependence. Some investigators attributed the deviation to several factors, such as acceleration and deceleration of tidal currents, density stratification, the transport of material as bed-load, etc. (Collins et al., 1998). Among these factors, the effects of density gradient, including longitudinal density gradients and vertical stratification on the velocity profiles, have been recognized by many investigators as a significant factor in the action of current circulation and salt water intrusion in the field of estuarine and coastal physical oceanography (Soulsby and Wainwright, 1987; Hetland and Geyer, 2004). Early studies of the well-mixed estuaries assumed that the flow was influenced persistently by the horizontal density gradient, namely the baroclinic effect, which brought the near-bed seawater

* Corresponding author. Key Laboratory of Virtual Geographic Environment (Ministry of Education), Nanjing Normal University, Nanjing, 210023, PR China. E-mail address: [email protected] (Z. Zhang). http://dx.doi.org/10.1016/j.ecss.2016.10.025 0272-7714/© 2016 Elsevier Ltd. All rights reserved.

landward and counteracted some near-bed runoff seaward. The well-known estuarine circulation theory was proposed and analyzed by Pritchard (1952, 1954), Hansen and Rattray (1965), and Chatwin (1976), with the assumption that the geometry was idealized and the mixing was constant throughout the estuary. The formula of estuarine circulation was later used widely to interpret observations in the well mixed and partially mixed estuaries by many scholars (Geyer, 1997; Geyer et al., 2000; Vinita et al., 2015). Their theory assumed that, if the section averaged salinity increased linearly along the estuary, then both the gravitational and the stratification circulation were constant along the estuary. Prandle (2004) established theoretical solutions for the vertical profiles of tidal currents and residual currents under a linear approximation to the familiar quadratic bed friction law. His results indicated that neglecting the longitudinal density gradient was entirely valid for vertical structure of tidal currents but incorrect for that of the residual current because the density gradient introduced small, but significant, tidally averaged residual currents. Different from Hansen and Rattray (1965) and Prandle (2004), who assumed a constant eddy viscosity, Burchard and Hetland, (2010) presented a log-linear expression for tidal-mean currents, assuming a more realistic parabolic eddy viscosity, showing that the intensity of the

164

Z. Zhang et al. / Estuarine, Coastal and Shelf Science 183 (2016) 163e178

gravitational circulation scaled with the horizontal Richardson number. MacCready (2004) concluded that all these theories relied on knowledge of the vertical turbulent momentum flux and the Hansen and Rattray (1965) solution always became diffusiondominated near the mouth. Recent studies concentrated on internal tidal asymmetry (called tidal straining), which was thought to be a major mechanism to form estuarine circulation (Simpson et al., 2002, 2005). On the other hand, the investigation into the effects of density stratification on the velocity profiles was initiated by scholars and engineers outside the field of estuarine and coastal physical oceanography. Numerous experimental results from atmospheric fluid dynamics introduced the Monin-Obukov length to account for the bottom boundary layer under the impact of stable stratification, such as a boundary layer cooled intensely from below (Monin and Yaglom, 1971). These ideas have been extended to density stratification brought by salt water intrusion, thermohaline stratification, or suspended sediment in marine bottom boundary layers through some revisions. For example, Anwar (1983) and Sanford et al. (1991) applied log-linear profiles to explain the stratification in the Great Ouse and Chesapeake Bay estuaries. Adams and Weatherly (1981) examined the stably stratified oceanic bottom boundary layer using the Mellor-Yamada turbulence closure (MY) model (Mellor and Yamada, 1982) and suggested an expression for the near-bottom velocity profile, which linked the reduction of the Karman coefficient to the constant Richardson's number throughout depth. Friedrichs and Wright (1997) applied the improved log-quadratic velocity profile to obtain realistic estimates of the bottom stresses, which was usually overestimated by the classical log profiles. Herrmann and Madsen (2007) examined the sedimentinduced stratification effects on the velocity profiles and sediment concentration distribution in a steady and uniform turbulent flow. However, the above studies neglected the horizontal density gradient effect on velocity profiles. In numerous estuaries, the salinity distribution changes alternatively between being well-mixed and being partially-stratified within tide cycles. In addition, the regime of the estuarine mixing may vary with the different locations of the estuary according to the distance from the river mouths. As a result, the influence of horizontal density gradient and stratification (by vertical density gradient) may coexist and interact in the hydrodynamic mechanics of the estuarine circulation. Unfortunately, most of the former theoretical formulae concentrated on either the horizontal density gradient (e.g., Hansen and Rattray (1965); Burchard and Hetland, (2010)) or stratification (e.g., Friedrichs and Wright (1997); Herrmann and Madsen (2007)). The theoretical expression of the velocity profile explicitly including both factors is rarely presented. Although an increasing number of numerical simulations can take multiple factors into account and obtain agreeable results with the measurements (Mazumder and Ghoshal, 2006), they are incapable of differentiating the mechanism of the horizontal density gradient and the stratification influence as explicitly as the theoretical formula do. Moreover, some boundary layer parameters such as friction velocity, roughness length and drag coefficient in estuarine environments, estimated based on theoretical velocity profiles, are required for: (i) defining flow conditions at the sediment-water interface; and (ii) predicting sediment and contaminant transport rates (Collions et al., 1998). In this paper, we investigate more closely into the impact of the horizontal density gradient on unstratified and stably stratified boundary flow through new analytic solutions of shearing stress distribution, viscosity and velocity profiles. Thus, we have restricted our attention to well-mixed and periodically stratified estuaries, neglecting the earth's rotation, wind stress and advection. For the flow at peak flooding or ebbing time, the tidal

acceleration is weak and can be ignored. The tidal straining induced by tidal asymmetry is also excluded for the instantaneous flow. The paper is organized as follows: First, the shearing stress profile is derived, including the interaction between the horizontal density gradient and the bottom friction. Second, the analytic velocity profiles considering the horizontal density gradient and stratification are proposed according to the Prandtl's mixing length and Monin-Obukhow's stratification theory. Third, a series of ideal numerical simulations using vertical parameterization through Mellor-Yamada turbulence closure (MY) model (Mellor and Yamada, 1982) is performed to validate the theoretical results. The fitting results indicate that, although some differences appear between the MY model and the Prandtl's mixing model, the analytic velocity profile fits the numerical results satisfactorily. Moreover, the new-proposed velocity profiles can be applied to determine the valid range of the log-fit with enough accuracy for estimating the friction velocity in baroclinic flows. While the log-fit is not accurate enough, the estimating error can be given.

2. Theory 2.1. The vertical shearing stress The momentum equation regardless of earth's rotation along the x-coordinate is

vu vuw 1 vP 1 vt þ þ  ¼0 vt vx r0 vx r0 vz

(1)

with the velocity u and w along the x-coordinate and z-coordinate, the pressure P, the shear stress t along the x-coordinate and the averaged water density r0. According to the static pressure assumption, the pressure is given by

1 vP v2 g ¼g þ r0 vx vx r0

Z2 z

vr dz vx

(2)

where 2 denotes the water elevation,g denotes the gravity acceleration and r is the water density at a certain depth. Substituting (2) into (1) and neglecting the non-linear advection term and local acceleration term, one can obtain

1 vt v2 g ¼g þ r0 vz vx r0

Z2 z

vr dz vx

(3)

By integrating Eq. (3) in conjunction with the boundary condition and eliminating the barotropic pressure by Eq. (1), the solution of the shear stress is obtained as

t ¼ tb ð1  sÞ  0:5gD2

 vr  2 s s vx

(4a)

or

  t ¼ ð1  sÞ  Rix s2  s tb

(4b)

with the bottom shear stress tb, the dimensionless vertical height s ¼ hþz D , the static water depth h and the instantaneous water depth D. The modified horizontal Richardson number is defined as

Z. Zhang et al. / Estuarine, Coastal and Shelf Science 183 (2016) 163e178

Rix ¼

gD2 vr 2r0 u2* vx

(5)

which physically means the relative contribution of the horizontal density gradient force to the bottom friction force. The larger is the value of Rix , the more important role is taken by the horizontal buoyancy gradient over the bottom friction on boundary flow. Burchard and Hetland, (2010) proposed the similar form of Riwt x , in which he replaced the friction velocity with the current velocity to allow easier comparison of the theory with the field observations. Fig. 1 indicates that the shearing stress profiles deviate from the linear distribution under the effect of the horizontal buoyancy gradient. When the flow follows the density gradient (Rix>0), the shearing stress curves upward. When the flow is against the density gradient (Rix<0), the shearing stress curves downward. In the latter case the shearing stress profiles may have a zero point at a certain depth, corresponding to the extreme point in the velocity profiles. Assuming the constant viscosity Nz, we can derive the velocity profile from Eq. (4) as follows:

u¼

  tb D  gD3 vr  2 3s  2s3 2s  s2 þ 2r0 Nz 12r0 Nz vx

Given the runoff flux ur ¼

u¼

gD3

vr 

48r0 Nz vx

R1 0

(6)

uds, Eq. (6) can be written as



 3   8s3 þ 15s2  6s þ ur 2s  s2 2

(7)

which is exactly the classical estuarine circulation presented by Hansen and Rattray (1965). Burchard and Hetland, (2010) recently proposed a new velocity profile based on a more realistic parabolic and constant-in-time eddy viscosity profile given by

      gD2 vr 1 z þ z0 z ur z ln þ ln  u¼ H z0 2kru* vx 2N z0 N in which N ¼ D1

(8)

  RD zþz0 ln dz z0 z0

2.2. The eddy viscosity based on Prandtl's mixing-length theory The assumption of a constant viscosity is an acceptable approximation at the high or low water slack. Therefore Eq. (7) was often verified by the data in estuaries with negligible residual runoff (urz0). However, it is a rough approximation for the peak

165

flooding or ebbing flow in estuaries. Though the parabolic eddy viscosity is more realistic than the constant one, it is experiential and short of theory foundation. Thus, the Prandtl's mixing length model is applied to obtain the valid eddy viscosity profiles. The turbulent shear stress is given by





vu vu t ¼ rl2   vz

 (9)

vz

where l is the mixing length of the eddy, representing the transverse distance of fluid elements exchanging momentum. In the unstratified boundary flow, l can be described by

l ¼ kzð1  sÞ0:5

(10)

with the van Karman number k¼0.4. Substitute Eqs. (9) and (10) into Eq. (4) and simplify the equation, then we get

vu u* pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ 1 þ Rix s vz kz

ðRix  1Þ

8 > u* pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > > 1 þ Rix s < vu kz ¼ > u* pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi vz > > : 1  Rix s kz

s < sR

(11a)

ðRix <  1Þ

(11b)

s  sR

where sR ¼ Ri1x is the height where the shearing stress is zero when Rix is less than 1.    The eddy viscosity y ¼ l2 vu vz  can be expressed by

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

y ¼ ku* Dsð1  sÞ 1 þ Rix s ðRix  1Þ

y¼

(12a)

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

ku* Dsð1  sÞp1 þ Rix s ffi s < sR ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ku* Dsð1  sÞ 1  Rix s s  sR

ðRix <  1Þ

(12b)

It shows that when Rix<1, the viscosity profiles estimated by the mixing-length theory are characterized by the B-shape distribution: Besides at the surface and the bottom, an additional zero point appears at a certain height. The maximum viscosity first reduces with the decrease of Rix until Rix¼2, then increases rapidly when Rix<2. For deriving simple velocity profiles, (Eq. (12b)) can be fitted by the parabolic function as

y ¼ ku* Dsð1  sÞFðRix Þ ðRix <  1Þ

FðRix Þ ¼

(12c)

0:212Rix þ 0:956 ð  2 < Rix <  1Þ 0:0152Rix 2  0:3808Rix  0:1688 ðRix <  2Þ

A comparison of Eqs. (12b) and (12c) is shown in Fig. 2. The approximate expression estimates the same maximum and similar distribution of viscosity with that directly deduced by the mixinglength theory. The validity of the approximation in the derived velocity profile will be tested in Section 4. Eq. (12c) can be regarded as a modification in baroclinic flow to the empirical parabolic viscosity in neutral flow. In stratified boundary flow, according to the theoretical and observational investigations by Monin and Yaglom (1971), and Turner (1973), l should be replaced by l' in the form of

l0 ¼ l Fig. 1. The shearing stresses profiles for different values of Rix.

. z 1þa L

and viscosity is modified to

(13)

166

Z. Zhang et al. / Estuarine, Coastal and Shelf Science 183 (2016) 163e178

Fig. 2. The dimensionless viscosity profiles for different values of Rix. The real line is for Eq. (12) (b) and the dash line is for Eq. (12) (c).

y0 ¼

y 1 þ a sLD

(14)

where L is the Monin-Obukov length which varies on the gradient g vr

Richardson number Riz ¼   vz2 and az5 (Webb, 1970; Kundu, r

vu vz

2.3. The velocity profile

k

! pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ Rix s  1 þ 2 1 þ Rix s þ C ln pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ Rix s þ 1

(15a)

For stratified flow:

u¼

u*

k

u*

kFðRix Þ

½ln s þ Rix s þ C

(16a)

u¼

   u* aD Ri ln s þ Rix s þ s þ x s2 þ C L kFðRix Þ 2

(16b)

with the integration constant C decided by the fit data. With the definition of relative roughness height s0, all above velocity profiles, taking Eq. (15b) for example, can be written in another form of

If Rix1, direct integration of (11a) yields. For unstratified flow:

u*

u¼

For stratified flow:

1990) is a parameter estimated by observations from stably stratified boundary layers. Similar form with (14) was applied by Styles and Glenn (2000) in stratified flow.

u¼

shape viscosity Eq. (12b) results in the following velocity profiles: For unstratified flow:

! pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 1 þ Rix s  1 bu* D þ 2 1 þ Rix s þ ð1 þ Rix sÞ2 þ C ln pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Rix L 1 þ Rix s þ 1 (15b)

where b ¼ 23ak ¼ 8:33 if Rix<1, substituting the parabolic viscosity Eq. (12c) for the B-

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ Rix s  1 1 þ Rix s0 þ 1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u* 4 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  þ 2 1 þ Rix s ln pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u¼ k 1 þ Rix s þ 1 1 þ Rix s0  1 3 i pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 3 bu* D h ð1 þ Rix sÞ2  ð1 þ Rix s0 Þ2  2 1 þ Rix s0 5 þ Rix L (17) Some special cases for (17) are listed as follows:

Z. Zhang et al. / Estuarine, Coastal and Shelf Science 183 (2016) 163e178

Rix /0

u¼

u*

k



 ln

s aD þ s0 L

(18)

which was proposed and applied by Webb (1970) and Kundu (1990) in describing the stratified flow.

Rix /0 and L/ þ ∞

u¼

u*

k

ln

s s0

where l is the turbulence scale, and q is a scaling velocity for the large energy-containing eddies (q2/2 is the turbulent kinetic en~ ¼ 1 þ E ðl=kLÞ is the proximate wall function and ergy). W 2 L1¼(2z)1þ(Dz)1. E1¼1.8 E2¼1.33 are the constants. After derivation of q2 and l, the eddy viscosity is written as

yt ¼ qlSm (19)

which is the classical log velocity profile for the neutral flow.

167

(25)

0:42753:354Gh Where Sm ¼ ð134:676G is the stability h Þð11:672Gh Þ 2 Gh ¼ ql2 rg r 0 The vertical eddy diffusion coefficient is written as

function,

2.4. Scaling analyses

Kq ¼ 0:2ql

Some assumptions are necessary for simplifying the hydrodynamic equations so that an analytic solution can be derived. Nevertheless, the applicable extent is restricted upon these assumptions. Here we focus on the estuarine environment and the followings listed are the scales typical of an estuary:

For solving the hydrodynamic and turbulence closure equations numerically, the following boundary condition should be supplemented as follows

y

5

Estuarine length: l ~10 m Width: d ~103 m Depth: h ~102 m Flow longitudinal speed: V ~1 m s1 Wind speed: W ~1 m s1 Coriolis parameter: F ~104 s1



From the above scales, the aspect ratio (d¼l/d¼100) is large, so the earth rotation effect and lateral movement can be neglected and the flow can be consider one-dimensional. The Rossby number (R0¼V/ Fl¼0.1) is small so that the advection momentum term can be neglected as long as the mild-slope topography dominates in major estuarine areas. Moreover, the statistic data show that the windy days (with the wind power being stronger than Level 6) cover less than 50 days per year. Thus, the wind shear stress can be neglected for most of the calm days. Anyway, the assumptions are reasonable for usage of the new velocity profiles in estuarine environments. 3. Numerical model The vertical structure of the flow can be obtained using the vertical two-dimensional Navier Stokes equations and the continuity equation:

vDu vu v2 þ þ ¼0 vx vs vt vu vu vu v2 1 þu þw ¼ g þ vt vx vz vx r0

(20) Zz 0

vr v2 u dz þ y 2 vx vz

v2 vUD ¼0 þ vx vt

 2  vu k ¼ u2b ¼ u2* vz z¼0 lnðzb =z0 Þ   q2 ; q2 l 

z¼0

  2=3 ¼ B1 u2* ; 0

(27)

(28)

where B1¼16.6 is one of the turbulence closure constant. A finite difference approximation is used to discretize the governing equations and the boundary conditions. A staggered grid mesh system with a set of uniform rectangular M  N cells in the xand s-directions is employed. The velocity is defined at the side of the cell, and the scalar data are defined at the center of the cell. Using the semi-implicit scheme, one time step can be divided into four substeps. First, through the discretized forms of (21) and (27), the horizontal velocity is expressed as a function of water elevation; second, substituting the function into the discretized form of (22), a positive definite matrix system with unknown water elevations is obtained. The unique solution of the matrix system can be efficiently determined by preconditioned conjugate gradient iterations. Third, the horizontal velocity can be obtained by the solved elevation field, and the vertical velocity can be solved through the discretized form of (20). Four, the MY turbulence closure equations (23) and (24) are solved for the scaling velocity and the turbulence scale in conjunction with Eqs. (26) and (28), and then the eddy viscosity is obtained through (25). For details of the semi-implicit scheme, interesting readers can refer to the works by Casulli, 1990; Casulli and Cattani, 1993.

4. Model experiments and comparison

(21) 4.1. Model setup

(22)

RD where U ¼ D1 0 udz denotes the vertically averaged velocity. The MY turbulence closure model is

    vq2 D v Kq vq2 2yt vu 2 2g vr 2Dq3  ¼0  þ yt þ vt vs D vs D vz r0 vs B1 l

(26)

(23)

"   #   vq2 lD v Kq vq2 l yt vu 2 g vr Dq3 ~   E1 l þ yt  W¼0 vt vs D vs D vz r 0 vs B1 l (24)

Numerical experiments were performed using the model described in Section 3, and the results were compared with the theoretical results derived in Section 2. To investigate the distinct influences of the horizontal and vertical density gradients on the velocity profiles, various environmental density field scenarios (including four horizontal density gradients combined with two stratifications) are considered. Among these scenarios, No. 1 and No. 2 are the cases in neutral current to validate the numerical results. Nos. 3 to 6 are the cases in baroclinic currents under the downstream horizontal density gradient (following the currents). Nos. 7 to 10 are similar, but with the upstream density gradient (opposing the currents). Nos. 11 to 14 are the cases in the stratified flow under the impact of the horizontal density gradient. The computational domain is set to be a long channel with dimensions of 5000 m in length and 10 m in depth, with 10 layers in

168

Z. Zhang et al. / Estuarine, Coastal and Shelf Science 183 (2016) 163e178

the vertical direction. The numerical model starts from zero velocity and ramps up to the full values after the initial steps. The flux inlet ur is tuned to gradually approximate the designated values with the time-marching steps after the oscillation waves propagate throughout the computational domain. 4.2. Comparison of the shearing stress Eq. (4) demonstrates that the distribution of the shearing stress is determined by the bottom stress tb and Rix, in which the latter parameter causes deviation from the linear distribution. The

curvature of the shearing stress also becomes higher with the increase of Rix. To validate Eq. (4), a comparison between Eq. (4) and the numerical results is performed. The shearing stress under the impact of the downstream horizontal density gradient curves upward (shown in Fig. 3). Compared with Nos. 3 and 4, the shearing stresses of Nos. 5 and 6 show higher curvature as a result of the much larger Rix caused by the smaller flux and bottom friction stress. The shearing stress under the impact of the upstream density gradient (shown in Fig. 4) curves downward, opposite to that with the downstream density gradient. The shearing stresses of Nos. 9 and 10 have higher curvature than

Fig. 3. The shearing stresses profiles given different the downstream density gradients. Nos.1 and 2 are linear distributions in neutral currents for reference. The theoretical results are depicted in real lines and the numerical results are in scatters for comparison. The pattern of lines is followed by the following figures.

Z. Zhang et al. / Estuarine, Coastal and Shelf Science 183 (2016) 163e178

169

Fig. 4. The shearing stresses profiles given different the upstream density gradients. The parameters are listed in Table 1 for Nos. 7 to 9.

Table 1 The experiment parameters (runoff, horizontal density gradient, vertical density gradient).

Table 2 The fitting parameters (friction stress, horizontal Richardson number and the Monin-Obukov length) for the eddy viscosity profiles.

No.

ur (m s1)

vr/vx(kg m4)

vr/vz(kg m4)

No.

tb (N m2)

Rix

L (m)

1 2 3 4 5 6 7 8 9 10 11 12 13 14

0.45 0.19 0.45 0.45 0.19 0.19 0.45 0.45 0.19 0.19 0.45 0.45 0.45 0.45

0 0 1.54e-3 4.62e-3 1.54e-3 4.62e-3 1.54e-3 4.62e-3 1.54e-3 4.62e-3 4.62e-3 4.62e-3 4.62e-3 4.62e-3

0 0 0 0 0 0 0 0 0 0 3.85e-1 7.70e-1 3.85e-1 7.70e-1

1 2 3 4 5 6 7 8 9 10 11 12 13 14

0.55 0.11 0.51 0.48 0.08 0.05 0.63 0.77 0.14 0.22 0.36 0.17 0.73 0.67

0 0 1.5 4.8 9.6 46.2 1.2 3.0 5.5 10.5 6.4 13.6 3.1 3.3

e e e e e e e e e e 21 7 9 5

170

Z. Zhang et al. / Estuarine, Coastal and Shelf Science 183 (2016) 163e178

Fig. 5. The vertical distribution of viscosity given different downstream density gradients. The parameters are listed in Tables 1 and 2 for Nos. 1 to 6.

those of Nos. 7 and 8 due to the smaller flux and larger Rix. Note that the shearing stress under the upstream density gradient may intersect with the s-axis at a certain depth between 0 and 1. According to the relationship between the shearing stress and velocity, this point corresponds to the maximum extreme point in the velocity profile at the same depth. Another concern is the relationship between the bottom stress tb and Rix listed in Table 2. The relationship shows that the bottom stress decreases with the growth of the downstream horizontal density gradient at the same flux (Nos. 1, 3, and 4, and Nos. 2, 5, and 6). However, it increases with the growth of the upstream horizontal density gradient with the same flux (Nos. 7 & 8 and Nos. 9 & 10).

4.3. Comparison of the eddy viscosity profiles Conventionally, the eddy viscosity is considered to be a constant for acquiring the velocity profiles in estuarine circulation. A constant eddy viscosity is an acceptable approximation at the high or low water slack but not in peak tidal flow. Therefore, Eq. (12) is introduced to describe the eddy viscosity profile under the impact of a constant horizontal density gradient. Eq. (12) indicates that besides u*, Rix is another vital parameter to determine the distribution of the eddy viscosity, which makes the biggest difference to the empirical parabolic viscosity. Fig. 5 shows that the viscosity profiles are in accordance with the parabolic distribution when no horizontal density gradient exists (see Nos. 1 and 2). However, for the flows under the impact of

Z. Zhang et al. / Estuarine, Coastal and Shelf Science 183 (2016) 163e178

171

Fig. 6. The vertical distribution of the viscosity given different the downstream density gradients. The parameters are listed in Tables 1 and 2 for Nos. 7 to 10.

the downstream density gradient (Nos. 3, 4, 5, and 6), the distribution of the eddy viscosity is contradictory against the parabolic distribution because it increases with the growth of Rix and the reduction of tb. If the eddy viscosity followed a parabolic eddy viscosity, as Winterwerp et al., 2006 applied, the magnitude of the viscosity should be only proportional to the friction velocity u* and hence should decrease with the growth of Rix, resulting from the reduction of bottom stress tb listed in Table 2. Therefore, it is not correct to apply the same parabolic distribution in baraclinic flow as that in barotropic flow. The viscosity under the impact of the upstream density gradient is shown in Fig. 6 (Nos. 7, 8, 9, and 10), in which the B-shape

viscosity distributions are verified by the numerical results. Through the vertical transport of turbulence kinetic energy and scale, the MY turbulence obtains smooth results without a zero point appearing in the Prandtl's mixing length model. Overall, the viscosity profiles obtained using the Prandtl's mixing length model are consistent with those using the MY turbulence model, with a slight overestimation or underestimation at the extreme point. For stratified flows (Nos. 11, 12, 13, and 14), the Monin-Obukov length L must be introduced through Eqs. (13) and (14) to account for the turbulence-restraining effect by vertical buoyancy. Using a reasonable L (shown in Table 2), the theoretical viscosity is consistent with the numerical results shown in Fig. 7.

172

Z. Zhang et al. / Estuarine, Coastal and Shelf Science 183 (2016) 163e178

Fig. 7. Comparison of the vertical distribution of the viscosity in stratified (Nos. 11 to 14) and unstratified flow (Nos. 4 and 8) with ur¼0.45 m s1.

4.4. Comparison of the velocity profiles Eqs. (15) and (16) demonstrate that the velocity profiles in baroclinic flows under the impact of the density gradient deviates from the classical log profile in neutral flows. For the unstratified currents, Rix is the only factor to cause deviation, and for the stratified currents, Rix and L both cause deviation. Fig. 8 shows the velocity profiles in unstratified flows under the impact of the downstream density gradient by the numerical model (shown in scatters) and the fits using Eq. (15a) (shown in lines). The fitting parameters are listed in Table 3. The fitting profiles agree

well with the numerical results, and the estimates of u* and Rix via the least-square-fit are close to the designated values in Table 2, but not exactly the same. The discrepancy may result from the different estimates of viscosity between the Prandtl's mixing length model and the MY turbulence closure model. This comparisons between Nos. 1, 3 and 4, and Nos. 2, 5 and 6 show that under the same flux, the friction velocity u* reduces with the increase of Rix, causing the near-bottom current to weaken and the remote-bottom current to strengthen. This trend is enhanced as Rix increases. Fig. 9 shows the velocity profiles in unstratified flows under the impact of the upstream density gradient. Although an

Z. Zhang et al. / Estuarine, Coastal and Shelf Science 183 (2016) 163e178

173

Fig. 8. The comparison of the velocity profiles with different the downstream density gradients. The parameters are listed in Tables 1 and 3 for Nos. 1 to 6.

approximation is taken to the viscosity for obtaining a simple velocity profile formula of Eq. (16a), the agreement between the formula and the numerical results is still very good. Differing from the log-linear profile suggested by Burchard and Hetland (2010) and Winterwerp et al., 2006, Eq. (16a) takes into account the effect of Rix on viscosity via the function of F(Rix). Thus, the more correct boundary layer parameters are obtained in Table 3. Fig. 10 shows the velocity profiles in stratified flow under the impact of the density gradient. For comparison, the unstratified

flow results (No. 4 and No. 8) are also depicted. The results indicate that the stratification greatly strengthens the vertical velocity gradient compared to the that in unstratified flows under the same horizontal density gradient. By introducing and tuning the Monin-Obukov length L, the theoretical formula (15b) and (16b) can well fit the numerical results. Moreover, the estimate of L (shown in Table 3) through the velocity-profile-fit are close to that through the viscosity-profile-fit (shown in Table 2).

174

Z. Zhang et al. / Estuarine, Coastal and Shelf Science 183 (2016) 163e178

Table 3 The fitting parameters (friction velocity, horizontal Richardson number and the Monin-Obukov length) for the velocity profiles. No.

u* (m s1)

Rix

C

L (m)

1 2 3 4 5 6 7 8 9 10 11 12 13 14

0.023 0.010 0.022 0.022 0.009 0.007 0.025 0.028 0.012 0.015 0.019 0.013 0.024 0.021

0 0 2.5 8.3 10.1 43.2 1.2 3.1 5.5 11.4 8.3 11.0 3.0 3.1

0.52 0.21 0.40 0.30 0.09 0.04 0.60 0.66 0.26 0.29 0.18 0.11 0.70 0.68

e e e e e e e e e e 18 10 11 6

5. Discussion The classical log velocity profile is widely used in neutral flows due to its convenience (Zhang et al., 2016). However, that profile may introduce enormous errors in the estimate of friction velocity and roughness length when applied in baroclinic flows with a buoyancy gradient. So it is significant for log-fit to present the conditions under which it is permissible to neglect the buoyancy gradient effect. If it is permissible, the classical log profile is suggested for convenience. Otherwise, the resulting errors in estimate of the friction velocity should be given. This following example shows how to utilize the velocity profile (15a) for a given downstream density gradient to estimate the valid range of the log-fit estimate of the friction velocity within an acceptable margin of error. The velocity profiles can be written as the following general form

u¼

u*

k

XþC

(29) 6. Summary and conclusions

where X is a function of s, Rix and L, and C is a constant to be determined.

X ¼ X1 ¼ ln s

(30)

For neutral flow

For unstratified baroclinic flow X ¼ X2 ! pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ Rix s  1 þ 2 1 þ Rix s ¼ ln pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ Rix s þ 1

(31)

Given ui is the velocity magnitude at depth si, and Xi is the corresponding value according to si, then through the least-squarefit the ratio of the friction velocity estimated by the log law and that by the presented formula (15a) is written as

h¼

u*LOG ¼ u*NEW

where u*LOG is the velocity friction obtained by the log-fit and u*NEW is obtained by Eq. (15a). It shows that h depends on Rix and ui. Here we designate the average velocity ur¼0.19 m s1 and the friction velocity u*¼0.010 m s1, then the velocity profiles can be obtained through Eq. (15a) corresponding to Rix¼ 1, 5, 10, 30 respectively. Under extreme condition of Rix/0, Eq. (15a) approaches the log profile and h/1. Other values of h can be obtained through Eq. (32). The velocity profiles corresponding to Rix and h are shown in Fig. 11. The numerical results of Nos.2, 5 and 6 are also depicted in scatter plot form for comparison. Under the impact of the downstream density gradient, the log profile overestimates u* with increasing errors with the growth of Rix. For example, h¼3.03 means that the estimated u* by log-fit is more than three times of the designated one under the condition of Rix¼30. For practical use, the horizontal density gradient can be neglected if the condition 0.85h1.15 is satisfied, corresponding to less than 15% errors in the estimate of u*. As shown in Fig. 11, Rix¼1 corresponds to h¼1.14, which means that the errors in the estimate of u* by log-fit is about 14%. Therefore, the log profile can be used under the condition of Rix1 with the average flow ur¼0.19 m s1, according to 15% error threshold in u*. It is also shown in Fig. 11 that the simulated velocity profiles of No. 5 and No. 6 are outside of the range of 15% errors by log-fit with the value of Rix much more than 1. The error for No. 5 is about 95% and for No. 6 is more than 200%. In practical cases, the range of parameters typical of an estuary is listed as follows: ur¼0.10.5 m s1; u*¼0.010.1 m s1; vr 4 kg m4; D¼510 m. According to the definition vx ¼ ð1  5Þ  10 of Rix, the magnitude range of Rix is less than 5 in most cases. If the averaged flow velocity is near 0.19 m s1, the upper threshold of density gradient vvxr for log-fit is about 2104 kg m4, corresponding to Rix¼1. The upper threshold of density gradient increase with the growth of ur and u*. Hence, if given similar salinity conditions, the estuary with the larger flux can obtain more accurate estimate of u* by log-fit.

N N

PN

PN

PN

PN

PN

PN

i¼1 X1i ui  i¼1 X2i ui 

i¼1 X1i

i¼1 X2i

!

i¼1 ui

i¼1 ui

N !$ N

Mass transport and the structure of the water column in an estuary, especially in terms of the salinity flux, sediment, and contaminant transport, are a key concern for many researchers, given their importance in numerical modelling, for example. In order to estimate the flux and mass transport with enough accuracy, the structure and characteristics of the estuarine boundary layer should be well understood. The structure of the estuarine boundary layer is influenced by various processes and multiple dynamic factors, among which the density gradient plays an indispensable role in the circulation formulation. However, independent of tidal acceleration, wind forcing, bathymetry and Coriolis force, the effect of the density gradient on velocity profiles is still ubiquitous, although its contribution to estuarine circulation has been identified for a long time. The first contribution of the work is the derivation of the shear

PN

PN

PN

PN

PN

PN

2 i¼1 X2i  2 i¼1 X1i 

i¼1 X2i

i¼1 X1i

!

i¼1 X2i

i¼1 X1i

!

(32)

Z. Zhang et al. / Estuarine, Coastal and Shelf Science 183 (2016) 163e178

175

Fig. 9. The comparison of the velocity profiles with different upstream density gradients. The parameters are listed in Tables 1 and 3 for Nos. 7 to 10.

stress distribution theoretically through simplified hydrodynamic equation by neglecting acceleration term, advection term and wind shearing stress. Differing from that in steady barotropic flow, the shearing stress distribution in baroclinic flow deviates from the linear distribution but satisfies the parabolic distribution, which is upward for the downstream density gradient and downward for the upstream density gradient. To assess the deviation, the horizontal Richardson number, denoted by Rix, is introduced in the paper. The deviation is enhanced with the growth of Rix and vice versa with the decrease of Rix.

Another contribution is the derivation of the velocity profiles for the baroclinic flow under the impact of the density gradient. Applying the state-of-the-art two-equation turbulence closure model, we can obtain realistic eddy viscosity profiles and the velocity profiles in the numerical model. However, there is hardly a theoretical model to describe the velocity profile under the impact of density gradient due to the difficulty in describing the vertical eddy viscosity theoretically. To address this issue, the Prandtl's mixing length model is applied to obtain the eddy viscosity for balancing accuracy and complexity. Through

176

Z. Zhang et al. / Estuarine, Coastal and Shelf Science 183 (2016) 163e178

Fig. 10. The comparison of the velocity profiles in stratified flow (Nos. 11 to 14). The velocity profiles in unstratified flow (Nos. 4 and 8) are also depicted for reference.

comparison with the two-equation turbulence closure model results, the Prandtl's mixing length model obtains similar results with the numerical results. More complicated than the empirically parabolic distribution, the distribution of the eddy viscosity by the Prandtl's mixing length model can be divided into two situations: one is under the impact of upstream density gradient, and the other is under the impact of downstream density gradient. For the former situation, the parabolic distribution should be multiplied by a correction factor, as in Eq. (12a). But for the latter situation, a more complex B-shape viscosity profile was derived in Eq. (12b), with an additional zero-viscosity point in a

depth between s¼0 and s¼1. For simplicity, an amended parabolic distribution by multiplying a correction function, as in Eq. (12c), was raised in place of Eq. (12b) to obtain the velocity profiles. The comparison in tests demonstrates a good agreement between the theoretical and numerical velocity profiles. Thus, this result indicates that the velocity profiles are not as sensitive as the viscosity and the fit errors in viscosity do not have a significant effect on the velocity profiles with proper parameters. For the stably stratified flow, the Monin-Obukov length L is introduced, following the work of some scholars (Anwar, 1983; Sanford et al., 1991). The difference is that those scholars focused

Z. Zhang et al. / Estuarine, Coastal and Shelf Science 183 (2016) 163e178

177

Fig. 11. The log-fit error for different Rix values under the flux ur¼0.19 m s1. Several kinds of lines represent the velocity profiles by Eq. (15a) with different h and Rix. The scatters represent the numerical results under the same condition.

on vertical buoyancy and did not take into account the horizontal density gradient. Here, we concentrate on the influence of the stratification on the flow with a horizontal density gradient. The theoretical and numerical results both demonstrate that the viscosity and the bottom shearing stress are reduced by stratification because the vertical buoyancy restrains the turbulence and vertical mixing. Therefore, the main effect of stratification on the velocity profiles is that the deviation from the log profiles by the horizontal density gradient is significantly enhanced. The presented formula including the effect of horizontal density and stratification satisfactorily fits the numerical results. In numerous estuarine and coastal cases, the log-fit correlation of the velocity is acceptable, but the estimates of the friction velocity may deviate greatly from the proper values (Collins et al., 1998). The formula described here provides an approach for evaluating the errors of friction velocity caused by the classical log-fit in baroclinic flows, thereby defining a valid range of conditions within a reasonable margin of around 15% relative error, when a log-fit is applied. Acknowledgements We would like to thank Nanjing Hydraulic Research Institute for the technical support. The study was financially supported by the National Natural Science Foundation of China (Grant No. 41306078; 41301414; 51309159) and by Geographic Information research project of JPBSMG (Grant No. JSCHKY201504). References Adams, C.E., Weatherly, G.L., 1981. Some effects of suspended sediment stratification on an oceanic bottom boundary layer. J. Geophys. Res. 86, 4161e4172. Anwar, H.O., 1983. Turbulence measurements in stratified and well-mixed estuarine flows. Estuar. Coast. Shelf Sci. 17, 243e260. Burchard, H., Hetland, R.D., 2010. Quantifying the contributions of tidal straining and gravitational circulation to residual circulation in periodically stratified tidal estuaries. J. Phys. Oceanogr. 40, 1243e1262. Casulli, V., 1990. Semi-implicit finite difference methods for the two-dimensional shallow water equations. J. Comput. Phys. 86, 56e74. Casulli, V., Cattani, E., 1993. Stability, accuracy and efficiency of a semi-implicit

method for three-dimensional shallow water flow. Comput. Math. Appl. 27 (4), 99e112. Chatwin, P.C., 1976. Some remarks on the maintenance of the salinity distribution in estuaries. Estuar. Coast. Shelf Sci. 4, 555e566. Collins, M.B., Ke, X., Gao, S., 1998. Tidally-induced flow structure over intertidal flats. Estuarine. Coast. Shelf Sci. 46, 233e250. Friedrichs, C.T., Wright, L.D., 1997. Sensitivity of bottom stress and bottom roughness estimates to density stratification, Eckernforde Bay, southern Baltic Sea. J. Geophys. Res. 102 (C3), 5721e5732. Geyer, W.R., 1997. Influence of wind on dynamics and flushing of shallow estuaries. Estuar. Coast. Shelf Sci. 44, 713e722. Geyer, W.R., Trowbridge, J.H., Bowen, M.M., 2000. The dynamics of a partially mixed estuary. J. Phys. Oceanogr. 30, 2035e2048. Hansen, D.V., Rattray, M., 1965. Gravitational circulation in straits and estuaries. J. Mar. Res. 23, 104e122. Herrmann, M.J., Madsen, O.S., 2007. Effect of stratification due to suspended sand on velocity and concentration distribution in unidirectional flows. J. Geophys. Res. 112, C02006. http://dx.doi.org/10.1029/2006JC003569. Hetland, R.D., Geyer, W.R., 2004. An idealized study of the structure of long, partially mixed estuaries. J. Phys. Oceanogr. 34, 2677e2691. Kundu, P.K., 1990. Turbulence. chap. 12. In: Fluid Mechanics. Academic, San Diego, Calif, pp. 416e473. MacCready, P., 2004. Toward a unified theory of tidally-averaged estuarine salinity structure. Estuaries 27, 561e570. Mazumder, B.S., Ghoshal, K., 2006. Velocity and concentration profiles in uniform sediment-laden flow. Appl. Math. Model. 30, 164e176. Mellor, G.L., Yamada, T., 1982. Development of a turbulence closure model for geophysical fluid problems. Rev. Geophys. Space Phys. 20, 851e875. Monin, A.S., Yaglom, A.M., 1971. Application of dimensionless reasoning to turbulence in a stratified medium. In: Statistical Fluid Mechanics, Vol,1. MIT Press, Cambridge, Mass, pp. 425e442 section 7.2. Prandle, D., 2004. Saline intrusion in partially mixed estuaries. Estuar. Coast. Shelf Sci. 59, 385e397. Pritchard, D.W., 1952. Salinity distribution and circulation in the Chesapeake Bay estuarine system. J. Mar. Res. 11, 106e123. Pritchard, D.W., 1954. A study of the salt balance in a coastal plain estuary. J. Mar. Res. 13, 133e144. Sanford, L.P., Panageotou, W., Halka, J.P., 1991. Tidal resuspension of sediments in northern Chesapeake Bay. Mar. Geol. 97, 87e10. Simpson, J.H., Burchard, H., Fisher, N.R., Rippeth, T.P., 2002. The semidiurnal cycle of dissipation in a ROFI: model-measurement comparisons. Cont. Shelf Resour. 22, 1615e1628. Simpson, J.H., Williams, E., Brasseur, L.H., Brubaker, J.M., 2005. The impact of tidal straining on the cycle of turbulence in a partially stratified estuary. Cont. Shelf Resour. 25, 51e64. Soulsby, R.L., Wainwright, B.L.S.A., 1987. A criterion for the effect of suspended sediment on near-bottom velocity profiles. J. Hydraulic Res. 25 (3), 341e356. Styles, R., Glenn, S.M., 2000. Modelling stratified wave and current bottom boundary layers on the continental shelf. J. Geophys. Res. 105 (C10),

178

Z. Zhang et al. / Estuarine, Coastal and Shelf Science 183 (2016) 163e178

24119e24139. Turner, J.S., 1973. Turbulent shear flows in a stratified fluid. chap.5. In: Buoyancy Effects in Fluids. Cambridge University Press, New York, pp. 127e165. Vinita, J., Shivaprasad, A., Revichandran, C., Manoj, N.T., Muraleedharan, K.R., Binzy, J., 2015. Salinity response to seasonal runoff in a complex estuarine system (Cochin estuary, west coast of India). J. Coast. Res. 31 (4), 869e878.

Webb, E.K., 1970. Profile relationships: the log-linear range, and extension into strong stability. Q. J. R. Meteorol. Soc. 96, 67e90. Winterwerp, J.C., Wang, Z.B., van der Kaaij, T., et al., 2006. Flow velocity profiles in the Lower Scheldt estuary. Ocean. Dyn. 56, 284e294. Zhang, Z., Song, Z.Y., Zhang, D., Guo, F., Hu, D., 2016. Application of power law to the currents in a tidal channel. J. Oceanogr. 72, 589e600.