How Do Fractures Influence Hyporheic Exchange in Sedimentary Rock Riverbeds?

Bedrock rivers occur where surface water flows along an exposed bedrock surface. Tinkler and Wohl (1998) define bedrock rivers as reaches where a greater proportion (>50%) of the riverbed is an exposed rock surface or is covered by a thin alluvial veneer of weathered bedrock fragments that can be easily mobilized and removed during high flows. The associated flow patterns and hydraulics of bedrock rivers are strongly influenced by their underlying bedrock geometry and inherent resistance to erosion.

natural ecosystem functions) unique to bedrock river environments. Here, we are interested in understanding how a discrete fracture network and porous rock matrix will influence hyporheic exchange pathways in a sedimentary bedrock river environment and whether these conditions will play a role in attenuating aquatic pollutants within a catchment.
To address this knowledge gap, we investigate the nature of GW-SW exchange in a sedimentary bedrock river at the bedform-scale (centimeter to meter). We do so using a field-informed numerical modeling approach to simulate bedform-scale hyporheic exchange, capitalizing on direct field measurements of hydraulic properties (in both the bedrock and river) at a previously studied bedrock river field site located along the Eramosa River, Ontario, Canada (Capes et al., 2018;Kennedy, 2017;Steelman et al., 2017). Our scale of investigation is consistent with a discrete fracture-matrix (DFM) characterization approach (Parker et al., 2012;Walton et al., 2017), which explicitly represents the fractures as distinct features with variable apertures, length, and rock matrix properties contributing to groundwater flow and solute transport (Chapman et al., 2013;MacQuarrie & Mayer, 2005). By numerically modeling flow and transport, our aim is to gain insights about the potential for hyporheic exchange in fractured bedrock rivers, as well as the physical properties and conditions controlling its spatial extent and residence time distribution.
Although the same hydrostatic and hydrodynamic drivers for hyporheic exchange within fluvial river settings will drive hyporheic exchange in bedrock rivers , differences between fracture and matrix properties within fractured rock environments are likely to result in variable spatial and temporal scales of water and solute exchange across the riverbed interface (Boano et al., 2014). Finding a way to simulate hyporheic exchange in bedrock rivers, while incorporating our current knowledge of surface flow hydraulics and groundwater flow through a DFM system, will help identify the physical conditions and mechanisms controlling the spatial and temporal scales of hyporheic exchange in these environments.
Dual porosity systems like fractured sedimentary rock can be represented as an equivalent porous medium (EPM) for flow, such that hydraulic properties are integrated using an appropriate weighting scheme. In this case, the porosity, tortuosity, and mineral surface area of a fracture-matrix system are converted to an apparent parameter set (MacQuarrie & Mayer, 2005). An EPM approach may be reasonable for bulk groundwater flow in environments with pervasive fractures that form well-interconnected networks at a macro-scale. However, it is less ideal for simulating transport because it neglects the nonlinear effects of diffusive interactions between the more transmissive fractures and less transmissive rock matrix (Chapman et al., 2014). Furthermore, EPM approaches have been shown to underestimate the natural spatial variability or heterogeneity inherent to discretely fractured porous rock systems (Cacas et al., 1990). MacQuarrie and Mayer (2005) have stated that solute travel times (especially nonreactive ones) can be significantly underestimated when using an EPM approach due to the averaging of fracture and matrix properties.
This paper aims to evaluate factors controlling hyporheic exchange pathway distributions and residence times along a fractured sedimentary bedrock river system by developing a stochastic DFM framework informed by field site data sets to capture the subsurface properties and conditions within a 2D numerical groundwater model with a boundary condition that approximates (empirically) the effect of surface water CHOW ET AL.  flow dynamics. We propose a set of motivating questions: can hyporheic exchange (i.e., closed flow cell with water entering and exiting the streambed) occur within bedrock rivers at the bedform-scale? If yes, then how does the fracture network within a low permeability porous rock matrix affect the residence time distribution of hyporheic exchange? Finally, what is the spatial extent and residence time of hyporheic exchange within a sedimentary bedrock river at the bedform-scale?
The goal of this study is to understand how fracture network heterogeneity influences the longitudinal extent of hyporheic exchange flow paths and residence time distributions commonly associated with a fracture-dominated bedrock river flow system, informed by a well-characterized field condition. In order to simulate this in a physically meaningful way, we first constrain our model using geological and hydrogeological data collected at the Eramosa Bedrock River Site (EBRS) (Capes et al., 2018;Kennedy, 2017;Steelman et al., 2017). Second, we develop a stochastic framework to simulate flow and transport for a range of discrete fracture network conditions with homogeneous and isotropic rock matrix properties.

Background and Previous Studies
A field investigation between a bedrock river, Twenty Mile Creek, and a local aquifer in Smithville, Ontario, Canada by Oxtobee and Novakowski (2002) provided some early insights into GW-SW interactions in bedrock rivers. Their investigation revealed that groundwater discharge was extremely limited within Twenty Mile Creek, with more than 95% of the groundwater underflowing the creek during baseflow conditions. Their results indicated that groundwater discharge into Twenty Mile Creek was primarily through discrete fractures and that poor vertical fracture connectivity was responsible for the limited GW-SW exchange. A subsequent numerical model by Oxtobee and Novakowski (2003) indicated that groundwater discharge (or recharge) to the bedrock river depended on the aperture size of the discharging features and the vertical hydraulic head distribution within the fracture network. Although their studies were limited to groundwater discharge to surface water, Novakowski (2002, 2003) emphasized the importance of characterizing fracture network properties when simulating interactions in fractured riverbed environments.

The Eramosa River
The Eramosa River is a major tributary of the Speed River that is located within the Grand River Watershed, Ontario, Canada ( Figure 2). This river is predominantly underlain by a bedrock surface that is variably connected to two regional bedrock aquifers: the unconfined Guelph Formation (Fm) and the confined Gasport Fm (AquaResource Inc., 2010;Brunton, 2009). These bedrock aquifers consist of densely fractured dolostone with variable occurrence of dissolution-enhanced karst features (Kunert & Coniglio, 2002;Kunert et al., 1998) resulting in high permeability and storage potential. Although these regional bedrock aquifers are known to contribute to surface water flow along the Speed and Eramosa Rivers (Cole et al., 2009), the potential effects of increased groundwater demand on GW-SW interactions along the adjacent bedrock river has only been recently explored (e.g., Capes et al., 2018;Kennedy, 2017;Steelman et al., 2017).

Geological Setting
East of Guelph, the Eramosa River incises the Vinemont Member (Mb) of the Eramosa Fm by 2-3 m, exposing a network of vertical and horizontal fractures with little to no loose sediments (i.e., fluvial deposits) along the riverbed (Kennedy, 2017). Underlying the Eramosa Fm is approximately 3 m of the cherty, marble-like Goat Island Fm (Steelman et al., 2017), which is underlain by more than 15 m of Gasport Fm (Brunton, 2009) that serves as the regional water supply aquifer. Regionally, the Eramosa Fm is up to 11 m thick and is composed predominantly of fractured, organic-rich shaley dolostone, and acts as a discontinuous aquitard between the overlying Guelph and underlying Goat Island Fms (Cole et al., 2009).

Site Description and Previous Field Investigations
In this study, we construct a numerical GW-SW model based on previous field investigations completed along the Eramosa River in eastern Guelph ( Figure 2) (Capes et al., 2018;Kennedy, 2017;Steelman et al., 2017). The EBRS has been characterized in terms of fracture frequency, orientation, and lengths CHOW ET AL.

10.1029/2020WR028476
from three vertical and three angled, continuously cored holes using core logs and acoustic televiewer geophysical logs, as well as 2D mapping of fractures exposed on the bedrock riverbed. These data support our stochastic modeling framework by constraining the variability in the fracture lengths, orientations, thus network connectivity and hydraulic properties. Further details about the EBRS fracture properties and stochastic framework are provided in the Supporting Information.
The EBRS encompasses a riffle-pool sequence approximately 61 m in length, along a channel meander of the Eramosa River with minimal loose sediment ( Figure 3). It is important to note that the words "riffle" and "pool" are borrowed from the geomorphology and river morphology literature, which depict rivers with CHOW ET AL.  The Eramosa River flows in a south-easterly direction and encompasses a riffle-pool-riffle sequence (dark blue) within a river meander (modified from Kennedy, 2017 andSteelman et al., 2017). Core holes SCA1 and SCA2 are angled beneath the river (green shading). The light blue areas within the EBRS are rubble-free zones where surficial fracture trace lengths were mapped and reported by Kennedy (2017). unconsolidated sedimentary features. For convenience, we have adopted these terms to mean riffle-like and pool-like with respect to the bedform geometry and surface water flow features due to their geometric similarities ( Figure 1). We anticipate different residence times for the different material types (i.e., fractured rock and alluvial sediment). Our study focuses on the hyporheic exchange occurring within the riverbed with an exposed fractured sedimentary rock surface ( Figure 3b). Kennedy (2017) developed a conceptual site model of the EBRS using a DFM approach (Parker et al., 2012), which uses several high-resolution measurement methods in rock cored holes. These methods include bedrock seepage meters (BSMs) specifically designed to assess GW-SW exchange in bedrock rivers. These BSMs are constructed as 0.10 m diameter, and 0.30 m deep cored holes below the riverbed with mechanically inflated rubber seals to isolate the shallow groundwater zone beneath the riverbed to measure flow along bedding plane fractures and to obtain vertical hydraulic gradients across the riverbed interface (Kennedy, 2017). A total of 24 BSMs were installed at the site to a maximum depth of 25 cm within the pool. Hydraulic heads were measured for each BSM over a 28day period beginning on July 7, 2014 and ending on August 18, 2014. The conceptual site model was further informed by a series of continuously cored boreholes, along with geophysical and hydrogeological logging to map the fracture networks beneath the riverbed, and outcrop mapping of exposed fractures along the streambed.
Although Kennedy (2017) noted local upward gradients across the pool, and thus potential for groundwater discharge, it was unclear whether there was bedform-scale hyporheic exchange (i.e., closed loops of downwelling and upwelling where surface water enters and leaves the bedrock flow system). Based on time-lapse electrical resistivity measurements within the pool and up-stream riffle, Steelman et al. (2017) noted seasonal transience in the electrical conductivity of the upper 2-3 m of the riverbed suggesting a potential zone of GW-SW interaction. However, it was unclear whether these seasonal resistivity dynamics were predominantly associated with changes in electrolytic concentration (i.e., variable groundwater/surface water mixing) or porewater temperature. Temporal analysis of thermal logs collected across the borehole network at the site by Capes et al. (2018), revealed spatial and temporal variations in groundwater temperature associated with flow in discrete fractures. Although these thermal measurements delineated active groundwater flow zones deeper in the bedrock, the presence of surface casing within the upper 2 m of rock limited their utility in assessing fluid dynamics near the riverbed interface.

Site Conceptual Model and Modeling Approach
A major issue with 3D simulations of flow and transport, that incorporate the details of discrete fracture networks, is their high computational cost. As a possible remedy, Chapman et al. (2014) proposed a combined equivalent porous media-discrete fracture network (EPM-DFN) approach to simulate transport in sedimentary bedrock aquifers. The EPM-DFN approach involves using a 3D EPM model that incorporates macro-complexity (e.g., hydrogeological units) to capture bulk groundwater flow, then to identify the predominant plume flow path using forward particle tracking from the source zone. Once the plume flow path is identified, a 2D vertical cross-section with discrete fractures is used to model solute transport.
To determine the GW-SW exchange potential in a fractured river system, we model flow and transport within a 2D DFM cross section that is consistent with the anticipated scale of hydraulic interactions (  groundwater flow direction is assumed to be the same as the river flow direction. Therefore, we only use the latter half of the EPM-DFN approach (i.e., the 2D vertical cross section) proposed by Chapman et al. (2014). The 2D DFM steady-state flow and transient transport model was developed using HydroGeoSphere™ (Aquanty Inc., 2015).
Our conceptual model of GW-SW flow at the EBRS is shown in Figure 4. The model comprises of a 2D cross section parallel to the thalweg along a 61 m reach of a riffle-pool sequence. Our conceptualization addresses potential bedform-scale hyporheic exchange between the Vinemont Mb and the surface water of the Eramosa River. The top of the model domain follows the elevation of the riverbed along the thalweg based on riverbed elevation measurements (Figure 5a; Kennedy, 2017). The model extends downward to 303.0 meters above sea level (masl), which is the average elevation between the Vinemount Mb and the underlying Goat Island Fm.

Model Boundary Conditions
Hyporheic exchange at the bedform-scale is predominantly driven by hydrodynamic pressure variations for river reaches with a relatively gentle slope (Boano et al., 2014). We used a simplified approximation to represent hyporheic exchange across a riffle-pool sequence (e.g., Cardenas, 2008;Käser et al., 2013;Pryshlak et al., 2015;Stonedahl et al., 2010), whereby the total hydrodynamic pressure head h total along the riverbed is treated as a constant head boundary condition approximated by the modified representation of the basic pumping model for a single bedform, which in our case is a riffle-pool sequence (Elliott & Brooks, 1997;Fehlman, 1985;Stonedahl et al., 2013): where s h is the hydrostatic head component and  m h is the half-amplitude of the head variation given by:   Steelman et al., 2017). Hyporheic exchange is controlled by the open fractures, which are more transmissive than the surrounding rock matrix. Open fractures will be hydraulically connected to the porous matrix, and thus, support potential solute exchange (i.e., by diffusion) between the two flow regimes.
shows the resulting boundary condition. The hydrodynamic head low within the first pool aligns with the maximum upward vertical head gradient observed in the field (Kennedy, 2017). This boundary condition will drive upwelling within the first pool, which is consistent with field observations of fluid potential.
The hydrostatic head component (h s ) was set to a constant value of 311.022 m, which represents the average river stage elevation. A constant stage elevation is deemed reasonable given the relatively short and gentle slope of the EBRS (Figure 3). Thus, by assuming the water level fluctuations over this reach are negligible, we are able to examine the impact of hydrodynamic head variations on hyporheic exchange at the bedform-scale, independent of larger-scale influences from groundwater flow.
Steady-state hydraulic head conditions are used to represent long-term conditions at the EBRS, negating the need for specific storage parameters.
Hydraulic packer testing at the EBRS indicated that the Goat Island Fm had a bulk average    10 9.4 10 m / s K (Kennedy, 2017). Therefore, we treat the Goat Island Fm as an aquitard and assign a no-flow boundary to the bottom of the model domain. We assigned the left and right sides of the model domain with a constant head boundary of 311.022 m, which represents the average river stage elevation over the entire length of the EBRS based on an average flow rate of 1.5 m 3 /s (Kennedy, 2017). Further details on the model discretization can be found in the Supporting Information.
CHOW ET AL.   (Kennedy, 2017). The hydrodynamic head component is calculated using the basic pumping model for a single bed form (Elliott & Brooks, 1997;Fehlman, 1985) and added to a hydrostatic head value of 311.022 masl, which represents the mean river stage elevation.

Stochastic Framework for Discrete Fracture Generation
Our DFM model is evaluated within a stochastic framework, which means the subsurface is represented by an ensemble of 100 DFM realizations. Each DFM realization is randomly generated based on statistical distributions consistent with observed spatial and depth-dependent (i.e., geomechanical layering) fracture properties obtained from the field site, which is then used to simulate a steady-state head distribution and the transient transport of a conservative solute tracer within the bedrock riverbed. Probability distribution functions of hyporheic exchange residence times and spatial extents are then produced from the ensemble of solute transport simulations. This results in ranges of bedform-scale hyporheic exchange fluxes and residence times as opposed to a single deterministic result. The field data collected at the EBRS and their subsequent reduction to input parameters and boundary conditions for our 2D DFM stochastic model is described in the Supporting Information.

Transport Setup
To assess the hyporheic exchange potential, a conservative tracer was released for each DFM realization and tracked for 100 years. This allowed us to determine the residence time and spatial extent of various hyporheic exchange pathways. The conservative tracer was applied for each DFM realization over a 1-day period to the left-most surficial vertical fracture located between 15.0 and 19.4 m of the model domain ( Figure 6). Here, 15.0 m (approximately 25% of the longitudinal domain length) was chosen as the left-most extent for the tracer placement to avoid potential boundary effects, while 19.4 m was chosen as the right-most extent because it is where the elevation begins to decrease (i.e., where the pool begins). Releasing the conservative tracer in the left-most vertical fracture between 15.0 and 19.4 m meant that we were injecting along the largest gradient between the riffle (downwelling) and the pool (upwelling), and thus, simulating a dominant exchange pathway.
To assess the residence time distribution of the instantaneous 1-day tracer injection across the riffle-pool hyporheic exchange pathway, the total mass balance of the conservative tracer within the model domain was analyzed with respect to time for each DFM realization. The total tracer mass within the model domain increases over the injection time t 0 , until reaching a maximum mass M max (Figure 7a). Afterward, the mass within the model domain will decrease as the conservative tracer leaves the model domain through the top or side boundaries of the model (Figure 7b). Solute mass exiting the top of the model is conceptualized as bedform-scale or riffle-pool-scale hyporheic exchange, while the mass leaving through the sides of the model is conceptualized as larger-scale (e.g., potential meander-scale) hyporheic exchange beyond the extent of our model domain. To assess the solute exchange capacity between the high-velocity fracture network and the lower-permeable matrix, a cumulative distribution function (CDF) M norm (t) of the residence times CHOW ET AL.   Figure 6 shows the steady-state head distribution from a single DFM realization that closely resembles the average (in terms of residence time distribution) of the model ensemble. The area of highest vertical hydraulic head within the riverbed resides within the first riffle. Areas of lowest vertical hydraulic head reside within the first and second pool. This indicates that downwelling will occur in the first riffle, while upwelling will occur in the first and second pool. Figure 8 shows the transport of conservative tracer viewed over four orders of magnitude concentrations at select time steps (1 year, 25 years, 50 years, and 100 years). A constant initial concentration c 0 = 1,000 mg/L (1 kg/m 3 ) was applied over a single 1-day time-step (t 0 ) within the upper-most node of a vertical fracture that intersected the top surface of the model. Conceptually, this represents a 1-day pulse of conservative tracer injected within a vertical fracture that is hydraulically connected to the river. Figure 8a represents a DFM scenario based on the field data collected at the EBRS which includes a matrix porosity of 4%; Figure 8b depicts the same fracture scenario with a near-zero matrix porosity of 0.4%. A matrix porosity of 0.4% represents the upper limit for crystalline rock (Skagius & Neretnieks, 1986), which is generally considered to have negligible matrix porosity for decadal to century scale transport simulations. This low porosity scenario illustrates the impact of matrix diffusion on the distribution of solute mass within our sedimentary fractured riverbed.

Transport of a Conservative Tracer
Although much of the mass remains within the footprint of the first riffle for both scenarios, some of the mass reaches the first pool and discharges along a series of vertical fractures between 22 and 28 m. However, the solute extends slightly farther downgradient in the lower matrix porosity case (Figure 8b). At later times (i.e., 25, 50, and 100 years), connected fracture flow paths at depth (between 303.5 and 305.5 masl) deliver small concentrations (∼1 μg/L) of tracer mass to a zone of lower hydraulic head within the reach of the second pool, between 41 and 45 m. Although higher matrix porosity appears to slow the migration of solute within the riverbed the spatial extent of solute mass in the riverbed are similar in scale.

Hyporheic Exchange Residence Times
Conservative tracer mass was tracked until 95% of M max exited the domain. Most of the transport simulations meet this threshold without reaching an excessively long simulation time (>200 years). Therefore, all realizations for the residence time analysis have 5% of tracer mass remaining in the domain. Figure 9a shows the resulting ensemble of bedform-scale hyporheic exchange residence time CDFs. The DFM realization (no. 39) that produced the residence time CDF closest to the median (thick blue line in Figure 9a    was chosen as the most representative DFM realization for our fractured rock riverbed (shown in Figures 6  and 8). These results indicate that the residence times for hyporheic exchange of a conservative tracer in a fractured sedimentary bedrock riverbed with matrix porosity n = 0.04 range from hours to decades. The average residence time (i.e., 50% CDF) for bedform-scale hyporheic exchange in our sedimentary fractured bedrock river is 13 years. The interquartile range for the average residence time is between 4 and 35 years (Figure 9a), which indicates that the potential variability in residence time is large due to the inherent variability in the fracture network connectivity and properties (e.g., fracture spacing, length, orientation, and aperture). Once plotted on a logarithmic time scale (Figure 9b), the median residence time CDF (no. 39 with n = 0.04) reveals an inflection point indicating two distinct pulses of solute exiting the riverbed. Approximately 50% of the mass that enters the riverbed appears to discharge within 5 years and 5 m downstream from the injection site. The remaining mass discharges at a slower rate within the first pool between 20 and 30 m downstream from the injection site, taking on average 100 years before reaching the 95% threshold.
The effect of matrix porosity (n) on solute residence time is shown by the case of lowering the porosity from 0.04 to 0.004 for a single realization (red line in Figure 9b), which effectively depicts a fractured and weathered crystalline rock. Lowering the matrix porosity resulted in considerably lower residence times due to a reduction in mass transport rates into the matrix by reducing the effective diffusion coefficient and cross-sectional area for diffusion. The reduced matrix porosity also lowers the volume for temporary solute storage available in the immobile domain, which leaves more solute tracer mass within the higher velocity fractures. Lowering the matrix porosity enables solute mass to travel more quickly through the fractures without uptake and release from the matrix, which ultimately increases the mass discharge rate within a hyporheic flow cell. Nevertheless, a portion of the mass back diffuses into the fractures over an extended time period. By lowering matrix porosity, the average residence time for the solute was reduced from 13 to 1.5 years and the tailing (in terms of the 95% mass discharge) was reduced from 112 to 12 years.
To further illustrate the impact of fractures on hyporheic flow, we performed a single simulation with a homogeneous isotropic porosity of 0.3 and K = 1 × 10 −3 m/s to represent a sandy alluvial riverbed (green line in Figure 9b). This scenario does not include the contribution of discrete fractures within a porous rock matrix, and thus does not include interactions between matrix diffusion and fracture transport. The absence of these interactions dramatically reduces the residence time, from an average of 13-0.015 years (∼5.4 days), and the tailing (in terms of the 95% mass discharge) from 112 to 0.25 years (∼90 days). The reduced solute residence time for the homogeneous porous media scenario is the result of a higher bulk K-value and the absence of dual porosity interactions. Figure S7 in the Supporting Information depicts the transport of a conservative tracer at select time steps for the homogeneous porous media case. This porous media simulation shows that >85% of the tracer mass will exit the riverbed within the first 0.03 years (∼11 days) after injection approximately 3 m downstream from the injection site, followed by the remaining mass that exits the riverbed from the first pool approximately 25 m downstream after 0.1 years (∼36 days) as it approaches the 95% threshold. The marked difference in the timing and magnitude of the initial pulse of mass exiting the riverbed relative to the fractured rock scenarios indicates that discrete fractures can have a dominant influence on hyporheic exchange pathways.

Spatial Extent of Hyporheic Exchange in a Fractured Rock Riverbed
To evaluate the spatial extent of hyporheic exchange, we analyze the spatial distribution of normalized mass leaving the domain (Figure 10). Over the 100 DFM realizations, the average total mass that entered the model domain after the 1-day pulse injection of 1 kg/m 3 was M max = 1.1 × 10 −4 kg. Thus, 0.0001% is considered an estimate of the fraction of river flow that downwelled into the individual fracture. The blue line in Figure 10 shows the normalized mass exiting the riverbed (i.e., model domain top surface) averaged over 100 DFM realizations. The area under the blue line sums up to 93% of the injected mass, showing that most of the tracer eventually exited through the riverbed. Of the remaining mass, only 2% exited through the right side of the model domain and 5% remained within the model domain after 100 years. Most of the tracer mass exits the riverbed within the vicinity of the injection interval indicating that the predominant hyporheic exchange pathway occurs within a ∼4 m reach of the river. The next most common upwelling zone occurs within the first pool between 20 and 30 m, which coincides with the total hydraulic head low and the maximum observed upward vertical hydraulic gradient. These proximal and distal discharge zones are consistent with the two pulses of mass discharge highlighted by the residence time CDF (Figure 9b). Only minute quantities (<0.1%) of tracer mass exits the riverbed within the second pool between 41 and 45 m. Our analysis of the spatial distribution of the standard deviation among the realizations ( Figure S6), indicates that the greatest variability in mass discharge will occur where the majority of the mass leaves the domain and decreases farther down gradient away from the injection point.

Discussion
Previous studies at the EBRS have indicated the potential for GW-SW exchange across the fractured bedrock river based on fluid potential measurements (Kennedy, 2017), seasonal temperature dynamics within the upper few meters of riverbed observed using surface geophysics (Steelman et al., 2017), and along discrete fractures intersecting boreholes underlying the river (Capes et al., 2018). To further our understanding of potential hyporheic exchange in a sedimentary bedrock river, we produced a probabilistic distribution of potential hyporheic exchange extents and residence time distributions based on parameters measured from the field site and a physically based numerical GW-SW model.
Our results indicate that the spatial extent of bedform-scale hyporheic exchange will be largely controlled by the coincidence of fractures connected to the river surface and the hydrodynamic head gradients. Changes in matrix porosity (e.g., 4% in a sedimentary rock vs. 0.4% for a crystalline-type [i.e., granite] rock) will have a marked impact on the average residence time for bedform-scale hyporheic exchange. An average residence time of 13 years can be expected at our sedimentary bedrock river field site, while much shorter times (1.5 years) may be encountered at an equivalent fractured crystalline bedrock riverbed. These fracture rock scenarios represent considerably longer residence times compared to those reported for rivers composed of unconsolidated sediments (Boano et al., 2014).
Our stochastic simulations indicate that the average residence time may not be a good proxy for the potential natural attenuation capacity of fractured sedimentary bedrock rivers due to more skewed distributions in mass discharge along the riverbed and longer hyporheic exchange residence times compared to rivers composed of unconsolidated sediments. In addition, matrix diffusion becomes the dominant process in fractured sedimentary bedrock environments, leading to significant retardation and attenuation of surface water solutes delivered into a bedrock riverbed through hyporheic exchange along the fracture networks.
Given the delivery of solute into a fractured porous rock, we expect that interconnected fractures and dissolution features would be the only oxygen-rich portions of the riverbed, while the rock matrix porewater between the fractures would be more reducing. This could have significant implications for the spatial CHOW ET AL. distribution of biogeochemical processes, hence mobility of redox sensitive solutes, in a fractured bedrock river. In situations with high matrix storage without reactions, eventual back-diffusion of contaminants from the rock matrix into the fracture porosity would be a source of long-term trace contaminant concentrations, leading to persistent contamination problems that would be difficult to remediate. In contrast, a crystalline bedrock river with much lower diffusion rates due to much smaller matrix porosities would not possess these strong diffusion influences, and thus would have shorter solute residence times.
Nonetheless, a direct comparison of the current stochastic DFM model with a porous media scenario might provide additional insight into effective dispersion, sorption, and anisotropy values. Approximating fractured rock systems with an appropriately informed EPM model requires inclusion of nonlinear effects with time and distance scaled to the advection-diffusion interaction specific to the geologic conditions. If properly informed, these EPM simulations could be useful for more regional-scale water resource studies, bridging the gap that currently exists with computationally demanding DFM models. Additional analysis for future research could be a more in-depth sensitivity analysis of the various stochastic DFM model parameters to identify which fracture property dominates the simulated hyporheic exchange response (i.e., residence time distribution in space and time, tracer extents, peak flux magnitude and position, and concentration attenuation timelines for various solute types). This could include stochastifying parameters that were assumed to be constant (e.g., mechanical layer thickness, dip angle, longitudinal and transverse dispersivities, matrix porosity, tortuosity, diffusion coefficients).

Conclusions
This study evaluated the occurrence and spatial extent of hyporheic exchange within a sedimentary bedrock riverbed by quantifying the advection-diffusion interaction using a 2D discrete-fracture matrix numerical model based on parameters and conditions from a field research site. This was achieved by developing a 2D numerical groundwater model for flow and transport using boundary conditions that approximate (empirically) the effect of observed surface water flow dynamics, and the elevation of the bedrock riverbed along the thalweg with fracture and matrix parameters informed by field measurements. Flow was simulated within a stochastic framework of orthogonal discrete fractures in a porous sedimentary rock matrix, followed by transport simulations to quantify residence times of a conservative solute tracer released at a fracture intersecting the riverbed. Through this methodology we were able to evaluate the potential hyporheic exchange behavior along a well-characterized reach of a sedimentary bedrock river and the influence of matrix diffusion by comparison to a low matrix porosity case. Our study addressed three research questions: 1. Can hyporheic exchange occur within sedimentary bedrock rivers at the bedform-scale? Our model results show that indeed closed loops of hyporheic exchange are possible at the bedform-scale within bedrock rivers, and that the coincidence of vertically and horizontally connected fractures to the riverbed surface and hydrodynamic head gradients will largely control the longitudinal extents of hyporheic exchange pathways within a fractured bedrock river setting. 2. How does the fracture network within a porous, low permeability rock matrix affect the residence time distributions of hyporheic exchange? Bedform-scale hyporheic exchange in a sedimentary bedrock river will be orders of magnitude slower compared to a riverbed composed of unconsolidated sediments (i.e., alluvium) or low matrix porosity rocks such as granites. A fractured porous sedimentary rock will exhibit a significant delay in the residence time of dissolved groundwater solutes due to diffusion into the low-permeable rock matrix, effectively creating a reservoir of mass for delayed exchange back to the more permeable fracture network. This produces a long tailing with respect to the solute residence time and possibly accentuating the influences of hyporheic exchange processes in sedimentary bedrock river systems. As a result, the average residence time may not be a good indicator of hyporheic exchange given that a significant portion of solute will ultimately back diffuse from the rock matrix over an extended period. This delay in solute travel times within the bedform will be much less pronounced in rock with low matrix porosities (i.e., crystalline rocks) since there would be minimal, if any, matrix diffusion; however, abrupt fracture terminations and limited interconnectivity can contribute to zones with lower flow rates and longer residence times for some zones within the riverbed. Future research could investigate alternative metrics to describe hyporheic exchange potential in sedimentary bedrock rivers based on the timing of peak fluxes.
3. What is the extent and residence time of hyporheic exchange in a fractured sedimentary bedrock river?
Our results indicate that most of the conservative tracer mass would be recovered within a 15 m reach of a riffle-pool sequence. The average residence time from 100 DFM realizations is approximately 13 years with a considerably large uncertainty (i.e., tens of years) mainly due to inherent variability in fracture properties (e.g., fracture spacing, length, orientation, and aperture). While residence time in a fractured rock riverbed can drop by an order of magnitude with a similarly scaled reduction in matrix porosity, it is the fractures that enable deeper penetration of solutes into the riverbed, thus increasing the residence time relative to an unconsolidated fluvial system. Discrete fracture networks will influence the timing and magnitude of mass flux to a river and may contribute to a more variable mass discharge over time.
These field-based numerical results represent an important step toward the design of a conservative tracer (i.e., bromide) experiment to directly monitor hyporheic exchange in a fractured sedimentary bedrock river that could also be controlled by local groundwater flow gradients. Although we focus on hydrodynamic pressure variations at a local-scale, our numerical results indicate that bedform-scale hyporheic exchange pathways will largely be on the order of a few to several meters in extent along the fastest flowing sections of a river (i.e., the thalweg) and that a conservative tracer experiment would require a multi-year monitoring strategy.

Data Availability Statement
Original data sets are available through Kennedy (2017), Steelman et al. (2017), and Capes et al. (2018). Reduced data sets used in this research are provided in Tables S1-S7 of the Supporting Information.