Interactive comment on “ Modelling hyporheic processes for regulated rivers under transient hydrological and hydrogeological conditions ” by D .

Understanding the effects of major hydrogeological controls on hyporheic exchange and bank storage is essential for river water management, groundwater abstraction, restoration and ecosystem sustainability. Analytical models cannot adequately represent complex settings with, for example, transient boundary conditions, varying geometry of surface water–groundwater interface, unsaturated and overland flow, etc. To understand the influence of parameters such as (1) sloping river banks, (2) varying hydraulic conductivity of the riverbed and (3) different river discharge wave scenarios on hyporheic exchange characteristics such as (a) bank storage, (b) return flows and (c) residence time, a 2-D hydrogeological conceptual model and, subsequently, an adequate numerical model were developed. The numerical model was calibrated against observations in the aquifer adjacent to the hydropower-regulated Lule River, northern Sweden, which has predominantly diurnal discharge fluctuations during summer and long-lasting discharge peaks during autumn and winter. Modelling results revealed that bank storage increased with river wave amplitude, wave duration and smaller slope of the river bank, while maximum exchange flux decreased with wave duration. When a homogeneous clogging layer covered the entire river–aquifer interface, hydraulic conductivity positively affected bank storage. The presence of a clogging layer with hydraulic conductivity −1 significantly reduced the exchange flows and virtually eliminated bank storage. The bank storage return/fill time ratio was positively related to wave amplitude and the hydraulic conductivity of the interface and negatively to wave duration and bank slope. Discharge oscillations with short duration and small amplitude decreased bank storage and, therefore, the hyporheic exchange, which has implications for solute fluxes, redox conditions and the potential of riverbeds as fish-spawning locations. Based on these results, river regulation strategies can be improved by considering the effect of certain wave event configurations on hyporheic exchange to ensure harmonious hydrogeochemical functioning of the river–aquifer interfaces and related ecosystems.


Introduction
Surface water-groundwater interfaces have recently received growing research interest (Sophocleous, 2002) and have become the focus of multiple water resources management policies (Klöve et al., 2011).The hyporheic zone that harbours river-aquifer interactions plays a key role in riverine and riparian ecosystem functioning (e.g.Krause et al., 2011).The reactive nature of this zone maintains exchange and transformation of solutes along the pathways between surface water and groundwater.The hyporheic zone is an appreciated habitat for hyporheos, microorganisms and bacteria occupying the space below and along the river channel (Boulton et al., 1998).Besides having negative impacts on the ecosystem of the zone itself, with changes in hyporheic water composition (Calles et al., 2007;Siergieiev et al., 2014c), alteration of the hyporheic functionality due to surface water-aquifer disconnection can also severely modify neighbouring ecosystems.As a result, riparian zones and adjacent wetlands may experience changes in groundwater table, water and nutrient fluxes.Furthermore, restricted hyporheic exchange limits mobilisation of solutes from the ri-Published by Copernicus Publications on behalf of the European Geosciences Union.

D. Siergieiev et al.: Modelling hyporheic processes for regulated rivers
parian zone and their fluxes into the river during key hydrological events as a result of river-aquifer disconnection (Burt and Pinay, 2005), which can in turn affect surface water quality (Valett et al., 1996).
A large number of rivers worldwide are obstructed by dams (Nilsson et al., 2005) and are therefore subject to artificial discharge fluctuations.These fluctuations stress hyporheic exchange flows, which often results in degradation of the river-aquifer continuum.A major impact is the alteration of river sediment transport and, subsequently, increased colmation of the riverbed (Brunke and Gonser, 1997;Blaschke et al., 2003), which deteriorates river-aquifer hydraulic connectivity (Burt and Pinay, 2005) and controls functional changes in the hyporheic zone (Siergieiev et al., 2014c).In addition, inundation of the river banks by construction of reservoirs changes the shape of the river-aquifer interface, successively affecting the hyporheic exchange (Doble et al., 2012a).
Therefore, it is necessary to understand the hydrogeological functioning of this interface under artificial conditions, such as hydropower-regulated rivers, in order to incorporate this knowledge into water resource management and thereby improve the functional behaviour of the hyporheic zone by optimising river discharge strategies (e.g.Hanrahan, 2008).
The hyporheic zone size, bank storage volume and bank fluxes vary with river stage fluctuations (amplitude, duration) and river bank conditions (slope, hydraulic conductivity).Todd (1955) provided a first theoretical analysis of floodinduced bank storage.This work served as the foundation for later analytical solutions and numerical simulations of floodplain hydrology.The dynamics of bank storage were later estimated using analytical models (e.g.Cooper and Rorabaugh, 1963) based on the following simplifications: single flood wave, homogeneous aquifer and a fully penetrating vertical river bank.During base flow conditions, hydrological riveraquifer interactions can be described by precipitation-runoff models (Butturini et al., 2002).However, these often neglect the distributed effects of e.g.unsaturated zone processes or topography, resulting in residual unexplained variability in bank storage.Chen and Chen (2003) used a numerical approach to simulate the bank storage response to changes in river stage and riverbed hydraulic conductivity for a partially penetrating river with vertical banks in a fully saturated aquifer.They pointed out the importance of subsurface anisotropy, which governs the directions of hyporheic zone development.Simultaneous consideration of seepage and a variably saturated aquifer showed that unsaturated zone processes have a pronounced effect on bank storage (Li et al., 2008;Doble et al., 2012a).Inclusion of the unsaturated zone in bank storage simulations decreased the modelled storage and improved the return flows (Doble et al., 2012b).Furthermore, models that consider vertical river banks for sloping banks under-estimate bank storage (Doble et al., 2012a).Understanding of bank storage processes can improve hyporheic ecotone and river-aquifer continuum concepts, with positive implications for e.g.base flow separation techniques (Mc-Callum et al., 2010) and ecosystem sustainability (Schneider et al., 2011).
For several seasons, hyporheic exchange was studied in the hydropower regulated Lule River, northern Sweden (Siergieiev et al., 2014a, c).Low hydraulic conductivity of the riverbed and daily varying river discharge have resulted in depleted hyporheic exchange flows across the river-aquifer interface (Siergieiev et al., 2014a).Deteriorated river water quality as a result of regulation (Smedberg et al., 2009;Siergieiev et al., 2014b) may partly depend on suppression of hyporheic processes due to regulation (Valett et al., 1996).Improved understanding of the major hydrogeological controls of hyporheic exchange has legacy effects on understanding geochemical fluxes between surface water and groundwater (Fritz and Arntzen, 2007) and can provide a platform for implementation of environmental flows and improved management and ecological status of regulated rivers.
The aim of this study was therefore to provide a set of scenarios with variable river discharge schemes (wave duration and amplitude), river bank slope and riverbed hydraulic conductivity, in order to investigate the effects of these parameters on fluxes across the river-aquifer interface, bank storage volume, fill/return time ratio and residence time.A realistic case study was used to set up a conceptual and a numerical model and justify the use of the method for further scenario simulations.The major difference of this work to earlier studies is an attempt to demonstrate possible limits of the surface water-groundwater models and the data necessity.Furthermore, application of numerical surface water-groundwater exchange modelling in a highly dynamic environment, common for many other regulated rivers in the world, was evaluated.This kind of modelling can be used in different applications such as river restoration, implementation of environmental flows, integration with regional models etc.The topic is also of interest from the perspective of improvement of ecological status of impacted watercourses as requested by Water Framework Directive.Understanding and coupling of river-aquifer functioning to the transfer of nutrients across their interface and the impact on sediment transport and biological migration are essential for improved ecological status of regulated rivers, which are today severely impacted by disruption of hydrological connectivity and related biogeochemical consequences.

Site description
The measurement profile orthogonal to the river included an observation station in the river and two groundwater wells (Fig. 1) with hourly registration of water level during 2010-2011 and every 15 min during 2012.The riverbed at the site slopes gently towards the middle of the channel and is composed of silty-clayey material with layers of highly conductive stratum and laterally spread sparse patches of sand and gravel.The former floodplain is limited by an earth embankment from the river side and stretches for over 150 m where it meets the foot of the hillslope.The depth to the bedrock varies from none to 60 m with an average depth of around 10 m and thinnest overburden towards the hillslope.The flat topography of the floodplain and the hillslope orientation which is normal to the river suggest orthogonal groundwater flow towards the channel.Orientation of the groundwater well profile was chosen accordingly.
3 Materials and methods

Data collection
The soil was visually inspected during installation of groundwater wells.Samples were collected at 0.3 m interval and sieving analysis was performed on three samples from each of the two locations.Unsaturated flow parameters were estimated on these six selected samples using pressure pot experiments to obtain water-holding characteristics (Ehlert, 2014).
To assess saturated hydraulic conductivity, repeatable slug tests (three in each well) using both falling and rising hydraulic head were carried out in the wells, while a direct push piezometer (two repeatable tests at two locations; 1 and 2.5 m from the shoreline) using a falling head was performed at the riverbed at 0.1 m depth interval down to 0.7 m depth (Siergieiev et al., 2014a).A harmonic mean over the top 0.3 m least permeable layers was used as hydraulic conductivity of the clogging layer in the model.

Conceptual model
The data collected at the site did not allow development of a highly distributed model.Parameter values obtained in the field and in the laboratory were therefore averaged for the saturated and unsaturated zone and the clogging layer.The vertical 2-D conceptual model considered a homogeneous aquifer, partially penetrating river, sloping banks, unsaturated zone and clogging layer that covered the entire riveraquifer interface (Fig. 2).The Dirichlet boundary condition at the riverside was varied according to the measured water level time series.The top of the model was represented by a constant flux boundary to consider recharge, which was assumed to be 50 % of annual precipitation (Lemmelä, 1990) of 470 mm yr −1 , resulting in 6.5 × 10 −4 m d −1 recharge.No flow boundaries were assigned to the remaining borders.The distance to the right boundary was set to ensure that influences of river stage fluctuations did not reach this boundary for the longest fluctuation period.A first approximation of the distance of influence based on an analytical solution (Sawyer et al., 2009) assumed negligible riverbed resistance and is therefore questionable in the present case (e.g.Singh, 2004).This estimate of the maximum extent (180 m for a 1-month fluctuation period) was further tested using the numerical model.
The following assumptions were used in the model: -Two-dimensional model space  -The same unsaturated parameters for the entire model domain -Surface flow resistance and hydraulic effects of the river processes due to variable discharge were neglected -Viscosity effects (temperature and solute concentration differences between the river and the aquifer) were neglected.

Numerical model
The numerical modelling code FEFLOW 6.2 (Diersch, 2014) was used to simulate variably saturated flow during riveraquifer interaction by solving Richards' equation using the PARDISO solver (Schenk and Gärtner, 2004).The model domain was discretised using the triangle mesh generator.The time step was set to 30 min intervals to keep the computational time within reasonable limits and was increased to 1 hour during the validation run.To enable inverse parameter estimation for the van Genuchten model, a plug-in for coupling FePEST (graphical user interface for PEST by Doherty et al., 2011) with FEFLOW was developed (Ehlert, 2014).

Model calibration
The model was sequentially calibrated against measurements in L5 and L25 collected during June-October 2012.First, only hydraulic conductivity, specific storage and porosity were calibrated, followed by unsaturated van Genuchten parameters and maximum and residual saturation.Finally, all parameters were calibrated together.To track improvement of the fit between modelled and observed hydraulic head, the regression coefficient (R 2 ), Nash-Sutcliffe index (NS) and root mean square error (RMSE) were calculated for each calibration run.
The conceptual understanding of hydrogeological processes is often erroneous in terms of boundary and initial conditions (Bredehoeft, 2005).While the initial state of models is often calibrated, the importance of other conceptualisation aspects seems to be rarely verified, even though their ef- fects have been debated by different authors (e.g.Refsgaard et al., 2006).To test our assumptions on hydrogeological conditions in the area, a sensitivity analysis for model boundary conditions was carried out.The distance from the river to the aquifer boundary was varied from the reference distance (200 m) to 50, 500, 1000 and 5000 m, keeping the recharge at 50 % of annual precipitation.Afterwards, the recharge rate was changed from the assumed 50 to 30 and 70 %, keeping the distance from the river to the aquifer boundary at 200 m.

Modelling scenarios
Based on the calibrated model domain, the effect of multiple hydrogeological parameters on hyporheic exchange was evaluated, varying one parameter at a time.Artificial river stage variations were applied according to the distribution of commonly observed amplitudes and durations during 2012 (Fig. 3).The head boundary on the river side was varied as a cosine-shaped wave between t = 0 and t = t , with ampli- where h is the hydraulic head (m), t is time (h), t is the duration of the stage oscillation (h), h 0 is the head at t = 0, and h max is the maximum head (at t = t /2).All scenarios used a single wave event and were terminated after steadystate conditions were reached.
The sensitivity of the model to various scenarios was evaluated using the flux across the river-aquifer interface, bank storage and the ratio between the time to fill and to empty the bank storage.A reference simulation was based on the calibrated model (hydraulic conductivity of the aquifer 2.14 m d −1 and of the clogging layer 0.01 m d −1 ), the actual conditions at the site (bank slope 10 • ) and an input wave with 3 h duration and 0.4 m amplitude.Other scenarios included varying hydraulic conductivity of the clogging layer (0.001, 0.01, 0.1, 1 m d −1 ), river bank slope (5, 10, 15, 30, 45 • ), wave amplitude (0.03, 0.1, 0.4, 0.7, 1.0 m) and duration (3, 6, 12, 24, 168 h) (Fig. 4).The initial conditions were generated by running a transient simulation with constant hydraulic head at the riverside, no flow at the aquifer side and constant distributed diffuse recharge rate at the top for the time sufficient to recreate steady-state conditions.The flux into the aquifer was considered to be positive.The exchange flux per metre riverbed width (m 2 d −1 ) was identified as the difference between inflow and outflow rates across the river boundary.The bank storage (m 2 ) was the cumulative exchange flux multiplied by the time step assuming the flux being either positive or negative for in-and outflow from the model, respectively, while bank storage always positive.The fill time was the time required for the bank storage to reach its maximum, whereas the return time was the time between the maximum and zero bank storage on the falling limb.Residence time was calculated as the sum of the fill and return times.and the river (below) compared with the simulated results using calibrated parameters (see Table 1).

Model calibration
Sequential calibration using FePEST yielded minor underestimation of hydraulic head in the beginning of the simulation and minor over-estimation in later parts (Fig. 5).Generally, the fit was better for the well closer to the river, i.e.L5.The resulting R 2 , NS, RMSE were 0.95, 0.89, 0.06 m, respectively, for observation well L5 and 0.89, 0.77, 0.08 m, respectively, for well L25.The calibrated model (Table 1) was validated using observation data from 2010 with R 2 , NS, RMSE of 0.85, 0.97, 0.09 m, respectively, for well L5 and 0.74, −17.85, 0.43 m, respectively, for well L25 (data not shown).
Numerical simulation results verified that an aquifer boundary located 200 m away from the river was out of reach of influences induced by the river stage fluctuations used here.Sensitivity analysis of the boundary conditions resulted in substantial over-estimation of measured hydraulic head for a model domain size of 5000 m.The fit improved slightly for a domain size of 500 and 1000 m compared with the 200 m long domain.However, the larger model domains were discarded to keep the computation time reasonable.A model domain of 50 m tended to under-estimate the measured hydraulic head for both L5 and L25.The highest and lowest recharge rates over-and under-estimated the observed data, respectively.Therefore, the recharge rate taken as 50 % of precipitation was the most suitable solution for the present case.Overall, observation well L25 was more affected by changes in the recharge than well L5, indicating a strong influence of groundwater gradient on L25 and of the river boundary on L5.However, it was recognised here that the final influence of recharge and distance to the boundary on calibration results is a combined effect rather than a single effect of one of these.

Scenarios
The simulated scenarios of varying river bank slope, clogging layer hydraulic conductivity and input wave amplitude and duration were compared based on their effect on the resulting exchange fluxes, bank storage and residence time.

Exchange fluxes
There were variations in exchange fluxes across the riveraquifer interface as a result of varying forcing parameters (Fig. 6).The maximum exchange flux decreased with river bank slope (Fig. 6a).A change in the bank slope from 5 to 10 • caused a similar decrease in the exchange flux as a change from 10 to 45 • .The exchange flux increased with hydraulic conductivity of the interface (Fig. 6b).However, for the scenario with K = 0.001 m d −1 , the fluxes were always directed towards the river, indicating no bank storage and thus no hyporheic exchange.The increase in exchange flux was logarithmically proportional to the change in hydraulic conductivity, e.g. a 100-fold increase in hydraulic conductivity generated only a 7-fold rise in exchange flux.
The wave amplitude was related positively to the exchange flux across the river-aquifer interface (Fig. 6c).There were no positive (towards the aquifer) fluxes for amplitude 0.1 m and lower.Every further 0.3 m increase in amplitude resulted in a 0.4 m 2 d −1 increase in the maximum simulated fluxes.The maximum exchange flux was higher for shorter wave duration times, e.g.0.40 m 2 d −1 for a 3 h wave and 0.09 m 2 d −1 for a 168 h wave (Fig. 6d).

Bank storage
The effects of bank slope, hydraulic conductivity of the clogging layer and input wave amplitude and duration on bank storage volume were plotted as a function of time (Fig. 7).Bank storage increased with lower bank slope (Fig. 7a).For example, an almost 5-fold increase in bank slope from approx.10 to 45 • reduced bank storage by less than 50 %.Meanwhile, a decrease in bank slope from 10 to 5 slope for river-aquifer exchange.Overall, the bank storage for 5 • was 5-fold higher than that for 45 • .A 10-fold increase in hydraulic conductivity of the riverbed from the reference scenario (bank storage = 0.02 m 2 ) improved bank storage by 500 % (0.10 m 2 ), which further increased to 0.16 m 2 with another 10-fold increase in hydraulic conductivity.As was the case for exchange flux, virtually no bank storage occurred for the scenarios with 0.03 and 0.1 m wave amplitude (Fig. 7c), which were the most common amplitudes at the observation site (Fig. 3).An approximately 60 % rise in wave amplitude (from 0.4 to 0.7 m) resulted in a 30 % increase in maximum bank storage and a 50 % increase in maximum exchange flux.Duration of the river stage oscillation also positively affected bank storage (Fig. 7d).A 56-fold rise in duration (from 3 to 168 h) resulted in a 7-fold increase in maximum bank storage (from 0.02 to 0.14 m 2 ).

Residence time
The timing of bank storage (residence time and return/fill ratio) was examined under different modelling scenarios (Fig. 8).The residence time and return/fill time ratio decreased with increasing bank slope (Fig. 8a).The return time always exceeded the fill time except for the slopes above 30 • , which indicated t R /t F = 1.Increased hydraulic conductivity of the river-aquifer interface increased the return time, which positively affected the overall residence time of river water in the subsurface (Fig. 8b).Nonetheless, the residence time was highest (6.5 h) for the scenario with hydraulic conductivity of the interface of 0.1 m d −1 , marginally exceeding the scenario with the highest hydraulic conductivity (6 h).Return time increased with rising wave amplitude, as did the residence time, ranging from 0 h for the smallest wave to 9.6 h for the largest (Fig. 8c).The ratio between the change in res- idence time and the change in amplitude varied between 0.3 and 0.5 and was higher at lower amplitudes.The return time of bank storage was longer than the fill time for the waves with duration below 24 h and decreased with wave duration (Fig. 8d).The return/fill time ratio decreased by almost twothirds from the shortest wave duration (1.7) to the longest (0.6), whereas the residence time increased by more than one order of magnitude from the shortest wave to the longest (from 4.8 to 96 h).

Conceptual and numerical models and calibration
The model calibration resulted in hydraulic conductivity one order of magnitude higher than estimated via field tests.This difference is attributable to the limitations of slug tests, which provide point data for saturated hydraulic conductivity around the well filter.According to the observed soil profile (Fig. 1), grain size decreased with depth and can be a reason for lower measured and higher calibrated hydraulic conductivity of the aquifer.It was beyond the scope of this work to analyse whether the calibrated parameter set converged around a local or global optimum.However, there was a tendency for over-estimation of both saturated hydraulic conductivity and effective porosity (Table 1).Hydraulic conductivity and porosity are related through the hydraulic diffusivity term that controls the connectivity of a high-permeability flow path (Knudby and Carrera, 2006).Hydraulic diffusivity remains virtually the same with proportional change in both hydraulic conductivity and porosity.Therefore, even if

D. Siergieiev et al.: Modelling hyporheic processes for regulated rivers
over-estimation of both parameters by model calibration took place, this would have had a limited effect on the results.Although the validation results showed a good fit for hydraulic head at observation well L5, major under-prediction of hydraulic head at L25 was observed.This could be related to the substantially different precipitation patterns during the 2 years (2012 was used for calibration and 2010 for validation) and the fact that L5 is more influenced by the river and L25 by the aquifer.Recharge on top of the model, used for both calibration and validation runs, was assumed to be constant and equal to 50 % of mean annual precipitation.In the following sections, the implications of the modelling results for hyporheic exchange and the limitations due to the assumptions used are discussed.

General reflections
The implications of sloping banks for numerical modelling have been discussed previously, e.g. by Doble et al. (2012a).In terms of functioning of the river-aquifer interface, low bank slope angles increase the contact area between river and aquifer and therefore result in enlarged volume of bank storage and size of the hyporheic zone.This also positively affects the residence time, which is primarily governed by the penetration distance of river water and the return time necessary for it to discharge back into the river.Enhanced hyporheic exchange due to river bank slope would be true for many regulated rivers, as construction of reservoirs is associated with river floodplain inundation, and therefore with the formation of gently sloping banks along the channel.Consequently, lower bank slope angles have the potential to improve e.g.spawning conditions and species richness (Hanrahan, 2008).However, this appears not to be the case in several regulated temperate and boreal rivers due to the effect of colmation (Brunke and Gonser, 1997;Blaschke et al., 2003;Calles et al., 2007;Siergieiev et al., 2014a).In the present study, bank storage decreased with decreasing hydraulic conductivity of the river-aquifer interface (Fig. 7b).Therefore, in rivers with clear interstices, flat river banks contribute greatly to an increase in bank storage and hyporheic exchange, as demonstrated by the simulations, whereas this effect is hampered in rivers with a clogged riverbed.
River wave duration and amplitude were positively related to bank storage (Fig. 7c, d), but only amplitude positively affected exchange flux (Fig. 6c).As opposed to bank storage, exchange flux is more dependent on soil properties than on input wave configuration.This is supported by the fact that the peak exchange flux decreased with the prolonged wave duration (Fig. 6d), due to smaller hydraulic gradients at the river-aquifer interface, whereas the maximum bank storage increased (Fig. 7d).
A linear relationship between maximum bank storage and the product of wave amplitude and period has been reported previously by Todd (1955).However, this was only valid for a fully saturated homogeneous aquifer adjacent to a fully penetrating river.Using the results of the modelling scenarios in the present study, it was possible to show that there is a relationship between the ratio of wave duration/amplitude and bank storage or residence time for waves with amplitude exceeding 0.1 m (Fig. 9).This indicates that there is an optimal wave configuration (duration and amplitude) for every specific set of hydrogeological conditions that accounts for the highest bank storage and can potentially improve hyporheic exchange and minimise energy losses in hydropower-regulated rivers.
The hysteresis patterns observed for the t R /t F ratio for different modelling scenarios illustrate that the process of filling the pores of an aquifer is different from that of draining them and depends on hydraulic gradient and river-aquifer contact area (Fig. 8).The former is a function of the river wave configuration, while the latter depends on the river bank slope.The contribution of bank storage to river runoff is complex and of high importance in catchment hydrology (Harr, 1977;Turton et al., 1992;McGlynn et al., 2004).The results presented here suggest that with decreasing bank slope, the contribution of bank storage to the river extends in time, prolonging the falling limb of the river hydrograph.The same effect occurs with rising amplitude, which generates steeper hydraulic gradients across the river-aquifer interface.However, it requires less time to return the bank storage to the river with prolonged wave duration.A wave duration exceeding 24 h indicates faster return than the time required to fill the soil moisture deficit.These modelling results were obtained for a one-time wave event and no repeated wetting process was simulated.However, it is known that the hysteresis pattern can change direction over time (McGuire and McDonnell, 2010), which implies that the patterns observed here may differ for initially wet soil.Generally, a simple estimate of water residence time included in our scenarios is a useful proxy for hyporheic geochemical processing.
Investigation of theoretical scenarios creates a platform for development of more environmental friendly river regulation strategies as requested by Water Framework Directive.

Site-specific implications
Hyporheic exchange at the observation site was mainly characterised by hindered water flow across the river-aquifer interface and had residence time sufficient to establish suboxic conditions in the subsurface due to a number of reasons.The observed hydrograph for the Lule River was mainly dominated by short-term regulation with daily discharge peaks during July-early August (first 40 days of the simulation) and by long-term regulation with extended discharge waves during late August-October.A geochemical investigation of the hyporheic zone at the site revealed a basically suboxic environment, with elevated dissolved concentrations of Fe, Mn, NH 4 and organic carbon (Siergieiev et al., 2014c).These conditions suggest that the area along the banks and below the bed at the site experiences deficiencies in river water intrusion, which are primarily caused by the river discharge and the clogging layer.Using nitrogen as an example, the hyporheic zone is a nitrate source at low residence time and a nitrate sink at high residence time (Zarnetske et al., 2011).Consequently, biogeochemical activity in the hyporheic zone is controlled by exchange fluxes and bank storage (Gu et al., 2012).The combination of a rapidly rising discharge limb with long duration time favours intensive intrusion of river water into the subsurface, transfer of oxygen and dissolved organic carbon, and therefore promotes nitrification.Note that wave amplitude had a higher influence on the maximum flux across the river-aquifer interface, whereas wave duration affected total bank storage, i.e. the subsurface volume available for hyporheic exchange.This is explained by a steep hydraulic gradient across the river-aquifer interface and thus increased exchange flows due to a rapid rise in river discharge.To provide stable conditions for these ecologically important flows, an extended discharge wave with a sharp rising limb is required.At the Lule River study site, bank hyporheic exchange was triggered by 40 % of the wave events during 2012.These conditions had a potential to reset pore water geochemistry but only a part of it might satisfy sufficient residence time for reactions to occur and thus guarantee an effective biogeochemical exchange.The validity of this relationship requires further testing by e.g. a sediment transport survey, among other techniques, which can form the basis for implementation of environmental flows in restoration programmes (Schneider et al., 2011).
It is not only the hyporheic zone intimately connected to the river that can be affected by fluctuating river water stages, but also the distant groundwater.The simulation results indicated that groundwater head was affected by pressure propagation beyond observation well L25 (25 m distance to the river).This can have an impact on oxidation-reduction conditions in the aquifer due to changes in the redox potential during wetting and drying cycles of the soil (Reddy and Patrick, 1975;Cavanaugh et al., 2006).The relationship between the depth to the groundwater and groundwater composition in observation well L25, sampled during the pe-riod May-October 2011, was investigated.Based on nine water quality samples and principal component analysis, the first two significant components explained 77 % of the data variance, indicating a positive correlation between depth to groundwater and NO 3 concentration and a negative correlation between depth to groundwater and Fe and Al concentration.A positive correlation between Mn and alkalinity and depth to groundwater explained only 13 % of the data variance and P showed no relationship.For the significant correlations, a rising groundwater level promoted a more reduced environment and was associated with higher Fe and lower NO 3 concentrations.This suggests that transient changes in river water stages in response to hydropower management can force time-dependent alterations in groundwater quality, with further potential impacts on riparian soils.

Limitations
The assumptions made in this modelling study resulted in the following limitations: -Because of the vertical 2-D conceptualisation perpendicular to the river, longitudinal fluxes parallel to the river were neglected.A 3-D model is required for proper consideration of these processes.Bates et al. (2000) argued that the contribution of the longitudinal component is most important at the beginning and end of an event, implying confidence about timing but not about the absolute value of computed peak bank storage and fluxes using a 2-D approach.
-In aquifers with a clogging layer, hydraulic pressure propagation will always be ahead of water flow that follows oscillations at the river-aquifer interface (Welch et al., 2014).Assuming homogeneous subsurface media, solute travel time may exceed that of the pressure, resulting in over-estimated bank storage.In addition, it was assumed that all return flow came from bank storage, even though it contains a mixture of old water from the unsaturated zone and groundwater (Burt and Pinay, 2005).This is crucial for chemical fluxes through the hyporheic zone (McDonnell, 1990)  for chemical hydrograph separation (McCallum et al., 2010).Because the response of solute fluxes to bank storage is dependent on heterogeneity, verification of the fluxes obtained by pressure propagation using measurements of solute concentrations, e.g.electrical conductivity (Welch et al., 2014), or measurements of temperature (Anibas et al., 2012) may be required.
-Lateral variability in riverbed hydraulic conductivity at the site was simplified by implementing a continuous low hydraulic conductivity layer.In field settings, however, a riverbed with variable sediment composition is much more likely (Hancock and Boulton, 2005;Siergieiev et al., 2014a), which suggests that hyporheic exchange seeks more conductive patches.This assumption is likely to result in under-estimated hyporheic exchange (Kalbus et al., 2009) and partially compensate for using homogeneous media (see above).
-Hydraulic conductivity of the riverbed affected the initial distribution of hydraulic gradients.However, the difference in the initial conditions was negligible (4 % between the extreme scenarios) compared with the differences in bank storage caused by the presence of a clogging layer with variable hydraulic conductivity.Therefore, possible effects on bank storage can be ignored.
-In order to avoid over-parameterisation of the model and to limit the effect of sparse information availability regarding regional groundwater gradients, precipitation and evapotranspiration were approximated using a constant recharge flux term.
-The hydraulic effects of the river flow were not included in the model.
-Viscosity effects on hydraulic conductivity (Ma and Zheng, 2010) were excluded due to the small temperature difference between river and groundwater (max.10 • C) at the site.Ehlert (2014) has shown that a 24 % increase in hydraulic conductivity is possible due to this temperature difference.However, field measurements indicated solely conductive heat transport (Siergieiev et al., 2014a), due to attenuation of the advective-dispersive heat transfer by the clogging layer.

Conclusions
Bank hyporheic exchange was simulated using a field case scenario in an alluvial aquifer adjacent to the hydropowerregulated Lule River.The modelling showed that ecosystem requirements in terms of river-aquifer exchange flux are satisfied during 40 % of all wave events during the studied year.Discharge waves with longer duration and increased amplitude are essential for this site to improve hydrological exchange across the river-aquifer interface and bank hyporheic water quality.The combination of realistic and theoretical models improved current process understanding of hyporheic exchange in free-flowing and regulated rivers.Hypothetical scenarios included variable river discharge wave (duration and amplitude), river bank slope and hydraulic conductivity of the river-aquifer interface.Bank storage increased with lower bank slope, indicating the necessity of correct data on geometry of the river-aquifer interface when modelling surface water-groundwater interactions.Hydraulic conductivity of the riverbed positively affected bank storage.However, the influence on the residence time was not always consistent.Higher amplitude and longer wave duration increased bank storage, although larger maximum fluxes were observed for shorter waves at a given amplitude.There will be always a unique relationship between bank storage or residence time and the duration/amplitude wave ratio, which depends on the hydrogeological conditions.Hence, hyporheic exchange suppressed by colmation processes or flow manipulation can be improved by periodically releasing river discharge waves that are optimised for the specific river reach.This type of modelling offers a platform for understanding the transfer of nutrients across river-aquifer interfaces in natural and hydropower-impacted catchments.It further demonstrates the importance of conceptualisation (e.g.process consideration, boundary conditions, data necessity) and the limitations of numerical models.

Figure 1 .
Figure 1.Location of the observation site and cross section of the aquifer with groundwater wells, soil depth profiles and extent of the hyporheic zone.Numbers in well names indicate distance to the mean shoreline in metres.

Figure 2 .
Figure 2. Conceptual model of the site showing boundary conditions (BC), clogging layer and observation points.

Figure 3 .
Figure 3. Summary matrix of wave amplitude (m) and duration (h) as a fraction of all observed wave events.

Figure 5 .
Figure5.Hydraulic head at observation wells L5 (above) and L25 and the river (below) compared with the simulated results using calibrated parameters (see Table1).

Figure 6 .
Figure 6.Exchange flux for different bank slope (a), hydraulic conductivity (b), wave amplitude (c) and duration (d) scenarios, with the reference case (red) and insets of the input wave pulse.Note different scale.

Figure 7 .
Figure 7. Bank storage for different bank slope (a), hydraulic conductivity (b), wave amplitude (c) and duration (d) scenarios, with the reference case (red) and insets of the input wave pulse.Note different scale.

Figure 8 .
Figure 8. Return/fill time ratio (t R /t F ) and residence time for different bank slope (a), hydraulic conductivity (b), wave amplitude (c) and duration (d) scenarios, with the reference case (red) and insets of the input wave pulse.

Figure 9 .
Figure 9. Duration/amplitude ratio in relation to bank storage, residence time and return/fill time ratio (t R /t F ) for all waves exceeding 0.1 m amplitude.

Table 1 .
Model parameters measured and calibrated using FePEST.
* Mean of all measurements ± standard deviation.