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 profile, 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 specific 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 influence 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 influence 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 large-scale and small-scale eddies using a specified 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

and

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 f_i=(0, 2Ωcos(ϕ), 2Ωsin(ϕ)), u_{g, k} represents the geostrophic wind speed (as we only discuss the influence of topography without the geostrophic wind speed, u_{g, k}=0), ρ is the density of seawater, p is the hydrostatic pressure, is the subgrid-scale TKE, g is the acceleration due to gravity, μ/ρ=υ is the kinematic viscosity, μ is the dynamic viscosity, ε is the Levi-Civita symbol, T is the absolute temperature, Sa is the practical salinity, and Ψ_T and Ψ_Sa are the source/sink of T and Sa, respectively. The quotation marks indicate the subgrid-scale values, and the overbar represents the average value.

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

where ψ^*≪ψ₀, 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

Since

and

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 classified as 'resolved scales', while scales smaller than the cutoff scale Δx are filtered, and classified 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 K_m and K_h, and assumes that the secondary moments of the mean quantities are proportional to their gradients by

and

where K_m is the subgrid-scale eddy viscosity and K_h is the local subgrid-scale eddy diffusivity of heat. The viscosities can be parameterized as

and

where c_m=0.1, Δ=, where Δx, Δy and Δz are the grid intervals in the x, y and z directions, respectively, l is the subgrid-scale mixing length, which is determined by the height z (the distance from the wall when the topography is used), Δ, and the stratification, where

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 .

The lateral (streamwise and spanwise) boundary conditions are cyclic, with Dirichlet (inflow) and radiation (outflow) 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 δ=2a/(λ/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 described by the function , where is the height of the hill. The low hill topography is shown in Fig. 2, and is described by the same function, where . Calhoun and Street (2001) showed that, when the amplitude of the wavy boundary surface is high, a recirculation occurs in the trough, which decreases with the amplitude a until disappearing. Accordingly, we simulate the flow over two different topographies to further investigate the recirculation region.

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).

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 h_i=H (Figs. 1 and 2). The TKE generated by the influence 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 u_w^* corresponding to the water-side friction velocity, and calculated as

Here, u_a^* is the atmosphere-side friction velocity, ρ_a is the air density, and τ is the shear stress, which can be calculated using the wind speed 10 m above the sea surface (u₁₀) and the drag coefficient C_d,

where C_d is obtained from Wu (1980) as

Here we focus on the influence 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 u_w^* can normalize all velocity components.

The influence 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 profiles 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 profiles after 1 200 s. The profiles 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.

In Fig. 3, each profile 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 profiles of the v-component of the velocity at sections s1–s5 after 1 200 s, indicating similar profiles at all five sections at this time, and thus a lack of variation of v in the y-direction.

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.

3.3 w-component of the velocity

Figure 5 shows the profiles 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 is 0–0.5 m/s, indicating ascending flow in the area above the lee side. Figure 5 indicates that the profiles 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 profiles 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 subgrid-scale 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 simplification, we set the value of ρ_θ as constant to 1.025 7×10³ kg/m³. On the right-hand side of Eq.24, the first term represents the mean velocity transport T_m, 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 T_m is a significant part of the TKE budget, our results show that the order of magnitude of T_m is less than 10^-7 m²/s³, while the minimum value of the other terms is 10^-7 m²/s³. Therefore, T_m has 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 T_m, the production S, the transport terms P and T, and the turbulent dissipation D.

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, T_m, 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.

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 influence 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 simplified the filtered equations. Their results regarding the velocity field are similar to ours, although our analysis includes the influence 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 influence 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 influence 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 confirmed 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.