High-Performance Simulation of Dynamic Hydrologic Exchange and Implications for Surrogate Flow and Reactive Transport Modeling in a Large River Corridor

Hydrologic exchange flows (HEFs) have environmental significance in riverine ecosystems. Key river channel factors that influence the spatial and temporal variations of HEFs include river stage, riverbed morphology, and riverbed hydraulic conductivity. However, their impacts on HEFs were often evaluated independently or on small scales. In this study, we numerically evaluated the combined interactions of these factors on HEFs using a high-performance simulator, PFLOTRAN, for subsurface flow and transport. The model covers 51 square kilometers of a selected river corridor with large sinuosity along the Hanford Reach of the Columbia River in Washington, US. Three years of spatially distributed hourly river stages were applied to the riverbed. Compared to the simulation when riverbed heterogeneity is not ignored, the simulation using homogeneous riverbed conductivity underestimated HEFs, especially upwelling from lateral features, and overestimated the mean residence times derived from particle tracking. To derive a surrogate model for the river corridor, we amended the widely used transient storage model (TSM) for riverine solute study at reach scale with reactions. By treating the whole river corridor as a batch reactor, the temporal changes in the exchange rate coefficient for the TSM were derived from the dynamic residence time estimated from the hourly PFLOTRAN results. The TSM results were evaluated against the effective concentrations in the hyporheic zone calculated from the PFLOTRAN simulations. Our results show that there is potential to parameterize surrogate models such as TSM amended with biogeochemical reactions while incorporating small-scale process understandings and the signature of time-varying streamflow to advance the mechanistic understanding of river corridor processes at reach to watershed scales. However, the assumption of a well-mixed storage zone for TSM should be revisited when redox-sensitive reactions in the storage zones play important roles in river corridor functioning.


INTRODUCTION
The concept of river corridors that integrate in-stream transport with lateral and vertical exchange across interfaces (Harvey and Gooseff, 2015) has been increasingly recognized due to their valuable ecological functions. The surrounding sediments along the river corridor, or the hyporheic zones (HZ), directly relate with the stream surface water in terms of activity and hydrologic connectivity across microscale, streamreach scale, and catchment scale (Boulton et al., 1998). Bidirectional mixing, or hydrologic exchange flow, at the interface of stream and HZs is important in influencing stream water quality, water temperature, fish conservation, microbial communities, dissolved organic matter processing, and groundwater contamination (Boulton et al., 1998;Edwardson et al., 2003;DiStefano et al., 2009;Febria et al., 2010;Wexler et al., 2011;Dickson et al., 2012;Kawanishi et al., 2017;Haddou et al., 2018;Harjung et al., 2018;Lewandowski et al., 2019).
River corridors are heterogeneous in space, non-stationary, and involve multiscale feedbacks (Ward and Packman, 2019). Previous studies suggested that research on river corridor processes has suffered from the lack of the long-term (weeks to years) investigations on the dynamic driving forces (e.g., streamflow, evapotranspiration) of hyporheic processes across different scales (Boano et al., 2007;Helton et al., 2012;Larsen et al., 2014;Anibas et al., 2018;Galloway et al., 2019;Zheng et al., 2019). A recent kilometer scale simulation of a major gravel bed river corridor along the Columbia River by Zachara et al. (2020) suggested that hyporheic zones must be frequently monitored under different flow conditions for water quality concerns because of the complicated multidirectional flow behaviors. Studies accounting for the dynamic hydrologic processes on hyporheic exchange and biogeochemical reactions have started to emerge in recent years (Larsen et al., 2014;Rahimi et al., 2015;Gomez-Velez et al., 2017;Dwivedi et al., 2018;Liu and Chui, 2018;Singh et al., 2019;Zheng et al., 2019;Kruegler et al., 2020). However, these studies mainly focused on the spatial and temporal patterns of HEFs and biogeochemistry for synthetic individual bedform or two-dimensional cross-sections of single geomorphic features along the river corridor. Scaling these observed complex patterns from small scale localized sampling to reach scale and basin scale investigation is challenging as they are influenced by geomorphic and hydrologic dynamics (Helton et al., 2012), channel morphology, bed roughness, and permeability (Triska et al., 1989;Hassan et al., 2015), hydrologic connectivity (Datry et al., 2008), and vegetation feedbacks (Magliozzi et al., 2019).
On the other hand, reach scale river corridor exchanges are often quantified using the transient storage models (TSM) (e.g., Bencala, 1983;Harvey et al., 1996;Runkel et al., 1998;Briggs et al., 2009). Despite its simplicity, TSM has been extensively used to study riverine solute transport since its first introduction fifty years ago. TSM is a simple expansion of the one-dimensional instream transport equation to include bilateral exchange processes between the main-channel and the well-mixed transient storage zones (Bencala, 1983). This model assumes that the amounts of water transferred between the stream and subsurface are insignificant at the reach scale and the water transfer only influences solute transport (Bencala et al., 2011). The exchange process is often represented by a single-rate mass transfer. In the past, the hyporheic exchange rates were obtained from the calibration of field tracer experiments as they are impractical to measure directly over large areas (Boano et al., 2014). Reach scale tracer investigations can also be effective in providing bulk reactivity for a particular biogeochemical constituent of interest (Harvey and Fuller, 1998;Jonsson et al., 2004;Harvey et al., 2013;Lemke et al., 2013), but they rely on data fitting from the observations and mechanistic process understanding cannot be gained (Boano et al., 2014). There were a few examples of TSM applications to reactive tracers (Bencala, 1983;Harvey and Fuller, 1998;Gooseff et al., 2005;Knapp and Cirpka, 2017;Kelleher et al., 2019), but they only focused on a single reactive component.
A recent review article by Ward and Packman (2019) calls for a framework and theoretical models of river corridor exchange to accurately predict the hydrologic exchange in the river corridor at the scales of stream reaches and entire networks. Lewandowski et al. (2019) also highlighted the research focus on HEF processes that link hydrological, biogeochemical and ecological processes. In this research, we investigate how mechanistic understanding from high resolution three-dimensional (3D) simulations of flow and reactive transport induced by multiple years of transient stream discharge in a large gravel bed river corridor can be used to inform the development of reach scale reduced order models, such as the transient storage models. We use synthetic data generated by high-resolution 3D physically-based simulations to derive and test a surrogate TSM model and identify future development needs. Given the physical basis of the highresolution model it is reasonable to assume that it is a plausible representation of a realistic system, although we recognize that it is not a precise representation of any particular real system. While TSM models can be parameterized and tested using field tracer tests, the size and flow rate of the Columbia River renders such tests impractical in this case. Therefore, testing the performance of a simpler empirical model (the TSM) against the high-resolution model outputs is one reasonable way to assess potential strengths and weaknesses of the TSM for this large system. The results of the 3D simulations were used to parameterize the exchange coefficients of the TSM and were compared with the TSM simulation results, for both conservative and reactive transport, to identify TSM model weaknesses and future model development needs. We hypothesize that transient streamflow and riverbed heterogeneity have significant effects on the zonation of HZs, residence time of HEFs within HZs, and biogeochemical reaction rate constants at the reach scale.

Study Area
The study river corridor is located at the 100-H area along the Hanford Reach (an eight-order stream) of the Columbia River in the southeastern Washington State ( Figure 1A). The Hanford Reach provides a refuge for native fish, riverine wildlife and plants (Rickard and Gray, 1995). The bottom substrate of the river consists mainly of cobbles and boulders, and the flow rate in the Hanford Reach is regulated by upstream dams in response to hydroelectric power demand (Dauble et al., 1989), ranging between 1,900 and 5,600 m 3 /s depending on season (Zachara et al., 2020). The water level fluctuates up to 2 m daily due to discharge variation from the Priest Rapids Dam (Arntzen et al., 2006). At our selected location, the river meanders north and then southeast (Figure 1). The unconfined aquifer of the river corridor consists primarily of the Hanford Formation composed of Pleistocene-age flood deposits and the underlying late-Miocene to late-Pliocene Ringold Formation formed from river and lake sediments (Thorne et al., 2006;Shuai et al., 2019). The Ringold Formation in the 100-H area is broadly divided into Ringold E unit, Ringold Taylor Flats unit, and Ringold Lower Mud unit. The geologic unit of the simulated domain is shown in Figure 1B and their hydrologic properties are summarized in Table 1. The study domain occupies an area of 50.864 km 2 . The domain is discretized into 748 × 680 × 30 (15,096,000) elements of 10 m resolution in the horizontal direction and 1 m in the vertical direction. The bottom elevation is 90 m above sea level. The riverbed is delineated using a combination of LiDAR and other measured bathymetric elevation (Coleman et al., 2010).

Subsurface Flow and Reactive Transport Simulation
To investigate the hydrologic exchange fluxes at the riverbed driven by multi-year stream discharge, we used the state-of-theart massively parallel subsurface flow and reactive transport code, PFLOTRAN. PFLOTRAN is an open source code, which solves a system of non-linear partial differential equations to describe multiphase, multicomponent and multiscale reactive flow and transport in porous materials (Hammond et al., 2014). For this study, the Richards Equation was solved for the variably saturated flow with the van Genuchten (Van Genuchten, 1980) and Burdine (Burdine, 1953) constitutive relations for saturation and relative permeability, respectively. The simulation period is from 2013 to 2015, which covers two normal flow years and a low flow year in 2015.

Initial and Boundary Conditions for Flow
The initial pressure of the modeled domain was interpolated from the 2012 results of a coarser reach scale simulation in Shuai et al. (2019). The bottom boundary of the domain was prescribed as no flow condition as it was constrained by the low permeable Ringold units. The top boundary was also set as no-flow as the surface recharge is negligible in the semiarid climate zone (Rockhold et al., 1995). The inland boundaries of the domain were prescribed as hydrostatic conditions, using transient pressures interpolated from the reach scale simulation in Shuai et al. (2019). To represent a low permeable thin alluvium layer at the sediment-water interface, a heterogeneousconductance type boundary condition was imposed on each grid cell of the riverbed. The conductance type condition is similar to a seepage face, which dampens the effect of river stage fluctuations (Hammond and Lichtner, 2010). The flow rate for River boundary condition 2.58 × 10 −4 8.87 × 10 −6 3.59 × 10 −4 1.0 × 10 −5 this type of boundary is proportional to the head difference, with conductance as the constant of proportionality. Hourly river stages simulated using the Modular Aquatic Simulation System in two dimensions (MASS2) simulator (Perkins and Richmond, 2004), as described by Niehus et al. (2014), were assigned to each river boundary cell.

Reactive Transport
To evaluate the effect of dynamic HEFs on reaction rates at the reach scale, we added a conservative tracer and the following two-step denitrification reactions and an oxidative respiration reaction based on previous studies by Song et al. (2018) at a nearby area in the reactive transport simulation: Cell decay or lysis of biomass  and turnover back to CH 2 O was also considered in the model. Reaction rate formulations follow a modified dual-Monod type. The maximum uptake rate of O 2 , NO 3 − , and NO 2 − are 84.78, 28.26, and 23.28 d −1 , respectively. The initial and boundary conditions are shown in Table 2. NH 4 + was assumed to be unlimited. A conservative tracer with a concentration of 10 −3 mol/L was assigned at the river boundary to indicate the concentration of river water. The four lateral boundaries were treated as tracer free.

Case Scenarios
We ran simulations for five case scenarios: • Case 5. Flow and reactive transport for year 2014 using the end of the simulation in year 2013 of Case 1 as the initial condition for flow, and setting initial and boundary conditions of chemical components to those in Table 2.
Cases 1 and 2 were used to examine the effect of dynamic discharge and riverbed heterogeneity on HEFs. It is challenging to quantify the spatial distribution of conductance over the reach scale. A recent effort by Hou et al. (2019) provided a map of the riverbed sediments along the Hanford Reach based on the integration of diverse observations with numerical simulations of river hydrodynamics using facies classifications. This map enabled us to assign unique conductance for each of the four facies identified in the study area (Figure 2) to evaluate how significant the HEFs can be affected by riverbed heterogeneity. The conductance is 1.07 × 10 −13 , 9.11 × 10 −13 , 3.97 × 10 −12 , and 8.01 × 10 −12 m, for facies 1, 2, 3, and 4, respectively. Facies 1 to 4 correspond to increasing mean shear stresses from river hydrodynamics simulation (Hou et al., 2019). The homogeneous conductance is the geometric mean of these heterogeneous values. Cases 1 and 2 were used to define the hyporheic zones in the river corridor. Results from Cases 3, 4, and 5 were used for reduced order model development for reach-scale simulation. We generated synthetic reach-scale datasets of concentrations in the hyporheic zones for tracer and reactive components from cases 3 to 5 treating the whole simulation domain as a single river corridor of a reach. All simulations were run on the EMSL (the William R. Wiley Environmental Molecular Sciences Laboratory) cascade high performance computing cluster using 1,596 process cores. Hourly outputs of around 9.6 terabytes for each year were recorded for analysis and post-processing.

Reduced-Order Model
The transient storage model (TSM) (Bencala, 1983), a simplified reduced order model, has been extensively used to quantify hyporheic exchange flows based on tracer injection experiments (Bencala et al., 2011). The partial differential formulation of the model is shown in the following equations (Harvey et al., 1996;Boano et al., 2014;Larsen et al., 2014):  (Harvey et al., 1996;Larsen et al., 2014). The storage zone size and the residence time are the two key parameters required to run the TSM model. In order to use the conservative tracer simulation to provide meaningful parameter estimates for the reduced order model, massively parallel forward particle tracking algorithm was used to find the residence time distribution and flow paths of the hydrologic exchange flows. Particles were released from 100,000 randomly sampled locations on the river boundary at 1,000 time points randomly selected between October 2013 and September 2014.
We assumed a single transient storage zone or hyporheic zone (HZ) for the whole river corridor. After 1-year spinup, the volume of the storage zone was defined as the region with 95% of the stream water, or tracer concentration of 0.0095 mol/L, following the geochemical definition in the literature (Triska et al., 1989;Singh et al., 2019) using the results of Cases 1 and 2 in 2014. The total mass of solutes of Case 3 and Case 4 within the whole transient storage zone volume defined in Case 1 and Case 2 were calculated, respectively. The total mass in each case was then used to calculate the hourly whole reach HZ concentration (mass divided by HZ volume) response. The same procedure was performed for the reactive components, i.e., CH 2 O, NO 3 − , and O 2 in Case 5.
We amended the TSM model with reactions and the following equation was used to calculate the hyporheic zone reactive component concentrations at the reach scale by treating it as a batch reactor. Concentrations of reactive species are fixed in the river (i.e., Equation 4 is not solved) as shown in Table 2.
where R i is the reaction rate of reactions shown in Equations (1-3), and µ i is the stoichiometric coefficient of a species in reaction i, which is negative for reactants and positive for products. For example, the above equation for C 5 H 7 O 2 N and NO 2 − can be written as Equation (7) and Equation (8), respectively.
Species concentrations predicted by the reaction amended TSM represented by Equation (6) for CH 2 O, NO 3 − , and O 2 , will be compared with the average concentrations within the HZ calculated from the PFLOTRAN simulation. Equation 6 is solved using the PFLOTRAN batch reaction module.
Global sensitivity test using the Fourier Amplitude Sensitivity Test (FAST) method will be performed for the TSM as we expect there will be disagreement between the two models for reactive species.

Hydrologic Exchange Flow
We summed up the lateral and vertical hydrologic exchange flow rates for Case 1 and Case 2 after 1-year spinup. The maximum lateral flow is about a tenth of the vertical flow (Figures 3A,B). From the hourly PFLOTRAN output during the 2-year simulation, we found that the river gains water from the lateral flow 83.8% of the time for Case 1 and 85.5% of the time for Case 2. The river gains water from the vertical flow 54.7% of the time for Case 1 and 49.8% of the time for Case 2. While the percentages of times for upwelling (denoted as positive value) and downwelling do not make a big difference for the two cases, lateral flows contributed more to the river water for Case 2 with a mean of 0.46 m 3 /s due to the heterogenous conductance compared to Case 1 with a mean of 0.21 m 3 /s ( Figure 3A). On the other hand, vertical flow for Case 1 contributed more to the river, with a mean flow rate to the river of 0.47 m 3 /s compared to 0.1 m 3 /s for Case 2. Vertical downwelling flow rates were stronger than upwelling for both Case 1 and Case 2 ( Figure 3B). There was a seasonal decrease of lateral flow to the river ( Figure 3C). Vertical upwelling was strong in the fall, and downwelling is dominant in the wet season ( Figure 3D). Compared to the Columbia River discharge (1,900-5,600 m 3 /s), the hydrologic exchange flow has little effect on river hydrodynamics and transport.

Volume of the Hyporheic Zone
The volume of the hyporheic zone (HZ) for Case 1 and Case 2 was calculated by summing up the pore volumes of those grid cells with tracer concentration equals to or > 0.0095 mol/L (i.e., ≥ 95% river water). Hourly HZ volumes varied between 5.7 × 10 6 and 1.23 × 10 7 m 3 for Case 1, and between 6.11 × 10 6 and 1.27 × 10 7 m 3 for Case 2. The distributions of the HZ volume displayed pronounced multimodal character as shown in the bean plot created using the R package "beanplot" (Figure 4A). The difference of mean HZ volume between the two cases was small. The seasonal expansion and contraction of HZ follow the dynamics of stream discharge ( Figure 4B). HZ volume for  Case 2 was larger than Case 1 when streamflow increased, and smaller when streamflow decreased due to the high riverbed conductance of Facies 4. There was a temporal lag of the peak HZ volume compared to the stream discharge due to the variation in residence time of the hydrologic exchange fluxes.

Reach-Scale Residence Time
A mean residence time of all released particles at a selected time was calculated by weighting the residence time of each particle with its releasing flux. A total of 1,000 reach-scale mean residence times were calculated with the 100 million particles. Compared to Case 1, the longest residence time in Case 2 was shorter than the long-term average of Case 1 (Figure 5). Heterogeneous conductance decreased the long-term mean residence time by half, with a mean of 1,639.35 h for Case 1, and 782.89 h for Case 2. Case 1 residence time had a larger variability than Case 2, ranging from 259 to 2,570 h, compared to 291 to 1,520 h for Case 2. In general, the residence time is positively related to stream discharge ( Figure 5B).

Model Validation With Conservative Tracer
Summing up the volume of the grid cells identified as hyporheic zone in year 2014 for Case 1 and Case 2 at every time step, we obtained a reach-scale hyporheic zone volume for  each Case. The mass of conservative tracer simulated in year 2014 from Case 3 and Case 4 in those hyporheic zone grid cells identified in Case 1 were then summed up. Dividing the total mass of conservative tracer by the reach-scale hyporheic zone volume, we obtained a synthetic hourly response of the conservative tracer concentration in the hyporheic zone driven by the dynamic stream discharge. The exchange rate coefficient in Equation (6) is now time-varying. At every time step, the residence time was interpolated from those estimated using particles released during 2014 for both Cases 1 and 2. Our simulation showed that the transient storage model (TSM) without reactions could explain the synthetic tracer data with an R-squared of 0.97 for Case 1 and 0.95 for Case 2 (Figure 6). However, TSM underpredicted the concentration for Case 1, and overpredicted for Case 2 compared to the synthetic PFLOTRAN results.

Model Validation With Reactive Components
As shown in Figure 5, the minimum reach-scale residence time was > 200 h. Because the time scales of the reactions shown in section Reactive transport were < 1 h, the reaction Damköhler number, Da was >> 1. This means hyporheic exchange flow paths spent too much time in storage, and they were transport limited . Case 5 was run for year 2014 using the end of 2013 simulation of Case 1 as the initial condition for flow and chemical initial and boundary conditions as shown in Table 2. We summed up the mass of CH 2 O, NO 3 , and O 2 for Frontiers in Water | www.frontiersin.org the hyporheic grid cells identified in Case 1 for 2014, and then normalized to the hyporheic volume calculated for Case 1 at each time step of year 2014 to create a temporal synthetic reach-scale chemical concentration data set. Running the reaction amended TSM using the same exchange rate coefficients as for the tracer validation in Case 3 and the same reaction kinetics in Case 5, the concentration profiles for CH 2 O, NO 3 , and O 2 are shown in Figure 7.
R-squared for CH 2 O, NO 3 , and O 2 are 0.19, 0.99, and 0.13, respectively. The TSM model has a very good performance for NO 3 except for at low concentrations after 2,000 h (Figure 7b), but biased high for CH 2 O (Figure 7a) and O 2 (Figure 7c) compared to the synthetic data. However, the high variability of CH 2 O exhibited in the synthetic data was simulated using TSM. The high variability was mainly caused by the balance between the production by biomass decay and exchange, and the consumption by the respiration reactions (Figure 7d). The concentration change with time for NO 3 and O 2 are relatively smooth. The hyporheic exchange contributed to the HZ concentrations of CH 2 O (Figure 7d), and O 2 (Figure 7f) during the simulation, but was a sink for HZ NO 3 till 2,000 h (Figure 7e). NO 3 consumption was dominated by hyporheic exchange up to 2,000 h. The TSM model performance for NO 3 at low concentrations started to deteriorate when production due to exchange balances the reaction consumption.

Global Sensitivity Analysis
Apart from the residence times, there are 8 reaction kinetic parameters involved in the TSM with reactions. To evaluate whether there exists a single set of rate parameters that can explain the results obtained from the PFLOTRAN simulations using the TSM, we ran simulations using a range of rate parameters previously reported in the literature. Parameter importance was evaluated to identify which rate parameters we should pay our attention to for TSM parameterization. Rate parameters were sampled from 6.24 to 240 (1/d) for the maximum microbial reaction rate of nitrate (k NO3 ), 0.001-0.05 (mmol/L) for the Monod half-saturation constant of nitrate (K a,NO3 ), 0.0833-0.2583 (mmol/L) for the half-saturation constant of CH 2 O (K d,CH2O ), 4.8-139.2 (1/d) for the maximum microbial reaction rate of aerobic respiration (k O2 ), 0.001-0.0969 (mmol/L) for the half-saturation constant of O 2 (K a,O2 ) based on parameter bounds in Zarnetske et al. (2011) surveyed from the literature. The lower bound for the maximum microbial reaction rate of nitrite was set to one-tenth of that in Killingstad et al. (2002), and the rest of rate constant bounds for nitrite were the same as nitrate reduction reaction. We also assumed 0.242-3.1536 (1/d) for the biomass decay rate. One thousand six hundred samples were generated using the Fourier Amplitude Sensitivity Test (FAST) method implemented in the Python sensitivity analysis library (SALib) (Herman and Usher, 2017).  The TSM results from the bulk of the simulations did not fit the PLOTRAN results well (Figure 9). To check if there is a common set of effective rate parameters for all three reactants, we narrowed down to 5% of the model runs that can explain the observed synthetic concentrations for each reactive component. Plotting the parameters, we were not able to find a common set of important parameters for all of them (Figure 10). The parameters also varied with time. For instance, to explain the CH 2 O concentration, k NO3 is > 100 (1/d) during the earlier time and smaller than 50 (1/d) between 1,000 and 2,000 h ( Figure 10A). We only showed the parameter space for k NO3 as an example. It is challenging to find a constant value for the other parameters too.

DISCUSSION
Do Spatial and Temporal Heterogeneity Matter for Reach-Scale Modeling?
There have been an increasing number of numerical studies emphasizing the impact of spatially and temporally variable conditions on the hyporheic exchange within river corridors. However, none were carried out at the reach-scale accounting for both the long-term dynamics in driving force and the smallscale hydromorphic heterogeneity. This study investigated the response of river corridor hydrologic exchange as a result of time varying streamflow and riverbed heterogeneity for a reach in an eighth-order stream. Riverbed conductance heterogeneity increased the mean lateral hydrologic exchange flow into the river by a factor of 2.2 and decreased the vertical exchange flow into the river by a factor of 4.7. On the other hand, dynamic streamflow changed the dominant downwelling in the wet season to upwelling dominant in the dry season for the vertical flow ( Figure 3D). As the upwelling becomes stronger in the dry season (Figure 3D), the hyporheic volume is significantly reduced (Figure 4B), which echoes the findings in Boano et al. (2008) that the upwelling of subsurface water can significantly reduce the volume of hyporheic sediments that receive water from the stream. HEF variability caused by streamflow is larger than that by riverbed conductance for this case study.
Hyporheic volume defined using the geochemical definition showed that the largest expansion at the peak streamflow can be twice as big as that during the low streamflow. The discrepancy in hyporheic volume due to riverbed conductance heterogeneity is not significant compared to that caused by streamflow variability. It is consistent with the conclusion in Sawyer and Cardenas, 2009 that permeability heterogeneity in riverbed sediment is more important in changing the residence times than the volume of the hyporheic zone. Shorter mean residence time can be caused by riverbed heterogeneity compared to the homogenous simulation, which has also been previously reported (Salehin et al., 2004;Tonina et al., 2016). The riverbed conductance heterogeneity and streamflow variability both have impacts on the reach scale residence times. The difference caused by conductance heterogeneity can be comparable to the difference caused by streamflow variability.
Can a Transient Storage Model Be Developed for Reach-Scale Modeling?
Though not perfect, we showed that it is possible to derive the TSM based on the residence time and synthetic dataset simulated using high resolution flow and reactive transport. We developed a simplified transient storage model (TSM) integrating the streamflow dynamics and riverbed heterogeneity that have been reflected in the residence time to simulate reach-scale processes. The TSM model was agreeable simulating the synthetic reach-scale tracer data. But the discrepancy due to heterogeneity was exaggerated due to the inaccuracy predicting the synthetic observation. The mismatch between the synthetic observations and TSM may be due to the definition of hyporheic zone or some HEF pathways that contributed significantly to the tracer concentrations were not accounted for when deriving the mean residence times. The high bias in O 2 and CH 2 O coincided with the low bias of tracer, which may be caused by the cumulative effect of reaction and transport that cannot be simply represented by using the same reaction parameters apart from the aforementioned possible causes.
When including reactions in the TSM, there is a surprisingly good agreement between the TSM and the synthetic observation of nitrate. TSM results shows that the nitrate loss due to the biological reactions varied from 34 to 66% of the total loss combined with hyporheic exchange. For this simulation, nitrate from the hyporheic zone was mainly a source to the river water until HZ nitrate concentration is lower than river concentration. The high bias of CH 2 O and O 2 concentrations suggest that the agreement between the TSM and the observations for nitrate might give us a bulk nitrate uptake rate, but it is not enough to explain the other two reactive components. Because the tracer concentration was being underpredicted by the TSM as shown in Figure 6, the high bias of CH 2 O and O 2 could not be caused by the hyporheic exchange, but effective reaction rates. The results also indicate that the well-mixed assumption in the HZs for TSM should be revisited, especially for redox-sensitive reactions. Adding one-dimensional flow and reactive transport in the HZ might improve the model simulation.

What Did We Learn From the Reduced Order Model?
One of the questions asked in a recent commentary by Knapp and Kelleher (2020) was whether the processes estimated from conceptual parameters are interpretable for TSM. We showed that the TSM model parameters can be interpretable. Our TSM parameters carry the signatures of dynamics streamflow and subsurface heterogeneity. The simplified TSM can assist in uncertainty analysis efficiently to understand the controlling parameters for reach-scale system behavior. Sensitivity analysis on reaction rate parameters showed that they can be time dependent, which is caused by inaccurate conceptualization, such as well-mixed storage zones. The maximum change in cumulative consumption of nitrate is only 2.5% if a constant mean residence time was used for TSM because of the long mean residence time. This suggested that the dynamic exchange rate coefficient was not necessary for this study case. However, the difference between the simulations using a constant vs. time varying residence times became larger if we decreased the residence times by a factor of 10, reaching 7.8%.
High resolution models can give more mechanistic details of a study system, but it is not realistic for every river corridor as it is quite expensive, especially for uncertainty analysis. Combined with TSM models amended with biogeochemical reactions and machine learning (future work) from the high resolution model simulation results, it is possible to better constrain the TSM models and improve the mechanistic understanding of the river corridors and better inform field investigations so that the vulnerability of watersheds to water quality changes due to a variety of stressors can be identified and protected (Ranalli and Macalady, 2010).

CONCLUSIONS
Using high performance flow and reactive transport simulations and a massively parallel particle tracking method, we calculated the dynamic hydrologic exchange fluxes and residence time distributions driven by long-term time-varying stream discharges within a large river corridor, located along an eighth-order stream. While it is possible to interpret the synthetic tracer data at the reach-scale using a transient storage model integrating the signatures of dynamics stream discharge and subsurface heterogeneity, there are still challenges representing multiple reactive components using simplified TSM that's amended with reactions. A bulk reaction rate for a single reactive component may be far from reflecting the true mechanisms in the system of interest. A transferrable, novel multiscale modeling framework integrating multiscale characteristics from high resolution simulations and modified TSMs, e.g., accounting for multirate mass transfer, dynamic reaction kinetics, or depthresolved flow and transport, is in need in order to improve mechanistic predictions of water quality in streams and map the watersheds of vulnerability for restoration.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
YF did model simulations and wrote the paper. XS did particle tracking. HR, WP, PS, and MR developed the dataset for model boundary conditions. All co-authors contributed to the writing and model discussion of the paper.