Effects of a two-equation turbulence model on the simulation of the internal lee waves

Xuhui Hu , Xuefeng Zhng , , Lei Lin , Linxin Zhng , Shifn Wng

a School of Marine Science and Technology, Tianjin University, Tianjin, China

b College of Ocean Science and Engineering, Shandong University of Science and Technology, Qingdao, China

c Key Laboratory of Marine Environmental Information Technology, National Marine Data and Information Service, Ministry of Natural Resources, Tianjin, China

Keywords:Non-hydrostatic MERF Internal wave Turbulent mixing

A B S T R A C T In order to understand the wave—turbulence interaction under non-hydrostatic conditions to prepare future advanced very-high-resolution ocean reanalysis data, an σ-coordinate ocean model —namely, the Marine Environment Research and Forecasting (MERF) model —with an idealized supercritical slope topography is applied to conduct a series of high-resolution numerical experiments with and without the non-hydrostatic approximation.The popular Mellor—Yamada two-equation turbulence model (MY2.5) is enclosed in MERF to validate its effect on small-scale internal lee waves. Instantaneous results show that the internal lee-wave processes are relaxed through employment of the MY2.5 scheme, whether or not in the non-hydrostatic model. Time averaged results suggest the influences of the vertical mixing parameterization scheme on the numerical results are more dominant than the non-hydrostatic/hydrostatic selection for the large-scale dynamic process. Besides, diagnostic analyses of the energy budget show that the spread of internal lee waves at the slope is dramatically suppressed by the vertical turbulent mixing, indicating more tidal energy is able to be converted into the irreversible mixing when the two-equation turbulence model is employed.

1. Introduction

Internal waves are a kind of dynamic process existing universally in the ocean. The energy of large and mesoscale processes can be transferred to small-scale ones through the nonlinear wave—wave interaction associated with internal waves, and eventually dissipates in the form of turbulent mixing ( Orr and Mignerey, 2003). The safety of marine engineering and underwater moving objects is also dramatically affected by the vertical fluctuation and horizontal shear produced by internal waves in the ocean ( Osborne and Burch, 1980 ; Du and Fang, 2003 ).With the continuous improvement in computing power and more and more urgent needs for small-scale signals in marine environmental security, non-hydrostatic ocean models have been developed over the last decade ( Liu et al., 2016 ; Wu et al., 2018 ), and thus numerical simulation has become a common method to study internal waves ( Lai et al., 2010 ;Zhang et al., 2011 ).

To gain insight into internal waves, Xing and Davies (2006) conducted a series of numerical simulation experiments on tidally induced internal lee waves using a non-hydrostatic model in cross-sectional form with idealized topography representing a sill. Berntsen et al. (2008) used the BOM (Bergen Ocean Model) to investigate the generation and propagation of tidally forced internal waves near the sill, and further analyzed the sensitivity of numerical results to horizontal grid size. Xing and Davies (2009a) found that the existence of the second sill, which is close to the first one, led to standing internal wave generation with associated changes both in wave spectrum and Richardson number in the inter sill region. Subsequently, Xing and Davies (2009b) studied the influence of bottom friction upon unsteady lee-wave generation. Both short- and longer-term integration with a number of buoyancy frequencies confirmed that increasing bottom friction modified the flow and unsteady lee-wave distribution on the downstream side of a sill. Further, the role of the bottom boundary layer upon lee-wave generation in sill regions was investigated in deeper water and shallower water,separately ( Xing and Davies, 2010 ). For a given tidal forcing, as the sill depth increased, the lee-wave response and vertical mixing decreased( Xing and Davies, 2011 ).

Although the overall performance of numerical methods is relatively high, it is undeniable that there are still some differences between simulation and real-world results. In response to ocean turbulent mixing, a fundamental problem in ocean modeling, parameterization schemes of vertical mixing have played key roles in determining model performance( Ma et al., 2014 ; Zhu and Zhang, 2018 ). Li et al. (2016) compared the simulation capabilities of three vertical mixing schemes —namely, KPP(K-Profile Parameterization), KT (Kraus and Turner), and MY2.5 (Mellor and Yamada) —in tropical and mid—high latitude regions. At the climate scale, the overall simulation effect of the KPP scheme was the best in the tropical sea area. In middle and high latitudes, there was little difference between the KPP scheme and MY2.5 scheme, while the KT scheme was more suitable over those regions. The performance of vertical mixing parameterization schemes is not only affected by geographical locations,but also changes with different physical phenomena. Although the average flow fields in the upper ocean simulated by the above three mixing schemes were similar to OFES (OGCM for the Earth Simulator) data,slight mismatches were found in flow direction and amplitude ( Li et al.,2017 ).

The turbulence closure scheme based on the Reynolds averaging method has still been a frequently used method incorporated into a non-hydrostatic model when the small-scale dynamical processes are resolved, since LES (large eddy simulation) is too insensitive for a fine-scale simulation with a fine resolution

O

(1)m . Despite the earlier studies on non-hydrostatic models and vertical turbulent mixing mentioned above, the effects of turbulence closure schemes on simulation of small-scale internal waves have not yet been thoroughly investigated with non-hydrostatic models. Since the Mellor—Yamada closure scheme( Mellor and Yamada, 1982 ) is one of the most popular turbulent closure models enclosed in many ocean models, including non-hydrostatic ones, this study uses the Marine Environment Research and Forecasting(MERF) model ( Liu et al., 2016 ; Lin and Liu, 2019a , 2019b ), a nonhydrostatic model, to simulate tidal-induced internal waves to explore the effects of the MY2.5 scheme on the numerical results, in order to provide a reference for further small-scale reanalysis research with nonhydrostatic models.

This paper is organized as follows. The setup of the numerical experiments and calculation of diagnostic variables are given in Section 2 .Section 3 reports the sensitivity of the results to vertical mixing parameterization schemes. A summary and general discussion are provided in Section 4.

2. Numerical experiments

2.1. Model and experimental setup

The ocean model used in this study is MERF, which is introduced in detail in Liu et al. (2016) . Cartesian coordinates (

σ

-coordinates) are used in the horizontal (vertical) direction. (The

σ

-coordinates may produce pressure gradient errors in regions of steep topography and coarse resolution, while the experiments here are conducted with sufficiently small grid size, which can effectively depress the negative impact). The model variables are discretized on a C-grid. Vertically, the conversion from

z

-coordinates to

σ

-coordinate follows the rule

σ

= (

z

?

η

)∕(

H

+

η

) , where

η

is the surface elevation and

H

is the bottom depth. The total pressure is decomposed into two parts: hydrostatic and dynamic. The PIFD (par-tially implicit finite difference) scheme proposed by Liu et al. (2016) is used to calculate the dynamic pressure. The solution to the surface elevation is mainly based on

θ

-method implicit calculation ( Casulli and Cheng, 1992 ). MERF has been tested by five benchmark cases ( Liu et al.,2016 ; Lin and Liu, 2019a ).This study relies on a tidally induced internal lee wave experiment.The initial conditions are set according to Xing and Davies (2006) and Berntsen et al. (2009) . The depth at the western boundary is constant —namely,

h

= 50 m —and the sill depth

h

= 15 m. On the eastern side,

h

= 100 m. The initial temperature decreases linearly from the sea surface to the sea floor. The salinity is spatially uniform and constant. The initial buoyancy frequency is

N

= 0.01

s

in all experiments.Horizontally, the applied grid sizes obey a 10 m equidistant distribution. Vertically, there are 40 uniform layers. An

M

tide (principal lunar semidiurnal tidal constituent) is imposed at the western open boundary of the domain with a velocity of

u

= 0.3 sin(1.40526 ×10

t

) m s.The time step is set to 0.2 s, and the model is integrated for 2 days (four tidal periods).The horizontal viscosity (diffusivity)

A

M

= 10ms(

A

=10ms) is in all studies, which is consistent with Liu et al. (2016) .The detailed settings of the vertical mixing coefficient are shown in Table 1 .

Table 1 List of experiments.

2.2. Model diagnostics

Kang and Fringer (2012) proposed a theoretical framework that is now widely used in diagnostic calculation ( Chen et al., 2013 ; Han and Eden, 2019 ) to analyze internal tide energetics, including the nonlinear and non-hydrostatic effects, based on which the energy analysis in this study is performed. The velocity is divided into barotropic and baroclinic parts as follows:

u

=

U

+

u

and

w

=

W

+

w

. The barotropic velocities are defined as

where

η

is the surface elevation and

H

is the bottom depth. The total density is given by

The baroclinic kinetic energy

E

is computed in each time step according to

The APE

E

p is computed in each time step according to

3. Results

3.1. Results at special moments

3.1.1.

Non-hydrostatic

solution

In all experiments ( Table 1 ), the sill depth

h

is initially set to 15 m,and the buoyancy frequency

N

is 0.01 swith a tidal current forcing of

U

= 0.3 m s. This means that the maximum Froude number

F

=

U

∕

N

h

>

1 , indicating the characteristic of supercritical flow. At

t

= 1/8

T

(where

T

is the

M

tidal period), there is a steep vertical displacement of temperature contours (

x

= 500 m) on the lee side of the sill ( Fig. 1 (a)). This phenomenon can be explained by the formation of a hydraulic jump related to the change between the top of sill and the lee side of the sill in Froude number ( Xing and Davies, 2007 ). The speed distribution at this moment ( Fig. 1 (a)) can also illustrate the change well.The horizontal velocity

u

at the top of sill is the maximum, which is about 0.7 m s. However, the across-sill velocity on the right-hand side of the sill is less than 0.3 m s. Corresponding to the downward movement of the isotherm, the vertical velocity

w

in this area is also changed( Fig. 1 (a)), with the upwelling and downwelling alternating. At the same time, the positive value of dynamic pressure

QP

is next to the negative value on the lee side of the sill ( Fig. 1 (a)), meaning the mixing will occur.

Fig. 1. (a) Distribution of the temperature field (units: °C; contour interval: 0.25°C), horizontal velocity u (units: m s ? 1 ; contour interval: 0.05 m s ? 1 ), vertical velocity w (units: m s ? 1 ; contour interval: 0.02 m s ? 1 ), and dynamic pressure QP (units: N; contour interval: 0.001 N) at t = 1/8 T , where T is tidal period, computed with the non-hydrostatic model and constant vertical mixing coefficients (CONH). (b) As (a) but at t = 1/4 T . (c) As in (a) but at t = T . (d) As in (a) but just for temperature(°C), horizontal velocity u and vertical velocity w , computed with the hydrostatic model and constant vertical mixing coefficients at t = 1/4 T (COH).

At the time of maximum inflow (

t

= 1/4

T

), internal lee waves are generated on the lee side of the sill. The contours of the temperature field ( Fig. 1 (b)) from near the surface to the middle (

z

= ? 50 m) show the obvious vortex shape, while the temperature contours at depth are a bit like sine curves. The distribution of the horizontal velocity

u

near surface ( Fig. 1 (b)) corresponds to the isotherms’ status with the wavelength in the order of 100 m, which is the same as lee waves. The signs of internal waves can still be seen from the contours of dynamic pressure and vertical velocity on the lee side of the sill. Overall, the intensity of each physical quantity discussed above at 1/4

T

is much stronger than that at 1/8

T

.After a tidal cycle, the lee waves are almost invisible. The velocity components (

u

,

w

) and the dynamic pressure

QP

in the study area are close to zero ( Fig. 1 (c)). The distribution of temperature is not the same as the initial state. Large areas of well-mixed water appear in the upper layer on the right-hand side of the sill, as they do in the bottom layer ( Fig. 1 (c)). The increase in well-mixed areas is related to a series of physical processes that promote the density overturning to strengthen the shear-driven and buoyancy-driven turbulent mixing,such as the generation of the hydraulic jump, lee waves, and shear instability.When the constant vertical mixing coefficients (vertical viscosity

K

and vertical diffusivity

K

) are replaced by the variable coefficients from the MY2.5 scheme, something different happens ( Fig. 2 ). Although the hydraulic jump phenomenon also occurs at

t

= 1/8

T

in this case( Fig. 2 (a)), the change trend of the isotherm at

x

= 500 m is a bit smoother and the temperature in most areas is 0～1.3 °C lower than the case of constant coefficients ( Figs. 1 (a) and 2(a)). The most obvious difference is at the time of maximum inflow (

t

= 1/4

T

). No sign of lee-wave generation is seen at this special moment, except there are some faint ripples at the depth of the temperature field. The wave processes are smoothed by the viscosity and diffusivity terms ( Berntsen et al., 2008 ).Moreover, a surface temperature front advected by the tidal current occurs at about

x

= 1700 m rather than the spiral isotherm near the surface shown in Fig. 1 (b). Different from the case of constant coefficients, the temperature field distribution on the left of the sill at

T

( Fig. 2 (c)) indicates that the area has been well mixed. The MY2.5 scheme removes the physically realistic density inversion that can occur in a non-hydrostatic model. In general, the intensity of each variable in the case computed with the non-hydrostatic model and MY2.5 scheme (MYNH) at these three special moments ( Fig. 2 ) is weaker than the counterpart in the case of constant coefficients ( Fig. 1 ).

3.1.2.

Hydrostatic

solution

As described in Section 2 , the settings for the case computed with the hydrostatic model and constant vertical mixing coefficients (COH)and the case computed with the hydrostatic model and MY2.5 scheme(MYH) are the same as before, except that the hydrostatic model is used for the simulation. In the hydrostatic models, the steepening of the internal waves and eventual overturning of the density surface is prevented by the generation of rapid artificial convective mixing ( Xing and Davies, 2007 ). Large areas of well-mixed signals are found in the upper and middle regions (top 60 m) on the lee side of the sill in the hydrostatic model test with the constant vertical mixing coefficients( Fig. 1 (d)). In the deep area of the temperature field, sine ripples also exist, but the amplitude is not as large as those in the non-hydrostatic model case ( Fig. 1 (b,d)). The distribution of the

u

velocity field on the right of the sill ( Fig. 1 (d)) is smaller than that in the non-hydrostatic test( Fig. 1 (b)). From the vertical velocity field distribution ( Fig. 1 (d)), it can be seen that the wavelength generated from the lee waves in this case is longer, while the distance traveled downstream is farther compared with the result of the non-hydrostatic model with the same configurations ( Fig. 1 (b)).The differences between the hydrostatic and non-hydrostatic cases discussed above also occur in the tests associated with the MY2.5 scheme. The surface temperature front in the hydrostatic model( Fig. 2 (d)) is shallower than the one in the non-hydrostatic model( Fig. 2 (b)). As in the cases with constant mixing coefficients, the amplitude of the deep isotherm is smaller than that using the non-hydrostatic model ( Fig. 2 (b,d)). The horizontal velocity

u

( Fig. 2 (d)) is close to the counterpart found in the non-hydrostatic model ( Fig. 2 (b)). However,in the tests using the MY2.5 scheme, the vertical velocity from the hydrostatic model is generally larger than the one in the non-hydrostatic model, indicating the vertical mixing intensity in the hydrostatic model is stronger.

According to the above analysis, it is found that the small-scale dynamic process can be obviously suppressed by the MY2.5 scheme in both the hydrostatic and the non-hydrostatic model. In order to eliminate the instantaneous effects, a comparison of the results from the time average is performed.

3.2. Results of the time average

The relative spatial distributions of these physical quantities averaged to

t

=

T

are demonstrated in Fig. 3 . As shown in Fig. 3 (a), the average temperature field in the first tidal period in the case computed with the non-hydrostatic model and constant vertical mixing coefficients(CONH) is uniformly distributed. The small-scale signals at the instants disappear. Compared with the case of constant vertical mixing coeffi-cients, the simulated temperature in the MY2.5 case exhibits obvious differences. The results in MYNH are generally cold, especially in the middle layer on the left side of the sill and the upper-middle layer on the right side from about

z

= ? 10 m to

z

= ? 50 m. In addition to cold deviations, there is also a warm bias in the surface and the top of sill ( Fig. 3 (b)). In the difference between COH and CONH ( Fig. 3 (c)),the bias in the right-hand central part of the sill is negligible, which is

O

( 102 )C . The warm bias, however, remains in the surface. The differences between MYH and CONH ( Fig. 3 (d)) are similar to those between MYNH and CONH ( Fig. 3 (b)). Only some slight cold biases arise in the middle and bottom areas on the right side of the sill (not shown). This means that the influence of vertical mixing parameterization schemes is greater than that of model selection.Among the four physical quantities (horizontal velocity

u

and vertical velocity

w

not shown) averaged to

t

=

T

, the case of constant vertical mixing coefficients and the MY2.5 case have the most obvious difference in dynamic pressure simulation ( Fig. 3 (e,f)). It is worth noting that the dynamic pressure biases between MYNH and CONH and the reference dynamic pressure (CONH) are of a similar spatial pattern but with opposite signs. In addition, the absolute value of the corresponding position is almost the same, further revealing that the dynamic pressure field signal obtained by the MY2.5 scheme is weaker and close to zero.

Fig. 2. (a) Distribution of the temperature field (units: °C; contour interval: 0.25°C), horizontal velocity u (units: m s ? 1 ; contour interval: 0.05 m s ? 1 ), vertical velocity w (units: m s ? 1 ; contour interval: 0.02 m s ? 1 ), and dynamic pressure QP (units: N; contour interval: 0.001 N) at t = 1/8 T , where T is tidal period, computed with the non-hydrostatic model and MY2.5 (MYNH). (b) As in (a) but at t = 1/4 T. (c) As in (a) but at t = T . (d) As in (a) but just for temperature (°C), horizontal velocity u and vertical velocity w , but computed with the hydrostatic model and MY2.5 (MYH).

3.3. Analysis

To complement the previous comparisons for the constant vertical mixing coefficients and the MY2.5 scheme, the baroclinic kinetic energy and APE in the area near the sill are shown in Fig. 4 . As a whole,the baroclinic kinetic energy and the APE both peak at the time of maximum inflow (

t

= 1/4

T

), and then decrease when the flow reverses at its maximum (

t

= 3/4

T

) in those four different tests. Among the four special moments, the energy is the least at

t

= 1/8

T

. This trend shows that the generation of internal waves affects the distribution of energy in the ocean. As can be seen from Fig. 4 (a,b), in both the non-hydrostatic and hydrostatic model, baroclinic kinetic energy obtained from the MY2.5 scheme is smaller than that from the constant vertical mixing coeffi-cients scheme at

t

= 3/4

T

, the difference of which is about 7 ×10J mand 9 ×10J m, respectively. Compared with the baroclinic kinetic energy, the APE in Fig. 4 has the same magnitude, indicating that APE is also important in energy analysis. It refers to the small active portion of potential energy that is available to be converted into kinetic energy, which can ultimately contribute to mixing ( Chen et al., 2013 ).The distribution of APE at different times from CONH and MYNH shown in Fig. 4 (c) implies that more tidal energy is transferred to the mixing in the MY2.5 scheme with the non-hydrostatic model case compared with the constant mixing coefficient test ( Berntsen et al., 2008 ).

Fig. 3. (a) Distribution of the temperature field (units: °C; contour interval: 0.25°C) computed with the non-hydrostatic model and constant vertical mixing coefficients(CONH) averaged to t = T . (b) Temperature differences (units: °C; contour interval: 0.1°C) between MYNH and CONH (MYNH minus CONH), averaged to t = T . (c)As in (b) but between COH and CONH (COH minus CONH. (d) As in (b) but between the MYH and CONH (MYH minus CONH). (e) As in (a) but for dynamic pressure QP (units: N; contour interval: 0.0001 N). (f) As in (b) but for dynamic pressure QP (units: N; contour interval: 0.0001 N).

Fig. 4. (a) Baroclinic kinetic energy of four special moments with the non-hydrostatic model. (b) As in (a) but with the hydrostatic model. (c) Available potential energy of four special moments with the non-hydrostatic model. (d) As in (c) but with the hydrostatic model. The results for the constant vertical mixing coefficients experiments are given with solid lines and circles. The results for the MY2.5 scheme are given with dashed lines and circles.

4. Conclusion

Based on the MERF model with the non-hydrostatic scheme, the interaction of the tidal-induced internal waves and vertical turbulent mixing are numerically explored. Cases with the hydrostatic condition and constant mixing coefficients serve as a benchmark to analyze the scenario for the tidally induced internal lee waves. The conclusions are as follows.

Firstly, characteristics of the internal lee waves can be well embodied using the given constant mixing coefficients, consistent with previous research ( Xing and Davies, 2007 ). However, the small-scale dynamic process is easily impaired by the MY2.5 scheme, either in the hydrostatic case or in the non-hydrostatic one. The propagation of internal lee waves can be obviously decayed by the traditional MY2.5 scheme.

Secondly, numerical cases with and without the non-hydrostatic scheme, as well as with and without the turbulent parameterization scheme, are discussed from the perspective of the time average. With the MY2.5 scheme, the dynamic pressure associated with the nonhydrostatic process tends towards zero. The difference in the cases among the temperature illustrate that when the large-scale, timeaveraged processes are concerned, the effects of the vertical turbulent mixing schemes on the dynamic and thermodynamic processes are more prominent than the hydrostatic/non-hydrostatic alternative.

Finally, energy budget results confirm the view that the simulation with the MY2.5 scheme results in more tidal energy being converted into irreversible mixing. Therefore, the turbulent mixing from the MY2.5 scheme is more intensive than that from the constant vertical mixing coefficients scheme. The former restricts the baroclinic kinetic energy and enhances the APE. This will provide some insights into prospective small-scale ocean reanalyses where the cascade of energy happens from the grid-resolved scale to the sub-grid scale.

The horizontal grid size is set to 10 m in this study, which is a very fine level and thus does not cause serious pressure gradient errors. In practice, it is still inevitable that

σ

-coordinates will be used in regions of steep topography and coarse resolution. In order to improve the simulation accuracy of non-hydrostatic ocean models, the version of MERF based on generalized vertical coordinates is in the development stage,which will further reduce the error.

Declaration of Competing Interest

No potential conflict of interest was reported by the authors.

Funding

This study was supported by the National Key Research and Development Program of China [grant number 2016YFC1401800 ], the National Programme on Global Change and Air—Sea Interaction of China [grant number GASI-IPOVAI-04 ], and the National Natural Science Foundation of China [grant numbers 41876014 and 41606039 ].

Acknowledgments

MERF has been upgraded to MERF-NH 2.2, primarily because of the adoption of the Mellor—Yamada turbulence closure scheme (MY2.5).