Sensitivities of Bottom Stress Estimation to Sediment Stratification in a Tidal Coastal Bottom Boundary Layer

The bottom friction velocity (U*), which controls seabed erosion and deposition, plays a critical role in sediment transport in tidal coastal bottom boundary layers. Approaches have been proposed to calculate U*, including the log profile (LP) estimation, the direct covariance (COV) measurement, and the turbulent kinetic energy (TKE) method. However, the LP method assumes homogeneous flow and the effects of stratification need to be taken into account. Here, field investigations of hydrodynamics and sediment dynamics were carried out on the Jiangsu Coast, China. Two acoustic Doppler velocimeters (ADV) velocity measurements at 0.2 and 1 m above the seabed have been used to estimate U*, based on the aforementioned three methods. The COV and TKE methods provided reasonable estimations of U*, while a pronounced overestimation was identified when using the LP method. This overestimation can be attributed to the stratification effects associated with the vertical suspended sediment concentration (SSC) gradient near the bottom. Then, three models were utilized to correct the overestimation, in which the gradient/flux Richardson number was modified with empirical constants α, β, and A to parameterize the stratification effects in the logarithmic velocity distribution. The values of α, β, and A derived from the observation are smaller than the results from previous investigations. These modified logarithmic velocity distribution models can be applied in numerical simulations when sediment stratification is important.


Introduction
Coastal boundary layer processes (e.g., near-bottom blending, energy and momentum transformation, and sediment resuspension and transport) play a fundamental role in the nature of marine sediment dynamics [1][2][3][4]. Bottom friction velocity (U * ) and/or current-induced bottom shear stress, which largely control seabed erosion and deposition, play a critical role in sediment resuspension and transport in tidal coastal bottom boundary layers [1,5].
Several methods to calculate U * are summarized in [6], including the log profile estimate (LP), the direct covariance (COV) measurement, and the turbulent kinetic energy (TKE) method, which are widely used [7][8][9][10][11][12]. The LP method is based on the van Karmen-Prandtl boundary model, and the velocity profile can be fitted by a logarithmic curve to derive U * . The shear stress is directly determined

Study Area
The Jiangsu Coast is adjacent to the South Yellow Sea. Tides in the study area are irregular semidiurnal and mesotidal, with a mean tidal range of 2-4 m, and the tidal bottom boundary layer is well developed [20]. The offshore area is dominated by longshore southward flood currents and northward ebb currents because of the southward propagation of the tidal wave. The local significant wave height and average wave period range from 0.15-2.22 m and 2.06-6.82 s, respectively [21]. Suspended sediment concentration (SSC) gradients (>0.3 kg m −4 ) appear in the local bottom boundary layers and the maximum SSC can be greater than 5 kg m −3 , suggesting a distinct sediment stratification [22]. Thus, the study area is suitable to investigate the sensitivities of bottom stress estimation to sediment stratification in a tidal coastal bottom boundary layer. The seabed sediments are characterized by silty sand in the study area [23].

Field Observations
A bottom-boundary-layer instrumented tripod was deployed to measure the hydrodynamics and sediment dynamics at an anchor station (31.617 • N, 120.654 • E) located outside the Jiangsu Coast, China, during January 18-19, 2015 ( Figure 1). The tidal average water depth at the site is about 7 m [22]. The station is located far from rivers with significant freshwater discharge (the annual average discharge is 66.6 m 3 s −1 ) [22], and the well-mixed water caused by strong tides has negligible vertical salinity variation.
Two acoustic Doppler velocimeters (ADV) were mounted at 0.2 and 1 m above seabed (asb), respectively, to measure three-dimensional and high-frequency (16 Hz) velocities for a period of 128 s in each burst interval of 15 min. At the same heights, i.e., 0.2 and 1 m asb, two optical backscatter sensors (OBS) also measured water turbidities at a sampling frequency of 1 Hz for a period of 30 s in each burst interval of 3 min. Besides, a Seabird Electronic SBE26 wave and tide recorder was installed at 1 m asb, measuring water depth and wave height and period at 4 Hz for 300 s within a 1 hour burst interval. The deployment period was 25 h for continuous observation, from 13:00 on January 18 to 14:00 on January 19. The deployment schemes are summarized in Table 1.
Turbidities from the optical backscatter sensors (OBS) in nephelometric turbidity units (NTU) were calibrated to suspended sediment concentration (SSC) (kg m −3 ) in the laboratory using the in situ suspended sediment samples. Seabed surface sediment samples, collected at the observation station, were stored in plastic bags and brought to the laboratory for grain size analysis using the laser diffraction method (Mastersize2000 granulometer, Malvern Ltd., Malvern, UK). situ suspended sediment samples. Seabed surface sediment samples, collected at the observation station, were stored in plastic bags and brought to the laboratory for grain size analysis using the laser diffraction method (Mastersize2000 granulometer, Malvern Ltd, Malvern, UK).

Figure 1.
Map of the study area. The black star denotes the location of the anchor station, the arrows denote the flood and ebb directions of the rectilinear tidal currents, and the positive direction is 10° anticlockwise from the north (ebb direction). Unfortunately, the update data are insufficient to identify depth contour lines in the study area. The LP method relies on the classical "law of the wall," which is expressed by the Karman-Prandtl equation for fully rough turbulent boundary layers: Figure 1. Map of the study area. The black star denotes the location of the anchor station, the arrows denote the flood and ebb directions of the rectilinear tidal currents, and the positive direction is 10 • anticlockwise from the north (ebb direction). Unfortunately, the update data are insufficient to identify depth contour lines in the study area.

LP Method
The LP method relies on the classical "law of the wall," which is expressed by the Karman-Prandtl equation for fully rough turbulent boundary layers: where U is the current velocity at a certain height, U * is bottom friction velocity, Z 0 is bed roughness length, and Z is the height above the bed. The Karman constant κ is approximately 0.41. In this study, we only considered the tidal-induced current bottom boundary layer, in which the thickness of the log layer reached 10 1 m [1,5]. Z 0 and U * can be obtained by fitting ln(Z) and U for a velocity profile in the log layer: However, in this study, two ADV velocity measurements at 0.2 and 1 m asb, which both fall in the tidal-induced log layer, could be used for calculation. You [24] extended the classical Karman-Prandtl equation to estimate the bed roughness length from mean velocities measured at two levels near the seabed over a certain period (e.g., several tidal cycles). Because Z 0 is primarily controlled by the bed sediment grain size and bedform geometry [25,26], Z 0 was assumed to be constant for several days, for example. The Karman-Prandtl formula can be converted to: where U 2 and U 1 are the burst averaged (the 128 s average velocity in each burst interval of 15 min) current speeds at the elevation Z 2 and Z 1 above seabed (asb), respectively, K is the slope with the least-squares fitting method from time series of U 1 and U 2 , and Z 0 is the averaged bed roughness length during the observation period. U * can be then expressed as: In this way, the measured U 1 at Z 1 and U 2 at Z 2 here were used to calculate U *.

COV Method
The current-related bed shear stress is determined from the Reynolds stress method [6,27,28]: where τ zx and τ zy are the Reynolds stress components along and across the main flow direction, respectively, ρ w is the density of seawater, u', v', and w' are the instantaneous turbulence components of the horizontal (along and across the main flow direction) and vertical velocities in a certain height (Z) above the bed, respectively, which can be described as the difference of the measured velocity and mean components according to the Reynolds decomposition (i.e., u = u + u ), and the overbar denotes the burst average. Thus: where τ z is the Reynolds stress at a certain height (Z) above the bed. The measured local shear stress is corrected for the effect of water depth variation [29,30]: where τ b is bottom shear stress at the bed, h is water depth, and z is referenced to the seabed.
Thus, both the calculated time series of Reynolds stress at 0.2 and 1 m asb in Equations (6)-(8) can be used to calculate time series of τ b , expressed as τ b0.2 and τ b1 , respectively. Their arithmetic mean is then: Finally, U * can be calculated as:

TKE Method
Previous studies have found that the mean ratio of shear stress to TKE is constant [31,32]. The shear stress can be calculated through TKE: where E is TKE, u', v', and w' are the turbulence components of the horizontal (along and across the mean flow direction) and vertical velocities, respectively. C is a proportionality constant representing the ratio of shear stress to TKE. Soulsby and Dyer [31] adopted C to be~0.20, Huntley [33] took C at 0.19, and Kim et al. [6] assumed C to be~0.21. Here, C was chosen to be~0.20 (more on this later). Like the COV method, U * can be calculated from Equations (9)-(11).

Stratification Effect
To predict the log velocity distribution in sediment stratified flows, Karman's constant (κ) is usually modified by introducing the Richardson number, which represents the impact of stratification. Thus, the velocity profile reflecting the stratification effects can be expressed as: where κ is the modified Karman constant. Fischer et al. [17] modified κ to: In Smith and McLean [18], κ was corrected to: Adams and Weatherly [19] adjusted the constant as: where α, β and A are empirical constants. Fischer et al. [17] suggested that α is equal to 3.33, and thereafter Hamblin [34] used this value in their numerical models. β was proposed to be 4.7 ± 0.5 by Businger et al. [35], which was adopted in the modeling work of Smith and McLean [18]. The value of A was suggested to be 5.5 [19], which was then applied in the models of [36]. R f and R i are the flux Richardson number and gradient Richardson number, respectively, which can be derived from the level-2 closure scheme [37]: For the case of stratification induced by a suspended sediment gradient, the density of the seawater with suspended sediment (ρ) can be calculated by: where ρ s is the specific density of sediment (2650 kg m −3 in this study) and C is the concentration of suspended sediment. The measured velocity and SSC at 0.2 and 1.0 m asb were used to calculate R f and R i at the mean location between the two sensor locations, 0.6 m asb. According to Equations (1), (15)- (17), LP method results can be modified by Equations (21)-(23) to account for the effects of stratification: where U * m is the modified friction velocity, reflecting the effects of stratification. During the observation, the time series of the modified U * m (Equations (21)- (23)) was fitted linearly with that of U * estimated by the COV method, forcing the intercept to be zero. When the slope in the regression is equal to 1, reflecting a statistical equivalence of U * m and U * estimated by the COV method, the empirical constants α, β, and A can be determined.

Observation
During the 25-hour bottom tripod measurement, the water depth varied from 5.19 to 8.42 m, and the average depth was 7.15 m. This area was characterized by weak wave activity, the significant wave height ranging from 0.08 to 0.14 m, and the mean period varying from 2.17 to 2.95 s ( Figure 2). The tidal currents were rectilinear, and velocities in the main direction at 0.2 and 1 m asb had considerable tidal variations, ranging from −0.52 to 0.49 m s −1 and from −0.72 to 0.64 m s −1 , respectively, the positive direction being 10 • counterclockwise from the north. The strong tides and weak waves denote the tidal-dominated hydrodynamic environment during the observation. SSC at 0.2 and 1 m asb changed from 0.94 to 1.40 kg m −3 and from 0.77 to 1.25 kg m −3 , and their mean values were 1.09 and 0.97 kg m −3 , respectively. Accordingly, SSC gradients near the seabed were fairly large (0.07-0.45 kg m −4 , with an average of 0.14 kg m −4 ).
The Reynolds stress (calculated from the COV method) in the main flow direction at 0.2 and 1 m asb varied from −1.23 to 1.32 Pa and from −1.15 to 1.29 Pa, respectively, and the positive direction is the same as the current velocity. The time series of Reynolds shear stress at 0.2 m asb were slightly larger than at 1 m asb, because the Reynolds stress increase linearly downward from the surface to the seabed (Equation (9)). The gradient Richardson number and the flux Richardson number were 0.01-14.79 and 0.01-0.25, respectively, indicating pronounced sediment stratification (Figure 2). The area is characterized by sandy bed with a mean grain size of 3.9 ϕ, and the contents of sand, silt, and clay were 65.50%, 32.35%, and 2.15%, respectively. The threshold friction velocity at the seabed is 0.011 m s −1 .

Estimation of Friction Velocity
All the three methods (Section 3.2) were used to estimate the bottom shear stress. Firstly, the LP method was applied. The burst averaged current speeds at 0.2 (U 0.2 ) and 1 (U 1 ) m asb were fitted, which shows a linear relationship (R 2 = 0.98) (Figure 3). The average ratio of U 1 to U 0.2 was 1.32, and thus the averaged Z 0 was 1.34 × 10 −3 m based on Equation (3). The friction velocity U *0.2 and U *1 were calculated by U 0.2 and U 1 , respectively, using Equation (5), and their time series are approximately coincident (Figure 4a). The 1:1 relationship between U *0.2 and U *1 indicates that U * estimated by Equation (5) is irrelevant to the elevation of current velocity (Figure 4b).
The relationship between the bottom shear stress estimated from the COV method (τ COV ) and the turbulent kinetic energy (TKE) is shown in Figure 5. They are highly correlated (R 2 = 0.96), and the ratio between τ COV and TKE is 0.20, which is in agreement with the value of the constant C in Equation (12) (the ratio of the shear stress estimated from the TKE method and TKE) used in this study. Therefore, the τ COV is consistent with the bottom shear stress estimated from the TKE method.  (Figure 4a). The 1:1 relationship between U*0.2 and U*1 indicates that U* estimated by Equation (5) is irrelevant to the elevation of current velocity (Figure 4b).
The relationship between the bottom shear stress estimated from the COV method (τCOV) and the turbulent kinetic energy (TKE) is shown in Figure 5. They are highly correlated (R 2 = 0.96), and the ratio between τCOV and TKE is 0.20, which is in agreement with the value of the constant C in Equation (12) (the ratio of the shear stress estimated from the TKE method and TKE) used in this study. Therefore, the τCOV is consistent with the bottom shear stress estimated from the TKE method.         The time series of U* estimated from the method of LP (based on the tidal averaged Z0 of 1.34 × 10 −3 m and the burst averaged velocity measured at 1 m asb), COV (based on the high-frequency velocity measured at 0.2 and 1 m asb), and TKE (based on the high-frequency velocity measured at 0.2 and 1 m asb) are shown in Figure 6a. The COV and TKE methods have consistent results, and the LP method significantly overestimated U*. The time series of U * estimated from the method of LP (based on the tidal averaged Z 0 of 1.34 × 10 −3 m and the burst averaged velocity measured at 1 m asb), COV (based on the high-frequency velocity measured at 0.2 and 1 m asb), and TKE (based on the high-frequency velocity measured at 0.2 and 1 m asb) are shown in Figure 6a. The COV and TKE methods have consistent results, and the LP method significantly overestimated U * .

Correction of the LP Method with Stratification Effects
The linear fit shows that U * from the LP method is 19% larger than U * from the COV method (Figure 7a), suggesting a pronounced overestimation by the LP method. The overestimate of U * from the LP method was generally consistent with the effects of sediment stratification [6,8].
In order to correct the overestimated U * from the LP method, the three methods described in Section 3.3 were applied, in which the Richardson number (Figure 2f,g) and empirical constants (α, β, and A) were adjusted to represent the stratification effects. Figure 7b-d show the adjustments of the U * from the LP method to that from the COV method (1:1 regression lines). The empirical constants α in Fischer et al.'s [17] model (Equation (21)) was estimated to be 3.16, while the empirical constants β in Smith and McLean' s [18] model (Equation (22)) was 2.99, and the empirical constants A in Adams and Weatherly's [19] model (Equation (23)) was 3.98. The time series of the corrected estimations of U * from Equations (21)-(23) are shown in Figure 6b, suggesting reasonable agreement with the U * from the COV method.

Correction of the LP Method with Stratification Effects
The linear fit shows that U* from the LP method is 19% larger than U* from the COV method (Figure 7a), suggesting a pronounced overestimation by the LP method. The overestimate of U* from the LP method was generally consistent with the effects of sediment stratification [6,8].
In order to correct the overestimated U* from the LP method, the three methods described in Section 3.3 were applied, in which the Richardson number (Figure 2f,g) and empirical constants (α, β, and A) were adjusted to represent the stratification effects. Figure 7b-d show the adjustments of the U* from the LP method to that from the COV method (1:1 regression lines). The empirical constants α in Fischer et al.'s [17] model (Equation (21)) was estimated to be 3.16, while the empirical constants β in Smith and McLean' s [18] model (Equation (22)) was 2.99, and the empirical constants A in Adams and Weatherly's [19] model (Equation (23)) was 3.98. The time series of the corrected estimations of U* from Equations (21)(22)(23) are shown in Figure 6b, suggesting reasonable agreement with the U* from the COV method.

Discussion
The COV method is suggested to provide unbiased estimates of U* [6,8,38]. The ratio of U* estimated from the COV method and TKE is 0.20, suggesting the TKE method performs well with the constant C = 0.20. This result is consistent with the findings of Soulsby and Dyer [31].
Bed roughness height Z0 is controlled by the inherent characteristics of the seabed (e.g. sediment grain size and bedform geometry) [13,39], and these characteristics can be assumed to remain essentially unchanged, thus Z0 should remain relatively constant throughout the observation period [25,40,41]. For a sandy bed, including the present study area, a typical value of 0.4 × 10 −3 m is suggested [1,13]. However, the Z0 derived from Equations (3) and (4) was 1.34 × 10 −3 m, significantly larger than the typical value. This is caused by stratification associated with the SSC gradient near the bottom.
The failure of the LP method to predict U* in a tidal bottom boundary layer is most likely due to the density stratification. Temperature, salinity, and SSC gradient near the bottom are the most direct mechanisms to cause density stratification. However, in the shallow well-mixed coastal water,

Discussion
The COV method is suggested to provide unbiased estimates of U * [6,8,38] . The ratio of U * estimated from the COV method and TKE is 0.20, suggesting the TKE method performs well with the constant C = 0.20. This result is consistent with the findings of Soulsby and Dyer [31].
Bed roughness height Z 0 is controlled by the inherent characteristics of the seabed (e.g., sediment grain size and bedform geometry) [13,39], and these characteristics can be assumed to remain essentially unchanged, thus Z 0 should remain relatively constant throughout the observation period [25,40,41]. For a sandy bed, including the present study area, a typical value of 0.4 × 10 −3 m is suggested [1,13]. However, the Z 0 derived from Equations (3) and (4) was 1.34 × 10 −3 m, significantly larger than the typical value. This is caused by stratification associated with the SSC gradient near the bottom.
The failure of the LP method to predict U * in a tidal bottom boundary layer is most likely due to the density stratification. Temperature, salinity, and SSC gradient near the bottom are the most direct mechanisms to cause density stratification. However, in the shallow well-mixed coastal water, temperature and salinity gradients are too small to exert impacts on stratification [42]. Thus, in this study, stratification arising from the pronounced SSC gradient near the bottom (Figure 2) is the main reason for the overestimation of U * from the LP method.
Three modified logarithmic velocity distribution models (Equations (15)- (17)) can be well applied to correct the overestimation with the present observational data, in which the gradient/flux Richardson number is modified with empirical constants α, β, and A to parameterize the stratification effects. The value of α in this study (Figure 7b) is 3.16, which is slightly smaller than that of 3.33 demonstrated by Fischer et al. [17]. The constant β is 2.99 (Figure 7c), which is~2/3 of the original value suggested by Businger et al. [35] and Smith and McLean [18]. However, the value selected for A is 4.7 [43], 5.2 [44], 5.5 [19], and 6.8 to 14.7 [45], while the present result is 3.98, which is significantly smaller compared to previous research.
The modified logarithmic velocity distribution models (Equations (15)- (17)) have been widely used in numerical modeling studies [34,36,46,47]. The present field observations provide new evidence for a slightly lower value of α, and pronounced lower values for β and A. The reason for these differences is not clear. It is hypothesized that the Fischer et al.'s [17] model with the gradient Richardson number (R i ) is more reliable. This problem will be investigated in future studies. At the present, it is suggested that field calibrations are essential for numerical models.

Conclusions
Two ADV velocities at 0.2 and 1 m asb observed in the tidal bottom boundary layers outside the Jiangsu Coast, China, were analyzed with the LP, COV, and TKE methods. The COV and TKE methods provided a reasonable estimation of U * , and pronounced overestimation was identified with the LP method. This overestimation can be attributed to the stratification effects associated with the vertical SSC gradient near the bottom. The modified logarithmic velocity distribution models could correct the overestimation with the empirical constant α, β, and A. The values of α, β, and A derived from the observation data were 3.16, 2.99, and 3.98, respectively, which are smaller than previous estimates.