## Abstract

Gravity currents are important in many fields, including the estuarine sciences, meteorology and hydraulic engineering. The NHWAVE (non-hydrostatic wave) model was applied to simulate the detailed interface structure between a lock-release gravity current and the ambient fluid. The simulated structures, including the front height, front position and velocity of the current, are consistent with the results of laboratory experiments. However, the internal structure of the current is different from that revealed by previous research. The Kelvin–Helmholtz phenomenon in the interface and the interface vortices were successfully captured by the NHWAVE model. The difference in velocity between the front and rear vortices leads to entrainment, further causing changes in the shapes and amount of vortices. Flow field results obtained by the NHWAVE model reveal the existence of a significant circular flow, as well as some small eddies within it. The significant circular flow supports the forward movement of the current, whereas the small eddies reflect interface vortices. In contrast, hydrostatic simulation with the same model settings fails to capture the vortices. This research shows that the NHWAVE model performs better than a hydrostatic model when simulating the Kelvin–Helmholtz instability phenomenon and vortex entrainment in a lock-release gravity current.

## INTRODUCTION

A gravity current is basically a horizontal flow driven by a difference in density between two fluids. Many natural phenomena, including the salt wedge in an estuary, a snow slide, a sand storm and density currents in a reservoir are closely related to gravity currents. The bottom outlet of a reservoir is dependent on the front location of a gravity current with sand (Hu *et al.* 2012) and when mixing caused by salinity in an estuary occurs the Kelvin–Helmholtz instability can be observed (Pu *et al.* 2015). Gravity currents are important in many fields which are important to the economy of China, including the estuarine sciences, meteorology and hydraulic engineering (Fan 2011).

Methods of studying gravity currents include theoretical derivation, laboratory measurement and numerical simulation. In theoretical studies, researchers mostly investigate gravity currents from a hydro-mechanical perspective and propose that a current undergoes three phases, namely, a slumping phase, an inertial phase and a viscous phase. For example, Huppert & Simpson (1980) investigated by dynamic analysis the relationship between the propagation distance and duration of each phase. A primary means of laboratory observation of gravity currents is flume experiments, which are usually performed in combination with particle image velocimetry (PIV). Hacker *et al.* (1996) used this experimental technique to determine the density structure of a gravity current and identified the importance of mixing for the movement of the current. They also postulated that by breaking Kelvin–Helmholtz waves a stratified fluid was produced in the rear of the current. Peng *et al.* (2015) recorded the flow structures of a gravity current using high-speed PIV, and found that the streamwise velocity profiles would oscillate when the slumping phase developed into the inertial phase. He *et al.* (2016) employed PIV to study the velocity field and vortex field of a gravity current. Their results showed that the Kelvin–Helmholtz instability phenomenon is one of the factors causing vortices. Since a gravity current involves a nonlinear, complex motion process, numerical simulation is often necessary when studying it in detail. On the basis of experimental measurements, Chen & Lee (2001) used a two-equation turbulence model to simulate the flow and mixing of a current. They found that the renormalization group (RNG) model for Reynolds-stress closure simulated the transitional and local turbulence effectively. Zhang *et al.* (2004) studied the head structure of a lock-release gravity current by using the RNG model and generated a velocity vector diagram, the distribution of turbulent kinetic energy and concentration contours.

However, turbulence in the interface is a complex element of hydromechanics, the simulation of internal mixing or entrainment is hard to achieve and the Kelvin–Helmholtz phenomenon is difficult to capture minutely. There are several major simulation methods that can be used to study the turbulence such as direct numerical simulation and large-eddy simulation (Ishihara *et al.* 2009; Remmler & Hickel 2013; Wei *et al.* 2015). Notably, due to the development of methods of dividing pressure into dynamic pressure and hydrostatic pressure, non-hydrostatic models are able to simulate turbulent phenomena well. In the Kelvin–Helmholtz phenomenon of a gravity current there is an intensive mixing between two fluids with different densities (Samothrakis & Cotel 2006). Intensive mixing causes a violent vertical motion of fluids, which is regarded as a high dynamic pressure in this process. In turn, the dynamic pressure has an effect on the vertical velocity. In such numerical studies, the pressure calculation term should be divided into static pressure and dynamic pressure with the dynamic pressure receiving special attention.

At present, the frequently-used non-hydrostatic models include NHWAVE, SWASH, SUNTANS, etc. Some researchers have introduced dynamic pressure into a hydrostatic model, and the hydrostatic models are thus changed into non-hydrostatic models, such as FVCOM-NH, ROMS, etc. These non-hydrostatic models have been widely used for modelling of estuarine and coastal hydrodynamics. For example, a model-splitting method has been developed to simulate breaking waves (Zhang 2017). Likewise, the simulation of a lock-release gravity current also requires the application of a non-hydrostatic model to progress the research of such currents.

The NHWAVE (non-hydrostatic wave) model has been successfully applied to the study of many non-hydrostatic processes, such as the bubble flow in a surf zone, the dispersion of a surface wave, a tsunami wave generated by deformable submarine landslide (Ma *et al.* 2011, 2012, 2013), etc. Consequently, in this paper, the NHWAVE model has been applied to simulate a lock-release gravity current.

## METHODS

*t*are given by: where ,

*h*is water depth, is surface elevation, is the velocity in the direction of Cartesian coordinates , is the vertical velocity in the -coordinate, is the water density calculated by temperature and salinity in this numerical experiment,

*p*is the dynamic pressure, the hydrostatic pressure is , and are turbulent diffusion terms.

Developers have incorporated a nonlinear model into NHWAVE to estimate the turbulent kinematic viscosity. A combined finite-volume and finite-difference scheme with a Godunov-type method was applied in the model and the computational variables were re-planned. Velocity was defined at the cell centre and dynamic pressure *p* was defined at the vertical cell face. To obtain second-order temporal accuracy, a two-stage second-order nonlinear strong stability preserving (SSP) Runge–Kutta scheme was adopted for time-stepping. The time-step was adaptive during the simulation, according to the Courant–Friedrichs–Lewy (CFL) criterion. The Courant number was taken to be 0.5 to ensure accuracy and stability in the current model.

*p*solved, the non-hydrostatic velocities at each stage can be updated in Equation (4).

## RESULTS AND DISCUSSION

### Model setup

In this paper, the numerical experiment to analyze the propagation of a gravity current is the same as the seventh laboratory experiment described in Chen & Lee (2001), on the basis that instantaneous frames of this case showed, in some detail, the start of the collapse of the dense fluid when the lock was released and the current shapes which evolved (unfortunately the quality of these frames was unsuitable for reproduction in this paper). In the numerical experiment, the horizontal length is 1.316 metres and the depth is 0.1 metres. There are 526 grids in the horizontal () direction and 40 grids in the vertical direction (i.e., horizontal grid size is and vertical size is 40 layers in the -coordinate). The length of the lock is = 0.06 metres. For the region where , the initial salinity is . For the region where , the salinity is zero. Therefore, there is a density difference between the lock and the region to its right (). All boundary conditions in the model were set as free-slip boundary conditions, and then a baroclinic calculation was carried out.

### Model validation

Numerical results were compared with those reported by Chen & Lee (2001). The contour line was taken as the outermost gravity current structure to analyze the current propagation. The forefront of the contour line was taken as the front position and velocity at this position was the measured front velocity. The front height was measured as the height of the top of this line. In the numerical simulation, the front of the gravity current is formed about 3 seconds after the lock has been released, and its height is half of the water depth. The initial front shape is shown in Figure 1(a), where *h* is the water depth and the PSU color bar is the salinity distribution legend (note that the following figures without a color bar use the same color bar as Figure 1).

As the propagation duration extends from 3 seconds to 6 seconds, the front external shape remains almost the same, and the front height is about half of the water depth. When the current has propagated for more than 6 seconds, the front height fluctuates and the maximum value is 0.0639 metres (Figure 1(b)).

During propagation, the gravity current undergoes two distinct phases: an initial slumping phase and a self-similar phase. The term ‘self-similar phase’ means that the viscous and inertial phase are contained in this phase. The front position and velocity are different functions of time in different phases (Figure 2). In the initial slumping phase, the position of the front is linearly related to , and the function is expressed as , where the coefficent 0.4746 is the Froude number of the front. This coefficient is slightly greater than the values measured in the laboratory (0.43, 0.45 or 0.46) which are cited by Chen & Lee (2001). This may be because the free-slip boundary conditions in the numerical model lead to lower damping. In the initial slumping phase, the front velocity remains almost constant (). In the self-similar phase, the position of the front in the experiment is related to raised to the power of 2/3, whereas the function derived from the numerical model is described as . This power is slightly greater than that determined in previous experiments. In the self-similar phase, the front velocity in the experiment is related to raised to the power of −1/3, whereas the function derived from the numerical model is expressed as . This power is also slightly greater than that obtained in previous experiments. Greater values of the exponents are associated with the free-slip boundary conditions.

## RESULTS

When the current has propagated for more than 6 seconds, the front shape starts to change markedly. As shown in Figure 3, the Kelvin–Helmholtz phenomenon manifests in the interface, and major vortices occur in the front of the current and their shape and size change with time. While the forefront vortex maintains its shape well, the rear vortices are intermixed.

Mixing vortices continue to advance and entrain the front vortex. The entrainment occurring in the back of the front of the current is captured by NHWAVE. The rear-most vortex is diving under the one in front of it and the latter is then forced to lift, and thus entrainment happens (see Figure 3(a)–3(c)). Entrainment of different vortices leads to the formation of a bigger vortex, which continues to proceed but does not catch up with the front-most vortex (see Figure 3(d)). In this process, the front height of the current fluctuates, and about 15 seconds after release of the lock, the front height remains unchanged and is half of the depth (see Figure 1(b)).

Figure 4 shows that there is a significant circular flow in the whole front of the current, and there are also several small eddies within the circular flow. Note that the trend of the height and the shape of the circular flows is consistent with that of the vortices described above (see Figure 3(c) and 3(d)). These flow fields provide hydrodynamic evidence of the occurrence of the vortices and of entrainment.

If the calculation of dynamics pressure is not performed, then the simulation results of a hydrostatic model can then be obtained as plotted in Figure 5. As shown in Figure 5(a) and 5(b), the front height of the current decreases after it has propagated for 10 seconds. In addition, there is vertical stratification of salinity, and the entrainment flow does not occur (see Figure 5(c) and 5(d)). In comparison with the non-hydrostatic results, the hydrostatic simulation fails to disclose the internal structure of the gravity current.

To study the effects of lock length and density differences on the propagation of a gravity current, other numerical experiments were undertaken including lock lengths of = 0.10 metres and = 0.13 metres, the other initial salinities including , which match the conditions also tested by Chen & Lee (2001).

The slumping time, when = 0.13 metres (), is about 4.6 seconds whereas when = 0.10 metres (), it is about 3.8 seconds. Both these slumping times are greater than the slumping time for the shortest lock length shown in Figure 1. Likewise the speed of the current increases with increasing lock length, namely, 1.0 metres (), 1.12 metres () and 1.18 metres () for lock lengths of = 0.06 metres, 0.10 metres and 0.13 metres respectively (see Figure 3(d)). That is because the larger lock length increases the current mass, which provides more power for the current to propagate.

A greater initial density difference between the current and ambient fluid also creates greater momentum for propagation. As shown in Figure 6 ( = 0.06 metres) and Figure 3(d), 18 seconds after lock release, the current formed under conditions with a greater density difference propagated further, with a velocity corresponding to . When the current has propagated the complete length of the horizontal domain.

## CONCLUSIONS

In this paper, a NHWAVE (non-hydrostatic wave) model has been used to study the propagation of a lock-release gravity current and to analyze the Kelvin–Helmholtz phenomenon between the current and ambient fluid. The major innovation of the paper is the simulation of interface vortices. Simulation results of the front height, position and the changes in current velocity with time are consistent with laboratory experimental results. Vortices occur between the current and ambient fluid. After the current has propagated for a while, these vortices entrain each other. Flow field results predicted by the NHWAVE model reveal the occurrence of a significant circular flow as well as small eddies within it. The circular flow allows the forward movement of fluid and the small eddies reflect the existence of interface vortices. In contrast, a hydrostatic simulation fails to capture the vortices and the front height of the current. Numerical results also show that the slumping time and propagation of current increase with increasing lock length. A larger difference of density also increases the momentum for the current, which increases its speed of propagation.

These simulation results advance the theoretical understanding of gravity currents and especially the interface entrainment. More detailed understanding of the structure of gravity currents will be advanced by the improvement of observation technology and computer technology. This research can also provide valuable reference for the study of sewage entrainment between water bodies of different density.

## ACKNOWLEDGEMENTS

The authors are very grateful to Professor Fengyan Shi of the University of Delaware for the use of the software NHWAVE. The authors also express appreciation to the anonymous reviewers for their valuable comments during the review process of the manuscript.