Turbulent flow dynamics and mass transport processes in a natural surface storage zone using field data and numerical simulations

In this work, the turbulent flow dynamics and mass transport mechanisms in a natural SSZis analyzed. The study site is a river reach of the Lluta River, located in northern Chile in a high-altitude Andean environment known as the Altiplano (∼ 4,000 masl) The large-scale turbulent coherent structures are characterized using field measurements and 3D numerical simulations. The detailed topography was measured through DGPS and digital image processing while the surface velocity field, through the LSPIV technique. Regarding the field data, numerical simulations were performed using a DES turbulence model coupled with a 3D passive scalar transport model for Re = 45,800. The coherent structure dynamics in the shear layer was identified as the main mechanism that drives the mass and momentum transport processes between the SSZ and the main channel. Also, the 2D vortical structures of the mean flow are analyzed within the lateral cavity, since they have a strong influence in mass transport, increasing mean residence times due to their lower velocities and longer exchange timescales. Finally, the performance of two simplified transport models is analyzed to represent the mass transport dynamics at larger scales.


Introduction
In general, stream channels can be separated into advective (main channel) and non-advective (transient storage) zones [5].Transient or surface storage zones (SSZ) are characterized by low velocities compared to the flow in the main channel, and long residence times that favor the deposition of contaminants, nutrient uptake and interactions with reactive sediments.For its role in mass transport processes, biogeochemical reactions and the habitats formation for aquatic communities, lateral cavities are one of the most important SSZ in riverine environments.The dynamics of the turbulent flows in these zones is very complex, typically characterized by a shear layer that induces a recirculating area, with multiple large-scale coherent structures of different temporal and spatial scales.Flow dynamics between the main channel and the SSZ has been widely studied both numerically and experimentally [e.g., [6][7][8][9][10][11][12].Nevertheless, limited works have determined and analyzed the turbulent flow features of natural SSZ with high Reynolds numbers [4], where the interaction between turbulent flow, solute transport and intricate stream bathymetry represent a very complex problem [13].
The objective of this work is to further the undestanding of the transport processes and the flow patterns developed in a SSZ in natural conditions.In this way, this work seeks to determine and analyze the main mechanisms through the turbulent flow leads the transport dyamics of mass and momentum between the main channel and the lateral cavity.Finally, this work seeks to elucidate some insights about the performance of some empirical relations to determine key indicators of mass transport dynamics.

Study site and field characterization
The study site is a river reach located approximately 1000 m downstream of the confluence between Azufre and Caracarani Rivers.There, a lateral cavity is formed in the downstream side of a boulder located in the left riverbank.The presence of a smaller boulder in the stream, combined with the riverbank boulder generates a narrowing section in the flow.The interaction between the flow and the boulders generate a combination of large scale structures that emerge from the narrowing section and are advected by the flow.A complex surface flow was observed in the field, characterized by a large recirculation area and the presence of a sink point located in the leftbank downstream the narrowing section, where the flow is divided in opposite directions: one part enters the cavity and the other one moves onto the main channel in streamwise direction (see Figure 1.a).Two methodologies were applied in order to obtain field data: Large-Scale Particle Image Velocimetry (LSPIV) to obtain detailed surface mean velocity field, and image analysis combined with DGPS survey to obtain three-dimensional bed topography.LSPIV is a technique, which by means of statistical analysis of consecutives images, seeks to determine the surface velocity field by the identification of patterns formed by the particles (most of the time, pieces of paper) that are moved by the flow.Particle displacements are calculated by computing two-dimensional correlation (or cross-correlation) on pairs of orthorectified images, separated by a given ∆t [14] and filtered by a correlation threshold.For this study, seven movies, of approximately 2 minutes each, were recorded in 1 hour's period with a sample rate of 30 frames per second (∆t = 33.3ms).To obtain the detailed bathymetry for numerical simulations, a bed topography survey was carried out.A set of photos of the dry bed were taken and ten points were obtained using a Digital Global Positioning System (DGPS) device.The field data was processed through Structure from Motion (SfM) photogrammetry technique using the software Photoscan.In Figure 1.c it is possible to see the Digital Elevation Model (DEM) obtained as result of image analysis process.Also, the discharge and water depth were measured using a handheld acoustic Doppler velocimeter (Flowtracker Handheld ADV, Sontek, San Diego, CA, USA) in a gauging section located at 2 meters downstream from the cavity.The discharge in the main channel was Q = 0.2683 m 3 /s.The mean water depth of the gauging section was h = 0.097 m with a channel width of b = 5.88 m.Thus, the bulk velocity in the main channel is therefore equal to U b = Q/bh = 0.472 m/s.With this information, and considering that in the main channel b >> h, the bulk Reynolds number (U b h/ν) is about 45800 and the Froude number (U b / gh) corresponds to 0.48 (subcritical flow).

Three-dimensional numerical simulations and computational details
Three-dimensional numerical simulations were carried out to study some important features about turbulent flow dynamics in the shear layer and within the cavity, as well as dynamics of mass transport between the main channel and the SSZ.The governing equations for the resolved flow are the 3-D unsteady incompressible Reynolds -averaged Navier -Stokes equations (URANS) for the conservation of mass and momentum.For the mass transport, the volume fraction of conservative solute (C) is modeled through an advection-diffusion equation coupled to the flow solver.In tensor notation these equation can be written as follows: For the turbulence model, a hybrid URANS/LES formulation known as Detached-Eddy Simulation (DES) was appliedd [15,16].This approach is based on the one-equation eddyviscosity model closure for the URANS equations [17].For details of the model and the flow solver description, refer to [18].Computational domain was obtained using DEM of the bathymetry generated from image analysis described above.The computational domain was discretized in two zones using domain decomposition with overset grids.The first zone corresponds to the main channel and the second zone corresponds to the lateral cavity, with 5.4x10 6 and 1.8x10 6 nodes respectively.

Simplified transport models
For understanding the mass transport dynamics and the exchange processes between the main channel and the cavity, it is also important to know how the concentration of contaminants in the cavity evolves over time.[19] proposed the following one-dimensional transport model for surface storage zones: Assuming a first order exchange rate for all of these processes the conservation equations for each of these regions can be written as follows: In the equations, C 1 stands for to the concentration in the outer region of the cavity, C 2 is the concentration of the inner region in the cavity, β is the exchange rate between the inner and outer region, α is the exchange rate from the outer region to the main channel, φ = V 1 /V 2 is the ratio of the volumes of the two regions, V 1 is the volume of the outer vortex and V 2 is the volume of the inner region.The advantage of using a two-storage model is the capability to represent in a coupled way different exchange rates timescales in the cavity.This is an important feature of the model because the exchange timescales of the inner region are greatly slower compared with the outer region.

Velocity Field and mean flow patterns description
Some of the most important flow structures can be analyzed trough the mean flow field.The presence of 2D vortical structures in the SSZ defines zones of the cavity with low velocities within the recirculation region which are particularly important for the mass transport between the main channel and the SSZ and its residence time in the cavity [20,21].These flow patterns can be detected by the averaging of the filtered set of photos through the LSPIV technique.In Figure 3 the most important flow patterns identified in this study case can be seen.At point p1, the primary gyre of the recirculation region of the cavity can be observed; at point p2, the upwelling structure which evidences the three-dimensional behavior of the flow in the cavity is seen; at point p3, the impingement point of the cavity is clearly identified, at the flow bifurcation at the trailing edge.Finally, a secondary gyre can be identified, which appears as a counter-rotating structure associated to the primary gyre recirculation pattern.

Shear layer and mass exchange dynamics
Mass transport dynamics obtained from the numerical simulations revealed the importance of the shear layer as the main mechanism that drives the mass exchange between the main channel and the cavity.This phenomenon is produced by the large-scale coherent structures that are continuously shed from the mixing layer and advected along the interface, which produce periodic alternation of positive and negative fluctuations of transverse velocity [11].
Results also exhibit great differences in C evolution for different zones within the cavity (see Figure 4.a).These differences are mainly related to the mean flow vortical structures, where the zone influenced by the primary gyre shows C evolution rates much lower than the rest of the cavity, generating the existence of two storage zones interacting between them.The difference in exchange timescales in both zones can be observed more clearly in Figure 4.a.v.
In this graphic, the instantaneous time series of computed concentration evolution is plotted at two points within the cavity.The first one represents a point located at the center of primary gyre, while the second one represents a point located near the reattachment point at the end of the shear layer.As discussed before, there are differences in time scales of concentration decaying in different zones of the cavity related to the flow mechanisms that lead the fluid motion on them.

Simplified mass transport models performance
The mass exchange was quantified by computing the non-dimensional exchange coefficient proposed by [22].The velocity exchange scale was obtained using the approach followed by [23] and later by [24]: The value for k is obtained as the ratio between velocity exchange scale and bulk velocity, k = E/U b .In this case, the mass exchange coefficient yields a value of k = 0.0298, which agrees well with the values reported in literature for cavity flows.
Finally, the time evolution of mean or space-averaged concentration within the cavity was computed In Figure 4.b, two transient storage models were adjusted: single transient storage model (STSM) and two transient storage model (TTSM).The results show that the time evolution of mean concentration within the cavity exhibit a non-ideal behavior with a change of the decaying rate slope approximately between 20 seconds and 40 seconds.This slope change is associated to the time when the contaminant in the zone surrounding to the primary gyre has been transported significantly and the contaminant existing in the influence zone of primary gyre, which has exchange rates considerably lower than the surrounding one, begins to be transported at a higher rate.The models adjusted show that despite the STSM begins with a good agreement with the simulated data, as time progresses its deviation increases.Since its formulation, this model is also not able to represent the change of decaying rate slope.For its part, TTSM exhibits a good fit both at the beginning and at the end of the simulation, capturing, in a global sense, the change of decaying rate slope of mean concentration.Nevertheless, it is not capable to represent adequately the detail of concentration evolution in the middle period.Despite this, TTSM shows a much better performance than STSM with a RMSE considerably slower (0.017 in the TTSM versus 0.0556 in the STSM).

Conclusions
Through the combination of different measurement techniques, it was possible to obtain detailed information about the hydrodynamic features of the flow in natural conditions.The topographical survey, using photogrammetry, provided a high-resolution bathymetry of the bottom in the zone of interest.This technique was widely useful in the complementation of other techniques, specially providing the necessary information to generate the grid used in  the numerical simulations of the flow.For its part, the LSPIV technique provided the necessary information of the mean flow for the identification of the main coherent structures which govern some of the most important transport processes between the main channel and the cavity, and the exchange timescales of solutes within it.The numerical simulations allow to integrate spatial and temporal scales at different resolutions.This provides the necessary information for a more detailed analysis of the flow dynamics.Particularly, numerical simulations allow the identification and analysis of the spatial distribution and the evolution of the most important coherent structures that lead the mass and momentum transport between the main channel and the cavity.Thanks to the visualization of solute concentration evolution, it was possible to identify the main mechanisms through which the shear layer dominates the exchange processes.The unsteady behavior of the shear layer, specially at the trailing edge of the cavity, generates high transport rates in its influence zone.In contrast, the center of the gyres exhibits the lowest transport rates, representing regions with longer residence times.These differences in transport rates can be seen at higher spatial scales in the time evolution of the mean concentration in the whole cavity, where the TTSM represents a better adjustment than STSM for the mass exchange processes description.
In future research, these 1D descriptions will be improved by calculating the parameters of the models from geometric characteristics of the SSZ and flow parameters.Another feature that also will be investigated is the efficiency of contaminant transport and downstream effects of SSZ in series along the river channel, as well as the consequences of contaminant transport produced by the SSZ during extreme flood events

Figure 1 .
Figure 1.a) Schematic representation of the main surface flow feature of the SSZ studied.Letter "A" will represent the riverbank boulder and letter "B" will represent the in-stream boulder hereinafter for all figures.b) Orthorectified image obtained from image analysis with Photoscan.c) DEM generated from the orthorectified image.The resolution of each cell is approximately 1 cm

Figure 2 .
Figure 2. a) Top-view of the computational domain used for numerical simulations.b) Bottom overset grid.The zones of interest as the interface between the main channel and the cavity, or the zone where the reatachment produces, have finner grid distribution

Figure 3 .
Figure 3. Mean flow field.A, B, C and D represent the main flow structures identified by LSPIV.a) Orthorectified image of the study zone.b) Mean flow field at the water surface obtained from LSPIV measurements.The R 2 threshold used to filter correlation between consecutve images was 0.85

Figure 4 .
Figure 4. Simplified transport models comparison for mean concentration evolution within the cavity