Lagrangian and Eulerian perspectives of turbulent transport mechanisms in a lateral cavity

Abstract

The dynamics of turbulent flows past lateral cavities is relevant for multiple environmental applications. In rivers and coastal environments, these lateral recirculating regions constitute surface storage zones, where large-scale turbulent coherent structures control the transport and fate of contaminants. Mass transport in these flows is typically represented by one-dimensional first-order equations that predict the evolution of the spatially integrated concentration between the cavity and the main channel. These models, however, cannot represent the long-term evolution of the concentration or incorporate memory effects induced by turbulence. In this investigation, we carry out large-eddy simulations (LES) of the open-channel flow with a lateral square cavity of Mignot et al. (Phys. Fluids, vol. 28, issue 4, 2016, 045104). The model is coupled with an advection–diffusion equation and a Lagrangian particle model to investigate the transport mechanisms in the cavity and across the interface. From the simulations we provide quantitative comparisons of the physical processes from both perspectives, and investigate the effects of turbulent coherent structures on residence times and trajectories from finite-time Lyapunov exponents. From the Lagrangian results, we identify general spatial distributions of time scales in the cavity associated with the dynamics of coherent structures, providing new insights into the mechanisms that drive the global transport. We also show that an upscaled model informed by LES and based on a fractional derivative captures the evolution of concentration, and the exchange between the cavity and the main channel, providing accurate predictions of mass transport and reproducing the temporal dependence observed at larger scales.

1. Introduction

Lateral recirculating regions play a fundamental role in the transport of mass and momentum in fluvial environments and in open-channel flows. Intense shear at the interface of the main turbulent channel with an adjacent volume of water creates vortex shedding and flow separation at the side, which modify the instantaneous local velocities and pressure fields (Jackson et al. 2012, 2013). In natural environments these regions are called in-stream surface storage zones (SSZs), which are typically formed by logs, boulders, vegetation or channel curvature that reshape the bank morphology. The SSZs provide habitat diversity by creating regions with slower velocities, and they also control the storage and release of contaminants along the river, affecting water quality, nutrient cycles and the fate of pollutants originated by acid-mine drainage (DeAngelis et al. 1995; Fernald, Wigington & Landers 2001; Sukhodolov et al. 2009; O'Connor, Hondzo & Harvey 2010; Sandoval et al. 2019). River restoration strategies often involve lateral recirculating regions, designed to take advantage of the hydrodynamic environment or to mitigate the impacts of pollutants in the river, modifying the channel morphology to increase residence times, and promote biogeochemical processes and exchange with subsurface flow (Mueller Price, Baker & Bledsoe 2016; Juez et al. 2018b).

Contaminant transport in these lateral recirculating regions is typically analysed from a global perspective, representing the concentration of passive scalars by one-dimensional (1-D) models (O'Connor et al. 2010; Jackson et al. 2012; Knapp & Kelleher 2020). Since these approaches are based on the gradient diffusion hypothesis and an exchange coefficient or by defining multiple exchange regions with more parameters, they have a limited accuracy in flows dominated by large-scale and highly 3-D turbulent coherent structures (Khosronejad et al. 2016).

To reduce the effects of specific arbitrary geometries and external factors, recent investigations have focused on rectangular geometries, i.e. lateral cavities for open-channel flows (e.g. Mignot et al. 2016, 2017; Sanjou, Okamoto & Nezu 2018; Ouro, Juez & Franca 2020). The flow past cavities has been extensively studied for multiple hydrodynamic and aerodynamic applications (Rockwell & Naudascher 1979; Liu & Katz 2013; Tuna & Rockwell 2014; Karimpour, Wang & Chu 2021), including processes such as noise generation and flow-induced vibrations. The experimental analysis of mass transport, however, has remained limited to dye experiments from an Eulerian perspective. Measurements of the concentration of passive scalars on planes inside the cavity or global measurements across the water depth have characterized the mass exchange through optical techniques, providing quantitative estimations of the transport across the interface between the cavity and the main channel (Uijttewaal, Lehmann & van Mazijk 2001; Weitbrecht, Socolofsky & Jirka 2008; Sanjou & Nezu 2013; Mignot et al. 2017).

The recent work of Engelen et al. (2021) has been the first experimental investigation to analyse the exchange process from a Lagrangian perspective. They used 3-D particle tracking velocimetry to capture the trajectories of microspheres slightly denser than water, and classified their paths at the interface to compute the global exchange coefficient. The observations revealed an intricate Lagrangian dynamics, with particles zigzagging and repeatedly crossing the interface, being either entrained or ejected into the main channel. Quantifying the transport of particles for these measurements required a statistical approach to interpret the results from Lagrangian motion with caution, as some of the particles were lost during the tracking observations. Nevertheless, these results offered valuable insights into the impact of small-scale processes on mass exchange over time.

A comprehensive understanding of exchange processes driven by turbulent coherent structures is therefore critical for upscaling the results to the larger-scale 1-D models, and accurately quantifying transport for environmental or restoration applications. Previous studies have found that the initial mass exchange seems to follow the conventional first-order 1-D equation to represent the initial evolution of the concentration (Uijttewaal et al. 2001; Sandoval et al. 2019). However, as time progresses and smaller concentrations are exchanged, the global model can significantly overpredict the mass transport. This is especially problematic for contaminants that are toxic at low concentrations, as they may become trapped for longer times in the cavity, changing the exposure periods and the spatial distribution in the environment.

High-resolution numerical simulations can characterize the transport mechanisms by capturing the dynamics of the turbulent coherent structures, since the collective effects of turbulence on mass exchange in the cavity and at the interface with the main channel are perceived as memory effects in upscaled 1-D representations. In this investigation we perform large-eddy simulations (LES) of the lateral square cavity experiments of Mignot et al. (2016) to study the leading mechanisms of mass transport. The main objective of this research is to understand the transport dynamics from Eulerian and Lagrangian perspectives, providing physical insights into the effects of small-scale processes on global transport. The present work focuses on quantifying residence times and characterizing regions of the flow that give rise to emergent large-scale behaviour with longer temporal dependence by studying the concentration of a passive scalar or the trajectories of particles, developing models that consider the influence of subgrid scales. Through this analysis, different time scales observed in mass transport are identified and linked to the complex spatial distribution of residence time statistics and the divergence of trajectories in the cavity, quantified with finite-time Lyapunov exponents.

The paper is structured as follows: in § 2 the LES numerical model and the Eulerian and Lagrangian mass transport equations to study the exchange between the cavity and the main channel are described. Subsequently, in § 3 the results of the time-averaged flow field and instantaneous features are presented and compared with the experimental measurements. In § 4 the main characteristics of mass transport from both Lagrangian and Eulerian perspectives are discussed, analysing the small-scale effects on global transport represented from emergent large-scale relations. The statistics of particle exchange between the cavity and the main channel are explored in § 5, by studying the distribution of residence times inside the cavity as a function of the initial conditions, and expanding the analysis to quantify the separation of trajectories from the finite-time Lyapunov exponent in § 6. The conclusions in § 7 contain a summary of our findings and ideas for future research.

2. Numerical framework for LES and mass transport

2.1. Governing equations of the flow

The governing equations for LES are the 3-D spatially filtered unsteady, incompressible Navier–Stokes equations. They are expressed in non-dimensional form and Cartesian coordinates as follows:

(2.1)\begin{gather} \frac{\partial \bar{u}_i}{\partial x_i} = 0, \end{gather}

(2.2)\begin{gather} \frac{\partial \bar{u}_i}{\partial t}+u_j\frac{\partial \bar{u}_i}{\partial x_j}={-}\frac{\partial \bar{P}}{\partial x_i} + \frac{1}{Re} \frac{\partial^2 \bar{u}_i}{\partial x_j\partial x_j} - \frac{\partial \tau_{ij}}{\partial x_j}, \end{gather}

where $\bar {u}_i$ are the filtered components of velocity with the bar denoting the grid filter, $\bar {P}$ is the filtered pressure, $Re$ is the Reynolds number and $\tau _{ij}$ are the components of the subgrid-scale (SGS) stress tensor. In the above equations all quantities have been made non-dimensional using the flow depth $h$ and the bulk velocity in the main channel $U_b$ as characteristic length and velocity scales, respectively. Consequently, the Reynolds number is defined as $Re = U_bh/\nu$, where $\nu$ is the kinematic viscosity of the fluid.

The SGS tensor is defined as $\tau _{ij} = \overline {u_i u_j} - \bar {u}_i \bar {u}_j = -2 \nu _{t} \bar {S}_{ij}$, with $\bar {S}_{ij} = ({\partial \bar {u}_i}/{\partial x_j} + {\partial \bar {u}_j}/{\partial x_i})/2$ being the resolved-scale strain-rate tensor. The eddy viscosity is calculated as $\nu _t = C_v \bar {\varDelta } \sqrt {k^{sgs}}$ based on a one-equation LES turbulence model, the locally dynamic kinetic energy model developed by Kim & Menon (1999), where $k^{sgs} = (\overline {u_k u_k} - \bar {u}_k \bar {u}_k)/2$. The model considers the filter size $\bar {\varDelta }$ from the grid, solving a transport equation for the SGS turbulent kinetic energy $k^{sgs}$

(2.3)\begin{equation} \frac{\partial k^{sgs}}{\partial t}+\frac{\partial }{\partial x_j}(\bar{u}_j k^{sgs})={-}\tau^{sgs}_{ij}\bar{S}_{ij}+\frac{\partial}{\partial x_j}\left[\left(\frac{1}{Re}+\nu_t\right)\frac{\partial k^{sgs}}{\partial x_j}\right]-C_\epsilon\frac{(k^{sgs})^{3/2}}{\bar{\varDelta}}. \end{equation}

In the above equation, the three terms on the right-hand side correspond to the production, diffusion and dissipation of $k^{sgs}$. The dynamic coefficients $C_\epsilon$ and $C_v$ are employed to determine the SGS dissipation rate $\epsilon _{ij}$ and the turbulent viscosity $\nu _t$, respectively, and they are computed from the resolved flow field considering a test filter size $\hat {\varDelta } = 2 \bar {\varDelta }$. The coefficient $C_v$ is obtained by applying the least-square method suggested by Lilly (1992) and $C_\epsilon$ based on the SGS dissipation $\epsilon ^{sgs}$

(2.4)\begin{gather} C_v = \frac{1}{2} \frac{L_{ij}\sigma_{ij}}{\sigma_{lm}\sigma_{lm}}, \end{gather}

(2.5)\begin{gather} C_\epsilon = \left(\frac{1}{Re} + \nu_t\right) \frac{\hat{\varDelta}}{{(k^{test})}^{3/2}} \left[\widehat{\frac{\partial \bar{u}_i}{\partial x_j}\frac{\partial \bar{u}_i}{\partial x_j}} - \frac{\partial \hat{\bar{u}}_i}{\partial x_j}\frac{\partial \hat{\bar{u}}_i}{\partial x_j}\right], \end{gather}

where $L_{ij} = \widehat {\bar {u}_i \bar {u}_j} - \widehat {\bar {u}_i}\widehat {\bar {u}_j}$ is the Leonard stress tensor, $\sigma _{ij} = \hat {\varDelta } \ \hat {\bar {S}}_{ij} \sqrt {k^{test}}$ and $k^{test} = (\widehat {\bar {u}_k \bar {u}_k} - \hat {\bar {u}}_k\hat {\bar {u}}_k)/2$.

2.2. Eulerian mass transport model

The solute transport is based on an advection–diffusion equation for a passive scalar, to predict the transport of the volumetric concentration

(2.6)\begin{equation} \frac{\partial \bar{C}}{\partial t}+\bar{u}_j\frac{\partial \bar{C}}{\partial x_j}=\frac{\partial}{\partial x_j} \left[\left(D+D_t\right) \frac{\partial\bar{C}}{\partial x_j}\right], \end{equation}

where $\bar {C}$ is the filtered concentration, $\bar {u}_j$ are the resolved velocity components of the flow and $D$ and $D_t$ are the non-dimensionalized molecular and turbulent diffusion coefficients, respectively. The diffusion coefficients are defined from the molecular Schmidt number $Sc=Re^{-1}/D$ equal to 100, and a turbulent Schmidt number $Sc_t=\nu _t/D_t$ equal to 1, which represent correctly the transport of passive scalars in water at ambient temperature (Gualtieri et al. 2017).

2.3. Lagrangian model for fluid particles

In the Lagrangian approach, the path of discrete point elements of fluid without inertia are tracked in time. These simulation of particle trajectories can represent small neutrally buoyant particles associated with nutrients and minerals, pollutants such as microplastics or suspended fine sediment particles with low inertia and small Stokes number (Escauriaza & Sotiropoulos 2009; Abarca et al. 2017).

To compute the fluid velocity at the particle location, the fourth-order partial Hermite interpolation is used, which has been successfully employed in previous investigations of particle simulations in turbulent flows (e.g. Balachandar & Maxey 1989; Choi, Yeo & Lee 2004; Polanco 2019). To incorporate the effects of the SGS velocities, a stochastic equation is used to obtain the local flow velocity considering the effects of unresolved scales (Haworth & Pope 1986; Pope 1994; Gicquel et al. 2002; Fede et al. 2006; Berrouk et al. 2007; Marchioli 2017). In this case, the SGS particle model of Vinkovic et al. (2006) is implemented to obtain a velocity fluctuation ($v_i'$) that is added to the filtered flow field at the particle position ($\overline {v_i}$). The motion of each particle is therefore computed from the following set of equations:

(2.7)\begin{gather} \frac{{\rm d}\kern0.7pt x_i}{{\rm d} t} = \bar{v}_i + v_i', \end{gather}

(2.8)\begin{gather} {\rm d} v_i' = \left( -\frac{1}{T_L} + \frac{1}{2 k^{sgs}} \frac{{\rm d} k^{sgs}}{{\rm d} t} \right) v_i'\, {\rm d}t + \sqrt{\frac{4 k^{sgs}}{3 T_L}} \,{\rm d} \eta_{i}(t). \end{gather}

Equations (2.7) and (2.8) correspond to the particle trajectory and the stochastic SGS model based on a Langevin equation (Vinkovic et al. 2006), respectively. The first term in (2.8) is the deterministic memory of the velocity fluctuation, and the second is a zero-mean Wiener noise scaled by the SGS turbulent kinetic energy, where the time scale is equal to $T_L = {4 k^{sgs}}/{3 C_0 \epsilon^{sgs}}$, $\epsilon^{sgs} = {C_{1} k^{{sgs}^{3/2}}}/{\tilde {\triangle }}$ and ${\rm d} \eta _{i} = \eta \sqrt {t}$ ($\eta \sim N(0,1)$). The values of the model constants $C_0$ and $C_{1}$ are typically assumed as $C_0=2.1$ and $C_1=1$ (Haworth & Pope 1986; Pope 1994; Gicquel et al. 2002; Fede et al. 2006; Berrouk et al. 2007).

2.4. Computational details

The 3-D cavity flow configuration studied by Mignot et al. (2016) was simulated by considering the entire experimental set-up and the same Reynolds number equal to $Re=11\,667$. The rectangular channel length is 4.9 m long and $B=0.3$ m wide.

The computational domain considers the entire experimental channel, spanning $70h \times 4.29h \times h$ in the streamwise ($x$), spanwise ($y$) and vertical ($z$) directions, respectively, as shown in figure 1. The lateral square cavity is located at the centre of the total length with sides of dimensions $4.29h$.

Figure 1. Geometry of the computational domain (not to scale). The grid comprises the entire length of the experimental set-up.

The structured numerical grid, shown in figure 2, consists of $981 \times 281 \times 101$ points in the $x$, $y$ and $z$ directions, corresponding to a total of approximately 27 million nodes with a resolution at solid boundaries of $\triangle x_{i}^+ \le 1$. No-slip boundary conditions are applied at solid surfaces, and the free surface is represented as symmetry boundary due to the small Froude number of the flow ($Fr = 0.2$), which represents values observed in low-gradient streams, where no significant deformation of the free surface is observed. At the inflow of the computational domain that coincides with the channel entrance, we prescribe a uniform bulk flow to which we add a synthetic unsteady stochastic inlet using the random flow generator of Smirnov, Shi & Celik (2001), representing the turbulence of the experimental inlet in Mignot et al. (2016). The entire channel is simulated to make sure that the development of the boundary layer is reproduced, and that the outlet has no influence on the dynamics observed in the cavity section. The minimum, maximum and mean grid resolutions in wall units are presented in table 1. The maximum spacings in $x$, $y$ and $z$ correspond to the regions at the outlet, main channel and mid-channel, respectively. In table 1, the shear-layer region is centred at the interface, considering $20\,\%$ of the cavity volume, and a streamwise extension slightly longer than the cavity.

Figure 2. Details of the computational grid in the cavity region. The $X$ axis corresponds to the flow direction in the main channel.

Table 1. Maximum, minimum and mean spatial resolutions at different regions in wall units.

The filtered Navier–Stokes equations in generalized curvilinear coordinates are solved using a dual time-stepping artificial compressibility iteration scheme, employing second-order-accurate finite-volume method on a non-staggered computational grid. The discrete equations are advanced in time by adopting the pressure-based implicit preconditioner of Sotiropoulos & Constantinescu (1997), enhanced with local time stepping and V-cycle multigrid acceleration. The one-equation turbulence model is integrated with second-order schemes, advancing the solution in pseudo-time using the standard alternate direction implicit scheme. Applications and the performance of this model have been tested and discussed in great detail in a series of previous papers (Paik, Escauriaza & Sotiropoulos 2007, 2010; Escauriaza & Sotiropoulos 2011a,b,c; Link et al. 2012; Gajardo, Escauriaza & Ingram 2019; Gotelli et al. 2019; Sandoval et al. 2019; Soto-Rivas, Richter & Escauriaza 2019; Sandoval et al. 2021), in which the accuracy of the methods has been demonstrated by qualitative and quantitative comparisons with available experimental data, in terms of mean flow quantities and turbulence statistics.

The Eulerian mass transport equation (2.6) is solved using the algorithm implemented in Sandoval et al. (2019), advancing in time by using the dual-time-stepping method. The advective term of the transport equation is discretized using the second-order-accurate upwind scheme with a flux limiter (Leonard 1991) to resolve steep gradients and avoid unphysical oscillations, while central differencing is employed for the diffusive flux.

The particle trajectories in the Lagrangian model are obtained from (2.7) integrated in time, using a third-order Runge–Kutta scheme by coupling the LES model with the open-source code LIGGGHTS (Kloss et al. 2012) for particle simulations, as in the recent work of González et al. (2017) and Escauriaza et al. (2023). The stochastic velocity equation (2.8) is solved by an Euler–Maruyama scheme (Gicquel et al. 2002; Fede et al. 2006). All the simulations for the flow and mass transport are carried out using a non-dimensional physical time step of $\triangle t=0.01$ (0.0042 s), which is scaled with the water depth $h$, and the bulk velocity $U_b$.

The difference between Eulerian and Lagrangian approaches is due to the effect of diffusive transport in the continuum simulations, which has been also observed in simulations of mass transport in the blood flow of aneurysms (Reza & Arzani 2019) and in porous media (Henri & Diamantopoulos 2022). The numerical scheme of Leonard (1991) implemented in the solution of (2.6) has demonstrated its accuracy in capturing advective transport in open-channel flows for passive and active contaminants (Lin & Falconer 1997; Gross, Koseff & Monismith 1999; Devkota & Imberger 2009; Paik, Eghbalzadeh & Sotiropoulos 2009), and in LES of scalar transport (Sharan, Matheou & Dimotakis 2018). To verify the Eulerian results presented in this work, we carried out the procedure outlined by Muppidi & Mahesh (2008), performing a separate calculation with a significantly smaller time step. Muppidi & Mahesh (2008) commented on the potential influence of numerical dissipation on mass transport, and integrated the advection–diffusion equation in direct numerical simulations using a smaller time step to ensure the correct representation of advective mass fluxes at high resolutions. The results we obtained here with a higher temporal resolution, reducing the time step by one order of magnitude ($\triangle t=0.001$ or 0.00042 s), showed no changes on the instantaneous mass transport and statistics that we report on this investigation.

3. Flow in the lateral cavity

In this section we present briefly the time-averaged and instantaneous flow field obtained from the LES. We compare the computed results with the experimental measurements of Mignot et al. (2016), to show the good agreement of the simulations with the time-resolved 2-D particle-image velocimetry (PIV) observations reported in a horizontal plane at $z/h=0.7$. In all the results the velocity and time are made non-dimensional using the bulk velocity and the water depth.

3.1. Mean flow

The time-averaged flow field in the cavity is mostly two-dimensional, with a simple structure that is illustrated in the different planes depicted in figure 3. The flow exhibits the following features: (i) a shear layer with a strong transverse gradient of streamwise velocity and converging streamlines at the interface; (ii) a large vortical structure induced by the flow in the main channel, which occupies most of the cavity volume; (iii) a time-averaged inflow that enters the cavity from the trailing corner at the downstream boundary and close to the bed and a time-averaged outflow that occurs close to the upstream boundary mainly at the middle of the channel; (iv) secondary small corner vortices generated by the interaction of the rotating flow with the walls. In figure 3 we can also observe the disparity of velocity magnitudes inside the cavity, with velocities at the centre of the vortex that are one order of magnitude smaller compared with the shear layer, inducing a significant variability of time scales associated with these structures in the lateral region, as observed in previous investigations (Sanjou & Nezu 2013; Drost et al. 2014).

Figure 3. Non-dimensional time-averaged velocity magnitude and streamlines. (a) Experimental horizontal plane at $z/h=0.7$ for comparison. (b) Transverse plane at the interface showing the mean transverse velocity and where the time-averaged inflow and outflow occurs in the cavity; (c) and (d) correspond to the horizontal planes where most of the inflow and outflow occur, respectively.

The measured and calculated streamwise velocity profiles and gradients at the horizontal plane $z/h=0.7$ depicted in figure 4, show that the simulations are in very good agreement with the measurements, reproducing the velocity magnitude across the interface and the spreading of the time-averaged shear layer in the presence of the adverse pressure gradient imposed by the impinging point at the trailing edge of the cavity, which increases the thickness of the shear layer downstream.

Figure 4. (a) Time-averaged streamwise non-dimensional velocity profiles at the interface between the cavity and the main channel. (b) Non-dimensional streamwise velocity gradient. The solid lines correspond to the present simulations, and the red circles to the experimental data of Mignot et al. (2016).

3.2. Instantaneous flow field

To identify the coherent dynamics of the flow in the cavity and at the interface with the main channel, we visualize the instantaneous resolved vertical vorticity distributions in figure 5, and the 3-D instantaneous $q$-criterion iso-surfaces in figure 6.

Figure 5. Sequence of non-dimensional resolved vertical vorticity contours in the horizontal plane $z/h=0.7$. The dashed line aids in visualizing the precise location of the vortex cores with respect to the cavity zone. The images are separated by 0.86 s , equivalent to 2 non-dimensional times ($2\times h/U_b$). This temporal interval corresponds to half of the period of the lead frequency of vortex shedding.

Figure 6. Three-dimensional instantaneous evolution of the shear layer visualized with $q$-isosurfaces ($q=1.0$). The 3-D images are coloured by non-dimensional pressure, separated by a temporal interval equal to half the period of the vortex shedding. The locations of prominent vertical-axis and longitudinal vortices are indicated in the panels. Note that, to identify the instantaneous low-pressure cores, each panel has a different scale.

The images are snapshots of vorticity near the interface separated by 200 time steps (0.84 s), and the first and last encompass an entire period of the shear-layer vortex shedding, as described below. The sequence for the horizontal plane at $z/h=0.7$ shows the general characteristics of the instantaneous flow driven by a shear-layer instability that develops at the interface, which encounters an adverse pressure gradient produced by the flow impingement on the downstream corner of the cavity interface. Vortex shedding is initiated at the leading edge of the interface, generating a sequence of vertically aligned vortices that travel and impinge on the trailing corner. The instantaneous shear layer consists of up to two large counter-clockwise vortical structures that are continuously shed in a seemingly periodic manner, with associated smaller clockwise structures. These vortices divide when they impact the downstream corner, producing a wall jet inside the cavity that induces the rotation in the enclosed volume, in agreement with what is observed in the experiments of Mignot et al. (2016). The interface plane is marked with a dashed line in the figure to show that the core of shear-layer vortices are sometimes encountered inside the cavity, while outside excursions of the vortices can also occupy a section of the main channel in a low-frequency flapping motion of these structures across the interface.

The visualization of the flow in figure 6 also shows that longitudinal coherent structures engage in complex interactions with the large-scale vortices inside and outside the cavity, with a highly 3-D structure of the flow in the interface and cavity regions illustrated by the $q$ isosurfaces. The $q$-criterion (Hunt, Wray & Moin 1988) is defined in regions of the flow where the difference between the norm of the resolved rate of rotation ($\varOmega _{ij}$) and the resolved rate of strain tensors $S_{ij}$ is positive. This scalar quantity, obtained as $q=0.5(\varOmega _{ij} \varOmega _{ij} - S_{ij}S_{ij})>0$, can help identify vortical structures in regions where local rotation rate is larger than the strain rate, even though it has limitations where there is high shear and both quantities are very large.

The images and animations in the supplementary material show a very rich dynamics, in which columnar vortices of the shear layer that emerge at the interface are elongated and connected by longitudinal worm-like vortical streaks, impinging at the trailing edge of the cavity. The figure sequence helps to clarify the mechanisms that generate the rotation inside the lateral cavity, and how the wall jet produced at the interior downstream wall is heavily influenced by the periodic vortex shedding and internal feedback with the shear layer, which experiences a continuous flapping modulated by the low-frequency rotation inside the cavity (see similar processes described by Liu & Katz 2013).

These comparisons demonstrate the ability of the LES model to resolve the large-scale flow patterns, and leading mechanisms of momentum transport captured in previous experimental analyses (Liu & Katz 2013; Mignot et al. 2016). In figure 7, we also show the power spectra for time series of the three velocity components at two representative points in the domain: at the interface of the cavity with the main channel, and close to the end boundary near the opposite wall. The analysis shows that the model reproduces the leading frequencies measured in the experiment of Mignot et al. (2016), with a vortex-shedding frequency of the shear layer equal to 0.58 Hz, which corresponds to a period of four dimensionless times (the Strouhal number, $St_r = fh/U_b$, is equal to 0.25). The two frequency peaks in figure 7(a) show two harmonics ( $f_1/f_2 = 0.5$), confirming the strong influence of the shear layer in the slower and the less energetic flow dynamics inside the cavity, dominated by a single large vortical structure. It is worth noting that the spectrum inside the cavity, as depicted in figure 7(b), corroborates the lower energy levels in the slow recirculation flow near the opposite wall, in contrast to the turbulent shear layer at the interface.

Figure 7. Three-component velocity spectra at two different positions of the lateral cavity flow. (a) At the interface; and (b) inside the cavity, near the end boundary. Curves in blue, red and black correspond to spectra of streamwise, transverse and vertical velocity components, respectively.

The $-3$ slope observed at low frequencies could be due to the fact that the largest turbulent scales at the interface and inside the cavity are mostly two-dimensional, with vertically oriented axes of structures confined by the depth, creating a range of turbulent scales for which the 3-D energy transfer is suppressed (Jirka 2001). Analyses of similar flows and detailed discussions of this phenomenon are discussed in Uijttewaal & Booij (2000) and Proust & Nikora (2020).

4. Mass transport

Since the LES shows good qualitative and quantitative agreement with the experimental measurements, we use the numerical simulations to explore how the turbulent structures at the interface and inside the lateral cavity influence the mass exchange with the main channel. The main objective is to provide insights on how the exchange processes driven by turbulent coherent structures have a global effect on the transport of contaminants in open-channel flows and other environmental applications.

Here, we perform an analysis of mass transport from Eulerian and Lagrangian perspectives, starting from a converged instantaneous solution of the flow. For the continuum approach, we fill the cavity with a passive scalar assuming volumetric concentration $C=1.0$ in the entire lateral volume and solve (2.6), as described in the previous section. For the Lagrangian approach, we uniformly distribute $1.08\times 10^5$ passive particles inside the cavity and integrate their trajectories from the instantaneous flow field, using (2.7) and (2.8) to compute transport statistics. The findings presented here have also been tested with particle resolutions smaller by an order of magnitude, yielding consistent and identical results.

We analyse the time series of mass flux across the interface and study the implications of the small-scale processes on the global mass exchange and transport, which can be quantified from large-scale representations typically used in the study of contaminant transport in fluvial systems (O'Connor et al. 2010; Sandoval et al. 2019).

The evolution of transport for Eulerian and Lagrangian cases is depicted in figure 8, showing how the shear-layer vortices emerging at the interface entrap the contaminant in the lateral cavity, transporting mass downstream in the main channel. The time sequence of concentration contours in the Eulerian case, and particles coloured by their velocity magnitude in the Lagrangian simulations, can help visualize the leading role of these large-scale vortices, which also entrain fluid from the main channel into the recirculating region as the cavity is emptied.

Figure 8. Plan view of transport evolution from the lateral cavity in the Lagrangian (a–c) and Eulerian (d–f) simulations. The sequences show the influence of the shear layer on the mass exchange toward the main channel. Both sequences correspond to one shear-layer period (1.72 s, or 4 non-dimensional times), but they are displayed at different physical times to enhance clarity of the transport processes in each case. Additional details can be observed in the supplementary movies of the paper.

To provide quantitative insights into the global decrease of the concentration in the cavity we study the flux across the interface, which is the large-scale response of the system to the overall effects of the turbulent coherent structures on mass transport, between the channel and the lateral storage zone. For the Lagrangian data we computed the flux as the difference in the number of particles inside the cavity at two consecutive times divided by the total number of particles and the time window. For the Eulerian data we considered the difference in the volume-averaged concentration divided by the time. The time series of mass flux in figure 9, show that the LES captures the periodicity and transport fluctuations across the interface driven by the shear layer. The filtered data at larger time scales show that the leading frequency of transport observed in these plots coincides with the periodic velocities induced by the the vortex shedding at the interface. The fluxes at the interface also show negative excursions, or negative transport, due to the re-entrance of mass in the cavity, which is a phenomenon clearly captured by the Lagrangian model with particle trajectories that perform a complex zigzagging at the interface, as will be subsequently shown (see also Engelen et al. 2021).

Figure 9. Time series of Lagrangian and Eulerian mass flux across the cavity interface. The mass exchange reveals a periodicity linked to the vortex-shedding frequency. The figure at the right shows the more fluctuating Lagrangian flux with larger excursions.

To evaluate the mass transport at a global scale, lateral recirculating regions are typically analysed by using transient storage models. A 1-D transport equation for the averaged concentration represents the exchange with the main channel as proportional to the concentration gradient. Through the instantaneous information provided by LES, we can understand in detail the driving mechanisms for the mass exchange and improve the large-scale formulations computing the instantaneous volume-averaged concentration inside the cavity. The global large-scale model represents the time evolution of the concentration in the cavity as follows:

(4.1)\begin{equation} \frac{{\rm d}\langle C\rangle}{{\rm d} t} = \frac{Q_{mc}}{\forall_c} (C_m - \langle C\rangle), \end{equation}

where $\langle C\rangle$ and $C_m$ are the spatially averaged bulk concentration in the cavity and in the main channel respectively, $Q_{mc}$ is the volumetric flux from the main channel to the cavity and $\forall _c$ is the cavity volume. The instantaneous flux across the interface is computed from an exchange velocity across the interface $E$, and is expressed in terms of the interface width $W$

(4.2)\begin{equation} \frac{{\rm d} \langle C\rangle}{{\rm d} t} = \frac{E}{W} (C_m - \langle C\rangle). \end{equation}

The exchange velocity is estimated as proportional to the bulk velocity in the main channel $E=kU_b$, where $k$ is known as the dimensionless mass exchange coefficient. The model for the cavity concentration is written as

(4.3)\begin{equation} \frac{{\rm d} \langle C\rangle}{{\rm d} t} = \frac{k U_b}{W} (C_m - \langle C\rangle). \end{equation}

For the cavity that is initially full, the analytical solution of this equation shows that the concentration decays exponentially, such that

(4.4)\begin{equation} \langle C\rangle/C_0 = e^{{-}t/\tau}, \end{equation}

where $C_0$ is the initial concentration inside the cavity and corresponds to the maximum value, and $\tau$ is the mean residence time, related to the mass exchange coefficient $k$ as follows:

(4.5)\begin{equation} \tau = \frac{W}{k U_b}. \end{equation}

The magnitude of the exchange is typically estimated from experiments, using dye releases or measuring the transverse velocity at the interface, fitting the solution to the data (Knapp & Kelleher 2020). In our numerical simulations we can directly compute the contaminant mass remaining in the cavity, either from the instantaneous concentration field or from the particles in the Lagrangian model.

Figure 10 shows (4.4) with the flux computed from LES in a semilog plot, where the exchange coefficients and time scales associated with this model are related to the slope of the curves. As expected, the Eulerian and Lagrangian results yield different exchange coefficients, due to the mixing effects of the continuum approach. The contaminant is diluted in fluid from the main channel that has been entrained in the cavity by the shear layer. Clear fluid from upstream enters the cavity, facilitating the mixing with the contaminant by advective and diffusive processes. In the Lagrangian model, however, there is no molecular or turbulent mixing. This effect was verified in our calculation by changing the number of particles and by performing an additional Eulerian calculation at a higher temporal resolution, to discard the potential influence of the numerical discretization of the advective terms.

Figure 10. Comparison for the evolution of averaged concentration from LES, and the large-scale model with a fractional derivative of order $\alpha =0.96$. For Eulerian and Lagrangian cases the global model captures the evolution of the instantaneous mass inside the cavity, representing the memory effects induced by the collective interactions of turbulent coherent structures.

The numerical results show that the exponential decay of this first-order model separates from the actual decay when there is still $40\,\%$ of the initial total mass in the cavity. The exchange coefficient $k$ is estimated from the instantaneous average concentration in the Eulerian case, whereas, for the Lagrangian data, the value of $k$ is obtained by counting the fraction of the initial number of particles that remains inside the cavity over time.

The non-dimensional global residence times obtained from Lagrangian and Eulerian data are $\tau _L=102$ and $\tau _E=130$, respectively. This statistic represents the average time that particles or the passive scalar remain inside the cavity, if we assume that transport is proportional to the concentration difference between the two zones. The mean residence times obtained here are similar to the values obtained by previous studies for cavity flows at different aspect ratios $L/W$ (Jackson et al. 2012; Sanjou & Nezu 2013; Drost et al. 2014). The mass exchange coefficients from Lagrangian and Eulerian data are $k_L=0.042$ and $k_E=0.033$, respectively, and are also consistent with data in the literature for lateral cavity flows (Weitbrecht & Jirka 2001; Sanjou & Nezu 2013; Mignot et al. 2016, 2017; Sandoval et al. 2019).

From the work of Sandoval et al. (2019) on a natural lateral cavity, the global model based on a time-fractional differential equation can be used to better describe the evolution of mean concentration in the cavity, which captures the complex emergent dynamics that is the consequence of the interactions of turbulent coherent structures inside the recirculating region. We maintain the same exchange coefficients $k_L=0.042$ and $k_E=0.033$, solving the following equation:

(4.6)\begin{equation} \frac{{\rm d}^\alpha}{{\rm d}t^\alpha}\langle C\rangle={-}\frac{kU_b}{W}\frac{h_m}{h}\left(\langle C\rangle-C_m\right), \end{equation}

where $\alpha$ corresponds to the order of the fractional derivative. The solution of this equation corresponds to a Mittag–Leffler function, which generalizes the exponential solution.

In figure 10 we show that, at short time scales, the first-order linear model represents correctly the initial emptying process that transitions to a heavy-tail behaviour with an apparent change of the decay rate at approximately $t=120$ for $(C=0.4)$. The order of the fractional derivative $\alpha =0.96$ reflects the memory effects induced by the collective dynamics of the turbulent flow. Stronger memory effects have been also observed in natural surface storage zones, enhanced by complex bathymetries that induce upwelling events and large 3-D flow features as discussed in Sandoval et al. (2019).

All these characteristics observed in the global dynamics of mass transport, such as the periodicity of the fluxes and the evolution of the cavity concentrations that departs from classical gradient diffusion transport, point to the conclusion that multiple time scales influence the exchange with the main channel.

Therefore, memory effects generate a wider distribution of residence times in the cavity. We investigate their distribution by constructing their histogram, as shown in figure 11. The power-law probability density function associated with the memory effects revealed by the fractional model, with the tails decaying with a $-1.96$ slope in a semilog plot, is linked to the exponent of the fractional derivative and represents better the extreme values of the residence time distributions.

Figure 11. Histogram of Lagrangian residence times in the entire cavity volume, showing the theoretical decay given by the exponential function and the equivalent power-law distribution from the fractional model, which can represent better the extreme values of the distribution.

From the LES results and particle trajectories, we can further analyse the effects of coherent structures on global statistics and identify the topology of the flow in the cavity that explains the global observations in experiments and simulations. In the following sections we study in detail the spatial distribution of the residence times inside the lateral recirculating region, and compute statistics that are closely linked to particle trajectories across the interface.

5. Spatial residence time statistics

The analysis of the global transport reveals a wide distribution of residence times imposed by the turbulent coherent structures of the flow. Three-dimensional maps of the residence time distribution inside the cavity are generated by the following procedure: (a) particles are randomly distributed occupying the entire cavity volume; (b) the Lagrangian simulation is performed, tracking each individual particle position and recording the time they remain inside the cavity (residence time); and (c) after the simulation is finished, the particle initial position is coloured by the residence time obtained in the previous step. Therefore, by analysing how long particles remain inside the cavity and relating it to their initial positions, we can visualize the spatial distribution of these time scales and identify sections with similar residence times occupying distinct volumes inside the cavity.

In figure 12(a) we plot three different planes of the 3-D map of residence times, coloured by their magnitude at the initial location of each particle. The horizontal plane $x\unicode{x2013}y$, and vertical $y\unicode{x2013}z$ and $x\unicode{x2013}z$ planes are organized from top to bottom, showing the large disparity of time scales in the lateral recirculating region, which seem to be organized according to the large-scale vortices observed in the simulations. Two relevant aspects of the Lagrangian mass transport can be visualized: (i) there are global patterns of residence times that appear in the cavity volume, with smaller values near the interface and upstream wall, and a region occupied by longer time scales in the vicinity of the central vortex core; and (ii) inside these global sections that gather trajectories with similar residence times, there is an irregular distribution of scales and local fluctuations induced by turbulence, which are embedded in most of the cavity volume. The blue regions emerge as the clear imprint of shear-layer vortices, with short residence times at the locations close to the interface and return leg of the central core rotation, where these vortices can easily extract a large number of particles from the cavity.

Figure 12. Spatial distribution of particle residence times according to their initial positions. Slices correspond to $x\unicode{x2013}y$, $y\unicode{x2013}z$ and $x\unicode{x2013}z$ planes from top to bottom. (a) Non-averaged values of initial position of particles coloured by their residence time, showing regions with groups of particles with similar time scales. (b) Values averaged over the first four dimensionless times (equivalent to one shear-layer period).

To reduce the grainy structure produced by local fluctuations, we perform a period-average filtering, based on the scale of the shear layer $T=4.0$, where these specific global regions can be better identified. The plots in figure 12(b) show the same planes but computed with the period-averaged residence time, capturing a clearer spatial distribution of the large-scale global patterns. Regions of different time scales separated by two orders of magnitude are observed, the largest is associated with the main vortex core in the central region, with residence times above 120; and the second at the interface and upstream boundary influenced by the shear layer, with residence times lower than 40. All other regions have residence times between these two values (40–120). This finding is similar to the multiple-region models proposed by Drost et al. (2014) and Sandoval et al. (2019), which consist of two exchanges zones corresponding to the perimeter of the main eddy and the core of the recirculation region. Our findings, however, differ slightly from those of Drost et al. (2014) and correspond to two regions with completely different time scales: the vortex core time scale ($t_{vc}=150$) and the interface and upstream boundary influenced by the shear-layer time scale ($t_{sl}=4$). We found that all other areas inside the cavity are influenced by more than one time scale and mainly by the eddy turnover time scale ($t_{ed}=50$). Our findings provide a more complete understanding on the identification of regions that trap and release mass in lateral recirculating regions and are in agreement with Juez et al. (2018a), who observed that sediment particles tended to flow into the core of the main recirculating eddy and settle there (as also discussed in Ouro et al. 2020).

The fluctuations and spatial variability observed inside the cavity, as depicted in the fine-grained figure 12(a), indicate that particles starting very close to each other can have significantly different residence times, highlighting the complex Lagrangian dynamics of the flow driven by the interplay between the central slow vortex core motion, and particles approaching the interface that can rapidly leave the cavity or remain inside for arbitrarily long times. A closer analysis of the particle trajectories shows that streamwise vortices between pairs of the columnar vortices that are part of the shear layer induce a vertical motion that reintroduces particles inside the lateral cavity. The sequence of images in figure 13 shows the instantaneous positions coloured by residence time and vertical velocity, in which particles at the interface that approach the region between two vertical axis vortices are re-entrained and remain inside for a longer period. This can be identified by particle pockets of green and yellow residence times at the shear layer. Animations of the flow show that these particles oscillate at the interface, which was described as particle zigzagging by Engelen et al. (2021), leaving and re-entering the cavity one or more times. Re-entrained particles by vertical velocity events of longitudinal vortices circle around the central vortex at least once more, adding to their residence time a turnover time scale on average every cycle they remain inside. The vortex dynamics inside the cavity influences the approach trajectory of particles approaching the interface. When they reach the dark blue zones in figure 12, they will likely exit the cavity immediately.

Figure 13. Sequence of Lagrangian transport across the interface, with particles coloured by residence time (a–c) and vertical velocity (d–f). The re-entrainment of particles at the interface can be identified by particle pockets that sit between two vortices of the shear layer, which are advected back in the cavity by vertical motions. This 3-D effect of streamwise vortices is the leading re-entrainment mechanism.

To summarize the variability of particle residence times we compute the spatial distribution of the standard deviation of local residence time at each point (in units of time) and the non-dimensional coefficient of variation, corresponding to the ratio between the standard deviation and the mean residence time. Since these statistics exhibit a less pronounced variability in horizontal planes, in figure 14 we show depth-averaged magnitudes. The spatial distribution observed in the standard deviation of residence times is influenced by the complex instantaneous dynamics of the shear layer, and connects to the physical mechanisms observed in figure 13. Large and small values of standard deviation are related to regions of small and large residence times, respectively. The red region of large standard deviation at the left in figure 14 coincides with an internal shear layer produced by the vortex core in the centre of the cavity and the return flow near the upstream wall, which divides particles that exit or remain inside, where the latter circle the large-scale vortex one or more times. The coefficient of variation distribution shows the influence of the region between vortices, where there is a large disparity of residence times with respect to the average due to the periodic re-entrainment. Due to the small residence times in these zones, the large coefficient of variation implies that particles reaching these positions oscillate at the interface one or more times to finally exit or remain inside the cavity, and the largest values are explained by this zigzagging at the interface. Particles starting in these regions have residence times lower than 30 or 40, which is less than the time it takes to a particle that start close to the interface and circling the cavity before reaching the interface again.

Figure 14. Depth-integrated statistics of the residence times in the cavity volume. (a) Standard deviation; and (b) coefficient of variation.

6. Finite-time Lyapunov exponents

The spatial distribution of residence times generated by the coherent structures reveals regions built by particles trajectories that organize the flow into ordered patterns, accumulating trajectories with similar time scales. To formalize the detection of these volumes inside the cavity and at the interface, we study the range of time scales that emerge from the Lagrangian transport, calculating the 3-D finite-time Lyapunov exponents (FTLE) of particle trajectories (Haller 2015).

The FTLE are used to measure the rate of separation of neighbouring trajectories over a specific time interval. They are commonly used to identify Lagrangian coherent structures in the flow and to provide a quantitative evaluation of the time scales embedded in the particle transport, which characterize the flow complexity. The technique consists of computing the local flow deformation gradient, through the calculation of the largest stretching of trajectories that begin from a close initial condition, and assumes that the growth of the separation evolves in an exponential fashion, integrated forward in time.

The magnitude of this separation is used to identify flow structures associated with time scales of trajectories. Larger values of FTLE indicate the rate at which nearby particles move apart. We initially consider the uniform spatial distribution of particles inside the cavity at $t=0$, and a varying final integration time $t=T$. The local deformation gradient is a $3\times 3$ matrix at each node of the Eulerian grid based on the derivatives of particle positions. The magnitude of the deformation is the largest eigenvalue $\lambda _{max}$ of the Lagrangian strain tensor $\pmb {\varLambda }$, computed over the time interval

(6.1)\begin{equation} \sigma_{t}^T(x,y,z) = \frac{1}{T} \ln \sqrt{\lambda_{max}(\pmb{\varLambda})}, \end{equation}

where $\sigma _{t}^T(x,y,z)$ is the FTLE at the position $(x,y,z)$ in the time interval $0< t< T$.

In figure 15 we show the 3-D FTLE for a dimensionless time $T=30$, where we can identify three regions inside the cavity. The largest values of FTLE equal to or larger than $\sigma _{t}^{30}=0.17$ can be observed at the interface and in the upstream boundary in blue, whereas the smallest values that are one order of magnitude smaller are encountered at the centre of the cavity in the red region, $\sigma _{t}^{30}=0.08$, which is consistent with the analysis of residence times. The largest values of FTLE in the blue region at the interface and upstream wall, correspond to zones where a significant number of particles leave the cavity and others remain inside, increasing the separation distance. The smallest values of FTLE at the centre of the cavity are due to the longer time scales within the central vortex core, where particles remain for longer periods and their separation evolves slowly. The computation of FTLE for these longer time scales shows that, in the dynamically rich turbulent conditions of the flow, all particle trajectories ultimately diverge from all the regions identified inside the cavity, as they have positive repulsive values of the $\sigma _{t}^{30}$ exponents. The FTLE regions, depicted with different colours in figure 15, provide a simplified skeleton of the overall Lagrangian dynamics inside the cavity that can be used for choosing different zones in multiple-region models that solve a different 1-D mass transport equation for each of these regions (Drost et al. 2014; Sandoval et al. 2019).

Figure 15. Three-dimensional FTLE show regions of the cavity associated with different time scales for an observation period of a non-dimensional time equal to $T=30$.

Similarly, in figure 16 we show horizontal $x\unicode{x2013}y$ planes at middle depth of FTLE for different values of the integral scale $T=1,5,10,30$. In the context of the present simulations, it is important to note that $T=30$ was shown to be better for visualizing the regions of different time scales. Since all particles escape the cavity, the diagrams saturate for larger values of $T$ as particles separate faster when they are transported in the main channel. In figure 16 we show that the largest FTLE values are related to the evolution of the shear layer at the interface from $T=1$ to $T=10$, where the central region between vortices pairs identified in the residence time animation can be distinguished. At $T=30$ not only does the interface have larger values of FTLE but also the upstream boundary of the cavity, which corresponds to the faster-moving region or outer flow that envelops the central core. The FTLE analysis delineates clearly the three regions in the cavity: the shear layer, the intermediate zone of rotating particles and the core with an exponent one order of magnitude smaller. In this last region, around ${35}\,\%$ of the total particles remain when the global concentration is $C=0.4$, occupying the central ${25}\,\%$ of the cavity volume, underscoring its effects on the suppression of global transport for long times.

Figure 16. Horizontal $x\unicode{x2013}y$ planes of repulsive FTLE at different final times $T=1,5,10,30$ at middle depth $z/h=0.5$; (a) $T=1$, (b) $T=5$, (c) $T=10$, (d) $T=30$.

7. Conclusions and future work

We performed LES of the flow past an open-channel lateral square cavity to investigate the mass transport mechanisms using Eulerian and Lagrangian approaches. The model reproduces quantitatively the mean flow and instantaneous features observed in the experimental measurements, with the predominant role of the shear layer in the flow dynamics, shedding periodic vortices that impinge on the downstream corner of the cavity, inducing the recirculating flow inside the volume.

The mass transport between the lateral cavity and the main channel was studied by solving an advection–diffusion equation for the continuum analysis, and the transport of passive particles through the implementation of coupled equations for the particle paths and the instantaneous velocity, incorporating the effects of unresolved turbulent scales with the addition of a stochastic term. In both Eulerian and Lagrangian simulations, the total mass is initially located inside the cavity, to track the evolution of the exchange with the main channel and the temporal mass decay inside the cavity volume.

Simulations showed that the interplay of the shear layer with the central core region is the dominant transport mechanism between the cavity and the main channel, as the central vortex inside the volume retains mass for longer periods of time, on average, while the faster vortices of the shear layer extract mass continuously to the main channel. Even though this interaction is a general description of the mass exchange mechanisms, specific details of transport have implications at larger scales, as revealed by the simulation results.

The high-resolution mass transport of the LES model provided new insights into the exchange process from a global perspective, identifying small-scale processes that influence global transport and informing large-scale models that are commonly implemented to predict contaminant transport in environmental flows. The spatially integrated data of the mass inside the cavity showed that the classical approach with a first-order equation fails to represent the transport evolution for averaged concentrations below 0.4. The upscaled fractional equation model captured better the emergent behaviour of mass transport, in which the exchange coefficient $k$ is maintained for the continuum and particle descriptions, while the order of the time derivative represents the temporal dependence of the mass exchange. This fractional model represented better the statistical distribution of long residence times, and it showed the same exponent for Eulerian and Lagrangian transport, $\alpha =0.96$, exhibiting consistent memory effects induced by the combined influence of the interacting turbulent coherent structures of the flow for both transport approaches.

The spatial distribution of residence times for Lagrangian transport in the cavity revealed a complex fluctuating 3-D field, where initial conditions that exhibit similar residence times are grouped in specific regions inside the cavity volume. Mean residence times are longer for particles that start near the central vortex core with slower flow velocities, while particles that only remain in the cavity for shorter times typically start close to the interface or the upstream boundary. The complex instantaneous dynamics of the internal interaction of the shear layer with the central core produces an enveloping region with a high variability of residence times. Particles in this region that do not leave the cavity immediately go around the central vortex and face the interface again, repeating the process after a time scale equal to the turnover time of this large vortical structure. The analysis of individual trajectories also showed that the model captures the remarkable behaviour of particles at the interface described by Engelen et al. (2021), with particles that exit the cavity and are re-entrained, zigzagging at the interface due to the dynamics of the shear layer. The simulations results show that most of the re-entrainment events occur by streamwise vortices that produce vertical velocity events and strong inflow to the cavity of fluid that is encountered between the main columnar vortices of the shear layer.

The FTLE technique was also used to further investigate the mixing regions inside the cavity, quantifying the separation rate of trajectories starting from all the regions in the volume. Through the FTLE tested for different integration periods, we identified three regions with separation exponents of different orders of magnitude: the slowest central core where particles remain trapped for longer periods of time, the enveloping region around the cavity and the shear layer at the interface. Calculations of the FTLE for different time scales were also capable of capturing the re-entrainment region between two subsequent vortices of the shear layer. This analysis links the spatial distribution of the FTLE with the memory effects revealed by the fractional transport equation, as ${35}\,\%$ of remaining particles are located in the central region when the global particle concentration is equal to 0.4. The simplified dynamics depicted by the FTLE analysis with different zones of similar time scales inside the cavity is important for 1-D multiple-region models of mass transport.

The Lagrangian results presented here contribute to identifying critical zones of mass storage and release, and the general spatial distributions of time scales in the cavity associated with the dynamics of coherent structures, providing new insights into the mechanisms that drive the global transport. This Lagrangian framework has the advantage of generating complete information of the particle dynamics without the risk of losing trajectories, as can occur experimentally with particle tracking velocimetry (PTV) (Engelen et al. 2021), improving the reliability and sensitivity analysis of mass transport, and complementing previous experimental and field observations.

Future research will focus on extending the analysis of the Lagrangian patterns to open-channel lateral cavities. We will seek to connect the topology of time scales to global parameters such as the Froude number (Karimpour et al. 2021), the Reynolds number and the aspect ratio of the cavity via numerical simulations. Exploring this multi-dimensional phase space will provide a better understanding of the description of global transport characteristics and potentially propose modifications to the geometry of the lateral recirculating regions that can optimize or control residence time statistics. We will also simulate the transport of fine sediments or inertial particles as contaminants that can be deposited or resuspended, investigating the effects of the processes studied in this work for different geometries and flow parameters at laboratory and field scales.