1. Introduction
Gravity currents are flows that develop when there is a horizontal density gradient in the presence of a gravitational field. This density difference causes the heavier fluid (gravity current) to intrude beneath the lighter one (ambient fluid) with a mainly horizontal displacement (Huppert and Simpson, 1980; Linden, 2011). Gravity currents appear in a variety of natural and man-made phenomena at different scales (Simpson, 1999; Ungarish, 2009; Wu and Dai, 2020). Examples include thermally driven winds in the hillside of the mountains, sea breezes, oceanic fronts, squall lines in thunderstorms, pyroclastic flows, and the cold air from outside entering a warm house when the door is opened.
When a fluid is in a gravitational field, buoyancy forces appear whenever there are density variations. Particularly, when density varies in the horizontal direction, a gravity current (GC) will always result (Linden, 2012). The GC acquires an acceleration called reduced gravity, g', which depends on the density contrast between the GC (with density ρc) and the ambient fluid around it (with density ρa),
where
where u̅ is the mean velocity and h is the height of the GC. Additionally, the Reynolds number,
where L is a characteristic length and v is the kinematic viscosity, a key parameter for studying turbulent flows. Finally, when investigating the mixing between the gravity currents and the ambient fluid, the Richardson number is crucial. This dimensionless parameter expresses the ratio of the buoyancy term to the flow shear term, which can be estimated as:
where u is the velocity of the flow and ρ its density. The Richardson number in gravity currents is further analyzed in this study because the entrainment, E, depends on this dimensionless parameter (Ellison and Turner, 1959). Because they influence the dynamics of the atmosphere and the ocean and even have industrial applications (Fleischmann et al., 1994), gravity currents have been the subject of numerous studies. Many of them have focused on experimental gravity currents propagating over smooth flat bottoms or smooth slopes using the typical lock-exchange configuration (Mitsudera and Baines, 1992; Simpson, 1999; Reuten et al., 2010; Linden, 2012; Shin et al., 2004; Adduce et al., 2012; Lowe et al., 2002; Ungarish, 2016; Balasubramanian and Zhong, 2018; Mukherjee and Balasubramanian, 2020, 2021). Nevertheless, in nature, gravity currents propagate over complex topographies. In particular, in this study, we are interested in gravity currents that may simulate the thermally driven winds propagating over a topographic slope.
There has been a growing interest in understanding the impact of bottom roughness (Nogueira et al., 2013, 2014; Cenedese et al., 2016, 2018; Tanino et al., 2005; Zhang and Nepf, 2008; Bhaganagar and Pillalamarri, 2017; Wilson et al., 2017; Jiang and Liu, 2018; Xu et al., 2019; Zhou and Venayagamoorthy, 2020; Maggi et al., 2022) and isolated obstacles (Tokyay and Constantinescu. 2015, de Falco et al. 2021) on gravity current propagation. In particular, the presence of rectangular obstacles was shown to increase entrainment as the GC propagates downstream (Wilson et al., 2017). On the other hand, different bed roughnesses were tested by Nogueira et al. (2013), and they found that diverse granular media play an important role in the current dynamics, mainly because the extra drag at the bed decreases the front velocity (Nogueira et al., 2013). Bottom roughness has been experimentally idealized by different arrangements of rigid cylinders by Cenedese et al. (2016, 2018) and Maggi et al. (2022). The former measured density fields and found entrainment enhancement by different processes depending on the configuration of the cylinders. In contrast, the latter found that the relative height of the cylinders to the gravity current depth controls the structure and propagation speed of the current with higher roughness elements generating strong recirculation areas between cylinders that interact and modify the overlying current. These investigations have shown that idealized topography affects the gravity current profile (Tanino et al., 2005) and the mixing mechanisms (Cenedese et al., 2018). To quantify entrainment, previous authors (Cenedese and Adduce., 2008; Nogueira et al., 2013; Ottolenghi et al., 2016; Wilson et al., 2017; Ottolenghi et al., 2017; Balasubramanian and Zhong, 2018) have utilized the Morton-Taylor-Turner (MTT) entrainment hypothesis (Ellison and Turner, 1959):
where U is a characteristic velocity and w e is the entrainment velocity, perpendicular to U. As the GC propagates, the velocity of inflow into the turbulent region can be assumed proportional to the velocity scale of the turbulent layer, and the constant of proportionality is E. Ellison and Turner (1959) estimated E for neutral and traveling downslope jets; they showed that this coefficient ranges between 0.01 and 0.1. Other studies have considered these values for studying gravity currents (Turner, 1986; Princevac et al., 2008; Manins and Sawford, 1979). However, the MTT hypothesis may underestimate the entrainment rate in atmospheric flows (Princevac et al., 2008) because natural flows occur at much higher Reynolds numbers (Re≈107) than the ones in a laboratory tank (Re≈103). Experimental and numerical gravity current studies have extended the range of entrainment values between 10-4 and 10-1 (Ottolenghi et al., 2016; Cenedese and Adduce, 2008).
To estimate entrainment using equation 5, the velocities U and w e must be quantified. While it is relatively simple to calculate U, parametrizing the entrainment velocity w e is not easy. Some authors have evaluated w e as the flow of ambient fluid crossing the interface between the dense current and the ambient fluid around it per unit area. This idea was proposed by Cenedese and Adduce (2008), and it has been widely implemented to quantify entrainment (Nogueira et al., 2013; Ottolenghi et al., 2016, 2017; Wilson et al., 2017).
The present study investigates a synthetic topography’s influence on gravity currents’ dynamics. The entrainment coefficient is estimated from the velocity fields, which were measured experimentally using Particle Image Velocimetry (PIV). Utilizing these measured velocity fields, the entrainment is studied by the MTT assumption (equation 5) and with an analysis presented by Princevac et al. (2008), which will be further explained in the Methodology. We are particularly interested in the dependence of ntrainment on bottom topography. To get insight into this dependence, the experimental gravity currents were allowed to propagate over different sections of a topographic surface as similar as possible to a mountain slope near the Valley of Mexico. Additionally, qualitative visualizations are performed to better characterize the flow propagation and the spatial and temporal variations of the mixing. These visualizations illustrate the mixing while gravity currents propagate over the topographic slope.
The particular situation of gravity currents propagating downslope on a topographic surface is of interest in the Valley of Mexico because thermally-driven winds are known to propagate in that region (Doran et al., 1998; Fast and Zhong, 1998). Nevertheless, identifying these flows in the atmosphere is hardly achieved with the velocity fields because those fields are not as easily measured in the valley as they can be measured in a laboratory tank or a numerical simulation. Therefore, we are also interested in finding a dynamic condition to help us identify gravity currents in the atmosphere. We also utilized the measured velocity fields to estimate instantaneous pressure fields and obtain pressure time series over synthetic stations along the topographic surface that is shown in figure 1. These pressure time series allowed us to relate fast variations in the pressure to the arrival of gravity currents to synthetic stations. The measurement of atmospheric pressure in stations around a valley is much simpler and cheaper than the velocity field measurements. Therefore, such a pressure condition may help detect atmospheric gravity currents.
2. Methodology
2.1 Experimental set-up
Experimental lock-release gravity currents are produced in a tank of rectangular cross-section (see figure 2). The density contrasts between the two liquids are created using saline solutions. The resulting densities were measured with an Anton Paar densitometer. Each GC develops over a flat section (15 cm long) followed by a slope of 5º. A synthetic topography was printed over an acrylic slab following the “Tláloc’’ volcano hillside profile. The manufacturing process was possible using a computer numeric control (CNC) machine.
The synthetic topography on the slope is scaled with the one from the “Tláloc” volcano near the Valley of Mexico.
Because the horizontal scale is much larger than the vertical scale, we split the entire slope into three sections. We consider eight cases described in Table I. For comparison purposes, the base case is the GC propagating over a flat slope. The rest of the seven cases are performed using a synthetic topography. The entire section is divided into three subsections (Figure 1 b). We explore different combinations of parameters, such as g' and the varying topography over the slope. Firstly, we investigate the propagation of gravity currents with initial g' = 14 cm s-2 over each of the topographic and flat slope sections. Secondly, we study gravity currents varying g' and propagating over the lower section of the topographic slope. Two smaller and two larger g' than14 cm s-2 are considered.
Run | Bottom | Γ | g'(cms-2) | U(cms-1) | Re u | Fr u | Ri | E(mean) |
F01 | flat (no topography) | 0.014 | 14 | 6.23 | 18615 | 1.36 | 0.34 | 0.0406 |
TU01 | upper (topography) | 0.014 | 14 | 3.11 | 9292 | 0.68 | 0.25 | 0.0459 |
TI01 | intermediate (topography) | 0.014 | 14 | 6.70 | 20020 | 1.46 | 0.19 | 0.0430 |
TL01 | lower (topography) | 0.014 | 14 | 4.51 | 13476 | 0.98 | 0.38 | 0.0458 |
TL02 | lower (topography) | 0.008 | 8 | 4.25 | 12699 | 1.23 | 0.15 | 0.0457 |
TL03 | lower (topography) | 0.011 | 11 | 4.70 | 14044 | 1.16 | 0.14 | 0.0440 |
TL04 | lower (topography) | 0.018 | 18 | 5.10 | 15239 | 0.98 | 0.21 | 0.0451 |
TL05 | lower (topography) | 0.024 | 24 | 5.73 | 17121 | 0.95 | 0.22 | 0.0446 |
We utilized two different experimental techniques to study the GCs. The first is a simple qualitative visualization with a shadowgraph, and the second uses Particle Image Velocimetry (PIV). The shadowgraph visualization is implemented by displaying a white light source from behind the tank and placing a screen opposite it. A camera then records the pattern of lights and shadows on the screen.
A two-dimensional particle image velocimetry technique or 2D-2C PIV was employed to obtain the velocity fields of the gravity current. Illumination, capturing, and processing equipment were necessary. A Dantec Dynamics System was used to perform the measurements. The illumination system consisted of a doubled pulsed Nd:YAG Litron laser, illuminating the test section from below the tank using special optics that transformed a cylindrical beam into a light sheet. The fluid was previously seeded with neutrally buoyant silver-coated microspheres of glass. A Phantom CMOS Speed Sense M320 camera was used to capture the images and was placed facing perpendicular to the laser sheet. A 40 mm Nikon lens was also employed.
The camera has a resolution of 1920 × 1200 pixels; however, for these experiments, the region of interest was 1858 by 755 pixels. The laser and camera shots were synchronized and controlled by the Dynamic Studio software version 4.0.37 from Danted Dynamics. The velocity vectors were also obtained using this software. An image masking procedure was implemented to isolate the fluid region where vectors were computed. An adaptive cross-correlation algorithm was used to obtain the instantaneous velocity fields. Different criteria and filters were employed to eliminate spurious vectors and improve velocity fields. Typically, 2960 independent vectors were obtained with a sampling frequency of 55 Hz and a spatial resolution of 3.85 mm. The region of interest corresponds approximately to an area of 325.6 × 132.4 mm.
2.2 Entrainment estimations
To use the MTT assumption (equation 5) for the estimation of E, we consider the analysis presented by Manins and Sawford (1979) and Princevac et al. (2008). The Princevac analysis considers a downslope gravity current with a transformation of the coordinate system (s,n) such that s and n are the tangent and normal directions to the slope, respectively. For each gravity current, the height is represented with two values: H and h. The former (H) is the distance from the solid boundary to the interface between the current and its surroundings, taking into account the vortex layer; while h does not take this layer into account. The schematic representation for these variables is shown in figure 3. The vortex layer is the region where the entrainment occurs between the gravity current and its surroundings. The scale of the velocity U, parallel to the slope, is calculated from the PIV estimated velocity fields, considering the equation shown by Princevac et al. (2008):
where u is calculated from the measured velocity fields, it is noteworthy to mention that the velocities are measured inside the gravity current domain at each time step in the direction of the slope. The characteristic velocity for entrainment, w e , is rather difficult to be estimated from the experimental data (Wilson et al., 2017). Here, we developed a MATLAB code to estimate this w e by directly using velocity fields and not from the area considerations proposed by previous authors (Ottolenghi et al., 2016; Nogueira et al., 2013; Wilson et al., 2017). To do so, we calculate the maximum values of vorticity in the s, n coordinates; then, we consider a small region around that maximum for each s and define that range as the vortex layer (H − h). This small region was estimated by the 5 pixels around the maximum vorticity pixel. Then, we look at the velocities inside this region (H − h) and the n component inside that vortex layer, which are considered the entrainment velocities w e . For each instantaneous velocity field measurement, we estimate E using the mean values of U and w e inside the gravity current for that velocity field. Then, for each gravity current, we get an entrainment time series (equation 5).
2.3 Pressure time series
Pressure fields are estimated from the velocity field measurements. An algorithm developed by Dabiri et al. (2014) that integrates the pressure gradient term from the Navier-Stokes equations is implemented for our study. This algorithm is contained in a MATLAB toolbox called queen 2.0 Pressure Calculator (Dabiri et al., 2014). From our instantaneous 2D velocity field estimates, we compute the corresponding instantaneous pressure fields, which are later used to generate the pressure time series at the stations shown in figure 1 (c). Once having those instantaneous pressure fields, we define synthetic stations along the topographic slope and calculate the pressure time series in each station. We obtain a set of six pressure time series for each gravity current propagating along the topographic slope.
3. Results
3.1 Velocity Fields
As the gravity currents propagated over the flat and topographic slopes, velocity fields were obtained with the PIV technique. Figure 4 shows these fields for three of the cases of gravity currents: F01 (flat slope), TU01 (upper topographic slope), and TL01 (lower topographic section). All these three cases correspond to g' = 14 cm s-2; however, the velocity fields of these gravity currents show different structures, elucidating that the topography clearly impacts their propagation. Over the flat slope, higher velocities are reached, and they are uniform regardless of GC propagation. Moreover, the velocity reaches higher values as the GC evolves, showing a more rapid displacement at 5.45 and 7.27 s. Over the upper part of the topographic slope (TU01), as the head of the gravity current displaces over an irregular slope, higher velocities are reached in some zones of the GC compared with other zones; this becomes particularly clear at 5.45 and 7.27 s. Over the lower part of the topographic slope (TL01), non-uniform velocities also develop within the GC, and a clear vortex mixing layer is identified in the inter-phase between the GC and the ambient fluid.
3.2 Entrainment estimations
Entrainment is estimated for each GC as it propagates over the different sections of the slope. We compare the mean entrainment of each run and investigate the cases that resulted in higher entrainment and, therefore, more mixing. Entrainment time series are obtained to investigate how mixing processes occur while the gravity currents propagate. For each velocity field, E is calculated with the mean values of U and w e ; this last quantity is estimated from the velocities in the vortex layer following the procedure explained in section 2.2. We consider a dimensionless time t = t / t f , where t is the current time, and t f is the total time lapse. Figure 5 shows the entrainment E for two cases, with g' = 14 cm s-2. One corresponds to a GC over the flat slope; the other is a selected GC in a topographic section (runs F01 and TL01). The rest of the results are summarized in figure 6.
The entrainment time series of figure 5 shows that E rapidly fluctuates in time and reaches its highest values in the first half of the time series. During this lapse, the head of the gravity current (and the vortex behind this head) are in the domain of study, i.e., in the section of the flow being recorded by the camera. A decreasing tendency appears at the end of the time series when the main vortices that enhance mixing are already out of the domain of study. As the GC flows over the slope, it dissipates energy over the rough topography and, therefore, also loses entrainment potential.
Entrainment values are small (on average 4.59 × 10-2 and 4.15 × 10-2 for TL01 and F01, respectively), and due to the turbulent nature of the vortex region defined for the entrainment calculation, it shows a large variability. Entrainment variability is slightly more considerable for the smooth topography experiment (F01), with variations representing 13.9% of the mean compared to 12.4% for the synthetic topography experiment (TL01). Amplitude spectra for these time series (Figure S3 in the supplementary material) show that energy is distributed among fewer frequencies in F01 than in TL01, which suggests that there is less variability of scales for the smooth case compared to the rough topography case.
Autocorrelation plots (Figure S4) show a weak positive autocorrelation for both time series, with correlation values of 0.2 or less for TL01 and between 0.4 and 0.2 for F01, although correlation values for the smooth case (F01) almost double those for the synthetic topography (TL01) for the first seven lag values (between 0.18 and 1.26 s). Autocorrelation decreases as the lag increases for the smooth case, which reflects the fact that E has a linear trend and declines in time (Fig. 5a). This is not the case for the rough topography GC, for which correlations are immediately smaller and increase noisily to a maximum correlation of 0.3 at lag 12.
We estimate the Richardson number (equation 4) for each instantaneous velocity field. We use fixed values of d ρ / dz corresponding to Table I. A mean Richardson number (Ri) is calculated for each GC run. Figure 6 shows the relation between E and Ri. All the data points are added in figure 6, comparing all the experiments in Table I. A general trend indicates that the Richardson numbers are smaller for the synthetic topography, conversely with the flat cases with higher values. The lower the Ri, the higher the mixing; therefore, it appears that the topographic irregularities on the slope contribute to enhancing mixing.
Figure 6 shows that the lowest values of entrainment correspond to the flat slope (run F01, purple squares in the figure). Although this run has the highest velocity (and therefore the highest Reynolds number), the entrainment remains the lowest (E mean = 0.0406), and the Richardson number in this gravity current reaches high values compared with the Ri corresponding to the GCs over the topographic slope. The intermediate section of the topographic bottom is the flattest one of this topography (see figure 1). In this case, E is also low (E mean = 0.0430) compared with the values of E reached in the upper (E mean = 0.0459) and lower (E mean = 0.0458) sections of the topographic slope. These different values of E, corresponding to all the experimentally studied gravity currents, are presented in figure 6. The increase in irregularities along the topographic surface appears to be related to higher entrainment values. Although the highest velocities (and highest Re) correspond to the flat slope, it suggests that the Re and the irregularities on the bottom contribute to the mixing between the gravity currents and the ambient fluid. The relation between the initial density contrast and the entrainment values is insufficient to elucidate a clear relationship between the g' and the entrainment. The highest values of E correspond to the gravity current with g' = 18 cm s-2 propagating over the lower section of the topography; in that case, E ≈ 0.049. Nevertheless, smaller and larger density contrasts (g' = 11 cm s-2 and g' = 24 cm s-2) are used in the same section with a very similar mean result for entrainment values (E ≈ 0.045 for both values of g'), as can be seen in figure 6 and on Table I.
3.3 Pressure time series
The analysis of pressure time series allowed us to conceive a method for identifying gravity currents in the slope of mountains. We compute the time series from the pressure field estimates at the stations shown in figure 1 (c). Six pressure time series are obtained for each gravity current run in the lower section of the topographic profile. These six series correspond to each synthetic station shown in figure 1 (c). One set of time series (corresponding to the run TL01) is shown in figure 7, where higher frequency pressure fluctuations start successively as the GC reaches the stations. The GC arrived at station 1 before we started recording because pressure fluctuations have started since the beginning of the time series on that station. However, from stations 2 to 6, almost no pressure variations are seen at the beginning of the time series. Higher frequency variations are seen at successive times: At ≈ 0.8 s in station 2, at ≈ 1.6 s in station 3, at ≈ 2.8 s in station 4, at ≈ 3.8 s in station 5, and at ≈ 4.6 s in station 6. Similar behavior is observed in the pressure time series for all the runs. Additionally, we perform numerical test simulations, shown in the supplementary material, that corroborate this behavior qualitatively but support the existence of surface pressure variations as the gravity currents flow.
In addition, figure 7 shows similar pressure variations just before 4 seconds in stations 2, 3, and 4. This fast and high amplitude (≈ 0.5 Ba) pressure variation could be explained in terms of a hydraulic jump. The average Froude number in this particular gravity current (TL01), Fr =0.98, is very close to the critical value (Fr =1), and the pressure signal is similar to the one presented by previous authors in the suitable conditions for hydraulic jumps to develop (De Padova et al., 2018). However, further investigation is necessary to rigorously confirm the existence of a hydraulic jump in gravity current TL01.
More significant variability and energy in stations 2, 3, and 4 can also be confirmed from their standard deviations (0.19, 0. 23, 0.19 Ba, respectively) compared to stations 1, 5, and 6 (0.08, 0.12, 0.11 Ba, respectively), and greater power spectral density in amplitude spectra (Fig. S5 in supplementary material).
3.3 Qualitative flow visualization
The GCs are visualized to qualitatively characterize their shape and observe the structures that lead to mixing. This section presents visualizations of GCs with the lowest (TL02) and highest (TL05) density contrasts. The images in this subsection belong to runs TU01, TL02, and TL05. Figure 8 (a) shows a snapshot of a GC corresponding to run TU01. This photograph was taken using red light and dark ink. Clearly, a Kelvin-Helmholtz instability is formed at the interface. These instabilities are associated with mixing enhancement between the GC and the ambient fluid (Strang and Fernando, 2001). Shadowgraph visualizations are shown in Figures 8 (b) and (c), where the main structures of the gravity current’s head, tail, and vortex can be seen behind them. Two sets of snapshots, corresponding to gravity currents with different g' in the lower section of the topographic slope, are shown. The set of snapshots of Figure 8 (b) corresponds to the largest g' = 24 cm s-2 (run TL05), and the right panels (c) correspond to the smallest density contrast: g' = 8 cm s-2 (run TL02); thus, a fainter shadow and light contrast can be seen when compared with the left panels. This behavior can be attributed to the fact that, in the shadowgraph technique, smaller contrasts of lights and shadows correspond to smaller density gradients (Panigrahi and Muralidhar, 2012).
In addition, in the visualization corresponding to run TL02, the vortices have more defined boundaries with sharp edges. In comparison, more mixing appears to form in run TL05, where the vortices are more diffuse. Additionally, we are interested in visualizing the propagation of the gravity currents as they interact with the irregularities of the topography. For both currents, when the head of these currents is traveling upwards in a topographical irregularity (Figure 8 b1, c1, and b5), it appears to be more defined, and fewer vortices are visible. When the currents travel downwards over a topographical irregularity, more vortices and contrasts of shadows and lights result (Figure 8 b3, b4, c3, and c4), suggesting that mixing may increase.
4. Conclusions
In this study, we investigated experimentally downsloping gravity currents over varying synthetic topography, along which we identified and characterized the enhanced mixing regions by using the temporal and spatial variability of the entrainment coefficient, E. The synthetic topography is a simplified scaled profile analogous to one on the hillside of a real mountain, which, to the best of our knowledge, has yet to be experimentally studied with the PIV techniques presented in this investigation. We used the entrainment coefficient as a proxy for the development of the mixing between the gravity current and the ambient fluid based on previous definitions adapted for our PIV experiment measurements and a non-flat slope (see Methods). We perform several experiments with different density contrasts, and we find that the entrainment consistently increases as the roughness of the topographic bottom increases. The lowest values of E are found for the flat slope experiment (F01) (E mean = 0.041) and for the intermediate (and flattest) section of the synthetic hillside-like-topography (E mean = 0.043). For a given density contrast (g' = 14 cm s-2), the highest values of E correspond to the upper and lower sections of the topography (E mean = 0.0459 and E mean = 0.458 respectively), corresponding to the regions where the roughness increases. While we observe a dependence of E with the slope roughness, we did not find any clear dependence between the density contrast and the entrainment values. The dependence between E and the roughness may also be qualitatively identified in the shadowgraph images in figure 8, where more mixing occurs when the gravity currents propagate over the irregularities of the slope.
From the velocity fields obtained with the PIV, we derived the pressure fields, which represent our second quantitative analysis. For the six synthetic stations on the lower section of the topographic slope, we analyze the pressure time series showing that the arrival of the gravity current can be observed. High-frequency variations in pressure arise as the current propagates over the synthetic stations, suggesting that this analysis can be implemented for real atmospheric gravity currents at particular locations of mountain slopes. Gravity currents in the atmosphere usually have much more complex dynamics than in experiments. However, these experimental models help understand the most basic underlying dynamics of real geophysical flows. The present study contributes to this understanding by considering a more realistic synthetic topography and its influence on mixing processes. In addition, the signal of the pressure time series we found is a novel aspect in analyzing experimental gravity currents. Its applicability in thermally driven winds on the slope of mountains should be further investigated in future research.