2 Numerical model, test cases and model validation

The governing Navier–Stokes and density transport equations are solved in non-dimensional form with the channel height, $H$ , as the spatial scale and the buoyancy velocity, $u_{b}=\sqrt{g^{\prime }H}$ , as the velocity scale. The reduced gravity is $g^{\prime }=g(\unicode[STIX]{x1D70C}_{max}-\unicode[STIX]{x1D70C}_{min})/\unicode[STIX]{x1D70C}_{max}$ , where $g$ is the gravitational acceleration and $\unicode[STIX]{x1D70C}_{max}$ is the density of the colder fluid (with temperature $T_{min}^{\prime }$ ) situated in the region $x/H<0$ before the lock gate is removed. The other half of the channel ( $x/H>0$ ) initially contains warmer fluid with temperature $T_{max}^{\prime }$ and density $\unicode[STIX]{x1D70C}_{min}$ . The non-dimensional density is defined as $C=(T_{max}^{\prime }-T^{\prime })/(T_{max}^{\prime }-T_{min}^{\prime })$ , where $T^{\prime }$ is the dimensional temperature.

The finite-volume viscous solver (Pierce & Moin Reference Pierce and Moin2001; Chang, Constantinescu & Park Reference Chang, Constantinescu and Park2006) used to perform the simulations advances the governing (filtered) Navier–Stokes equations in time using a semi-implicit iterative method. The Boussinesq approximation is employed to account for stratification effects. The pressure-Poisson equation is solved using multigrid. The conservative form of the non-dimensional Navier–Stokes equations is integrated on non-uniform Cartesian meshes. All operators in the momentum and pressure equations are discretized using second-order central schemes. The algorithm is second order in time. Discrete energy conservation ensures robustness at relatively high Reynolds numbers despite using strictly non-dissipative (central) schemes to discretise the Navier–Stokes equations. A standard advection–diffusion equation is solved for the non-dimensional density. The QUICK scheme is used to discretise the convective term in the non-dimensional density equation. The two parameters in the non-dimensional governing equations are the channel Reynolds number, $Re=u_{b}H/\unicode[STIX]{x1D708}$ and the molecular Prandtl number, $Pr=\unicode[STIX]{x1D708}/\unicode[STIX]{x1D705}$ , in which $\unicode[STIX]{x1D705}$ is the molecular diffusivity. The subgrid-scale viscosity and the subgrid-scale diffusivity in the filtered non-dimensional momentum and temperature equations are calculated using the dynamic Smagorinsky model (Pierce & Moin Reference Pierce and Moin2001, Reference Pierce and Moin2004) based on the resolved velocity and non-dimensional density fields at each time instant. The model correctly predicts negligible values of the subgrid-scale viscosity and subgrid-scale diffusivity in regions where the flow in weakly or non-turbulent, even if the shear is significant, which is essential for accurate simulation of gravity current flows. The code was previously validated against flow over a bottom 2-D cavity initially filled with non-buoyant or buoyant pollutant (Chang et al. Reference Chang, Constantinescu and Park2006; Chang, Constantinescu & Park Reference Chang, Constantinescu and Park2007), constant-density flow in a channel containing an array of cylinders (Chang & Constantinescu Reference Chang and Constantinescu2015), lock-exchange currents propagating over a flat smooth surface (Ooi, Constantinescu & Weber Reference Ooi, Constantinescu and Weber2007a ; Ooi et al. Reference Ooi, Constantinescu and Weber2009), bottom gravity currents propagating over an isolated circular/rectangular cylinder situated at a small distance from the channel bottom (Gonzalez-Juez et al. Reference Gonzalez-Juez, Meiburg, Tokyay and Constantinescu2010), gravity currents propagating over a triangular dam (Tokyay & Constantinescu Reference Tokyay and Constantinescu2015) and gravity currents propagating over an array of bottom-mounted obstacles (Tokyay, Constantinescu & Meiburg Reference Tokyay, Constantinescu and Meiburg2011, Reference Tokyay, Constantinescu and Meiburg2012, Reference Tokyay, Constantinescu and Meiburg2014). Pierce & Moin (Reference Pierce and Moin2001, Reference Pierce and Moin2004) discuss detailed validation of the model for reacting turbulent flow simulations in which the transport equations are solved for both conserved (same as the equation solved for the non-dimensional density in the present study) and non-conserved scalars. Given that our level of mesh refinement is similar to the one used by Pierce & Moin (Reference Pierce and Moin2001, Reference Pierce and Moin2004), we expect the model will accurately predict mixing. Additional proof of the capability of the model to predict mixing for gravity current flows is given by the satisfactory comparison between the shape of the current predicted by 3-D LES and that obtained from video recordings of the evolution of the current in laboratory experiments (Ooi et al. Reference Ooi, Constantinescu and Weber2007a ; Ooi, Constantinescu & Weber Reference Ooi, Constantinescu and Weber2007b ; Ooi et al. Reference Ooi, Constantinescu and Weber2009; Tokyay & Constantinescu Reference Tokyay and Constantinescu2015).

The surface porous layer of height $h_{1}$ present in the left half of the channel $(-L<x<0)$ was populated with square cylinders of side length, $D$ (figure 1). Their axes were aligned with the spanwise $(z)$ direction. This is the main difference between the present numerical study and previous experimental studies of lock-exchange flows in vegetated channels (e.g. Zhang & Nepf Reference Zhang and Nepf2011), in which the vegetated layer was modelled with vertical circular cylinders. The cylinder shape and orientation was chosen to reduce the computational cost of the 3-D LES. Specifically, the mesh resolution required to numerically resolve the flow past the square cylinders is lower than that required for circular cylinders. This allowed us to conduct full 3-D sufficiently well-resolved simulations in long channels at relatively high Reynolds numbers. The main effect of using square cylinders is that for the same flow and geometrical conditions ( $h_{1}/H,D/H,\unicode[STIX]{x1D719}$ ) square cylinders produce somewhat higher total drag. The main effect of using horizontal cylinders instead of vertical cylinders is related to the orientation of the eddies shed in the separated shear layers of the cylinders. While horizontal cylinders shed primarily horizontal vorticity, which is oriented in the same direction as the baroclinic vorticity generated by the density differences and the vorticity generated in the boundary layer forming at the bottom edge of the porous layer, vertical cylinders shed primarily vertical vorticity. This is expected to generate some differences in the vortical interactions observed for vertical and horizontal cylinders, which should also have some effect on mixing. However, in both configurations the axes of the cylinders are oriented perpendicular to the flow generated by the advancing gravity current. For such cases, the cylinder drag coefficient is approximately the same, thus the capacity of the cylinders to retard the flow should be fairly similar for cases with identical ‘bulk’ descriptors of the porous medium ( $ah_{1}$ and $\unicode[STIX]{x1D719}$ ). This should result in a similar level of energy associated with the eddies generated by the cylinders inside the free surface current in the two cases.

To mimic in a more realistic way the relative positions of plant stems in a real patch of vegetation and to reduce the regularity of the wake-to-cylinder interactions inside the porous layer, a random 2-D displacement of approximately $0.5D$ – $1.0D$ was applied to the original staggered arrangement of the cylinders within the porous layer in the $x$ – $y$ plane. The number of solid cylinders placed within the porous layer varied between $N=118$ and $N=314$ (table 1). The surface of the square cylinders, the horizontal channel bottom and the lateral boundaries situated at $|x|=L$ were treated as no-slip smooth surfaces (zero velocity). The top boundary was treated as a free-slip boundary with zero normal velocity. This is standard numerical treatment for simulations in which the top boundary is characterized by fairly small deformations of the free surface during the duration of the flow (Rodi, Constantinescu & Stoesser Reference Rodi, Constantinescu and Stoesser2013). This assumption is consistent with what was observed in the experiments of Zhang & Nepf (Reference Zhang and Nepf2011). Before the lock gate was released at $t=0$ , the non-dimensional density was $C=1$ on the left side of the domain ( $x/H<0$ ) containing the heavier, colder fluid and $C=0$ on the right side containing lighter, warmer fluid. The surface-normal concentration gradient was set to zero at all no-slip and free-slip boundaries (i.e. no-flux boundary conditions for $C$ ). The molecular Prandtl number was equal to 7.

The Reynolds number in all but one of the simulations matched the value ( $Re=5758$ ) used in the experimental study of Zhang & Nepf (Reference Zhang and Nepf2011) conducted with $H=15$  cm and $u_{b}\approx 0.038~\text{m}~\text{s}^{-1}$ . This Reynolds number falls on the low end of the expected field conditions in fresh-water vegetated areas ( $5000<Re<80\,000$ corresponding to $H=0.5{-}2$  m, temperature differences of $1{-}2\,^{\circ }\text{C}$ and $u_{b}\approx 0.01{-}0.04~\text{m}~\text{s}^{-1}$ ). For deeper channels, the Reynolds number can be over 80 000. To investigate Reynolds number scale effects and to understand how representative the low Reynolds number laboratory results are to the full range of field conditions, an additional simulation ( $\unicode[STIX]{x1D719}=8\,\%$ , $h_{1}/H=0.28$ ) was performed at $Re=500\,000$ . This Reynolds number was sufficiently high that the current evolution was close to the inviscid limit and within the typical range of values ( $Re=10^{5}{-}10^{6}$ ) used to numerically study viscous effects on the evolution of lock-exchange gravity currents (e.g. see Ooi et al. Reference Ooi, Constantinescu and Weber2009; Ozan et al. Reference Ozan, Constantinescu and Hogg2015).

Table 1.

Main parameters of the simulations.

The ratio $D/H$ was the same (0.035) as the one used by Zhang & Nepf (Reference Zhang and Nepf2011). The length of the channel was $2L=36H$ , which is more than twice the value ( $2L=14H$ ) used in the experiments. Simulations were conducted with $\unicode[STIX]{x1D719}=8\,\%$ , 16 % and 28 % (see table 1), which fall within the range of values observed in the field ( $\unicode[STIX]{x1D719}=0.01{-}0.45$ , Zhang & Nepf Reference Zhang and Nepf2011; Downing-Kunz & Stacey Reference Downing-Kunz and Stacey2012). The range of values of the non-dimensional frontal area per unit volume, $aH$ , was 2.0–6.4, while $ah_{1}$ varied between 0.23 and 2.24. Finally, simulations were conducted with porous layer heights $h_{1}/H=0.1$ , 0.28 and 0.47. For field conditions, one expects $0.1<h_{1}/H<0.8$ (Zhang & Nepf Reference Zhang and Nepf2011).

The width of the computational domain was $W=H$ and the flow was assumed to be periodic in the spanwise direction. In the $Re=5758$ simulations, the grid contained $7200\times 192\times 48$ nodes in the streamwise, spanwise and vertical directions, respectively. In the simulation conducted with $Re=500\,000$ , the grid was finer $(9000\times 288\times 64)$ because of the need to sufficiently resolve the near-wall flow. The time step was $0.001t_{0}$ , where the non-dimensional time scale was $t_{0}=H/u_{b}$ .

Figure 2 compares the velocity profiles measured (red curves) in the experiments of Zhang & Nepf (Reference Zhang and Nepf2011) to those predicted (black curves) by the simulations for two cases. In the simulations, the velocity profiles were averaged over a streamwise distance of $H/2$ to remove the dependence of the velocity profile from the exact position of the cylinders within the porous layer in a certain vertical section. One should mention that the simulations were conducted to square horizontal cylinders and with a slightly lower solid volume fraction compared to the corresponding experiment conducted with circular vertical cylinders. In the case with a low porous layer thickness ( $h_{1}/H\approx 0.1$ ) the agreement is satisfactory inside the bottom current. Some disagreement is observed inside the intrusion layer, where the peak velocity is situated slightly closer to the top boundary in the simulation. The agreement between simulation and experiment is poorer in the case with a high porous layer thickness ( $h_{1}/H=0.28$ ). However, this simulation correctly predicts the value and location of the peak velocity in the intrusion layer. The disagreement between the simulation and experimental results in figure 2 can be partly attributed to the difference in cylinder geometry. The simulations we run were limited to rectangular cylinders, which produce a higher drag than the circular cylinders used in the experiments. We tried to compensate for this by choosing a simulation with a slightly lower solid volume fraction than the experiment so that the total drag force induced by the obstacles present inside the porous layer was as close as possible. However, this adjustment could not be exact, since the drag force was not known a priori, and this may now explain the disagreement between simulation and observation in figure 2.

Figure 2.

Comparison between velocity profiles predicted by the experiments of Zhang & Nepf (Reference Zhang and Nepf2011) and present simulations ( $Re=5758$ , square cylinders). Results in the figure are compared for two sets of cases: (1) LES: $\unicode[STIX]{x1D719}=8\,\%$ , $h_{1}/H=0.28$ versus exp.: $\unicode[STIX]{x1D719}=10\,\%$ , $h_{1}/H=0.28$ ; (2) LES: $\unicode[STIX]{x1D719}=8\,\%$ , $h_{1}/H=0.1$ versus exp.: $\unicode[STIX]{x1D719}=10\,\%$ , $h_{1}/H=0.13$ . The streamwise velocity profiles are non-dimensionalized with the discharge per unit width of the free-surface current in the cross-section, $q_{c}$ .

3 Structure of the free-surface gravity current

3.1 Density and vorticity fields

The two main geometrical parameters that determine the structure of the surface current are the solid volume fraction, $\unicode[STIX]{x1D719}$ , and the height of the porous layer, $h_{1}/H$ . Their influence is discussed based on simulations conducted with $Re=5758$ (table 1). First, consider the simulations conducted with a porous layer depth $h_{1}/H=0.28$ . In these cases the surface currents are thicker than the porous layer, and their thickness increases with increasing $\unicode[STIX]{x1D719}$ (e.g. from $0.6H$ for $\unicode[STIX]{x1D719}=8\,\%$ in figure 3 a to approximately $0.75H$ for $\unicode[STIX]{x1D719}=24\,\%$ in figure 3 b). This is because as $\unicode[STIX]{x1D719}$ increases, the drag force per unit length of porous layer increases, and more of the lighter fluid entering the left side of the channel is diverted into the intrusion layer close to $x/H=0$ . As no cylinders are present inside the intrusion layer ( $\unicode[STIX]{x1D719}=0\,\%$ ), the mean streamwise velocity within the intrusion layer is larger than the velocity inside the porous layer, at least away from the head region of the free-surface current (see also discussion of figure 5 a,c).

Figure 3.

Distribution of the non-dimensional density, $C$ , within and around the head of the surface current in the simulations with: (a) $\unicode[STIX]{x1D719}=8\,\%$ , $h_{1}/H=0.28$ , $Re=5758$ ; (b)  $\unicode[STIX]{x1D719}=24\,\%$ , $h_{1}/H=0.28$ , $Re=5758$ ; (c) $\unicode[STIX]{x1D719}=8\,\%$ , $h_{1}/H=0.1$ , $Re=5758$ ; (d) $\unicode[STIX]{x1D719}=8\,\%$ , $h_{1}/H=0.47$ , $Re=5758$ ; (e) $\unicode[STIX]{x1D719}=8\,\%$ , $h_{1}/H=0.28$ , $Re=500\,000$ . The solid line shows the $C=0.99$ isocontour. Results are compared when $|x_{f}/H|\approx 15$ . The arrows point toward intrusions of lighter fluid into heavier fluid occurring at the front generated by flow acceleration as lighter fluid from behind the front is pushed into the space between two neighbouring cylinders.

Figure 4.

Distribution of the spanwise vorticity magnitude, $\unicode[STIX]{x1D714}_{z}(H/u_{b})$ , within and around the head of the surface current in the simulations with: (a) $\unicode[STIX]{x1D719}=8\,\%$ , $h_{1}/H=0.28$ , $Re=5758$ ; (b) $\unicode[STIX]{x1D719}=24\,\%$ , $h_{1}/H=0.28$ , $Re=5758$ ; (c) $\unicode[STIX]{x1D719}=8\,\%$ , $h_{1}/H=0.1$ , $Re=5758$ ; (d)  $\unicode[STIX]{x1D719}=8\,\%$ , $h_{1}/H=0.47$ , $Re=5758$ ; (e) $\unicode[STIX]{x1D719}=8\,\%$ , $h_{1}/H=0.28$ , $Re=500\,000$ . The solid line shows the $C=0.99$ isocontour. Results are compared when $|x_{f}/H|\approx 15$ .

Depending on the porous layer solid volume fraction, two scenarios are possible near the front of the current. First, for high $\unicode[STIX]{x1D719}$ values, the front advances faster in the intrusion layer than in the porous layer (see supplementary movie 1 available athttps://doi.org/10.1017/jfm.2016.698). This case is illustrated in figure 3(b) for $\unicode[STIX]{x1D719}=24\,\%$ and also occurs for $\unicode[STIX]{x1D719}=16\,\%$ (not shown in figure 3). Specifically, at the time shown in figure 3(b), the warm-water front has only advanced to $x/H=-13.8$ in the porous layer, but has advanced to $x/H=-14.8$ in the intrusion layer. As both fronts progress, the distance between them remains bounded (never larger than $4h_{1}$ ). This is because an unstable stratification is generated at the head of the intrusion layer as the lighter fluid within the porous layer advances ahead of the heavier fluid in the porous layer, carrying lighter fluid beneath a region of heavier fluid within the porous layer. As a result of this unstable density gradient (see region with $13.8<|x/H|<14.5$ in figure 3 b), lighter fluid from the intrusion layer penetrates the porous layer and mixes with the slower advancing fluid inside the porous layer. The non-dimensional density and spanwise vorticity magnitude contours in figures 3(b) and 4(b) show that intrusions of mixed fluid with a density close to $C\approx 0.5$ occur over a distance of about $2H(7h_{1})$ behind the front in the porous layer ( $11.8<|x/H|<13.8$ ). The jet-like intrusions correspond to regions of high vorticity magnitude oriented upwards, starting at the bottom of the porous layer in figure 4(b) (e.g. from $|x/H|=11.8$ until close to the front). So, the mass exchange between the porous and intrusion layers is not restricted to the region between the streamwise positions of the front of the free-surface current in the two layers (e.g. the front is close to $|x/H|=14.8$ inside the intrusion layer and $|x/H|=13.8$ inside the porous layer in figure 3 b).

Figure 5.

Vertical profiles of the streamwise velocity, $u/u_{b}$ , at different streamwise locations. (a) $\unicode[STIX]{x1D719}=24\,\%$ , $h_{1}/H=0.28$ , $Re=5758$ ; (b) $\unicode[STIX]{x1D719}=8\,\%$ , $h_{1}/H=0.47$ , $Re=5758$ ; (c) $\unicode[STIX]{x1D719}=8\,\%$ , $h_{1}/H=0.28$ , $Re=5758$ ; (d) $Re=5758$ versus $Re=500\,000$ for $\unicode[STIX]{x1D719}=8\,\%$ , $h_{1}/H=0.28$ . Results are compared when $|x_{f}|/H\approx 18$ . The height of the intrusion layer, $h_{2}$ , is determined by the bottom of the porous layer and the vertical location where $u/u_{b}=0$ .

Second, for low $\unicode[STIX]{x1D719}$ values, the fronts in the porous and intrusion layers advance at the same rate (e.g. $\unicode[STIX]{x1D719}=8\,\%$ in figure 3 a and movie 1). The velocity difference between the porous and intrusion layers is still high (e.g. see figure 5 c), but vertical exchange of fluid occurs more rapidly, keeping the two fronts aligned. The impact of this vertical mixing is evident in the uniform density (uniform colour) between the porous and intrusion layers in figure 3(a). In this case, most of the fluid present close to the front of the surface current in the porous layer originates in the porous layer. This can be inferred from the vorticity magnitude fields that show that the separated shear layers (SSLs) forming behind the individual cylinders are close to parallel to the streamwise direction, even for the cylinders situated close to the front (e.g. see figure 4 a,c). The detachment of energetic eddies from the SSLs is the main mechanism for mixing within the porous layer in the cases with low $\unicode[STIX]{x1D719}$ values. This mixing mechanism is not present in the cases with high $\unicode[STIX]{x1D719}$ values, for which no strong unsteady SSLs form on the sides of the individual cylinders (e.g. figure 4 b). The main reason is that for high $\unicode[STIX]{x1D719}$ values, most of the flow leaves the porous layer close to the lock gate position due to the large drag induced by the cylinders and re-enters the porous layer close to the front (figure 1 b). This is confirmed in a quantitative way by examining the streamwise variation of the vertical flux of fluid entering/leaving the porous layer through its bottom boundary, $Q_{v}(x)$ , shown in figure 6 (see detailed discussion in § 3.2). Similar to cases with high $\unicode[STIX]{x1D719}$ values, lighter fluid from the intrusion layer is also advected into the porous layer in the cases with low $\unicode[STIX]{x1D719}$ values. This can be inferred by looking at the orientation of the SSLs, which is parallel to the ‘mean’ incoming flow approaching each cylinder. Moreover, the vorticity magnitude inside the SSLs of each cylinder is proportional to the ‘mean’ velocity magnitude in the incoming flow approaching that cylinder. In the cases with low $\unicode[STIX]{x1D719}$ values, the orientation of the SSLs is generally not horizontal for the cylinders situated at the bottom of the porous layer (e.g. see figure 4 a, $11<|x/H|<14$ , where the SSLs are oriented at an angle of $0^{\circ }$ – $30^{\circ }$ with respect to the horizontal direction and pointing toward the free surface, which indicates flow entering the porous layer).

Figure 6.

Volumetric streamwise flux discharge within the porous layer, $Q_{t}(x)/Q_{0}$ , net vertical flux at the interface between the top and intrusion layers, $Q_{v}(x)/Q_{0}$ , and total streamwise flux within the free-surface current, $Q_{ti}(x)/Q_{0}$ ( $Q_{0}=u_{b}HW$ ) (a) $\unicode[STIX]{x1D719}=24\,\%$ , $h_{1}/H=0.28$ , $Re=5758$ ; (b) $\unicode[STIX]{x1D719}=8\,\%$ , $h_{1}/H=0.47$ , $Re=5758$ ; (c) $\unicode[STIX]{x1D719}=8\,\%$ , $h_{1}/H=0.28$ , $Re=5758$ . Results are compared when $|x_{f}|/H\approx 18$ .

The influence of the porous layer depth, $h_{1}$ , on the flow structure was investigated for a relatively low solid volume fraction ( $\unicode[STIX]{x1D719}=8\,\%$ , see movie 2). Figure 3(a,c,d) compares the flow structure among cases with different $h_{1}/H$ ( $=$ 0.28, 0.1 and 0.47, respectively) when the front is situated close to $|x/H|=14.8$ . For $h_{1}/H=0.1$ and 0.28, the front has a fairly uniform distribution over depth, i.e. the front has a blunt nose. However, the structure of the surface current front changes for higher $h_{1}/H$ values, for which most of the current in contained within the porous layer, and the height of the intrusion layer, $h_{2}$ , is small ( $h_{2}\ll h_{1}$ ). For such cases (e.g. $h_{1}/H=0.47$ in figure 3 d), the front position varies linearly over most of the porous layer. This is similar to the shape observed for currents propagating in a fully vegetated channel (Tanino et al. Reference Tanino, Nepf and Kulis2005; Ozan et al. Reference Ozan, Constantinescu and Hogg2015). Because of the linear distribution of the front, for cases with fairly low $\unicode[STIX]{x1D719}$ and high $h_{1}/H$ , the front in the intrusion layer is situated behind the most advanced part of the front within the porous layer. The strength of the SSLs and the coherence of the eddies shed from the cylinders both decrease with increasing $h_{1}/H$ (compare figure 4 a,d). For example, for $h_{1}/H=0.47$ (figure 4 d), the SSLs are too weak to generate large-scale wake eddies even for the cylinders situated close to the front. While the total streamwise flux discharge within the porous and intrusion layers decays significantly with increasing $h_{1}/H$ , the streamwise flux discharge within the porous layer, $Q_{t}$ , is much less dependent on $h_{1}/H$ (e.g. see discussion of figure 6 b,c in § 3.2). This means that the mean streamwise velocity within the porous layer, $Q_{t}/h_{1}$ , decreases with increasing $h_{1}/H$ , which explains the aforementioned reduction in the strength of the SSLs with increasing $h_{1}/H$ .

Increasing the Reynolds number from $Re=5758$ to $Re=500\,000$ decreases the size of the large-scale eddies, as shown by the finer-scale variation in vorticity in figure 4(e) (high $Re$ ) compared to figure 4(a) (low $Re$ ). In addition, the strength of the wake-to-cylinder interactions was greater for the higher $Re$ case (see movie 3). Both these effects are a result of vortex stretching phenomena becoming more intense with the increase in the Reynolds number. Another consequence of increasing the Reynolds number is greater mixing, which results is a more diffuse density interface (see figure 3 a,e) and a slight increase of the height of the surface current (see also figure 5 d, where the bottom of the free-surface current corresponds to the vertical location where the streamwise velocity is equal to zero).

4 Mixing induced by free-surface current

The presence of cylinders in the channel has a large influence on entrainment and mixing associated with the propagation of the surface current. Strong mixing occurs in the region situated close to the front and in regions where energetic interfacial Kelvin–Helmholtz (KH) billows are present. Though the discussion is this section strictly applies to the case of channels with horizontal cylinders, it is expected that the effects of varying the main non-dimensional geometrical parameters on mixing will be qualitatively similar for channels with vertical and horizontal cylinders.

Figure 7.

Evolution of the normalized volume of mixed fluid, $V^{\prime }$ , with the front position, $|x_{f}|/H$ . The volume $V^{\prime }$ is defined as the fluid with $0.01<C<0.99$ inside the left side of the channel ( $x/H<0$ ). (a) Effect of the solid volume fraction in the $Re=5758$ , $h_{1}/H=0.28$ simulations with $\unicode[STIX]{x1D719}=8\,\%$ (light blue line), $\unicode[STIX]{x1D719}=16\,\%$ (blue line) and $\unicode[STIX]{x1D719}=24\,\%$ (black line); (b) effect of $h_{1}/H$ in the ( $Re=5758$ , $\unicode[STIX]{x1D719}=8\,\%$ ) simulations with $h_{1}/H=0.1$ (orange line), $h_{1}/H=0.28$ (pink line) and $h_{1}/H=0.47$ (green line); (c) effect of Reynolds number in the ( $\unicode[STIX]{x1D719}=8\,\%$ , $h_{1}/H=0.28$ ) simulations with $Re=5758$ (light blue line), $Re=500\,000$ (red line).

A way to quantify mixing and entrainment in a flow that initially contains two unmixed regions (e.g. Necker et al. Reference Necker, Härtel, Kleiser and Meiburg2005) is to estimate the variation with time of the volume containing fluid with a density falling between the initial densities of the lighter and heavier fluids. After non-dimensionalizing with the initial volume of lighter fluid, $V_{0}$ , the normalized volume containing mixed fluid is $V^{\prime }(x_{f}/H)=(1/V_{0})\int \unicode[STIX]{x1D6FE}\,\text{d}V$ , in which $\unicode[STIX]{x1D6FE}=1$ if $0.01<C<0.99$ and $\unicode[STIX]{x1D6FE}=0$ otherwise, which defines mixed fluid as fluid for which the non-dimensional density differs from the initial values by at least 1 %. The integration is done over the domain $x/H<0$ and only over the part of the domain containing fluid. Initially, the channel contains only unmixed fluid with $C=0$ or $C=1$ , so $V^{\prime }(0)=0$ .

The temporal evolution of the normalized volume of mixed fluid, $V^{\prime }$ , is a function of $\unicode[STIX]{x1D719}$ (figure 7), because $\unicode[STIX]{x1D719}$ determines the relative strength of the main mixing mechanisms. For example, consider the $Re=5758$ simulations with $h_{1}/H=0.28$ (figure 7 a). As $\unicode[STIX]{x1D719}$ increases, there is a decrease in the coherence of the KH billows along the bottom of the intrusion layer (figure 4 a,b). In addition, for higher $\unicode[STIX]{x1D719}$ there is a longer region of unstably stratified fluid in between the streamwise positions of the front of the free-surface current within the porous and intrusion layers (figure 3 b). As $\unicode[STIX]{x1D719}$ increases, the volume of mixed fluid at the head is larger (e.g. compare the spatial extent of regions with $0.2<C<0.8$ in figure 3 a,b showing the head region; see also figure 8 a,b). This increase in mixing is mainly because of intrusions of lighter fluid penetrating into the porous layer (see figures 3 b and 4 b for $11.8<|x/H|<14.5$ ). Moreover, animations of the non-dimensional density fields show that the interaction of the front with the cylinders increases the volume of mixed fluid via engulfment of heavier fluid from behind the cylinders by the advancing front of the surface current (see movie 1). Cases with a higher $\unicode[STIX]{x1D719}$ have a proportionally larger number of cylinders per unit streamwise length, because the cylinder size is constant. Assuming that the amount of higher-density fluid engulfed behind each cylinder is independent of $\unicode[STIX]{x1D719}$ , this mixing mechanism will become more important as $\unicode[STIX]{x1D719}$ increases further. The changes in mixing regime with increasing $\unicode[STIX]{x1D719}$ can now be used to explain the dependence of $V^{\prime }$ on $\unicode[STIX]{x1D719}$ in figure 7(a).

Figure 8.

Distribution of the non-dimensional density, $C$ , in the left side of the channel $(x/H<0)$ showing regions containing mixed fluid in the simulations with: (a) $\unicode[STIX]{x1D719}=8\,\%$ , $h_{1}/H=0.28$ , $Re=5758$ ; (b) $\unicode[STIX]{x1D719}=24\,\%$ , $h_{1}/H=0.28$ , $Re=5758$ ; (c) $\unicode[STIX]{x1D719}=8\,\%$ , $h_{1}/H=0.1$ , $Re=5758$ ; (d) $\unicode[STIX]{x1D719}=8\,\%$ , $h_{1}/H=0.47$ , $Re=5758$ ; (e) $\unicode[STIX]{x1D719}=8\,\%$ , $h_{1}/H=0.28$ , $Re=500\,000$ . Results are compared when $|x_{f}|/H\approx 15$ . Regions containing fluid with $C>0.99$ and $C<0.01$ were blanked out.

Figure 9.

Front position as a function of time for the surface current (UL) and the bottom-propagating current (LL). (a) Effect of the solid volume fraction in the $Re=5758$ , $h_{1}/H=0.28$ simulations with $\unicode[STIX]{x1D719}=8\,\%$ (light blue line), $\unicode[STIX]{x1D719}=16\,\%$ (blue line) and $\unicode[STIX]{x1D719}=24\,\%$ (black line); (b) effect of $h_{1}/H$ in the ( $Re=5758$ , $\unicode[STIX]{x1D719}=8\,\%$ ) simulations with $h_{1}/H=0.1$ (orange line), $h_{1}/H=0.28$ (light blue line) and $h_{1}/H=0.47$ (green line); (c) effect of Reynolds number in the ( $\unicode[STIX]{x1D719}=8\,\%$ , $h_{1}/H=0.28$ ) simulations with $Re=5758$ (light blue line), $Re=500\,000$ (red line).The black dashed straight lines show the range over which $x_{f}/H\sim (t/t_{0})^{\unicode[STIX]{x1D6FC}}$ , where $\unicode[STIX]{x1D6FC}$ is the power law coefficient, also given in the figure.

The rate of growth of the normalized volume of mixed fluid during the initial stages of the propagation of the surface current is fairly independent of $\unicode[STIX]{x1D719}$ for $\unicode[STIX]{x1D719}>16\,\%$ (figure 7 a), because of the following compensating effects. As $\unicode[STIX]{x1D719}$ increases, the interfacial KH billows become weaker, such that mixing associated with these billows is reduced. However, as $\unicode[STIX]{x1D719}$ increases the set-up of unstable stratification at the front of the free-surface current and associated mixing increases. In the three cases compared in figure 7(a), the free-surface current enters the drag-dominated regime ( $x_{f}\sim t^{0.75}$ ) when $|x_{f}/H|\approx 6$ for $\unicode[STIX]{x1D719}=8\,\%$ and $|x_{f}/H|\approx 7{-}8$ for $\unicode[STIX]{x1D719}=16\,\%$ and 24 % (see discussion of figure 9 a). That the transition to the drag-dominated regime occurs at a greater distance with higher $\unicode[STIX]{x1D719}$ is somewhat counter-intuitive, i.e. with a higher $\unicode[STIX]{x1D719}$ one expects a higher drag in the porous layer and a more rapid transition to a drag-dominated flow. Indeed, this occurs in a fully vegetated channel for which the flow is largely one-dimensional (figure 1 a), and for which the transition distance is proportional to $\unicode[STIX]{x1D719}^{-1}$ (Tanino et al. Reference Tanino, Nepf and Kulis2005). However, in a partially vegetated channel (figure 1 b), the flow is strongly two-dimensional, especially at high $\unicode[STIX]{x1D719}$ . Most of the flow entering the porous layer at $x=0$ quickly leaves it (blue arrows in figure 1 b) and re-enters the porous layer close to the front (green arrows in figure 1 b). As a result, the velocity within the porous layer is negligible over most of the length of the current (figure 6 a), and the advance of the front in the porous layer is not associated with the horizontal current in this layer. Instead, the advance of the front in the porous layer is associated with the infiltration of lighter fluid from the intrusion layer beneath it (figure 1 b). Because the lighter fluid moves faster in the intrusion layer, a zone of unstable stratification is created that drives the vertical migration of lighter fluid from the intrusion layer into the porous layer (figure 3 b). It is this migration of lighter fluid from the intrusion layer that appears as the lighter fluid ‘advancing’ inside the porous layer. As $\unicode[STIX]{x1D719}$ increases, this exchange slows, and it takes more time and greater streamwise distance for the advance of the front in the porous layer to reach an equilibrium with the advance of the front in the intrusion layer. The transition to the drag-dominated regime occurs when the distance between the front in the porous layer and the front in the intrusion layer is close to constant. The time and distance travelled by the front of the free surface current inside the intrusion layer until the transition to the drag-dominated regime occurs increases with increasing $\unicode[STIX]{x1D719}$ . Once the transition to the drag-dominated regime occurs, the rate of growth of $V^{\prime }$ decreases. As the transition between the fast and the slow growth regimes occurs at a larger $|x_{f}/H|$ as $\unicode[STIX]{x1D719}$ increases, $V^{\prime }$ nearly doubles its value at $|x_{f}/H|=15$ between $\unicode[STIX]{x1D719}=8\,\%$ and 24 %.

The effect of increasing $h_{1}/H$ for the cases with $Re=5758$ and $\unicode[STIX]{x1D719}=8\,\%$ is a slight monotonic increase of $V^{\prime }$ with decreasing $h_{1}/H$ (at a given $|x_{f}/H|$ ) during the later stages of propagation of the current (figure 7 b). Examination of the non-dimensional density contours in figure 8(a,c,d) shows that the contribution to $V^{\prime }$ due to mixing by the KH billows is monotonically decreasing with increasing $h_{1}/H$ . Meanwhile, a larger part of the surface current is contained within the porous layer. Mixing induced by the interaction of the head with the cylinders grows monotonically with the increase in $h_{1}/H$ despite the fact that the wakes of the cylinders are more energetic in cases with a low $h_{1}/H$ (see figure 4 a,c,d).

Comparison of the non-dimensional density fields show that the effect of increasing the Reynolds number in the case with $\unicode[STIX]{x1D719}=8\,\%$ and $h_{1}/H=0.28$ is felt starting with $|x_{f}/H|>4$ (figure 7 c). For small $|x_{f}/H|$ , the coherence of the KH billows is slightly higher in the lower Reynolds number simulation, while mixing behind the front is higher in the higher Reynolds number simulation (see movie 3). For large $|x_{f}/H|$ , mixing generated near the front becomes significantly larger in the higher Reynolds number simulation (figure 8 a,e). This is because the size and penetration length of intrusions of lower-density fluid into the low-velocity higher-density fluid situated close to the front of the surface current increase with increasing $Re$ (see arrows at the front of the free-surface current pointing toward such intrusions in figure 3 a,e). The intrusions occur as lighter fluid moving with velocities close to the front velocity is locally accelerated as it is pushed in the opening between two cylinders. As a result, it can penetrate into a region of higher-density fluid and create unstable stratification. Moreover, most of the increase of the normalized volume of mixed fluid, $V^{\prime }$ , is due to the sharp increase in the number of highly 3-D smaller eddies in the wakes of the cylinders and within the intrusion layer (figure 4 a,e). These eddies are very effective in promoting mixing.

5 Front velocity

For comparison, we quickly review the fully vegetated case, for which the evolutions of the surface and bottom-propagating currents are close to anti-symmetric, and for sufficiently large values of $\unicode[STIX]{x1D719}$ and $Re$ , the two currents transition to a quadratic-drag-dominated regime where $x_{f}/H\sim (t/t_{0})^{3/4}$ (see Ozan et al. Reference Ozan, Constantinescu and Hogg2015). The value of the exponent ( $\unicode[STIX]{x1D6FC}=3/4$ ) is close, but not identical, to the shallow-water theory predictions ( $\unicode[STIX]{x1D6FC}=2/3$ ). The difference in the temporal evolution of the front position was attributed to the mixing occurring close to the front that is not accounted for by theory. Moreover, during the same regime the non-dimensional discharge of the top current at the lock gate, $q(x=0,t)/q_{0}$ , varies proportionally to $(t/t_{0})^{-1/3}$ , in excellent agreement with theory ( $q_{0}=u_{b}WH/2$ ).

In this study, the front position is defined using the non-dimensional concentration contours, specifically the streamwise position at which the spanwise-averaged value of $C$ was equal to 0.99 (surface current) or 0.01 (bottom current). In the present configuration, the bottom current advances in a region with $\unicode[STIX]{x1D719}=0\,\%$ , for which one might expect that the front velocity would remain constant with time. In contrast, the surface current advances, at least in part, within the porous layer. Thus, for sufficiently large $\unicode[STIX]{x1D719}$ and/or $x_{f}$ , one expects that the surface current will start decelerating due to the porous layer drag. These trends are confirmed by the temporal variations of the front position plotted in figure 9 in log–log scale. Because of continuity, the decay of the surface current discharge at $x/H=0$ should eventually result in a decay of the bottom current front velocity as well.

For $h_{1}/H=0.28$ and $\unicode[STIX]{x1D719}\geqslant 8\,\%$ (figure 9 a), the gravity current transitions to a quadratic-drag-dominated regime with $x_{f}/H\sim (t/t_{0})^{3/4}$ . Despite the relatively low value of the Reynolds number, the surface current does not transition to the linear drag-dominated regime within the simulation time. Meanwhile, the discharge at $x/H=0$ varies proportionally to $(t/t_{0})^{-1/3}$ (figure 10). This behaviour is identical to that observed in the fully vegetated case (Ozan et al. Reference Ozan, Constantinescu and Hogg2015). One suspects again that mixing at the front, which is stronger than that in the fully vegetated case, is responsible for the larger rate of increase of $|x_{f}|$ with time with respect to the self-similar theoretical solution. At any given time, the distance between the lock gate and the front of the surface current decreases with the increase in $\unicode[STIX]{x1D719}$ , though the difference between the three cases is small. In all three simulations the bottom current transitions to a slumping phase with a front velocity of $0.41u_{b}$ . While the slumping phase is observed until the end of the simulation in the $\unicode[STIX]{x1D719}=8\,\%$ case (linear variation of $x_{f}/H$ with $t/t_{0}$ ), results for the $\unicode[STIX]{x1D719}=16\,\%$ and $\unicode[STIX]{x1D719}=24\,\%$ cases show that the front velocity of the bottom current starts decaying slightly with time for $|x_{f}/H|>10$ (e.g. in the $\unicode[STIX]{x1D719}=24\,\%$ case, the front velocity decays from $0.41u_{b}$ at $|x_{f}/H|=10$ to $0.36u_{b}$ at $|x_{f}/H|=16$ ).

Figure 10.

Temporal evolution of the discharge at the lock gate position, $q/q_{0}(q_{0}=u_{b}WH/2)$ , for the free-surface current in the simulations with $\unicode[STIX]{x1D719}=8\,\%$ , $h_{1}/H=0.28$ , $Re=5758$ (light blue line), $\unicode[STIX]{x1D719}=16\,\%$ , $h_{1}/H=0.28$ , $Re=5758$ (blue line), $\unicode[STIX]{x1D719}=24\,\%$ , $h_{1}/H=0.28$ , $Re=5758$ (black line), $\unicode[STIX]{x1D719}=8\,\%$ , $h_{1}/H=0.1$ , $Re=5758$ (orange line) and $\unicode[STIX]{x1D719}=8\,\%$ , $h_{1}/H=0.47$ , $Re=5758$ (green line). The black dashed straight lines show the range over which $q/q_{0}\sim (t/t_{0})^{\unicode[STIX]{x1D6FE}}$ , where $\unicode[STIX]{x1D6FE}$ is the power law coefficient also given in the figure.

For the same $\unicode[STIX]{x1D719}$ , the advance of the surface current differs with $h_{1}/H$ during the later stages of propagation. For sufficiently high values of $h_{1}/H$ (e.g. 0.28 and 0.47, see figure 9 b) the surface current transitions to the drag-dominated regime, with $\unicode[STIX]{x1D6FC}=0.75$ . The transition happens sooner with larger $h_{1}/H$ . In contrast, for $h_{1}/H=0.1$ the front velocity starts to decay with time, but does not reach the $\unicode[STIX]{x1D6FC}=3/4$ regime. A best fit gives $\unicode[STIX]{x1D6FC}=0.95$ . For small values of $h_{1}/H$ , most of the surface current is contained within the intrusion layer and the cylinders within the porous layer act as large-scale roughness on the top of the current. Most probably, inertia and drag are both important and comparable to the buoyancy forces, such that one expects a value of $\unicode[STIX]{x1D6FC}$ between 1, corresponding to no additional drag forces, i.e. $\unicode[STIX]{x1D719}=0\,\%$ , and $3/4$ , corresponding to porous layer drag forces much larger than the inertial forces. It is also likely that in a longer channel one will observe a gradual decrease of $\unicode[STIX]{x1D6FC}$ toward $3/4$ , even for $h_{1}/H=0.1$ . Consistent with the variation of $\unicode[STIX]{x1D6FC}$ in the other cases in which $h_{1}$ was not much smaller than the height of the intrusion layer, $h_{2}$ , the discharge at $x/H=0$ varies proportionally to $(t/t_{0})^{-1/3}$ in the simulations with $h_{1}/H=0.28$ and 0.47. However, $q/q_{0}\approx (t/t_{0})^{-1/5}$ in the simulation with $h_{1}/H=0.1$ (figure 10). In all three simulations, the bottom-propagating current transitions to a slumping phase, but the constant front velocity over this regime decreases slightly with increasing $h_{1}/H$ . In the simulation with $h_{1}/H=0.47$ , the front velocity of the bottom current starts decaying with time for $|x_{f}/H|>12$ , which means the current leaves the slumping phase due to reduced discharge at $x/H=0$ .

The main effect of increasing the Reynolds number is to slightly increase the front velocity (e.g. by about 8–10 %). This is due to the increase of the front velocity with the Reynolds number for the part of the front situated within the intrusion layer ( $\unicode[STIX]{x1D719}=0\,\%$ ). However, once the transition to the drag-dominated regime is completed and drag forces largely dominate inertia forces, the difference between the front positions in the two cases remains close to constant in time (figure 9 c). For both cases, $\unicode[STIX]{x1D6FC}=3/4$ . The front velocity of the bottom current during the slumping phase increases from $0.41u_{b}$ to $0.46u_{b}$ , which is close to the inviscid limit ( $0.5u_{b}$ ) expected for lock-exchange currents propagating over a flat surface.

6 Flushing induced by the surface current

Convective flow driven by differences in water temperature plays an important role in the transport of nutrients and other substances between regions with and without vegetation (James & Barko Reference James and Barko1991; James, Barko & Eakin Reference James, Barko and Eakin1994; Zhang & Nepf Reference Zhang and Nepf2008). For floating vegetation, the residence time within the root (porous) layer is a critical parameter that can control the rate of nutrient uptake by roots, and the growth rate of floating vegetation. For example, when the biomass of hyacinth gets too high, which would be associated with a denser root zone, the growth rate declines (Imaoka & Teranishi Reference Imaoka and Teranishi1988). This might be explained by a decreased flushing rate associated with the denser root zone, which would limit the supply of nutrients. As the root zone becomes denser, the flushing rate declines, increasing the root zone residence time. If the root zone residence time becomes much longer that the uptake time scale, the nutrient uptake becomes supply limited, and overall uptake rate to the plant declines. Similarly, for water quality applications, the residence time of water within the porous (root) layer may control the overall nutrient removal from the water by both the roots and the biofilms existing on the root surfaces.

The residence time within the porous layer can be inferred from the time history of the volume of water with $C\leqslant 0.99$ advancing through the porous layer, $V$ . As the present investigation only considered cases in which the domain width, $W$ , and porous layer height, $h_{1}$ , were constant, we consider the variable $V/(Wh_{1})$ , which has the dimension of a length, or its non-dimensional equivalent $V/(HWh_{1})$ . Its slope characterizes the rate of streamwise advance (velocity) of lighter water through the porous layer, e.g. the rate at which new (lighter) water passes through the root layer (figure 11). Importantly, this quantity allows a more accurate estimation of the flushing time than the advection inferred from the velocity field within the porous layer, which would yield very long residence time, due to the near zero mean velocity in the porous layer in many cases (figure 5).

Figure 11.

Evolution of the normalized volume of heavier fluid removed by the surface current from the porous layer, $V/(h_{1}W)$ . The volume $V$ corresponds to the region with $C<0.99$ inside the porous layer ( $x/H<0$ , $1-h_{1}<y/H<1$ ). (a) Effect of solid volume fraction in the $Re=5758$ , $h_{1}/H=0.28$ simulations with $\unicode[STIX]{x1D719}=8\,\%$ (light blue line), $\unicode[STIX]{x1D719}=16\,\%$ (blue line) and $\unicode[STIX]{x1D719}=24\,\%$ (black line); (b) effect of $h_{1}/H$ in the ( $Re=5758$ , $\unicode[STIX]{x1D719}=8\,\%$ ) simulations with $h_{1}/H=0.1$ (orange line), $h_{1}/H=0.28$ (light blue line) and $h_{1}/H=0.47$ (green line); (c) effect of Reynolds number in the ( $\unicode[STIX]{x1D719}=8\,\%$ , $h_{1}/H=0.28$ ) simulations with $Re=5758$ (light blue line), $Re=500\,000$ (red line).

The rate of increase of the normalized volume of lighter water advancing through the porous layer, $V/(HWh_{1})$ , i.e. the flushing rate, increases with decreasing $\unicode[STIX]{x1D719}$ (figure 11 a) and decreasing $h_{1}/H$ (figure 11 b). For relatively small values of $\unicode[STIX]{x1D719}$ and sufficiently low values of $h_{1}/H$ , $V/(HWh_{1})$ increases at close to linear rate with time (constant flushing rate), with $V\sim t^{\unicode[STIX]{x1D6FD}}$ , and $\unicode[STIX]{x1D6FD}=0.95$ for the $h_{1}/H=0.1$ , $\unicode[STIX]{x1D719}=8\,\%$ case in figure 11(b). This result is not surprising given the close to linear increase of the front position with time ( $\unicode[STIX]{x1D6FC}=0.95$ ) for cases with a very thin porous layer ( $h_{1}/H=0.1$ , $\unicode[STIX]{x1D719}=8\,\%$ case in figure 9 b). For all other cases, the increase of $V/(HWh_{1})$ with time is less than linear and $\unicode[STIX]{x1D6FD}$ is slightly greater than $\unicode[STIX]{x1D6FC}$ ( $=$ 0.75). Only for the case where $h_{1}/H$ is sufficiently high ( $h_{1}/H=0.47~\unicode[STIX]{x1D719}=8\,\%$ ), such that most of the surface current is contained within the porous layer, $\unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D6FC}=0.75$ . In the case with the largest $\unicode[STIX]{x1D719}$ (e.g. $\unicode[STIX]{x1D719}=24\,\%$ , $h_{1}/H=0.28$ case in figure 11 a), the temporal variation of the rate of increase of $V/(HWh_{1})$ with time is non-monotonic. This is mainly because of the strong temporal variation in the distance between the front of the free-surface current in the porous layer and the front of the same current in the intrusion layer. The intrusions of lighter fluid into the porous layer in between the positions of the front of the free-surface current in the porous layer and in the intrusion layer are also contributing to the non-monotonic variation of the rate of increase of $V/(HWh_{1})$ . During the time intervals when the distance between the front positions in the two layers is relatively large, the rate of increase of $V/(HWh_{1})$ is high ( $\unicode[STIX]{x1D6FD}\approx 1$ ). Finally, the effect of increasing the Reynolds number is felt only after the end of the transition to the drag-dominated regime (figure 11 c). As expected, $V/(HWh_{1})$ increases slightly faster in the higher Reynolds number case.

As the discussion of the structure of the flow in §§ 3.1 and 3.2 has shown, the mean streamwise velocity inside the porous layer is small, suggesting that the renewal of water in the porous (root) layer would be minimal, especially for porous layers with large $\unicode[STIX]{x1D719}$ . However, the numerical results clearly show that the intrusion of less dense water from beneath the porous layer leads to significant convective mixing of fluid between the porous layer and the intrusion layer. This observation reveals the important, and previously unidentified, role of convective mixing and associated vertical exchange in setting the residence time of the root/porous layer. For example, for $h_{1}/H=0.28$ and $\unicode[STIX]{x1D719}=8\,\%$ the residence time estimated from the velocity within the root layer is 6 times the actual flushing time estimated through the concentration fields (figure 11). This discrepancy increases with increasing solid volume fraction, such that for $\unicode[STIX]{x1D719}=24\,\%$ the discrepancy is two orders of magnitude.

The biological uptake by roots and biofilms takes place only in the root/porous layer. Using information from figure 11 on the variation of the normalized volume of lighter water advancing through the porous layer, $V/(HWh_{1})$ , before the end of the transition to the drag-dominated regime ( $t/t_{0}<30$ for most cases, see figures 9 and 10) together with the value of $\unicode[STIX]{x1D6FD}$ one can estimate this variable at any time. This allows calculation of the length of the vegetated region that can be flushed in a typical day in real systems, $L_{d}$ . Ultsch (Reference Ultsch1973) found that the average temperature difference between floating hyacinth and open water is about $1^{\circ }\text{C}$ and persists for about 6 h each day (see figure 1–4 in Ultsch (Reference Ultsch1973). This temperature difference corresponds to a density difference of $0.25~\text{Kg}~\text{m}^{-3}$ . For typical applications, the mean water depth in the pond, $H$ , is about 1 m, the buoyancy velocity $u_{b}$ is $0.05~\text{m}~\text{s}^{-1}$ and the time scale $t_{0}$ ( $=H/u_{b}$ ) is 20 s. If one assumes that the relative height of the root layer is $h_{1}/H=0.28$ , then the value of $V/(HWh_{1})$ after 6 h ( $t/t_{0}\approx 1100$ ) is 95 for $\unicode[STIX]{x1D719}=8\,\%$ and 60 for $\unicode[STIX]{x1D719}=24\,\%$ . These values were obtained using the values of $V/(HWh_{1})$ from the corresponding curves in figure 11(a) at $t/t_{0}=60$ and then the assuming a power law variation $(V\sim t^{\unicode[STIX]{x1D6FD}})$ with the estimated values of $\unicode[STIX]{x1D6FD}$ ( ${\approx}$ 0.78 for both $\unicode[STIX]{x1D719}=8\,\%$ and $\unicode[STIX]{x1D719}=24\,\%$ cases with $h_{1}/H=0.28$ ) for $60<t/t_{0}<1100$ . Then $L_{d}=V(t=1,100t_{0}=6~\text{h})/(Wh_{1})$ . Using $H=1$  m, one obtains $L_{d}=95$  m for $\unicode[STIX]{x1D719}=8\,\%$ and $L_{d}=60$  m for $\unicode[STIX]{x1D719}=24\,\%$ .

Ultsch’s (Reference Ultsch1973) measurements also suggest that the thermal exchange is only driven during heating, and during the night-time cooling there is negligible temperature difference. So, one could reasonably estimate a cumulative flushing over several days just from the day-time exchange for cases where the total length of the vegetated region, $L_{v}$ , (which might be of the order of tens to hundreds of meters in the field or treatment wetlands) is larger than $L_{d}$ . The root layer residence time can then be estimated as $T_{R}=Lv/L_{d}$ (days). Assuming a streamwise extent of about 50 m for the floating vegetation mats (Hill, Webb & Smith Reference Hill, Webb and Smith1987), our estimate for typical thermal conditions given above indicates that the thermal exchange can fully flush the root layer in one day. Wang & Sample (Reference Wang and Sample2014) described the uptake of phosphorus and nitrogen by floating mats of soft-stem bulrush (Schoenoplectus tabernaemontani) and pickerel weed (Pontederia cordata L.) by a first-order reaction with first-order rate constants of $k=0.135\pm 0.001~\text{day}^{-1}$ (phosphorus) and $0.072\pm 0.05~\text{day}^{-1}$ (nitrogen). These rates correspond to uptake reaction time scales, $T_{k}=1/k$ , of approximately 7–14 days. Since in our example $T_{R}<T_{k}$ , the flushing can significantly enhance the overall rate of nutrient removal by maintaining higher concentration of nutrient with the root zone. Though in our example the exchange flow can fully flush the area each day that may not be the case when the average temperature difference is smaller.

7 Summary and conclusions

Highly resolved, 3-D LES were used to study the effects of a surface porous layer present on one side of a lock gate on the long-term evolution of the lock-exchange flow. This problem was motivated by thermally driven exchange between open water and a region containing floating vegetation. For all cases there is a net mass exchange within the surface current, between the porous layer and the intrusion layer beneath. Near the lock gate, fluid is advected downward from the porous layer into the intrusion layer, while the opposite happens close to the front (figure 1). Moreover, for sufficiently high values of $\unicode[STIX]{x1D719}$ , most of the flow leaves the porous layer close to origin (leading edge of the porous layer) and a region of negligible velocity within the porous layer forms until the head region is approached. For high porous layer heights ( $h_{1}/H>0.4$ ), the structure of the surface current resembles that observed in the fully vegetated case (e.g. the interface inside the porous region varies linearly with the distance from the front, starting some distance away from the front).

For sufficiently high values of the $\unicode[STIX]{x1D719}$ and $h_{1}/H$ , the surface current transitions to a quadratic-drag-dominated regime in which the front velocity is proportional to $t^{-1/4}$ , and the discharge at the lock gate position decays proportional to $t^{-1/3}$ . Both results are similar to those obtained by Ozan et al. (Reference Ozan, Constantinescu and Hogg2015) for the fully vegetated case. The decreasing discharge within the surface current also caused the velocity of the bottom current, advancing into a region with zero $\unicode[STIX]{x1D719}$ , to decay with time. Even in cases for which a region with negligible flow velocities developed within the porous layer and the flow in the two layers was far from horizontal, the front velocity decayed proportional to $t^{-1/4}$ , provided that $h_{1}/H$ was sufficiently high. For cases with a small $h_{1}/H$ , most of the surface current propagates through the intrusion layer, and the decay of the front velocity is weaker (smaller power law coefficient) compared to cases with deeper $h_{1}/H$ .

The simulations revealed the details of an unstable stratified region that develops at the interface between the porous and intrusion layers due to the difference in the velocity within the two layers. Especially for cases with a large $\unicode[STIX]{x1D719}$ , the removal of the heavier fluid from the porous layer is driven by convective exchange with lighter fluid originating in the intrusion layer (figure 3). As a result, 2-D effects are important and the flow within the porous and intrusion layer is far from being unidirectional. This means that the utility of shallow-water theory to predict the front velocity and other flow variables for the present three-layer problem is questionable.

The numerical results illustrate how the root zone is flushed by a combination of the advancing intrusion layer and the vertical exchange between the intrusion layer and the porous layer. The residence time of the root zone was estimated by tracking the volume of lighter fluid inside the porous layer, and in most cases this true residence was much shorter (in some cases by two orders of magnitude) than what would be estimated from the velocity within the root zone, highlighting the importance of the vertical exchange. Estimates of root zone residence time in field conditions suggests that the flushing of the root zone by thermally driven exchange currents can significantly enhance to the uptake of nutrients by the root system.