Modeling storm events to investigate the influence of the stream‐catchment interface zone on stream biogeochemistry

We formulate a new mixing model to explore hydrological and chemical conditions under which the interface between the stream and catchment interface (SCI) influences the release of reactive solutes into stream water during storms. Physically, the SCI corresponds to the hyporheic/riparian sediments. In the new model this interface is coupled through a bidirectional water exchange to the conventional two components mixing model. Simulations show that the influence of the SCI on stream solute dynamics during storms is detectable when the runoff event is dominated by the infiltrated groundwater component that flows through the SCI before entering the stream and when the flux of solutes released from SCI sediments is similar to, or higher than, the solute flux carried by the groundwater. Dissolved organic carbon (DOC) and nitrate data from two small Mediterranean streams obtained during storms are compared to results from simulations using the new model to discern the circumstances under which the SCI is likely to control the dynamics of reactive solutes in streams. The simulations and the comparisons with empirical data suggest that the new mixing model may be especially appropriate for streams in which the periodic, or persistent, abrupt changes in the level of riparian groundwater exert hydrologic control on flux of biologically reactive fluxes between the riparian/hyporheic compartment and the stream water.


Introduction
[2] Recently, stream hydrology and biogeochemistry studies have focused on modeling approaches to determine the potential influence of the stream-catchment interface (SCI) on the dynamics of dissolved nutrients [Runkel et al., 2003]. Physically, these zones are interfaces between the stream and its hyporheic and riparian sediments with a continuous exchange of water and dissolved inorganic and organic solutes [Butturini et al., 2003;Jones and Holmes, 1996]. The stream-groundwater interface has been studied primarily at the hyporheic sediment scale and at the stream reach scale [Boulton et al., 1998]. However, recent new field studies have dealt with the influence of SCI on nitrate chemistry in streams at the catchment scale [Vidon and Hill, 2004;Wigington et al., 2003].
[3] In catchment-scale field studies, the mass balance mixing model (hereafter, Mix mod ) is widely used to evaluate the contribution of different flow components to streamflow during storms and, consequently, to discern the sources of each solute [Christophersen et al., 1990]. The flow components that feed streams during a storm are associated with the rapid hillslope water (overland flow, subsurface flow through soil macropores, throughfall, and direct channel interception) and with the slow infiltrated water (groundwater) [Chanat and Hornberger, 2003]. The slow subsurface water input through soil has been identified as a third important source of runoff in several humid catchments [McHale et al., 2002;Hornberger et al., 2001;Rice and Hornberger, 1998]. According to Evans and Davies [1998], the solute concentration-discharge hysteresis (c-q hysteresis) observed during a storm can be explained using the Mix mod model with two or three input components assuming a constant solute concentration in the input components. This conventional Mix mod model is a simplification of complex hydrological processes generated at the catchment scale. For instance, the assumption of constant solute concentration in time and space of the input components is arbitrary in most cases. This assumption adds uncertainty both to the interpretation of the c-q hysteresis [Bonell, 1998] and to quantitative estimates of the relative contribution of input components to the total stream runoff [Bishop et al., 2004;Chanat et al., 2002;Joerin et al., 2002;Hornberger et al., 2001].
[4] Historically, hydrological influence of the SCI has been ignored in the classic mass balance mixing model. In this model, the stream is viewed essentially as a conduit receiving water and solutes from the catchment [Bencala, 1993]. However, field observations suggest that the SCI can be relevant in regulating stream runoff in [McGlynn and McDonnell, 2003;. Recently, Chanat and Hornberger [2003] have integrated the hydrological storage effect of SCI into a two-component mixing model. In this model, the two input components had mixed within the SCI before entering the stream and water flow was assumed to be unidirectional from SCI to stream. Both modeling and field studies at the reach-scale have demonstrated that storm events cause dramatic changes in the stream catchment interface [Serrano and Workmann, 1998].
Briefly, storms provoke in the SCI (1) abrupt changes in groundwater levels, with consequent saturation/inundation of the riparian soil, (2) occurrence of reverse fluxes in the subsurface SCI in response to changes in the hydrologic gradient, and (3) rapid expansion and shrinkage of the boundary of the SCI in response to flooding and drying [Butturini et al., 2003]. These hydrologic processes may control the availability of reactive solutes, such as nitrate and dissolved organic carbon (DOC), in small streams. DOC and NO 3 are typically flushed during storms. Pulses in response to storms have been attributed to solute release from catchment organic soils generating higher solute concentrations in the hillslope water than in groundwater [Sickman et al., 2003;Carey, 2003;Brown et al., 1999;Hornberger et al., 1994]. Under these particular conditions, DOC has been occasionally used as a conservative solute quantifying the contribution and timing of hillslope water during storm events [Katsuyama and Ohte, 2002;Brown et al., 1999]. Alternatively, other studies have pointed out the relevance of soils/sediments in the stream-catchment interfaces important sources of solute components during storms [Bechtold et al., 2003;Creed and Band, 1998;Fiebig et al., 1990]. In this paper we examined how concentrations of reactive solutes in stream water change when the SCI compartment, a source of solutes, is coupled with the conventional mass balance mixing model.
[5] The main objectives of this study are (1) to generate a new mass balance mixing model that includes the streamcatchment interface compartment and the release of reactive solutes from sediments to the infiltrated water (hereafter, SCI-Mix mod ), (2) to explore, through a Monte Carlo approach, hydrological and chemical conditions which help discern the influence of the SCI, separate from catchment inputs, on stream reactive solutes during storms, and (3) to compare the outputs of the Mix mod and the SCI-Mix mod models using DOC and nitrate data obtained from two Mediterranean streams Sabater, 2000, 2002;Bernal et al., 2002]. The influence of the SCI on stream solute dynamics was separated from that of the catchment by comparing the main features of the c-q hystereses from simulations with both models. The latter assumes solutes behave conservatively. We also compare the general patterns of scatterplots of solute concentration versus discharge obtained as outputs from the two models and from field data of DOC and nitrate concentrations obtained in two Mediterranean streams.

Model Formulation: Conventional Mixing Model (Mix mod )
[6] The mixing of waters and solutes from different input sources during a storm episode can be simulated by a mass balance approach (Figure 1a): where Ws and Ms are the water volume (m 3 ) and the solute mass (mol) in the stream, respectively. Q(t) and M(t) are the a b c Figure 1. Schematic representation (a) of the conventional mixing model, Mix mod, and (b and c) of the modified mixing model integrating the stream-catchment interface, SCI-Mix mod , and solute release from the sediment compartment. For symbols, see details in the text.
input waters (m 3 /s) and the solute input fluxes (mol/s), respectively. The sub indices ''H'' and ''G'' indicate the rapid Hillslope and the slow Groundwater flow inputs, respectively [Chanat and Hornberger, 2003]. The index ''t'' indicates that input fluxes change temporally during a storm episode. The terms Qout(t) and Mout(t) are output water and solute flux, respectively, from the stream compartment.
2.2. Model Formulation: Accounting for the Stream Catchment Interface in the Mixing Model (SCI-Mix mod ) [7] To include the interface between the stream and catchment in the Mix mod model (equation (1)), the stream was divided into two subcompartments. These were connected by bidirectional water exchange between the surface stream channel and the interstitial SCI ( Figure 1b). The interstitial water is in close contact with the hyporheicriparian sediments (hereafter SCI sed ) ( Figure 1c): with the following initial and boundary conditions: where W SC and M SC are the water volume (m 3 ) and solute mass (mol) in the interstitial SCI. WS(0) and WSC(0) are the initial water volumes (m 3 ) under basal discharge conditions in stream and SCI zone, respectively. Q basal is the discharge at base flow (m 3 /s). Q HPeak and Q GPeak are the peak discharges (m 3 /s) of hillslope and groundwater respectively. In this model, the hillslope water (Q H (t)) flows directly into the stream compartment, while the groundwater (Q G (t)) flows into the stream through the SCI. a S and a SC (1/s) are the first-order exchange rates between the stream and the interstitial waters. The first-order exchange rates describing bidirectional water fluxes between stream and SCI (a S , a SC ) imply instantaneous mixing [Wörman et al., 2002]. The formulation of a S and a SC into the SCI-Mix mod is similar to that appearing in the advection-dispersion-transient storage model [Stream Solute Workshop, 1990]. To apply this model, it is recommended to work under base flow discharge conditions (for a remarkable exception, see Runkel et al. [1998]). Therefore the water masses are stationary (dW S /dt = 0 and dW SC /dt = 0 in equations 2a and 2b), and water exchanges are at equilibrium (i.e., a S W S = a SC W SC ). The SCI-Mix mod model, in contrast, assumes that, during storm events, water fluxes between the SCI and the stream are not in balance (i.e., a S W S 6 ¼ a SC W SC ) because the water masses of both change continuously.
[8] M S (0) and M SC (0) are the initial solute masses under base flow conditions. C S (0) and C SC (0) are the initial solute concentrations in stream and SCI compartments (mol/m 3 ) respectively. M H (t) and M G (t) are the hillslope and groundwater solute fluxes (mol/s). M HPeak and M GPeak are the peak solute input fluxes: where C H and C G are the solute concentrations (mol/m 3 ) of hillslope and groundwater inputs, respectively.
[9] In the SCI-Mix mod model, we assume that solute release predominates over solute immobilization during storms. Immobilization is considered negligible and k rel (1/s) gives the first-order release coefficient of the solute S (mol) adsorbed onto sediments [Schnabel and Fitting, 1988;Lookman et al., 1995]. The first-order release rate represents a simplification of adsorption-desorption kinetics at the sediment-water interface [Selim and Amacher, 1997]. However, the SCI compartment is a dynamic open system with continuous water flux in which the solute released from the sediment is continuously removed from the interstitial water. Studies using experimental soil columns with a continuous water flow have made a similar assumption, thus estimating a release rate which is independent of solute concentrations [Sparks et al., 1980].
[10] In order to account for the vertical heterogeneity in solute content, the compartment SCI sed was divided into ''n'' cells (DSCI sed (i), where 1 i n), with each ''i'' cell contains an initial solute mass (DS i (0)). Solute concentration is assumed to be higher in the upper soil layer and to decline through the soil profile [Bishop et al., 2004;Creed et al., 1996]. This approach allowed us to introduce solute gradients down the sediment profile. Using a similar approach, Hornberger et al. [1994] simulated the DOC flushing in an alpine stream.
[11] During a storm event, the water mass present in the SCI varies continuously (i.e., dW SC /dt 6 ¼ 0 in equation (2b)). Because of this, when dW SC /dt > 0 at an initial stage, the sediment in the SCI (SCI sed ) gradually saturates and subsequently, when dW SC /dt becomes less than 0, the SCI sed dries. Consequently, the total solute mass that is susceptible to release from SCI sed at each time ''t'' (S(t)), is a function of (1) the initial content of solute mass (DS i (0)) in each cell, DSCI sed (i), (2) the volume of interstitial water (W SC ), (3) the release rate k rel , and (4) the time of saturation of each cell.
[12] The number of cells ''n'' is identical to the time ''T'' required for maximal saturation of the SCI sed compartment (i.e., when the volume of interstitial water W SC is maximal). Thus a cell ''i'' is saturated (or dried) at each time ''t''. Formally, the total solute mass at each time (S(t)) is estimated as follows: where S(T) represents the total solute mass available in the SCI sed at time ''T''.
[13] Figure 2 shows the dynamics of S, with water mass in the SCI varying in time, obtained by using equation (4) under three different scenarios of initial sediment profile of solute masses (lineal, random and exponential). A gradual increase is observed in the solute mass (S) adsorbed onto sediments during the saturation of the SCI sed (i.e., when t < T). This general increase in S reflects the accumulation of DS i (0) elements, which is implicit in equation (4). On the other hand, the time of saturation of each cell, coupled with the value of k rel , determines the mass of solute remaining adsorbed onto the sediments when SCI sed is drying. Thus the relationship S versus W SC shows a hysteretic pattern.
[14] To compare results obtained with the SCI-Mix mod and the Mix mod models, it was important that the solute release flux (dS/dt) from the SCI sediments could be easily compared to that of the groundwater input (M G ) flowing through the SCI media. Therefore we defined as an input parameter the maximum solute release flux occurring at t = T: {dS/dt} T . In the present study, {dS/dt} T was related to the solute mass peak of the infiltrated input water (M GPeak ). If {dS/dt} T = 10%, that meant that {dS/dt} T = 0.1M GPeak .
[15] In order to consider the dynamics of hydrologic conditions, we assumed that a S and a SC were functions of the water volumes in the stream and SCI compartments (W S , W SC ), following the equations: where a S0 and a SC0 (1/s) are the exchange rates between stream and SCI under basal conditions, when water fluxes between the stream and the SCI compartments are assumed to be at equilibrium (equation (5c)). b (1/s) is the maximum exchange rate. According to equation (4), a S and a SC range between:

C-q Hysteresis Descriptors and Model Simulations
[16] To compare the c-q hystereses simulated with the SCI-Mix mod and the Mix mod models, the two models were run simultaneously assuming identical peaks of water inputs (Q HPeak and Q GPeak ), identical solute concentrations in water inputs (C H and C G ), and identical initial water mass in the stream channel (WS (0)).
[17] The c-q hystereses obtained with the two models were compared using three descriptors summarizing their main features (Figure 3). The similarity of general trends in c-q hystereses obtained from the alternative models was estimated by the Pearson correlation coefficient (r) between the simulated solute concentrations. When r < 0 the two c-q Figure 2. Profiles of the initial distribution of solute mass in the stream catchment porous media of the stream catchment zone (SCI sed ), used in the Monte Carlo simulations, and associated evolution of the relative total solute mass (S* = S (t) /S (T) ) with respect to the water mass in the stream-catchment interface (W* SC = W SC (t)/W SC (T)). These profiles represent the example at which k rel = 0.0005 1/s, W S (0)/W SC (0) = 25, and Q HPeak /Q GPeak = 1.
hystereses followed opposite trends. The differences in the areas (A) and rotational patterns of the c-q hystereses were estimated as the product: where A SCI /A Mix is the ratio between the c-q hysteresis area obtained with the SCI-Mix mod model (A SCI , mol/s) and that obtained using the mixing model (A Mix , mol/s). The areas were estimated after standardizing discharges and solute concentrations to a unity scale. The term Rp summarizes the rotational pattern of the c-q hystereses. If the c-q hystereses obtained with the two models had the same rotational pattern (i.e., both were either clockwise or counter clockwise), then Rp = 1; otherwise Rp = À1. The rotational pattern of a c-q hysteresis may be ambiguous and the presence of a node in the middle of the c-q was considered to indicate an ambiguous rotational pattern, therefore Rp = 0. For instance, if DRp = À0.25, the c-q hysteresis obtained with the SCI-Mix mod model would be 25% sharper than that obtained with the Mix mod model and would have an opposite rotational pattern.
[18] The indices r and DRp summarized differences in the shape of the c-q hystereses obtained using the alternative models. Divergences between the two models were considered substantial when both indices were negative.
[19] The ratio between the absolute range in concentration observed with the SCI-Mix mod model ((C max À C min ) SCI ) and that observed with the Mix mod model ((C max À C min ) Mix ) summarized the differences in stream concentration changes estimated using each model: [20] When DC rel $ 1, the c-q hystereses in both models showed similar absolute ranges in concentration. For this study, we considered that differences in stream concentrations between the two models were substantial when DC rel > 2.
[21] The concentration versus discharge plots (hereafter referred to as scatterplots) obtained from using the SCI-Mix mod and the Mix mod models were drawn by extracting values of solute concentration at maximum discharge at each run of the two models. The general patterns of these scatterplots were compared with DOC and nitrate scatterplots observed in two Mediterranean streams.
[22] In order to explore how input parameters drove model outputs (i.e., the c-q hysteresis descriptors), the Mix mod and SCI-Mix mod models were run simultaneously 10 4 times using a Monte Carlo approach. In each simulation, we randomly assigned values to the six input param- by giving values between 0.001 and 2 times a nominal value arbitrarily selected for each parameter (Table 1).
[23] The volume of stream water is reportedly larger than that of the transient storage water [Butturini and Sabater, 1999]. Therefore the nominal value of the ratio W S (0)/ W SC (0) was taken as 25.
[24] For those parameters that were included in both models (Q HPeak /Q GPeak and C H /C G ), the nominal values were unity. Thus the simulations were uniformly distributed within the Q HPeak /Q GPeak versus C H /C G plane. Ranges for the input parameters were selected with an aim to covering wide spectra of hypothetical conditions. The range of the Q HPeak /Q GPeak ratio included more extreme conditions than those reported in studies from humid catchments, where groundwater and subsurface inputs usually predominate over total storm runoff [Bazemore et al., 1994;Brown et al., 1999;Hornberger et al., 2001]. The magnitudes of the hydrological (W S(0) /W SC(0) , b) and chemical (k rel , {dS/dt} T ) parameters specifically introduced in the SCI-Mix mod model are not known, therefore the Monte Carlo approach avoided selecting a subjective combination of values for these Figure 3. Schematic description of the c-q hysteresis descriptors (r, DR p , and DC rel ). C* and Q* are concentration and discharge standardized between 0 and 1. The crossed circle indicates a node. parameters [Kalos and Whitlock, 1986]. The distribution of initial solute mass in the SCI compartment (SDS i (0)) was selected randomly among the three hypothetical distributions shown in Figure 2.
[25] After conducting the Monte Carlo simulation and calculating output descriptors, we estimated the probability of obtaining divergences between the two models in the Q HPeak /Q GPeak versus C H /C G plane. In this plane, four regions can be identified (Figure 4). When Q HPeak /Q GPeak > 1 (regions A and B), the peak of hillslope water is larger than that of the groundwater. When C H /C G > 1 (regions A and D), the concentration of the hillslope water is higher than that of groundwater. Therefore the region A represents the case in which the water flux peak and the solute concentration are larger in the hillslope water than in the groundwater (Q HPeak /Q GPeak > 1 and C H /C G > 1). Region C represents the opposite condition.

Study Sites Description and Stream Water DOC and Nitrate Data
[26] The chemical data are from two Mediterranean streams in northeastern Spain called Fuirosos and Riera Major with similar total area (15.5 and 16 km 2 , respectively) and with negligible human activity. Both streams drain forested catchments and saturated areas in the riparian zone and hillslope are nonexistent.
[27] In Riera Major, the discharge is permanent and base flows range between 25 and 60 L s À1 in summer and winter, respectively [Butturini and Sabater, 2000]. The hyporheic zone is rather shallow (between 0 and 1 m depth) with rapid water exchange between stream and hyporheic waters [Butturini and Sabater, 1999]. The quantity of particulate organic matter in the streambed ranged between 0 and 60 g DW m À2 , with a peak during the short, leaf fall period, in October [Romaní, 1997].
[28] In Fuirosos, discharge is intermittent, with a long, summer dry period followed by a short but intense stream recharge period in late summer-early fall and a subsequent humid period lasting until late spring [Butturini et al., 2003]. The hyporheic-riparian compartment in Fuirosos is much deeper than in Riera Major and the granitic bedrock is located up to 3 -6 m below the streambed and the riparian soil [Butturini et al., 2003]. During the stream recharge period, the stream water can infiltrate 10 m into the riparian zone [Butturini et al., 2003]. In summer, the riparian tree community suffers intense hydric stress, causing a marked input of up to 450 g DW m À2 of detritus which accumulates on the streambed and stream margins [Sabater et al., 2001].  [29] DOC and nitrate concentrations typically increased during storm events in both streams. In Riera Major, the DOC and NO 3 -N concentrations ranged between 1 and 8 ppm and 0.15 and 2.7 ppm respectively. In groundwater, DOC and NO 3 -N concentrations were in the same range as that measured in stream water Sabater, 2000, 2002]. In Fuirosos, DOC and NO 3 -N concentrations in stream water ranged between 3 and 20 ppm and 0.01 and 3 ppm respectively [Bernal et al., , 2004. In ground-water, NO 3 -N concentration was virtually nil, while DOC concentrations were typically less than 3 ppm. On the other hand, DOC and NO 3 -N concentrations in hillslope water flowing through organic forest soil during rains averaged 35 ppm and 3 ppm, respectively (unpublished data from authors).

Monte Carlo Simulations
[30] Figure 5 shows a visual example of the differences in the c-q hystereses obtained with the conventional Mix mod model and with the SCI-Mix mod model under the assumption of a relatively low solute input from the SCI interface ({dS/dt} T = 50%). In Figures 5a -5c the ratio of concentrations of the input water components (C H /C G ) was set at 0.5, while the ratio Q HPeak /Q GPeak was set at 0.5, 1 and 2, respectively. Under these assumptions, differences between the two models in the c-q hystereses were more evident when the groundwater was larger than the hillslope water (i.e., Q HPeak /Q GPeak < 1, Figure 5a). On the other hand, in Figures 5d-5f, the ratio Q HPeak /Q GPeak was fixed at 0.5, while the ratio C H /C G was at 0.5, 1 and 2, respectively. In this case, when the solute concentration in hillslope water was larger than in the groundwater (i.e., C H /C G > 1, Figure 5f) the c-q hystereses that were obtained with the two models, differ in their concentration changes, but not in their shapes and rotational patterns.
[31] The Q HPeak /Q GPeak versus C H /C G plot represents the values of the three c-q hysteresis descriptors obtained during the Monte Carlo experiment, summarizing the differences between the two models ( Figure 6). The probability of obtaining different results by using the Mix mod model or the SCI-Mix mod model was not uniformly distributed within the Q HPeak /Q GPeak versus C H /C G plane. The two models followed opposite global c-q hysteresis trends (r < 0) in 20% of cases, and 58% of these cases were located in region C, while none were found in regions A and D. The agreement between the two models (r > 0.5) was more probable in region A (Figure 6a).
[32] DRp showed a similar distribution to that of r ( Figure 6b). The two models had opposite rotational patterns (DRp < 0) in 11% of cases and 57% of these cases were in region C. Again, in regions A and D the likelihood of finding DRp < 0 was nil. Nevertheless, c-q hystereses with ambiguous rotational patterns (DRp = 0) amounted to 31% of the total number of cases and these were uniformly scattered within the Q HPeak /Q GPeak versus C H /C G plane. Approximately 72% of the values of the ratio A SCI /A Mix were lower than 0.75, indicating that the SCI-Mix mod model favored sharper c-q hystereses than the Mix mod model did. The 90% of c-q hystereses with ambiguous rotational pattern were associated with values of A SCI /A Mix < 0.5. The probability of simultaneously observing DRp < 0 and r < 0 was 9% of the total number of cases and 65% of these cases were in region C.
[33] Substantial differences in the ranges of simulated stream concentrations using the two models (DC rel > 2) were observed in 31% of cases, and 77% of these were in regions C and D. High DC rel values (DC rel > 20) were found almost exclusively in the region C. With the increase in the proportion of hillslope water with respect to groundwater entering the stream (i.e., when Q HPeak /Q GPeak > 1), the DC rel > 2 values emerged around the condition C H /C G $ 1 (Figure 6c).
[34] The magnitude of the SCI solute flux ({dS/dt} T ) was a key variable modifying the c-q hysteresis shapes and driving the probability of obtaining divergent results between both models. Figure 7 represents the distribution, within the Q HPeak /Q GPeak versus C H /C G plane, of cases when DC rel > 2, and both r and DR p were negative under different SCI solute flux ranges. Increases in {dS/dt} T clearly increased the chance of observing differences between the two models: the probability of observing DC rel > 2 increased from 8% to 53%, with most cases located in regions C and D, while the probability of simultaneously observing DRp < 0 and r < 0 increased from 3 to 14%, with cases being located in regions B and C.

From Simulations to Empirical Data: Comparison of Scatterplots
[35] The Monte Carlo simulations showed that divergences between the two models was most likely when (1) Q HPeak /Q GPeak < 1 (i.e., regions C and D) and (2) {dS/dt} T > 100%. Specifically, in region C (C H /C G < 1 and Q HPeak /Q GPeak < 1), the three c-q hysteresis descriptors helped to detect divergences between the two models (i.e., DC rel > 2; DR p < 0 and r < 0) while in region D (C H /C G > 1 and Q HPeak /Q GPeak < 1), solely the descriptor DC rel identified disagreements between the two models. This implies that, when solute concentrations in the rapid hillslope water is higher than in groundwater (i.e., C H /C G > 1), the analysis of the dispersion of solute concentrations with respect to discharge (i.e., the scatterplot) is essential for comparing the SCI-Mix mod outputs with the empirical data. Under this condition, the scatterplot obtained from the Mix mod model showed a general positive trend, with a gradual spreading of concentration data points with increases of discharge (Figure 8a, left). On the other hand, this characteristic pattern did not emerge with the SCI-Mix mod model because of solute peaks at low discharges (Figure 8a, right).
[36] The scatterplots for Riera Major show the highest DOC and NO 3 concentrations at moderately high discharges (i.e., between 300 and 400 L s À1 ). At lower discharges, their patterns are similar to those expected with the conventional Mix mod model, while at discharges higher than 400 L s À1 , the DOC and NO 3 dispersions decreased markedly (Figure 8b, left). In Fuirosos, the DOC and NO 3 scatterplots obtained from field data differed greatly from plots expected by using the Mix mod model. Indeed, the highest dispersions occurred at low discharges (Figure 8b, right). The DOC variability decreased for discharges higher than 150 L s À1 , while the variability of NO 3 remained relatively high over the entire discharge range. As for DOC, the solute peaks were frequently associated to the discharge rising limb of storm events monitored in early autumn, during the stream recharge period. Later, solute concentrations decreased rapidly, generating a c-q hysteresis with a marked clockwise rotational pattern. These DOC-q hystereses are very different from hystereses generated by storms of similar magnitude monitored during the successive humid period. Then, DOC concentrations are lower and hystereses show sharp areas and unclear rotational patterns (Figure 9a).

Discussion
[37] This study describes a new mass balance model that incorporates the chemical and hydrological properties of the interface between the stream and the catchment. Much of the current interest in the stream-catchment interface is related to findings that the riparian soil/sediment may mitigate diffuse nitrate pollution from shallow groundwater [Vidon and Hill, 2004]. Hill [1996] pointed out that most uncertainties about the effective role of riparian zones in nitrate removal derive from insufficient attention to hydrology. In fact, this conclusion is from studies performed essentially at reach scales and under base flow discharge . There have also been studies analyzing the biogeochemical behavior of the riparian zone during storms, which have reported nitrate release rather than removal [Konohira et al., 2001;Butturini et al., 2003]. Potential influences and effects of the interface compartment over stream biogeochemistry are increasingly recognized. Yet the conventional mixing model ignores its hydrological and chemical effects. This leads to uncertainties in interpretations of response during storm events, of both reactive (this study) and conservative solutes [e.g., Chanat and Hornberger, 2003]. Thus research on the role of the stream-catchment interface has stimulated reviews of the classical mixing model formulation.
[38] The results of the Monte Carlo simulations reported here reveal that the probability of discerning the influence of the SCI on stream carbon and nitrogen flux during storms increases when the storm runoff is dominated by a groundwater component flowing through the SCI before entering the stream. Under the common situation where the concentration of reactive solutes (such as DOC and NO 3 ) in the hillslope water exceeding that in the groundwater (i.e., C H /C G > 1) [Sickman et al., 2003;Brown et al., 1999;Hornberger et al., 1994], no clear differences between models (SCI-Mix mod and the Mix mod ) were apparent in the rotational patterns of the c-q hystereses. Evans and Davies [1998] recognized the problematic interpretation of flushing when C H /C G . Indeed, conflicting interpretations can be found in the literature. For instance, clockwise DOC and NO 3 -q hystereses were observed in small humid catchments and have been associated with solute flushing in the upper organic-rich soil horizon [Creed and Band, 1998;Figure 8. Solute concentration versus discharge scatterplots (a) obtained from the Monte Carlo simulations, under the condition C H /C G > 1, by using the Mix mod and SCI-Mix mod models and (b) obtained from field observations of NO 3 and DOC concentrations in two pristine Mediterranean streams. C* and Q* in Figure 8a are concentrations and discharges standardized between 0 and 1. Each dot in Figure 8a is a run of each of the two models. Hornberger et al., 1994]. In contrast, Meyer and Tate [1983] have suggested that the clockwise DOC-q hysteresis could be caused by solute mobilization from organic matter stored in the streambed and hyporheic sediments. Figures 9b and 9c illustrate the ambiguity in the interpretation of clockwise c-q hystereses monitored in Fuirosos. Figure 9b shows that the sharp DOC-q hysteresis monitored in this stream during the humid period (23 October 2000) can fit reasonably well the conventional Mix mod model under the assumption that: C H /C G = 35/4, Q HPeak /Q GPeak = 0.5, and that peaks of the two components coincided in time. In this example, C G = 4 ppm is the DOC concentration at basal discharge when streamflow is feed by groundwater, while C H = 35 ppm is the estimated DOC concentration in hillslope water [unpublished data from authors]. If we assume that values of C H /C G , Q HPeak /Q GPeak and the hydrograph template of the storm event occurring in the stream discharge period (19 September 2000) are identical, we see that the SCI-Mix mod model captures reasonably well the abrupt DOC peak occurring during the discharge rising limb and that this model generates a marked clockwise DOC-q hysteresis when {dS/dt} T = 200% and solute adsorbed on sediment is randomly distributed through the sediment profile ( Figure 9b). Nevertheless, the conventional Mix mod model can simulate a clockwise DOC-q hysteresis using identical C H /C G and Q HPeak /Q GPeak values and anticipate the peak of hillslope water with respect to the groundwater (Figure 9c). The c-q hysteresis obtained with SCI-Mix mod simulated better the global trend of DOC-q hysteresis than the hysteresis obtained using the conventional Mix mod model. However, the c-q hystereses obtained by using the two models have the same rotational patterns and consequently follow similar overall trends. Neither of these features are a proof of the relevance of SCI for stream chemistry. To evaluate the importance of including the stream-catchment interface in models describing actual c-q hystereses of reactive solutes it is essential to integrate the information obtained from c-q hystereses monitored during different hydrologic periods with information obtained from more specific studies focused on the hydrology and biogeochemistry of the SCI compartment. In Fuirosos, the stream recharge period is characterized by extended water exchange between the stream and the surrounding riparian/hyporheic compartment [Butturini et al., 2003] and the streambed is covered by large amounts of organic matter accumulated during the previous dry period [Sabater et al., 2001]. This information specific to Fuirosos suggests that these typically marked clockwise DOC-q hystereses, with high DOC concentrations, which occur during the stream discharge period, might reflect the immediate DOC flushing from the organic matter stored on the streambed rather than an input from water flowing over the forest hillslopes.
[39] Moreover, the Monte Carlo simulations suggested that, if detailed information about the hydrology and biogeochemistry of the SCI compartment is unavailable, an analysis of the general patterns of the scatterplots of solute concentration versus discharge may be useful as a preliminary tool for discerning whether solute dynamics in a stream have a SCI control or not. DOC and NO 3 scatterplots from Riera Major, or from other small catchments in humid regions [Brown et al., 1999;Creed et al., 1996], fit into that expected from the Mix mod model, suggesting that the  1). In this simulation the solute was concentration changed randomly through the sediment profile (see Figure 3). (c) Comparison of the DOC-q hysteresis monitored on 19 September 2000 with that simulated (shaded dashed line) with the conventional Mix mod model (Q HPeak /Q GPeak = 0.5, C H /C G = 9). Arrows indicate the rotational pattern of the c-q hysteresis. The insets in Figures 9b and 9c show the hydrograph template used in the simulations (the solid line is Q H , and the dashed line is Q G ). dynamics of these solutes are essentially driven by the catchment, as opposed to the SCI. Moreover, the DOC and NO 3 -q hystereses that were reported for Riera Major are typically counterclockwise . These rotational patterns corroborated the hypothesis that the SCI does not play a relevant role in the dynamics of these solutes during storms.
[40] Conversely in Fuirosos, the SCI-Mix mod model can explain better than the conventional model the high DOC and NO 3 concentration variability observed at low discharges and suggests that the SCI may have exerted significant control on DOC and NO 3 .

Conclusion
[41] We have presented a mass balance model (SCI-Mix mod ) which integrates the chemical and hydrological properties of the stream-catchment interface with the conventional two components mixing model. Field and theoretical approaches provide evidence of the influence of the hyporheic and riparian zones on stream biogeochemistry at the stream reach scale (for a review, see Runkel et al. [2003]). Simulations with the SCI-Mix mod and the comparison with empirical data from two Mediterranean streams with contrasting DOC and NO 3 concentration patterns, shows that studies of the riparian/hyporheic zone will benefit from information about storm runoff events of low magnitude dominated by the groundwater flow component. The SCI-Mix mod permitted an analysis of the relevance of this interface on stream chemistry from a catchment-scale perspective. This model may be especially appropriate for streams in which the periodic, or persistent, abrupt changes in the level of the riparian groundwater exert a hydrologic control on chemical connections between the riparian/hyporheic zone and the stream water.