Effects of Successive Peak Flow Events on Hyporheic Exchange and Residence Times

Hyporheic exchange is a crucial control of the type and rates of streambed biogeochemical processes, including metabolism, respiration, nutrient turnover, and the transformation of pollutants. Previous work has shown that increasing discharge during an individual peak flow event strengthens biogeochemical turnover by enhancing the exchange of water and dissolved solutes. However, due to the nonsteady nature of the exchange process, successive peak flow events do not exhibit proportional variations in residence time and turnover, and in some cases, can reduce the hyporheic zones' biogeochemical potential. Here, we used a process‐based model to explore the role of successive peak flow events on the flow and transport characteristics of bedform‐induced hyporheic exchange. We conducted a systematic analysis of the impacts of the events' magnitude, duration, and time between peaks in the hyporheic zone's fluxes, penetration, and residence times. The relative contribution of each event to the transport of solutes across the sediment‐water interface was inferred from transport simulations of a conservative solute. In addition to temporal variations in the hyporheic flow field, our results demonstrate that the separation between two events determines the temporal evolution of residence time and that event time lags longer than the memory of the system result in successive events that can be treated independently. This study highlights the importance of discharge variability in the dynamics of hyporheic exchange and its potential implications for biogeochemical transformations and fate of contaminants along river corridors.


Introduction
River discharge and stage are characterized by significant temporal variability and frequent peak flow events (Bernard-Jannin et al., 2016;Gomez-Velez et al., 2017;Krause et al., 2011;Mojarrad et al., 2019;Salazar et al., 2014;Sawyer et al., 2009;Singh et al., 2019;Trauth & Fleckenstein, 2017;Vivoni et al., 2006;Wu et al., 2018). This variability results from natural drivers (snow and precipitation events) and anthropogenic interference (such as artificial inputs from waste water treatment plants or dam operations). Perturbations in river stage and discharge-alter the pressure distributions at the interface between the water column and the streambed, the main driver for hyporheic exchange, and thus lead to exchange fluxes several orders of magnitude higher than the ones under the base flow conditions (Gomez-Velez et al., 2017;Singh et al., 2019;Wu et al., 2018). In addition to river stage fluctuations and pressure variations, continuous exchange of water, solutes, and energy between the water column and the hyporheic zone is controlled by interactions between geomorphological settings, channel gradient, hydraulic conductivity, sediment heterogeneity, and spatial variability in heads at the sediment-water interface (SWI) and preferential flow paths Gomez-Velez & Harvey, 2014;Lotts & Hester, 2020;Marzadri et al., 2016;Menichino & Hester, 2015;Tonina & Buffington, 2011). Interactions between these driving mechanisms determine the hydrodynamics and residence times within the hyporheic zone.
Successive peak flow events can cause variability in physiochemical characteristics of streambed fluids (Bruno et al., 2009;Fritz & Arntzen, 2007;Hinton et al., 1997;Inamdar et al., 2004;Krause et al., 2011). The main impacts of peak flow events relate to the activation of deeper subsurface flow paths, enhanced reactivity owing to the faster transport of solutes into the reactive zones and contaminant attenuation. For example, when denitrification in anaerobic zones of the hyporheic zone is limited by the supply of labile dissolved organic carbon (DOC) (Zarnetske et al., 2011), peak flow events can facilitate higher and faster transport of DOC into the streambed and thus enhance nitrogen cycling (Hinton et al., 1997;Inamdar et al., 2004). Fritz and Arntzen (2007) found that the higher influx of river water during increased water levels resulted in lower uranium concentrations due to dilution effects. Bruno et al. (2009) showed that peak flow event-induced disturbance can have negative effects for hyporheic invertebrates because sharp increase in discharge can lead to scouring, disturbing the benthic invertebrates and salmonid eggs, and eventually impacting overall ecosystem functioning.
Previous modeling studies have analyzed the impact of dynamic stream flow on groundwater-surface water exchange flow processes and hyporheic exchange (Boano et al., 2013;Dudley-Southern & Binley, 2015;Malzone, Anseeuw, et al., 2016;McCallum & Shanafield, 2016;Schmadel et al., 2016;Ward et al., 2013Ward et al., , 2018, and for this, primarily focused on the effect of individual discharge events on hyporheic exchange processes (Gomez-Velez et al., 2017;Liang et al., 2018;Singh et al., 2019;Trauth & Fleckenstein, 2017;Wu et al., 2018). These studies led to the identification of dominant drivers and controls of hyporheic exchange flows during transient stream flow conditions. For example, Harvey et al. (2012) found that larger-magnitude and longer-lasting events have higher potential to activate hyporheic flow paths and engage in solute and fine particulate transport. The modeling of dynamic hyporheic exchange fluxes by  showed that the volume of the hyporheic zone and the hyporheic exchange fluxes are modulated by annual and storm-induced groundwater fluctuations. Schmadel et al. (2016) investigated the importance of diel hydrologic fluctuation and controls such as hillslope lag, amplitude of the hillslope, and cross-valley and down-valley slopes on hyporheic flow path and residence times. McCallum and Shanafield (2016) found alterations in the residence time distributions of bank inflows and outflows for different flow events, and Gomez-Velez et al. (2017) explored the role of high-discharge events on the spatial and temporal evolution of river bank storage and sinuosity-driven hyporheic exchange. In addition, a field-based study incorporating successive storm events by Dudley-Southern and Binley (2015) showed that such events could lead to a reversal of the vertical hydraulic gradient and can enhance mixing up to 30 cm within the streambed. Similar outcomes were also found by Krause and Bronstert (2007) and Munz et al. (2011). To the authors' knowledge, none of the previous modeling studies incorporated the effects of successive hydrological events as they commonly occur in natural settings (e.g., driven by storm events, cyclic behavior of wastewater treatment plants, and dam operations); therefore, mechanistic understanding of the impacts of successive peak flow events remains elusive. Moreover, the effect of the interaction between two successive events on water and solute transport within the hyporheic zone remains largely unaccounted in these studies.
The goals of this multiparametric study are to compare and contrast the hydrodynamics and solute transport characteristics of hyporheic exchange driven by a single peak flow event and two consecutive flood events. This is done by systematically modeling the flow of water, transport of a conservative tracer, and the evolution of mean residence times. We provide a systematic approach to decipher the potential impacts of successive peak flow events on hyporheic exchange flows, using reduced complexity models of idealized, uniform, and single types of bedform-induced hyporheic exchange. The term reduced complexity indicates that the model formulation assumes and captures first-order drivers and controls of the exchange process, ignoring some of the higher-order complexities such as heterogeneity in the sediments . These analyses allow us to develop a comprehensive understanding from a large number of simulations enabled by the reduced computational effort. We evaluate the impacts of two successive peak flow events, the first event is referred to as antecedent peak flow event, and the second one as the subsequent peak flow event, by studying the effects of three different parameters, namely (i) the time lag between two peak flow events (t lag ); (ii) the magnitude ( H p 1 − H p 2 ); and (iii) the duration of subsequent flow event (t d1 − t d2 ) in relation to the preceding event. This study provides important mechanistic understanding of the flow and transport behavior in dynamically forced hyporheic zones,with relevance for biogeochemical cycling at the interface of surface and groundwater systems. Our results can therefore be used as a guide for designing field experiments and the interpretation of tracer studies in dynamic flow systems, as well as interpretation of hyporheic zone turnover efficiencies.

Methodology
A process-based model was developed in this study and used for performing model simulations following an approach of six major steps (1) setup of simulation framework including streambed sediment geometry and parameterization of properties; (2) generation of successive peak flow event stage hydrograph; (3) inclusion of flow in porous media model; (4) solute transport model to track the response of hyporheic exchange and dynamic hyporheic zone development; (5) implementation and simulation of residence time model; and (6) simulation of numerical breakthrough curves and development of model scenarios.

Conceptual Model
We use a simplified conceptualization of a streambed-river interface to systematically explore the impact of successive peak flow events on bedform-driven hyporheic exchange (Elliott & Brooks, 1997;Singh et al., 2019;Stonedahl et al., 2010). For this purpose, we implement detailed flow, transport, and residence time models in the model domain. The modeling domain (Ω) represents homogeneous and isotropic stream sediments. It is bounded at the top by the SWI (∂Ω SWI ), which is assumed sinusoidal Z SW I ðxÞ ¼ The numerical solution of the proposed model is implemented into COMSOL Multiphysics for simulation. Mesh-independent solutions are achieved with a resultant mesh of approximately 40,000 triangular elements with telescopic refinement closer to the sediment-water interface. This refinement is required in order to capture the effect of local, fast-flowing hyporheic circulation cells and have accurate flux integrals along the boundaries.

Successive Peak Flow Events Hydrograph Generation
Individual peak flow pulses or successions thereof are used as simplification to represent the dynamic nature of river discharge. The deterministic stage hydrograph is modeled with an asymmetric curve as proposed by Cooper and Rorabaugh (1963): where H p1 , t p1 , and t d1 describe the antecedent peak flow event and H p2 , t p2 , t d2 , and t lag describe the subsequent peak flow event. Note that for all the explored scenarios, antecedent peak flow event commence at t = 0. For convenience, from here t lag is used to define the lag between two peak flow events and the time at which subsequent event commences, t d2 is used to describe the time at which duration of the subsequent event ends (i.e., at t lag + t d2 ), and t p2 for the time at which peak for the subsequent event occurs (i.e., at t lag + t p2 ).
With the estimated time-varying river stage (H sup (t)), we describe the pressure distribution at the SWI (∂Ω SWI ). For simplicity, we use an expression for prescribed head distribution that assumes a linear combination of head fluctuations induced by large-and small-scale bed topography (Stonedahl et al., 2010;Wörman et al., 2006): where S is channel slope, H sup (t) [L] is the time-varying river stage, Z SWI (x) is the function describing the bed topography, and h d (t) is the intensity of the dynamic head fluctuations (Elliott & Brooks, 1997), where the mean velocity is estimated with the Chezy equation for a rectangular channel as U s (t) = M −1 H sup (t) 2/3 S 1/2 with M as the Manning coefficient [L −1/3 T] (Dingman, 2009). Notice that the pressure distribution at the sediment-water interface is a function of both space x and time t, where the temporal fluctuations are induced by the peak flow events (see Equation 3).

Flow Model
Flow within the domain is driven by pressure gradients at the sediment-water interface (∂Ω SWI ). Neglecting the storage term, a reasonable assumption for submerged channel sediments, flow within the domain is described by the following version of the groundwater flow equation and Darcy's law Assuming that bedforms repeat periodically along the channel, we implemented a periodic boundary condition for the lateral boundaries (∂Ω u and ∂Ω d ; ). Under neutral groundwater conditions (i.e., without gaining and losing groundwater conditions), the only groundwater flow constraining the hyporheic zone is the ambient groundwater flow driven by the channel gradient (i.e., horizontal underflow component), and therefore no flow is assumed at the lower model boundary (∂Ω b ). The depth of this boundary (d b ) was selected to 10.1029/2020WR027113 Water Resources Research minimize boundary effects. Finally, the solution under steady state (i.e., baseflow conditions) is used as the initial condition for the transient simulations (i.e., during the peak flow event). This method of calculating pressure distribution at the SWI reproduces reasonable observations. It also allows the exploration of large number of scenarios with fewer complexities when implemented in the model, and with reduced computational demand.

Solute Transport Model and Delineation of the Hyporheic Zone
The advection-dispersion equation (ADE) is used to model the transport of conservative solutes within the streambed sediments , and D = {D ij } is the dispersion-diffusion tensor defined as (Bear, 1972) with α T (0.05 m) and α L d(α T /10) the transverse and longitudinal dispersivities [L] (Gomez-Velez & Harvey, 2014;Gomez-Velez et al., 2017), D m the effective molecular self-diffusion coefficient, ξ m = θ −1/3 is the fluid tortuosity (defined here with the Millington and Quirk model, Millington & Quirk, 1961), and δ ij is the Kronecker delta function.
Modeling the transport of a conservative tracer allows us to explore the mixing and extent of the hyporheic zone. We assume that the concentration of the tracer in the stream water column is C s , and therefore a prescribed boundary condition C(x, t) = C s is used along the SWI's inflow areas (∂Ω IN ). Outflow areas (∂Ω OUT ) along the SWI are advective boundaries where n·(D ∇C) = 0. Lateral boundaries (∂Ω u and ∂Ω d ) are periodic boundaries C(x = −L, y) = C(x = 2L, y) and the bottom boundary (∂Ω b ) is a no-flow boundary n·(q C−D ∇C) = 0. An initial condition for the concentration field is obtained from a steady-state simulation of the transport model (Equation 6) under baseflow conditions (i.e., H s = H 0 ). In this case, the hyporheic zone is defined as the zone with at least 90% of the pore water originated from the stream (i.e., C ≥ 0.9C s ). This definition is similar to the one proposed by Triska et al. (1989), (Gomez-Velez et al., 2017). Throughout the manuscript, we refer to this definition as the biogeochemical definition of the hyporheic zone.

Residence Time Model
The hyporheic zone residence time describes the time that water and solutes are exposed to the stream sediment biogeochemical conditions. Here, we evaluate the impacts of transient flow, driven by successive peak flow event, on the first moment of the hyporheic zone's residence time distribution. To this end, we use the approach outlined in Gomez-Velez et al. (2012), Gomez-Velez and Wilson (2013), and Gomez-Velez et al. (2017) where the moment of the residence time distribution is described by an ADE of the form 10.1029/2020WR027113 (Bear, 1972), and a 1 (

Water Resources Research
Initial and boundary conditions are defined following the approach of Gomez-Velez and Wilson (2013) and Gomez-Velez et al. (2017). Similar to the conservative transport model, the initial distribution of the first moment of the residence time distribution (mean residence time) was estimated under steady baseflow conditions.

Numerical Breakthrough Curves and Scenarios
To model the breakthrough curves, first, the advection-dispersion equation is implemented on multiple conservative tracers. Multiple conservative tracers are deployed to potentially estimate relative contributions by not only antecedent and subsequent events individually but also when the events are superimposed. We define the time of the occurrence of the subsequent event (t lag ) by the equation where η = 1/4, 2/4, 3/4, 4/4, 5/4, and 6/4. Depending on the values of η describing the time lag between ). c 1 , c 2 , and c 3 are the conservative tracers made active depending on the occurrence of the subsequent event. Cases I (top) and II (middle) constitute the scenarios with superimposition or overlapping of two events whereas Case III (bottom) represents scenarios with no overlapping of the two peak flow events.

10.1029/2020WR027113
Water Resources Research both events, tracking of the mixing between the two events would need deployment of two or three conservative tracers. Following are the formulated cases (See Figure 2): To track the mixing for Case I, two tracers are injected, namely c 1 and c 2 . The tracers are described by the piece-wise function where c 1 = c s (where we assume c s is the concentration of the conservative tracer in the water column) when the antecedent event is active or the rise in stage is caused by the antecedent event, c 2 = c s when both the events are active, that is, the stage rise is due to combined effect of the two events and finally we use c 3 = c s for stage rise due to the subsequent event only (see Figure 2 for further clarification).
For the tracer entering the sediment-water interface, that is, ∂Ω IN , the total mass that entered the system is estimated as where ∂Ω IN is the inflow boundary at the sediment-water interface. The mass leaving the hyporheic zone is given by Finally, the cumulative mass (or recovery) gives the breakthrough curve (BTC) estimated by the following mathematical statement Figure 3. Depiction of dimensionless stage hydrographs. The antecedent peak flow event is of higher magnitude and is identical for all the scenarios. For (a), i.e., Category 1, the subsequent peak flow event is of shorter magnitude and occurs at t lag (H p 1 > H p 2 and t d1 > t d2 ). As shown, scenarios with η = 1/4, 2/4, and 3/4 the two events are superimposed. For the scenarios with η = 4/4, 5/4, and 6/4, the subsequent event occurs after the duration of the antecedent event. For (b), that is, Category 2, H p 1 ¼ H p 2 and t d1 > t d2 and for (c), that is, Category 3, H p 1 > H p 2 and t d1 ¼ t d2 .

Water Resources Research
All the scenarios used for the simulations have the same duration (t d1 ) and magnitude of the antecedent event (H p 1 ; see Figure 3). Time to peak of the antecedent event t p 1 ¼ Sk 1 × t d1 and time to peak for the subsequent event t p 2 ¼ Sk 1 × t d2 , where Sk 1 = 1/8. Category 1 represents the scenarios with magnitude of the subsequent event (H p 2 ) less than the magnitude of the antecedent event (H p 1 ) and duration of the subsequent event (t d2 ) is less than the duration of the antecedent event (t d1 ) (Figure 3a). Category 2 (Figure 3b) corresponds to scenarios where the magnitude of the two events is the same (H p 1 ¼ H p 2 ), and the duration of the subsequent event is less than the duration of the antecedent event ( t d2 <t d1 ). Finally, Category 3 ( Figure 3c) corresponds to scenarios where the magnitude of the subsequent event is lower than the one of the antecedent event (H p 2 < H p 1 ), and the duration of events is the same (t d2 ¼ t d1 ). (Figure 3c).

Hyporheic Flow Field and Extent
The peak flow event induced pressure distributions along the sediment-water interface vary dynamically. The temporal evolution of flow fields are illustrated in Figure 4. Time t ≤ 0 represents the base flow conditions and time t ≤ t d2 represents the dynamics during the superimposed event for the scenario η = 3/4. The hyporheic zone dynamically expands and contracts during these events, with maximum expansion at t ¼ t p 1 due to the large magnitude of the antecedent event. Also, the hyporheic zone expands at t ¼ t p 2 ; however, the magnitude of the expansion is relatively less than that compared to t ¼ t p 1 due to the smaller magnitude of the subsequent event in the illustrated example. The magnitude and duration of the events determine the dynamic changes in the flow field.
We also use the biogeochemical definition to define and illustrate the hyporheic zone (see the black line in Figure 4). This is numerically achieved by introduction of a conservative tracer and tracking its transportation in the streambed. The contour line delineating the biogeochemically defined hyporheic zone with at least 90% of the water originating from the stream indicates that both hyporheic zone shape and size vary considerably for different times. However, these patterns tend to be more stable than compared to the variations

Water Resources Research
in the flow field. Such behavior can be explained by the flow field instantaneously propagating with pressure fluctuations at the sediment-water interface. On the other hand, the hyporheic zone extent based on biogeochemical definition takes into consideration transport and retention processes within the hyporheic zone. Observations demonstrate the variability in the hyporheic zone extent showing considerable change for t ¼ t p 1 , followed by t ¼ t p 2 , both linked to maximum rise in stage.
Another characteristic behavior demonstrated by the model results is the horizontal and longitudinal movement of stagnation zones (in each subplot of Figure 4) which would also be observed in the case of single peak flow event . These are the regions or locations where the absolute Darcy velocity is zero or near-zero to form stagnant zones. The presence of stagnation zones is dependent on the existence of the counter-directional local flow systems (Jiang et al., 2011). Dynamic stagnation zones are usually characterized as reactive hot spots, playing a primary role in types and rates biogeochemical transformations, such as the reduction of nitrates and nitrites in the nitrogen cycle Pinay et al., 2015). Furthermore, the dynamic variations in the flow field and the stagnation points bears relevance for fate and transport of contaminants in river-corridors.

Flux-Weighted Mean Residence Time
A representative value of residence times for the hyporheic exchange process is estimated by flux-weighting the modeled mean residence time along the sections of the SWI discharging hyporheic water into the stream (described in the residence time model). The mean residence time, μ corresponds to the first central moment, a 1 . Here we estimate flux-weighted mean residence time, μ * τ relative to baseflow conditions. Flux-weighted mean residence time ( μ * τ ) slightly increases at first with the antecedent peak flow event ( Figures 5, 6 and 7). This behavior is caused by the first high flush in all scenarios expanding the hyporheic zones, penetrating deeper into the sediments and discharging older hyporheic waters, reaching its maximum before an event's time to peak. Increased mean residence time of water discharged from the hyporheic zones for a short period of time is followed by progressive discharge of younger hyporheic waters due to activation of short flow paths closer to the sediment-water interface. This is in line with findings of Singh et al. (2019) and Gomez-Velez et al. (2017).
The antecedent peak flow event is followed by the low-magnitude (Figures 5 and 7) and high-magnitude subsequent peak flow event ( Figure 6) which is superimposed for the cases with η = 1/4, 2/4, and 3/4 (first row of

Water Resources Research
Figure 5 and 6) and separated for η = 4/4, 5/4, and 6/4 (second row of Figure 5). This event disturbs the system causing variability, preventing the system to attain base flow conditions for a certain period of time depending on the values on η, magnitude, and duration of the subsequent peak flow event. This leads to discharge of relatively older water (compared to reference conditions; see gray lines in Figures 5, 6, and 7) for a short duration which eventually starts declining shortly after the second peak is attained.
Variations in different scenarios is observed for different values of η. The results indicate that consecutively increasing the separation between the time-to-peak of the two events, t lag , exhibits higher variability of μ * τ Figure 6. Temporal evolution of flux-weighted residence time μ * τ of hyporheic waters leaving the sediment-water interface, as a function of dimensionless time. Category 2 scenarios are shown where H p 1 ¼ H p 2 and t d1 > t d2 . Vertical red lines and black lines correspond to time-to-peak and duration of the events, respectively. Different colors represents different η values. Black horizontal dotted line represents the baseflow conditions where μ * τ ¼ 1. Black solid curve correspond to the shape of the hydrograph and gray line to the μ * τ for the reference scenario. Figure 7. Temporal evolution of flux-weighted residence time μ * τ of hyporheic waters leaving the sediment-water interface, as a function of dimensionless time. Category 3 scenarios are shown where H p 1 > H p 2 and t d1 ¼ t d2 . Vertical red lines and black lines correspond to time-to-peak and duration of the events, respectively. Different colors represents different η values. Black horizontal dotted line represents the baseflow conditions where μ * τ ¼ 1. The black solid curve corresponds to the shape of the hydrograph and gray line to the μ * τ for the reference scenario.

Water Resources Research
(increase in the difference between the troughs in mean residence time, that is, the variability in mean residence time from one event to the next). The hyporheic waters discharged from the sediment-water interface during successive peak flow events is older, compared to the μ * τ for reference scenario (i.e., when the subsequent peak flow event does not exist) for the scenarios with η = 3/4, 4/4, 5/4, and 6/4. On the other hand, two events with a short interval of time in between causes lesser variability primarily because of activation of shorter subsurface flow paths with high flow velocities, leading to mixing of younger water in the hyporheic zone and decreasing μ * τ . Increased duration of the subsequent event deviates the system from the reference scenario for a longer period of time (Figures 5 and  7). Moreover, increased magnitude of the subsequent peak flow event (i.e., when H p 1 ¼ H p 2 ; Figure 6) leads to discharge of younger hyporheic water into the channel (i.e., water with lower values of μ * τ ) after the time-to-peak of the event, when compared to the scenarios with H p 1 > H p 2 (Figures 5  and 7).
Further analysis of maxima of the first central moment of flux-weighted mean residence time for the second peak, μ max shows interesting results ( Figure 8). When two peak flow events occur with shorter separation time (t lag ), mean age of the hyporheic waters leaving the streambed is relatively younger. Consecutively, increasing the value of t lag , μ max increases. When the η is 5/4 and 6/4, the maxima of central moment (μ max ) slowly reaches the plateau where for any further increase μ max become constant. This is an indication that the two events become independent with no memory effects from the antecedent peak flow event that is carried forward to the subsequent event and hence any response of the system is solely due to the subsequent event. Furthermore, subsequent peak flow event with larger magnitude increased the peaked mean RT (μ max ) (red curve in Figure 8). Also note that for η = 1/4, scenarios associated with Values are plotted as a function of time, that is, the time when second (i.e., subsequent) peak flow event commences (t lag ). Colors and symbols correspond to different values of η. Gray, red, and black curves represent Categories 1, 2, and 3 scenarios, respectively.

Figure 9.
Modeled breakthrough curve as a function of dimensionless time. Category 1 scenarios are shown where H p 1 > H p 2 and t d1 > t d2 . Solid colored lines correspond to c 1 (i.e., when the first event is active), dashed colored lines for c 2 (i.e., when both the events are active at the same time), and dotted lines is for c 3 (i.e., when only the second event is active). Vertical red lines and black lines correspond to time-to-peak and duration of the events, respectively. Gray line represents breakthrough curves for the reference scenario. Different colors represent different η values.

Water Resources Research
Categories 1 and 3 do not exhibit a μ max (see black and gray curves in Figure 8; instead, there is a slight deviations from attaining reference scenario conditions.

Numerical Breakthrough Curves
Conservative solute tracer retention was determined by estimating and comparing the integrated mass fluxes for different case-defined concentrations tracers. Mass entering and leaving the system was computed by multiplying the discharge by the concentration of conservative solute tracer entering the streambed and leaving the hyporheic zones (i.e., C ≥ 0.9C s ), respectively. Breakthrough curves of tracers were determined by comparing the integrated mass fluxes leaving the system to the total mass injected.
Observations from the numerical tracer tests and its response to successive peak flow events provided insights into the complex retention and transport behavior of conservative solute tracers within a sediment domain. For the scenario with η = 1/4 in Figure 9, the breakthrough curve for c 1 is steeper, due to higher mass transfer rate owing to shorter lag between the two events. As the concentration of c 2 emerges, the c 1 is pressed out of the hyporheic zones at a faster rate. The c 2 tracer corresponds to the concentrations which are the combination or mixing of antecedent and the subsequent peak flow events which exhibit broader breakthrough curve for η = 1/4. Another spike in c 1 indicates activation of the antecedent peak flow event after cessation of the subsequent peak flow event. Longest mass transfer, that is, flatter breakthrough curve for the c 1 is observed for the scenarios with η = 4/4, 5/4, and 6/4, primarily due to noninterference and nonsuperimposition of subsequent peak flow event. Similar behavior is observed for the subsequent peak flow event, that is, for c 3 . Shape of the breakthrough curve is dependent on the flow rate and diffusion characteristics. When the lag between the two events is shorter, the solutes are relatively less retained in the system and transported at a faster rate, potentially from the shorter flow paths. In addition, it is also observed that as t lag increases, the antecedent event's behavior becomes similar to the reference concentration which represents the scenario with a single peak flow event, that is, c 1 for η = 6/4 nearly coincided with reference concentration (black solid curve). Furthermore, increase in the magnitude of the subsequent peak flow event causes the concentration to increase at a faster rate (see Figure 10). Category 2 scenarios are shown where H p 1 ¼ H p 2 and t d1 > t d2 . Solid colored lines correspond to c 1 (i.e., when the first event is active), dashed colored lines for c 2 (i.e when both the events are active at the same time), and dotted lines is for c 3 (i.e., when only the second event is active). Vertical red lines and black lines correspond to time-to-peak and duration of the events, respectively. Gray line represents breakthrough curves for the reference scenario. Different colors represent different η values.

Water Resources Research
A prominent observation in breakthrough curves for scenarios with longer duration of the subsequent event ( Figure 11) is that the curves are much flatter when compared to scenarios with shorter duration. This can be explained by the slower mass transfer due to the longer spread of the event. Moreover, for η = 1/4, c 3 tracer is also injected as the longer duration of the subsequent event causes partial over lapping of the events and not complete as for the same condition in Category 1 scenario (Figure 9).

Solute Transport Retention and Dynamics During the Successive Peak Flow Events
Observations of the flux-weighted mean age and breakthrough curves unraveled interesting results detailing the impact of successive peak flow events on water and solutes' residence times and transport. The increased transient pressure gradient due to the peak flow events induces hyporheic flows that are primarily an instantaneous response to changes in the river stage. This instantaneous increase enlarges the size of hyporheic zones, elongates the subsurface flow paths, and increases hyporheic fluxes. Sudden rises in stage induce newly activated deeper and elongated subsurface flow paths that increase the transport time scale for the water and solutes. The results indicated that the duration of the impact in the residence times caused by the antecedent event is highly determined by the flow characteristics of the subsequent peak flow event. For example, if the subsequent peak flow event with a shorter magnitude occurs immediately after the antecedent event, there is only minor deviation in the change in residence times. The intensity, time of occurrence, and duration of the subsequent peak flow event determine the magnitude of variations in the residence times over the course of events. Primarily, when the subsequent peak flow event is of a large magnitude, the changes can be considerable and can exhibit long-term change in the system. As also illustrated by Harvey et al. (2012), discharge with longer duration and larger magnitude will potentially completely flush deeper and longer flow paths. It can furthermore be useful to consider the modulating role of groundwater fluxes in tandem with the effects of sequential peak flow events. The role of groundwater upwelling and downwelling has been documented in the study by Wu et al. (2018) for single peak flow events, where the authors found that the regional groundwater flow and geomorphological settings greatly modulate the temporal evolution of bedform-induced hyporheic responses driven by a single peak flow event. We expect that groundwater fluxes will have a similar modulating effect when a sequence of flood events is considered. Figure 11. Modeled breakthrough curve as a function of dimensionless time. Category 3 scenarios are shown where H p 1 > H p 2 and t d1 ¼ t d2 . Solid colored lines correspond to c 1 (i.e., when the first event is active), dashed colored lines for c 2 (i.e., when both the events are active at the same time), and dotted lines is for c 3 (i.e., when only the second event is active). Vertical red lines and black lines correspond to time-to-peak and duration of the events, respectively. Gray line represents breakthrough curves for the reference scenario. Different colors represent different η values.

Water Resources Research
For example, the subsequent event is likely to counteract the modulation of groundwater fluxes, changing the gradients adjacent to the streams, and hence reducing residence times and changing the patterns of solute transport. Future work will address this aspect in more detail.
Solute tracers results demonstrated that the peak flow events can lead to prolonged storage and delayed release of solute tracers, for the scenarios where the two peak flow events are not superimposed and occur for a longer duration of time. This is mainly because of the slower mass transfer due to the longer spread of the event. This is also indicated in the study by Harvey et al. (2012). Our results indicate that time lag between two peak flow events (t lag ) is a crucial determinant of water and solutes retention time in the hyporheic zones. It is observed that if the time of the occurrence of the subsequent event is longer than the memory of the system from the antecedent event, the response of the successive events can be treated independently. However, other geomorphological parameters such as the channel gradient, bedform amplitude, wavelength (Elliott & Brooks, 1997;Singh et al., 2019;Wu et al., 2018) and flood characteristics like intensity, duration, and skewness of the event also determine the transport and retention behavior of water and solutes.

Biogeochemical and Ecological Implications
Herein we use numerical tracer tests to investigate the impacts of successive peak flow events on water and solute retention and transport in sediment domain. This research mainly focuses on variations in hyporheic flux velocities, residence times, and transport of conservative solutes under dynamic flow conditions.Our results indicate that high-flow events caused an expansion of the hyporheic zones and increase hyporheic fluxes, indicating that events have the potential to enhance downwelling of surface water rich in oxygen, dissolved organic matter, and other nutrients, being delivered into the hyporheic zone at high concentrations into greater depths and larger streambed areas (Gomez-Velez et al., 2017;Hester & Doyle, 2008;Singh et al., 2019). Findings by Drummond et al. (2017) suggested that storm events can mobilize or remobilize solutes and fine particles within the streambed and therefore further fuel biogeochemical transformations. That said, successive events occurring within a short duration of time, that is, when system is unable to recover from antecedent event may not enhance stream metabolism, owing to the larger amount of water pressed from shorter subsurface flow paths.
Sequential events have the potential to induce pockets with distinct chemical composition and steep reactive gradients within the hyporheic zone. For example, surface water solutes from the antecedent peak flow event can be pushed at greater depths when the subsequent event commences. Then, a large proportion of the dissolved solutes from the subsequent event are flushed first through shorter and shallow subsurface flow paths while the solutes from the antecedent event remain within the streambed sediments. This mechanism can lead to biogeochemical layering and significantly affect the transport of dissolved organic carbon, dissolved oxygen, contaminants, and microbes.
During successive peak flow events, the altered hydraulic conditions affect the solute transport and retention. Increased residence times with longer separation between the two events suggests that high discharge events when occurred after a certain period of time has the potential to positively contribute to processes such as nitrification and denitrification. These results suggest that depending on the magnitude of the events, increased hyporheic flow paths due to repeated peak flow events may enhance the stream metabolism. The results presented by Trauth and Fleckenstein (2017) demonstrated that a single discharge event increased the reactive efficiency of the hyporheic zone. Our findings highlight the effects of consecutive flow perturbations in the hydrodynamics and transport characteristics of the hyporheic exchange process. The compounding nature of these perturbations can have significant implications for the reactive efficiency of the hyporheic zone, and therefore, it should be considered in rivers exposed to natural and artificial (e.g., dam operations and wastewater discharge plants) fluctuation of flow As the biogeochemical time scales for oxygen consumption are important determinant for shifting the respiration from aerobic to anaerobic conditions in the streambed, elongated subsurface flow paths can promote biogeochemical transformations in deeper parts of the streambed (Gu et al., 2008;Hinton et al., 1997;Inamdar et al., 2004;Kennedy et al., 2009;Krause et al., 2013;Malcolm et al., 2004). In addition, transience variability in hydraulic conditions and stream chemistry can have cascading influence on diversity and productivity of hyporheic organisms (Bruno et al., 2009). 10.1029/2020WR027113

Conclusions
Research results of this study highlight the relevance of considering the impacts of successive peak flow events on hyporheic exchange flow, transport, and retention times of water and conservative solutes within the streambed. Systematic hydrograph scenarios with variations in time lag between the two flow events, magnitude and duration of the events were simulated with a two-dimensional model for sinusoidal bedforms and analyzed for flow, transport and residence time distributions. The time lags between simulated event scenarios ranged from complete overlapping, that is, superimposition of two events, partial overlapping to scenarios with no overlapping, that is, representing repeated peak flow events.
During peak flow events, instantaneous changes in the temporally variable hyporheic flow field were observed due to rapid pressure wave propagation within the streambed. However, the actual increase in hyporheic flow as defined based on the geochemical definition of hyporheic zone reacted generally less sudden. The results demonstrate that longer time lags between successive peak flow events cause higher temporal variability in the flux-weighted mean residence times, with a higher potential to discharge older hyporheic pore water compared to the reference scenario. Furthermore, with increasing time lag between two flow events, the maxima of the first central moment, that is, mean age slowly reaches the plateau where for any further increase, the value remains constant. Simulated breakthrough curves depicted generally reduced gradient for scenarios with longer time lags between events and enhanced duration of the subsequent event. Variability in event magnitude also had an impact on the shape of the breakthrough curves. Higher magnitude of the subsequent event caused the breakthrough curves to become steeper, when compared to low magnitude and longer duration of the subsequent event. It is observed that successive events can be treated independently if the occurrence of subsequent event is longer than the memory of the system from the antecedent event. The model-based analysis of simulated breakthrough curves of conservative tracer for successive peak flow events highlights the impact of the type and duration of peak flow successions on residence time distributions and solute transport that are relevant for biogeochemical turnover in hyporheic zones.

Data Availability Statement
All data required to reproduce the figures in this paper are available on the data repository of the University of Birmingham's library.