Observations of large-scale coherent structures in gravity currents: implications for flow dynamics

Density driven flows, also known as gravity currents, comprise a head, body, and tail. Yet whilst the body typically forms the largest part of such flows, its structure remains poorly understood. In this work, experimental data gathered using particle image velocimetry enables the instantaneous, whole-field dynamics of constant-influx solute-based gravity currents to be resolved. While averaged turbulent kinetic energy profiles are comparable to previous work, the instantaneous data sets reveal significant temporal variation, with velocity measurements indicating large-scale wave-like motions within the body. Spectral analysis and dynamic mode decomposition, of streamwise and vertical velocity, are used to identify the frequencies and structures of the dominant motions within the flow. By considering an idealised theoretical density profile, it is suggested that these structures may be internal gravity waves that form a critical layer within the flow located at the height of the maximum internal velocity. Irreversible internal wave breaking that has been postulated to occur at this critical layer suggests formation of internal eddy transport barriers, demonstrating that new dynamic models of turbulent mixing in gravity currents are needed.


Introduction
Gravity currents, also known as density currents, are a common class of geophysical flow that occur in many natural and man-made environments (Ungarish 2009;Simpson 1997). They are of particular relevance to the study of atmospheres and oceans, with examples including thunderstorm outflows, and sediment transport in lakes and oceans (Parsons and Garcıa 1998;Britter and Linden 1980;Simpson 1997;Bonnecaze et al. 1993;Talling 2014). Owing to their prevalence, and the fact that they are the primary mechanism for the transport of sediment, solutes, and heat in oceans (Dorrell et al. 2019;Talling 2014), extensive research has been conducted to establish the structure and dynamics of gravity currents through both experimental (Ellison and Turner 1959;Hacker et al. 1996;Kneller et al. 1999;Hallworth et al. 1996;Middleton 1966;Gray et al. 2005) and numerical investigations (Cantero et al. 2007;Hogg et al. 2016;Meiburg et al. 2015;Özgökmen et al. 2004;Stacey and Bowen 1988).
Experimentally, gravity currents can be categorised by mode of generation. Two often-discussed examples are constant-flux and constant-volume flows. There are significant differences in structure between the two flow types (Nogueira et al. 2014;Gerber et al. 2010). In a constantflux current, the continuous replenishment of dense fluid results in the bulk of the head and body of the current remaining undiluted (Hallworth et al. 1996). This category of flow has a prolonged body, which is assumed to be statistically steady with only small variations in characteristic variables, such as velocity and density (Gerber et al. 2010;Kneller and Buckee 2000). In a constant-volume flow, this body section is much shorter with a proportionally far more prominent head section. Existing experimental research has primarily considered constant-volume type flows, and examined the structure of the head of the flow over periods ranging from 10 s to a few minutes (Hacker et al. 1996;Hallworth et al. 1996;Middleton 1966). However, in oceanic gravity currents the head typically is not present, or forms only a small portion of the flow. Such currents are either quasi-permanent (for example the flow resulting from the connection between the Mediterranean Sea and the Black Sea (Sumner et al. 2014)) or have been observed to persist for several hours, or even days, and as a result the flow is assumed to be predominantly statistically steady (Peakall and Sumner 2015;Parsons et al. 2007;Azpiroz-Zabala et al. 2017;Khripounoff et al. 2003;Simpson 1997). However, despite the body of gravity currents normally forming the bulk of the flow (Azpiroz-Zabala et al. 2017;Sumner et al. 2014;Özsoy et al. 2001), the structure of turbulence within the body remains poorly understood (Wells and Dorrell 2021).
There are two primary mixing processes that occur in the head of the flow (Hacker et al. 1996;Simpson 1997;Kneller and Buckee 2000;Ungarish 2009): vortices that form as a result of Kelvin-Helmholtz instabilities generated by shear between the current and ambient fluid; and ambient fluid incorporated into the current as the raised nose of the current propagates over a no-slip surface, leading to a three-dimensional lobe and cleft structure (Simpson 1997). Existing experimental work has primarily considered lock-exchange constant-volume type flows with short statistically steady sections, or are based on at-a-point or profile measurements of constant-flux flows which limit the analysis techniques available (Buckee et al. 2001;Gray et al. 2006;Islam and Imran 2010;Kneller et al. 1999;Cossu and Wells 2012;Buckee et al. 2001;Davarpanah Jazi et al. 2020). Properties such as instantaneous and time-averaged Reynolds stress and turbulent kinetic energy have been measured in several previous works (Buckee et al. 2001;Islam and Imran 2010;Cossu and Wells 2012;Gray et al. 2006). For example, Buckee et al. (2001) examined a time-averaged profile of turbulent kinetic energy, and identified shear from the mean flow as the primary source of turbulence in the body of gravity currents, and Kneller et al. (1999) used instantaneous velocity fluctuations and Reynolds stresses to suggest the presence of eddies with size on the same order as the height of the body.
The currently accepted structure of gravity current body velocity and density profiles are shown in Fig. 1a (Kneller and Buckee 2000;Buckee et al. 2001;Sequeiros et al. 2010;Dorrell et al. 2019;Kneller et al. 2016;Abad et al. 2011;Sequeiros et al. 2010). The body is generally assumed to be two-dimensional and quasi-steady, with structure that can be divided into upper and lower layers by the height of the velocity maximum (Gray et al. 2006;Islam and Imran 2010;Kneller and Buckee 2000). The upper layer structure is determined by density stratification and shear with ambient fluid, and the shape appears comparable with that of a wall bounded jet. Unlike the jet, however, the lower layer of the gravity current can be approximated as similar to an open-channel flow (Dorrell et al. 2019).
There has recently been increased interest in understanding the instantaneous structure of the body, with several works highlighting the importance of instantaneous wholefield measurements and in particular of identifying and quantifying large-scale coherent structures (Lefauve et al. 2018;Odier et al. 2009Odier et al. , 2012Odier et al. , 2014Krug et al. 2013Krug et al. , 2015. Interfacial instabilities, in the form of Holmboe waves, between the flow and the ambient, have been observed in a gravity current in an inclined duct (Lefauve et al. 2018). However, this work considered an exchange-type flow, with equivalent velocity magnitude in the current and ambient fluids (resulting in a significant increase in shear compared with currents with little ambient flow). Mixing in the body has been investigated through simultaneous velocity and density measurements by Odier et al. (2009Odier et al. ( , 2012Odier et al. ( , 2014 and Krug et al. (2013Krug et al. ( , 2015, though despite being described as statistically steady gravity current flows all of these works consider very similar domains with flows that expand downstream and current heights that scale with the outlet size. As the domain used was not long enough to allow a transition to turbulence driven by the instabilities expected in gravity current flows, an array of active grids was employed to enhance turbulence near the inlet. Further, some of these works have low temporal resolution (Odier et al. 2009(Odier et al. , 2012(Odier et al. , 2014. Moreover, others only considered a small volume limited to the current/ambient interface ( 4 × 4 × 2 cm from a domain 200 × 50 × 50 cm) (Krug et al. 2015). In the same experimental domain, Neamtu-Halic et al.(2019) identified vortical Lagrangian coherent structures in the flow body capable of affecting the height of the turbulent/non-turbulent interface, though measurements were again limited to the current-ambient interface (this work considering 4 regions of 9 × 9 × 4 cm) precluding identification of structures lower in the flow. As a result of the limitations of these works, significant questions regarding the structure of the gravity current body remain.
Recent field measurements and laboratory experiments have suggested that the current gravity current model may need to be revised to a dynamic version considering forcing of flow-scale turbulent structures (Dorrell et al. 2018;Kostaschuk et al. 2018;Best et al. 2005). Flow-scale mixing by periodic internal gravity waves has been postulated to explain the structure of field-scale gravity currents in data collected from the body of a natural saline gravity current (see Fig. 1b) (Dorrell et al. 2019). A thorough understanding of the structure of the body of gravity currents is critical for accurate predictions of flow duration and interaction with the surroundings (Azpiroz-Zabala et al. 2017;Kneller et al. 1999).
In this paper, particle image velocimetry (PIV) is used to generate non-intrusive whole-field measurements of the instantaneous velocity structure of constant-influx, solute-based gravity currents. Quantification of the turbulence structure within the pseudo-steady body is used to improve existing understanding of gravity current flows. Therefore, it enables a far more detailed analysis of the nature of turbulence and flow structure within the body of a density current. Specifically, the key aims are to assess: (1) whether the pseudo-steady body of gravity currents can be described by flow-scale structures, (2) how these structures change with flow Reynolds number, (3) what these structures imply for our existing understanding of gravity currents, and (4) how such structures interact with the environment. The experimental setup is described in Sect. 2, followed by discussion of the turbulence structure in Sects. 3 and 4.

Methodology
The constant-flux gravity current experiments are conducted in a tank 0.1 m wide, 0.2 m deep and 2 m long (schematic shown in Fig. 2). The system is designed such that fluid leaves the outlet at the same rate that it is pumped in through the inlet. Raised sections, added at either end, prevent air entrainment through the inlet or outlet reaching the measurement region. The bed-slope ( ) for these experiments is set to 0.1 • . A large sump at the outlet prolongs flow duration by slowing the rate of current fluid pollution into the ambient fluid. The tank is initially filled with ambient solution, a 6% by mass solution of glycerol. The dense fluid, a 6% by mass solution of potassium dihydrogen phosphate (KDP), is pumped in at a constant rate through the inlet using a positive-displacement gear pump to provide a steady inflow with an inverter to control the flow rate. A coarse mesh with holes of diameter 7.8 mm is fitted over the inlet to provide a homogeneous inflow. Before entering the tank, the dense fluid first passes through a bubble trap, removing any air entrained by the gear pump. The bubble trap consists of a 1 m long, 0.1 m Fig. 1 a The currently accepted idealised structure of gravity current body velocity and density profiles (Kneller and Buckee 2000;Davarpanah Jazi et al. 2020;Altinakar et al. 1996;Garcia 1994;Abad et al. 2011;Sequeiros et al. 2010) and b a postulated flow structure from Dorrell et al. (2019) based on field scale gravity current measurements and comparison with zonal jet flows (Dritschel and Scott 2011). In (b), the coherent structure associated with large scale mixing is equivalent to a wave depending on frame of reference. The presence of dispersive waves leads to momentum transport due to anti-diffusive mixing and radiation stresses (Dorrell et al. 2019). Internal waves break close to the critical layer, leading to deposition of angular momentum and flow acceleration diameter cylinder filled with dense fluid. The dense fluid is pumped in to the top of this cylinder and out at the bottom, removing any small bubbles. The top and back of the tank are covered in black aluminium polyethylene composite panels to improve the image quality and contain the laser light.

Refractive index matching
Two 150 L mixing tanks are used to mix the ambient and dense fluids, which have a ≈ 3% density difference (see Table 1). These two fluids, as well as a mixture of the two, are refractive index matched, as required for PIV (Alahyari and Longmire 1994). While other pairs of fluids also meet the required criteria, such as solutions of ethyl alcohol and sodium chloride (Hannoun et al. 1988) or epsom salts and sugar (Mcdougall 1979), as in Alahyari and Longmire (1994) glycerol and KDP were chosen for their lack of volatility, comparatively low cost, and ease of handling. The fluid concentration and refractive index matching is tested using both a Reichert AR200 digital refractometer and an Anton Paar DMA TM 35 Basic density meter as well as by monitoring the temperature (which can have a significant impact on refractive index). The refractive index of each fluid is required to be equal to the reference value in Table 1 to the precision of the refractometer (5 significant figures) and constant across 3 consecutive readings at least 5 minutes apart. The density of the fluids is allowed to vary from the reference values due to temperature variation, with readings in the ranges 1012.9 ± 0.2 kg m −3 for the glycerol solution and 1043.0 ± 0.5 kg m −3 for the KDP. A single pair of RI matched fluids is chosen to closely and consistently match viscosity (see Table 1). The refractive index matching requirements therefore set the density difference between the current and ambient fluids. The fluid pair selected restricts the experiments to sub-critical flow ( Fr < 1 ). The range of experimental conditions possible is defined by the minimum and maximum discharge of steady, pulseless flow.

The PIV system
Planar particle image velocimetry (PIV) is used to generate measurements of downstream and vertical velocities (Raffel et al. 2018;Adrian and Westerweel 2011). Silver-coated hollow glass spheres are used as PIV seeding particles, with mean diameter of 10 m and density of 1400 kg m −3 , at concentrations of 0.0015 g L −1 and 0.0014 g L −1 for the current and ambient fluids respectively. The Stokes velocity of these seeding particles, and their relaxation time, r = d 2 p p 18 ≈ 6.50 × 10 −9 s, suggest that they follow the flow sufficiently for them to be suitable for use as PIV seeding particles.
A central vertical plane, parallel to the flow direction, is illuminated using a 532 nm Nd:YAG laser (with maximum energy 50 mJ) pulsed at 100 Hz. The images are captured using a DANTEC Dynamics SpeedSense camera with a Zeiss ZF.2 50mm f/1.4 lens with aperture set to f/2.0 that captures approximately 0.3 m of flow horizontally and 0.18 m vertically (see Fig. 2). The images are captured at a rate of 50 Hz in single-frame mode, and processed using DAN-TEC Dynamic Studio version 6.4 adaptive PIV. The resulting velocity field is mapped to a grid with spacing 3 × 3 mm, with temporal separation 0.02 s.  Table 1 Details of the density, , kinematic viscosity, , and refractive index, n, of 6% by mass solutions of the ambient (glycerol) and dense (potassium dihydrogen phosphate) solutes in tap water at 20 • C, from Haynes (2014) Solute Glycerol (ambient fluid) 1012.0 1.14 × 10 −6 1.3400 Potassium dihydrogen phosphate (current fluid) 1041.4 1.09 × 10 −6 1.3400 120 Page 6 of 18

The experimental cases
A series of experiments are conducted to establish the effect of influx (and therefore Reynolds number) on gravity current dynamics, and in particular the turbulence structure of the body of the gravity current. The influx rates for this set of experiments are shown in Table 2, along with characteristic length L c , velocity U c and time t c scales, and Reynolds ( Re = U c L c ∕ ) and densimetric Froude ( where is the kinematic viscosity, and g ′ is the reduced gravity. Here the characteristic scales are defined using averaged profiles. Case 7 is not considered further due to lower data quality compared with Cases 1 to 6, however plots including this case are presented in the supplementary material for completeness. In order to calculate time averages of the body data ( U , V ), the location of the body must be identified. Here, body data are defined by measuring the time taken for the current front to travel across the measurement region, and then waiting that length of time again before averaging across all downstream locations and time. For all cases, this definition results in consistent downstream velocity averages, whether averaging over 5 s or 20 s of data, suggesting that the data are approximately quasi-steady.
The characteristic velocity scale is taken to be the maximum value in the vertical profile of downstream velocity averaged over downstream position and body timesteps, U c = U max . The length scale is chosen to be the Ellison and Turner integral length scale (Ellison and Turner 1959), defined as where L is the height of the tank and Ū is the mean velocity relative to that in the ambient. The Ellison and Turner length scale is observed to be different from the current height (where the downstream velocity averaged over time U = 0 ). The characteristic time scale is the ratio of the two scales, t c = L c ∕U c . The kinematic viscosities of the dense , and ambient fluids are very similar (see Table 1) and for these calculations we use the viscosity of the dense fluid. These characteristic scales lead to Reynolds and Froude numbers that are output parameters, and therefore a doubling of influx does not result in a doubling of Reynolds number. The experiments captured between 35 and 55 s of data, including both the head and body of the flow, dependent on the rate of pollution of the ambient as a result of mixing. This flow duration is shorter for the higher influx cases owing to the higher rate of mixing, and the smaller ratio of ambient fluid to dense fluid.

Results
As seen in some previous works (Kostaschuk et al. 2018;Best et al. 2005), the height of the velocity maximum oscillates over time for influx values greater than 0.11 L s −1 (see Fig. 3), which may be linked to cross-stream flow, the presence of low frequency waves, and enhanced turbulent mixing (Dorrell et al. 2018). Mean velocity profiles collapse by shifting velocities by dividing by U c , and the vertical location by subtracting the averaged height of the velocity maximum, and then dividing by L c , For all cases, this normalisation unexpectedly results in average non-dimensional current height ≈ 1 . The downstream velocity averages for data in the body of each case are shown in Fig. 4a, with averaged (dU * ∕dt) and dU * ∕dY * for two cases in Fig. 4b and 4c. Average flow in the ambient is negative with magnitude O(10% → 20%) of the downstream flow, resulting in a slight increase in shear compared with a current with no ambient flow. This is comparable to flows considered in previous works (Kneller et al. 1999;Gray et al. 2005;Cossu and Wells 2012) and some real-world Ellison and Turner integral length scale (Ellison and Turner 1959),  (Sumner et al. 2014). The addition of lines illustrating the average height of the velocity maximum, the height of minimum dU * ∕dY * , and current height h, defined as the point where the average downstream velocity changes from positive to negative, allows closer comparison between the two cases. Both cases have local maxima in dU * ∕dt * (Fig. 4b) that coincide with the height of the local minimum in dU * ∕dY * (Fig. 4c), suggesting flow acceleration at this height. The higher influx case also has a local maximum in acceleration at the average height of the velocity maximum, though the lower influx case is decelerating at this height. This acceleration suggests that the body is not statistically steady as typically assumed (and as used to identify body data earlier in this paper). While increasing with increased Reynolds number, in the flows considered here the magnitude of the acceleration is small relative to the mean flow such that averages over the time windows considered are constant to leading order. However, the effect of such acceleration on long duration gravity currents remains an open question, beyond the scope of this paper. Subtracting the averaged profiles in Fig. 4a from the instantaneous data in the body gives the downstream and vertical fluctuations from the mean ( U * � = U * − U * , V * � = V * − V * ). As in previous work (Buckee et al. 2001), cross stream velocity fluctuations W * � are assumed to be small compared with U * � and V * � and the averaged normalised turbulent kinetic energy is defined as k * = 0.5(U * �2 + V * �2 ) . There are significant similarities between the profiles of time-averaged normalised turbulent kinetic energy presented in Fig. 4d and those presented in previous studies (Buckee et al. 2001;Islam and Imran 2010;Cossu and Wells 2012;Gray et al. 2006), specifically a local minimum just above the velocity maximum and a local maximum between the velocity maximum and current height, here close to the height of maximum negative shear. The maximum value of k * is highly dependent on the time range chosen for averaging (suggesting significant temporal variation) though the profile shape remains similar. All cases here have a maximum value of k * located below the velocity maximum.
In order to examine the structure within the body, the instantaneous velocity data and k * = 0.5(U * �2 + V * �2 ) are presented. Two representative cases that illustrate the changes in structure with Reynolds number are considered -a low influx case with Q = 0.09 L s −1 , and a high influx case with Q = 0.16 L s −1 . Fig. 5 shows instantaneous downstream and vertical velocity overlaid with velocity streamlines at a central non-dimensional downstream location, i.e. fixed X * (where X * = X∕L c and X is downstream location), as a function of Y * within the body over time. This central location is defined as the centre of the images, with the exact location varying slightly from case to case but always ≈ 1.5 m from the inlet. In both cases, these plots illustrate the instability of the interface between the current and ambient fluids. Additionally, the streamlines highlight the presence of wave-like motions within the body of the flow. These flow features are also evident in the supplementary video, depicting velocity components overlaid with velocity streamlines in the whole measurement region over time. Figure 6 shows similar plots of vertical velocity fluctuations over a shorter time range to allow closer inspection of the flow. In the low influx cases there is a regular pattern within the current body (Fig. 6a). As Reynolds number increases, this pattern becomes less regular and higher frequency. Additionally, as influx increases, a second regular structure appears at the height of the current with a frequency that decreases with increasing Reynolds number (Fig. 6b). Figure 7 shows k * at a central downstream location for the same two cases. These plots suggest significant temporal variation within the body. For the lowest influx cases, this takes the form of intermittent peaks at the height of maximum negative shear (marked with a dashed line). As influx increases, intermittent peaks begin to appear at the height of the velocity maximum. By Q = 0.16 L s −1 , these intermittent peaks are also present at the upper interface. Time-averaging    Fig. 4a) at a central downstream location (defined as in Fig. 3) turbulence statistics may therefore discard data important for the understanding and modelling of turbulence structure. As V * is small compared with U * , U * 2 max is representative of the kinetic energy in the mean flow. As the magnitudes of k * and k * are small ( ≪ 1 ), the energy contained within turbulent fluctuations is small compared with that in the mean flow.
Having established the presence of spatio-temporal structure in the body, suitable techniques to analyse this structure are required. Some motions with significant impact on the flow can be identified purely by examination of the plots of V * � . Figure 6 shows motions with period ≈ 3 s above the velocity maximum in the 0.09 L s −1 case, and motions with a period ≈ 5 s on the upper interface in the 0.16 L s −1 case. It is therefore expected that motions with frequency ≈ 0.33 Hz and ≈ 0.2 Hz respectively have a significant impact on the flow. For a more detailed analysis of these motions and their structure, we here employ discrete Fourier transforms and Dynamic Mode Decomposition (DMD). These techniques will be used in combination, with the aim of providing a thorough understanding of the flow structure, and will be applied to the data in dimensional form.

Fast Fourier transform
A discrete Fourier transform in time (Briggs et al. 1995;Iftekharuddin and Awwal 2012) is used to decompose a signal of N data snapshots that are equally spaced over time into the frequencies that make up that signal, →̂ where ̂ is the data transformed into the frequency domain. Here, a central downstream location is selected and a fast Fourier transform over time performed using the MATLAB fft function (MATLAB 2020) on the body data (defined as the data used to calculate the averages shown in Fig. 4a). Figure 8 shows the frequencies that dominate the flow (where significant frequencies are determined by overlaps between peaks in the Fourier transform data and peaks in dynamic mode amplitude, shown in Fig. 9). For Q = 0.09 L s −1 , the FFT identifies key frequencies of interest: 0.35 Hz, 0.50 Hz, 0.65 Hz, and 0.85 Hz. As the frequency of these modes increases, the height of the relevant mode gradually moves from the height of maximum negative shear to slightly above. There is far more variation in the FFT data for Q = 0.16 L s −1 . A mode with frequency 0.2 Hz is present at the height of the current. Compared with the lower influx case, other key modes are more difficult to identify. There is an indication of modes at 0.35 Hz and 0.50 Hz at the height of maximum negative shear, along with a variety of higher frequency modes with smaller magnitude at the same height and that of the height of maximum downstream velocity. However, this analysis does not show us the structure of these motions.

Dynamic mode decomposition
To consider the structure of the motions, DMD is utilised to decompose the body data (defined as in the FFT analysis) into modes in the form of waves with particular frequencies (Schmid 2010;Tu et al. 2014;Kou and Zhang 2017). A linear relationship is assumed between N data snapshots 1∶N = { 1 , 2 , ..., N } separated by time t, As the relationship is linear, the eigenvalues of A contain the dynamical characteristics of the system (Tu et al. 2014). However, the dimensionality of A is very large, so for practicality and accuracy in this work a similar matrix to A , Ã , is constructed, with fewer dimensions (Kou and Zhang 2017). This can be done using singular value decomposition, in this case the MATLAB svd function (MATLAB 2020), which is a generalisation of eigendecomposition to a non-square matrix. This decomposes a data matrix of size p × q into 2 unitary matrices, C and D , with sizes p × p and q × q respectively, and a diagonal matrix, , of size p × q containing the 'singular values' where H indicates a Hermitian transpose. As they are similar matrices, the eigenvalues and eigenvectors of Ã , j and j , are a subset of those of A , and contain the same dynamical information. These eigenvalues and eigenvectors are used to calculate the dynamic modes j = C j , the modal angular frequencies j = Im(log( j ))∕ t (and corresponding modal physical frequencies f j = j ∕2 ), growth rates g j = Re(log( j ))∕ t , and amplitudes a j = 1∕||F(∶, j)|| (where F = D −1 ). Each mode is a wave of the form j = a j e (g j +i j )t (Schmid 2011;Tu et al. 2014;Kou and Zhang 2017;Richecoeur et al. 2012). Figure 9 shows the amplitudes of the dynamic modes for the Q = 0.09 L s −1 and Q = 0.16 L s −1 cases. These mode amplitudes align well with the FFT analysis data, highlighting significant modes with similar frequencies in the same locations within the flow. Figure 10a and 10b show the structures of those modes. As indicated by the FFT data, the vertical location of these waves varies between the height of the current and the height of the velocity maximum. For the Q = 0.09 L s −1 case, the waves are at the height of maximum negative shear. For Q = 0.16 L s −1 , there are also waves at the current height and the velocity maximum.
Similar modes can be identified across all cases, and their frequencies and phase speeds compared. Figure 11a shows the wavelength of particular modes plotted against Reynolds number. These wavelengths were estimated by inspection of the mode velocity and vorticity plots (as in Fig. 10). The cases with Re < 3250 contain modes with a smaller range of frequencies than the higher Re cases. The frequencies of most modes stay roughly constant, with only small increases as Re increases. The frequency of those modes present at both high and low Re however experiences a significant increase when Re increases beyond ≈ 3750 . Figure 11b shows the wavelengths of those modes against frequency, and demonstrates two distinct categories of wave with the modes collapsing onto two separate lines. When plotting the phase speed, c = f , of the modes against Re as in Fig. 11c, those modes on the lower line in Fig. 11b all have c ≈ 0.025 m s −1 . Those modes on the upper line have a wider spread of phase speeds, between 0.06 m s −1 and 0.1 m s −1 . Those cases with Re < 3250 exclusively have modes with c ≈ 0.025 m s −1 , while the higher Re cases contain signs of internal waves with a wider spread of phase speeds.
To demonstrate whether the observed waves are due to buoyancy, inspection of the density profile would be advantageous. This data is not available from the PIV experiments, and instead simplifying assumptions are employed to obtain a heuristic estimate of the Brunt-Väisälä buoyancy frequency, N, which is the upper limit on angular frequency of internal waves due to buoyancy gradient (Sutherland 2010), Fig. 9 Plots of DMD mode amplitude against mode frequency for cases with (left) 0.09 L s −1 and (right) 0.16 L s −1 influx, with circles indicating the modes plotted in Figure 10a and 10b. These amplitudes were calculated based on DMD of 20 s body data for these two cases (as defined in Figure 4a) where is the average density profile, and 0 is taken to be the mean of the glycerine and KDP densities. An excess density distribution similar to that in Fig. 1 for solute based flows is assumed based on profiles obtained through threedimensional direct numerical simulation, with constant excess density ( e = − a , where e is the excess density, and a = 1012 kg m −3 is the density of the ambient fluid) both above the current height and below the velocity maximum and a linear distribution between the two. Above the current, the excess density is taken to be e = 0 . Below the velocity maximum, the excess density is estimated by requiring conservation of density flux between the inlet and the data. The inlet excess density flux ( F I ) is estimated by taking the product of the fluid influx, Q, and the excess density of the KDP. If the velocity and density profiles are assumed to be constant in the cross-stream direction, the excess density flux from the data ( F e ) may be estimated by (where W is the width of the tank). The excess density below the velocity maximum is estimated for each case by requiring that F I = F e (Table 3), and hence an approximate density profile for each case is established. As the maximum velocity will slow towards the side walls, these are likely underestimates of the maximum density within the body. Figure 12 shows the percentage difference between our observed wave frequencies and the estimated buoyancy frequencies after applying the Doppler shift due to the mean flow (Sutherland 1999), where N DS is the frequency measured by a stationary observer, and k x the downstream wavenumber, as it has been assumed that the waves propagate purely downstream. As seen in Fig. 12 many of the observed waves are within 25% of the estimated values. Given that the vertical location and wavenumber of the wave, and therefore the mean flow speed and hence Doppler shift, are approximate, that any wave propagation in the cross-stream or vertical directions would act to decrease the observed frequencies, and that the exact density profile is unknown, all frequencies observed here are the right order of magnitude and are considered sufficiently close to indicate that these may be internal waves due to buoyancy. The greater difference between the observed waves and the estimated values included in Fig. 12b may be attributed to larger uncertainty in the height of the wave compared to those in Fig. 12a, or could be the result of some more complex three-dimensional structure. Figure 13 demonstrates that the phase speed of some of these potential gravity waves is approximately equal to the mean downstream flow speed at the wave height. This suggests the presence of critical layers, defined to be where U = c (Maslowe 1986), within the flow.

Flow dynamics
Fourier transforms and dynamic mode decomposition have been used to identify the most energetic motions in   Table 3 Details of the inlet excess density flux F I calculated by taking the product of the fluid influx and the excess density of the KDP, and the estimated maximum excess density within the body calculated by requiring F I = F e (where F e is defined in (7) the body of gravity current flows. This analysis has been used to identify potential flow-scale internal gravity waves centred on the height of maximum negative shear, and the velocity maximum. As Reynolds number increases, the frequencies of the dominant waves within the flow change.
The modes at the height of maximum negative shear become less significant, and the wavelengths of the most energetic modes decrease and their frequencies increase. At this point we also find the height of the velocity maximum starts to vary, and modes on the velocity maximum become significant in terms of flow dynamics. It has been shown that for some of these waves the wave phase speed is approximately equal to the mean flow speed, indicating a potential critical layer within the gravity current body. The presence of internal waves in the gravity current body has been postulated by Dorrell et al. (2019), who suggest that the gravity current body has a structure similar to that of a zonal jet (Fig. 1b) (Dritschel and Scott 2011;Rossby and Zhang 2001;Maxworthy 1984;Bower and Hogg 1996). In zonal jets, the breaking of dispersive internal waves near a critical layer results in self-organisation of the flow and net momentum transport towards the jet core (Dritschel and McIntyre 2008;Dorrell et al. 2019;Bühler 2014;Dritschel and Scott 2011). Close to critical layers, breaking waves homogenise potential vorticity and steepen the potential vorticity gradient (Dorrell et al. 2019). Unless they have sufficient strength, this steep gradient is difficult for eddies to penetrate. The gradient therefore acts as a barrier to mixing between the flow above and below the critical layer, preventing dilution of the lower part of the flow and sharpening the density profile. Thus, the presence of a critical layer may result in the maintenance of far stronger density stratification than predicted by existing models (Booker and Bretherton 1967;Maslowe 1986;Dritschel and McIntyre 2008;Dorrell et al. 2019). This would lead to faster flow velocities and longer flow durations than expected from the current theory. Internal waves absorbed at the critical layer transfer horizontal momentum into the mean flow (Booker and Bretherton 1967;Maslowe 1986;Thorpe 1975), increasing the mean downstream velocity. This would imply that the gravity current body may not be statistically steady as typically assumed (Gray et al. 2006;Kneller and Buckee 2000;Islam and Imran 2010).
Indeed, acceleration of the flow at the height of the potential internal waves has been identified in this work.  As this acceleration is largely in the upper part of the flow, this acceleration would change the profile of dU/dt, with the maximum downwards shear moving further from the velocity maximum and closer to the current height over time. The presence of these waves thus may help to explain discrepancies between data from a real-world flow in the Black Sea (Dorrell et al. 2019) and predictions from traditional models of the gravity current body. As indications of waves have here been identified even in flows over smooth surfaces, and a rough surface would act to increase the prevalence of gravity waves (Aguilar and Sutherland 2006;Sarkar and Scotti 2017), their existence may be an inherent characteristic of gravity current body flow. This reinforces the need for a dynamic model of body flow, allowing for sharpening of density and velocity profiles.
The aspect ratio of the flows in this work is comparable to those in previous works quantifying the velocity structure (Cossu and Wells 2012;Gray et al. 2005;Islam and Imran 2010); however sidewall drag may affect the flow. As Reynolds number increases the effect of side walls is expected to be reduced, yet a greater range of frequencies with high amplitude dynamic modes have been identified in the higher Reynolds number cases. Therefore, the identified structures are likely not dominantly a result of sidewall interactions. However, future work considering a wider domain would strengthen the conclusions presented in this work. Similarly, identification of internal gravity waves from the presented data would be bolstered by the addition of density measurements. While beyond the scope of this work, the addition of simultaneous, instantaneous, non-intrusive density measurements conducted using a technique such as laser induced fluorescence (such as that used by Ermanyuk et al. (2017)) or a background oriented schlieren method (as described by Verso and Liberzon (2015)), or consideration of an alternative fluid pairing with a different expected buoyancy frequency, would allow conclusive identification of the observed structures as internal gravity waves.
As the Reynolds numbers of the flows considered in this work are significantly lower than those of real-world flows, it is important to consider trends in the strength of identified modes as Reynolds number increases. The magnitudes of the vertical velocity fluctuations are equivalent in the two cases (Fig. 6), while the amplitudes of the dynamic modes are significantly higher in the higher Reynolds number case (Fig. 9). Similarly the magnitudes of the Fourier transform of downstream velocity, and of the downstream velocity component of the dynamic modes, are equivalent in the two cases (Figs. 8 and 10). However, as Reynolds number increases the magnitudes of the Fourier transform of vertical velocity, and of the vertical velocity component of the dynamic modes, is smaller than the lower Reynolds number case. This may be a result of increased cross-stream velocity associated with these modes. Therefore, in order to identify a clear trend in mode strength with increased Reynolds number, further investigation including this cross-stream velocity component is needed.

Implications for natural flows
Real-world thermohaline gravity currents (Ivanov et al. 2004;Legg et al. 2009), for example the flow at the Strait of Gibraltar, are a crucial and common class of geophysical flow responsible for driving oceanic circulation. Furthermore, comparison can be made between the structures of thermohaline flows and of sediment-driven turbidity currents (Kneller and Buckee 2000;Garcia 1994;Moodie 2002). It is, therefore, interesting to consider how applicable the results presented in this work are to such sediment-laden flows. It has previously been claimed that solute-based flows are dynamically similar to fine-grained conservative gravity current flows (Kneller and Buckee 2000;Cossu and Wells 2012). While the structure of coarse-grained non-conservative gravity currents may differ in some respects (Stacey and Bowen 1988;Kneller and Buckee 2000;Hogg et al. 2005;Cossu and Wells 2012;Wells and Dorrell 2021), for example an expected decrease in current height, increase in downstream velocity, and greater stratification of the density profile in the lower part of the current, the potential for internal gravity waves to form requires only that stable stratification of density exists (Staquet and Sommeria 2002). The formation of a critical layer requires only that the phase speed of these waves equals the flow speed (Maslowe 1986). It is, therefore, reasonable to suppose that similar structures may be present in sediment-laden flows also.

Conclusions
Large-scale coherent structures are identified within the pseudo-steady body of laboratory-scale constant-influx solute-based gravity currents. Inspection of the velocity field reveals wave-like structures in the body, prompting the use of dynamic mode decomposition to characterise the dominant motions. The presence of these structures during flow over a smooth surface suggests that they are an inherent characteristic of the body flow due to instability at the upper interface. Estimation of the Brunt-Väisälä buoyancy frequency, using an idealised expected density profile, suggests that these structures may be internal gravity waves. Observations indicate that, due to internal velocity variation, these waves may support development of a critical layer within the flow, which generates a barrier to mixing; prevents the dilution of the lower part of the current; and changes the expected structure of the density and velocity profiles. The presence of internal gravity waves suggests that the body may not be statistically steady as typically assumed. Instead, internal waves, which can transfer horizontal momentum within a flow, and deposit it at critical layers, are likely critical to explain gravity current dynamics.