Development of gravity currents on slopes under different interfacial instability conditions

We present experimental results on the development of gravity currents moving onto sloping boundaries with slope angles $\unicode[STIX]{x1D703}=7^{\circ }$ , $10^{\circ }$ and $15^{\circ }$ . Different regimes of flow development are observed depending on the slope angle and on the initial velocity and density profiles, characterized by the Richardson number $J_{i}=\unicode[STIX]{x1D6FF}_{i}{g_{0}}^{\prime }/\unicode[STIX]{x0394}u_{i}^{2}$ , where $\unicode[STIX]{x1D6FF}_{i}$ , $\unicode[STIX]{x0394}u_{i}$ and $g_{0}^{\prime }$ are, respectively, the velocity interface thickness, the maximum velocity difference and reduced gravity at the beginning of the slope. For $J_{i}>0.7$ and the larger slope angle, the flow strongly accelerates, reaches a maximum at the beginning of the Kelvin–Helmholtz instability, then decelerates and re-accelerates again. For $0.3<J_{i}<0.6$ , instability occurs earlier and velocity oscillations are less. When $J_{i}\leqslant 0.3$ the increase in velocity is smooth. The magnitude of velocity oscillation depends on the combined effect of $J_{i}$ and slope angle, expressed by an overall acceleration parameter $\overline{T_{a}}=(\unicode[STIX]{x1D6FF}_{i}/U_{i})((U_{c}-U_{i})/x_{c})$ , which, to first order, is given by $J_{i}\sin \unicode[STIX]{x1D703}$ , where $U_{c}$ and $x_{c}$ are, respectively, the velocity and position at instability onset. The velocity increases smoothly up to an equilibrium state when $\overline{T_{a}}\leqslant 0.06$ and exhibits an irregular behaviour at larger values of $\overline{T_{a}}$ . The critical Richardson number $J_{c}$ decreases with increasing $J_{i}$ (increasing $\unicode[STIX]{x1D6FF}_{i}/h_{i}$ ) which is due to wall effects and $\unicode[STIX]{x1D6FF}/\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D70C}}\neq 1$ . After the beginning of Kelvin–Helmholtz instability, entrainment rates are close to those of a mixing layer, decreasing to values of a gravity current after the mixing layer reaches the boundary. It is shown here that the interfacial instability during current development affects the bottom shear stress which can reach values of $c_{D}\approx 0.03$ regardless of initial conditions. By solving numerically the depth integrated governing equations, the gravity flow velocity, depth and buoyant acceleration in the flow direction can be well predicted for all the performed experiments over the full measurement domain. The numerical results for the experiments with $J_{i}>0.3$ predict that the current requires a distance of at least $x_{n}\approx 40h_{i}$ to reach a normal state of constant velocity, which is much larger than the distance $x_{n}\approx 10h_{i}$ required in the case of a current with $J_{i}\leqslant 0.3$ that is commonly assumed for downslope currents.

or interleaves at the level of neutral buoyancy. Along its path, ambient water is entrained, diluting the current and thus changing its dynamics and the properties of the water masses. A quantification of these dilution processes is relevant for example to understand the transport of dense water from the Antarctic shelf toward the abyssal plain, which is part of the meridional overturning circulation (Baines & Condie 1998). Gravity currents are also important in morphodynamics, causing sediment erosion, sediment transport and creating bed forms (Garcia & Parker 1993;Sequeiros et al. 2010;Zordan et al. 2018). A knowledge of the bottom drag of developing as well as equilibrium state currents is therefore essential for understanding these transport processes.
Most laboratory experiments considered gravity currents on horizontal or constant, small slope boundaries, using lock exchange, finite volume releases or constant supply conditions (Simpson & Britter 1979;Britter & Linden 1980;Simpson 1982;Altinakar, Graf & Hopfinger 1990;Cenedese & Adduce 2010;Ungarish 2011). Studies of constant supply gravity currents on sloping boundaries, with slope angles varying over a large range, have focused on the normal, constant velocity state of the current (Ellison & Turner 1959;Britter & Linden 1980). This state is considered to be reached at a distance from the supply of ≈10h i (Odier, Chen & Ecke 2014), where here h i is the initial depth of the current. The flow up to this equilibrium state and the dependency of the current development on initial conditions and slope angle, including the spatial development of the current due to topography changes, has received little attentions. Pawlak & Armi (2000) considered the development region of an arrested wedge flow over a sill on slopes up to 14.5 • and found entrainment rates in the initial region two to three times larger than those known for a developed current on these slope angles. Similar results have been obtained by Odier et al. (2014) for a constant inflow light gravity current moving up a sloping boundary of 10 • . It may be expected that the current development will depend on initial interfacial instability conditions in addition to the slope. However, until the paper of Negretti, Flòr & Hopfinger (2017), this has been ignored. Further, excepting Negretti et al. (2017), no simultaneous measurements of the bottom drag and entrainment in gravity currents on slopes considering their mutual interaction have been conducted. Gravity currents may often be in a developing state and the properties of the developed current will depend on its history i.e. interfacial entrainment and bottom drag during flow development. Field measurement by Van Haren et al. (2014), in the Romanche fracture channel, clearly indicate that different initial conditions at the beginning of the slope drastically influence the further development of the current. In a recent paper Zordan, Schleiss & Franca (2019) also point to the importance of initial conditions of turbidity current development. Negretti et al. (2017) considered only flow development of large interfacial Richardson number before moving on to concave or linear slopes. The resulting down-slope current development exhibited a behaviour completely different from that of low initial Richardson number, considered by Pawlak & Armi (2000). Three distinct development regions have been identified: (i) nearly free fall acceleration, (ii) beginning of Kelvin-Helmholtz instability (KHI) with high entrainment and related high bottom drag, (iii) subsequent collapse of KH billows and readjustment of the current to a force equilibrium. Under these conditions gravity currents, even at x ≈ 30h i , (distance of experimental observation) were far from equilibrium, i.e. a normal state, constant velocity current. Similar observations have also been reported by Kostaschuk et al. (2018), which include data from numerical and laboratory studies and from field observations. In nature, slope changes are frequent and consequently, gravity currents may often not be in a normal state, contrary to what is generally assumed to be the case (Dallimore, Imberger & Ishikawa 2001;Danabasoglu, Large & Briegleb 2010). There is therefore sufficient motivation to clarify the effects of the flow conditions of the current at the beginning of the slope on its further development on the slope.
We present here experimental results demonstrating the dependency of initial conditions on the current development. Different flow development regimes are observed that depend on the value of an overall acceleration parameter (introduced here) which also determines the distance required to reach an equilibrium, constant velocity regime. The slope angles considered, i.e. 7 • to 15 • , are relevant for ocean applications (see also Pawlak & Armi 2000) and with these small slopes we cover already the current behaviour characteristic of large slopes, expressed by a relatively large value of acceleration parameter, i.e. large velocity oscillations with large entrainment and bottom drag; when this parameter is small the approach of an equilibrium state is smooth (no velocity oscillations). We present measurements of both bottom drag and interfacial entrainment, whereas most previous studies focus on one only, because these two processes are often considered to be uncorrelated. We show here that the bottom drag is correlated with interfacial instability (interfacial entrainment) and can thus take values much larger than estimated from boundary layer considerations. For equilibrium currents a correlation between bottom drag and entrainment has been considered by Dallimore et al. (2001). Using the estimated entrainment and drag coefficient (E and c D ) we solve the governing equations reporting a good agreement with the experimental data and allowing us to predict the distance required for the current to reach the equilibrium.
The paper is organized as follows: in § 2 we present the experimental set-up, and the measurement techniques. In § 3 the initial conditions are defined and discussed. Section 4 includes the derivation of the governing equations and in § 5 we classify the different flow regimes depending on the characteristic parameters of the spatially evolving flow. Using the entrainment and drag estimated in § 6, we compare the experimental data to the numerical solution of the governing equations in § 7. Section 8 summarizes the results and includes concluding remarks.

Experimental set-up
A schematic of the experimental set-up is shown in figure 1. The buoyancy driven flow has been generated using saline solutions injected at the bottom of one channel end at constant flow rate Q = qb, via a pump, from an external 800 l reservoir. The first portion of the channel is elevated by 30 cm, and is horizontal of total length of L = 230 cm and width b = 25 cm. The end of the horizontal portion is followed by a sloping boundary of slope angle θ on which the flow accelerates. Once the pump was switched on, the outlet at the opposing end of the tank was opened in order to enable the flow to evacuate to prevent a return flow and ensure a constant total water depth H 0 . The whole tank downstream of the gate was filled with fresh water with density ρ a to a level of H 0 = 18 ± 0.5 cm from the bottom of the elevated channel and the portion upstream of the gate was filled with salt water. The density ρ s of the saline solution was measured just before each experiment with a densimeter (DMA TM 35, Anton Paar   figure 1 for the notations) of the particle image velocimetry (PIV) velocity measurements. The flow rate imposed by the pump was Q p = 0.4 l s −1 , the total water depth in the channel H 0 = (18 ± 0.5) cm and the reduced gravity g 0 = g((ρ s − ρ a )/ρ a ) = 4.3 ± 0.1 cm s −2 were kept equal in all the experiments. D represents the distance of the gate from the beginning of the slope. The gate opening was h 0 = 6 cm and was completely removed in D0 experiments (arrested wedge condition, AW). The Richardson number is J i = g 0 δ i /( u i ) 2 . h i and U i are the current height and average velocity at 5 cm upstream of the beginning of the slope and given by the integrals (4.2.a,b) (see figure 1). The overall and integral scale Reynolds numbers are Re = uh/ν ∼ 3000, Re L = u 2 √ 15/(ν ) ∼ 60, respectively. The field of view (FOV) was between −10 cm < x < 76 cm. In the experiments marked with an asterisk * also runs with fluorescent dye (Rhodamine 6G) were conducted. The distance x c , measured from the beginning of the slope, corresponds to beginning of KHI. of up to 10 % of the dense current velocity. In addition, experiments with complete removal of the gate at D = 5 cm (experiment D0) were run and were thus similar to those of Pawlak & Armi (2000), with h 0 ≈ H 0 /2. The return flow was in these experiments was larger, approximately 30 % of the dense current velocity. The slope angles were varied at 7 • , 10 • and 15 • . The depth integrated initial velocity U i , the integral current height h i (see (4.2.a,b)) and the interface thickness δ i , given in table 1, have been determined from the PIV experimental data at a distance of 5 cm upstream of the beginning of the slope (x = 0 cm) (see figure 1). The corresponding hydraulic initial conditions are expressed in terms of an initial Froude number Fr 2 = U i 3 /B 0 = 0.94 ± 0.07, where B 0 = g 0 q is the initial buoyancy flux, q = U i h i is the initial flow rate per unit width and g 0 = ((ρ s − ρ a )/ρ a )g is the initial reduced gravity.

Measurement technique
The velocities were determined using the optical, non-intrusive experimental technique of PIV. The PIV set-up consisted of a light source, light sheet optics, seeding particles, a camera and a PC equipped with a frame grabber and image acquisition software. Polyamide particles (PA12, Vestosint 1101 white) with a diameter in the range d p = 100-250 µm with a specific density of ρ p = 1.060 g cm −3 were added in both sides of the channel containing fresh and salt water as tracer material for the PIV measurements. It turned out that an appropriate particle density in the water body could be achieved by adding ≈0.04 g particles per litre. Because of the small density difference between the two reservoirs and the given density distribution of these specific particles, it was possible to use these as seeding material for both the fresh water and the salty water. With the estimated particle time response τ p = d 2 p ρ p /18µ ≈ 4.8 × 10 −4 and the Kolmogorov time scale τ K = (ν/ ) 0.5 ≈ 0.04, with ≈ u 3 /h ≈ 0.01 cm 2 s −3 , the Stokes number was estimated to be St = τ p /τ K 1 thus particles are considered to follow the fluid streamlines closely.
A 6 Watt argon-ion laser (Coherent) operating in multimode (λ b = 488 nm, λ g = 514 nm) has been used as a continuous light source. The beam was transmitted through a fibre optic cable to a line generator with spherical lenses (OZ Optics Ltd., Nepean, Ontario). The generated laser sheet on the slope had a length of approximately 1 m and a width of 5 mm and was positioned in the middle of the channel.
Images of 86 cm × 65 cm were captured with a CCD camera (FlowMaster3, 14 bit, 1600 × 1200 pixels) at a frame rate of 23.22 Hz and an exposure time of 9000 µs. A wide angle lens (SIGMA AF EX 1.8/24 DG Macro AF for Nikon) was used at a distance of 2 m to the field of view, leading to a spatial resolution of 0.540 mm pixel −1 . Time series of 300 s leading to 7000 raw images were stored in real time on a raid system. The raw images were then processed using a cross-correlation PIV algorithm to compute the velocity fields, each from two consecutive raw images, using the software package DaVis (LaVision). An adaptive multipass routine was used, starting with an interrogation window of 64 × 64 pixels and a final window size of 16 × 16 pixels with 50 % overlap. Each vector of the resulting vector field represents an area of approximately 0.43 cm × 0.43 cm. The velocity vectors were post-processed using a local median filter and checked by the distinctiveness of the highest correlation peak. Given the velocities encountered in the experiments, the experimental error in the instantaneous velocity is estimated to be approximately 3 %. Further error estimations from the experiments are given in appendix A.
For dye visualizations Rhodamine 6G was added to the salt water of the gravity current with a concentration of less than 20 µg l −1 , and they were used to estimate the density profiles and some averaged values of the density field by normalizing locally each vertical section with the maximal value.

Initial conditions and parameters
The initial density interface is strongly affected by the gate position and configuration. Figure 2   showing time sequences of vertical sections at x = −5 cm highlighting interfacial density fluctuations for experiments D2-7 (a), D1-7 (d) and D0-10 (g). It is seen that the density interface fluctuations increase as the release distance (distance of gate before the beginning of the slope) is reduced: in experiments D2-7 and D1-7, the dye interface indicates the existence of Holmboe waves before the beginning of the slope, decaying in conditions D2, while for experiment D0-10 KHI appears close to the beginning of the slope. The initial vertical profile of the streamwise velocity is also strongly affected by the gate position and configuration. Figure 2(b,e,h) shows time-averaged vertical velocity profiles, with the confidence interval on the time average (grey shadow), along with a fit with the hyperbolic tangent function (black symbols). While for experiments D1-7 (e) and D0-10 (h) the fit reproduces well the experimental profile, in experiment D2-7 (b) the fit fails in the upper part of the profile and close to the bottom. The discrepancy at the bottom is due to the development of a boundary layer on the relatively long horizontal boundary and in the outer part viscous diffusion and some mixing by Holmboe waves affect the velocity profile.
Holmboe instability (HI) is expected to exist when the velocity shear layer thickness δ is larger than the thickness of the density interface δ ρ and Hazel (1972) demonstrated that HI exists at any weak shear provided δ/δ ρ > 2. The Holmboe instability is characterized as two sets of waves with one layer cusping into the upper layer and the other into the lower layer, and propagating in the opposite direction. The cuspids can be seen in the insets of figure 2(a,d). Theoretical studies by Haigh & Lawrence (1999) of salinity stratified flows showed that when the bulk Richardson number J, defined as is less than 0.071, KHI dominates, even when δ/δ ρ > 2, in agreement with the experimental findings by Koop & Browand (1979). Zhu & Lawrence (2001) observed symmetric HI in the subcritical region of an exchange flow over a sill, which decayed when J > 0.65-0.75 depending on Reynolds number. Similar results have been reported by Strang & Fernando (2001) and Hogg & Ivey (2003).
To evaluate the Richardson number J, it is necessary to measure the shear interface thickness δ. Since in the D2 experiments the velocity interface cannot be fitted by a tanh distribution, we determined δ by using the same method as in Pawlak & Armi (2000), i.e. the best linear fit through the normalized velocity profile between the 15 % and 85 % values, where u is the velocity normalized by the local velocity difference u = u(z)| max − u(z)| min . The initial values δ i are given in table 1 together with the initial Richardson number values. In experiment D2, J i ≈ 0.72, the Holmboe waves are of much lower amplitude than in experiment D1 where J i ≈ 0.53 (figure 2a,d,g). This agrees with the observations by Zhu & Lawrence (2001), i.e. Holmboe waves decay when J > 0.7. Figure 2(c,f,i) shows the change in shear layer thickness on the slope defined as the 15 %-85 % bounds of the maximum velocity u m (black continuous lines) in the downstream direction along with the Reynolds stress u w , which will be discussed in § 6.2.

Governing equations
Following the procedure of Ellison & Turner (1959), Negretti et al. (2017) derived the governing equations of the depth integrated velocity and density on slopes including angle variations. The assumptions made were two-dimensional stationary flow with the boundary layer and Boussinesq approximations. The x-momentum, the mass and buoyancy flux conservation equations are where c D = τ w /ρU 2 and E are respectively the bottom drag and entrainment coefficients, and τ w the bottom shear stress, θ the slope angle, q = Uh is the flow rate per unit width and B = g Uh the buoyancy flux, which is equal to the buoyancy flux supplied at the gate, B 0 . The parameters S 1 and S 2 are the empirical shape factors which will be discussed in more detail in appendix B. The depth integrated velocity U, the corresponding depth h and the depth integrated reduced gravity g are quantities defined making use of the similarity hypothesis, which implies the following relations: that depend on x only. We note that the system of differential equations (4.1) reduces to a shallow water model (Ungarish 2011) when fixing S 1 = S 2 = 1, which is equivalent to considering 'top-hat' density and velocity profiles. All dimensional variables (x, h, u, q and g ), denoted below by·, are made dimensionless using the characteristic length and velocity scales L = (q 2 0 (4.3) Equations (4.1) represent a system of three differential equations with the unknown variables U, h, g and can be combined and rewritten to finally give Note that c D is of order E and cannot be neglected (Negretti et al. 2017;Kostaschuk et al. 2018). As outlined in Negretti et al. (2017), when the Richardson number at the beginning of the slope is large, J i ≈ 1, the interface is KH stable at the beginning of the slope and in the initial acceleration region, up to x c , E = c D ≈ 0. In this case (4.4) can be simplified to which has the implicit solutioñ The integration constant C 0 = (1 + 2S 1 cos θ ) is determined imposing the initial conditionŨ =Ũ i , forX = 0 and assuming the hydraulic condition at the beginning of the slope for the Froude number Fr = 1. We point out that (4.6) differs from the free fall solution presented in Negretti et al. (2017) by the pressure term multiplied by S 1 . The different theoretical solutions and their comparison with experiments will be presented and discussed in § 7.

Flow development regimes
In this section time-averaged quantities are presented. For each experiment, all values have been averaged in time by taking into account only the steady layer flow (after the current head has passed). A moving averaging filter of spatial width L is applied and data are sub-sampled with the same spatial scale to improve the clarity of the figures. A detailed estimation of the errors is given in appendix A.
Dimensionless time-averaged velocities U/U i and current heights h/h i as a function of x/h i are presented in figures 3(a,c) and 3(b,d), respectively for experiments at θ = 15 • (a,b) and 7 • (c,d). Different symbols and grey scales correspond to the release distances (D2, A ; D1, 6 ; D0, @ ). The grey shades in figure 3(a,b) represent the confidence interval of the time-averaged velocities. Figure 3(a,c) shows that in experiments D2 and D1, a strong acceleration of the current occurs until a certain distance, x c . At this position, in D2 ( A ) experiments, the current velocity reaches a maximum and then decreases rapidly, whereas in D1 ( 6 ) experiments, an abrupt change in velocity increase occurs at x c . In D0 ( @ ) conditions, only a small change in velocity increase at x c is observed. This critical position x c corresponds to beginning of KHI, as already observed by Negretti et al. (2017), and is marked by vertical arrows in figure 3(a,c). They vary from x c /h i ≈ 6-8 (15 and 7 • respectively) in D2 experiments, to x c /h i ≈ 4 in D1 experiments and to x c /h i ≈ 2 in D0 experiments. Similar considerations hold for the current height h (figure 3b,d): a strong decrease until a minimum at x c /h i is observed in experiments D2 and D1, while in the D0 configuration the depth of the current remains nearly constant over the full measured domain.
From these results it is clear that the conditions at the beginning of the slope strongly affect the development of the flow on the slope. As already mentioned in § 3, the current released far upstream from the slope (D2 experiments) has a large interfacial Richardson number J i ≈ 0.67-0.77 (cf. table 1) at the beginning of the slope, implying that the Holmboe waves, that are present some distance downstream of the gate, are decaying (Zhu & Lawrence 2001). When the gate is closer to the beginning of the slope, the initial Richardson number is J i = 0.53 (D1 experiments) and the interface remains Holmboe unstable at the beginning of the slope. In the D0 experiment, J i ≈ 0.3 and the current is KH unstable shortly after the beginning of the slope. Also, the effect of the slope is evident. In general, larger slopes generate larger changes in velocity and depth of the current, while on the smallest slope, 7 • , the change in velocity is smoother. Figure 4 gives the variation of J as a function of x/h i for the experiments of slopes 15 • (a) and 7 • (b). Again, the different symbols correspond to the different release distances (D2, A ; D1, 6 ; D0, @ ). It is seen that the larger is J i , the lower is the critical value of the Richardson number: J c ≈ 0.24 for D0, ≈0.2 for D1 and ≈0.13 for D2.
Onset of KHI depends on the ratio δ/δ ρ and on the wall proximity, expressed by (δ/h) c . As has been shown by Hazel (1972), the critical Richardson number J c decreases with increasing ratio δ/h, as evident from figure 4(c) from our experiments. In the D2 experiments, J i is large (weak Holmboe wave mixing) and (δ/h) c ≈ 1, hence a low value of J c and the current will accelerate over a longer distance to reach this value. In the D1 experiment (δ/h) c < 1 and Holmboe waves are more intense and can mix the density interface so that δ/δ ρ −→ 1, hence a larger value of J c is expected. In the D0 experiment J c ≈ 0.24 because δ/δ ρ ≈ 1 and (δ/h) c ≈ 0.5. Figure 5 shows the normalized velocity profiles u/U i before the beginning of the slope at x = −5 cm (a) and at the location of minimum Richardson number, at x c 1.0 for slope angle 15 • (b). At the beginning, the velocity profiles are similar (with the exception of the D0 experiment, which has a larger return flow-lock exchange effect), but at the location of minimum Richardson number, the velocity maximum of D2 experiment is much closer to the bottom boundary. In experiments D1 and D0 the thickness of uniform velocity is much larger. The variation of the ratio δ/h is shown in figure 6 for experiments with 15 • (a) and 7 • (b). In the D0 and D1 experiments δ/h 1, while in experiments D2, δ/h ≈ 1. In the latter, the KH billows are thus larger and closer to the wall, causing boundary layer separation and large c D , as has been shown by Negretti et al. (2017). Negretti et al. (2017) explained the current development on slopes by introducing an overall acceleration parameter T a = (h i /x c )( U/U m ), with the velocity difference U from the beginning of the slope to x c . They suggested that when T a > 0.05 velocity oscillations exist whereas when T a < 0.05 the velocity increase is regular, approaching an equilibrium state within approximately 10h 0 . In the more general case, where the initial conditions are varied (which was not the case in Negretti et al. 2017), a physically more appropriate definition of an overall acceleration parameter is This acceleration parameter can be expressed in terms of the initial Richardson number and slope angle by relating U c to the critical Richardson number J c and getting x c from assuming free fall up to x c . This gives  where C is a function of the ratio of u m /U and of (u m − u min )/U. In experiments D1 and D2, C ≈ 1.1 and in D0, with larger back flow, C ≈ 1.3, giving respectively T a ≈ 0.8J i sin θ and ≈1.1J i sin θ. A similar dimensionless number is defined in Negretti, Socolofsky & Jirka (2008), using a temporal linear stability, normal modes, analysis of a two-layer stratified exchange flow. There it is shown that the growth rate of the instability depends on the acceleration parameter. Figure 7 shows T a as a function of J i sin θ for all the performed experiments, including those of Negretti et al. (2017) and Pawlak & Armi (2000). The pre-factor 0.8 shown in the figure by the dashed line compares reasonably well with the experimental value determined from (5.1). Although the definition of T a differs from that of Negretti et al. (2017), we see that transition from no velocity oscillations to oscillations of increasing amplitude with increasing T a is at T a ≈ 0.06. In experiment D2-15 for instance, with T a ≈ 0.12 velocity oscillations are large. When T a 0.05, D0, the increase in velocity up to an equilibrium state is smooth.
6. Entrainment and bottom drag 6.1. Determination of the entrainment coefficient The entrainment rate is defined as the volume of ambient fluid engulfed into the current along its descent. It can obtained from the discretization of (4.1b), i.e. taking the difference between the input and the output fluxes in a control volume with a longitudinal width x, here taken equal to h i . The entrainment coefficient is thus (6.1) Figure 8(a,b) shows the variation in volume flux as a function of x/h i . As is evident from the data, a direct calculation of the gradient would give large variations. A primary source of error in these estimates of the volume flux gradient results from the slow variation at the origin of the instabilities. In order to reduce the variations, the gradient is determined using a moving averaging filter of spatial width L = L and sub-sampling the data with the same spatial scale. The dashed lines in figure 8(a,b) at x/h i > 10 indicate the expected entrainment rate of an equilibrium current. It is seen that in figure 8(a) there is good agreement whereas in the experiments with a slope of 7 • in figure 8(b) there is a relatively large variation because values are very small.
The entrainment coefficients obtained using (6.1) are shown in figure 8(c,d). In experiments D0 ( @ ) the highest entrainment coefficient is near the beginning of the slope, corresponding to the beginning of KHI, it then decreases along the slope down to a value of approximately 0.02 on slope 15 • , whereas in the experiments on slope of 7 • , E q decreases to approximately 0.005. This is in agreement with the behaviour observed by Pawlak & Armi (2000). While in experiments D2 and D1, the final values are similar, the highest values are observed shortly after KHI. These highest values of entrainment rates are close to 0.065 ± 0.005 for larger slope angles and 0.04 ± 0.1 for the smallest slope angles, which are values similar to what is observed in free shear layers.
Another estimate of entrainment is obtained from the velocity w h normal to the interface as proposed by Morton & Turner (1956). The z-integrated continuity equation To minimize errors, we estimated the velocity w h at the interface from the PIV velocity data by averaging three values of the normal velocity around the averaged value of h for each x-position and normalizing with the depth integrated velocity at the same location. We also compared the values obtained directly from averaged velocity fields and interface positions and from instantaneous velocity fields and interface positions averaged after. The differences between the two methods is small. Details on the experimental errors relative to the entrainment are given in appendix A. The entrainment rates E w are shown in figure 8(  differences in E q and E w . In general, the entrainment E w is larger than E q and after the first acceleration region and establishment of the shear layer, the two differ by a factor of two at x/h i > 10: E w ≈ 0.05 and E q ≈ 0.02 for the larger slopes and E w ≈ 0.02 and E q ≈ 0.01 for θ = 7 • . These values are similar to those obtained by Odier et al. (2014) (cf. their figure 13): E q ≈ 0.02 and E w ≈ 0.04. In figure 9 these differences are highlighted by plotting E q (dashed lines), E w (continuous lines) and T a = (h/U) dU/dx (grey dash-dotted lines) for experiments D2-15 (a) and D0-7 (b). After the acceleration region, in which entrainment is small or nearly inexistent in the absence of interfacial instabilities, in the D2-15 experiment both E q and E w increase and reach a maximum value of ≈0.07 and then E q decreases rapidly down to 0.02, while E w ≈ 0.05. In the experiment D0-7, we observe a discrepancy at the beginning, with E q > E w until the beginning of KHI, while the discrepancy between the two definitions is smaller further downstream in this case. The grey dash-dotted lines in both figures represent the acceleration parameter T a = (h/U)dU/dx, which can be seen as part of the entrainment, E q = T a + dh/dx, as outlined in Negretti et al. (2017). In experiment D2-15 (figure 9a), T a is large at the beginning, is negative after the beginning of KHI to become weakly positive at x/h i ≈ 12. For experiment D0-7 (figure 9b) we see that for 0 < x/h i < 4, T a ≈ E q , so that the current depth is nearly constant, and then decreases toward zero. Since at slope angles θ 15 equilibrium entrainment rates are very small, errors in measuring entrainment rates can be large, particularly in E q . Nevertheless, entrainment coefficients determined from the change in volume flux E q tend toward the values known for equilibrium flow conditions as seen in figure 10 by comparing (a) and (b) for experiments D2 and D0, respectively. Red arrows highlight the path toward the equilibrium values for the two release cases: it is direct in the D0 experiments with a smooth increase, while for the D2 experiments there is also an oscillation in the Richardson number before reaching a constant equilibrium value in figure 4(a). This behaviour is consistent with the variation with downstream distance of the Froude number Fr (figure 10c), which also tends toward the equilibrium values expected for the slopes considered. The expected equilibrium is represented in figure 10(a,b)   the empirical relations proposed by Turner (1986) and Cenedese & Adduce (2010) deduced from gravity currents at equilibrium; E w , in contrast, overestimates the values known for equilibrium flow conditions. Odier et al. (2014) suggest taking an average between E w and E q for the equilibrium entrainment value, which is hardly satisfying. Krug et al. (2013) also find E w = 0.04 for their gravity current on a slope of 10 • , and suggest that this is the correct entrainment coefficient at equilibrium state of the current. The difference with respect to Ellison & Turner (1959) is attributed to larger Reynolds numbers in their experiments. Visibly, this is not the case. In the present experiments we also obtain E w ≈ 0.04, but E q ≈ 0.02 when the equilibrium state is approached, in agreement with Pawlak & Armi (2000), van Reeuwijk, Holzner & Caulfield (2019), Odier et al. (2014) and even Krug et al. (2013). Indeed, from their figure 12, we obtain E q = dh/dx close to 0.02 in the equilibrium state at a distance ten times the inlet height, i.e. at x > 52 cm. In the present experiments E w ≈ E q at the beginning of KH instability when the mixing length m is expected to be large (large scale overturning), but when the equilibrium state is approached the mixing length decreases and E q < E w . Krug et al. (2013) present measurements of W over a large range of z and it is seen that W/U = 0.04 at the interface and drops to zero at z ≈ 0.3h away from it. This is indicative of a small mixing length; in Krug et al. (2013) m ≈ h/10. For the foregoing reasons we take E q as the physically correct entrainment coefficient.

Bottom drag coefficient
In viscous flows, c D is inversely proportional to the Reynolds number, as has been shown by Cenedese et al. (2004) for a gravity current of Re < 300. In the present experiments Re > 1400 (for D2-15, Re = Uh/ν ≈ 1700 at the beginning of the slope and is Re ≈ 2800 at x/h i = 15). Hence, the boundary layer region is turbulent after onset of KH instability. In this case we use known values and expressions for a turbulent boundary region. For a wall jet (that has a velocity profile similar to a gravity current) of Reynolds number Re m = u m z m /ν ≈ 10 3 , where z m is the height of  maximal velocity, u m ≈ 8 cm s −1 the maximal velocity and ν ≈ 10 −2 cm 2 s −1 the viscosity, c D ≈ 5 × 10 −3 (George et al. 2000). For a turbulent boundary layer c D is given by the empirical relation c D = 0.0113Re m −0.178 , which gives a value close to the bottom drag coefficient for a wall jet at the same Reynolds number. These bottom drag coefficients are small compared with the large values determined for gravity currents on slopes in the present experiments (see figure 11 and also 12) and those of Negretti et al. (2017). This is due to the interaction of KH billows with the boundary region and not due to boundary layer turbulence. Therefore, the known Reynolds number dependency of the boundary layer drag coefficient does likely not apply for gravity currents on slopes, whereas for gravity currents on horizontal boundaries, the classical boundary layer values apply with the drag coefficient tending asymptotically towards 10 −3 , when Reynolds numbers are very large (Cossu & Wells 2012).
Since accurate velocity measurements are difficult close to the bottom boundary, we estimated c D using an average of the Reynolds stresses u w and the local averaged velocity where δ B is the distance from the bottom boundary where the local velocity becomes maximum. Figure 11 shows  Negretti et al. (2017). The values are respectively c D ≈ 0.03 and 0.015 for the 15 • slope and 7 • slope. In order to highlight the relation between current acceleration, entrainment and bottom drag, we display in figure 12 the variation with downstream distance of the three quantities for experiments on slope θ = 15 • . Figure 12(a) suggests that the strongest decrease in velocity after its maximum value takes place when the maxima of entrainment and bottom friction are correlated (D2 experiment). In the D2 experiment the current accelerates when both entrainment and bottom drag are small. At maximum velocity, entrainment and drag start to increase rapidly and the current velocity starts to decrease to then increase again after maximum values of E w and c D . This figure also shows the correlation between entrainment and bottom drag. In D1 and D0 experiments (figure 12b), the rate of velocity increase changes when E w reaches a maximum and c D starts to increase. Here, there is a shift between maximum entrainment and maximum bottom drag and this is larger in the D0 experiment ( figure 12c).
In D2 experiments, the bottom drag c D increases when the entrainment coefficient E q increases. The link between the upper and the lower part of the gravity current in the D2 configuration is evident in figure 2: we observe that the normalized Reynolds stress u w /U 2 , at the interface and at the bottom, increases at the same location (with opposite signs). The spatial correlation coefficient that we estimated between these two terms is above 0.75 for the D2 configuration and below 0.25 for the D0 configuration. Dallimore et al. (2001) use the turbulent kinetic energy budget of Sherman, Imberger & Corcos (1978) for a uniform flow to link entrainment to the bottom drag. This correlation does not apply in the present strongly developing flow. In general c D = f (E q , δ/h, Fr, τ ν ), where τ ν is the viscous shear stress. Since E q depends on Fr and τ ν is negligible, the important parameters are E q and the wall proximity, expressed by δ/h, which introduces a shift in c D with respect to E q when δ/h < 1 at instability onset. Maximum c D occurs when the interfacial layer approaches the wall, i.e. when δ/h 1. The D2 experiments (figure 12d) indicate a linear relation between E q and c D of the form with B ≈ 0.6 and the initial value E i ≈ 0.01. In figure 12(d) the bottom drag coefficient c D is plotted as a function of E q and compared with the empirical expression (6.4). The experimental entrainment coefficients and bottom drag are used in the following section to solve numerically the integrated equation (7.1).

Comparison with numerical solutions of the governing equations
Equation (4.4) has been solved numerically to obtain the three unknown variables U, h, g using a Runge-Kutta scheme of order 4 (RK4) with a step dx = 10 −2 mm, using the experimental entrainment E = E q (x) and bottom friction c D (x). The shape factors were obtained from the 2.5 layer model derived in appendix B, and are S 1 = 0.5 and S 2 = 0.6. Neglecting the spatial variations of S 1 , S 2 and S 3 , the numerical solution has been obtained solving the following equation: dŨ dx = S 2Ũ 2 sin θ −Ũ 5 (c D + E) − S 1 E cos θŨ/2 − (q/R)(S 3 EŨ 4 + S 1 sin θŨ/2) q(U 3 − S 1 cos θ ) .  figure 13(c,f ) have been calculated indirectly asg = B 0 /(Uh) since the downstream spatial variation of the density fields could not be extrapolated from the fluorescent dye visualizations. The numerical solution of (7.1) is represented by continuous lines, black lines using E q and grey lines using E w . The free fall solution proposed in Negretti et al. (2017) is represented by dashed lines and the implicit shallow water (SW) solution (4.6) with S 1 = S 2 = 1 by dotted lines. We tested the numerical solution for all experiments and found a very good agreement, the best agreement for all variablesŨ, h andg using the numerical solution of (7.1), E q and c D from figure 11, over the full measured domain. The SW solution (4.6) reproduces fairly well the data forũ m and badly those relative to h in the first accelerating region, in which entrainment and bottom friction  can be neglected. It is curious to see how the free fall relation proposed by Negretti et al. (2017) reproduces well the data forŨ in the same region, neglecting the pressure term with respect to the SW model: the positive contribution given by the pressure term is compensated by the shape factor S 1 . Also, the critical distance x c is well detected for both the velocity and the current depth. We tested solving (7.1) initialized using the constant entrainment coefficient of Turner (1973) (not shown) and constant c D = 10 −3 and observed a bad agreement with the experimental data. Sensitivity tests with different shape factors show that the numerical solution forŨ, h andg strongly depends on these parameters, which justifies the importance in their proper estimation as we made in appendix B.
We can use the numerical solution to predict the distance needed to reach the equilibrium between the gravity and the drag terms x n where acceleration ceases and the flow attains a constant velocity. The dimensionless equilibrium velocityŨ n defined in Negretti et al. (2017) from their (2.11) and here obtained from (4.4) is where Ri = B 0 cos θ/U 3 . The right-hand side term of (7.2) is deduced from the differential equation for Ri in Ellison & Turner (1959) (their equation (13)) by assuming dRi/dx = 0.  In figure 14 the solutions of (7.1) are shown for experiments D2-15 (a) and D1-7 (b). The equilibrium values, given by (7.1) are indicated by the horizontal red lines, for experiment D2-15: Ri = 0.25, E = 0.02 and c D = 3 × 10 −3 (higher value ofŨ n , dash-dotted line) and c D = 3 × 10 −2 (lower value ofŨ n , dashed line), as extrapolated from figures 8 and 11. It is interesting to point out that the two relations proposed by Ellison & Turner (1959) for E(θ) and E(Ri) coincide only if assuming a large c D : the dashed line in figure 14(a) is equivalent using E(θ ) and E(Ri) only if a large value of c D = 3 × 10 −2 is used. Figure 14(b) shows the velocity variation for D2-7, taking in this case E = 0.015 and c D = 0.005 (cf. figures 8 and 11). Figure 14 indicates that the gravity current needs a distance of 40x to reach an equilibrium state, which is much larger than the commonly assumed distance of 10x. No extrapolation is needed for the D0 experiment because, as seen in figure 12(c),Ũ has practically reached its equilibrium value atx ≈ 10.

Further discussion and conclusions
The four important conclusions that emerge are: (i) the development of gravity currents on slopes can be regular or highly irregular depending on the interfacial instability conditions at the beginning of the slope and on slope angle and this is expressed by an overall acceleration parameter; (ii) the bottom drag coefficient is large and cannot be neglected even in an equilibrium current and, more importantly, it is related to the interfacial entrainment rate; (iii) the entrainment coefficients determined from the change in volume flux E q = (dq/dx)/U on the one hand and from velocity measurements at the outer interface of the current E w = −w h /U on the other, are different when flow equilibrium is approached, as already observed by Odier et al. (2014). While Odier et al. (2014) suggest taking an average between E w and E q for the equilibrium entrainment value, we give arguments, based on turbulent shear stress measurements, that the equilibrium value is E q , in accordance with Pawlak & Armi (2000); (iv), when the initial interfacial Richardson number and slope angles are such that T a is large, an equilibrium state is reached at a distance of approximately 50h 0 (instead the usually assumed distance of 10h 0 ).
The current development is shown to depend on slope angle and on the interfacial Richardson number at the beginning of the slope J i = δ i g 0 / u 2 i , where δ i , u i and g 0 are, respectively, the velocity interface thickness, the maximum velocity difference and reduced gravity at slope origin. When J i is small, here J i < 0.3, as is the case in outflow over a sill (Pawlak & Armi 2000) or inflow at the beginning of the slope (here D0 experiment), an equilibrium flow of constant velocity is reached at the expected downstream distance of approximately 10h i , where h i is the current height at the beginning of the slope estimated with (4.2.a,b). For this reason, the development phase is often neglected as in Dallimore et al. (2001) or Odier et al. (2014). However, as has already been demonstrated by Negretti et al. (2017), the flow development can be entirely different, approaching equilibrium much further downstream, when the initial interfacial Richardson number and slope angle are large. The flow development is shown to depend on an overall acceleration parameter T a = (δ i /U i )(U c −U i /x c ), which, to first order, is given by J i sin θ, where U i is the velocity at the beginning of the slope and U c and x c are, respectively, the velocity and position at instability onset. When T a < 0.06 the velocity increases steadily toward an equilibrium state, whereas at large values of T a the development is highly irregular with large velocity variation.
To summarize the current development, there are three different regimes depending on the overall acceleration parameter. When T a is large, T a ≈ 0.12, as in the D2-15 experiments, the current strongly accelerates up to x = x c , reaching a maximum of velocity at the beginning of KHI, then decelerates and finally re-accelerates again. In this case the velocity oscillations are large, due to the large scale, intense vortical structures at the start of KHI that cause sudden large entrainment and boundary layer separation because (δ/h) c 1, and thus spatially correlated high values of the bottom drag coefficient (see § 6). In the second regime, 0.06 < T a < 0.1, corresponding to the D1 experiments, the current undergoes an abrupt change in acceleration at x = x c . In this case, the interface is thinner, 0.3 < J i < 0.6 and δ/δ ρ ≈ 1, so that KHI occurs at higher J and thus at a shorter distance from the beginning of the slope. The third regime is characterized by T a < 0.06, here corresponding to the D0 experiments, and represents the case were a smooth monotonic increase of the velocity occurs, with an initial J i 0.3, and KHI is already present near the beginning of the slope. In this regime, an equilibrium, with constant current velocity, is generally reached within approximately 10h i .
Numerical integration of the depth integrated momentum equation, using experimental values of E and c D , reproduces well the current development and clarifies the validity of free fall assumption versus a shallow water model. Furthermore, the numerical solution has been used to predict the distance x n required to reach the equilibrium or normal state with an asymptotic constant velocity, assuming constant entrainment and bottom friction coefficients, as reported experimentally after re-adjustment of the current. When T a > 0.06 this distance is x n ≈ 40h 0 , and thus much larger than that predicted when assuming a current in the D0 regime with an unstable interface close to the beginning of the slope (J i < 0.3).
The entrainment coefficient has been determined from the change in volume flux and from a direct measurement of the velocity normal to the current velocity. The overall behaviour is similar, however, when the normal state is approached, with a difference of approximately a factor of two. This difference has also been measured by Odier et al. (2014), who chose the larger value obtained by the direct method. At slope angles θ 15, equilibrium entrainment rates are very small so that errors in the measurement of entrainment rates are large. Nevertheless, entrainment coefficients E q determined from the change in volume flux tend toward the values known for equilibrium flow conditions. These values have therefore been used in the numerical extension of the experimental results to predict the current equilibrium.
After beginning of the KHI at J c < 0.25, observed entrainment rates are close to those of a mixing layer, decreasing to values of a gravity current after the mixing layer has reached the bottom boundary. This large entrainment rate affects the bottom shear stress with the bottom drag coefficient reaching, in all cases, O(10 −2 ). When δ/h 1 at J c , the bottom drag is directly correlated with entrainment, whereas when δ/h < 1 at J c there is a lag in bottom drag, reaching its maximum value when the shear layer has grown down to the boundary. This large values of bottom drag during current development are of importance for applications where slope changes are frequent and, consequently, gravity currents may often not be in an equilibrium state. Also, large bottom drag is indicative of the capability of the current to pick up sediment. Gravity currents are often modelled by shallow water equations (Ungarish 2011;Wirth 2011) that assume a constant layer-averaged density and velocity and which lead to an overestimation of the driving force due to the pressure gradient by 33 %, and up to 50 % depending on the chosen approximated shapes of the profiles (Pokrajac, Venuleo & Franca 2017). Figure 15 displays time-averaged density (dashed lines) and velocity (continuous lines) profiles for the experiment D2-7 at x = −5 cm (a) and at x = 50 cm (b). The density and velocity profiles are close to the 'top-hat' shape upstream of the slope (figure 15a) as long as the interface is stable and little mixing at the interface takes place, whereas when the shear layer develops a triangular shape is approached for both density and velocity profile. For the experiments D0 the initial velocity and density profiles are closer to the 'top-hat' profile. After the beginning of the KHI, both profiles are very similar for all experiments and close to the triangular profile as highlighted in figure 15(b).
The shape factors introduced by Ellison & Turner (1959) allow us to take into account the deviation of the instantaneous velocity and density profiles from the simplified rectangular (top-hat) profile by assuming the self-similarity of these profiles, and hence to correct the shallow water solution.
The integral velocity U and current depth h of the current are deduced from the experimental observations by integrating the velocity profile as in (4.2c,d) over a depth h v at which the velocity has become a percentage of the maximum velocity u m . In this paper this percentage is set to 8 % based on the experimental resolution of the velocity data. Figure 16(a,b) displays the different depths: the integral velocity layer depth h (grey dashed line), the integral density layer depth h ρ (black dotted line) and the constant maximum density layer η B (black continuous line) for D2-7 (a) and D1-7 (b). We see that the depth of the density layer h ρ (dashed line in figure 16a,b) behaves similarly to h, while the bottom constant density layer η B (continuous line in figure 16a,b) always decreases in the downstream direction to reach a thickness of a few centimetres only, thus suggesting that a triangular shape is approached moving downstream.  FIGURE 16. Integral current depth h (grey dashed lines), constant density layer η B (black continuous lines), density profile height h ρ (black dashed-dotted lines), for D2-7 (a) and D1-7 (b). Shape factors S 1 (c,d) and S 2 (e,f ) for D2-7 (c,e) and D1-7 (d,f ). 'Top-hat' density profiles with a total thickness h ρ − η I /2 (black dotted lines), triangular density profiles (black lines), no assumption on the velocity profile, assuming both triangular density and velocity profiles (black dashed lines) and using the 2.5 layer model (continuous grey lines).
Using (4.2c,d) and assuming a 'top-hat' density profile with a constant reduced gravity g through the depth h ρ , we get the following relations for S 1 and S 2 : The corresponding shape factors assuming a triangular shape of the density profile are obtained multiplying the above values of S 1 and S 2 by 1/3 and 1/2, respectively. The shape factors estimated using 'top-hat' and triangular shapes of the density profile and no assumption on the velocity profile are illustrated in figure 16(c-f ) using dash-dotted and continuous lines, respectively. In the 'top-hat' case, the factors are both of O(1) over the full distance considered with standard deviations of 0.12 and 0.05 for S 1 and S 2 , respectively, larger for S 1 because of the squared relation in (B 1). The constant value around one for this shape confirms that the pressure and gravity terms are overestimated in this case. Using triangular shapes, the values start from 0.3-0.5 for S 1 and 0.55-0.7 for S 2 atx = 0 and increase in downstream direction to reach finally values of approximately 0.7 for S 1 and 0.75-0.9 for S 2 atx = 20 for both experiments, in accord with the observed development of the velocity and density profiles of the gravity flow.
So far, no assumption has been made on the shape of the velocity profile. If it is assumed, for example, that the velocity and density profiles are both triangular with a height h ρ = h v = (4/3)h, the shape factors S 1 and S 2 become S 1 = 16 27 g 0 g and S 2 = 2 3 g 0 g , (B 2a,b) and S 3 can be estimated by assuming u = u m (1 − z/h v ) from which follows (cf. equations (4.2c,d)): U = 2/3u m , h = 3/4h v and S 3 = 2/3. The shape factors obtained using triangular profiles for both the density and velocity profiles are represented in figure 16(c-f ) by the dashed lines. They differ from the continuous lines, where no assumption on the velocity profile has been made, in the first developing region (0 <x < 5), but vary similarly forx > 5. An alternative way to estimate the shape factors, which is more sophisticated than considering a two-layer system, as is the case in the SW model, is to take into account an additional layer of thickness δ I , as depicted in figure 17. The top layer is assumed quiescent and the shear layer velocity is supposed to be half of the velocity of the bottom layer, but the model can be easily extended to any other ratio between the velocities in the shear and bottom layer. Similarly, the density is assumed to be half of the density of the bottom layer in the sheared velocity layer, and constant, but the approach can be adapted to other interface configurations. Neglecting the drag terms we have, whereũ m = u m /U corresponds to the dimensionless maximum free-stream velocity, h B = h B /L is the dimensionless bottom layer thickness,q B = q B /q 0 =ũ mhB andq I = q I /q 0 =ũ mδI /2. The unknowns areũ m ,h B andδ I . Equation (B 3a) can be analytically solved to give the implicit solution, u 2 m + 2 h B +δ I 2 cos θ = 2x sin θ + ζ , with ζ =ũ 2 m0 + 2 h B0 +δ I0 2 cos θ , (B 4a,b) whereũ m0 andh B0 are estimated using the initial condition F 2 = 1 and are linked tõ U 0 andh 0 through Turner's integrals (4.2.a,b), u m0 = χ 0Ũ0 , (B 5a) h B0 +δ I0 /2 = χ −1 0h 0 , (B 5b) The variable χ 0 links Turner's variables (Ũ,h) to the free-stream layer variables (ũ m ,h B ). Combining (B 5b) and (B 5c) we get χ 0 = 2(h 0 − h 2 0 −h 0δI0 )/δ I0 . The shape factors can now be obtained as S 1 = χ −2 (g 0 /g ) and S 2 = χ −1 (g 0 /g ) and are represented in figure 16(c-f ) by the grey continuous lines. We see that S 1 = 0.56-0.87 and S 2 = 0.75-0.93. The so obtained values will be finally used to solve the system (4.1) since they are the most representative, considering both the shear layer thickness and the bottom layer, and are also shown to produce better agreement when comparing the numerical solution with the experimental data as in § 7.