Evaluation of entrainment velocity induced by wind stress in a two-layer system

This study aims to reveal the factors most useful for evaluating the influence of wind on entrainment from the lower to upper layer in a two-layer system. Lake Abashiri, which is a typical brackish lake, is chosen as a study area because a distinct two-layer system exists due to salt-wedge intrusion from the ocean. To evaluate entrainment velocity between two layers, a 3D numerical computation is applied, showing good agreement with field observations. Computations suggest that the entrainment velocity estimated using the Richardson number is smaller than that from field observations. Thus, other factors are computed suggesting that use of the Lake number is more effective in estimating entrainment velocity due to internal wave breaking and interfacial fluctuation than use of the Effective Wedderburn number.


INTRODUCTION
In lakes and enclosed bays stratification is likely during summer due to strong solar radiation and a supply of fresh water from rivers forming a warm, fresh layer of lower density at the surface (Nakayama et al., 2005;Antenucci et al., 2000). Since strong stratification suppresses mixing through the density interface between the upper and lower layers, residence times tend to increase. Longer residence times can result in anoxia or hypoxia in the lower layer due to the consumption of dissolved oxygen during degradation of accumulated particulate organic matter. Anoxia and hypoxia have been revealed to have significant effects on water quality and thus the ecological function of enclosed water bodies (Okada et al., 2007;Okada et al., 2009;Wolanski et al., 2006). As a result many previous studies have made an attempt to clarify the mechanisms for formation and deformation of the density interface (Antenucci et al., 2001;Antenucci et al., 2003;Hodges et al., 2000a). Mixing due to turbulent energy driven by wind and flow has been shown to influence formation of density interface using 3D numerical models. Modeling has also demonstrated that the breaking of internal waves over topography has an especially large influence on mixing at the density interface (Hodges et al., 2000b).
In lakes and the ocean it has been shown that climate change may influence water quality and ecological systems due to a variety of factors, such as changes in rainfall patterns, increases in air temperatures and sea-level rise. Although 3D numerical models have clarified the mechanisms of formation and breakdown of the density interface, it is difficult to extend 3D numerical computations for the analysis of long-term scenarios such as climate change. In contrast to 3D numerical models, the runtime cost of a 1D vertical numerical model is very cheap, which is useful for scenario analysis of climate change (Yamashiki et al., 2010;. As a 1D vertical numerical model can take into account stratification effects only, the effect of internal waves breaking over topography has to be modeled using entrainment velocity. Entrainment velocity is considered to be one of the most significant factors controlling the level of the density interface and stratification. Previous studies have revealed that internal wave breaking and interfacial fluctuation enhance entrainment velocity in a two-layer system, though entrainment velocity is generally estimated by using the Richardson number only (Nakayama and Imberger, 2010). Therefore, this study aims to identify the factors which determine the influence of internal wave breaking on entrainment in a twolayer system, and thereby enhance the potential applicability of 1D models for climate change analysis.

METHOD
A two-layer system is the most appropriate in which to easily evaluate the influence of internal wave breaking over topography on entrainment at the density interface. Lake Abashiri, located along the Okhotsk Sea in northern Hokkaido, Japan, was thus selected as the study site because of the distinct density interface induced by sea-water intrusion from the sea to the lake (Figure 1). Lake Abashiri is a typical brackish lake and has water quality issues, including anoxia. Nakayama et al. (2009) and Maruya et al. (2010) investigated the influence of wind, river inflow and tide on the deformation of the density interface, which demonstrated the robustness of a 3D numerical model in this system. We thus made an attempt to analyze the stratified flow field of October in 2006, in which the applicability of a 3D numerical model has already been confirmed. The total water input from rivers was 1.8 × 10 8 m 3 from the 7 th to 16 th of October, which almost equals the total volume of Lake Abashiri, 2.3 × 10 8 m 3 (Figure 2). The peak discharge occurred at 18:00 on the 8 th of October, at around 800 m 3 s −1 . As the maximum wind speed was high on the 8 th of October, at about 21 m s −1 , wind is expected to play a great role in deformation of the density interface. The 3D numerical model, Fantom3D, was applied to evaluate changes in the density interface from the 1 st to 30 th of October in 2006. Fantom3D is the object programming model for the analysis of environmental fluid mechanics (Shintani and Nakayma, 2009). Since the thickness of the density interface in Lake Abashiri is about 0.5 m, a mesh size decreasing from 1 m close to water surface to 0.1 m around density interface was used over the total depth of 17 m, resulting in a total of 100 vertical meshes. The size of horizontal meshes was 200 m in the x and y plane, giving a total number of meshes of 55 × 27 × 100. Parallel computing was employed which resulted in the computational domain being divided into 4 separate domains of 28 × 14 × 100 meshes. The computational time step was set to 40 sec. Meteorological boundary conditions were given using Automated Meteorological Data Acquisition System (AMeDAS) observational data, and river discharge was evaluated by combining field observations and theoretically derived values. To analyze the density interface regarding topography effect accurately and to prevent numerical diffusion from occurring, a z coordinate system was employed. Large eddy simulation was applied as the turbulent model, the advection momentum scheme was rational cubic-polynomial interpolation, the advection scalar scheme was ultimate quickest, and the advection turbulence scheme was upwind.
The validity of Fantom3D in Lake Abashiri from the 1 st to 30 th of October in 2006 was investigated using field observations from st.1, st.2 and st.3. At each station water temperature and salinity were measured vertically at an interval of 0.1 m every 1 hour ( Figure 1). As st.1 is located in the center of the lake, it can capture the amplitude of oscillation of the density interface due to internal Poincare waves. Significant amplitude of oscillation was not found to exist, which suggests that internal Poincare waves were not predominant during the analytical period. The density interface was revealed to decrease by about 2 m, from about 5 m to 7 m, which was confirmed in both the field observations and the numerical model (period (a) in Figure  3a and Figure 3b, please find the details in Shintani and Nakayama, 2009). As the average wind speed from 16:00 on the 7 th of October to 9:00 on the 9 th of October was high, at more than 10 m s −1 , wind is expected to have played a great role in decreasing the level of the density interface level and, correspondingly, increasing the entrainment velocity in this two-layer system. On the other hand, there is the potential for significant influence from internal Kelvin waves, with field observations and numerical modeling results showing the occurrence of large amplitude internal Kelvin waves at st.2 and st.3. Thus, comparison with field data demonstrates the robustness of the numerical computation model. It should be noted that the density interface is defined to be the level where salinity is 10 compared to a mean value of 0 and 21 in the upper and lower layers, respectively.
Numerical computations revealed that winds greater than 10 m s −1 drove upwelling of the lower layer water at 16:00 on the 7 th of October. This actually caused a 'blue tide' of hypoxic water in the shallows (period (a) in Figure 3e). At 1:00 on the 8 th of October, winds increased above 10 m s −1 , which caused upwelling and the lower layer water was found to rise to the water surface (S1. and S2a.). This upwelling induced vertical second mode internal waves, which is found from the thickening of the density interface after winds dropped below 10 m s −1 (S1. and S2b.). At 3:00 on the 8 th , the upwelling was again confirmed to be induced by strong winds greater than 10 m s −1 , and strong mixing over the slope was revealed to occur due to the breaking of the density interface (S1. and S2c.). In contrast to the upwelling, the interface was confirmed to be kept horizontal when wind speed is weak, 1.9 m s −1 , at 8:00 on the 7 th (S2d.). Therefore, these two strong wind events are considered to represent the dominant phenomena controlling the level of the density interface.

RESULTS
To evaluate the influence of interface fluctuation and internal wave breaking on the level of the density interface, entrainment velocity was computed during these two strong wind events from the change in the density interface level. Firstly, the volume of the upper layer was evaluated as the area where density is less than the mean density between the upper and lower layers. Then, the entrainment velocity was computed from the temporal change in the volume of the upper layer by being divided by area (Figure 4). In previous studies, entrainment velocity has been found to be a function of the Richardson number by using the upper layer mean velocity. Thus, entrainment velocity estimated from the Richardson number was computed by using equations (1) to (3) (Figure 4) (Price et al., 1980). The actual entrainment velocity was found to be much larger than that estimated using the Richardson number at around 2:00 and 13:00 on the 8 th of October (period (b) and (d) in Figure 4). (1) where, R i is the Richardson number, ε is the density difference ratio, g is the gravitational acceleration, l is the surface layer thickness, u 1 is the representative speed of surface layer, w e is the entrainment velocity, w h is the variation velocity of density interface, and h is the level of density interface.

DISCUSSION
Underestimation of the entrainment velocity using the Richardson number may indicate the occurrence of internal waves due to strong winds and resultant enhancement of entrainment velocity due to breaking of the internal waves over a slope. Therefore, it is necessary to introduce new factors which can allow evaluation of the entrainment velocity in addition to the Richardson number. We thus made an attempt to apply the Lake number and Effective Wedderburn number, which can demonstrate the stability of the stratified flow field ( Figure 5). The Lake number is defined by equation (4) (Imberger J., and J. C. Patterson, 1990). The Effective Wedderburn number is described in Shintani et al. (2010). (4) where, : density stratification of potential energy, h m is the surface layer thickness, ρ 0 is the density of water, u * is the wind friction velocity, H is the water depth, A s is the surface area of the lake A(H), and z g is the height to the center of volume of the lake.
Since both the Lake number and Effective Wedderburn number were less than 1 from 0:00 to 3:00 and from 11:00 to 14:00 on the 8 th of October, the stratified water body became unstable resulting in actual upwelling (period (b) and (d) in Figure 4 and period (b) and (d) in Figure 5). The unstable conditions caused instability of the density interface, and drove larger entrainment velocities (period (b) and (d) in Figure 4). Furthermore, the density interface was found to fluctuate greatly between these two entrainment velocity peaks, and internal waves were induced (period (c) in Figure 4 and period (c) in Figure 5 and square of S3.). The breaking of the internal waves is considered to have caused mixing over the slope and increased entrainment velocity in period (c) in Figure 4. In addition, as the Lake number was less than 1 from 21:00 on the 7 th to 4:00 on the 9 th of October, the stratified water body became unstable and internal waves were induced, which resulted in relatively larger entrainment velocities (period (e) in Figure 4).
Therefore, it was found that the Lake number and Effective Wedderburn number are useful for the evaluation of entrainment velocity caused by internal wave breaking and interface fluctuations, which cannot be evaluated by using the Richardson number. To clarify which parameter is more appropriate for representing entrainment velocity, entrainment velocity due to internal wave breaking was plotted against the Lake number, and the Effective Wedderburn number was also compared with entrainment velocity (Figure 6). Although there are some anomalies shown by the solid-line square, it was found that the larger the Lake number, the smaller the entrainment velocity due to internal wave breaking and interface fluctuation ( Figure  6a). In contrast, the Effective Wedderburn number does not seem to show a clear relationship with entrainment velocity due to internal wave breaking and interface fluctuation though a peak entrainment coefficient appears to exist around the Effective Wedderburn number, 1.0 (Figure 6b). The Lake number and Effective Wedderburn number represent instability of the stratified flow field and the potential for upwelling, respectively. Since internal wave breaking and interface fluctuation are likely to be induced by instability of the stratified flow field, the Lake number is considered to be more appropriate for the evaluation of entrainment velocity compared to the Effective Wedderburn number.

ACKNOWLEDGEMENT
This work has been supported by the Japan Society for the Promotion of Science.