2.1. The lock-exchange set-up

In this study, we aim to understand the behaviour of gravity currents in the presence of external forcing by an oscillatory uniform horizontal pressure gradient using a 2-D numerical set-up. We use the classical lock-exchange set-up (see a simplified sketch in figure 1b), which is commonly used in DNS studies (Härtel et al. Reference Härtel, Meiburg and Necker2000; Necker et al. Reference Necker, Härtel, Kleiser and Meiburg2002, Reference Necker, Härtel, Kleiser and Meiburg2005; Cantero et al. Reference Cantero, Balachandar, García and Ferry2006). Initially, there is heavy fluid (with density $\rho _1$) in one side of the domain (grey area in figure 1b) and there is light fluid (with density $\rho _0$) in the other side of the domain. Both fluid layers extend over the entire height $H$ of the channel. There is a gate between the heavy ($\rho _1$) and light ($\rho _0$) fluid, which is instantaneously removed at time $t=0$. A sketch of the gravity current interface separating the heavy and light fluid after a certain time from the removal of the gate is schematically shown with a black line in figure 1(b).

2.2. Governing equations

The equations of motion for fluid flow with density variations in the presence of gravity are the continuity, Navier–Stokes and density transport equations. For fluids with small density differences, the Boussinesq approximation can be applied to the governing equations. In this study, we are interested in gravity currents caused by small density differences of a few percent which safely fall within this approximation. Therefore, we apply this Boussinesq approximation simplifying the governing equations by keeping the density constant in the Navier–Stokes equation except in the buoyancy term. The resulting equations from this approximation are called the Boussinesq equations. These equations formulated using Einstein's notation in a Cartesian coordinate system ${x}_i=({x},{y})$, with ${x}$ and ${y}$ the streamwise and wall-normal directions, respectively, are

(2.1)\begin{gather} \frac{\partial {u}_j }{\partial {x}_j} = 0, \end{gather}

(2.2)\begin{gather} \frac{\partial {u}_i }{\partial t} + \frac{\partial ( {u}_i {u}_j)}{\partial x_j} =- \frac{1}{ {\rho }_0} \frac{\partial {p}}{\partial {x}_i} + {\nu} \frac{\partial^2 {u}_i }{\partial {x}_j \partial {x}_j} - \frac{{\tilde{\rho}}}{ {\rho }_0} {g} \delta_{i2} + {U}_0 {\omega} \cos{( {\omega} {t} )} \delta _{i1}, \end{gather}

(2.3)\begin{gather} \frac{ \partial {\rho}}{\partial t} + \frac{ \partial ( {\rho} {u}_j)}{\partial {x}_j} = {\alpha} \frac{ \partial ^2 {\rho}} {\partial {x}_j \partial {x}_j}. \end{gather}

Here, ${u}_i=({u},{v})$ denotes the Cartesian velocity components, ${p}$ the pressure, $g$ the gravitational acceleration acting in the negative $y$ direction, ${\rho _0}$ the reference density or the density of the light fluid, $\tilde {\rho }$ the total density, $\rho = \tilde {\rho } - \rho _0$ the density variable with respect to $\rho _0$, $\delta _{i1}$ and $\delta _{i2}$ are the Kronecker delta symbols, ${\nu }$ the kinematic viscosity, and ${\alpha }$ the molecular diffusivity of the density field.

The last term on the right-hand side of (2.2) represents a horizontal oscillating pressure gradient force given by

(2.4)\begin{equation} {-}\frac{1}{\rho_0} \frac{{{\rm d}} {{p_0}}}{{{\rm d}}{\kern0.9pt}{{x}}} = {U}_0 {\omega} \cos{( {\omega} {t}) }, \end{equation}

where $p_0$ represents an externally imposed pressure distribution, ${U}_0$ denotes a free-stream velocity amplitude, ${\omega }= 2 {\rm \pi}/T_{osc}$ is the frequency of oscillations (with $T_{osc}$ the period of oscillations) and $t$ represents time. This pressure gradient results in a sinusoidal free-stream velocity $U_0 (t) = U_0 \sin (\omega t)$, as shown in the inset of figure 2, where the phase of the oscillation, $\phi = 360^{\circ } (t/T_{osc})$, is shown in degrees. The horizontal free-stream velocity reaches the maximum velocity at $\phi = 90^{\circ }$ and $\phi = 270^{\circ }$, and the vertical distribution of the velocity for $\phi = 90^{\circ }$ is shown in figure 2 for different values of the oscillation period (quantified by $KC_b$ defined below).

Figure 2.

Horizontal velocity profile caused by an oscillatory horizontal pressure gradient as a function of the channel depth at $\phi =90^{\circ }$ (the inset shows that this is the instant when maximum positive sinusoidal free-stream velocity ${U}_0$ is reached). The Stokes boundary-layer thickness grows with increasing $KC_b$.

Dimensional quantities that define the flow are shown in table 1. Dimensional analysis results in four dimensionless numbers that govern the flow. The Reynolds number $Re$, defined as

(2.5)\begin{equation} Re = \frac{{u}_b {H}} {{\nu}}, \end{equation}

represents the ratio of the inertial to viscous forces. Here, $u_b$ is the buoyancy velocity,

(2.6)\begin{equation} {u}_b=\sqrt{{g} ' {H}}, \end{equation}

where ${g} '$ indicates the reduced gravitational acceleration defined as

(2.7)\begin{equation} {g} ' = {g} \frac{\Delta{\rho}} {{{\rho }}_0}, \end{equation}

with $\Delta {\rho } = \rho _1 - \rho _0 >0$ the density difference between the two fluids. The Schmidt number,

(2.8)\begin{equation} Sc = \frac{{\nu}} {{\alpha}}, \end{equation}

represents the ratio of kinematic viscosity of the fluid to molecular diffusivity of the density field. The Froude number $Fr$, defined as

(2.9)\begin{equation} Fr = \frac{{U}_0 } {{u}_b }, \end{equation}

is the ratio between the free-stream velocity amplitude of oscillation ${U}_0$ to the buoyancy velocity ${u}_b$. Lastly, we introduce the Keulegan–Carpenter number ($KC_b$), which represents the ratio between the oscillation period ($T_{osc}$) to the advection time scale associated with the propagation of the gravity current (${H}/{{u}_b}$). It is given as

(2.10)\begin{equation} KC_b = \frac{ {T}_{osc} {u}_b}{{H}}. \end{equation}

The Keulegan–Carpenter number $KC_b$ can also be interpreted as a dimensionless period of the oscillations of the externally applied flow field.

Table 1.

Dimensional quantities that govern the fluid flow for our set-up.

We will use the governing equations (2.1)–(2.3) in dimensionless form by using

(2.11a–e)\begin{equation} {u_i}^* = \frac{{u_i}}{{{u}_b} }, \quad {x_i}^* = \frac{{x_i}}{{H} }, \quad t^* = \frac{{t}}{{H/{u}_b} }, \quad p^* = \frac{{p}}{{\rho }_0 \, {{u}_b^2} }, \quad \rho ^* = \frac{ \tilde{\rho}} { \Delta {\rho} } \end{equation}

to make velocity, length, time, pressure and density dimensionless. Henceforth, we will use the dimensionless form of the equations and we eliminate the asterisk for convenience:

(2.12)\begin{gather} \frac{\partial {u}_j }{\partial {x}_j} = 0, \end{gather}

(2.13)\begin{gather} \frac{\partial {u}_i }{\partial t} + \frac{\partial ( {u}_i {u}_j)}{\partial x_j} =- \frac{\partial {p}}{\partial {x}_i} + \frac{1}{Re} \frac{\partial^2 {u}_i }{\partial {x}_j \partial {x}_j} - {\rho} \delta _{i2} + \frac{Fr}{KC_b} 2{\rm \pi} \cos \left( \frac{2{\rm \pi} t}{KC_b} \right) \delta _{i1}, \end{gather}

(2.14)\begin{gather} \frac{ \partial {\rho}}{\partial t} + \frac{ \partial ( {\rho} {u}_j)}{\partial {x}_j} = \frac{1}{Re Sc} \frac{ \partial ^2 {\rho}} {\partial {x}_j \partial {x}_j}. \end{gather}

The forcing of the flow above the flat no-slip bottom wall by the oscillatory horizontal pressure gradient introduces a shear layer (the Stokes boundary layer) adjacent to the bottom boundary and an oscillating uniform free-stream velocity in the bulk above this layer. For later use (§ 4), it is instructive to compare the evolution of a gravity current in each of the limiting cases ($KC_b\rightarrow 0$ and $KC_b\rightarrow \infty$) with the freely evolving gravity current. On the buoyancy time scale, the propagating gravity current experiences a uniform quasi-steady shear flow for $KC_b \rightarrow \infty$, while it does not, on average, experience a free-stream flow field for $KC_b \rightarrow 0$ since $T_{osc}\ll H/u_b$.

Let us first focus on the latter case, $KC_b\rightarrow 0$. For high-frequency oscillations, the shear layer generated by the applied oscillating pressure gradient becomes very thin (see figure 2) and the duration of a single oscillation cycle is very short compared with the typical advection time scale associated with the buoyancy-driven dynamics of the gravity current. We expect that these high-frequency oscillations hardly affect the evolution of the gravity current on the buoyancy time scale, leaving its gross features almost untouched. The externally imposed oscillatory flow can nevertheless significantly affect the fluid motion near the stagnation point at the front of the current. The impact on details of the front dynamics and density redistribution will then always be significantly different compared with the freely evolving case, as will be illustrated in § 4.

However, when $KC_b\rightarrow \infty$, the period of oscillations of the applied pressure gradient increases, and the Stokes boundary layer thickness grows and induces an ever thicker region with (quasi-steady) shear flow adjacent to the wall. The gravity current is then for a significant part embedded in and affected by this (quasi-steady) shear flow that, depending on the phase $\phi$ of the forcing cycle, has a free-stream velocity in the propagation direction of the gravity current ($\phi =90^{\circ }$) or in the opposite direction ($\phi =270^{\circ }$). Summarizing, the gross features of the dynamics of the freely evolving gravity current will show similarities with those of the externally forced cases with $KC_b$ small (high-frequency limit) and behaves completely differently compared with oscillatory-forced gravity currents with $KC_b$ very large (low-frequency limit). Visual evidence of these distinct behaviours for the low- and high-frequency ranges will be provided in § 4.

2.3. Parameter settings for our numerical studies

The dimensionless numbers governing freely evolving gravity currents are $Re$ and $Sc$. The total number of grid points of a DNS of such 2-D gravity currents, for fixed $Sc$ and assuming proper resolution of the small-scale turbulence, scales with $O(Re)$ (while for 3-D gravity currents, we expect a scaling of the number of grid points with $O(Re^{9/4})$). Restrictions on the time step typically provides an additional cost proportional with $\sqrt {Re}$. Therefore, DNS studies are mostly limited to laboratory scale $Re$ numbers, even more so for 3-D turbulent gravity currents. Nonetheless, the results from such DNS studies are used to understand a variety of physical processes in more detail and contribute to the development of large-scale models suitable for higher $Re$ applications. For our DNS, we choose a fixed value of $Re=3000$. This selection allows for a feasible parameter study while aligning with the order of magnitude employed in earlier DNS studies of gravity currents (Härtel et al. Reference Härtel, Meiburg and Necker2000; Necker et al. Reference Necker, Härtel, Kleiser and Meiburg2002, Reference Necker, Härtel, Kleiser and Meiburg2005; Cantero et al. Reference Cantero, Balachandar, García and Ferry2006).

The Schmidt number displays significant variation, across several orders of magnitude for different environmental flows, since the diffusivity of the scalar that causes the density variation can differ significantly. For saline water, $Sc \simeq 700$. In the case of thermal density differences in water, $Sc\simeq 7$, while in air, $Sc \simeq 0.7$ (which is also known as the Prandtl number). However, the lower range $Sc$ are mostly used for DNS studies of gravity currents to decrease the computational cost (typically $Sc= O(1)$). This is justified by Necker et al. (Reference Necker, Härtel, Kleiser and Meiburg2002) indicating that the integral properties of the current are independent of the actual $Sc$ value as long as it is not much smaller than one. Moreover, Birman, Martin & Meiburg (Reference Birman, Martin and Meiburg2005) stated that varying $Sc$ between 0.2 and 5 has an influence although it is small. For turbulent flows, the turbulent Schmidt number (the ratio of turbulent viscosity and turbulent diffusion coefficient) is typically around 0.7 to 1, while it increases for strongly stratified flows up to $Sc \simeq 5$. Ralston, Geyer & Lerczak (Reference Ralston, Geyer and Lerczak2008) found that the effective $Sc \simeq 2$ for a salt wedge in an estuary. Donzis et al. (Reference Donzis, Aditya, Sreenivasan and Yeung2014) showed, based on a high-resolution homogeneous isotropic turbulence simulation, that the turbulent Schmidt number is between 1 and 2 for Péclet number ($Pe = Re \, Sc$) higher than 300, which is consistent with the findings by Necker et al. (Reference Necker, Härtel, Kleiser and Meiburg2002). In our study, we use $Sc=5$ in the majority of our simulations, but we also investigate the effect of $Sc$ by varying it in the range 1–20 with the aim of understanding the effect of diffusion on our results.

The external forcing is characterized by the Froude and Keulegan–Carpenter numbers, with the free-stream velocity amplitude characterized by $Fr$ and the oscillation period by $KC_b$. The tidal Froude number, $Fr_T$, is the ratio between the maximum velocity in the tidal cycle and buoyancy velocity. It varies for different estuaries. In some estuaries, the tidal effect is smaller, and hence, $Fr_T$ becomes smaller. Typical values reported in the literature are $Fr_T \approx 0.35$ for the Mississippi River, $Fr_T \approx 0.4$ for the Connecticut River Estuary and $Fr_T \approx 0.45$ for the Tweed River (Ralston, Geyer & Lerczak Reference Ralston, Geyer and Lerczak2010). Meanwhile, $Fr_T$ is higher for some other estuaries. Examples are $Fr_T \approx 0.7$ for the Columbia River Estuary, the Fraser River and the Merrimack River Estuary, and $Fr_T \approx 0.8$ for the Snohomish River Estuary (Geyer et al. Reference Geyer, Lavery, Scully and Trowbridge2010). Therefore, we vary $Fr$ between 0.25 and 2 to understand the effect of the amplitude of the external oscillatory forcing on gravity currents.

Finally, we need to discuss an appropriate range of $KC_b$ for our numerical simulations. Based on previous DNS studies (Härtel et al. Reference Härtel, Meiburg and Necker2000; Necker et al. Reference Necker, Härtel, Kleiser and Meiburg2002, Reference Necker, Härtel, Kleiser and Meiburg2005; Cantero et al. Reference Cantero, Balachandar, García and Ferry2006) on canonical freely evolving gravity currents, distinct characteristics are observed. Specifically, at the leading edge of the gravity current, a pronounced front emerges. This front is characterized by a steep increase in the current height, followed by a subsequent decrease towards the body of the gravity current. The location where the maximum current height occurs at the front is referred to as the head position. However, in nature, gravity currents are not always freely evolving and can be subjected to different external forces including an external oscillatory forcing. Salt wedges in estuaries are influenced by surface waves which typically have periods of the order of seconds as well as tidal forcing with a period of 12.42 or 24.84 hours on average. Recently, the effect of surface waves on gravity currents has been studied by using LES and $k-\epsilon$ models and experiments, and it was demonstrated that the presence of oscillatory forcing modifies the (turbulent) mixing happening at the front of the fresh-water current near the surface (Stancanelli et al. Reference Stancanelli, Musumeci and Foti2018a,Reference Stancanelli, Musumeci and Fotib; Cui et al. Reference Cui, Kamath, Wang, Han and Bihs2022; Marino et al. Reference Marino, Stagnitti, Stancanelli, Musumeci and Foti2023). However, the effect of surface waves on the heavier current near the bottom was less significant. The density distribution of the heavy front was actually similar in the absence of surface waves. Measurements of the salt intrusion in the Rotterdam Waterway (de Nijs, Pietrzak & Winterwerp Reference de Nijs, Pietrzak and Winterwerp2011) indicate significant differences in the front characteristics of the heavy current compared with freely evolving gravity currents. The front characteristics are also different compared with those of gravity currents affected by surface waves (Stancanelli et al. Reference Stancanelli, Musumeci and Foti2018a,Reference Stancanelli, Musumeci and Fotib; Marino et al. Reference Marino, Stagnitti, Stancanelli, Musumeci and Foti2023). For tidal forcing in the Rotterdam Waterway, with $H\approx 15\ {\rm {m}}$ and $u_b=1.5\unicode{x2013}1.8\ {\rm m}\ {\rm s}^{-1}$ (depending on the salt concentration), the main tidal constituent (M2) with a period of 12.42 h typically gives $KC_b \simeq 5000$. However, computationally, it is currently not feasible to use $KC_b$ numbers as high as 5000, since such simulations would require a domain with a very large aspect ratio and they would need very long integration times. We will show later that such simulations are fortunately not needed. Since distinct characteristics are observed between freely evolving gravity currents (Härtel et al. Reference Härtel, Meiburg and Necker2000; Cantero et al. Reference Cantero, Balachandar, García and Ferry2006) and the salt wedge in the Rotterdam Waterway (de Nijs et al. Reference de Nijs, Pietrzak and Winterwerp2011), where the oscillation period is large, it is necessary to explore the effect of different periods by varying $KC_b$. We will start with smaller periods of oscillation and increase the oscillation period to reach $KC_b$ values which can show similar flow features as those observed in estuarine salt intrusions. For this purpose, we will use $KC_b$ values of 5, 10, 25, 50 and 100.

3.1. Numerical model

A variety of approaches are available for computational modelling of gravity currents, see, for example, the review by Meiburg, Radhakrishnan & Nasr-Azadani (Reference Meiburg, Radhakrishnan and Nasr-Azadani2015) and more recently the study by Van Reeuwijk, Holzner & Caulfield (Reference Van Reeuwijk, Holzner and Caulfield2019) using SPARKLE (Craske & Van Reeuwijk Reference Craske and Van Reeuwijk2015). We perform 2-D DNS to solve the governing equations (2.12)–(2.14) in non-dimensional form using the incompressible Navier–Stokes solver PARTIES (Biegert, Vowinckel & Meiburg Reference Biegert, Vowinckel and Meiburg2017a; Biegert et al. Reference Biegert, Vowinckel, Ouillon and Meiburg2017b). The code is based on a finite-difference approach, employing the fractional step method by Kim & Moin (Reference Kim and Moin1985), and uses a third-order explicit Runge–Kutta scheme with three substeps (Harten Reference Harten1997) to discretize the equation in time. The projection method (Chorin Reference Chorin1968) is used to ensure incompressibility, and the resulting Poisson equation is treated using a direct fast Fourier transform (FFT) solver at each Runge–Kutta substep (Biegert et al. Reference Biegert, Vowinckel and Meiburg2017a). The convective term in the momentum equation is explicitly solved, leading to the well-known Courant–Friedrichs–Lewy restriction. Spatial derivatives of the diffusion term in momentum equations are assessed using a second-order central finite differences scheme. A second-order upwind scheme is used to handle convective terms in both the momentum and transport equations. The diffusion term in both the momentum equation and the density transport equation is treated implicitly using second-order central differences in combination with a conjugate-gradient solver (Saad Reference Saad2003). Furthermore, the code has been parallelized using the MPI library to enhance performance. The code's accuracy and reliability have been rigorously validated through extensive testing, as demonstrated by Nasr-Azadani & Meiburg (Reference Nasr-Azadani and Meiburg2011).

3.2. Numerical set-up, boundary and initial conditions

Initially, the entire fluid is in a state of rest. At $t=0$, a (virtual) gate is removed and due to buoyancy, a gravity current develops. Additionally, an external forcing in the form of an oscillatory pressure gradient is imposed in the along-channel direction. To apply periodic boundary conditions in the along-channel direction for pressure, velocity and density, we propose a computational set-up as sketched in figure 3. The middle section of the domain contains heavier fluid with a density $\rho = 1$ (the grey area in figure 3), while the left- and right-hand sides of the domain contain light fluid with a density $\rho = 0$. Our primary focus will be on the developing gravity current in the right half of the domain (in particular, the region with $x\ge 0$). The phase of the free-stream velocity, resulting from the imposed oscillatory pressure gradient, is initialized at $\phi =0^{\circ }$ (see figure 2). On the left-hand side of the domain ($x<-170$), we effectively replicate (at $t=0$) the initial density distribution as on the right-hand side ($x>-170$), but the gravity current (moving to the left) experiences an effective free-stream velocity with a $180^{\circ }$ phase difference compared with the right-hand side. This deliberate configuration provides us with a valuable opportunity to contrast the effect of oscillatory forcing starting from different phases. In the subsequent discussions, we will refer to the right part of the channel, except when the left part is explicitly mentioned. A no-slip boundary condition at the bottom wall and a stress-free boundary condition at the top boundary of the channel (a rigid free surface) are employed. The rigid free surface does not always accurately represent the free-surface behaviour for gravity currents in laboratory studies or estuaries, but it is commonly used for high-resolution numerical simulations of gravity currents.

Figure 3.

Initial density field of the simulation set-up. The primary focus is the gravity current forming on the right-hand side of the domain ($x>-170$). A visual representation of the interface between heavy and light fluid for the gravity current, which will form at a certain time after removal of the gate, is shown with a solid line. The gravity current at the left-hand side experiences oscillatory forcing with a phase difference of $180^{\circ }$ compared with the right-hand side.

The length of the domain is selected such that the gravity current does not extend to the sides of the domain during the full simulation covering 400 dimensionless time units. To include a minimum of four oscillation cycles of the externally imposed pressure gradient for the lowest frequency case ($KC_b=100$), we need to choose the length of our channel equal to $L=680$ (and we use this for all the simulations). We use an equidistant Cartesian mesh with $\delta x = \delta y = l$. According to our convergence study, a non-dimensional grid size of $l=0.008$ for simulations with $Sc=1$ and $5$ is sufficient, while a grid size of $l=0.0025$ is necessary for simulations with $Sc=20$. Therefore, with our domain size of $680\times 1$ dimensionless units, we need approximately $10^7$ grid points for runs with $Sc=1$ and $5$, while we need approximately $10^8$ grid points for the simulations with $Sc=20$.

We performed our simulations on the Dutch Super Computer Snellius, using 128 cores. The simulations over 400 dimensionless time units using $l=0.008$ required approximately 9000 CPU hours (approximately 71 wall-clock hours). For simulations with $l=0.0025$, the computational cost is nearly 528 000 CPU hours. We will cover a wide range of the parameter space by varying the three dimensionless numbers $KC_b$, $Fr$ and $Sc$. In the next sections, we present results from 13 simulations whose parameters are given in table 2, and we discuss (first qualitatively and later on quantitatively) the effect of varying these three parameters on the evolution and oscillatory dynamics of the gravity current. (For fresh and saline water, with $g'\approx 10^{-2}\ {\rm g}$, we obtain the following typical lab scale values: $H\approx 5$ cm and $u_b\approx 7\ {\rm cm}\ {\rm s}^{-1}$ based on $Re=3000$.)

Table 2.

The parameters for the simulations in this study. $Re=3000$ is maintained for all simulations. In the next section, simulations 1 to 6 are compared to evaluate the effect of different forcing periods (by varying $KC_b$). To examine the impact of different ambient velocity amplitude ($Fr$), results from simulations 1, 5, 7, 8 and 9 are compared for $Sc=5$. Additionally, the influence of diffusivity is analysed by contrasting results from simulation 1, 5, 10, 11, 12 and 13, where three different $Sc$ numbers are employed for two different $KC_b$ numbers. The simulations 1, 10 and 12 represent freely evolving gravity currents.

4.1. Evolution of gravity currents during the forcing cycle

We start our exploration with an illustration of the gravity current evolution and the change in the (fine-scale) density structure during one complete forcing cycle. The redistribution of the density occurring at the front of the gravity current can be illustrated with snapshots of the density field at different phases of the forcing cycle (see also the animation of the density field evolution for $50\le t\le 200$ in supplementary movie 1). Here, we consider the simulation with $KC_b=50$ (figure 4). In the beginning of the second oscillation cycle ($\phi =0^{\circ }$, starting at $t=50$, see figure 4a), the gravity current propagates into the fresh-water region near the bottom of the channel. In the first half of the oscillation cycle ($0^{\circ } \le \phi \le 180^{\circ }$), the gravity current propagates in the same direction as the imposed oscillatory flow. For $0^{\circ } \le \phi \le 90^{\circ }$, the imposed ambient flow tends to accelerate the gravity current and to decelerate it when $90^{\circ } \le \phi \le 180^{\circ }$. The imposed flow has a height dependent velocity in the (Stokes) boundary layer, which promotes differential advection. This is clearly visible for $\phi =90^{\circ }$ in figure 4(b). Heavy fluid has moved over lighter fluid near the bottom wall. We call this the lifting of the gravity current, and the area underneath the heavy current the lifting area. During the deceleration phase ($90^{\circ } \le \phi \le 180^{\circ }$), the lifted part of the gravity current becomes unstable and starts to promote eddy formation, inducing enhanced (vertical) redistribution of the density field (see supplementary movie 1; $67\lesssim t\lesssim 77$ and $125\lesssim t\lesssim 135$). During the first half of the oscillation cycle, the thickness of the front of the gravity current is thinner (see figure 4a,b,e) compared with the freely evolving gravity current (shown in figure 4f).

In the second phase of the oscillation cycle ($180^{\circ } \le \phi \le 360^{\circ }$), the externally applied flow acts against the propagation direction of the gravity current. This situation no longer supports lifting. Since an unstable flow configuration already emerged during the decelerating phase ($90^{\circ } \le \phi \le 180^{\circ }$), the flow reversal provokes strong (vertical) density redistribution (figure 4c,d). As can be seen in supplementary movie 1, the heavy fluid sinks while the lighter fluid rises towards the upper interface of the gravity current. This lighter fluid subsequently moves further backwards towards the body of the gravity current due to relatively strong horizontal advection. This process takes place predominantly for $180^{\circ } \le \phi \le 270^{\circ }$ when the externally applied opposing flow accelerates to its maximum value at $\phi =270^{\circ }$ (figure 4d).

4.2. Coincidence of the imposed flow with the current propagation $(0^{\circ }\le \phi \le 180^{\circ })$

From the previous description, we derive two distinct behaviours regarding the (small-scale) redistribution of the density field in the gravity current depending on the relative direction of the external forcing with respect to the propagation direction of the gravity current. When the oscillatory forcing (thus the externally imposed flow) coincides with the propagation direction of the gravity current, the gravity current is subjected to differential advection. This effect is most pronounced for $\phi =90^{\circ }$, when the oscillatory forcing reaches a maximum value in the direction of the gravity current. The behaviour of the gravity current is completely different when it is subjected to opposing flow (with a maximum at $\phi =270^{\circ }$). This dissimilarity will be further explored in § 4.3, with density field snapshots at $\phi =360^{\circ }$ (or equivalently $\phi =0^{\circ }$) to illustrate the significant changes in the density distribution induced by the opposing flow. We shall first discuss the first half of the oscillation period.

Snapshots of the different density fields provide us with a qualitative picture of how differential advection affects the gravity current as a function of $KC_b$ (see also the animation of the density field evolution for $0\le t\le 200$ in supplementary movie 2). Before we proceed, it is useful to discuss the approach to present our data in this section. We compare the density fields at the same phase of forcing, here at $\phi =90^{\circ }$, for simulations covering all $KC_b$. Therefore, we first select a time instance where, for all considered $KC_b$, an oscillation cycle concludes, thus $\phi =360^{\circ }$. Subsequently, we select instances that coincide at $\phi =90^{\circ }$ for the different $KC_b$, during the preceding oscillation cycle. By choosing $t=100$, we compare then the density fields at $\phi =90^{\circ }$ for $KC_b=5$, 10, 25, 50 and 100 at $t=96.25$, 92.5, 81.25, 62.5 and 25, respectively. For the canonical lock-exchange gravity current (without external forcing), we use the density field at $t=96.25$.

The density field of the gravity current forming from the classical lock-exchange set-up (see figure 5a) illustrates the same characteristics as those observed in earlier DNS studies (Härtel et al. Reference Härtel, Meiburg and Necker2000; Necker et al. Reference Necker, Härtel, Kleiser and Meiburg2002, Reference Necker, Härtel, Kleiser and Meiburg2005; Cantero et al. Reference Cantero, Balachandar, García and Ferry2006; Anjum, McElwaine & Caulfield Reference Anjum, McElwaine and Caulfield2013). At the leading edge of the gravity current, the height of the current increases steeply to a local maximum and subsequently to a local minimum in the body of the gravity current. The density fields at $\phi =90^{\circ }$ for gravity currents in the presence of oscillatory forcing (figure 5b–f) show distinct shapes when compared with the freely evolving gravity current (figure 5a). Quite remarkably, the gravity currents for $KC_b\le 25$ have propagated a larger distance compared with the non-oscillatory case. Since we have a no-slip bottom, the magnitude of the externally imposed flow field must decrease towards the bottom and differential advection occurs. Heavier fluid farther from the bottom experiences a larger along-stream flow velocity than fluid near the bottom wall. Consequently, it moves over lighter fluid close to the bottom wall, particularly near the nose of the gravity current. This results in an apparent lifting of the gravity current, and the area beneath the gravity current is then called the lifting area. The lifting of the gravity current usually does not occur, or is almost negligible, for the freely evolving gravity current over a horizontal wall. The significant lifting is thus entirely due to the presence of oscillatory forcing and subsequent differential advection. We clearly observe the lifting and squishing of the gravity currents in figure 5(b–f). Lifting is weak for $KC_b=5$, but it increases with $KC_b$ because differential advection occurs for longer as the oscillation cycles become longer (thus larger $KC_b$). Additionally, the Stokes boundary-layer thickness becomes larger with increasing $KC_b$ (see figure 2). The layer of slowly moving light fluid becomes thicker, contributing to larger redistribution of light and heavy fluid beneath the heavy fluid at the front. For a quantitative discussion, see §§ 5.1 and 5.2.

In figure 6, we show the same density fields as in figure 5, but zoomed-in on the front of the gravity current. The horizontal and vertical length scales in the panels of figure 6 are identical (in contrast with figure 5) to appreciate the fine-scale structure within and beneath the gravity current, and the details of the front of the gravity current, including distinctly different (small-scale) redistribution of the density for the different cases.

In some of the snapshots shown in figure 6, in particular for $KC_b=100$ and to a lesser extent for $KC_b=50$, there seems to be evidence of the presence of a Rayleigh–Taylor instability (see also figure 7b,c,d for $KC_b=5$, 10 and 25, respectively). For a closer inspection, the reader is also referred to the supplementary movies. This is an intriguing observation that deserves further exploration in future studies of both 2-D and 3-D DNS of oscillatory-forced gravity currents.

4.3. Gravity current propagation with opposing ambient flow ($180^{\circ }\le \phi \le 360^{\circ }$)

We now consider the situation where the propagating gravity current is subjected to an opposing externally imposed flow. This occurs for $180^{\circ }\le \phi \le 360^{\circ }$. The maximum amplitude of the opposing flow occurs at $\phi =270^{\circ }$, but the influence of the opposing external flow is best appreciated when we consider its impact on advection and local redistribution of the density field at the end of the cycle with the opposing flow, i.e. at $\phi =360^{\circ }$. Hence, we compare the flow and density fields for each $KC_b$ always at the same instance after the release of the gravity current. The choice of the fixed times is dictated by the largest value of $KC_b$ and is typically $t=100$, 200, 300 or 400. We choose to illustrate advection and density redistribution with snapshots of the density field taken at $t=100$ for all cases (non-oscillating and $KC_b =$ 5, 10, 25, 50 and 100, respectively), see figure 7. These snapshots show a completely different structure compared with the case with a coinciding ambient flow.

Before proceeding with the qualitative comparison of the behaviour of oscillatory-forced gravity currents, it is necessary to address the representativeness of the density fields obtained during the relatively early stages after the release of the gravity current (for $t\le 100$). The freely evolving gravity current propagates in a quasi-steady fashion, only driven by buoyancy due to the density difference between light and heavy fluid. Considering snapshots of the density field for this case taken at different instances, here at $t=100$ (see figure 7a), 200, 300 and 400 (see figures S10, S11 and S12 of the Supplementary Material), it can be observed that the front shape remains similar for all times, while the heavy current (front and body) is diluted due to entrainment of light fluid (and scalar diffusion). For the same times, the front of the gravity current in the presence of oscillatory forcing shows distinctively different patterns of small-scale density redistributions depending on the particular value of $KC_b$, but for each value of $KC_b$, this pattern of density redistribution does not strongly change with increasing time (except dilution by also diffusive processes). Signatures of flow instabilities (shear-induced eddy formation) and density-current lifting is clearly visible for $KC_b=5$, 10 and 25 (see figure 7b–d for $t=100$ and also the snapshots shown in figures S1–S9 of the Supplementary Material for $t=200$, 300 and 400, respectively). The main features of the evolution of the gravity current for $KC_b=50$ and 100 remain visible up to $t=400$. From this comparison, we can conclude that for the purpose of this qualitative comparison, we can rely on the density fields evaluated in the oscillation cycle ending at $t=100$.

We first consider $KC_b\le 10$ (figure 7b,c). The externally induced oscillating flow acts against the propagation of the gravity current and could eliminate effects of density-current lifting by advection and enhanced (vertical) density redistributions. The duration of the opposing flow is apparently too short for these smaller values of $KC_b$ to achieve substantial (vertical) redistribution of density. The front of the gravity current thus always stays lifted and the interface (both beneath and above the current) is disturbed by eddies due to the presence of shear layers. Clearly, no full vertical density redistribution is achieved. Gravity currents in the presence of large-period oscillatory forcing (e.g. those with $KC_b=50$ and $100$, in figure 7e,f) experience strong changes in the density structure resulting in significant horizontal homogenization of the density field. Shear instabilities and eddies at the interface between heavy and light fluid have ample time to fully redistribute heavy and light fluid near this interface. These gravity currents are horizontally well homogenized at $t=100$, and their shapes are clearly different from the front shape of gravity currents with $KC_b\lesssim 25$ (figure 7b–d). Furthermore, the front of the gravity current has a clear inclined shape for $KC_b \ge 50$, which is similar to the inclined shape of salt intrusions at the ebb phase, as observed with salinity measurements in, for example, the Rotterdam Waterway reported by de Nijs et al. (Reference de Nijs, Pietrzak and Winterwerp2011). Therefore, we conclude that for $KC_b \gtrsim 50$, the behaviour of density redistributions at the front and the shape of the gravity current seem representative for salt intrusions as observed in environmental flows (usually in the presence of large-period forcing, typically with $KC_b \simeq 5000$). Lastly, for $KC_b=25$, see figure 7(d), the front is not very homogenized compared with for $KC_b=50$ and $100$, although the front clearly has a more inclined shape compared with the cases with $KC_b=5$ and 10.

4.4. The evolution of the local current height $h(x,t)$

The shape of the gravity current can be displayed and easily compared for different $KC_b$ by finding the vertical position of the interface as a function of the horizontal coordinate. Consider first the situation of a density current without lifting. The vertical position of the interface can be estimated by defining a local current height for which different methods have been proposed in earlier studies by, for example, Shin, Dalziel & Linden (Reference Shin, Dalziel and Linden2004) and Anjum et al. (Reference Anjum, McElwaine and Caulfield2013). As suggested by Shin et al. (Reference Shin, Dalziel and Linden2004), the dimensionless local current height $h_{s}(x,t)$ is defined by the vertical integration of the density,

(4.1)\begin{equation} h_s (x,t) = \int_{0}^{1} \rho(x,y,t) \,{{\rm d} y}. \end{equation}

This method is not suitable in the presence of significant vertical density redistributions as observed by Anjum et al. (Reference Anjum, McElwaine and Caulfield2013). They proposed to define a spatially and temporally varying depth of the gravity current, denoted by $h_a(x,t)$, using the available potential energy contained in the current. They derived $h_a(x,t)$ by first introducing the dimensionless centre of mass position $h_{cm}(x,t)$ of the heavy current in the vertical direction, which is defined as

(4.2)\begin{equation} h_{cm}(x,t) = \frac{\int_{0}^{1} \rho (x,y,t) y \, {{\rm d} y}}{\int_{0}^{1} \rho (x,y,t) \, {{\rm d} y}}. \end{equation}

This definition implicitly assumes that $h_{cm}(x,t)$ is not defined right of the front position $X_{fr}$ as then $h_s(x,t) = 0$. We determined the front position by approaching the gravity current from the right-hand side and finding the position $x=X_{fr}$ of the first occurrence where $h_s(x,t)=0.005$ (varying this threshold value to determine $X_{fr}$ between 0.01 and 0.0025 affects the results by less than $1\,\%$). This implies that the value of $h_{cm}(x,t)$ is always clipped to zero for $x\ge X_{fr}$.

The local current height $h_{a}(x,t)$ is then defined as twice the centre of mass height, $h_{a}(x,t)=2h_{cm}(x,t)$. However, gravity currents in the presence of an imposed oscillating flow field experience lifting, and therefore, $h_{a}(x,t)$ overestimates the local current height for the cases with lifting. To include this effect, we propose to quantify the current height $h(x,t)$ as the sum of the (vertical) centre of mass position $h_{cm}(x,t)$ of the heavy current and the local current height with respect to $h_{cm}(x,t)$ (which is not affected by the lifting height):

(4.3)\begin{equation} h (x,t) = h_{cm}(x,t) + 2 \frac{\int_{h_{cm}(x,t)}^{1} \rho (x,y,t) (y-h_{cm}(x,t)) \, {{\rm d} y}}{\int_{h_{cm}(x,t)}^{1} \rho (x,y,t) \, {{\rm d} y}}. \end{equation}

With our definition of $X_{fr}$ and the constraint on $h_{cm}(x,t)$ introduced above, we immediately see that $h(x,t)=0$ for $x\ge X_{fr}$.

In figure 8, we provide a comparison of the local current height according to the methods introduced by Shin et al. (Reference Shin, Dalziel and Linden2004) and Anjum et al. (Reference Anjum, McElwaine and Caulfield2013), given by $h_{s}(x,t)$ and $h_{a}(x,t)$, respectively, with the method introduced in this work for the local current height $h(x,t)$ at the front of a gravity current. It clearly shows the suitability of the modified measure introduced here (see the red line in figure 8).

The shape and evolution of gravity currents for a range of $KC_b$ are explored by visualizing the local current height $h(x,t)$, representing the interface between heavy and light fluid, and computed according to (4.3), at $t=100$, 200, 300 and 400 (figure 9). Since we are not interested in the local details and features of the interface, such as instantaneous variations due to the presence of small vortices at the interface between the heavy and light fluid, we use a moving average of the data from 10 neighbouring grid cells to find the average current height. From figure 9, it is clearly seen that, for $KC_b = 50$ and 100, the interface maintains an inclined shape at the end of each cycle. This is not the case for smaller $KC_b$. In particular, for the cases with $KC_b=5$ and 10, there is a steep increase in the current height at the front of the gravity current. This is similar to the characteristics of the front shape seen for the freely evolving gravity current.

4.5. Summary

The density field $\rho (x,y,t)$ and current height $h(x,t)$ for different simulations show that the gravity current is lifted in the first half of the oscillation cycle ($0^{\circ }\le \phi \le 180^{\circ }$). During this part of the oscillation cycle, the externally imposed flow field aligns with the propagation direction of the gravity current. In the second half of the oscillation cycle ($180^{\circ }\le \phi \le 360^{\circ }$), the unstable situation of heavy fluid being above light fluid causes enhanced (vertical) density redistributions within the gravity current. Enhanced density redistributions already start when the coinciding imposed flow starts to decelerate (for $\phi \ge 90^{\circ }$), and it is further strengthened when the imposed flow is opposite to the propagation direction of the gravity current.

From the comparison of both the density field of the gravity currents and the current heights for different $KC_b$ (and for the freely evolving gravity current), three regimes, each with different behaviour of (small-scale) density redistributions, are observed. The front of the gravity current exhibits instabilities and additional redistribution of the density in the presence of the externally imposed oscillating flow field with $KC_b \lesssim 10$. The gravity current is lifted in the first half of the oscillation cycle and remains lifted throughout the whole cycle since the oscillation period is too short compared with the required time for full reattachment of the heavy current to the bottom of the channel. As a consequence, the subsequent phase with coinciding imposed external flow promotes lifting of the current again before reattachment occurs. Since the duration of the mean shear caused by the imposed flow is shorter, the shape of the current front is similar to the front shape for a freely evolving gravity current. Gravity currents with $KC_b \gtrsim 50$ show density redistributions similar in shape as those observed for salt wedges at the end of the ebb phase in estuarine flow configurations. Such salt wedges also occur, for example, in the Rotterdam Waterway with $KC_b \simeq 5000$. The front of the gravity current is well homogenized and has an inclined shape for $KC_b \gtrsim 50$. Lastly, the cases with $10 \lesssim KC_b \lesssim 50$ show an intermediate behaviour, where the current front has a more inclined shape compared with a smaller period of oscillations. However, the heavy current does not have enough time to become horizontally well homogenized.

5.1. Impact of $KC_b$ on the lifting area and front shape

Observations from the evolving density field of several gravity currents, discussed in the previous section, uncovered a new phenomenon: the lifting of the gravity current. It happens primarily during the first half of the oscillation cycle ($0^{\circ } \le \phi \le 180^{\circ }$). The second half of the oscillation cycle, for $180^{\circ } \le \phi \le 360^{\circ }$, is inherently unstable, since the heavy fluid has moved over light fluid in combination with the externally induced opposing flow, and lifting is no longer supported. In this section, we introduce a measure for and quantify the lifting area $A_{l}$ for gravity currents as a function of $KC_b$, thus as a function of the oscillation period of the external forcing. Additionally, we quantify the redistribution of density as well.

We define the dimensionless lifting height $h_l(x,t)$ in a similar way as we determined the current height $h(x,t)$ in (4.3). More precisely, we quantify the lifting height as the difference of the (vertical) centre of mass position $h_{cm}(x,t)$ of the heavy current and the local current depth with respect to $h_{cm}(x,t)$:

(5.1)\begin{equation} h_l (x,t) = h_{cm}(x,t) - 2 \frac{\int_{0}^{h_{cm}(x,t)} \rho (x,y,t) (h_{cm}(x,t)-y) \, {{\rm d} y}}{\int_{0}^{h_{cm}(x,t)} \rho (x,y,t) \, {{\rm d} y}}. \end{equation}

We use the lifting height $h_l(x,t)$ to introduce a measure for the lifting area $A_l(t)$ beneath the gravity current as a function of time (note that $h_l(x,t)=0$ for $x>X_{fr}$), so that

(5.2)\begin{equation} A_l(t) = \int_{0}^{X_{fr}} h_l(x,t) \, {{\rm d}\kern0.06em x}. \end{equation}

The start of additional density transport due to the decelerating phase of the oscillatory forcing can typically be observed after reaching a maximum lifting area $A_{l,max}$. We found that the lifting area becomes maximum at $\phi \approx 200^{\circ }$, $220^{\circ }$, $240^{\circ }$, $240^{\circ }$, $250^{\circ }$ for $KC_b = 100$, 50, 25, 10 and 5, respectively, and this occurs during each oscillation cycle. Therefore, we determine the maximum lifting area observed during each cycle for the full simulation and average $A_{l,max}$ over all cycles, and denote it with $\langle A_{l,max}\rangle$. With the simulation taking 400 time units, we average over $400/KC_b$ occurrences. A similar procedure is applied for the minimum lifting area $A_{l,min}$, since it is an indicator of whether the unstable configuration (heavy fluid over light fluid) continues to exist during the full cycle (which would imply $\langle A_{l,min}\rangle >0$) or the frontal region reaches a stable and vertically homogenized state within an oscillation cycle. The latter condition would mean that $\langle A_{l,min}\rangle \approx 0$.

The (averaged) maximum lifting area $\langle A_{l,max}\rangle$, shown in figure 10(a), clearly increases with $KC_b$. For the freely evolving gravity current, the averaged maximum lifting area is very small, but not zero, due to a thin boundary layer above the no-slip wall. This very thin boundary layer contains light fluid left between the propagating gravity current and the no-slip wall, and is not related to the differential-advection mechanism discussed for the cases with $KC_b\gtrsim 5$. The data for the (averaged) minimum lifting area show that $\langle A_{l,min}\rangle >0$ for $KC_b=5, 10$ and $25$. Thus, a part of the lifting occurring during the first half of the forcing cycle survives during the second half of the forcing cycle. This is distinctively different for $KC_b=50$ and $100$, where the heavy current is already vertically homogenized before the next oscillation cycle starts.

Figure 10.

(a) Average maximum lifting area $\langle A_{l,max}\rangle$ (red symbols) and average minimum lifting area $\langle A_{l,min}\rangle$ (black symbols) as a function of $KC_b$ (and $Fr=1$). The freely evolving case is plotted on the vertical axis. (b) Scaling of the average maximum lifting area $\langle A_{l,max}\rangle$; the black dashed line indicates a scaling $\langle A_{l,max}\rangle \propto KC_b^{1.1\pm 0.1}$. The fit is based on nonlinear least squares from the data $KC_b\in$(5, 10, 25, 50, 100).

We take a closer quantitative look at the averaged maximum lifting area. Therefore, we replot the data for $\langle A_{l,max}\rangle$ versus $KC_b$ logarithmically in figure 10(b), and the following approximate scaling is found: $\langle A_{l,max}\rangle \simeq 1.62\times 10^{-2} KC_b^{1.1}$. It is, thus, close to linear, for the range $5\le KC_b\le 100$. This scaling behaviour is consistent with the assumption that the lifting phenomenon is due to differential advection with, in this case, $Fr$ constant (thus, the same amplitude $U_0$ of the externally imposed oscillating flow in all cases). During one half of the oscillation cycle, with the externally imposed flow aligned with the propagation direction of the gravity current, the displacement of the heavy fluid in the front of the gravity current due to the externally imposed flow scales with $KC_b$. In contrast, the light fluid near the bottom wall only moves marginally due to the no-slip condition. An approximation for the lifting area would then take the form: $\langle A_{l,max} \rangle \propto \tfrac {1}{2}U_0 KC_b$ (or $\propto Fr KC_b$), consistent with our observation.

To quantify density redistributions in the presence of lifting and squishing phenomena, we focus on the gravity current front. It can be deduced from figure 9 that within a certain distance $\Delta L$ (for example, $\Delta L = 15$) from the front of the gravity current, the current height is most dynamic and irregular. This allows us to illustrate and quantify the different density redistribution behaviours from this part of the gravity current. We define the extent of the gravity current front as the distance between the front position $X_{fr}$, defined in § 4.4, and $X_b=X_{fr} - \Delta L$. The precise choice of $\Delta L$ will be discussed later.

To illustrate the differences in density redistributions with varying $KC_b$, we look at the spatially averaged current thickness. This quantity serves as a geometrical indicator of how much entrainment and redistribution of the density field happens between the heavy and light fluid in the gravity current. For this purpose, we take the current area within the front ($X_b \le x \le X_{fr}$) and normalize it with $\Delta L = X_{fr}-X_b$ to find the average thickness of the current defined as

(5.3)\begin{equation} h_{t} (t) = \frac{1}{\Delta L}\int_{X_b}^{X_{fr}} [h(x,t) - h_l (x,t)] \,{{\rm d}\kern0.06em x}. \end{equation}

A further indicator quantifying density redistributions is the average density in the gravity current front, denoted by $\rho _a(t)$ and defined as

(5.4)\begin{equation} \rho_{a} (t) = \frac{1}{h_t(t)\Delta L}\int_{X_b}^{X_{fr}} \biggl(\int_{h_l(x,t)}^{h(x,t)}\rho (x,y,t) \,{{\rm d} y}\biggr) \,{{\rm d}\kern0.06em x}. \end{equation}

Both $h_t(t)$ and $\rho _a(t)$ are evaluated with a predetermined length of the front region, $\Delta L = 15$. The choice of $\Delta L = 15$ might seem subjective at first glance, but we will explain this choice below.

In figure 11(a), we show $h_t(t)$ for the freely evolving and the five different oscillatory-forced cases with $\Delta L = 15$. For all simulations, we show the results for $t=50, 100,\ldots, 400$. Similarly, we show $\rho _a(t)$ in figure 11(b) for the same times. The freely evolving case is shown by the black filled squares. The average thickness of the gravity current front decreases during the first part of the propagation phase (up to $t\approx 200$) before it tends to saturate at a constant value. The average density $\rho _a(t)$ in the gravity current front shows that, for increasing $KC_b$, the entrainment and redistribution of the density field in the gravity current front is enhanced compared with the freely evolving case, and we see that increasing $KC_b$ implies smaller values for $\rho _{a}(t)$. The curves shown in figure 11(b) suggest an approximately exponential decay of the average density, $\rho _a(t)\approx \exp (-t/\tau _{\rho })$, with $\tau _{\rho }$ a density relaxation time scale which depends on $KC_b$. We have computed $\ln \rho _a(t)$ for the range of $\Delta L$ between $5$ and $20$, and found that $\tau _{\rho }$ becomes independent of the choice of $\Delta L$ for $\Delta L\gtrsim 15$. This observation motivates our choice $\Delta L=15$ for this quantitative analysis. The inset in figure 11(b) shows $\tau _{\rho }$ for the range of $KC_b$ explored in this study. It is clear from the monotonic decrease of $\tau _{\rho }$ with $KC_b$ in the range $0\le KC_b\le 100$ that density homogenization occurs faster with increasing $KC_b$.

Figure 11.

(a) Average current thickness $h_t(t)$ of the gravity current front and (b) average density $\rho _a(t)$ in the gravity current front for the freely evolving gravity current and the five cases with an externally applied oscillating pressure gradient. For all cases, $Fr=1$. The horizontal extent of the gravity current front is taken as $\Delta L=15$. The inset in panel (a) shows the current thickness averaged over the times $200, 250, \ldots, 400$. Estimated error margins are included. The inset in panel (b) shows the density relaxation time $\tau _\rho$ for the same cases. For both insets, the values for the freely evolving case are on the vertical axis.

These conclusions can be further substantiated by considering the probability distribution of the density, $P(\rho )$, in the same area (defined by $X_b\le x\le X_{fr}$ and $h_l(x,t)\le y\le h(x,t)$) used to compute $\rho _a(t)$. The density probability distribution is normalized using $\rho _a(t)$ such that $\int _0^1 P(\rho )\,{\rm d}\rho =1$ for each set of $KC_b$, $Fr$ and $Sc$ values. In figure 12, we show $P(\rho )$ at $t=100$ (figure 12a) and $t=200$ (figure 12b) for the freely evolving and oscillatory-forced gravity currents (with $Fr=1$ and $Sc=5$). The characteristics of the density redistributions for $KC_b=5$ and 10 are similar to those of the freely evolving gravity current, and enhanced density homogenization is observed for $KC_b=50$ and 100. Initially, the characteristics of the density redistributions for the case with $KC_b=25$ is close to the freely evolving case, but at later times, it shows similar homogenization of the density field as observed for the higher $KC_b$ runs. For $t\gtrsim 300$, the density probability distribution tends to a continuous uniform distribution given by $P(\rho )=1/(\rho _{max}-\rho _{min})$ for $\rho _{min}\le \rho \le \rho _{max}$ and zero elsewhere. Typically, $\rho _{min}\approx 0.1$ and $\rho _{max}\approx 0.5\unicode{x2013}0.7$ (the upper bound for small $KC_b$ and lower bound for large $KC_b$). In this interval, density values are more or less equally probable.

Figure 12.

Density probability distribution $P(\rho )$ as a function of $KC_b$ (and $Fr=1$, $Sc=5$) taken at (a) $t= 100$ and (b) $t=200$.

The data for $h_t(t)$ support the observation from the previous section that we can distinguish three regimes. For $KC_b\lesssim 10$, the gravity current front remains always thicker than for the freely evolving case, while it is typically thinner for $KC_b\gtrsim 50$. The case $KC_b=25$ represents a transition regime with the thickness of the gravity current front being similar to that in the freely evolving gravity current. These observations suggest that, for $KC_b\lesssim 10$, entrainment and density redistributions mostly take place near the gravity current front, resulting predominantly in growth of its thickness. For larger $KC_b$, there is a strong interplay between entrainment and (vertical) density redistribution in the gravity current front, which promotes enhanced thickness, with strong horizontal advection. The latter mechanism removes parts of the top layer of the gravity current front to the rear, effectively resulting in a thinner but nevertheless lighter gravity current front.

We now consider average values of the current thickness, $\langle h_t\rangle =\tfrac {1}{5}\sum _{n=4}^8 h(t=50n)$, thus using only the values outside the transient regime occurring for $t \le 150$. In the inset of figure 11(a), we show $\langle h_t\rangle$ (with an estimated error margin) obtained for $\Delta L=15$. The results support the conclusion above that, for $KC_b \le 10$, the average thickness of the front is larger compared with the freely evolving case and the other way around for $KC_b\ge 50$.

5.2. Effect of velocity amplitude (or $Fr$) of the externally imposed oscillating flow

The overall behaviour of gravity currents with $KC_b \gtrsim 50$ show similarities with those of salt intrusions in estuaries. Therefore, we explored the effect of $Fr$ for simulations with $KC_b = 50$. It allows us to consider a substantial number of oscillation cycles (eight in total) to distinguish between transient effects (typically during the first three cycles) and quasi-steady oscillatory propagation of the gravity current during the last cycles. We consider $Fr = 0.25$, 0.50, 1 and 2 and the freely evolving case ($Fr=0$). This comparison is thus based on the simulations 1, 5, 7, 8 and 9 in table 2. We will incorporate results from simulations where the initial phase of forcing is $\phi _{init} = 0^{\circ }$ and $180^{\circ }$. An animation of the density field evolution for $50\le t\le 200$ is shown in supplementary movie 3.

In §§ 4 and 5.1, we have discussed the basic mechanisms that affect the lifting height, its temporal evolution, the maximum lifting area $A_{l,max}$ and the role of $KC_b$ (with $Fr$ kept fixed) on these quantities. In figure 13, we show $A_{l,max}$ at the end of each oscillation cycle for $Fr=0$, 0.25, 0.5, 1 and 2. From these data, it is clear that the first three cycles represent a transient period before the lifting area reaches an almost steady state. From the steady-state range, we get confirmation of our earlier hypothesized scaling of the maximum lifting area, $A_{l,max}\propto Fr KC_b$. This observation holds for both initialization phases, $\phi _{init}=0^{\circ }$ and $\phi _{init}=180^{\circ }$.

Figure 13.

(a) Maximum lifting area $A_{l,max}$ for each cycle of the oscillatory forcing of the different gravity currents under consideration ($KC_b=50$). Two initializations are considered: $\phi _{init}=0^{\circ }$ (solid line) and $\phi _{init}=180^{\circ }$ (dotted line). Inset shows the vertical profile of the horizontal velocity of the externally imposed flow field at $\phi = 90^{\circ }$. The results from the steady-state range support our hypothesis $A_{l,max} \propto Fr$ and hardly depend on $\phi _{init}$.

To understand the effect of $Fr$ on the distance travelled by the front of the gravity current, we investigate the front position $X_{fr}$, as defined in § 4.4, for different $Fr$ and compare it with the freely evolving case ($Fr=0$). The front positions $X_{fr}$ as a function of time are shown in figure 14(a). To distinguish the effect of externally imposed oscillatory forcing on the front position, we determine the front position relative to the freely evolving case (with front position $X_{fr,0}$), denoted by $\Delta X_{fr}= X_{fr}-X_{fr,0}$ (figure 14b). The first three oscillation cycles are affected by transient behaviour and are, hence, disregarded in further discussions. To quantify the effect of the oscillatory forcing on $X_{fr}$, we consider the following cycle-averaged quantity:

(5.5)\begin{equation} \langle \Delta X_{fr}\rangle_n = \sqrt{\frac{1}{KC_b}\int_{nKC_b}^{(n+1)KC_b}\Delta X_{fr}^2\,{\rm d}t}, \end{equation}

with $3 \le n \le 7$. It represents the averaged effect of the oscillatory forcing on $\Delta X_{fr}$ in each individual forcing cycle. The results for $\langle \Delta X_{fr}\rangle _{n}$ are shown in figure 14(c). We have computed the average $\langle \Delta X_{fr}\rangle =\tfrac {1}{2}(\langle \Delta X_{fr}\rangle _6+\langle \Delta X_{fr}\rangle _7)$ for the case with $\phi _{init}=0^{\circ }$. The results show that $\langle \Delta X_{fr}\rangle \propto Fr$ and advances the freely evolving case (see the solid curves in figure 14b,c). Furthermore, if we focus on simulations with varying $KC_b$ and keep $Fr=1$, we observe that $\langle \Delta X_{fr}\rangle \propto KC_b$. Carrying out the same approach for simulations with $\phi _{init} = 180^{\circ }$ shows again that $\langle \Delta X_{fr}\rangle \propto Fr$, but lags behind the freely evolving case (see the dotted curves in figure 14b,c). Our results provide strong indications that $\langle \Delta X_{fr}\rangle \propto Fr KC_b$, which makes sense intuitively, because $\langle \Delta X_{fr}\rangle$ will scale with the maximum velocity of the oscillation and the period of the oscillation. This suggests that differential advection is the key process determining the relative front position.

Figure 14.

Front position of the gravity current for simulations with varying $Fr$ with (a) $X_{fr}$, (b) front position for simulations with varying $Fr$ relative to the non-oscillating case, $\Delta X_{fr} = X_{fr}-X_{fr,0}$, and (c) period-averaged propagation speed $\langle \Delta X_{fr}\rangle _{n}$. The solid and dotted curves represent results with $\phi _{init}=0^{\circ }$ and $\phi _{init}=180^{\circ }$, respectively. The colour code indicates $Fr$.

The frontal region of the gravity current becomes thinner and more diluted in the presence of oscillatory forcing with large $KC_b$; see, for example, figures 7 and 11. To quantify the change in density within the current, we first introduce the total mass of the heavy fluid, $M_{Fr} (t) = h_t(t)\Delta L \rho _a(t)$, in the frontal region of the gravity current. Note that the dimensionless density of the light fluid is $\rho =0$ and of the heavy fluid $\rho =1$, and $0 \le \rho _a(t) \le 1$. The frontal region is defined by $X_b \le x \le X_{fr}$, with $X_b=X_{fr} - \Delta L$ (with $\Delta L = 15$), and $M_{Fr} (t)$ can be determined for each $Fr$. The relative mass difference in the frontal region compared with the freely evolving case is defined as $\Delta M_{Fr}(t)={[M_{Fr}(t)-M_0(t)]}/{M_{0}(t)}$, with $M_{0}(t)$ the total mass in the frontal region for the freely evolving gravity current. For $Fr\ge 0.5$, we find that the period-averaged value of $\Delta M_{Fr}(t)$ is negative, and for $Fr=0.25$, it is approximately zero. We start by considering $\phi _{init}= 0^{\circ }$. The period-averaged relative mass difference $\langle \Delta M_{Fr}\rangle _n$ in the frontal region is defined as

(5.6)\begin{equation} \langle \Delta M_{Fr}\rangle_n = \sqrt{\frac{1}{KC_b}\int_{nKC_b}^{(n+1)KC_b}\Delta M_{Fr}^2\,{\rm d}t}, \end{equation}

and is by definition positive. It increases significantly with increasing $Fr$, see figure 15(a), and also shows transient behaviour for $t \le 150$. Therefore, we focus on the average of $\langle \Delta M_{Fr}\rangle _n$ for $200 \le t \le 400$, denoted by $\langle \Delta M_{Fr}\rangle = \tfrac {1}{4}\sum _{n=4}^7\langle \Delta M_{Fr}\rangle _n$. Due to the fact that the period-averaged value of $\Delta M_{Fr}(t)$ is typically negative, there is no ambiguity in the interpretation of $\langle \Delta M_{Fr}\rangle$ (which is always positive). Figure 15(a) shows an increase of $\langle \Delta M_{Fr}$ with $Fr$. For $\phi _{init}=180^{\circ }$, very similar behaviour is found for $\langle \Delta M_{Fr}\rangle _n$ (see figure 15a).

Figure 15.

(a) Relative total mass difference $\langle \Delta M_{Fr}\rangle _n$ for different $Fr$ (and $KC_b=50$) in the gravity current front averaged for each cycle compared with the freely evolving case. (b) Relative total mass through the gate $\langle \Delta M_{l,Fr}\rangle _n$ averaged for each oscillation cycle for the four values of $Fr$ relative to the freely evolving case. Solid lines, $\phi _{init}=0^{\circ }$; and dotted lines, $\phi _{init}=180^{\circ }$. The frontal region is defined by $X_b \le x \le X_{fr}$, with $X_b=X_{fr} - \Delta L$ (with $\Delta L = 15$).

Also here, we can consider the probability distribution of the density, $P(\rho )$. In figure 16, we show $P(\rho )$ at $t=100$ (figure 16a) and $t=200$ (figure 16b) for the freely evolving and oscillatory-forced gravity currents with $KC_b=50$ and $Fr=0$, 0.25, 0.5, 1 and 2 (with $Sc=5$). We observe a gradual increase of the redistribution of the density with $Fr$ with a strong impact for $Fr\gtrsim 1$. This behaviour can also be observed qualitatively with animations of the density fields from simulations with varying $Fr$ (for $50\le t\le 200$) shown in supplementary movie 3.

Figure 16.

Density probability distribution $P(\rho )$ as a function of $Fr$ (and $KC_b=50$, $Sc=5$) taken at (a) $t=100$ and (b) $t=200$.

The decrease of the total mass in the frontal region supports our qualitative observations that the upper part of the gravity current front is pushed further backwards to its rear by the strong opposing flow (for $KC_b=50$). The net advection of heavy fluid towards the frontal region is thus, on average, less in the presence of oscillatory forcing, which reduces the gravity current thickness and its average density. However, its propagation properties are dominated by differential advection, as discussed in the previous paragraph.

To further understand the bulk density transport, we compare the mass of heavy fluid passing through the gate (at $x=0$) for gravity currents with $KC_b=50$ and $Fr=0.25$, 0.5, 1 and 2, and the freely evolving case ($Fr=0$). To find the total mass of heavy fluid that passed the gate position, we introduce $M_{l,Fr}(t)$, defined as

(5.7)\begin{equation} M_{l,Fr} (t) = \int_{0}^{({1}/{4})L}\int_0^1 \rho (x,y,t) \,{{\rm d}\kern0.06em x} \,{{\rm d} y}. \end{equation}

In a similar way as before, we introduce $\Delta M_{l,Fr}(t)= {[M_{l,Fr}(t)-M_{l,0}(t)]}/{M_{l,0}(t)}$, which also depends on $\phi _{init}$. For $\phi _{init}=0^{\circ }$, the period-averaged value of $\Delta M_{l,Fr}$ is positive, while it is negative for $\phi _{init}=180^{\circ }$ due to asymmetry of the forcing. The averaged total mass of heavy fluid passing through the gate is defined in a similar way as $\langle \Delta M_{Fr}\rangle$ and results shown in figure 15(b) suggest $\langle \Delta M_{l,Fr} \rangle \propto Fr$. For $\phi _{init}=0^{\circ }$, the total mass transport through the gate is always larger compared with the freely evolving case for $0.25\le Fr\le 2$. For $\phi _{init}=180^{\circ }$, the total mass transport through the gate is less for $0.25\le Fr\le 2$ compared with the freely evolving case (note that, according to our definition, $\langle \Delta M_{l,Fr} \rangle >0$). For both cases, the data show that $\langle \Delta M_{l,Fr}\rangle \propto Fr$, thus is connected with the excursion length of the oscillations and differential advection.

5.3. Effect of the Schmidt number on lifting area and propagation speed

The rate of molecular diffusivity changes according to the scalar that causes the density difference. Additionally, the effective Prandtl–Schmidt number also changes for different applications and conditions in estuaries with different stratification. Therefore, we are interested to see the effect of diffusivity on mass transport (of heavy fluid) and gravity current dynamics. For this purpose, we carried out simulations with $Sc=1$, 5 and 20, while we kept the other parameters fixed: $KC_b=50$, $Fr=1$ and $Re=3000$.

To show the effects of varying $Sc$ on the front position of the gravity current, we only focus on oscillatory-forced gravity currents with $\phi _{init}=0^{\circ }$. The cycle-averaged front position of the gravity current,

(5.8)\begin{equation} \langle \Delta X_{fr}\rangle_n = \sqrt{\frac{1}{KC_b}\int_{nKC_b}^{(n+1)KC_b} {[X_{fr,1} - X_{fr,0}]^2\,{\rm d}t}}, \end{equation}

represents an averaged propagation difference between the gravity current in the presence of oscillatory forcing (with front position $X_{fr,1}$) compared with the freely evolving gravity current (with front position $X_{fr,0}$) for varying $Sc$. The average $\langle \Delta X_{fr}\rangle$ over the last two forcing cycles ($300 \le t \le 400$) shows that the cycle-averaged front position of the gravity current advances the freely evolving case by a similar amount for the three values of $Sc$, see figure 17(a), and hardly influences the relative propagation of gravity currents. The maximum lifting area $A_{l,max}$ of the gravity currents, using the definition from (5.2), shows a similar trend for different $Sc$ (see figure 17b). They show a transient behaviour, where $A_{l,max}$ increases for $n \le 3$, before saturating. Since the molecular diffusion for $Sc=1$ is higher, the gravity current becomes more diluted at the front while differential advection occurs. Therefore, we expect that $\langle A_{l,max} \rangle$ will be somewhat smaller for $Sc=1$ compared with $Sc=5$ and 20.

Figure 17.

(a) Period-averaged front position of the oscillatory-forced gravity current with respect to propagation of the freely evolving gravity current, $\langle \Delta X_{fr}\rangle _{n}$, (b) maximum lifting area of gravity current, $A_{l,max}$, for simulations with different $Sc$ and (c) relative total mass difference $\langle \Delta M_{Sc} \rangle _n$ for different $Sc$ at the gravity current front for each cycle compared with the freely evolving case. The initial phase of the forcing is $\phi _{init} = 0^{\circ }$ and the frontal region is defined by $X_b = X_{fr} - \Delta L \le x \le X_{fr}$ (with $\Delta L = 15$). The colour indicates $Sc$.

The relative mass of heavy fluid in the frontal region is defined as $\Delta M_{Sc}(t)={[M_{Sc}(t)-M_{Sc,0}(t)]}/{M_{Sc,0}(t)}$, with $M_{Sc,0}(t)$ the total mass in the frontal region of the freely evolving gravity current for the Schmidt number under consideration. The cycle-averaged relative mass of heavy fluid in the frontal region for varying $Sc$ is expressed as

(5.9)\begin{equation} \langle \Delta M_{Sc}\rangle_n = \sqrt{\frac{1}{KC_b}\int_{nKC_b}^{(n+1)KC_b}\Delta M_{Sc}^2\,{\rm d}t}. \end{equation}

Also for the cases evaluated here, there is a reduction in the amount of mass in the frontal region, so there is no ambiguity in the meaning of $\langle \Delta M_{Sc}\rangle _n$. The cycle-averaged relative mass $\langle \Delta M_{Sc}\rangle _n$ shows a significant decrease for $Sc=1$, see figure 17(c), while those for $Sc=5$ and 20 saturate. This is also reflected in the probability distribution of the density $P(\rho )$ for different $Sc$, not shown. The averaged value of $\langle \Delta M_{Sc}\rangle _n$ for $300 \le t \le 400$ shows that the relative mass at the front region for each $Sc$ is significantly reduced, particularly for $Sc=5$ and 20. This indicates that the relative mass difference at the frontal region between the oscillating and freely evolving case is larger for higher $Sc$ values since the freely evolving gravity current is less diluted at the front for those cases.