A multirate mass transfer model to represent the interaction of multicomponent biogeochemical processes between surface water and hyporheic zones (SWAT-MRMT-R 1.0)

Surface water quality along river corridors can be modulated by hyporheic zones (HZs) that are ubiquitous and biogeochemically active. Watershed management practices often ignore the potentially important role of HZs as a natural reactor. To investigate the effect of hydrological exchange and biogeochemical processes on the fate of nutrients in surface water and HZs, a novel model, SWAT-MRMT-R, was developed coupling the Soil and Water Assessment Tool (SWAT) watershed model 5 and the reaction module from a flow and reactive transport code (PFLOTRAN). SWAT-MRMT-R simulates concurrent nonlinear multicomponent biogeochemical reactions in both the channel water and its surrounding HZs, connecting the channel water and HZs through hyporheic exchanges using multirate mass transfer (MRMT) representation. Within the model, HZs are conceptualized as transient storage zones with distinguished exchange rates and residence times. The biogeochemical processes within HZs are different from those in the channel water. Hyporheic exchanges are modeled as multiple first order mass 10 transfers between the channel water and HZs. As a numerical example, SWAT-MRMT-R is applied to the Hanford Reach of the Columbia River, a large river in the United States, focusing on nitrate dynamics in the channel water. Major nitrate contaminants entering the Hanford Reach include those from the legacy waste, irrigation return flows (irrigation water that is not consumed by crops and runs off as point sources to the stream), and groundwater seepage resulted from irrigated agriculture. A two-step reaction sequence for denitrification and an aerobic respiration reaction is assumed to represent the biogeochem15 ical transformations taking place within the HZs. The spatially variable hyporheic exchange rates and residence times in this example are estimated with the basin-scale Networks with Exchange and Subsurface Storage (NEXSS) model. Our simulation results show that 1) given a residence time distribution, how to approximate the exchange fluxes to HZs when using MRMT can significantly change the amount of nitrate consumption in HZs through denitrification; and 2) source locations of nitrate have different impact on surface water quality due to the spatially variable hyporheic exchanges. 20

Hyporheic exchange and biogeochemical processes within HZs are largely absent in reach-scale models (Helton et al., 2010). We need a modeling framework that integrates the physical and biogeochemical processes in the surface water and HZs in order to bridge the gap in reach-scale models and answer questions of when and where HZs become important. An 35 integrated hydrologic model that couples land surface processes with spatially-explicit three-dimensional groundwater flow processes, such as SWAT-MODFLOW (Bailey et al., 2016), can explicitly simulate hyporheic processes. However, it is often computationally demanding to capture the hyporheic zones that are highly variable in space and time, and requires more model parameters that are difficult to characterize (Painter, 2018). In reality, hydrologic exchange between surface water and groundwater was commonly modeled by the transient storage model to investigate the physical processes governing riverine 40 solute transport (Bencala and Walters, 1983;Briggs et al., 2009;Boano et al., 2014;Runkel et al., 2003;Haggerty and Gorelick, 1995;Zaramella et al., 2003) because it is simple, accessible, and adaptable (Knapp and Kelleher, 2020). The transient storage model is a lumped model representing first-order exchange between the main channel and well-mixed storage zones or dead zones, usually parameterized with conservative tracer test. The transport of solute between the storage zone and channel is simply determined by the solute concentration difference between the channel and storage zone and an exchange coefficient 45 (Bencala and Walters, 1983). It has been modified to include various zones with one exchange coefficient within each zone (Briggs et al., 2010;Neilson et al., 2010), and to simulate reactive solute within the storage zones by coupling with a chemical equilibrium model (Runkel et al., 1996); however, the transient storage model is not capable of representing the long solute tail observed in the stream (Cardenas, 2015;Painter, 2018). Coupled with a transient storage zone solute transport model, Ye et al. (2012) used a dynamic hydrologic network model to simulate dissolved nutrient retention processes during transient flow 50 events at the channel network scale. They assumed a constant exchange rate and simple uptake of nutrient in the storage zone, but the interplay between mass exchange and biogeochemical processes as well as hydrological and geological control on mass exchange rates were ignored.
In this study, to model the exchange between the channel and HZs, we extended the transient storage model by using the multirate mass transfer (MRMT) formulation incorporating a spectrum of transition times that is commonly used to simulate 55 the late-time solute and contaminant tails observed in porous media (e.g., Haggerty and Gorelick, 1995;Wang et al., 2005;Liu et al., 2008;Fernandez-Garcia and Sanchez-Vila, 2015). The MRMT formulation is a discrete equivalent to convolution-based representations to model mass exchange, which can be solved in Lagrangian domain (Silva et al., 2009). Painter (2018) recently generalized the convolution-based model to include multicomponent reactive transport with general nonlinear reactions and tested the model using a single reach, assuming a steady state hyporheic flow field. It has not yet been integrated in watershed 60 modeling. Here we generalized the MRMT model to include mass transfer between the stream channel and multiple storage zones according to a spectrum of rates within each storage zone. We developed an integrated model, referred to as SWAT-MRMT-R, based on two open-source codes to simulate in-stream biogeochemical processes, mass transfer between the stream and HZ, and biogeochemical processes in HZs. These two codes are the Soil and Water Assessment Tool (SWAT) watershed model  and the subsurface flow and reactive transport code, PFLOTRAN (Lichtner et al., 2017). We applied the integrated model to the Hanford reach of the Columbia river, a large river in the southern Washington State of the United States. At the simulated reach, major impacts to stream nutrients come from the contaminated groundwater, irrigation return flows, and groundwater seepage resulted from irrigated agriculture (Evans et al., 2000). This paper starts with an overview of the components in the SWAT-MRMT-R modeling framework, followed by description of each component (Sect. 2). The model is then applied at the Hanford reach using the spatially variable hyporheic exchange 70 rates and residence times estimated with the basin-scale Networks with Exchange and Subsurface Storage (NEXSS) model (Gomez-Velez and Harvey, 2014) (Sect. 3). We focus on the nitrate dynamics in the results (Sect. 4), varying modeled processes and parameters. Finally, model limitation and future work are discussed. Our modeling framework is built within the in-stream nutrient submodel in SWAT, considering distinguished biogeochemical processing in both the water column and hyporheic zones. The SWAT model is developed by the U.S. Department of Agriculture to predict the impact of land management on water, nutrients, and sediments of large ungauged basins  and it is widely used worldwide. A schematic description of the modeling framework is shown in Fig. 1. At the top of Fig. 1A, the processes pointed to by the blue arrow are reactions related to the water column. The processes in the orange 80 circle at the bottom of the figure are example reactions related to the hyporheic zones (shown as black rectangular prisms with orange streamlines) in the vertical and lateral directions. The orange streamlines represent the river water entering the vertical and lateral hyporheic zones and returning back to the river within the same reach over a distribution of travel times, i.e., each streamline carries a residence time (τ s ) and exchange flux (q s ). In this study, the residence time and exchange flux are estimated by NEXSS (Gomez-Velez and Harvey, 2014), which will be briefly described in Sect. 3.3. We use an MRMT 85 model to represent solute mass transfer dynamics between the channel and storage zones, capturing the effect of hyporheic exchange in the fate and transport of solutes along the river corridor. In this model, all of the reactions and exchange processes shown in Fig. 1A are simply conceptualized as a series of hyporheic batch reactors connected to a stream batch reactor through mass exchange as shown in Fig. 1B. The hyporheic zones do not communicate with each other, but only with the stream water.
The mass exchange between the hyporheic reactor and the stream reactor is provided by NEXSS through user input. In other 90 words, we only use the NEXSS output and it is not part of the model development. The batch reactors in Fig. 1B are simulated simultaneously with the implicit time stepping through the Newton Raphson method in batch mode (i.e., no transport) of PFLOTRAN (Lichtner et al. (2017). In this study, we only modified the in-stream nutrient transformation module in SWAT and did not modify the hydrologic flow and solute transport by taking advantage of the operator-splitting approach in which the solute transport and the reaction steps are solved separately in SWAT. The model in this study is not limited by the number 95 of storage zones. The user can provide the number of zones through the input file without having to modify the code. In the following, we describe each of these components and their coupling in detail.

Mass Exchange Processes between the Channel and Storage Zones
MRMT is commonly used to simulate the late-time solute and contaminant tails in porous media (e.g., Haggerty and Gorelick, 1995;Wang et al., 2005;Liu et al., 2008;Fernandez-Garcia and Sanchez-Vila, 2015). We use it here to simulate mass transfer 100 between the stream channel and multiple storage zones according to a spectrum of rates within each storage zone. The mass exchange rate of a solute component due to hyporheic exchange is:

Mass Balance Equations of the Coupled System
Considering the biogeochemical processes in the water column, storage zones, and hyporheic exchange, the final mass balance equations for dissolved components in the channel and storage zones are: The area of the storage zones are calculated as a function of each storage zone's water exchange flux (q s ) and storage residence time (τ s ):

120
Equations (1)-(4) are a parsimonious version of the general form of the MRMT model (e.g., Haggerty and Gorelick, 1995;Anderson and Phanikumar, 2011) where a discrete number of storage zones, each with different geometry, linear rate, and biogeochemical reactions, are considered.

Water Column Biogeochemical Processes
The in-stream submodel within SWAT is based on the QUAL2E (Brown and Barnwell, 1987) model. QUAL2E simulates major 125 interactions of the nutrient cycle, atmospheric aeration, algae production, benthic oxygen demand, and carbonaceous oxygen contribute to r m in Eq. (2).

Storage Zone Biogeochemical Processes
Although any type of biogeochemical process of interest can be simulated by the model, we focus on the following reactions for aerobic respiration and denitrifiction in each storage zone:

140
where the reactions R1, R2, and R3 are associated with the electron acceptors O 2 , N O − 3 , and N O − 2 , respectively, and f 1 , f 2 , and f 3 are the fractions of the electron equivalents in CH 2 O used for energy production in each reaction. We used the approach proposed by Song et al. (2017) and Song et al. (2018) to model the reaction rates for R1, R2, and R3 (denoted as r 1 , r 2 , and r 3 , which contribute to r ms in Eq. (3):

145
These reaction rates incorporate unregulated and regulated effects. In particular, the unregulated effect is represented by a Monod-type kinetics coefficient (r kin i ) of the form where the k i , K a,i , and K di denote the reaction rate [mol mol −1 d On the other hand, the regulated effect is dictated by the enzyme activity (e i ) and determined by the cybernetic control law (Young and Ramkrishna, 2007) as follows (Song et al., 2018): This approach associates a community's traits with functional enzymes without directly considering the dynamics of indi-155 vidual microbial species or their guilds that synthesize their enzymes (i.e., the entire microbial community is treated as a single organism) and allowed the estimation of reaction stoichiometries and rates for denitrification, as well as biomass degradation rate (Song et al., 2017). Rather than using empirical inhibition kinetics, the formulation in Eq. 8 enables the oxygen regulation of denitrification in response to environmental variation (Song et al., 2018). During each time step, the initial condition in the water column was calculated using the following equation : , and V f lowin is the volume of water entering the channel at the end of time step [L 3 ], calculated using the Muskingum routing method (Overton, 1966

SWAT Model Configuration
We set up the SWAT model by combining a series of geospatial data. We derived topography information from the U.S.
Land covers including shrubland, forestland, grassland, developed land, barren land, and cultivated land were derived from

Mass Exchange Parameters
Quantifying hyporheic exchange in streams is difficult. Exchange rate coefficients are often determined by tracer injection in the stream (e.g., Harvey et al., 2013;Gooseff et al., 2011;Laenen and Bencala, 2001). Assuming steady uniform flow, the river corridor hydrologic exchange model NEXSS in Gomez-Velez and Harvey (2014) can estimate hydrologic exchange flux rates and distribution of residence times that can be used to parameterize the exchange rate coefficients. NEXSS is a single 205 tool that consolidates multiple previously published geomorphic and hyporheic flow models for the analysis of exchange at the basin scale (Gomez-Velez et al., 2015). The river network is discretized into individual reaches of length L, where values of channel width, bankfull width, depth, bankfull depth, discharge, average channel velocity, median grain size, channel slope, sinuosity, and regional hydraulic head gradient in the x and y direction are prescribed. The model is applied in two steps: geomorphic characterization and hyporheic exchange modelling. For each channel reach, the hyporheic exchange modelling 210 predicts average hyporheic exchange flux per unit area of streambed and residence time distributions, as well as median residence time for vertical exchange beneath submerged bedforms (ripples, dunes and alternate bars) and lateral exchange outside the wetted channel through emergent alternate bars and meanders. To estimate the net amount of vertical exchange flux, vertical exchange is conceptualized as a three-dimensional process, assuming a homogeneous and isotropic porous medium. The boundary conditions of the system are: the top boundary at the sediment-water interface has a prescribed head distribution, 215 lateral boundaries (y = 0 or w bf (bankful channel width) are impervious (q y (x, y = 0, z) = q y (x, y = w bf , z) = 0), and the bottom boundary has a prescribed flux, which is defined as q u ≈ 0.57KJ y (Boano et al., 2009), where K is hydraulic conductivity, and J y is the mean head gradient across the alluvial valley. Basal flux through the dunes and ripples induced by the regional groundwater flow is defined as q b = KJ x (Elliott and Brooks, 1997), where J x is the mean head gradient along the valley in the downstream direction. The lateral exchange is conceptualized as steady and two-dimensional process, and the flow NEXSS reach is different from the discretization used in the stream network delineated for the SWAT simulation (Fig. 2b). To address this issue of mapping the NEXSS exchange metrics into the SWAT stream network, we implemented a procedure that matches channels and transfers parameters over a rectangular grid that covers the entire watershed. As the reaches from the two networks do not exactly overlap (see Fig. 2b), we used a grid size of 200 m to search for matching reaches. Reaches that 240 intersected with this grid and the length of each reach within each grid cell is calculated. For example, NEXSS reaches 386 and 380 in (Fig. 2b)  Mapping the NEXSS output to the simulated watershed, we found that the reaches along the Columbia River are characterized by vertical residence times ranging from 8 hours to 10 days, and the lateral residence times ranging from 17 days to 10 years (Fig. 3). The vertical and lateral residence times for the reaches highlighted in Fig. 2a are shown in Table 1 in the Supplement. Compared to other reaches in the watershed, the Columbia River is characterized by relatively larger exchange flow and shorter residence times in vertical storage zones (Fig. 3). Hydraulic conductivity, head gradient, stream geomorphol-255 ogy, and flow paths together determine the hyporheic residence times (Kasahara and Wondzell, 2003;Cardenas et al., 2008).
The vertical exchange residence times are short, probably due to their hydraulic and geomorphic conditions that favor ripple and dune formation (Gomez-Velez et al., 2015). For example, hyporheic flow around steps and riffles with coarse-textured alluvium can result in short residence time distributions for hyporheic exchange flows (Kasahara and Wondzell, 2003). Vertical exchange flow is dominant along the Columbia River.

270
Starting from the base case, where mass transfer and reactions in the HZ were not simulated, we added mass transfers using two storage zones to represent vertical and lateral exchanges, and within each zone only one exchange rate calculated from the NEXSS mean residence time was used. The mass transfer parameters were estimated by NEXSS using the long-term average flow conditions. Then reactions in the HZs were added to isolate the impact of physical and biological processes on nitrate level in the river in order to answer the questions of the role of HZs and biogeochemical processes in HZs to nitrate concentration 275 attenuation in the surface water. These three cases are referred to as BASE, MRMT, and MRMT+BGC, respectively. Additional cases with exchange and BGC include: 1) Replacing the residence time and exchange flux with those predicted by NEXSS using seasonal flow conditions. For this case, NEXSS model was run 12 times using mean monthly streamflow conditions. The exchange fluxes and residence times now change monthly instead of being constants throughout the whole simulation. 2) Replacing the single storage zone in vertical and lateral with sub-storage zones within a storage zone, assuming a distribution of residence time to evaluate whether the commonly used single transient storage model is sufficient to model nitrate attenuation in the river; 3) Pulsed increase of nitrate in the hyporheic zones to see how surface water quality will be affected when pulsed nitrate from the legacy waste intrude into the HZs that could be caused by a flow event. To represent multiple exchange rates according to a residence time distribution, the discrete form of the probability density function of residence time is written as where τ s,j is the residence time of the j-th sub-storage zone, P j is the probability of occurrence of the residence time in the j-th sub-storage zone such that j P j = 1 (11) τ s,j can be determined from the following equation using interpolation method: in which f j is defined as the fraction or the area under the curve of the distribution function that is between τ s −∆τ s and τ s +∆τ s of the sub-storage zone within the vertical or horizontal storage zones that has an average residence time τ s . Here we assume equal fraction for each sub-storage zone, i.e., a vertical or horizontal storage is evenly divided into N s sub-storage zones. To extract N s residence times for the sub-storage zones in a given HZ with mean residence time τ m , we use the exponential 305 distribution (P = 1 − e −τs/τm) . We have one distribution for the vertical HZ and one for the lateral HZ, with means equal to the vertical and lateral residence time calculated from NEXSS. N s discrete residence time values with their corresponding probability were then extracted from this distribution for both HZs. For example, if N s = 10, the fraction of each sub-storage zone is 0.1. The average residence time of each sub-storage corresponds to a value with probability less than or equal to 0. to the small exchange flow compared to the stream discharge (Fig. 4a). Vertical HZ reached dynamic steady state in the selected reach, but the concentration in the lateral HZ is still decreasing (Fig. 4b,c), indicating continuous losing of concentrations in the HZ due to slow mass transfer.

HZ Microbial Respiration
Including the respiration reactions in the simulation results in consumption of stream nitrate (green line in Fig. 4a) because the reaction time is faster compared to vertical HZ residence times in some reaches along the river, but on the similar order. Nitrate concentration in the vertical HZ at the selected reach (RCH93) drops due to the respiration reactions, creating a concentration 330 difference between the channel and HZs, further decreasing the concentration in the channel compared with case MRMT.
Lateral HZ is diffusion limited (slow exchange), hence not much dynamics in nitrate concentration (Fig. 4c).

Seasonal exchange flux and residence time
For this case, dynamic exchange parameters were generated using NEXSS driven by long-term averaged monthly stream discharges. During the simulation, exchange parameters were updated each month. Figure 6a,b,c shows the change of vertical about 12% change of exchange flux and residence time. Using the seasonal exchange fluxes does not have a significant effect on the stream water nitrate concentration (Fig. 6d,e,f).

Multiple sub-storage zones
Zhang and Baeumer (2007) showed that the long tail behavior in stream tracer can be approximated by a small number of mass 340 transfer coefficients. We started with 10 mass transfer coefficients or sub-storage zones for both the lateral and vertical HZs.
Using exponential residence time distribution and 10 sub-storage zones in both the lateral and vertical HZs, we had 10 discrete residence times with the mean residence time in the lateral and vertical HZs from NEXSS following the procedure described in Section 3.5. Given the total lateral and vertical HZ volumes calculated using the exchange flux and mean residence time from NEXSS, we assumed 1) the exchange flux from the NEXSS estimation is equally distributed in each sub-storage zone , and 2) larger fluxes are associated with longer residence times. The second assumption is based on the theory in Elliott and Brooks (1997) that exchange flowpaths with higher fluxes penetrate deeper in the streambed with homogeneous and isotropic material and are associated with longer residence times. Compared to the single-rate simulation, simulation with multiple exchange rates and uniform hyporheic exchange fluxes removes only 3% more of nitrate in the HZ of the Columbia River reaches due to denitrification (Fig. 7a), and 40% less nitrate removal in the HZs of the Columbia River reaches (Fig. 7b) when high exchange 350 flux is associated with long residence time. These results show that the assumption of the exchange flux associated with each sub-storage zone can have a significant effect on the estimation of HZ nitrate removal and there is a need for more research on the reliable estimation of these exchange fluxes.

Pulsed nitrate increase in HZs
Contaminant can enter the stream during the upwelling of nitrate-contaminated groundwater from the legacy waste site. We 355 changed the nitrate concentration to 0.3 mmol L −1 in the vertical HZ of RCH27 and RCH101 separately at the beginning of 2011. Because the residence time in HZ of RCH27 is 30 times faster than RCH101 (Table 1 in the Supplement), the stream and HZ nitrate concentration reached steady state instantaneously (Fig. 8a,b) and the high concentration pulse propagated downstream (Fig. 8c,d) because of the longer residence time downstream. Local perturbation at the outlet (RCH101) results in a small change in stream water concentration because of the longer residence time (Fig. 8c,d). At the outlet, it takes 112 360 days for the reach concentration to return to a pre-perturbed condition. It has been observed in the field that high stream nitrate concentration than those shown in the BASE case can occur. The field observation may be explained by the contaminant intrusion as tested here. The non-equilibrium behavior shown in this test indicates that dynamic mass contribution from the groundwater should be considered in this model as it affects concentrations in the HZs.

365
There are five major wasteways carrying return flow water into the Columbia River (Fig. 2a). Bureau of Reclamation, Pacific Northwest Region, has recorded nitrate loading from these wasteways. We generated an hourly nitrate loading (Fig. 9)    into RCH100 alone. At the end of the simulation, the total denitrification at RCH77 and RCH88 are 7000 kg and 8400 kg due 375 to nitrate loading from the irrigation flow only. They are 178% and 214% of that at RCH100 where total denitrification is 3920 kg (Fig 10). In other words, HZ denitrification can remove more nitrate at RCH77 and RCH88 even though the loadings and HZ volumes are smaller compared to those at RCH100. Compared to RCH77 and RCH88, nitrate loadings into reaches with longer residence times, such as RCH100 and RCH101, have more negative effects on river water quality.

380
Our study focuses on the effect of bedform-induced (vertical) and sinuosity-driven (lateral) hyporheic exchange derived under steady-state conditions on stream water quality in large rivers using the SWAT-MRMT-R model we developed. Our studied system is dominated by inlet conditions if there are no additional sources in the HZs.

Coupled effect of physical and biogeochemical processes
Coupling reactive transport in channel and HZs, as well as MRMT in a large river using the mean residence times and exchange 385 fluxes calculated from NEXSS, our simulations show that HZs can attenuate the peak nitrate concentrations in the stream with can play an important role in storage and biogeochemical transformation in certain reaches as the vertical HZ (Fig. 5), which echoes the findings in Gomez-Velez et al. (2015).

Effect of dynamic hyporheic exchange
The Hanford Reach of the gravel-bed Columbia River exhibits frequent hydrologic mixing and steep physicochemical gradients caused by the combined effects of large, long-term river stage fluctuations driven by snowmelt runoff and short-term 395 fluctuations driven by upstream flow regulation at the Priest Rapids Dam (Graham et al., 2017;Stegen et al., 2016;Song et al., 2018). Gravel-bed river corridors, which are common worldwide (Nilsson et al., 2005;Sawyer and Cardenas, 2009;Hauer et al., 2016), exhibit coarse-textured sediments with high hydraulic conductivity. Hydrologic exchange due to dam operation and high hydraulic conductivity can drive river water into and out of the river banks, creating high frequency lateral exchanges in these type of river corridors. The study in Zhou et al. (2017) indicated the exchange fluxes in the shallow water near the river the groundwater and river stage. This is not captured by the steady-state NEXSS model, which could result in an inaccurate estimation of the attenuation capability by the HZs and recovery time for the stream water quality.

Effect of multirate mass transfer using multiple rates
Short residence times could provide faster supply of DOC into the HZs. Including the short residence time can frequently 405 change the chemical signature of the stream water in the HZ and could be used to explain, for example, the response of microbial community composition in Stegen et al. (2016). The importance of the biogeochemical processes in HZs on surface water quality could be affected by the fraction of this short residence time if Damköhler number for oxygen (Da O 2 ) is greater than 1 (Zarnetske et al., 2012). Our simulation with multiple rates shows that channel hyporheic zone nitrate removal due to denitrification could increase or decrease compared to a single rate (Fig. 7) depending how exchange flux for each sub-storage 410 zone is approximated. MRMT simulation with the uniform exchange flux approximation makes almost no difference in HZ nitrate consumption compared to the single rate simulation, while much less (40%) nitrate consumption when large fluxes are assumed to be associated with longer residence times. It has yet to be verified by field studies.

Implications for source water quality treatment
Spatial perturbation in the HZs and including point sources of irrigation return flow show that high-concentration pulses of 415 nitrate may stay in the stream before being attenuated (Figs. 8 c,d). This has implications for the effort to treat water quality in groundwater and irrigation return flows before entering the reaches with longer residence times, particularly for streams impacted by agricultural activities.

Limitations and future development
There are limitations in the current model such as 1) the lack of nitrification processes, 2) source uncertainties, e.g., lack of 420 groundwater contribution to the HZs in the model; 3) uncertainty in the exchange rates and residence times calculated by NEXSS; 4) uncertainty in the rates for biogeochemical reactions in the HZs; and 5) lack of field observations for model comparison and application. For example, the shape of residence time distributions can be drastically changed due to heterogeneity in stream sediment hydraulic conductivity as shown in the study of Pryshlak et al. (2015). We only considered denitrification processes in our study case, which is reasonable in this study as the Damköhler number for oxygen (Da O 2 ) is > 1 based on 425 residence times calculated from NEXSS and oxygen uptake rate assumed. If Da O 2 becomes less than 1, nitrification processes which dominate short residence times (Zarnetske et al., 2012) should be included in the model to better quantify the role of HZs as a net sink of nitrate. These processes can be easily implemented in PFLOTRAN. Also, we only evaluated the denitrification process using a single set of reaction rate parameters in the river corridors. However, nutrient uptake in streams can be highly spatially variable (McClain et al., 2003). Our perturbation and return flow simulations show that source term can have a big im-430 pact on stream water quality. Effect of HZs on stream water quality can be compromised by upwelling of nitrate-contaminated groundwater through adding nitrate to the stream or altering features of the hyporheic zone (Azizian et al., 2017). These factors need to be evaluated and combined with field characterization and model development in future studies.
Field experiments in large rivers are challenging due to the large spatial extent, spatial heterogeneity, and high stream discharge. Decision has to be made as to how and where to sample across the watershed. One of the advantages of the model 435 developed in this study is that it can be used to cost effectively define the reaches that are biogeochemical hot spots using physically derived parameters before conducting field experiments. The model can then be combined with experiment data collected at those locations for a better mechanistic understanding of river corridor functioning, and model improvement.
As we mentioned earlier, this model is an extension of the transient storage model. The sub-storage zones are assumed to be well mixed, and they only communicate with the stream water, not with each other. This assumption is unlike the integrated 440 model such as SWAT-MODFLOW which can rigorously resolve the redox sensitivity at each location in the modeling domain.
Future model development will consider approximating each storage zone as an advection-dispersion system like the model in Painter (2018) so that redox zonation in the storage zones can be represented. Models of different conceptualizations will be evaluated on the effect of watershed biogeochemical processing function and whether model complexity is warranted for reach scale process understanding.

Conclusion
Based on widely used open-source models SWAT and PFLOTRAN, we developed an integrated model, SWAT-MRMT-R, to account for hyporheic exchange and biogeochemical processes within HZs that are often absent in reach scale modeling. In this model, the MRMT module can incorporate a spectrum of transition times to represent solute mass transfer dynamics between the channel and storage zones. The model is flexible such that different shape of residence time distributions (e.g., exponential, 450 power-law, or lognormal distribution) and biogeochemical reactions in stream water and HZs can be easily implemented.
Exchange parameters in our study were estimated from NEXSS. Although our demonstration of the model uses mass transfer parameters derived from steady state discharge conditions, it is applicable to dynamic conditions as well once the corresponding mass transfer parameters are available. An integrated model as developed in this study can serve as a base framework to answer questions related to the significance of HZs in surface water quality in higher order river networks where characterization of 455 HZs can be challenging. Our simulation results show that it is critical to better approximate the hydrologic exchange flux that enters into each sub-storage zone given a residence time distribution because it can have a significant impact on how much nitrate can be removed through denitrification in HZs. Effort are needed to treat water quality in groundwater and irrigation return flows before entering the reaches with longer residence times. It is worth noting that this model has a multitude of parameters at reach scale. For real-world, large-scale studies model validation at multiple points within the study area is