Modeling nutrient in-stream processes at the watershed scale using Nutrient Spiralling metrics

One of the fundamental problems of using large-scale biogeochemical models is the uncertainty involved in aggregating the components of ﬁne-scale deterministic models in watershed applications, and in extrapolating the results of ﬁeld-scale measurements to larger spatial scales. Although spatial or temporal lumping may reduce the problem, 5 information obtained during ﬁne-scale research may not apply to lumped categories. Thus, the use of knowledge gained through ﬁne-scale studies to predict coarse-scale phenomena is not straightforward. In this study, we used the nutrient uptake metrics deﬁned in the Nutrient Spiralling concept to formulate the equations governing total phosphorus in-stream fate in a watershed-scale biogeochemical model. The ratio- 10 nale of this approach relies on the fact that the working unit for the nutrient in-stream processes of most watershed-scale models is the reach, the same unit used in ﬁeld research based on the Nutrient Spiralling concept. Automatic calibration of the model using data from the study watershed conﬁrmed that the Nutrient Spiralling formulation is a convenient simpliﬁcation of the biogeochem- ical transformations involved in total phosphorus in-stream fate. Following calibration, the model was used as a heuristic tool in two ways. First, we compared the Nutrient Spiralling metrics obtained during calibration with results obtained during ﬁeld-based research in the study watershed. The simulated and measured metrics were similar, suggesting that information collected at the reach scale during research based on 20 the Nutrient Spiralling concept can be directly incorporated into models, without the problems associated with upscaling results from ﬁne-scale studies. Second, we used results from our model to examine some patterns observed in several reports on Nutrient Spiralling metrics measured in impaired streams. Although these two exercises involve circular reasoning and, consequently, cannot validate any hypothesis, this is a 25 powerful example of how models


Introduction
Excess human-induced nutrient loading into rivers has led to freshwater eutrophication (Vollenweider, 1968;Heaney et al., 1992;Reynolds, 1992) and degradation of coastal areas and resources on a global scale (Walsh, 1991;Alexander et al., 2000;McIsaac et al., 2001). Thus, cultural eutrophication assessment and control are important is- 5 sues facing natural resource managers, especially in watersheds with high human impact. Control measures are frequently based on bulk calculations of river nutrient loading (e.g., Marcé et al., 2004), on crude mass-balance approximations (Howarth et al., 1996;Jaworski et al., 1992), on the nutrient export coefficient methodology (Beaulac and Reckhow, 1982), or on several refinements derived from it (Johnes, 1996;Johnes watershed applications (Rastetter et al., 1992), and in extrapolating the results of fieldscale measurements to larger spatial scales. This is important because the ability to use the knowledge gained through fine-scale studies (e.g. nutrient uptake rate for different river producers communities, nutrient fate in the food web, and so on) to predict coarse-scale phenomena (e.g. the overall nutrient export from watersheds) is highly 5 desirable. However, incorporating interactions between many components in a bigscale model can be cumbersome, simply because the number of possible interactions may be very large (Beven, 1989). The usual strategy to avoid a model including precise formulations for each interaction (and thus the counting of thousands of parameters) is to lump components into aggregated units. But although lumping might reduce the 10 number of parameters to a few tens, we still cannot guarantee that the information obtained during fine-scale research will apply to lumped categories. The behavior of an aggregate is not necessarily equivalent to the sum of the behaviors of the fine-scale components from which it is constituted (O'Neill and Rust, 1979).
Considering nutrient fate modeling at the watershed scale, the processes involved 15 in in-stream biogeochemical transformations are a major source of uncertainty. The working unit for the nutrient in-stream processes of most watershed-scale models is the reach. Within this topological unit, several formulations for biogeochemical reactions are included depending on the model complexity (e.g. adsorption mechanisms, algae uptake, benthic release, decomposition). But frequently modelers only have lim- 20 ited field information to parameterize these processes, and when this information is available, it usually comes from fine-scale research. The problem is that model formulations and the processes described at the field and their scales are not necessarily equivalent, and frequently the incorporation of field information in the model is not straightforward.
pirically quantifiable at the reach scale, then we will be able to apply the field research to the model without the problems associated with upscaling results from fine-scale studies. In the case of nutrient fate in streams, the Nutrient Spiralling concept (Newbold et al., 1981) could be a convenient simplification of the nutrient biogeochemical transformations involved, because the nutrient spiralling metrics are empirically eval-5 uated at the reach scale (Stream Solute Workshop, 1990), the same topological unit used by most watershed-scale models. Within this framework, the fate of a molecule in a stream is described as a spiral length, which is the average distance a molecule travels to complete a cycle from the dissolved state in the water column, to a streambed compartment, and eventually back to the water column. The spiral length consists of 10 two parts: the uptake length (S w ), which is the distance traveled in dissolved form, and the turnover length, which is the distance traveled within the benthic compartment. Usually, S w is much longer than turnover length, and research based on the nutrient spiralling concept focuses on it. S w is evaluated at the reach scale, with nutrient enrichment experiments (Payn et al., 2005), following nutrient decay downstream from 15 a point-source (Martí et al., 2004), or with transport-based analysis (Runkel, 2007).
In this study, we explored the possibility of using the mathematical formulation of the Nutrient Spiralling concept to define the in-stream processes affecting total phosphorus concentration in a customary watershed-scale deterministic model. First, we manipulated the model source code to include the nutrient spiralling equations. Then, 20 we implemented the model for a real watershed, and let a calibration algorithm fit the model to observed data to assess the performance of the model. In a second step, we analyzed whether the final model structure (i.e., the value of the adjustable parameters) were a realistic representation of the natural system. This consisted in a comparison between the adjusted nutrient spiraling metrics in the model and values from 25 field-based research performed in the watershed under study and in other systems worldwide.
Finally, it is worth noting that implementing the model in this manner (i.e., fitting the model to data instead of incorporating data from field-based research into the model)

505
we could use the model as a heuristic platform, discussing some patterns observed in the Nutrient Spiralling metrics measured in streams worlwide in the light of the results of our model. Of course, this was circular reasoning than could not validate any hypothesis, but given that a model cannot be used for formal testing anyway (Oreskes et al., 1994), we considered this procedure a much more interesting exercise.

Study site
We explored the possibility of using the Nutrient Spiralling formulation for the in-stream modules of a watershed-scale model in the Ter River watershed (Spain), including all watercourses upstream from Sau Reservoir (Fig. 1). Thus, we considered 1380 km 2 10 of land with a mixture of land use and vegetation. The headwaters are located in the Pyrenees above 2000 m a.s.l., and run over igneous and metamorphic rocks covered by mountain shrub communities and alpine meadows. Downstream, the watercourses are surrounded by a mixture of conifer and deciduous forest, and sedimentary rocks become dominant. The Ter River then enters the alluvial agricultural plain (400 m a.s.l.) 15 where non-irrigated crops dominate the landscape. The main Ter River tributaries are the Fresser River in the Pyrenees, the Gurri River on the agricultural plain, and Riera Major in the Sau Reservoir basin.
The Ter River watershed includes several urban settlements, especially on the agricultural plain (100 000 inhabitants). Industrial activity is also important, with numerous 20 phosphorus point-sources (Fig. 1A) coming from industrial spills and effluents from wastewater treatment plants (WWTP). Additionally, pig farming is an increasing activity, generating large amounts of slurry that are directly applied to the nearby crops as a fertilizer, at a rate of 200 kg P ha −1 yr −1 (Consell Comarcal d'Osona, 2003). The median flow of the river at Roda de Ter (Fig. 1) is 10 m 3 s −1 , and total phosphorus (TP) 25 concentration frequently exceeds 0.2 mg P L −1 .

2.2 Modeling framework
The main target of the watershed-scale model was the prediction of daily TP river concentration at Roda de Ter (Fig. 1A). We used the Hydrological Simulation Program-Fortran (HSPF), a deterministic, semi-distributed model that simulates water routing in the watershed and water quality constituents (Bicknell et al., 2001). HSPF simu-5 lates streamflow using meteorological inputs and information on several terrain features (land use, slope. . . ), and it discriminates between surface and subsurface contributions to streams. As a semi-distributed application, HSPF splits the watershed into different sub-basins (e.g., Fig. 1A). Each sub-basin consists of a river reach, the terrain drained by it, and upstream and downstream reach boundaries to solve for lotic 10 transport across the watershed. Only limited, very rough spatial resolution is considered inside sub-basins, and explicit spatial relationships are present only in the form of reach boundaries. HSPF solves the hydrological and biogeochemical equations of the model inside sub-basins, and the resolution of each sub-basin is hierarchically sorted in order to adequately simulate mass and energy transport as water moves downstream 15 ( Fig. 2). Hydrology and river temperature have previously been simulated and validated in the Ter River watershed using HSPF on a daily and hourly time scale . Figure 3 shows the simulated daily river streamflow and temperature against observations at Roda de Ter for sampling dates when river TP 20 concentration was available. For simulations included in this study, we used the water routing and river temperature results from  and , respectively. We also refer the reader to  for the sub-basin delineation procedure and other details of the semi-distributed model. 25 This section describes how point and diffuse sources of TP to stream reaches were calculated for each subbasin defined in HSPF (Figs. 1A and 2). TP concentration and 507 water load information for point sources comes from ACA, and consisted of a georrefenced, heterogeneous database with very detailed data for some spills, and crude annual values for others. As a result, we decided to include in the model an adjustable multiplicative factor for WWTP inputs (C w ) and another for industrial spills (C i ), in order to correct for potential monotonous biases in the database (Table 1). Thus, the daily 5 TP load from point sources for a particular reach was the sum of all spills located in the corresponding subbasin times the correction factor. Note that the correction factor value was the same for all spills of the same kind (i.e., industrial or WWTP) throughout the watershed.

Point sources and diffuse inputs of phosphorus
Diffuse TP inputs into the watercourses were modeled using water routing results 10 from . Since we were mainly interested in the in-stream processes, and in order to keep the model structure as simple as possible, we calibrated the model against river TP data collected on sampling dates for which there was no surface runoff for at least seven days previously. Thus, we ignored TP transport in surface runoff. TP concentration in interflow and groundwater flow (diffuse sources in Fig. 2) 15 was modeled assuming a power dilution dynamics. We modified the HSPF code to include the following formulations where TP i and TP g are TP concentration (mg P L −1 ) in interflow and groundwater dis-20 charge, respectively. Q i and Q g are the interflow and groundwater discharge (mm) coming from the land drained by the reach. a i , a g , b i , and b g are adjustable parameters. Note that we did not consider spatial heterogeneity for these parameters (i.e., a different adjustable value for each sub-basin). Thus, they should be considered as averages for the entire watershed. However, as we will see later, river TP data for cali-25 bration of the model came from one sampling point. As a consequence, the optimized parameter values will more closely correspond to the situation around this sampling point, and they will be less reliable far from it. HSPF includes a module to simulate the biogeochemical transformations of TP inside river reaches (i.e., the in-stream processes, Fig. 2B). Several processes can be defined in this module, including assimilation/release by algae, adsorption/desorption mechanisms, sedimentation of particulate material, decomposition of organic mate-5 rials, among others (Bicknell et al., 2001). The main objective of this study was to explore the possibility of simplifying all these in-stream processes using an aggregate process: TP retention as defined by the Nutrient Spiralling concept. We modified the HSPF code to include formulations that follow. The in-stream TP fate was modeled as a first order decay following the Stream Solute

10
Workshop (1990) and can be conceptualized as where t is time (s), x is distance (m), Q is river discharge (m 3 s −1 ), A is river crosssectional area (m 2 ), and k c (s −1 ) is an overall uptake rate coefficient. Q i and Q g are as in Eq.
(1) but expressed in m 3 s −1 . The first term of the equation refers to advection, the 15 second to dispersion, and third and fourth to lateral subsurface inflows. In the context of the HSPF modeling framework, all these terms refer to TP inputs to the reach, and were solved as explained above. The last term in Eq.
(2) simulates solute transfers between water column and benthic compartment (this is what we considered in-stream processes in this paper). Of course 20 this represents an extremely simplified formulation, and must be interpreted as a net transport, because more complex settings account for independent dynamics of benthic release and concentration in one or more benthic compartments (Newbold et al., 1983). One important limitation of this formulation is that k c is a constant, and applying a single value in a system with varying water depth may be very unrealistic. A much 25 more convenient formulation of the last term in Eq. (2) where h is river depth. Obviously, from this we can establish v f =h×k c , which implies that v f is a scale free parameter (Stream Solute Workshop, 1990). We modified the 5 HSPF code to incorporate this formulation as the only modeled in-stream process, also including a built-in HSPF temperature correction factor. The final formulation of the in-stream processes was where TC is the temperature correction factor and T w ( • C) is river water temperature. 10 Thus, the in-stream module of the watershed-scale model only included two adjustable parameters (Table 1). v f is related to the Nutrient Spiralling metric S w through the following relationship where u is water velocity (m s −1 ). Since nutrient uptake experiments in rivers and 15 streams usually report S w values for representative reaches, we can calibrate the watershed model with observed data and compare the obtained S w with reported values from real systems (including data from the Ter River watershed). Regarding Eq. (4), we are assuming that areal uptake rate (U=v f ×TP) is independent of nutrient concentration. Although a Monod function relating U and nutrient concentra-20 tion is usually applied for this purpose, high TP concentrations in the Ter River watershed streams are well established in the asymptotic section of the relationship (Mulholland et al., 1990). Although this is not a realistic assumption for some pristine reaches of the Ter River headwaters and Riera Major, it probably applies to reaches around the 510 Introduction

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version Interactive Discussion TP sampling point, to which model calibration will be more sensitive. Nonetheless, the saturation assumption would not apply to the particulate fraction of TP, beacuse the most important retention process of this phosphorus fraction (sedimentation losses) does not exhibit saturation kinetics. However, particulate phosphorus accounted for a minor part of the TP pool at the sampling point (36% on average), and no conspicu-5 ous effect of TP on U was expected in most occasions. However, when the river carries a lot of particulate material (i.e., during floods) the model structure as defined so far could not be entirely appropiate. Still regarding Eq. (4), we are assuming a monotonous effect of temperature on solute transfer in the range of water temperatures measured in our streams.

10
As above, note that we did not consider spatial heterogeneity for the nutrient retention parameters (i.e., different adjustable values for each reach defined in the HSPF model). Thus, adjusted Nutrient Spiralling metrics reported in this study (v f and S w ) should be considered as averages for the entire watershed. As in the preceeding section, optimized parameter values will more closely correspond to the situation around 15 the TP sampling point, and they will be less reliable as we move upstream.

Calibration strategy
River TP concentration data for this study came from the Sau Reservoir long-term monitoring program, which includes a sampling point upstream of the reservoir at Roda de Ter (Fig. 1A). Sampling was weekly to monthly, from January 1999 to July 2004. 20 Samples were analyzed using the alkaline persulfate oxidation method (Grasshoff et al., 1983). Among available data, we only considered 106 river TP concentration values measured on sampling dates for which there was no surface runoff for at least seven days previously (see Sect. 2.3). These data was the basic data used for calibration and validation of the HSPF model. In addition, TP data from 14 sampling stations run by the 25 local water agency (Agència Catalana de l'Aigua, ACA) were used as a supplementary set for model verification (Fig. 1A). The amount of data from these stations was highly variable, and the reliability of many figures was dubious (e.g. precision only to the first 511 decimal place on most occasions). Thus, we did not consider this information adequate for model calibration.
We calibrated the 8 parameter-model (Table 1) using TP data collected from the  Roda de Ter sampling point from 1999 to 2002. TP data for the period 2003-2004 were left for the validation check and not used during calibration. However, since river 5 discharge used during calibration was a modeled variable, we corrected the possible effects of errors in discharge simulation on modeled TP values. TP concentration in the river at Roda de Ter followed a power dilution dynamics with discharge (TP=0.35×Discharge −0.36 , p<0.0001, n=106, r 2 =0.45). Therefore, any mismatch between observed and modeled discharge will have a profound effect on the calibration 10 process, especially at low discharges. To solve this problem, we performed calibration on a corrected TP observed series, using TP c =TP TP mod TP obs (6) where TP c is the corrected TP observed value. TP mod and TP obs are the TP values predicted by the above power regression using the modeled and the observed discharge, 15 respectively (Fig. 3A). The correcting quotient in Eq. (6) averaged 1.09 for all TP data used during calibration. Calibration was automatically done using the Shuffled Complex Evolution algorithm (SCE-UA), which was developed to deal with highly non-linear problems (Duan et al., 1992). From an initial population of randomly generated parameters, the algorithm 20 uses shuffling, competitive evolution, and random search to efficiently find the parameter set that minimizes an objective function (OF). In this case, the OF was the sum of the squared errors between model outcomes and corresponding TP c values. We performed the calibration run using SCE-UA as implemented in the PEST package (Doherty, 2003), with parameter bounds detailed in Table 1. 25 512 2.6 Model structure coherence In order to assess whether the final model structure was realistic, we compared the adjusted values of the nutrient spiraling metrics in the HSPF model with values from field-based research performed in the watershed under study and in other systems worldwide. The comparison with metrics measured in the Ter watershed was difficult, 5 because published field estimations of Nutrient Spiralling metrics from the Ter watershed mostly report data for pristine streams (Martí and Sabater, 1996;Butturini and Sabater, 1998), while the calibration of the HSPF model is based on data collected downstream a highly human impacted area. Thus, comparing retention metrics from these studies with the fitted metrics in our model could be misleading. Fortunately, 10 Martí et al. (2004) reported v f for two phosphorus retention experiments in a reach in the impaired Riera de Tona (Gurri River tributary, Fig. 1B), a location close to our sampling TP point.
We could take the comparison between modeled retention metrics and field-based estimations a step further. During recent years, researchers have accumulated data 15 that suggest nutrient enriched streams have lower retention efficiency (i.e., lower v f or higher S w ) than pristine streams (Doyle et al., 2003;Martí et al., 2004;Haggard et al., 2005;Merseburger et al., 2005;Gücker and Pusch, 2006;Ruggiero et al., 2006). To test how our model results fit into this picture, we collected S w results for phosphorus (for many studies v f results were not available) from pristine and nutrient enriched 20 streams. If the fitted S w in our model is a realistic approximation of the real value, it must resemble S w values measured in impaired streams. Note that collected results come from very heterogeneous field procedures (nutrient additions, nutrient decay downstream from a point source, isotopic tracers), and that they lump seasonal studies with one-measure data, and habitat specific experiments with whole stream 25 determinations. The most important implication is that while S w for pristine streams is usually assessed with nutrient enrichment experiments, thus reporting gross retention (Martí et al., 1997), most data from impaired streams comes from ambient nutrient de-513 cay experiments, which must be considered reporting net retention metrics. Obviously, our model estimates for the Ter watershed should be considered as a net retention. Finally, values from the literature are based on dissolved inorganic phosphorus retention while our model predicts TP. Although this could introduce some bias in the analysis, the low proportion of particulate phosphorus in this human impacted stream suggests 5 that the comparison between our results and the bibliographical values is acceptable.

Results
During HSPF calibration with SCE-UA, convergence to an optimized parameter set (see Table 1) was achieved after 7000 model runs. Factors for point source correction (C i and C w ) were adjusted to values different than one, suggesting that the available 10 database for point sources had significant biases. The TP load from WWTP seemed to be overestimated in the database, while the industrial spills were slightly underestimated. Applying C w and C i for the mean annual TP loads we obtained 19 000 kg P yr −1 from WWTP and 12 300 kg P yr −1 from industrial spills. Considering the diffuse TP inputs, the power function fitted for groundwater TP concentration had a very gentle 15 slope (b g , Table 1), implying that TP g was nearly a constant value in the range of Q g modeled in the Ter watershed (TP g around 0.06 mg P L −1 ). By contrast, the slope for the power relationship between TP i and Q i defined a clear dilution dynamics, with TP i concentration ranging from 0.6 to 0.04 mg P L −1 depending on Q i values. Using these power relationships with the time series of Q i and Q g we obtained mean annual TP 20 loads of 23 600 kg P yr −1 from groundwater discharge and 12 800 kg P yr −1 from interflow discharge. The mass transfer coefficient v f was optimized to a very low value (Table 1). On the other hand, the temperature correction factor (TC, Table 1) was adjusted to 1.06. Considering that mean daily river water temperature in the watershed ranges from 5 to 25 27 • C (Fig. 3), this implies that v f values were multiplied by a factor (Eq. 4) that ranged 514 from 0.4 to 1.3. Thus, actual v f values after temperature correction ranged between 5.6×10 −7 and 1.8×10 −6 m s −1 .
The fit between observed data and model outcomes at Roda de Ter was satisfactory (Fig. 4). The model explained 72% of variance in river TP c values during the calibration period (the contribution of the very high value during year 2000 was modest. 5 Without this point the explained variance amounted 69%). However, the model performed worse during high flow conditions (or low TP concentrations), as Fig. 5 clearly shows. This was most evident during the validation period, a very wet period (Fig. 3). In addition, the fit between median TP values coming from ACA stations and model results was good (Fig. 6), although ACA station 7 showed observed values that were 10 considerably higher than model outcomes.
From results found in the literature (Table 2), a clear power relationship could be established between S w values and discharge (Fig. 7). This relationship could be split differentiating pristine streams (1622 Q 0.65 , n=44, p<0.0001, r 2 =0.56) and data coming from nutrient-enriched streams (13 163 Q 0.51 , n=20, p<0.0097, r 2 =0.32). The 15 power relationship obtained by transforming the adjusted v f value for the Ter watershed to S w with Eq. (5) is the bold dotted line in Fig. 7, and corresponds to the equation 24 742 Q 0.77 . Note that for comparisons between the different power regressions, the adequate parameter is the intercept of the power regression, because the slope will depend on the geomorphologic traits of the rivers included in each relationship 20 (Stream Solute Workshop, 1990). Bearing this in mind, power regressions for the Ter River watershed and for impaired streams were similar, especially if we reevaluate the power regression for impaired streams discarding points labeled as j, r, and n in Fig. 7 (21 256 Q 0.49 , n=17, p<0.0001, r 2 =0.73, bold line in Fig. 7). The presence of these points, which represent very short phosphorus S w in nutrient enriched streams, should 25 be attributed to methodological constraints. Most of the nutrient retention experiments in impaired streams were measuring net retention. Since in impaired streams point sources and diffuse inputs can be inextricably linked (Merseburger et al., 2005), it is not easy to assign this low S w to the effect of actual in-stream processes or to lateral 515 inflows of nutrients by seepage.

Discussion
The low mass transfer coefficient v f optimized in our model is only comparable with values obtained in point-source impaired streams (Doyle et al., 2003;Martí et al., 2004). Values from pristine streams usually fall between 10 −3 and 10 −5 m s −1 (Doyle et al., 5 2003). Our low v f defines a watershed with watercourses with very low phosphorus retention capacity. Of course, this would probably hold in reaches around the sampling point at Roda de Ter, while in headwater streams the value will probably be underestimated. Thus, we must take this v f figure as a coarse-scale value. On the other hand, the significant dependence on water temperature suggested that v f for TP in this 10 watershed is controlled to some extent by biological activity. However, as an empirical correction factor, this could also reflect any seasonal process related to TP retention showing covariance with stream temperature. Thus, results from this study cannot be used to state that temperature is modulating TP retention. Concerning the model fit, it seemed that the model was missing some significant 15 effect at high flows, which could be attributed to physically-mediated higher retention during high flows not accounted for in our formulation, or to an overestimation of TP g during very wet periods. The particulate fraction of TP could play a role in these misfit situations, but data from this study did not allow an accurate assessment of this possibility. Low TP values modeled for ACA station 7 should be attributed to a missing 20 point source in the database upstream from this sampling point, considering that the adjusted v f value for the watershed represented a very low retention efficiency. However, despite these shortcomings, results from this study showed that the formulation on which the Nutrient Spiralling concept research is based is a good alternative for modelling the nutrient in-stream processes in a watershed-scale model. Even considering that we worked in a worst case scenario, in the sense that limited river TP concentration data were available to calibrate the model, model outcomes were satisfactory and adjusted parameter values realistic. These results pose the following 516 question: can we use field estimations of Nutrient Spiralling metrics to feed our model? Of course, the best method to test this possibility would be to measure S w (and then calculate v f with Eq. 5) in several reaches in the Ter watershed, and then compare this with our estimate. But this is beyond the scope of this work. However, the mean v f for two nutrient retention experiments in a reach in the impaired Riera de Tona (Gurri 5 River tributary, Fig. 1B) was 4.6×10 −6 m s −1 (Martí et al., 2004), which is an astonishingly similar figure compared to our adjusted reference value (Table 1). In fact, using Martí et al.'s empirical value in our model only caused a slight deviation in the model results (66% of TP explained variance compared to 72% with the optimized parameter). A more general test of the adequacy of the model structure is the comparison with 10 S w vs. streamflow power regressions based on data coming from impaired streams of the world. The dependence of S w on streamflow was already reported for phosphorus (Butturini and Sabater, 1998) and ammonia retention (Peterson et al., 2001) in pristine streams. Our fitted power relationship between S w and discharge in pristine streams slightly differed from the equation reported by Butturini and Sabater (1998), because 15 our database includes recent data. However, the most interesting fact in Fig. 7 was that a significant power relationship was also fitted with data coming from nutrient-enriched streams. Thus, the lack of relationship between phosphorus S w and discharge reported in impaired streams (Martí et al., 2004) can be attributed to a narrower discharge range in previous studies. In fact, the relationship between S w and discharge is highly plau-20 sible considering Eq. (5) (Stream Solute Workshop, 1990). The resemblance between the power relationship obtained by transforming the adjusted v f value for the Ter watershed and the obtained for impaired streams is notable, and suggest that the model structure used in our model is adequate and realistic. Results from Fig. 7 could be interpreted in two ways. First, retention efficiency greatly result for a nutrient enriched system. The second interpretation is as follows: due to this coincidence we can use data from empirical studies of nutrient retention to parameterize a watershed-scale model. Obviously, this is circular reasoning, and the model results cannot be used to state that we demonstrated what the above interpretations imply. But this is a nice example of how models can work as heuristic tools to compare 5 hypotheses and stimulate research (see Oreskes et al., 1994).

Conclusions
To conclude, we have demonstrated that a lumped, hardly parameterized formulation of the in-stream nutrient fate in rivers could be very efficient in a large-scale model, and that this opens the very interesting possibility of directly using data collected in the field 10 in large-scale applications. This avoids the exercise of upscaling fine-scale research results to parameterize do-everything models with many parameters, many of them finally adjusted to bibliographical values on most occasions. This reasoning should apply for any large scale model, in the sense that formulation of lumped processes must be prioritized, especially when information equivalent to those lumped processes can 15 be obtained in the field to directly parameterize the model. Of course, this is not a valid option if the detailed biogeochemical processes are research targets, or if we need explicit formulations of these processes to simulate complex biotic or abiotic interactions. However, the coarse-formulation approach should suffice in many modeling exercises. 130,[215][216][217][218][219][220][221][222][223][224][225][226][227]1994. 525 Marcé, R., Comerma, M., García, J. C., and Armengol, J.: A neuro-fuzzy modelling tool to estimate fluvial nutrient loads in watersheds under time-varying human impact, Limnol. 520 Intercept for TP vs. interflow discharge (Eq. 1) mg P L −1 3.5×10 −5 -0.38 0.002 b g Slope for TP vs. groundwater discharge (Eq. 1) mm −1 0-1.8 0.026 a g Intercept for TP vs. groundwater discharge (Eq. 1) mg P L −1 3.5×10 −5 -0.38 0.05 Point-sources correction C w Correction factor for TP load fom WWTP's -0-9 0.63 C i Correction factor for TP load from industrial spills -0-9 1.16 524 Pioneer Creek (USA) 0.0856 370 Davis and Minshall (1999) 5 Bear Brook (USA) 0.0145    Table 2. See the text for details on power regressions.