Inﬂuence of Streambed Heterogeneity on Hyporheic Flow and Sorptive Solute Transport

: The subsurface region where river water and groundwater actively mix (the hyporheic zone) plays an important role in conservative and reactive solute transport along rivers. Deposits of high-conductivity ( K ) sediments along rivers can strongly control hyporheic processes by channeling ﬂow along preferential ﬂow paths wherever they intersect the channel boundary. Our goal is to understand how sediment heterogeneity inﬂuences conservative and sorptive solute transport within hyporheic zones containing high- and low- K sediment facies types. The sedimentary architecture of high- K facies is modeled using commonly observed characteristics (e.g., volume proportion and mean length), and their spatial connectivity is quantiﬁed to evaluate its e ﬀ ect on hyporheic mixing dynamics. Numerical simulations incorporate physical and chemical heterogeneity by representing spatial variability in both K and in the sediment sorption distribution coe ﬃ cient ( K d ). Sediment heterogeneity signiﬁcantly enhances hyporheic exchange and skews solute breakthrough behavior, while in homogeneous sediments, interfacial ﬂux and solute transport are instead controlled by geomorphology and local-scale riverbed topographies. The hyporheic zone is compressed in sediments with high sorptive capacity, which limits solute interactions to only a small portion of the sedimentary architecture and thus increases retention. Our results have practical implications for groundwater quality, including remediation strategies for contaminants of emerging concern.


Introduction
River bedforms enhance the exchange of water and solutes between the surface water and shallow riverbed, known as hyporheic exchange. Hyporheic exchange is typically characterized by complex flow processes at nested spatial scales, ranging from ripples and dunes to pool-riffles and gravel bars to meanders [1][2][3]. The architectural and textural characteristics of sediment facies (i.e., three-dimensional (3D) sediment bodies [4]) across different scales control the distribution of both the physical (e.g., hydraulic conductivity (K), porosity, grain size, and bulk density) and chemical (e.g., mineralogical composition, equilibrium sorption distribution coefficient (K d ), and reactive surface area) properties of the subsurface [5][6][7][8][9][10][11][12][13][14][15]. These properties in turn influence fluid flow and solute mass transport. Sorption processes may also hinder solute migration and increase the interfacial area of concentration gradients, which indirectly enhances diffusion and mixing [16,17]. Roughly 88% of whole-stream ecosystem respiration has been attributed to processing within the hyporheic zone (HZ) [18].
Numerous studies have shown that sediment heterogeneity can have a significant impact on hyporheic exchange fluxes, residence times (RTs), and flow patterns [19][20][21]. These studies mainly represented riverbed heterogeneity using geostatistical approaches, which rely on generating spatial variability in conductivity without incorporating the underlying sedimentary architecture. Using smooth, continuous conductivity fields in flume and numerical experiments, Salehin [22] found that heterogeneity can create preferential flow paths and increase hyporheic exchange. Sawyer and Cardenas [23] showed that cross-bedded conductivity structures slightly increase the depth of hyporheic exchange, while Liu and Chui [24] demonstrated that heterogeneity decreases mean and median hyporheic RTs while increasing the range of the RTs. Abrupt, discontinuous conductivity transitions are common in fluvial deposits composed of multiple sediment facies types, however, with potentially dramatic effects on hyporheic exchange [25]. To more accurately replicate field conditions, prior studies have represented the riverbed as a bimodal distribution of high-and low-K sediments (e.g., [4,26,27]). Pryshlak et al. [26] showed that the more abrupt the K transitions, the more impactful they can be. Several prior studies [28][29][30] also show that relatively low-K zones can be as important as high-K zones in driving exchange and influencing the spatial patterns of exchange. Stonedahl et al. [31] explored the influence of discontinuous conductivity fields on hyporheic exchange and RT distributions in a sand and sandy gravel framework, and found that there is a strong correlation between the volumetric flow rate and the arrangement of high conductivity grid cells. Using numerical models, Zhou et al. [4] quantified the influence of high-K open-framework gravel (OFG) strata on hyporheic flow and showed that such strata have a greater influence than the channel morphology (also see Pryshlak et al. [26]). Thus, in strongly heterogeneous sediments, facies architecture is often more influential than channel morphology in predicting hyporheic exchange and RTs.
Hydraulically connected, high-K facies create preferential flow pathways that significantly influence fluid flow and solute transport processes [32][33][34]. In general, connectivity increases with the volume proportion of high-K facies [35][36][37]. For example, Kennedy et al. [38] found that 70% of nitrate flux through the hyporheic zone of a small stream in an agricultural watershed occurred across only 38% of the streambed area, related to the intersection of high-K cross-strata with the channel boundary. These findings echo those of Zhou et al. [4], who determined that the majority of hyporheic exchange occurs at locations of high-K OFG cross-strata. Zinn and Harvey [39] observed a fundamental shift in transport behavior toward rate-limited mass transfer (i.e., heavy tailing) as the connectivity of high-K strata increased. A number of previous studies have investigated the effects of bimodal sediment heterogeneity on hyporheic exchange rates and solute RTs (e.g., [4,26,27,31]), but the impacts of low connectivity and sorption distribution have not been determined.
Chemical heterogeneity may also influence RTs and biogeochemical reaction rates [40][41][42][43][44]. A number of reactive contaminants may be temporarily removed through adsorption to particulate matter but are eventually released back into the water column depending on the sorption characteristics of the sediments [45]. Peter et al. [46] found that sorption to immobile zones, biofilms, and organic matter is a controlling mechanism in emerging contaminant removal in the HZ and results in longer hyporheic retention times with improved contaminant attenuation. Both positive and negative cross-correlation between K and K d can significantly affect the values of the effective retardation factor and reactive solute transport in the temporal and spatial domains [47,48].
The goal of this study is to determine how the spatial organization of sedimentary facies types, their spatial connectivity, and the resulting heterogeneity in K and in K d influence hyporheic exchange and solute transport. Specifically, we focus on the differences between fully connected (i.e., percolating) high-K facies and disconnected (i.e., non-percolating) high-K facies. A series of numerical flow and mass transport models of bimodal sediments were developed within the HZ of a gravelly meander with local-scale topography in the riverbed and were used to understand the fate of conservative and sorptive solutes along hyporheic flow paths. The sorptive capacity of the high-K sediments was varied to investigate the influence of chemical heterogeneity on sorptive solute transport. Models of homogeneous sediments with the same bulk composition were compared to further understand how sediment heterogeneity influences mass transport.

Hyporheic Flow Model
The pressure over the streambed influences the water and solute exchange across the stream-sediment interface, which consequently governs biogeochemical transformations in the hyporheic zone. Surface water flow over complex streambed morphologies and its resulting pressure can be simulated using computational fluid dynamics (CFD) solutions [49][50][51][52][53][54]. Our numerical models follow the framework of Stonedahl et al. [3] and Zhou et al. [4]. Specifically, Stonedahl et al. [3] created a three-dimensional (3D), multi-scale model for hyporheic exchange within a meandering stream. They used bedform topography quantified from a flume experiment to calculate the pressure over the streambed interface with an approach adapted from empirical relationships defined in [55][56][57]. Pore water flow was simulated using the finite difference method in MODFLOW [58], with the calculation of bedform-induced head fluctuations over the streambed represented as the top boundary condition on the subsurface flow model. Stonedahl et al. [3] evaluated the influence of morphology by running simulations at varying scales and quantified differences in the interfacial flux. Zhou et al. [4] modified the finite-difference model in [3] by renormalizing it to the scale of a 100 m wide channel meander in the well-studied Sagavanirktok River, and included a representation of OFG cross-strata in the streambed.
Here, we build on the model from Zhou et al. [4]. Full details of the model are available in Zhou et al. [4], but below, we provide a brief summary of the key model features and details about new additions to the model. As in Zhou et al. [4], the model domain is 800 m long, 400 m wide, and 3 m deep (see Figure 1; also see Figure 3 in [4]). We used a channel width of 144 m (highlighted in blue in the conceptual domain shown in Figure 1), which reflects the average channel width in the Sagavanirktok River system [59]. The model was uniformly discretized with a 2 m horizontal and 0.05 m vertical resolution (4.8 × 10 6 MODFLOW grid cells). We made sure the grid was fine enough to effectively capture the size and connectivity of dipping stratasets. To keep the problem computationally tractable, only one meander from Stonedahl et al. [3] was included in the simulation. The boundary conditions and hydraulic head distribution from Zhou et al. [4] were adopted here. The prescribed head values along the top boundary were obtained from Stonedahl et al. [3], preserving the overall hydraulic gradient (1.3 × 10 −3 ) and local variations in the head associated with small-scale bedform topography and channel morphology. Note that the gradient matches the slope of the energy grade line in the gravelly rivers, as reported by [59]. The prescribed head was also applied to the upstream and downstream boundaries. The model base and sides were defined as no-flow boundaries.
Water 2020, 12, x FOR PEER REVIEW 3 of 15 transport. Models of homogeneous sediments with the same bulk composition were compared to further understand how sediment heterogeneity influences mass transport.

Hyporheic Flow Model
The pressure over the streambed influences the water and solute exchange across the streamsediment interface, which consequently governs biogeochemical transformations in the hyporheic zone. Surface water flow over complex streambed morphologies and its resulting pressure can be simulated using computational fluid dynamics (CFD) solutions [49][50][51][52][53][54]. Our numerical models follow the framework of Stonedahl et al. [3] and Zhou et al. [4]. Specifically, Stonedahl et al. [3] created a three-dimensional (3D), multi-scale model for hyporheic exchange within a meandering stream. They used bedform topography quantified from a flume experiment to calculate the pressure over the streambed interface with an approach adapted from empirical relationships defined in [55][56][57]. Pore water flow was simulated using the finite difference method in MODFLOW [58], with the calculation of bedform-induced head fluctuations over the streambed represented as the top boundary condition on the subsurface flow model. Stonedahl et al. [3] evaluated the influence of morphology by running simulations at varying scales and quantified differences in the interfacial flux. Zhou et al. [4] modified the finite-difference model in [3] by renormalizing it to the scale of a 100 m wide channel meander in the well-studied Sagavanirktok River, and included a representation of OFG cross-strata in the streambed.
Here, we build on the model from Zhou et al. [4]. Full details of the model are available in Zhou et al. [4], but below, we provide a brief summary of the key model features and details about new additions to the model. As in Zhou et al. [4], the model domain is 800 meters long, 400 meters wide, and 3 meters deep (see Figure 1; also see Figure 3 in [4]). We used a channel width of 144 m (highlighted in blue in the conceptual domain shown in Figure 1), which reflects the average channel width in the Sagavanirktok River system [59]. The model was uniformly discretized with a 2 m horizontal and 0.05 m vertical resolution (4.8 × 10 6 MODFLOW grid cells). We made sure the grid was fine enough to effectively capture the size and connectivity of dipping stratasets. To keep the problem computationally tractable, only one meander from Stonedahl et al. [3] was included in the simulation. The boundary conditions and hydraulic head distribution from Zhou et al. [4] were adopted here. The prescribed head values along the top boundary were obtained from Stonedahl et al. [3], preserving the overall hydraulic gradient (1.3 × 10 -3 ) and local variations in the head associated with small-scale bedform topography and channel morphology. Note that the gradient matches the slope of the energy grade line in the gravelly rivers, as reported by [59]. The prescribed head was also applied to the upstream and downstream boundaries. The model base and sides were defined as noflow boundaries.   As in Zhou et al. [4], we used the bimodal sediment heterogeneity distribution of [4], which was simulated using a transition-probability (i.e., the probability of transitioning from one facies to another)-based approach in T-PROGS [60] in two stages. First, a Markov Chain (MC) model and indicator simulation with quenching (the T-PROGS code [60]) were used to simulate the structure of high-and low-K facies. The MC model uses specified mean lengths and volume proportions of each facies to model the transition rates. The quenching step ensures that the specified volume proportions are honored for each modeled heterogeneity realization. Second, MODFLOW cells were populated with variable conductivity drawn from a statistical distribution for each of the two stratal facies types (see Figure 2, a description below, and also [4] for a detailed description).
(28% high-K facies type by riverbed volume) with hydraulically connected high-K clusters spanning the domain, and a corresponding homogeneous case. As mentioned earlier, facies connectivity is controlled by the volume proportion of high-K facies types, and clusters become fully connected when that proportion reaches between 20 and 25 percent for a given realization. Therefore, in our non-percolating case, the volume proportion of 10% is well below threshold. Both high and low volume proportions of high-K sediments have been observed in a range of fluvial deposits [61][62][63][64]. Here, high-K facies cells are considered fully connected (i.e., percolating) if a single cluster spans between any pair of opposing domain boundaries. Our goal is to compare cases with and without percolating clusters (i.e., full connectivity). A vertical anisotropy ratio in K of 0.1 was used for all facies in all scenarios. For the high-K facies, a random number generator was used to generate ln(K) values with a mean of −2.3025 (geometric mean of 0.1 m s -1 ) and variance of 1.0 [4]. Low-K facies were assigned a mean of −7.6009 (geometric mean of 5×10 -4 m s -1 ) and a variance of 1.0. Note that these are values are horizontal K. The frequency distributions of K in the heterogeneous cases are shown in Figure 2. An equivalent homogeneous model was developed for both scenarios using the geometric mean of hydraulic conductivity (2.2×10 −3 m s −1 and 9×10 −4 m s −1 for percolating and non-percolating sediments, respectively) ( Table 1).  Our model also contains additions to the model of Zhou et al. [4] relevant to studying the influence of heterogeneity on hyporheic exchange and the associated solute mixing. For example, two scenarios were developed to quantify the influence of sediment heterogeneity on hyporheic exchange and mass transport: a non-percolating case (10% high-K facies type by riverbed volume) with limited hydraulic connectivity and no spanning cluster (see Figure 3), and a corresponding homogeneous case. These cases complement the prior cases from Zhou et al. [4]: a percolating case (28% high-K facies type by riverbed volume) with hydraulically connected high-K clusters spanning the domain, and a corresponding homogeneous case. As mentioned earlier, facies connectivity is controlled by the volume proportion of high-K facies types, and clusters become fully connected when that proportion reaches between 20 and 25 percent for a given realization. Therefore, in our non-percolating case, the volume proportion of 10% is well below threshold. Both high and low volume proportions of high-K sediments have been observed in a range of fluvial deposits [61][62][63][64]. Here, high-K facies cells are considered fully connected (i.e., percolating) if a single cluster spans between any pair of opposing domain boundaries. Our goal is to compare cases with and without percolating clusters (i.e., full connectivity). A vertical anisotropy ratio in K of 0.1 was used for all facies in all scenarios. For the high-K facies, a random number generator was used to generate ln(K) values with a mean of −2.3025 (geometric mean of 0.1 m s −1 ) and variance of 1.0 [4]. Low-K facies were assigned a mean of −7.6009 (geometric mean of 5 × 10 −4 m s −1 ) and a variance of 1.0. Note that these are values are horizontal K. The frequency distributions of K in the heterogeneous cases are shown in Figure 2. An equivalent homogeneous model was developed for both scenarios using the geometric mean of hydraulic conductivity (2.2 × 10 −3 m s −1 and 9 × 10 −4 m s −1 for percolating and non-percolating sediments, respectively) ( Table 1). Water 2020, 12, x FOR PEER REVIEW 5 of 15

Solute Transport Model
Solute transport through the hyporheic zone was simulated using the modular multispecies transport model (MT3DMS) [65], a solute transport package coupled with MODFLOW. The pore water velocities are calculated from the hydraulic head distribution within the model domain, which are then used to solve the transport equation [44,66].
We followed the common approach of using a generic solute for all solute transport simulations (e.g., [23]). The initial solute concentration in the riverbed sediments was zero. Along the sedimentwater interface, solutes were released at downwelling locations across the streambed at a concentration of 1 mg L −1 (representing river water with a constant, normalized solute concentration). Upwelling cells were defined based on interfacial flux calculations in MODFLOW and were assigned as convective flux boundaries. Solute breakthrough curves were calculated by dividing the fluxweighted concentration integrated across all upwelling cells by the flux-weighted concentration integrated across all downwelling cells. No additional source or sink terms were specified.
A linear isotherm was used to represent sorption. Thus, a retardation factor (R) was used in the advection-dispersion equation (ADE).
where is the bulk density of the subsurface medium (mg m −1 ), n is porosity (), and is the equilibrium sorption distribution coefficient (m 3 mg −1 ). The x denotes the spatial variability of in the x, y, and z directions, which results in spatial variability in R.
More robust and mechanistic models (e.g., non-linear isotherms, Surface Complexation Models) could be used to model sorption (e.g., [67,68]), but our goal was to address the lack of understanding of how simple sorption processes affect solute transport within the hyporheic zone. Such knowledge is critical in designing remediation strategies for solutes that undergo sorption (e.g., perfluoroalkyl and polyfluoroalkyl substances (PFAS), 1,4-dioxane, perchloroethylene, trichloroethylene, and trace metals) as they are transported through porous media, particularly within the heterogeneous riverbed.
To better understand the impacts of heterogeneity in the equilibrium sorption distribution coefficient, three sorption scenarios were simulated in the present study: 1) a conservative case with no sorption (i.e., no Kd), 2) a low sorption case where the of high-K facies was ten-fold lower than

Solute Transport Model
Solute transport through the hyporheic zone was simulated using the modular multispecies transport model (MT3DMS) [65], a solute transport package coupled with MODFLOW. The pore water velocities are calculated from the hydraulic head distribution within the model domain, which are then used to solve the transport equation [44,66].
We followed the common approach of using a generic solute for all solute transport simulations (e.g., [23]). The initial solute concentration in the riverbed sediments was zero. Along the sediment-water interface, solutes were released at downwelling locations across the streambed at a concentration of 1 mg L −1 (representing river water with a constant, normalized solute concentration). Upwelling cells were defined based on interfacial flux calculations in MODFLOW and were assigned as convective flux boundaries. Solute breakthrough curves were calculated by dividing the flux-weighted concentration integrated across all upwelling cells by the flux-weighted concentration integrated across all downwelling cells. No additional source or sink terms were specified.
A linear isotherm was used to represent sorption. Thus, a retardation factor (R) was used in the advection-dispersion equation (ADE).
where ρ b is the bulk density of the subsurface medium (mg m −1 ), n is porosity (), and K d is the equilibrium sorption distribution coefficient (m 3 mg −1 ). The x denotes the spatial variability of K d in the x, y, and z directions, which results in spatial variability in R.
More robust and mechanistic models (e.g., non-linear isotherms, Surface Complexation Models) could be used to model sorption (e.g., [67,68]), but our goal was to address the lack of understanding Water 2020, 12, 1547 6 of 14 of how simple sorption processes affect solute transport within the hyporheic zone. Such knowledge is critical in designing remediation strategies for solutes that undergo sorption (e.g., perfluoroalkyl and polyfluoroalkyl substances (PFAS), 1,4-dioxane, perchloroethylene, trichloroethylene, and trace metals) as they are transported through porous media, particularly within the heterogeneous riverbed.
To better understand the impacts of heterogeneity in the equilibrium sorption distribution coefficient, three sorption scenarios were simulated in the present study: (1) a conservative case with no sorption (i.e., no K d ), (2) a low sorption case where the K d of high-K facies was ten-fold lower than the K d of the surrounding low-K sediment (1 m 3 mg −1 ) (0.1 K d hereafter), and (3) a high sorption case where the K d of high-K sediment was 10 times higher than the K d of the surrounding sediment (10 K d hereafter). In sorptive scenarios, the K d values within each facies were spatially uniform. Thus, in the 0.1 K d case, there would be less sorption along preferential solute transport pathways (i.e., high-K facies, especially when they are percolating). This is reversed in the 10 K d case. For equivalent homogeneous simulations corresponding to 28% high-K facies with 10 K d and 0.1 K d , average K d values of 3.52 m 3 mg −1 and 0.748 m 3 mg −1 were assigned, respectively. Similarly, for homogeneous simulations of 10% corresponding to high-K facies with 10 K d and 0.1 K d , the average K d values were 1.9 m 3 mg −1 and 0.91 m 3 mg −1 , respectively.

Results
Heterogeneity increased hyporheic exchange by nearly an order of magnitude relative to in homogeneous scenarios. The total exchange volumes through the percolating and non-percolating heterogeneous models were 22.55 m 3 s −1 and 9.81 m 3 s −1 , respectively, while the total exchange volumes through the associated homogeneous models were 3.33 m 3 s −1 and 1.31 m 3 s −1 , respectively ( Table 1). The average interfacial flow increased with connectivity, consistent with previous studies (e.g., [4,26]). The results for both percolating and non-percolating cases are consistent with prior work on the impacts of heterogeneity on the dynamics of hyporheic exchange (e.g., [19,22,23,[69][70][71][72][73]).
Sedimentary architecture strongly influences hyporheic exchange in heterogeneous simulations by governing both the location and magnitude of flux across the riverbed. Zones of upwelling and downwelling were co-located with regions where high-K facies intersected the sediment-water interface (Figure 4a,b). In percolating simulations (i.e., 28% high-K facies by volume), nearly 70% of hyporheic exchange occurred across just 30% of the channel boundary, coinciding with locations where high-K facies intersected the channel. Conversely, channel morphology predominantly controlled the distribution of hyporheic flowpaths in homogeneous sediments, and small-scale bedform topographies were locally important (Figure 4c,d). As a result, the magnitude of flux was greater in the homogenous model corresponding to a 28% proportion of high-K facies, but the flow patterns remained consistent.
Sorption strongly influenced solute transport through the hyporheic zone, but its effect was tempered by the flow variability associated with K within the sediments ( Figure 5). Less than 10% of the infiltrated solute mass returned to the river in both the heterogeneous and homogeneous models, with the rest advected across the downstream model boundary without returning to the river due to the natural hydraulic gradient. Slightly more mass returned back to river in heterogeneous cases (both with and without percolation), but the amount was generally lower in percolating sediments, which drive exchange and direct more flow away from the river. The conservative and low-sorption (0.1 K d ) scenarios were comparable in percolating and non-percolating conditions. However, the returned mass to the upwelling cells increased by nearly 30% with increasing K d in percolating sediments, and by approximately 8% in non-percolating sediments. Roughly 7% of the infiltrated solute mass returned to the river in homogeneous scenarios. Contrary to the trend in heterogeneous sediments, the amount of sorptive solute mass that returned to the river in homogeneous sediments decreased with increasing K d .  Sorption strongly influences the breakthrough curves and the paths of solute migration. In both percolating and non-percolating heterogeneous sediments, an increase in the sorption coefficient of high-K facies altered the distribution of solutes throughout the riverbed, compressing the solute mixing zone and skewing the breakthrough behavior. When sorption was low within the high-K facies, the breakthrough behavior of sorptive solutes was comparable to that of conservative solutes; sorptive solute migration through the hyporheic zone was lower in non-percolating sediments. Sorption significantly impedes solute transport in homogeneous cases; the homogenous model corresponding to 10% high-K sediments with high sorption had the longest breakthrough curves. Breakthrough times were generally longer in non-percolating cases, and earlier breakthrough behavior was associated with percolating sediments. The sorptive behavior of high-K sediments The sorption characteristics of high-K sediments changed solute breakthrough behavior (Figure 6). High-K facies in heterogeneous sediments exhibited a higher retardation factor, which skewed breakthrough patterns ( Figure 6) and significantly delayed solute transport in both percolating and non-percolating sediments, particularly when sorption was high. More conservative mass returned to the river in the high sorption (10 K d ) scenario relative to the low sorption (0.1 K d ) scenario, suggesting the compression of the solute mixing zone (Figure 7). As both hydraulic connectivity and the sorptive capacity of high-K facies increased within heterogeneous sediments, more solute mass returned to the river. In homogeneous sediments, sorption significantly delayed solute transport and resulted in longer solute breakthrough behavior compared to conservative scenarios, with the longest breakthrough curves observed when sorption was highest. The solute mixing zone was compressed in the high sorption (10 K d ) scenario relative to the conservative and low sorption (0.1 K d ) scenarios, which were comparable (Figure 7). Thus, the solute mixing zone is compressed as sorption increases, which causes solutes only to move along short, shallow flow paths. As a result, more solute mass returns to the river from the shallow hyporheic zone, which dominates mixing and is less affected by heterogeneity structure.
Water 2020, 12, x FOR PEER REVIEW 9 of 15 caused greater solute retention and a shallower solute mixing zone, which resulted in more mass exchange through the longer flow paths.  The sorptive potential of heterogeneous sediments could influence the subsurface transport of emerging contaminants, which are becoming increasingly abundant in the biosphere. Given their novel geochemical properties, detection is often difficult, and optimal remediation strategies remain uncertain. A necessary first step in addressing this growing issue is a fundamental understanding of

Conclusions
Strongly heterogeneous sediments with high sorptive capacity compress the solute mixing zone and return more solute mass to the river than homogeneous sediments, even when the bulk physical and chemical properties are equivalent. Permeable, hydraulically connected facies form preferential flow pathways that increase hyporheic exchange and, therefore, solute delivery but shorten solute breakthrough curves along shallow flow paths. These results suggest that contaminant remediation in aquatic sediments is often least efficient where loads are greatest, an unfortunate circumstance that would presumably limit the attenuation of emerging contaminant loads in the environment. However, high-K sediments often exhibit low sorptive capacity due to having less organic matter and smaller surface sorption sites, in which case, conservative solutes may be used to estimate RT distributions and the exchange flux of sorptive solutes to the river. Given the importance of physical and chemical heterogeneity on hyporheic exchange and solute transport, their associated processes should be considered when designing remediation strategies for geogenic and emerging contaminants.

Discussion
Hyporheic exchange is strongly influenced by sediment heterogeneity, even when high-K facies types are fully connected (i.e., non-percolating case). Heterogeneity produces irregularity in pore water flow relative to homogeneous sediments, which alters hyporheic exchange patterns. By contrast, the topography of the meander, the distribution of bedforms, and sediment conductivity govern the formation of hyporheic flowpaths through homogeneous sediments. Sediment heterogeneity also strongly affects the magnitude of hyporheic exchange and solute breakthrough behavior. Heterogeneity generally increases the complexity of hyporheic exchange paths and, depending on the conductivity structure, may also moderate the depth of mixing [23]. As a result, solutes move along shorter and faster preferential flow paths, and the overall penetration depth is often reduced [22]. In this study, the breakthrough curves for heterogeneous sediments were ten-fold shorter than even the shortest breakthrough curves in equivalent homogeneous sediments. These results reflect the findings of [74], who showed that the influence of sediment heterogeneity on RTs partly depends on the bimodal variance in K. The complex facies geometries associated with bimodal sediments also produces steep pore water head gradients, which result in flatter, more vertically compressed streamlines through the heterogeneous riverbed.
Sorption strongly influences the breakthrough curves and the paths of solute migration. In both percolating and non-percolating heterogeneous sediments, an increase in the sorption coefficient of high-K facies altered the distribution of solutes throughout the riverbed, compressing the solute mixing zone and skewing the breakthrough behavior. When sorption was low within the high-K facies, the breakthrough behavior of sorptive solutes was comparable to that of conservative solutes; sorptive solute migration through the hyporheic zone was lower in non-percolating sediments. Sorption significantly impedes solute transport in homogeneous cases; the homogenous model corresponding to 10% high-K sediments with high sorption had the longest breakthrough curves. Breakthrough times were generally longer in non-percolating cases, and earlier breakthrough behavior was associated with percolating sediments. The sorptive behavior of high-K sediments caused greater solute retention and a shallower solute mixing zone, which resulted in more mass exchange through the longer flow paths.
The sorptive potential of heterogeneous sediments could influence the subsurface transport of emerging contaminants, which are becoming increasingly abundant in the biosphere. Given their novel geochemical properties, detection is often difficult, and optimal remediation strategies remain uncertain. A necessary first step in addressing this growing issue is a fundamental understanding of how sorption capacity affects reactive solute transport through the hyporheic zone, which has the geomicrobiological potential for contaminant remediation. In our models, increased sorption in high-K sediments limited solute interactions to only a small subset of the sedimentary architecture as the solute mixing zone was compressed (Figure 7), and as a result, more solute mass was returned to the river instead of being lost to groundwater. This behavior could prove advantageous in the context of emerging contaminant transport, especially given the abundance of communities that rely on groundwater as a primary water source. High-K sediments do not often exhibit high sorptive capacity, however, but instead often contain less organic matter and smaller surface sorption sites, which decrease sorptive capacity [75][76][77][78]. As a result, sorptive solutes are more strongly controlled by the riverbed heterogeneity structure when high-K facies are more abundant and hydraulically connected. Our results indicate that in such conditions, conservative solutes may be used to estimate breakthrough behavior, RT distributions, and the exchange flux of sorptive solutes to the river. Note that RT distributions are likely predictable even though they were not assessed in this paper. The fate of sorptive solutes (i.e., emerging contaminants) lost to the groundwater cannot be approximated by conservative solutes, however, due to their non-ideal transport and mixing behavior (e.g., [12,13,79,80]), and thus a more comprehensive understanding of their transport in complex groundwater systems is needed.

Conclusions
Strongly heterogeneous sediments with high sorptive capacity compress the solute mixing zone and return more solute mass to the river than homogeneous sediments, even when the bulk physical and chemical properties are equivalent. Permeable, hydraulically connected facies form preferential flow pathways that increase hyporheic exchange and, therefore, solute delivery but shorten solute breakthrough curves along shallow flow paths. These results suggest that contaminant remediation in aquatic sediments is often least efficient where loads are greatest, an unfortunate circumstance that would presumably limit the attenuation of emerging contaminant loads in the environment. However, high-K sediments often exhibit low sorptive capacity due to having less organic matter and smaller surface sorption sites, in which case, conservative solutes may be used to estimate RT distributions and the exchange flux of sorptive solutes to the river. Given the importance of physical and chemical heterogeneity on hyporheic exchange and solute transport, their associated processes should be considered when designing remediation strategies for geogenic and emerging contaminants.