Subsurface flow mixing in coarse , braided river deposits

Introduction Conclusions References Tables Figures


Introduction
The subsurface heterogeneity at the 1 to 100 m scale can induce significant subsurface flow mixing that is relevant for aquifer remediation or drinking water extraction near a river or a contaminated area (e.g.Kitanidis, 1994;Mattle et al., 2001;Mays and Neupauer, 2012;Cirpka et al., 2015).Subsurface flow mixing is generally decomposed into an advective transport process combined with diffusion/dispersion (e.g.Mays and Neupauer, 2012).The advective transport process is best visualised with streamlines or streamtubes.Two-dimensional and three-dimensional flows exhibit a different streamline rearrangement when flowing through heterogeneities (Steward, 1998).Whereas two-dimensional, divergence-free flows locally deform the streamline geometry, three-dimensional, non-axisymetric flows permanently rearrange their streamtubes by redistributing the fluid within the subsurface (Steward, 1998;Janković et al., 2009).Janković et al. (2009) illustrated this difference by comparing two-dimensional and three-dimensional flows through an isolated, high-permeable subsurface structure, whose rotational axis was not aligned with the mean flow direction (i.e.non-axisymetric flows).For two-dimensional flows, the distance between the streamlines at a large distance upstream and downstream from the high-permeable structure remains the same.In contrast, the streamlines of threedimensional flows are permanently deformed downstream from the high-permeable subsurface structure resulting in a complex intertwining of streamlines.Janković et al. (2009) coined the phrase advective mixing to describe these phenomena.Cirpka et al. (2015) identified three advective mixing phenomena that enhance solute mixing: (i) streamline focusing/defocusing, (ii) depth-dependent streamline meandering (i.e.streamline deviation), and (iii) secondary motion consisting in persistent twisting, folding, and intertwining of streamlines.Chiogna et al. (2015) demonstrated the occurrence of macroscopic helical flow in subsurface flow simulations where the hydraulic conductivity field was heterogeneous and locally isotropic.Despite the importance of advective mixing in solute mixing processes that enhance diffusion/dispersion (Hemker et al., 2004;Janković et al., 2009;Cirpka et al., 2015;Ye et al., 2015), volumetric concentration measurements on the field cannot distinguish between advective mixing and dispersion/diffusion (Janković et al., 2009).Table 1.Hydraulic properties of the main sedimentary structures (after Jussel et al., 1994a).
Poorly sorted gravel Bimodal gravel Open-framework gravel Porosity 0.2 0.25 This study is part of a research project on the heterogeneity characterisation of coarse, braided river deposits on different scales.We focus on one important aspect of heterogeneity, namely its influence on advective mixing.Coarse, braided river deposits are highly heterogeneous in terms of hydraulic properties (e.g.Jussel et al., 1994a;Anderson et al., 1999;Lunt et al., 2004) and make up many groundwater reservoirs worldwide (Huggenberger and Aigner, 1999;Klingbeil et al., 1999;Bayer et al., 2011) and more than two-thirds of the exploited aquifers in Switzerland (Huggenberger, 1993).In this study the sedimentary heterogeneity is characterised following the hierarchy proposed by Huggenberger and Regli (2006).In order of increasing size, this hierarchy consists of sedimentary textures, sedimentary structures, and depositional elements.As schematically represented in Fig. 1, coarse, braided river deposits are composed of two main depositional elements, remnants of gravel sheets (Huber and Huggenberger, 2015) and trough fills with clear-cut erosional lower-bounding surfaces (e.g.Siegenthaler and Huggenberger, 1993;Jussel et al., 1994a;Beres et al., 1995Beres et al., , 1999;;Rauber et al., 1998;Stauffer and Rauber, 1998;Teutsch et al., 1998;Anderson et al., 1999;Klingbeil et al., 1999;Whittaker and Teutsch, 1999;Heinz and Aigner, 2003;Heinz et al., 2003;Huggenberger and Regli, 2006;Bayer et al., 2011).The sedimentary structure of the remnants of gravel sheets consists of horizontal to sub-horizontal layers with a poorly sorted gravel texture.The sedimentary structure of the fills generally consists of alternating openframework-bimodal gravel couplet cross-beds, although fills consisting of poorly sorted cross-beds or of interfingering cross-beds of poorly sorted gravel and sand are not uncommon (e.g.Siegenthaler and Huggenberger, 1993).Other less frequent sedimentary structures and depositional elements are described in the references above.Because the permeability contrast between the open-framework gravel texture and the other textures (bimodal gravel, poorly sorted gravel) is up to 3 orders of magnitude (e.g.Jussel et al., 1994a, Table 1), the spatial distribution of the open-framework gravel texture is expected to strongly influence the subsurface flow field and therefore to enhance advective mixing (Stauffer, 2007).
Based on observations of hydrofacies or sedimentary structures, several studies developed hydrogeological models of coarse, braided river deposits to investigate subsurface transport.Most of these studies assessed either macrodispersion processes (e.g.Jussel et al., 1994b;Stauffer and Rauber, 1998), sorption processes (e.g.Rauber et al., 1998;  Teutsch et al., 1998), or particle concentrations (e.g.Anderson et al., 1999;Heinz et al., 2003), mainly through the analysis of breakthrough curves.Stauffer (2007) modelled a trough fill of alternating open-framework-bimodal gravel couplets by a highly permeable rectangular cuboid with an anisotropic hydraulic conductivity tensor.He quantified the subsurface flow disturbance downstream of the cuboid embedded in a homogeneous background matrix as a function of the angle of anisotropy of the hydraulic conductivity tensor.He noticed that "the disturbance manifests itself by a distinct distortion of the streamtubes.Laterally, the influenced width is about 2.5 times the width of the [cuboid] for the considered case.Vertically, this influenced width makes up about 10 times the thickness of the [cuboid]" (Stauffer, 2007).
To the best of our knowledge the influence of trough fills on advective mixing has not been investigated with the exception of the work of Stauffer (2007) in which the complex trough fill structure was reduced to a simple cuboid with an homogeneous anisotropic conductivity.The present work aims to assess the influence of a geologically realistic representation of high-permeable trough fills on advective mixing.
The flow simulation is performed on a synthetic, conceptual model derived from ground-penetrating radar (GPR) data recorded over a small area (about 100 m × 50 m) of the river bed of the coarse, braided Tagliamento River (northeast Italy).First, the sedimentary structure of two overlapping trough fills is inferred from three GPR profiles, one 53 m long approximatively parallel to the main flow direction and two 7.5 and 10 m long approximatively perpendicular to the main flow direction.Simple geometric objects corresponding to each sedimentary structure are manually fitted to the in-terpreted GPR records.Then, a high-resolution, steady-state, three-dimensional groundwater model is set up based on hydraulic properties borrowed from the literature.Finally, advective mixing is investigated with particle tracking.

Ground-penetrating radar data acquisition
The objective of the project was to quantify the proportion of depositional elements in the sedimentary deposits.Because the erosional lower-bounding surfaces of trough-shaped depositional elements can be followed over large distances (> 25 m), 14 widely spaced GPR lines (about 25 m line spacing on average) were acquired in a 100 m × 200 m large area on the river bed of the coarse, braided Tagliamento River downstream from the Cimano bridge (46 • 12 37.945 N, 13 • 0 50.165E; WGS1984).The GPR data were recorded with a PulseEkko Pro GPR system (Sensors & Software Inc., Mississauga, Canada) with 100 MHz antennae.The nominal spatial resolution length of the 100 MHz antennae is of the order of 0.3 m (Bridge, 2009).The topography of the GPR profiles was surveyed with a Total Station.
The GPR data were processed as follows: -Time-zero adjustment.
-Direct current-offset (DC-offset) removal based on samples before time zero.
-Dewowing of each trace by removal of the trend estimated with a Hampel filter (Pearson, 2002).
-A spherical and exponential gain was applied to compensate for geometric spreading and attenuation (Kruse and Jol, 2003;Grimm et al., 2006).This gain preserves the relative amplitudes.
-Time-to-depth conversion with a constant velocity of 0.1 m ns −1 that leads to results that are sufficiently accurate for the purpose of this study.The velocity was estimated from previous common mid-point surveys recorded in the same area.

Ground-penetrating radar data interpretation
The interpretation of the GPR profiles is based on (i) continuity of the dominant reflectors within and between the profiles, (ii) differences of reflection patterns, and (iii) angular unconformity between the reflectors that can indicate an erosion surface or the superposition of two sedimentary structures with different sedimentary textures (Beres et al., 1995(Beres et al., , 1999)).Three GPR profiles image three relatively well-preserved, overlapping trough fill structures that are identified by their erosional lower-bounding surfaces.Figure 2 shows the three GPR profiles as well as their interpretation.The GPR data indicate that the trough fills are elongated in the main flow direction (i.e. the valley orientation) with cross-tangential reflector.The GPR profile "xline1" (perpendicular to the mean flow direction; Fig. 4a) displays asymmetrical circular-arced reflectors that are almost symmetrical on the profile "xline2".Most of the older trough (represented in green in Fig. 2) is eroded by the younger troughs (represented in blue and red in Fig. 2).

Subsurface structural modelling
The observed reflections are consistent with the results of many studies on coarse deposits that compared GPR reflections with sedimentological structures of outcrop exposures (e.g.Huggenberger, 1993;Bayer et al., 2011).Because only three GPR records image the trough fills, a conceptual representation of the sedimentary structure is needed to infer the three-dimensional structure of the imaged trough fills at a high resolution.The approach proposed by Siegenthaler and Huggenberger (1993) is adopted.Siegenthaler and Huggenberger (1993) hypothesised that trough fills originate from confluence scours that can migrate.Therefore, they sug-gested to simulate the internal structure of the trough fills based on geometric considerations, i.e. by several shifted half-ellipsoids representing the trough migration (see also Best and Rhoads, 2008).In this study, the trough fills are represented by truncated ellipsoids.The position and the size of several truncated ellipsoids was manually adjusted (i) to match the positions of the identified erosional lowerbounding surfaces and (ii) to respect the orientations of the internal structures of the trough fills that are visible on the GPR records.A top view of the resulting subsurface structural model is shown in Fig. 3.The GPR profiles are compared to vertical sections of the structural model as well as to vertical gravel pit exposures of coarse, braided river deposits located in northeast Switzerland (Fig. 4).

Hydrogeological model
The three-dimensional model grid has a size of 75 m×80 m× 9 m and a resolution of 0.5 m × 0.5 m × 0.1 m.The truncated ellipsoids are located between 0.6 and 3.1 m below the surface.Because of the close correspondence of the GPR reflection patterns and of the sorting process with the observations made by Siegenthaler and Huggenberger (1993), Huggenberger (1993), Beres et al. (1995Beres et al. ( , 1999)), and Heinz et al. (2003), we assume the hydraulic properties of the different types of gravel texture to be in the same order of mag- nitude as those estimated from measurements on disturbed and undisturbed samples in Quaternary coarse gravel deposits in northeast Switzerland (Jussel et al., 1994a).The hydraulic properties of the poorly sorted gravel (see Table 1) are attributed to the background matrix, while the hydraulic properties of the bimodal and open-framework gravel (Table 1) are alternatively assigned to the voxels located between two consecutive truncated ellipsoids, following the conceptual model shown in Fig. 1.For each voxel the hydraulic conductivities are drawn from log-normal distributions neglecting any spatial correlation (they are identically and independently distributed).The resulting conductivity field is displayed in Fig. 5a.The hydraulic conductivity tensors of the bimodal and open-framework gravel are both isotropic.A vertical anisotropy of the hydraulic conductivity (K h / K v = 6) is assigned to the poorly sorted gravel texture to reflect the layered structure that hinders vertical flow.
All model boundaries are set as a no-flow boundary with the exception of the inflow (x = 0 m) and outflow (x = 75 m) faces where constant head boundary conditions are specified (Fig. 5).The gradient between the inflow and the outflow model faces is 0.03 and corresponds to a locally large hydraulic gradient as found in situations where groundwatersurface water interactions occur.The saturated, steady-state subsurface flow simulation is performed with MODFLOW, the USGS three-dimensional finite-difference groundwater model (Harbaugh, 2005).

Advective mixing quantification
The advective flow is simulated with the particle-tracking code MODPATH (Pollock, 2012).One particle per cell is set on the model inflow face and the position of the particles travelling through the model is recorded.The resulting streamlines combined with a judicious colour scheme al-low for visualisation of advective mixing.Furthermore, we quantify advective mixing by evaluating (i) particle deviation, (ii) particle divergence, and (iii) particle intertwining between the inflow face and the outflow face.
The particle deviation ( ) is the transverse distance between the particle positions on the inflow face (y i , z i ) and on the outflow face (y o , z o ): For each cells of the outflow face, we compute the median particle deviation from all particles within the cell.
The particle divergence indicates how far a particle flowed away from its eight particle inflow neighbours.For each particle we compute the absolute difference between (i) the median distance between the particle and its eight neighbours on the inflow face and (ii) the median distance between the particle and its eight neighbours from the inflow face on the outflow face.
The neighbours of a particle on the inflow face can be different from the neighbours of the same particle on the outflow face.Therefore, the particle intertwining is estimated for each particle by the proportion of its four inflow that are still its neighbours on the outflow face.In order to include all neighbour particles, the neighbours on the outflow face are defined as the first-and second-order neighbours of the Delaunay triangles (i.e. the particles connected to the considered particle through an edge or two edges of the Delaunay triangles.3 Results and discussion

Hydraulic heads
Similar to a high-permeable homogeneous structure, the overlapping trough fills significantly influence the hydraulic head distribution -vertically (Fig. 5b) and horizontally (Fig. 6) -inducing an asymmetric flow focusing and defocusing (compare with Fig. 7). Figure 6 shows within longitudinal cross sections how the vertical distribution of the hydraulic heads is significantly influenced by the trough fills: the hydraulic gradient is oriented upward, toward the trough fills at their upstream end and downward, away from the trough fills at their downstream end.However, even in the centre of the model this pattern is never symmetric (Fig. 6b) because of (i) asymmetry of the internal structure of the trough fills and (ii) non-alignment of the trough fills with the mean flow direction.The asymmetry of the vertical hydraulic head distribution becomes more asymmetric close to the lateral model boundaries.The upward gradient upstream from the trough fills slowly disappears toward the right model boundary (looking downstream; Fig. 6a), while the downward gradient downstream from the trough fills slowly disappears toward the left model boundary (Fig. 6c).The hydraulic gradient within the trough fills is very small (about 0.002).The asymmetry of the three-dimensional hydraulic head distribution causes a permanent rearrangement of the streamlines.Therefore, in addition to a flow focusing and defocus-  ing effect, persistent streamline deformations and rearrangements are expected.

Particle tracking
Figure 8 shows the position of the particles on the model outflow face coloured by their initial y and z coordinates on the inflow face.The convex hull of the particles on the outflow face that flowed through the trough fills as well as the shape of the trough fills projected on the outflow face are also represented.The size of the projected trough fill shape and convex hull are about 38.5 m×2.2 m and 52.0 m×6.7 m, respectively.On the inflow face, the shape of the convex hull of the particles that flow through the trough fills (not shown) is up to a lateral shift of 8 m nearly identical to the convex hull shown in Fig. 8.This could indicate a similar flow focusing and defocusing effect combined with a lateral flow deviation.However, a notable particle deviation is clearly visible inside and outside the convex hull (see also Fig. 9).The median particle deviation is 4.0 m whereas the maximum is 28.1 m.The particle deviation outside the convex hull is very small Hydrol.Earth Syst. Sci., 20, 2035-2046, 2016 www.hydrol-earth-syst-sci.net/20/2035/2016/ with the exception of some particles below the convex hull (up to 12 m).Even if small, the particle deviation outside the convex hull is smoothly varying because these particles flowed through the low heterogeneous poorly sorted gravel.The largest particle deviations are observed within the convex hull.There, the particle deviations are irregular in amplitude and direction but still show an horizontal trend as expected from the orientation of the trough fills.Note that the asymmetry of the trough fills causes a partial, large-scale rotation of the particles.
The largest median distances between each particle and its eight inflow neighbours on the outflow face are found within the convex hull (Fig. 10a), where most of the particles lay at least 4 times farther away from their inflow neighbours as on the inflow face.The median distance between a particle and its eight neighbours is 0.1 m on the inflow face and less than 2 % of the particles are more than 10 m away from their neighbours.The largest distances are found in the central part of the convex hull that is associated to the two younger trough fills (trough fills 2 and 3 in Figs. 2 and 3).More than the half of the particles outside the convex hull lay closer to their inflow neighbours on the outflow face.The analysis of the remaining neighbours (Fig. 10b) attests for a strong particle intertwining as indicated by Fig. 10a.Indeed, about 70 % of the particles in the convex hull on the outflow face are no longer surrounded by their four initial neighbours from the inflow face.

Advective mixing mechanism
For clarity, Fig. 11 shows only a few particle paths that cross the trough fills.The particles upstream from the trough fills are attracted by the highly permeable layers of the open-framework gravel.Shortly before the particles enter the trough fills, some of them show a strongly curved path toward the trough fills.The particles that enter the open-framework gravel layers move horizontally within these layers until they dip upward.A closer look on Fig. 11 reveals series of sharp vertical zigzags of the particle paths, predominantly at the downstream end of the trough fills where the layers of openframework gravel dip upward.These zigzags occur where the particles tightly jump vertically between two adjacent layers of open-framework gravel.
Figure 12 displays an enlarged view of a vertical section of the model along the main flow direction that shows the layers of open-framework gravels as well as the vertical hydraulic head distribution.The arrows represent the volumetric flux (Darcy's flux) vectors projected on the vertical section for each cell of the open-framework layers.Note that the hydraulic conductivity tensor within the trough fills is isotropic.Therefore, the volumetric flux along each dimension of the Cartesian coordinate system is proportional to the hydraulic conductivity at the cell interface times the hydraulic gradient along the same dimension.Figure 12 reveals a complex spatial distribution of the volumetric flux that appears rather chaotic in the upward-dipping part.However, we observe that four of the upward-dipping layers of the openframework present a similar pattern: although very small in amplitude, the volumetric flux of the lower cells of these layers tend to point downward whereas in the upper cells the flux tend to point upward.The vertical position of the particles within the open-framework gravel layers is therefore critical because two closely spaced particles can flow in opposite direction.As a consequence, the volumetric flux pointing downward lets some of the ascending particles exit the trough fill earlier (see Fig. 11).In a similar way, two closely spaced particles do not enter the trough fills at the same position and therefore follow different paths within the trough fills.Small spatial variations of the volumetric flux (not only vertically but also horizontally) can drive the particles far away from each other (Fig. 11).This advective mixing illustrates the importance of the interplay between the hydraulic head field and the spatially distributed hydraulic conductivity that results in an heterogeneous volumetric flux distribution within the trough fills.
In consequence, the transport process through the trough fills can be viewed as a chaotic process where the particle positions on the outflow face depends on the initial particle positions on the inflow face (Neupauer et al., 2014).Note that the same effect is obtained with homogeneous hydraulic conductivity for each sedimentary texture.Spatial random hydraulic conductivity values increase advective mixing at a level that is negligible compared with the advective mixing resulting from the three-dimensional arrangement of the different textures.
Investigation on the influence of hydraulic gradient, hydraulic conductivity of the open-framework gravel, vertical anisotropy, and trough fill orientation on advective mixing showed the following: (i) the decrease of the hydraulic gradient significantly increases the lateral deviation of the particles; (ii) the extent of the convex hull of the particles that crossed the trough fills, the particle deviation, and mixing increase with increasing hydraulic conductivity of the openframework gravel; (iii) the vertical extent of the convex hull zone downstream of the trough fills as well as the vertical particle deviation are inversely proportional to the vertical anisotropy (K h / K v ) of the poorly sorted gravel texture (matrix) because a large vertical anisotropy of the poorly sorted gravel texture hampers vertical flow; (iv) the angle between the trough fills and the main flow direction plays an important role for the mixing processes.The width and height of the mixing zones negatively correlate when the orientation of the trough fills changes impacting significantly advective mixing.Furthermore, when the trough fills are aligned with the main flow direction, a partial, transverse rotation of the particles is observed within the convex hull.When the trough fills are perpendicular to the main flow direction, the advective mixing is the smallest.The largest convex hull, particle deviation and mixing are found when the trough fills form an 45 • angle with the main flow direction.

Discussion
Advective mixing is enhanced by the spatial distribution of trough fills in the sedimentary record and by the unsteady flow magnitude and direction.The advective mixing zones of closely spaced trough fills can interfere, resulting in a more complex pattern of subsurface flow.Under transient boundary conditions the mean flow direction and therefore the angle between the trough fills and the main flow direction change with time.In such a situation, the advective mixing zone as well as the flow patterns are expected to vary spatially and temporally, most likely leading to enhanced advective mixing.Because of this complexity, the present experiment is a starting point for further investigations into the influence of different proportions and types of trough fills on advective mixing in coarse fluvial aquifers on the 1 to 100 m scale.
In the presented synthetic model, the layers of poorly sorted gravel are modelled by an uniform, anisotropic matrix because the interface between the layers of poorly sorted gravel are barely identifiable on the GPR records.While the model set-up (isolated trough fills embedded in poorly sorted gravel) was observed in gravel quarries (e.g.Siegenthaler and Huggenberger, 1993), thin, finite layers of open-framework gravel can also be found within the layers of poorly sorted gravel (e.g.Huggenberger and Regli, 2006).However, the contribution of these thin, high-permeable structures to advective mixing is expected to be negligible compared to that of the trough fills.The hydraulic conductivity tensors of the bimodal and open-framework gravel are both isotropic.However, upon upscaling, the open-framework-bimodal gravel couplets lead to an anisotropic hydraulic-conductivity tensor (e.g.Jussel et al., 1994a;Stauffer, 2007) with the thicknessweighted arithmetic mean of the two conductivities within the layering and the thickness-weighted harmonic mean perpendicular to it.Therefore, on larger scales, the flow direction may not be parallel to the hydraulic head gradient.
Note that the use of an interpolation scheme is superfluous if densely sampled GPR data are available (e.g.pseudo three dimensional GPR survey) and the different sedimentary textures are well-resolved by GPR.

Conclusions
In this study, the hydraulic heterogeneity of coarse, braided river deposits is modelled through a simple geometrical model based on geological observations.The modelled trough fills (i) act as an attractor for the groundwater upstream of the trough fills, (ii) induce a significant intertwining of the streamlines that flow through it, resulting in strong advective mixing, and (iii) cause a strong horizontal streamline deviation that results in a partial, large-scale flow rotation.Furthermore, the anisotropy of the hydraulic conductivity of the poorly sorted gravel strongly influences vertical advective mixing whereas the orientation of the trough fills determine the flow patterns and therefore the degree of mixing.The advective mixing produced by the trough fills resembles a chaotic process that is very sensitive to the initial positions of the streamlines.While the emphasis is often put on the fast flow pathways and their connectivity, this study demonstrates the importance of the sedimentary structure of the whole geological fabrics in interaction with the hydraulic boundary conditions (see also Voss, 2011).
This study is only valid for the considered type of trough fills, i.e. trough fills consisting of layers of bimodal and openframework gravel, and for the proposed conceptual model.Trough fills consisting of cross-bedded poorly sorted gravel or of interfingering cross-beds are very likely to lead to different flow structures and therefore to different mixing patterns.The subsurface structure could more accurately be modelled with high-resolution GPR data, thereby making the use of the geometrical model unnecessary.
The study findings shed light on possible advective mixing in natural environments and indicate complex advective mixing in dynamic systems such as in systems characterised by significant groundwater-surface water interactions.A better understanding of the sedimentary structure can provide an additional support to the interpretation of the ecological processes in the hyporheic zone.

Figure 1 .
Figure 1.Simplified conceptual model of a single trough fill (with alternating open-framework-bimodal gravel couplets) embedded into layers of poorly sorted gravel.

FFigure 2 .
Figure 2. Fence diagram of the GPR data and their interpretation.The black arrows indicate the GPR survey direction.

Figure 3 .
Figure 3. Top view of the geometrical trough fill model (Coordinate system: WGS 1984, UTM Zone 33N).The trough fills are represented by green, blue, and red ellipses.The black lines indicate the position of the ground-penetrating radar profiles and the black arrows the GPR survey direction.

Figure 4 .
Figure 4. (a)-(c) Ground-penetrating radar data, sections of the geometric model and vertical outcrop exposures (northeast Switzerland) for comparison purposes.The trough fills are represented by green, blue, and red ellipses.The black arrows indicate the GPR survey direction

Figure 5 .
Figure 5. (a) Hydrogeological model set-up with spatial distribution of the hydraulic conductivity values.(b) Hydraulic head at the upper model boundary (top view, contour every 0.2 m).The blue arrows indicate the main flow direction.

Figure 6 .
Figure 6.Cross sections of the hydrogeological model along the x axis (see the coordinate system defined in Fig. 5a) with hydraulic head contours (every 0.2 m) superimposed on the hydraulic head values.The blue arrows indicate the main flow direction.The grey pixels correspond to the highly permeable layers of open-framework gravels.

Figure 7 .Figure 8 .
Figure 7. Particles coloured by their (a) y coordinate position and (b) z coordinate position on the inflow face.The blue arrows indicate the main flow direction.

Figure 9 .
Figure 9. Median particle deviation between the inflow face and the outflow face (computed vertically for every five cells) represented by arrows.The arrow length and colour correspond to the deviation magnitude.The black line represents the shape of the trough fills projected on the outflow face and the dashed, red line represents convex hull of the particles on the outflow face that flowed through the trough fills.The blue arrow indicates the main flow direction.

Figure 10 .
Figure 10.Particles on the model outflow face.(a) Median distance between each particle and its eight inflow-face neighbours computed on the outflow face.(b) For each particle on the outflow face, number of remaining neighbours from their four inflow-face neighbours.The black line represents the shape of the trough fills projected on the outflow face and the dashed line represents convex hull of the particles on the outflow face that flowed through the trough fills.The blue arrow indicates the main flow direction.

Figure 11 .
Figure 11.Selected particles coloured by their (a) y coordinate position and (b) z coordinate position on the inflow face.The blue arrow indicates the main flow direction.

Figure 12 .
Figure 12.Enlarged view of the vertical section of the hydrogeological model along the x axis with the hydraulic head contours (every 0.01 m) superimposed on the hydraulic head values.The grey rectangles represent the open-framework cells.The arrows correspond to the volumetric flux vectors projected on the model section; red indicates that the flux flows downward, blue upward.The large blue arrow on top indicates the main flow direction