Effect of surface waves on full depth lock exchange gravity currents hydrodynamics

Hydrodynamics of full depth lock exchange gravity currents is investigated by means of a numerical model which takes into account both surface wave propagation and gravity currents. The adopted modeling approach couples a Boussinesq-type of model for surface waves and a gravity current model for stratified flows. Such a model is calibrated by means of experimental data, obtained both in the absence and in the presence of surface regular waves. The model provides a fairly good prediction of the heavy front position in time in both conditions. Velocity fields, heavy and light front position, turbulence and vorticity distribution are also analyzed. In particular, results highlight the fact that the presence of the oscillatory motion causes a reduction of turbulence and of mixing between heavy and light fluids. related flow from surface wave motion, by separating the orbital velocity from the gravity current velocity. Such a model was derived under the assumption of inviscid flow and free slip condition at the bottom. The validation of the numerical model was carried out by means of the experimental data of Musumeci and Foti (2014). The heavy front propagation showed a delay for the wave case with respect to the still water case, in both experimental and numerical data. Moreover the behavior of the undulations predicted by the model at the interface between the two fluids, in the rear part of the heavy front, was qualitatively in agreement with the experimental observations, notwithstanding turbulent effects were neglected in the gravity current model. However, the assumption of inviscid flow and the absence of turbulence effects caused an overestimation of the heavy front velocity propagation computed by the numerical model. In the present work, the above limit is overcome by modeling turbulence by adding up a depth uniform eddy viscosity model and a Smagorinsky (1964) model. The first formulation takes into account large scale turbulence due to the gravity current propagation. The latter formulation simulates the effects of the small scale turbulence. In detail, Section 2 describes the proposed numerical approach, focusing on the turbulence model. Section 3 describes the model setup and the comparison with experimental data. In Section 4 the results of the application of the numerical model to lock exchange release are analyzed, by focusing on hydrodynamics both in the absence and in the presence of surface waves. Finally concluding remarks are summarized in Section 5. 2 ADOPTED NUMERICAL MODELLING The velocity field u is separated into two components, uB and ud, respectively related to surface waves and density variability in space. Such components are then submitted into the continuity, Reynolds averaged and density transport equations. The first component is assumed to be independent from density gradients and it is computed by a Boussinesq-type of model, proposed by Musumeci et al. (2005). In particular, in such a model the flow is assumed rotational after breaking and the governing equations are derived with no assumptions on the order of magnitude of nonlinear effects. Moreover, in such a model the velocity is influenced by the effects of vorticity due to breaking and the possibility to consider the rotationality of the flow allows to represent mixing which occurs in the surf zone in a more realistic way. The wave propagation model solves the surface elevation and the wave-related depth-averaged horizontal velocity, by integrating the depth-averaged continuity and momentum equations. At leading order the influence of buoyancy on surface wave propagation is neglected, therefore the depth-integrated wave numerical model is applied to calculate free surface elevation and depth averaged horizontal velocities. The horizontal component of orbital velocity field uB can be retrieved by applying the quadratic horizontal velocity profile, adopted by Musumeci et al. (2005). The vertical component of orbital velocity along the water column is obtained by applying the continuity equation. The derivation of equations for estimating gravity current velocity ud is obtained similarly to Viviano et al. (2014), with the reference Cartesian system (x,z) located at the still water level, as reported in Figure 1. The density ρ varies in the range [ρ 2 , ρ 1], with ρ 1>ρ 2; waves propagate in the positive x direction by causing a variation of free surface elevation ζ and of orbital velocity uB; the local water depth h is measured from the still water level.


INTRODUCTION
Buoyancy driven flows in the presence of surface waves are common in estuarine and coastal regions of natural, urban or industrial zones. A numerical model for the interaction of gravity currents and surface waves has been implemented here for the analysis of the hydrodynamics of a buoyancy-driven flow under monochromatic waves.
Starting from Benjamin (1968), the dynamics of gravity currents has been often investigated by considering the well-known lock-exchange problem (Nogueira et al., 2013(Nogueira et al., , 2014: two fluids having different density, initially separated by a vertical gate, move toward a more stable condition with the lighter fluid at the top and the denser fluid at bottom. However, the effect of surface waves on the propagation of gravity currents has not been extensively investigated. Ng and Fu (2002) developed an asymptotic theory for the spreading of a layer of heavy immiscible liquid at the bottom under the effect of both surface wave and gravity current. Such authors considered a multiple scale perturbation analysis under the assumption that surface waves and gravity currents act at different time scales, being the sea waves faster than gravity and streaming effects. Evolution equations of profile distribution of the heavy liquid are obtained as function of slowtime variables. They obtained that the streaming current, induced by surface waves, usually predominates and drives the heavy liquid to propagate together with surface waves. In the case of standing waves, the alternating component of the streaming current can induce the development of an undulating interface between the two fluids. The formation of such an undulation can dramatically change the local structure of mass transport near the bottom and it is related to the presence of circulating cells. The characteristics of such cells, above and below the interface, are dependent from the thickness of the dense layer. Robinson et al. (2013) studied the evolution of an heavy liquid at the bottom of surface waves by means of experiments carried out in a wave tank where a finite volume of saline liquid was released. They showed that, after the initial falling, the gravity currents propagates with two fronts moving in opposite directions. The overall length of heavy release was weakly affected by wave action, being function of buoyancy. On the contrary the position of the middle point and the shape of heavy liquid were significantly affected by surface waves: i) Stokes drift dominates for long waves and the middle point of heavy liquid is moved downstream; ii) Eulerianaveraged velocity is most important near the bed for short waves and the gravity currents moves upstream, against the versus of wave propagation. Moreover they found that the centre of the current is advected at a fraction of the mean ambient flow velocity, i.e. 0.63 of the mean Lagrangian velocity at the bed. This relationship was found for a broad range of inertia-dominated gravity currents, in the presence of shallow water waves.
More recently Viviano et al. (2014) proposed a numerical model for waves and gravity currents, which has been developed by decoupling buoyancy-Effect of surface waves on full depth lock exchange gravity currents hydrodynamics A. Viviano, R.E. Musumeci & E. Foti Dept. Civil Engineering and Architecture, University of Catania, Catania, Italy ABSTRACT: Hydrodynamics of full depth lock exchange gravity currents is investigated by means of a numerical model which takes into account both surface wave propagation and gravity currents. The adopted modeling approach couples a Boussinesq-type of model for surface waves and a gravity current model for stratified flows. Such a model is calibrated by means of experimental data, obtained both in the absence and in the presence of surface regular waves. The model provides a fairly good prediction of the heavy front position in time in both conditions. Velocity fields, heavy and light front position, turbulence and vorticity distribution are also analyzed. In particular, results highlight the fact that the presence of the oscillatory motion causes a reduction of turbulence and of mixing between heavy and light fluids. related flow from surface wave motion, by separating the orbital velocity from the gravity current velocity. Such a model was derived under the assumption of inviscid flow and free slip condition at the bottom. The validation of the numerical model was carried out by means of the experimental data of Musumeci and Foti (2014). The heavy front propagation showed a delay for the wave case with respect to the still water case, in both experimental and numerical data. Moreover the behavior of the undulations predicted by the model at the interface between the two fluids, in the rear part of the heavy front, was qualitatively in agreement with the experimental observations, notwithstanding turbulent effects were neglected in the gravity current model. However, the assumption of inviscid flow and the absence of turbulence effects caused an overestimation of the heavy front velocity propagation computed by the numerical model.
In the present work, the above limit is overcome by modeling turbulence by adding up a depth uniform eddy viscosity model and a Smagorinsky (1964) model. The first formulation takes into account large scale turbulence due to the gravity current propagation. The latter formulation simulates the effects of the small scale turbulence. In detail, Section 2 describes the proposed numerical approach, focusing on the turbulence model. Section 3 describes the model setup and the comparison with experimental data. In Section 4 the results of the application of the numerical model to lock exchange release are analyzed, by focusing on hydrodynamics both in the absence and in the presence of surface waves. Finally concluding remarks are summarized in Section 5.

ADOPTED NUMERICAL MODELLING
The velocity field u is separated into two components, u B and u d , respectively related to surface waves and density variability in space. Such components are then submitted into the continuity, Reynolds averaged and density transport equations. The first component is assumed to be independent from density gradients and it is computed by a Boussinesq-type of model, proposed by Musumeci et al. (2005). In particular, in such a model the flow is assumed rotational after breaking and the governing equations are derived with no assumptions on the order of magnitude of nonlinear effects. Moreover, in such a model the velocity is influenced by the effects of vorticity due to breaking and the possibility to consider the rotationality of the flow allows to represent mixing which occurs in the surf zone in a more realistic way. The wave propagation model solves the surface elevation and the wave-related depth-averaged horizontal velocity, by integrating the depth-averaged continuity and momentum equations.
At leading order the influence of buoyancy on surface wave propagation is neglected, therefore the depth-integrated wave numerical model is applied to calculate free surface elevation and depth averaged horizontal velocities. The horizontal component of orbital velocity field u B can be retrieved by applying the quadratic horizontal velocity profile, adopted by Musumeci et al. (2005). The vertical component of orbital velocity along the water column is obtained by applying the continuity equation.
The derivation of equations for estimating gravity current velocity u d is obtained similarly to Viviano et al. (2014), with the reference Cartesian system (x,z) located at the still water level, as reported in  Gravity current velocity u d is related to density variability by means of pressure variation due to density field: p d = p -p 0 , where p 0 is the reference pressure for non-stratified flow. Starting from Reynolds-averaged Navier-Stokes momentum equations along the vertical direction, applied to both the total velocity and the orbital velocity, it is possible to obtain the term p d at a generic height z by integrating numerically the following differential equation between z and free surface elevation ζ: in which ∆ρ = ρρ 0 is the local density variation which represents the gravity current forcing, ρ 0 is the density in non-stratified flow; τ xz and τ' xz are the viscous and turbulent shear stress, which can be expressed in terms of kinematic and eddy viscosity, respectively. Such stresses have been considered here, in order to allow a more realistic representation of the gravity current, with respect to the previous inviscid version of the numerical model, proposed by Viviano et al. (2014). In particular the eddy viscosity ν t , considered in the newer version of the model, is estimated by coupling a depth-uniform eddy viscosity with a Smagorinsky (1964) formulation: where h is the local water depth; g' is the reduced gravity, defined as g(ρ 1ρ 2 )/ ρ 2 , ρ 1 and ρ 2 are the density of heavy and light flow respectively; C u and C S are two calibration coefficients related to uniform and Smagorinsky eddy viscosity formulations, respectively.
Once the pressure related to density p d is computed from numerical integration of eq. (1), the two components of gravity currents (u d ,w d ) can be obtained by solving iteratively the gravity current component of the horizontal momentum equation and of the continuity equation, which can be expressed as follows: Moreover the estimation of density is achieved by solving the density transport equation which considers both convective and diffusive effects, in order to cope with the general case of miscible fluids: where k is the density diffusion coefficient. Such a coefficient is equals to zero for immiscible fluids, here it is assumed equal to 1.5·10 -9 m 2 /s for salt diffusion.
In the proposed formulation for gravity current modelling, free surface elevation related to waves ζ B is assumed to be one order of magnitude greater than the correspondent density related free surface ζ d .
Under such an assumption wave motion can reasonably be assumed independent from gravity currents.
However the evaluation of ζ d is necessary for the estimation of buoyancy effects, since the free surface elevation modifications due to density cause the backward motion of lighter front in the lockexchange problem. Such a variable can be estimated by solving the following equation obtained from continuity equation for the gravity currents (eq.4) along with the kinematic boundary conditions at bottom (z = -h) and free surface (z = ζ ): The time stepping of eqs. (3, 5, 6) adopts the Adams-Bashfort-Multon scheme, which is third order at predictor and fourth order at corrector step, consistently with the scheme used for the solution of the Boussinesq equations. The corrector step, which is accurate up to O(Δt 4 ), is repeated until a minimum relative error is reached. Such a scheme has been chosen for its good stability properties, as specified by Press et al. (1992). Moreover the free slip condition has been considered at the bottom.
The described solution methodology is valid also in the absence of surface waves. In such a case the wave propagation model is not executed and the orbital velocities are set equal to zero.

MODEL SETUP
The adopted numerical approach is applied and calibrated for a full depth lock exchange release in the presence of periodic waves, considering experimental data obtained at the Hydraulic Laboratory of Catania University, within a wave flume 0.5m wide, 0.7m high and 9m long (see Figure 2). The numerical tank has a similar geometry, with uniform still water depth h=0.2m. At x=0m an absorbinggenerating boundary condition is located which allows to propagate waves inside the domain and to absorb the reflected waves; at the onshore side of the domain a vertical wall with a sponge layer is placed which allows to absorb the incident wave energy, thus minimizing the reflected waves. The horizontal extension of numerical tank is larger than the length of the experimental wave flume, in order to limit the effects of numerical boundary conditions. In particular the numerical tank length is equal to 14 m. Its height is 50% greater than still water depth, in order to allow free surface variations. The wave numerical model is executed for 5 wave periods in the absence of density variations, in order to reach a stable wave distribution before the onset of lock exchange. Although different from the experimental procedure, this is required by the wave generation module, and it is numerically possible thanks to the separation between the wave module and the gravity current module. The opening of the gate corresponds to the starting of gravity current in the numerical model. Initial conditions of such a model is null density currents u d , free surface elevation correspondent with that obtained by the wave model and density distribution characterized by two adjacent zones, each of them having uniform density along the water depth. In particular density is equal to greatest value ρ 1 for x=0÷7m, i.e. on the left side of the domain; in the rest of the numerical tank (x=7÷14m) density get its lower value ρ 2 .
Spatial discretization of numerical model is different in horizontal and vertical direction, with ∆x=0.05m and ∆z=0.0075m. In such a way it is possible to obtain a fairly good simulation of vertical distribution of density related velocity with more than 25 grid points along the water column.
The calibration of the numerical model is performed by considering three experimental dataset. Such tests are subdivided in simple gravity current (S1) and gravity current plus surface waves (W1, W2), whose characteristics are reported in Table 1. Both surface wave heights (H) and periods (T) vary in the surface wave tests in order to study the effects of wave slope on lock exchange propagations having similar reduced gravity g'. The calibration procedure has been carried out by varying the coefficients introduced in the eddy viscosity formulation for turbulence shown in eq. (2). In all the test the calibration procedure has provided a fairly good agreement in the estimation of heavy front propagation, between the adopted model and the experiments (see for example tests W1 and S1 in Figure 3). In particular, eddy viscosity coefficient C u is quite constant, with values in the range 0.06÷0.08. Instead the Smagorinsky formulation coefficients C S shows a larger variation, being in the range 0.04÷0.25. The latter calibration coefficient may be related to the presence and the characteristics of waves, being inversely proportional to the wave steepness. Moreover, lowest values of turbulence coefficients have been introduced in the absence of surface waves (test S1 num. W1 exp. W1 num. S1 exp. S1 Figure 3. Tests W1 and S1: comparison of heavy front position in time between experimental (black crosses and red circle) and numerical results (solid lines).

ANALYSIS OF RESULTS
Since the calibration procedure has been successfully carried out, the numerical model outputs can be used to describe the hydrodynamics in more details. Indeed the experimental measurements give only the front position in time and do not furnish measurements on hydrodynamics mechanism related to the effects of surface waves on buoyancy currents. On the contrary, the model allows to obtain the full flow field. In particular it is possible to investigate the distribution of velocity, free surface elevation, turbulence and other dependant variables such as vorticity. In particular, the latter variable (i.e. vorticity) has been here deeply investigated, since it allows to clearly highlight the interaction between the two fluids. A first analysis involves the evolution of the gravity currents after the time t 0 . Such a trigger time of the lock exchange is related to the reduced gravity g' and to the water column (h), by considering the falling time from the height h/4, and it can be considered as the end of the initial accelerated stage of the density current: Since we focus on the inertia-dominated propagation of the front, the interval between t = 0 and t = t 0 is not taken into account here. A comparison of vorticity in the presence and in the absence of waves can be carried out by means of Figures 4 and 5, respectively. In both figures the velocity field and vorticity distribution is plotted at nine instants each wave period (i.e. 0.99 s for W1 test). It is possible to note that after three wave periods, i.e. from plot (d), the differences between wave conditions and absence of waves are quite appreciable with a predominance of positive clockwise vorticity along the water column in the gate section (x = 7 m), for the wave case. 6.5 6.6 6.7 6.8 6.9 7 7.1 7.2 7.3 7.4 7.5 -0.2 6.5 6.6 6.7 6.8 6.9 7 7.1 7.2 7.3 7.4 7.5 -0.  6.5 6.6 6.7 6.8 6.9 7 7.1 7.2 7.3 7.4 7.5 -0.2 -0.1 0 0.1 Figure 5. Test S1: vorticity distribution and velocity field at nine instants, starting from t 0 and separated by the same wave period considered in test W1. Thick line shows the interface between the heavy and lighter fluid. Gate is located at x = 7m.
ω [Hz] ω [Hz] The vorticity cells in the presence of waves (test W1) are wider but the maximum intensity of vorticity is lower with respect to the test S1 without waves. The turbulence is highly related to vorticity, thus it can be argued that the presence of waves causes a reduction of mixing intensity at the interface between heavy and light liquid. Such an effect is similar to what happens in other flows characterized by wave current interaction, as it has been described by Musumeci at al. (2006).
Moreover it is evident, from the comparison of Figures 4 and 5, that the presence of waves reduces the front velocity. Such a reduction can be related to the wide extension of positive vorticity cells, which extends over the entire water column, so reducing the forcing due to buoyancy.
The effect of varying wave steepness can be investigated by comparing Figures 6 and 7, in which the vorticity and velocity field are plotted each T/4. Such a representation allows to perform an intrawave analysis, from one crest to the following.
Considering the evolution of vorticity over a wave period, it can be observed that in both cases we have two initial high positive vorticity regions, which becomes three at the end of the wave angle. However, as the waves become steeper the intensity of vorticity is slightly increased while the vertical dimension of the vortices is decreased by the effects of the orbital motion.

CONCLUSIONS AND FUTURE WORKS
A numerical model has been implemented and applied coupling an existing wave propagation Boussinesq-type of model for breaking waves with a gravity current model. The latter model, which takes into account both uniform and Smagorinsky-type eddy viscosity, allow to give fairly good approximations of buoyancy hydrodynamics notwithstanding its simplicity and numerical efficiency.
The proposed numerical model for the simulation of gravity currents under surface waves has been calibrated considering experimental tests, with comparable reduced gravity. The calibration coefficients are represented by the constants in eddy viscosity formulations for Reynolds stresses.
The effect of the presence of waves has been studied by comparing the results of a surface wave test with the test carried out in the absence of surface waves.  Figure 7. Test W2: vorticity distribution (colored contours), interface between the fluids(thick line) and velocity field (vectors), at five intervals of the wave period: (a) t = 3.4s; (b) t = 3.4s + T/4; (c) t = 3.4s + T/2; (d)t=3.4s +3T/4; (e) t = 3.4s + T.
ω [Hz] ω [Hz] Such a calibrated model has been used as a tool for studying the hydrodynamics of gravity currents lock exchange, also in the presence of oscillatory motion. Such an analysis reveals that the vorticity is strongly dominated from buoyancy. However, the presence of waves causes appreciable modifications to the vorticity distribution which start to be significant about four wave periods after the initial transient. In particular a reduction of intensity of vorticity and an increase of vortex dimensions can be noted in the presence of surface waves. As a consequence, the presence of wave motion causes a reduction of mixing of gravity current.
The effect of wave steepness on buoyancy has been studied by comparing the vorticity under two wave tests having different wave period and height. Such an analysis has been carried out over a wave perios in order to focus on the waves phases. It has been shown that the increase of wave stepness causes a reduction of vortex dimensions and an increase of maximum absolute values of vorticity, so bringing the system toward a condition similar to the absence of waves. Such a result have to be investigated more in depth by extending the comparison with a wider range of wave conditions. Moreover the adopted turbulence formulation will be upgraded in the future by considering more complex k-ε or k-ω formulations.