Large-eddy simulation of the in fluence of a wavy lower boundary on the turbulence kinetic energy budget redistribution*

LU Zongze (陸宗澤) , FAN Wei (樊偉) , , LI Shuang (李爽) , GE Jianzhong (葛建忠)

1 Ocean College, Zhejiang University, Zhoushan 316000, China

2 State Key Laboratory of Estuarine and Coastal Research, East China Normal University, Shanghai 200062, China

Abstract Oceanic turbulence plays an important role in coastal flow. However, as the effect of an uneven lower boundary on the adjacent turbulence is still not well understood, we explore the mechanics of nearshore turbulence with a turbulence-resolving numerical model known as a large-eddy-simulation model for an idealized scenario in a coastal region for which the lower boundary is a solid sinusoidal wave.The numerical simulation demonstrates how the mechanical energy of the current is transferred into local turbulence mixing, and shows the changes in turbulent intensity over the continuous phase change of the lower topography. The strongest turbulent kinetic energy is concentrated above the trough of the wavy surface. The turbulence mixing is mainly generated by the shear forces; the magnitude of shear production has a local maximum over the crest of the seabed topography, and there is an asymmetry in the shear production between the leeward and windward slopes. The numerical results are consistent with results from laboratory experiments. Our analysis provides an important insight into the mechanism of turbulent kinetic energy production and development.

Keyword: large-eddy simulation; wavy lower boundary; oceanic turbulence; nearshore

1 INTRODUCTION

Coastal water is always well mixed compared with deep ocean water. While relatively well-mixed coastal water promotes the generation of coastal flow (Li et al., 2005), coastal mixing does not remove the vertical current pro file, implying that turbulence mixing still requires appropriate modeling. Models based on Reynolds-averaged Navier-Stokes (RANS) equations were used to describe coastal vertical mixing in early studies, which either parameterized the ocean turbulence as a bulk quantity (Large et al., 1994;Large and Gent, 1999; McWilliams and Sullivan,2000; Smyth et al., 2002; Wijesekera et al., 2003) or used speci fic parameterizations of terms of the differential equations (Mellor and Yamada, 1982;Wilcox, 1988; Umlauf and Burchard, 2003). While RANS models assume that the turbulence can be characterized by the background low-frequency flow,or by the degree of convective instability, this assumption results in excessive noise and does not capture the nature of the turbulence. In response, the coastal turbulence study of Zeng et al. (2008)proposed a turbulence-resolving model to correct the uncertainty of RANS models.

Large-eddy simulations (LES) have produced successful numerical solutions of coastal turbulence(Li et al., 2013, 2016; Walker et al., 2016). For example, for small-scale unstable flow above sand grains, Chang and Scotti (2004) compared a RANS model with an LES model to show similar results for only the vertical variation of the streamwise velocity component, with obvious differences in the vertical velocity component, turbulence kinetic energy (TKE)budget, and other higher-order quantities. The study concluded that the RANS model is not suitable for simulating suspended sediment transport, while the LES model gives more accurate results in agreement with the findings of experiments in a laboratory setting.

Large-eddy simulation is commonly used to study the flow over wavy lower boundaries (Broglia et al.,2003; Grigoriadis et al., 2012, 2013; Harris and Grilli,2012; Soldati and Marchioli, 2012) such as over sand ripples and sandbanks, which are very common in coastal waters and have a strong in fluence on sediment transport and wave energy dissipation.

Jackson and Hunt (1975) investigated the turbulent flow of air over hilly terrain in the laboratory, and showed that the flow acceleration reached a maximum above the hilltop, with the velocity over the top of the hill equal to the velocity at the upwind side of the hill at the same height. Henn and Sykes (1999) performed an LES of the flow over a rippled surface, and concluded that large spanwise fluctuations occur on the upslope boundary, and a detached shear layer exists in the lee of the crest, resulting in strong turbulence over the trough. Calhoun and Street (2001)also studied turbulent flow over a wavy surface, and showed that turbulence parameters above the valleys,such as turbulence intensity, turbulence transport and dissipation, and TKE, are larger than those above the hill tops. However, the mechanism behind the observed phenomena is still unclear.

Therefore, we use an LES model to explore the in fluence of a wavy seabed on the spatial distribution of the velocity and redistribution of the TKE budget,with Section 2 describing the LES model, and introducing the initial and boundary conditions, as well as important parameter settings. Section 3 presents an analysis of the simulation results, with a discussion and conclusions presented in Sections 4 and 5, respectively.

2 THEORY

2.1 Model description

The LES model used here is the PArallelized Large-Eddy Simulation Model (PALM) for Atmospheric and Oceanic Flows, which was developed by the Institute of Meteorology and Climatology of Leibniz University, Germany. In general, LES is based on the spatial average of the turbulent fluctuations, which are divided into largescale and small-scale eddies using a speci fied filter function. The large-scale eddies are then directly simulated, while the small-scale eddies are parameterized (Maronga et al., 2015). The basic governing equations for the LES model are given by

where t is time, and the coordinates ( i, j, k) represent the Cartesian ( x, y, z) directions, respectively. Two additional variables, the Earth’s rotational velocity Ω and the geographical latitude ?, describe the Coriolis parameter fi=(0, 2Ωcos( ?), 2Ωsin( ?)), ug,krepresents the geostrophic wind speed (as we only discuss the in fluence of topography without the geostrophic wind speed, ug,k=0), ρ is the density of seawater, p is the

In general, scalar quantities can be divided into resolved ψ0and unresolved scalars ψ*as

where ψ*＜ ＜ ψ0, with ψ represents absolute temperature T, pressure p and density ρ.

The equation of state

can be written as

where R=8.314 41± 0.000 26 [J/mol.K] is the molar gas constant, to give

Eq.1 can be written as

where δ is the Kronecker delta.

By applying a volume filter to the governing equations, those scales larger than the cutoff scale Δ x are retained and classi fied as ‘resolved scales’, while scales smaller than the cutoff scale Δ x are filtered, and classi fied as sub-grid scales. The filtered equation yields an additional term, which can be parameterized using an improved version of the 1.5-order closure from Moeng and Wyngaard (1988) and Saiki et al.(2000). Similar to the RANS model, PALM introduces eddy viscosities Kmand Kh, and assumes that the secondary moments of the mean quantities are proportional to their gradients by

where Kmis the subgrid-scale eddy viscosity and Khis the local subgrid-scale eddy diffusivity of heat. The viscosities can be parameterized as

Here, ρθis calculated by the method proposed by Jackett et al. (2006) using the polynomials determined by the variables Sa, θ and p (see Jackett et al., 2006,Table A2 ). Only the initial values of p are used here.

2.2 Initial conditions

The size of the domain is 10.0 m×10.0 m×10.0 m,with a mesh resolution of 0.1 m×0.1 m×0.1 m. The potential temperature at the sea surface is set to 275 K,the vertical gradient for -5 m ≤ z ≤ 0 is 1.5 K/100 m, and that for -10 m ≤ z ≤-5 m is 1.0 K/100 m. The salinity at the sea surface is 32.0, with vertical gradients for -5 m ≤z ≤ 0 and -10 m ≤ z ≤-5 m of 1.5/100 m and 1.0 /100 m,respectively. The u-component of the background velocity is set to 1 m/s, with the v-component set to zero. The u-component of the geostrophic velocity at the surface is 1 m/s, while the v-component at the surface is zero. The Coriolis force is negligible, and the roughness length is 0.1 m.

2.3 Boundary conditions

The upper boundary condition of the horizontal velocity components u and v follow the free-slip condition,

where u and v are the horizontal components of the water velocity at the surface. The lower boundary condition of u and v follow the no-slip condition,

The lower boundary condition of the perturbation pressure is set as p( k=0)=0, and the upper boundary condition of the perturbation pressure is set as

Fig.1 High hill topography ( H=5/2 m) at the seabed boundary

The lateral (streamwise and spanwise) boundary conditions are cyclic, with Dirichlet (in flow) and radiation (out flow) conditions allowed along the x- or y-direction.

If the topography is described by a linear function,a singularity results at the top of the hill. To avoid this, while simulating a more realistic ocean-bed topography, the seabed is described by trigonometric functions. Calhoun and Street (2001) modeled turbulent flow over a wavy surface and showed that a recirculation region develops with the amplitude of the trigonometric function. Here, we investigate two cases with different topographies δ=0.75 and δ=0.2,where δ=2 a/( λ/2) is the slope steepness, a is the amplitude of the wavy surface, and λ is the wavelength.The high hill topography is shown in Fig.1, and is

To discuss the TKE in different areas of the flow,and especially to compare the windward slope with the leeward slope of the hill (see Section 4.4), the hill topography is divided into five sections (s1, s2, s3, s4,s5) (Fig.1).

Fig.2 Low hill topography ( H=2/3 m) at the seabed boundary

The peak is located at (5 m, -7.5 m). The area is divided into five sections (s1–s5). Section s1 covers the area 1 m ≤ x ≤ 2 m, s2 covers the area 3 m ≤ x ≤ 4 m, s3 covers the area 4.5 m ≤ x ≤ 5.5 m, s4 covers the area 6 m ≤ x ≤ 7 m, and s5 covers the area 8 m ≤ x ≤ 9 m. The sections represent the first valley (s1), windward slope (s2), peak (s3), leeward slope (s), and second valley (s5).

According to the early work of Jackson and Hunt(1975), the area over a two-dimensional wavy lower boundary can be divided into two regions: an inner region and an outer region. Similarly, in our model,the area over a three-dimensional wavy seabed boundary is divided into an inner region and an outer region, with the two regions having a different dimensionless height H. Here, we estimate the inner region depth as hi= H (Figs.1 and 2). The TKE generated by the in fluence of the wave lower boundary cannot impact the outer region, while the lower hill produces the main turbulence shear in the inner region.

3 RESULT

According to the normalization method of Kantha and Clayson (2004), the u-, v-, and w-components of the velocity can be normalized by uw*corresponding to the water-side friction velocity, and calculated as

Here, ua*is the atmosphere-side friction velocity, ρa(bǔ)is the air density, and τ is the shear stress, which can be calculated using the wind speed 10 m above the sea surface ( u10) and the drag coefficient Cd,

Fig.3 Depth pro files of the high hill model showing the u-component of the velocity at sections s1–s5 in the y-direction after 1 200 s for (a) y=0.5 m, (b) y=3.5 m, (c) y=5.5 m, (d) y=7.5 m, and (e) y=10.0 m

where Cdis obtained from Wu (1980) as

Here we focus on the in fluence of the wavy seabed boundary. All velocity components have similar orders of magnitude because of the conservation of the whole water mass. Therefore, the water-side friction velocity uw*can normalize all velocity components.

The in fluence of the lower boundary on the turbulent statistics takes 1 200 s to reach a state of equilibrium.

3.1 u-component of the velocity

After initialization, the numerical model integrates for 7 200 s (2 h). In the y-direction, the pro files of the u-component of the velocity at s1–s5 are captured at y=0.5 m, 3.5 m, 5.5 m, 7.5 m, and 10 m, with Fig.3 showing these pro files after 1 200 s. The pro files at other times present almost identical values of the velocity, as well as velocity gradients, indicating that the model has attained a stable stage at this time.

Fig.4 Depth pro files of the high hill model showing the v-component of the velocity at sections s1–s5 in the y-direction after 1 200 s for (a) y=0, (b) y=3 m, (c) y=5 m, (d) y=7 m, and (e) y=10 m

In Fig.3, each pro file is divided into two parts. For-7 m ≤ z ≤ 0, has a maximum of about 1.5 m/s near z=-3.5 m, and decreases toward the upper surface, as well as at -7 m, which is the height around the peak of the hill. For z ≤-7 m, there is a countercurrent on the windward slope, with a maximum of about -0.5 m/s.On the leeward slope, the velocity is usually positive,with a maximum of about 0.5 m/s.

Furthermore, at z=-7 m, the vertical gradient of u is large, which implies that most of the TKE is generated in this area as a result of the shear induced by the topography. The area above the center of the trough is the maximum of the negative value, which implies high turbulence.

3.2 v-component of the velocity

Figure 4 shows the pro files of the v-component of the velocity at sections s1–s5 after 1 200 s, indicating similar pro files at all five sections at this time, and thus a lack of variation of v in the y-direction.

Fig.5 Depth pro files of the high hill model showing the w-component of velocity at sections s1–s5 in the y-direction after 1 200 s for (a) y=0.5 m, (b) y=3.5 m, (c) y=5.5 m, (d) y=7.5 m, and (e) y=10 m

At 1 200 s, there is a strong countercurrent on the windward slope of the middle hill (see Fig.4) with a maximum v of about -0.9 m/s, while v above the countercurrent is positive and has a maximum of about 0.9 m/s. At z=-7 m, TKE production is strong,especially over the windward slope where the two opposite groups touch. The countercurrent in the area above the trough is still strong, showing a similar distribution pattern to that of u. Thus, the TKE over the windward slope is still strong.

Figure 5 shows the pro files of the w-component of the velocity at sections s1–s5 after 1 200 s, which is negative above the windward slope with a maximum of about -0.9 m/s, implying the flow descends beyond the middle hill. The velocity above the leeward slope

3.3 w-component of the velocity

Fig.6 High hill model depth pro files of TKE budget parameters for 1 800 s, 3 600 s, 5 400 s and 7 200 s after model initialization

is 0–0.5 m/s, indicating ascending flow in the area above the lee side. Figure 5 indicates that the pro files of w in the y-direction are similar to each other, as well as to v in the y-direction, thus demonstrating the conservation of water mass, such that any volume of water flowing upward, must be replaced by water moving downward, resulting in convective motion.The pro files clearly show two counteracting groups at the hill top above the valley, indicating a strong TKE at this location, which agrees with the experimental results reported in Poggi et al. (2007).

3.4 TKE budget

3.4.1 Horizontal average over the total model domain Following Eq.1, the horizontally-averaged subgridscale TKE equation after tensor contraction is

where the quotation marks represent the subgrid-scale values and the overbars represent average values.Here, ρθcan be calculated from the state equation of seawater as proposed by Jackett et al. (2006), where ρθdepends on Sa, T and p (see Jackett et al., 2006,Table A2). For simpli fication, we set the value of ρθas constant to 1.025 7×103kg/m3. On the right-hand side of Eq.24, the first term represents the mean velocity transport Tm, the second term represents the shear production S, the third term represents the buoyancy production B, the fourth and the fifth terms represent the pressure and turbulent transports denoted P and T,respectively, the sixth term represents the molecular viscosity diffusion M, and the last term represents the turbulent dissipation D. The term on the left-hand side of the equation is denoted C, and represents the time derivative of the subgrid-scale TKE. Here, the buoyancy production B is zero because of the constant density ρθ, and the molecular diffusion M is several orders of magnitude smaller than the other terms on the right-hand side, so that it may be ignored. While Calhoun and Street (2001) showed that Tmis a signi ficant part of the TKE budget, our results show that the order of magnitude of Tmis less than 10-7m2/s3,while the minimum value of the other terms is 10-7m2/s3. Therefore, Tmhas a smaller effect in this simulation. Equation 24 implies that the rate of change over time of the subgrid-scale TKE C is equal to the sum of the mean velocity transport Tm, the production S, the transport terms P and T, and the turbulent dissipation D.

Fig.7 Low hill model depth pro files of the TKE budget parameters for 1 800 s, 3 600 s, 5 400 s, and 7 200 s

Fig.8 Depth pro file after 1 800 s comparing the shear production for different sections and in the total simulation domain

As shown in Fig.6, the subgrid-scale TKE remains almost constant in time, indicating that the rate of change of the subgrid-scale TKE with time is close to zero. Hence, the sum of the terms on the right-hand side of Eq.24 should also be close to zero. Indeed, in Fig.6, for each horizontal sea layer, the sum of S, Tm,P, T, and D is approximately zero.

Close to the sea surface and z=-7.5 m (the top of the hills), S and D change suddenly, while P and T do not. Changes in the area close to the sea surface result from various factors acting at the sea surface, such as the wind stress, so that the shear stress produces TKE.At the same time, the topographic changes in the area above the hilltops also cause a shear and TKE production. In the area above the hilltops, the maximum rate of change in TKE is at the peak of the hill. The redistribution of production and dissipation results in a variation of P and T, resulting in changes to the transfer of energy.

3.4.2 Shear production

As shown in Fig.8, the simulation domain is divided into five sections (s1, s2, s3, s4, s5), with s1 covering the area 1 m ≤ x ≤ 2 m, s2 the area 3 m ≤ x ≤ 4 m, s3 the area 4.5 m ≤ x ≤ 5.5 m, s4 the area 6 m ≤ x ≤ 7 m, and s5 the area 8 m ≤ x ≤ 9 m. The average is the horizontal average over each section, including the total domain average (s0).All the sections show a maximum at z=-7.5 m, with s3 having the highest value, indicating that the shear production is largest above the hilltop. Moreover, the shear production above the leeward slope is larger than that above the windward slope, indicating stronger turbulence above the leeward slope.

Fig.9 Depth pro file after 1 800 s

3.4.3 Dissipation

As shown in Fig.9, all the sections show peaks at z=-7.5 m, with s3 having the maximum value,indicating that the dissipation is largest above the hilltop. Thus, the turbulence dissipates rapidly in accordance with the slope of the topography.Moreover, the dissipation above the windward slope is smaller than that above the leeward slope, indicating reduced turbulence dissipation above the windward slope.

4 DISCUSSION

Analysis of the turbulent kinetic energy produced by the LES model PALM allowed the in fluence of a wavy lower boundary on the redistribution of the TKE budget to be investigated. One of the first studies to investigate in detail the turbulent air flow over a low hill in the atmospheric boundary layer, Jackson and Hunt (1975) found that the maximum velocity occurs mostly in the area above the top of the hill, with the velocity over the upwind slope of the hill at the same elevation almost equal to that over the top of the hill.Although Jackson and Hunt (1975) used a low hill as the topography and some of their conclusions are similar to those of our study, their conclusions only apply to the atmospheric boundary layer. In 1994, the Coastal and Hydraulics Laboratory (formerly the Coastal Engineering Research Center) of the US Army Engineer Waterways Experiment Station performed a variety of laboratory experiments including those with different topographies (http://www.frf.usace.army.mil/duck94/DUCK94.stm). Their experiment,‘Sediment Dynamics in the Nearshore Environment’,showed that (1) a wavy lower boundary can lead to higher water velocities, and (2) the velocity decays as the elevation rises. These conclusions are similar to our findings. Poggi et al. (2007) also investigated the distribution of the velocity field over a hilly surface,and simpli fied the filtered equations. Their results regarding the velocity field are similar to ours, although our analysis includes the in fluence of the wavy seabed boundary on the TKE budget redistribution.

We considered two different wavy boundaries with different relative hill height based on the steepness δ to distinguish the different hill heights ( δ=0.75 and 0.2). As the hill height decreases, the resolution of the calculation needs to increase and, therefore, the requirements of the computational performance and time are higher, leading to higher costs. We found that, for a sufficiently low hill height, the in fluence of the wavy lower boundary on the TKE budget redistribution only slightly varies. It is worth noting that, apart from the shape and height of the hills, the tilt of the hills and the distance between two adjacent hills also affects the TKE budget redistribution.

5 CONCLUSION

An LES model to simulate the in fluence of a wavy lower boundary allowed the analysis of the distribution of the velocity components and the redistribution of the TKE budget. The analysis of the distribution of the u-, v-, and w-components of the velocity shows that the flow over the center of the trough exhibits large velocity gradients, implying that the area over the center of the trough has a strong shear, leading to high TKE production.

The analysis of the redistribution of the TKE budget shows that, apart from the area close to the water surface where shear is produced by the complex sea-surface motion, the area above the seabed hilltops also produces turbulence. Our results show that shear production dominates the turbulent kinetic energy,although we have neglected buoyancy production here. The effect of the mean velocity transport is small. Furthermore, the shear production reaches a maximum above the hilltop and is larger above the leeward slope than above the windward slope.

Finally, the existence of a recirculation region(Calhoun and Street, 2001) over the trough is con firmed for different heights of the hill topography.Most of the TKE distribution occurs in the layer above the trough near the top of the hill, followed by in the recirculation region near the trough.

6 ACKNOWLEDGEMENT

We thank the two anonymous reviewers for their constructive comments.