Estimation of phosphorus export from a Mediterranean agricultural catchment with scarce data

Abstract

Knowledge of the amount of nutrient input to receiving waters by their catchments is often mandatory for water managers. In this study, we estimated the export of phosphorus from a Mediterranean catchment into a downstream reservoir under the constraint of limited sampling data availability. We calibrated the physically- based catchment-scale model ZIN-Sed 2D to a regionalized flow duration curve and scarce event mean sampling data of dissolved and particulate phosphorus. The model results were further confirmed by discussion of soft data from the literature. Our modelling approach differs from commonly applied models by using a time step of less than 1 day, and a Langmuir isotherm for phosphorus sorption instead of the linear isotherm. The successful model application showed that this combination of methods is applicable under data-scarce conditions and the long-term phosphorus export from the study site is below average for this region.

Editor D. Koutsoyiannis; Associate editor T. Okruszko

Citation Gassmann, M., Brito, D., and Olsson, O., 2013. Estimation of phosphorus export from a Mediterranean agricultural catchment with scarce data. Hydrological Sciences Journal, 59 (1), 221–233.

Résumé

La connaissance des quantités de nutriments apportés aux eaux réceptrices par leurs bassins versants est souvent obligatoire pour les gestionnaires de l’eau. Dans cette étude, nous avons évalué l’exportation de phosphore d’un bassin méditerranéen vers un réservoir situé en aval sous la contrainte de disponibilité limitée des données d’échantillonnage. Nous avons calé le modéle à base physique ZIN-Sed 2D à l’échelle du bassin versant à une courbe régionalisée des débits classés et à de rares données moyennes événementielles d’échantillonnage de phosphore particulaire et dissous. Les résultats du modéle ont de plus été confirmés par l’examen de données qualitatives issues de la littérature. Notre approche de modélisation différe des modéles communément appliqués, en utilisant un pas de temps infra-journalier, et une isotherme de Langmuir pour la sorption du phosphore au lieu de l’isotherme linéaire. L’application réussie du modéle a montré que cette combinaison de méthodes est applicable en condition de données limitées., L’exportation à long terme du phosphore du site d’étude se situe en dessous de la moyenne pour cette région.

Key words:

• phosphorus export modelling
• ZIN-Sed 2D model
• synthetic flow duration curve
• Langmuir isotherm
• water quality

Mots clefs:

• modélisation de l’exportation de phosphore
• modéle ZIN-Sed 2D
• courbe synthétique des débits classés
• isotherme de Langmuir
• qualité de l’eau

1 Introduction

The water quality of receiving waters such as lakes, reservoirs and lagoons is highly dependent on the quality of their recharging water (e.g. Velasco et al. 2006, Olsson et al. 2013). Thus, the solution to a water-quality problem of those receiving waters is often to be found in their catchments (Schindler 2006). In particular, non-point (diffuse) contaminant sources from agricultural areas account for the input of large amounts of nutrients into surface waters (e.g. Vanni et al. 2001). The input of phosphorus (P), for example, is critical as it is a minimum factor for primary production in the aquatic environment and thus a driver for algal blooms (Schindler 1974).

The most straightforward method for the determination of phosphorus export is direct sampling. However, not all catchments are well gauged and thus the long-term load of P and water being transferred from the catchment to the receiving water is largely unknown. A method to tackle the problem of ungauged catchments is the use of models for estimating P export. Generally, there are three types of these models in the literature: (a) Process-based models consider known P cycling processes and are usually linked to a hydrological model (Borah and Bera 2003). (b) Export coefficient models use coefficients of P release per area and time step for each land-use class to assess the total load at the catchment outlet with different complexity (Johnes 1996). (c) Empirical models try to establish relationships between sampled P concentrations and catchment variables like slope, catchment size or fraction of agricultural fields (e.g. Ekholm et al. 2000). While all three model types have been applied successfully in different catchments, only process-based models are applicable in ungauged catchments due to the amount of P sampling data required to establish other types of model (Sivapalan et al. 2003). However, even the best process-based models are not helpful without any information for calibration. One possibility is the calibration of the process-based model in a nearby catchment with similar characteristics, and its subsequent application in the catchment of interest (Parajuli et al. 2009). Another successfully applied approach is the regionalization of flow duration curves (FDCs) using many surrounding catchments (Yu and Yang 2000). If it is possible to obtain water samples in the ungauged basin, they may improve modelling results enormously, even if few in number (Seibert and Beven 2009). As shown by Seibert and McDonnell (2002) the use of soft data may help to lower the uncertainty of modelling results. Soft data are rather qualitative information about a certain environmental behaviour which cannot be expressed in exact numbers. Examples of soft data are a perceptual model for catchment behaviour found by experimental work, or the reasonability of calibrated parameter values. In particular, the conceptual correctness of a model can be evaluated by rejecting model set-ups which produce results contradicting an experimental perceptual model. In contrast to soft data, hard data contain information in the form of exact numbers, such as river discharge sampling.

There are many process-based dynamic models widely used in literature for P export from catchments: Annualized Agricultural Non-Point Source (AnnAGNPS; Bingner et al. 2011), Areal Nonpoint Source Watershed Environment Response Simulation (ANSWERS-2000; Bouraoui and Dillaha 2000), Groundwater Loading Effects of Agricultural Management Systems (GLEAMS; Leonard et al. 1987), Generalized Watershed Loading Function (AVGWLF; Evans et al. 2008), Hydrological Simulation Program Fortran (HSPF; Donigian et al. 1984), Soil and Water Assessment Tool (SWAT; Arnold et al. 1998) and Integrated Catchments model of Phosphorus dynamics (INCA-P; Wade et al. 2002). However, except for HSPF, all these models only run at a daily time step, but rainfall–runoff events in small catchments are usually shorter than 1 day. If only a few event samples are available and, thus, no appropriate daily mean concentration can be calculated, a daily modelling time step may overestimate or underestimate the model results.

The representation of P sorption is a central question for process-based modelling, as P is strongly adsorbed to soil particles (e.g. Shields et al. 2009). Smil (2000), for example, pointed out that more than 70% of the total phosphorus (TP) transported by rivers worldwide is bound to soil particles. A closer look to the equations behind the above-mentioned models reveals that they all use a linear sorption approach in their current version. However, many studies have found nonlinear sorption behaviour for P on soil particles. The Langmuir isotherm could often be fitted successfully to experimental data (e.g. House et al. 1995). In contrast to the linear isotherm, the Langmuir relationship has a maximum sorption capacity and thus more P remains in solution at high concentrations.

This paper presents a method for analysing the long-term export of phosphorus from Mediterranean catchments that is applicable under data-scarce conditions. It introduces the combination of using a method for discharge prediction in ungauged basins, scarce sampling data and a physically-based hydrologic model (ZIN Sed 2D) containing a nonlinear Langmuir isotherm for calculations of P sorption. In this manner, the short P release events in the river were simulated at sub-daily time steps. The results are discussed in the light of concepts and data from the literature (soft data validation). The Mediterranean catchment of the intermittent stream Enxoe, in Portugal, which discharges into a downstream reservoir that has suffered from cyanobacteria blooms in the past (Galvão et al. 2008) was chosen for the case study.

2 MATERIALS AND METHODS

2.1 Catchment and data availability

The Enxoe catchment is located in southeastern Portugal. It has an area of 61 km² and drains down to the Enxoe Reservoir which was designed as storage for drinking water supply (Fig. 1). However, the Enxoe Reservoir has suffered from cyanobacteria blooms in the last decade and thus its water could only be used for agricultural purposes during the toxic algal blooms. It was hypothesized, that the input of phosphorus (dissolved and particulate) from the agricultural area of the catchment stimulates the primary production in the reservoir and thus is responsible for algal bloom development (Coelho et al. 2008). In this study, only a part of the catchment down to the sampling station Enxoe Main (37 km²) is considered.

Fig. 1 The Enxoe catchment and location of the sampling sites.

The area has a Mediterranean climate with hot, dry summers and mild wet winters. Annual rainfall is about 493 mm (2001–2011) and average annual air temperature around 16.3°C. The distribution of soil types in the catchment based on grain-size is clay (27%) and (silty and sandy) clay loam (73%). Land use is dominated by extensive agriculture (95%) and forest (4%). The main agricultural crop is olive groves (50%). Only 1% of the catchment land area is urban. Although there is a wastewater treatment plant (WWTP) in the catchment, no significant anthropogenic influence is expected on water quality since the discharge of the WWTP is directed to outside the catchment.

In 2010–2011, automatic sampling for phosphorus and suspended sediment concentration (SSC) was conducted in the Enxoe catchment during eight flood events at Enxoe Main station (Fig. 1). For this purpose, an automatic water sampler (EcoTech Umwelt-Meßsysteme GmbH Bonn, Germany), with eight bottles (2 L), was installed. The sampler was located at a river reach where homogeneity of water movement was considered representative during all hydrological conditions. The beginning of sampling was triggered by a water level sensor as soon as the level changed by 10 cm on both the rising and falling limb of a flood wave. In addition, grab samples (2 L) were taken weekly at Vale do Vargo and Enxoe Tributary.

Water samples with volumes ranging from 250 to 1000 ml, depending on the suspended sediment load, were filtered in the laboratory to determine SSC using pre-weighed glass microfibre paper (GFF 0.75 µm). The sediments retained on the filter paper were oven dried at 50°C for 24 h. The filters were again weighed and SSC was calculated. Before analysing for total dissolved phosphorous (DP), each sample was filtered using a cellulose acetate membrane (0.45 µm). The TP was determined using unfiltered samples; TP and DP were quantified colorimetrically after sulphuric acid and nitric acid digestion by reacting with ammonium molybdate (Eaton et al. 1995). Particulate phosphorus (PP) was calculated as the difference between TP and DP concentrations.

Peak TP concentrations in the river were as high as 1–5 mg L^-1 with an event mean of 64% of PP and 36% of DP. However, no discharge measurements were conducted. Thus, for comparison of model results, sampled discharge values from surrounding catchments were used as described in the following sections. All discharge time series were obtained from the Sistema Nacional de Informação de Recursos Hídricos (http://www.snirh.pt), a freely accessible database of meteorological and hydrological data in Portugal. This database also provided hourly rainfall data for four stations surrounding the Enxoe catchment for the modelling in this study. Table 1 shows the time span of the available discharge and rainfall data of all neighbouring catchments to Enxoe. It also shows the long-term water balance of the catchments.

Table 1 Catchments neighbouring the Enxoe catchment, their available daily discharge data time span and water balance. d_Enxoe: distance to Enxoe catchment; P: precipitation; Q: mean discharge; C_Q: runoff fraction; ET_a: calculated storage change, including evapotranspiration; C_ET: ET_a as a fraction of Q. Values in brackets are standard deviations.

Catchment	Area (km^2)	dEnxoe (km)	Discharge data time span	P (mm)	Q (mm)	CQ (%)	ETa (mm)	CET (%)
Entradas	51	58	1971–2000	489 (162)	116 (134)	24	373 (85)	76
Limas	249	25	1982–1990	519 (202)	89 (101)	17	430 (114)	83
Monte da Arregota	97	38	1981–1990	529 (170)	95 (126)	18	434 (85)	82
Monte da Ponte	618	42	1971–1997	482 (163)	127 (135)	26	355 (74)	74

2.2 ZIN-Sed 2D model

The Zin-Sed 2D model is a further development of the ZIN-Sed model (Gassmann et al. 2012). In a fully distributed way, it is able to simulate hydrology, erosion/sediment transport and phosphorus release at the catchment scale. After erosion, sediment particles follow water flow until they either leave the catchment, or are deposited on the land surface or in the stream. Phosphorus is simulated as a dissolved and a sorbed substance. While PP follows sediment movement, DP is transported wherever water flows. The source and sink of DP is sorption to/from soil particles and fertilizer input. In this study, the model runs at an hourly time step in order to catch peaks of short events to be able to calculate representative event mean values.

2.2.1 Hydrological modules

Water infiltration through the soil surface is calculated by the Green and Ampt approach and (vertical and lateral) water movement in three soil layers by the Buckingham–Darcy equation. Unsaturated hydraulic conductivity is determined using the Mualem–van Genuchten equation (van Genuchten 1980). Additionally, an anisotropy factor is implemented which multiplies the hydraulic conductivity by a factor in the lateral direction. Water draining vertically out of the third soil layer is considered as loss to the system. The boundary condition for vertical drainage of water is the saturated hydraulic conductivity of the underlying bedrock. Overland flow is calculated from cell to cell in eight possible directions by a diffusive wave routing scheme solved by the Manning flow equation.

The Food and Agriculture Organization of the United Nations (FAO) crop method (Allen et al. 1998) was used in this study to generate potential evapotranspiration (ET) values. The ET is multiplied in the model by the relative filling of a soil layer between permanent wilting point (PWP) and field capacity (FK) in order to obtain actual ET (ET_a) values. The PWP and FK were determined as specific pressure heads of the water retention curve as given by van Genuchten (1980). The ET from soil is controlled by the rooting depth and depth of the soil layer from surface. The channel routing of ZIN-Sed 2D uses an explicit kinematic wave routing scheme as used, for example, in CASC2D (Rojas et al. 2003).

2.2.2 Erosion and sediment transport

Three erosion processes are modelled separately in ZIN-Sed 2D: interrill erosion in the field, rill erosion in the field and bed erosion in the river channel. Interrill or sheet erosion is calculated by a relationship between rainfall intensity and erosion (Hairsine and Rose 1991) as explained in Gassmann et al. (2012). Rill erosion follows a critical shear stress approach. For water flowing from cell to cell, shear stress is calculated from water density ρ, water velocity and the drag coefficient as (Schlichting 1979):

(1)

where C_d can be related to Manning’s roughness value, n, the depth of flowing water, h, and the gravity acceleration, g, by (DHI 2011):

(2)

As soon as the critical shear stress, τ_c,erosion is exceeded, erosion, E, starts with (Partheniades 1962):

(3)

The amount of erosion can be controlled by the erosion coefficient EC (g m⁻² s⁻¹). For the application to the field, equation (3) is multiplied by the Universal Soil Loss Equation (USLE) crop type factor, C_e (Wischmeier and Smith 1978).

Deposition, D, of suspended sediment occurs if the actual shear stress is below the critical value for deposition, τ_c,depo, influenced by suspended sediment concentration, , and particle fall velocity, (Krone 1962):

(4)

Maximum sediment transport capacity of flowing water is calculated by the approach of Govers (1990). Erosion and deposition in the channel are also solved by equations (1)–(4).

2.2.3 Phosphorus model

As a starting point for phosphorus calculations in ZIN-Sed 2D, the initial soil phosphorus concentration (mg kg^-1) has to be defined in the model. The soil phosphorus content in the soil is affected by sorption, drainage of the soil water and by fertilization. Erosion changes the total amount of P in soil, but not its soil concentration. Equations for uptake of phosphorus by plants, as well as transformation of organic P to mineral P and the degree-day plant growth model, were implemented in ZIN-Sed 2D as found in the SWAT model (Neitsch et al. 2005) and calculated daily.

The partitioning between dissolved (mg L^-1) and sorbed (mg g^-1) phosphorus is usually governed by an isotherm, i.e. the relationship between both. The quotient between sorbed and dissolved material is called the partitioning coefficient. As isotherms are generally nonlinear, the actual partitioning coefficient is calculated as the first derivative of the isotherm equation or the slope of the isotherm at a specific point. In this study, the Langmuir isotherm with the Langmuir constant K_L and the maximum sorption capacity c_max was used:

(5)

The sorbed equilibrium concentration τ_c,sorbed can be determined analytically if the total concentration, c_total, the suspended sediment concentration, SSC, and the dissolved equilibrium concentration, τ_c,solved, are known:

(6)

with

(7)

In order to describe the temporal delay of the sorption/desorption process (retardation), a pseudo first-order rate equation was used (Azizian 2004):

(8)

where c_e,sorbed is the sorbed concentration at equilibrium; r is the rate constant and dt is the time step. Adsorption and desorption velocities usually are different. Thus, the rate constant r generally has a different value for adsorption (r_ad) and desorption (r_de).

Compared to bulk soil, eroded sediment in overland flow is enriched with P (Sharpley 1995). This effect was explained by the fact that fine sediments are more readily eroded than coarse sediments and contain a larger maximum P transport capacity due to the larger surface. With increasing erosion, fine and coarse particles are eroded in the same amounts and thus the ratio between eroded sediment P content and bulk soil P content, the enrichment ratio (PER), decreases. This is expressed in the widely-used relationship between PER and the erosion amount SSY (CREAMS model, Menzel 1980), which was also used in the ZIN-Sed 2D model:

(9)

2.3 Parameterization and calibration

2.3.1 Hydrological calibration data

Since no discharge sampling data were available for the Enxoe catchment, a discharge regionalization approach developed for prediction of FDCs in ungauged basins was used. This approach consists of two steps:

1. a synthetic flow duration (SFDC) curve is derived from flow data of neighbouring catchments (Yu and Yang 2000), and

2. the hydrological model is calibrated using the SFDC.

The SFDC is determined by relating the catchment areas of surrounding catchments to quantiles of their FDC. Then, a linear regression (forced intercept 0.0) is fitted per quantile of discharge value dependent on catchment area. These regression equations are further used in combination with the area of the ungauged catchment in order to derive its FDC.

FDCs contain no information about the temporal distribution of discharge. Thus, we further compared seasonal discharge regimes of the neighbouring catchments to the modelled discharge regime using the principles of Pardé (1963). Finally, the water balance of the surrounding catchments and the modelled values were compared.

2.3.2 Hydrological module

Since ZIN-Sed 2D is a distributed physically-based model, it contains many parameters. In order to decrease the number of parameters for calibration, aggregate parameters (Yang et al. 2007) were used in this study in combination with literature-derived parameter values. Soil hydraulic parameters were obtained from pedotransfer functions used in the software ROSETTA (Schaap et al. 2001). As input to ROSETTA, grain-size distributions were taken from the soil map of Portugal. The crop coefficient for ET calculation per land use was taken from Allen et al. (1998). Although the crop coefficient usually changes according to crop growth stage, we kept it constant in time and space and calibrated it later. We assumed this simplification was applicable, as two thirds of the catchment were either covered by olive groves or agroforestry, and thus their crop coefficient was relatively stable all year. Moreover, crops were grown in the catchment only during winter time when ET was lowest and thus the potential error in ET_a least. Manning roughness values were taken from the literature for each land use and not calibrated any further, because the time step of the hydrological data for calibration (1 day) was in the order of mean event length. Thus, changing the roughness values would not have affected hydrological modelling much. Due to the lack of exact values, soil depth was assessed using values from the Harmonized World Soil Database (HWSD; FAO/IIASA/ISRIC/ISSCAS/JRC 2009) and distributed between 0.5 and 2.0 m, decreasing the depth with slope. Using aggregate parameters as explained above, the hydrological module in this study has five parameters left for calibration, the erosion module seven parameters and the phosphorus module six parameters (Table 2).

Table 2 Parameters, their units and calibrated values of the ZIN-Sed 2D model in this study. P: phosphorus.

Parameter	Unit	Description	Value
Kf1_m	–	Multiplier of hydraulic conductivity of soil layer 1	3
Kf23_m	–	Multiplier of hydraulic conductivity of soil layers 2 and 3	2
soilD_m	–	Multiplier of soil depth	0.4
Klat_m	–	Anisotropy factor in soil	5
Crop_m	–	Multiplier of FAO crop coefficient	0.9
τc,erosion	N m^-2	Critical shear stress for erosion in overland flow	0.15
τc,depo	N m^-2	Critical shear stress for deposition in overland flow	0.10
EC_OF	g m^-2 s^-1	Erosion coefficient in overland flow	0.003
VS_OF	mm s^-1	Settling velocity of sediments in overland flow	0.1
TC_OF	–	Transport capacity factor in overland flow	0.2
TC_River	–	Transport capacity factor in river	0.2
Dm	–	Multiplier of rainfall detachability factor	0.1
KL	L mg^-1	Parameter K of Langmuir isotherm	1.3
cmax	mg g^-1	Maximum sorbed concentration	3.5
rad	d^-1	Rate of adsorption	20
rde	d^-1	Rate of desorption	20
SPm	–	Multiplier of soil P content	4
b	–	Exponent of initial soil P profile	–0.5

2.3.3 Water-quality modules

In order to calibrate the water-quality modules, we compared sampled event mean values to model results for SSC, DP, PP and TP. Additionally, the DP model results at the outlet were compared to DP river concentrations during baseflow sampled in two tributaries (Enxoe Tributary and Vale do Vargo, Fig. 1). In doing so, especially the sorption behaviour of the soil storage of a catchment and, thus, the parameterization of the sorption isotherm, could be evaluated. Because mixing is a linear process and significant inflows at this short distance were not expected during baseflow conditions, the river concentration below the confluence of the two tributaries was assumed to be between the concentrations of the tributaries.

The main boundary condition for P release modelling from agricultural areas is soil P content. As a starting point for the soil concentration, the same values were used as assumed in the SWAT model: the initial TP concentration was set to 42.5 mg kg^-1 in pristine areas and 212.5 mg kg^-1 in agriculturally used areas. These values were multiplied by an aggregate factor in order to fit sampled concentrations. The vertical distribution of phosphorus in soil was initially calculated as a strongly decreasing (power function) with soil depth. This function was normalized in order to contain the given initial phosphorus value in the uppermost layer. The depth distribution function was parameterized by changing the exponent (Table 2).

2.4 Sensitivity analysis

Sensitivity analysis shows the ability of model parameters to alter model output variables (here discharge, SSC, DP and PP). All parameters in Table 2 were tested for their sensitivity regarding the calculation of annual yields of discharge (QY), suspended sediment (SSY), DP (DPY) and PP (PPY). Since this model has long run times, we decided to vary each parameter by ±50% and compare the results with each other and the baseline scenario, which was the set-up of the calibrated model run. The model was run for the hydrological year 2001/02 with an additional 6 months of warm-up phase, at a time step of 1 hour. Each sensitivity run took about 90 minutes and in total 37 runs were performed. Each parameter j was analysed for its local sensitivity SI in the parameter range Δp_j via n model runs. The model output variation ∂y_i was set in relation to the parameter variation ∂p_i, as used in Frey et al. (2009):

(10)

Results were normalized in order to compare the variation of SI values of discharge, sediment and TP release with each other.

3 RESULTS

3.1 Sensitivity analysis

Parameter sensitivity was analysed for all calibration parameters (Table 2) in relation to annual discharge, sediment, DP and PP release (Fig. 2). All the hydrological parameters, except Klat_m, were to some extent also changing water-quality variables. The most sensitive parameter for QY was the multiplier of the FAO crop ET coefficient. The multiplier of soil depth was the most sensitive parameter for PPY and SSY. Generally, PPY and SSY changed in a similar way to the reaction to changes of hydrological and sediment parameters. This shows the strong connection of both variables in the ZIN-Sed 2D model. The DPY changed most due to changes in the Langmuir isotherm parameters. The multiplier of soil phosphorus was also highly sensitive for phosphorus release. Only Dm, sediment sink velocity (VS_OF) and the adsorption rate coefficient r_ad were not sensitive in the analysed parameter ranges.

Fig. 2 Results of the local sensitivity analysis. Parameter abbreviations are explained in Table 2. The ZIN-Sed 2D modules (hydrology, erosion, phosphorus) of the parameters are indicated above the diagram. SSC: suspended sediment; PP: particulate phosphorus; and DP: dissolved phosphorus.

3.2 Hydrological modelling

The relationship between flow quantiles and catchment areas of the neighbour catchments of Enxoe had a coefficient of determination (R²) of R² > 0.90 for all but the 50th quantile (R² = 0.65) (Table 3). These correlations allowed for an assessment of the unknown FDC of the Enxoe catchment. The subsequent calibration of the ZIN-Sed 2D model resulted in very close modelled and synthetic FDCs with a R² of 1.00 (n = 13) and a Nash–Sutcliffe efficiency (Nash and Sutcliffe 1970) of 0.98 (Fig. 3). The long-term water balance (2001–2011) modelled by ZIN-Sed 2D brought an average annual discharge of about 110 mm year^-1 and a storage change of 383 mm year^-1, of which 330 mm were ET. A validation of temporal distribution of discharge can be drawn from Fig. 4; the modelled seasonal discharge regime follows the general trend of the sampled regimes of the catchments neighbouring Enxoe. More specifically, the coefficient of determination of modelled seasonal regime compared to median sampled regime was 0.92 (p < 0.05, n = 4).

Table 3 Quantiles (Q_x at x percent), regression slope (s), coefficient of determination (R²) and discharge values (Q_SFCD) of the synthetic flow duration curve. Regressions were made with n = 4 catchments (Table 1).

	R^2	s	QSFCD (m^3 s^-1)
Q0.05	0.96	0.34608	12.86
Q0.1	0.95	0.29550	10.98
Q0.5	0.96	0.14317	5.321
Q1.0	0.95	0.09815	3.648
Q2.5	0.96	0.03576	1.329
Q5.0	0.98	0.01325	0.492
Q10	1.00	0.00384	0.143
Q15	0.98	0.00142	0.053
Q20	0.98	0.00084	0.031
Q30	0.99	0.00028	0.010
Q40	0.92	0.00009	0.003
Q50	0.65	0.00002	0.001

Fig. 3 Synthetic and modelled flow duration curve of the Enxoe River.

Fig. 4 Seasonal discharge (Q) regimes as sampled in catchments neighbouring the Enxoe (ranges) and modelled for the Enxoe catchment by ZIN-Sed 2D (circles).

3.3 Water-quality modelling

The performance of the ZIN-Sed 2D model during rainfall–runoff events after calibration can be seen in Fig. 5 for water-quality parameters. Generally, all quality parameters could be reproduced with medium, but significant, coefficients of determination and low values of root mean square error (RMSE) over the whole concentration range. An outlier could be observed for PP and TP modelling at the upper edge of the concentration range.

Fig. 5 Modelled and sampled mean event concentrations of suspended sediment (SSC), dissolved (DP), sorbed (PP) and total (TP) phosphorus. The grey, dashed line denotes the 1:1 line. RMSE: root mean square error; R²: coefficient of determination with n = 8 values.

Figure 6 shows a comparison of model results for DP to a range of sampled river concentrations of two tributaries just above their confluence (compare Fig. 1). Of the modelled DP values, 39% were located within these boundaries and the general baseflow behaviour of DP could be reproduced. A significant correlation with a coefficient of determination of 0.45 (p <0.0001, n = 36) and RMSE of 0.07 mg L^-1 was found between modelled DP concentrations and the Enxoe tributary DP samples. No significant correlation could be found to Vale do Vargo tributary DP values (R² = 0.0, p < 0.9, n = 36, RMSE = 0.23). Especially at the beginning of the dry season (June 2011) there was a rise in DP concentrations sampled at Vale do Vargo tributary which could not be reproduced by the model. Since discharge is very low at these times, the total mass is also low compared to event release of phosphorus.

Fig. 6 Modelled and sampled baseflow concentrations of dissolved phosphorus (DP). The sampled DP range is the range between sampled values of the two tributaries.

The calibrated model was further used for long-term assessment of phosphorus and sediment export from the Enxoe catchment towards the Enxoe Reservoir for hydrological years 2001–2011 (Fig. 7). Within this time span there were high-flow as well as low-flow seasons. The fraction of average annual DP in discharge varied from 39% to 55% of TP. An annual average of 0.27 kg ha^-1 year^-1 was calculated in this decade for DP and 0.31 kg ha^-1 year^-1 for PP export. The TP export values altered between 0.002 and 1.66 kg ha^-1 year^-1. The mean annual sediment release was 0.27 t ha^-1 year^-1 (range 0.0–0.8 t ha^-1 year^-1) and the mean discharge 4.08 hm³ year^-1 (range 0.01–13.2 hm³ year^-1). In all years, 90% of annual TP was released during 8–17 (non-consecutive) days except in the dry year 2004/05 (90% TP during 2 days). The mean annual concentrations of TP were relatively stable ranging from 0.44 to 0.68 mg L^-1.

Fig. 7 Long-term assessment of mean annual suspended sediment (SSY), dissolved (DP) and particulate (PP) phosphorus release, discharge and total phosphorus (TP) concentrations from the Enxoe catchment.

4 DISCUSSION

The low data availability and thus lack of hard model validation data raise the question as to what extent the model results actually reflect the environmental situation in the Enxoe catchment. In the following, we discuss our results with respect to previously published work (soft data) in order to generate more reliability in our modelling results.

The functioning of the hydrological modules regarding calibration to the SFDC is a method that was proven before (Yu and Yang 2000) and successfully applied here. An additional confirmation of this was the fact that the long-term average water balance components (discharge and storage change) of the model results were in the same range as the water balance of neighbouring catchments. The discharge coefficient was found to be 18–26% in this region (Table 1). The model predictions gave a value of 22% for the Enxoe catchment. Furthermore, the modelled seasonal discharge regime is close to the sampled regimes of nearby catchments. Thus, the modelled total water volume per year, modelled temporal discharge distribution and the modelled probability of sampling a certain discharge value (FDC) were consistent and confirm the model results.

Correlations between sampled and modelled values of suspended sediment and phosphorus components were of medium range but significant (Fig. 5). Since calibration events covered the whole spectrum from low- to high-concentration events and throughout the year (except the dry summer months) there is a certain internal stability in ZIN-Sed 2D model predictions over a wide concentration range (0.3–2.2 mg TP/L) and under different hydro-meteorological conditions. Additionally, the general baseflow behaviour of DP could be reproduced by the model for the Enxoe tributary that was only influenced by agriculture. In the Vale do Vargo tributary, especially the rising P concentrations at the beginning of summer could not be reproduced. This was most likely due to small uncontrolled P discharges coming from the small settlement of Vale do Vargo. We assume that these inputs were small compared to the TP mass exported from the catchment but resulted in high concentrations when discharge was very low at the beginning of summer (Fig. 6). The overall successful P modelling confirmed that the Langmuir isotherm, mostly derived by small-scale laboratory investigations, can successfully be used at the catchment scale. The same result was recently also found by Rossi et al. (2012) applying a modified SWAT model to a small sub-catchment (0.07 km²).

Sensitivity analysis showed that parameters in the hydrological module were more sensitive to sediment release than the parameters in the sediment module. The same is true for PP and the phosphorus module. This underlines the importance of the underlying hydrological model for the calculation of substance release, which agrees with the findings of Cryer and Havens (1999). During a model sensitivity analysis in a Spanish catchment, Caille et al. (2012) recently found that topsoil P content was the most sensitive parameter for agricultural P release. Although using a different model, the high sensitivity of this parameter was confirmed in this study. Specifically, a variation of soil phosphorus content of 50% resulted in a 50–60% change of TP output values in this study. Torrent et al. (2007) underlined the importance of single events for TP loss in the Mediterranean. The analysis of our model results showed that, in all years, 90% of TP was released in less than 17 non-consecutive days. This testifies to the conceptual correctness of the phosphorus module parameterization of our model set-up.

Soil erosion and thus sediment export on the Iberian Peninsula is very heterogeneous (Cantón et al. 2011). Using a coarse resolution (1 km²), annual soil erosion by water was previously calculated in this region to be between 0 and 1 t ha^-1 year^-1 (Kirkby et al. 2008). Considering the gentle slopes of the Enxoe catchment, our calculated average sediment export of 0.3 t ha^-1 year^-1 was a reasonable value and confirms the consistency of the model set-up for erosion.

Recently, the average TP release of hydrological catchments into the Mediterranean Sea was calculated by Grizzetti et al. (2012) to be around 0.7 kg ha^-1 year^-1. This value was confirmed by model results in a large catchment in Spain (Caille et al. 2012) and sampling in a coastal catchment in France (Obermann et al. 2007). However, these findings included urban anthropogenic sources of phosphorus, such as industry and wastewater, which are not present in the Enxoe catchment. At a smaller scale, Torrent et al. (2007) sampled a mean annual release of about 1 kg ha^-1 year^-1. In other European agricultural catchments of similar size, a release of 0.1–4.7 kg ha^-1 year^-1 was found (Kronvang et al. 2007). In a study including 685 catchments worldwide, by far the majority of catchments had annual exports 0–1 kg ha^-1 year^-1, but values of up to 50 kg ha^-1 year^-1 were reported (Alvarez-Cobelas et al. 2009). Thus, the estimated long-term TP release of this study (0.58 kg ha^-1 year^-1) can be classified as below average for this region, and much below the highest worldwide values.

5 CONCLUSIONS

This study has examined the estimation of long-term phosphorus export in a case study of the Enxoe catchment. We calibrated the physically-based catchment model ZIN Sed 2D to a regionalized FDC and scarce water-quality sampling data. Modelling results were further confirmed by soft data obtained from the literature. Thus, we could show that it is possible to successfully combine a method for prediction in ungauged basins, sampling data and soft data for model calibration.

The availability of sampling data is often restricted to short periods and discontinuous time series under data scarcity. With the sub-daily time step used in this study, we were able to compare model results to short periods with available sampling data, which would have been difficult with a daily time step. These results let us suggest that the short time step is especially appropriate for data-scarce conditions. The results of this study also imply that the applied Langmuir isotherm for phosphorus sorption calculation is an alternative to using the linear isotherm at the catchment scale.

In the Enxoe Reservoir, eutrophication was a major problem in the last decade. The model results of this study showed that long-term average TP export from the Enxoe catchment was below average for this region, and low relative to the range worldwide. We suggest further research regarding phosphorus cycling processes in the Enxoe Reservoir using our results as boundary conditions for assessment of the eutrophication potential of TP discharged from the catchment.

Data scarcity is a problem that scientists and water managers face in many catchments worldwide. In particular, if receiving water bodies suffer from eutrophication, knowledge about nutrient input from the catchment is mandatory. The information in this paper brings new insights to possible options to tackle the problem of missing discharge data and scarce water-quality sampling data in Mediterranean catchments. It shows that very limited data can be sufficiently analysed with the proposed methodology.

Funding

This study was funded by the European Union through the Mediterranean Intermittent River Management (MIRAGE) project (ref: FP7 ENV 2007–1).