Modelling the hydrological behaviour of a coffee agroforestry basin in Costa Rica

The profitability of hydropower in Costa Rica is affected by soil erosion and sedimentation in dam reservoirs, which are in turn influenced by land use, infiltration and aquifer interactions with surface water. In order to foster the provision and payment for Hydrological Environmental Services (HES), a quantitative assessment of the impact of specific land uses on the functioning of drainage-basins is required. The present paper aims to study the water balance partitioning in a volcanic coffee agroforestry microbasin (1 km2, steep slopes) in Costa Rica, as a first step towards evaluating sediment or contaminant loads. The main hydrological processes were monitored during one year, using flume, eddy-covariance flux tower, soil water profiles and piezometers. A new Hydro-SVAT lumped model is proposed, that balances SVAT (Soil Vegetation Atmosphere Transfer) and basin-reservoir routines. The purpose of such a coupling was to achieve a trade-off between the expected performance of ecophysiological and hydrological models, which are often employed separately and at different spatial scales, either the plot or the basin. The calibration of the model to perform streamflow yielded a Nash-Sutcliffe (NS) coefficient equal to 0.89 for the year 2009, while the Correspondence to: F. Gómez-Delgado (fgomezd@ice.go.cr) validation of the water balance partitioning was consistent with the independent measurements of actual evapotranspiration (R2 = 0.79, energy balance closed independently), soil water content ( R2 = 0.35) and water table level ( R2 = 0.84). Eight months of data from 2010 were used to validate modelled streamflow, resulting in a NS = 0.75. An uncertainty analysis showed that the streamflow modelling was precise for nearly every time step, while a sensitivity analysis revealed which parameters mostly affected model precision, depending on the season. It was observed that 64% of the incident rainfall R flowed out of the basin as streamflow and 25% as evapotranspiration, while the remaining 11% is probably explained by deep percolation, measurement errors and/or inter-annual changes in soil and aquifer water stocks. The model indicated an interception loss equal to 4% of R, a surface runoff of 4% and an infiltration component of 92%. The modelled streamflow was constituted by 87% of baseflow originating from the aquifer, 7% of subsurface nonsaturated runoff and 6% of surface runoff. Given the low surface runoff observed under the current physical conditions (andisol) and management practices (no tillage, planted trees, bare soil kept by weeding), this agroforestry system on a volcanic soil demonstrated potential to provide valuable HES, such as a reduced superficial displacement-capacity for fertilizers, pesticides and sediments, as well as a streamflow regulation function provided by the highly efficient mechanisms Published by Copernicus Publications on behalf of the European Geosciences Union. 370 F. Ǵomez-Delgado et al.: Modelling hydrological behaviour of coffee agroforestry basin in Costa Rica of aquifer recharge and discharge. The proposed combination of experimentation and modelling across ecophysiological and hydrological approaches proved to be useful to account for the behaviour of a given basin, so that it can be applied to compare HES provision for different regions or management alternatives.


Introduction
The ability of ecosystems to infiltrate rainfall, sustain aquifers, and avoid erosion is a key determinant for the provision of hydrological environmental services (HES), especially in the humid tropics where surface fluxes can be very high (MEA, 2005). Woody plants and, in particular, agroforestry (AF) systems associating shade trees and perennial crops with deep root systems, are assumed to enhance these HES in comparison to traditional intensive cropping systems (Siles et al., 2010a;Ataroff and Monasterio, 1997;Vaast et al., 2005). Costa Rica is renowned as a promoter of HES by charging water users for the HES they receive from land owners (e.g. forest conservation), focusing on water quality (Pagiola, 2008). Hydropower producers, generating 78% of the total electricity consumption in Costa Rica during 2008 (ICE, 2009), are major HES payers. Coffee is one of the most traded agricultural commodities in the world employing 100 million people (Vega and Rosenquist, 2001). In Costa Rica, coffee accounted for 15% of the agricultural exports in 2008 and covered 2% of the territory (SEPSA, 2009). As coffee plantations are present in the main basins used for hydroelectric generation in Costa Rica, the possible trade-offs of the payment for HES from hydropower producers to coffee farmers become evident. Negotiation for these payments is facilitated between providers and purchasers when the service, and/or the impact of a given practice on the provision of the service, are clearly evaluated. However, links between land use, management practices (like tree planting) and hydrology in Costa Rica have not been thoroughly investigated by quantitative research (Anderson et al., 2006). There is a need of both, experimentation at the basin scale in order to evaluate the main hydrological processes, and of integrated modelling to understand the behaviour of all water compartments, including hidden ones (e.g. the aquifer).
The partitioning of the water balance (WB) is a prerequisite to evaluate HES such as infiltration, aquifer regulation capacity, erosion control and contaminants retention in coffee AF systems. Comprehensive WB studies at basin scale, including closure verification by independent methods, have been carried out in the developed world and for other land covers, like those reported by Roberts and Harding (1996), Dawes et al. (1997), Ceballos and Schnabel (1998), Wilson et al. (2001) and Maeda et al. (2006). Some experimental basins are located in the tropics, like those in Brazil, Costa Rica, Guadeloupe and Panamá (Charlier et al., 2008;Fujieda et al., 1997;Genereux et al., 2005;Kinner and Stallard, 2004), but no coffee AF basins have been equipped so far. Some reports are available for coffee AF systems but at the plot level and for some particular fluxes such as throughfall and stemflow (Siles et al., 2010b), tree and coffee transpiration (Dauzat et al., 2001;van Kanten and Vaast, 2006), surface runoff (Harmand et al., 2007), energy balance and latent heat flux (Gutiérrez et al., 1994). To our knowledge, there is no comprehensive study of the water balance partitioning of coffee AF systems at the basin level, including the behaviour of the aquifer.
Truly balanced combinations of hydrological and ecophysiological experiments and models remain scarce, although they intrinsically carry a more realistic and comprehensive representation of plant, soil and aquifer components at plot and basin scales. Most hydrological studies at basin scale use flumes for monitoring the streamflow and simply estimate evapotranspiration (ET), which prevents a true verification of the water balance closure.
As in the tropics we assumed that ET, including the reevaporation of intercepted water (R In ), is an important component of the water balance, even for precipitations around 3000 mm yr −1 , we decided to measure it directly by eddycovariance (Roupsard et al., 2006;Baldocchi and Meyers, 1998;Wilson et al., 2001), choosing a 0.9 km 2 micro-basin embedded in a very homogeneous coffee AF plantation. As an additional advantage, the eddy-covariance method can be validated itself by closing the energy balance (Falge et al., 2001).
Lumped, conceptual rainfall-streamflow models have been used in hydrology since the 1960s (e.g. Crawford and Linsley, 1966;Cormary and Guilbot, 1969;Duan et al., 1992;Bergström, 1995;Donigan et al., 1995;Havnø et al., 1995;Chahinian et al., 2005). These models consider the basin as an undivided entity, and use lumped values of input variables and parameters. For the most part (for a review, see Fleming, 1975;Singh, 1995), they have a conceptual structure based on the interaction between storage compartments, representing the different processes with mathematical functions to describe the fluxes between the compartments. Most hydrological models simplify the ET component based on reference ET routines (Allen et al., 1998) or using very empirical, non-validated models for actual ET. For instance, many models assume that actual evapotranspiration equals the reference ET, or a constant fraction of it, or a variable fraction that depends only on the soil water content in the root zone. However, improper parameterization of the crop coefficient may severely affect the parameterization of hydrological resistances and fluxes. In constrast, ecophysiological models may operate efficiently at plot level but miss the partitioning between lateral non-saturated (subsurface) runoff and vertical drainage, and the dynamics of water in aquifers and rivers. This is a major limitation for the assessment of HES, which is mainly desired at the basin scale.  In the present study we attempted to couple two lumped models into a new integrative one, chosen to be scalable and parsimonious: a basin reservoir model similar to the CREC model (Cormary and Guilbot, 1969) and employing the Diskin and Nazimov (1995) production function as proposed by Moussa et al. (2007a,b), and the SVAT model proposed by Granier et al. (1999). While the basin model was considered appropriate for its simplicity and capacity to support new routines, the SVAT model was chosen for its parsimony (three parameters in its basic formulation), its robustness (uses simple soil and stand data in order to produce model runs for many years, avoiding hydraulic parameters that are difficult to measure and scale up), its ability to quantify drought intensity and duration in forest stands, and for its successful past validation in various forest stands and climatic conditions, including tropical basins (Ruiz et al., 2010).
This paper aims to explain and model the hydrological behaviour of a coffee AF micro-basin in Costa Rica, assessing its infiltration capacity on andisols. The methodology consists of experimentation to assess the main water fluxes and modelling to reproduce the behaviour of the basin. First, we present the study site and the experimental design. Second, we develop a new lumped eco-hydrological model with balanced ecophysiological/hydrological modules (that we called Hydro-SVAT model). This model was tailored to the main hydrological processes that we recorded (streamflow, evapotranspiration, water content in the non-saturated zone and water table level) and that are described in the subsequent sections. Third, we propose a multi-variable calibration/validation strategy for the Hydro-SVAT model so we calibrate using the streamflow in 2009 and validate using the remaining three variables in 2009 and the streamflow in 2010. Fourth, we make an uncertainty analysis to produce a confidence interval around our modelled streamflow values, and a sensitivity analysis to assess from which parameters this uncertainty might come. Fifth, we attempt to simplify the initially proposed model based on the results of the sensitivity analysis. Finally, we discuss the main findings concerning the water balance in our experimental basin.

Location, soil and climate
The area of interest is located in Reventazón river basin, in the Central-Caribbean region of Costa Rica ( Fig. 1a and b). It lies on the slope of the Turrialba volcano (central volcanic mountain range of the country) and drains to the Caribbean Sea. The Aquiares coffee farm is one of the largest in Costa Rica (6.6 km 2 ), "Rainforest Alliance TM " certified, 15 km from CATIE (Centro Agronómico Tropical de Investigación y Enseñanza). Within the Aquiares farm, we selected the Mejías creek micro-basin (Fig. 1c) for the "Coffee-Flux" experiment. The basin is placed between the F. Gómez-Delgado et al.: Modelling hydrological behaviour of coffee agroforestry basin in Costa Rica coordinates −83 • 44 39 and −83 • 43 35 (West longitude), and between 9 • 56 8 and 9 • 56 35 (North latitude) and is homogeneously planted with coffee (Coffea arabica L., var Caturra) on bare soil, shaded by free-growing tall Erythrina poeppigiana trees. The initial planting density for coffee was 6300 plants ha −1 , with a current age >30 years, 20% canopy openness and 2.5 m canopy height. It is intensively managed and selectively pruned (20% per year, around March). Shade trees have a density of 12.8 trees ha −1 , with 12.3% canopy cover and 20 m canopy height. The experimental basin has an area of 0.9 km 2 , an elevation range from 1020 up to 1280 m a.s.l. and a mean slope of 20%. Permanent streams extend along 5.6 km, implying a drainage density of 6.2 km km −2 . The average slope of the main stream is 11%.
According to the classification by Mora-Chinchilla (2000), the experimental basin is located along a 1.3 km wide strip of volcanic avalanche deposits, characterized by chaotic deposits of blocks immersed in a matrix of medium-to-coarse sand, which is the product of the collapse of the southeastern slope of Turrialba volcano's ancient crater. The general classification given by the geological map of Costa Rica (MINAE-RECOPE, 1991) describes the general stratigraphy as shallow intrusive volcanic rocks, and the particular region as proximal facies of modern volcanic rocks (Quaternary), with presence of lava flows, agglomerates, lahars and ashes. Soils belong to the order of andisols according to the USDA soil taxonomy (USDA, 1999), which are soils developing from volcanic ejecta, under weathering and mineral transformation processes, very stable, with high organic matter content and biological activity and very large infiltration capacities.
According to Köppen-Geiger classification (Peel et al., 2007), the climate is tropical humid with no dry season and strongly influenced by the climatic conditions in the Caribbean hillside. The mean annual rainfall in the study region for the period 1973-2009 was estimated as 3014 mm at the Aquiares farm station (Fig. 2). At the experimental basin the rainfall in 2009 (3208 mm) was close to the annual mean, but showed a monthly mean deviation of ±100 mm around the historical regime. Mean monthly net radiation ranged in 2009 from 5.7 to 13.0 MJ m −2 d −1 , air temperature from 17.0 to 20.8 • C, relative humidity from 83 to 91%, windspeed at 2 m high from 0.4 to 1.6 m s −1 and Penman-Monteith reference evapotranspiration (Allen et al., 1998) from 1.7 to 3.8 mm d −1 .

Experimental setup
The "Coffee-Flux" experimental basin and instrument layout was designed to trace the main water balance components employing spatially representative methods (Fig. 1c). It is part of the FLUXNET network for the monitoring of greenhouse gases of terrestrial ecosystems. The hydrological measurements were recorded from December 2008 up to February 2010.  Rainfall and climate: rainfall was monitored at 3 m above ground in the middle of 3 transects of the basin, using three lab-intercalibrated ARG100 tipping-bucket (R. M. Young, MI, USA) connected to CR800 dataloggers (Campbell Scientific, Shepshed, UK), and integrated every 10 min. Other climate variables were logged on top of the eddy-flux tower with a CR1000, every 30 s, integrated half-hourly and using: Net radiation: NR-Lite (Kipp and Zonen, Delft, The Netherlands); PPFD: Sunshine sensor BF3 (Delta-T devices Ltd, UK); temperature and humidity: HMP45C in URS1 shelter (Campbell Scientific); wind-speed and direction: 03001 Wind Sentry (R. M. Young, MI, USA). The theoretical evapotranspiration from a wet grass placed under local climate conditions, ET 0 , was computed in accordance with Allen et al. (1998).
Streamflow: a long-throated steel flume (length: 3.9 m; width: 2.8 m; height: 1.2 m) was home-built to measure the streamflow at the outlet of the experimental basin, to record up to 3 m 3 s −1 , the maximum estimated discharge for the study period from an intensity-duration-frequency analysis. The flume was equipped with a PDCR-1830 pressure transducer (Campbell Scientific) to record water head at gauge point (30 s, 10 min integration), while the rating curve was calculated considering the geometric and hydraulic properties of the flume using Winflume software (Wahl et al., 2000). A validation of the rating curve was made successfully using the salt dilution method as well as a pygmy current meter.
Soil water content: a frequency-domain-reflectometry portable probe (FDR Diviner2000, Sentek Pty Ltd) was used to survey 20 access tubes distributed in the three study transects to provide the mean volumetric soil water in the basin. The sensor measures at 10 cm intervals, reaching a total depth of 1.6 m. A measurement campaign through the 20 sites was carried out every week. The sensors were calibrated by digging sampling pits in the vicinity of six test tubes, to obtain the actual volumetric soil water content from gravimetric content and dry bulk density.
Evapotranspiration: the actual evapotranspiration from the soil, coffee plants and shade trees was measured at reference height (26 m) on the eddy-covariance tower, similarly to Roupsard et al. (2006). . Raw data were collected and preprocessed by "Tourbillon" software (INRA-EPHYSE, Bordeaux, France) for a time-integration period of 300 s, then post-processed using EdiRe software (University of Edinburgh, UK) into half-hourly values and quality checked. A validation was made by direct comparison of the measured net radiation R n with the sum of sensible heat flux (H ) and latent heat flux (λE): at daily time step, this yielded H + λE = 0.92 R n (R 2 = 0.93) which was considered sufficiently accurate to assume that advection effects on λE could be neglected here. Due to lighting and sensor breakdown, 45 days of data were lost between July and August 2009. To gap-fill the missing period we used the Penman-Monteith model, whose canopy conductance was adjusted using measured values. Leaf Area Index (LAI): the coffee light transmittance was measured monthly in diffuse light conditions, for five rings at different zenital angles (LAI2000, Li-COR Corvallis, USA), along three 50 m-long transects through the flux tower plot, similarly to Roupsard et al. (2008). Effective coffee LAI, obtained from this light transmittance, was converted into actual LAI according to Nilson (1971), using a ratio of effective to actual LAI that was estimated from a dedicated calibration. The actual coffee LAI was measured directly on a small plot by counting total leaf number of 25 coffee plants, measuring leaf length and width every 20 leaves and using empirical relationships between leaf length and width and leaf area (LI-3100C, Li-COR) (R 2 > 0.95). On the same small plot, the effective LAI was measured with LAI2000. The ratio of effective to actual LAI was then calculated on this small plot (1.75) and was considered to be constant with time and space in the micro-basin, allowing the estimation of the actual LAI on the three LAI2000 transects. The LAI for shade trees was estimated using their crown cover projection (on average 12.3% over the whole basin) observed on a very high resolution panchromatic satellite image (World-View image, February 2008, 0.5 m resolution). As we did not have measurements of LAI for shade trees, we considered this LAI in the order of magnitude of coffee LAI on a crown-projected basis, and therefore we multiplied the actual coffee LAI measured on transects by 1.123, to estimate the ecosystem LAI (tree and coffee). In order to monitor the time-course of ecosystem LAI at the basin scale, we combined these ground measurements with time series of remotely-sensed images. We used time series of Normalized Difference Vegetation Index (NDVI) from Moderate Resolution Imaging Spectroradiometer (MODIS) data products MOD13Q1 and MYD13Q1 (16-Day composite data, 250 m resolution). NDVI is known to be correlated with the green LAI if it is low, for most ecosystems (Rouse et al., 1974). Twenty-three MODIS pixels covering the experimental basin were selected, and their NDVI time series were downloaded. We filtered the raw NDVI time series according to quality criterion given in the MODIS products, and we adjusted a smooth spline function on it as in Marsden et al. (2010). Then, a linear regression between the smoothed NDVI of the pixel including the flux tower and ground values of actual ecosystem LAI was calibrated (R 2 = 0.69). This regression was used on other pixels of the basin, and averaged to have the annual time-course of actual LAI.
Water table level: four piezometric wells measuring up to 4 m depth were built in the three main transects of study. They were equipped with pressure transducers (Mini-Divers, Schlumberger Water Services) that measure and record the water table level every 30 min.
Period of measurement, data gaps and gap-filling: the recording information is given in Table 1 for the five hydrologic variables. The frequency of measurement varies, but is finally calculated at the 30 min time step (except for soil water content that is a non-continuous measurement). When gaps are present in the measurements, a gap filling method was applied.

Hydro-SVAT lumped model
We designed a lumped, five-reservoir-layer eco-hydrological model in order to predict the water balance (WB) partitioning (stocks and fluxes) at the scale of the entire basin. It is based on the water balance models developed by Moussa et al. (2007a,b) and Granier et al. (1999), and built to reproduce the main hydrological processes measured at the experimental basin, which will be presented in Sect. 4. The model of Moussa et al. (2007a,b) works at the basin scale and simulates the ecosystem evapotranspiration rather roughly, while the one of Granier et al. (1999) works at the plot scale and totally ignores the lateral water fluxes through the soil and the role of the basin aquifers. The main novelties of the Hydro-SVAT model with respect to the model structure of Moussa et al. (2007a,b) are the inclusion of a land cover reservoir to separate the intercepted rainfall from the combined throughfall/stemflow component, and the partition of non-saturated soil into two reservoirs, one with and one without roots of plants and trees. The first of these innovations intends to take into account the non-negligible interception loss in coffee AF systems, as reported by Jiménez (1986), Harmand et al. (2007) and Siles (2007). The second new addition to the model is to better represent the water dynamics in the nonsaturated soil, given that only its upper layer will lose humidity by root extraction. The water balance model of Granier et al. (1999) is incorporated in this superficial reservoir but in a simplified form, so that both, the root distribution and the soil porosity, are homogeneous through the vertical, nonsaturated profile. Hence, the water content in this reservoir is the state variable linking our two parent models.
The modelling hypotheses governing the model architecture were: (a) the interception loss component is not negligible in the WB and is a function of rainfall intensity, (b) infiltration is a function of the soil water content in the nonsaturated reservoirs, (c) evapotranspiration is a significant component in the WB and is best described using a SVAT model that couples evapotranspiration to root water extraction from the soil, (d) the aquifer has a higher discharge rate above a threshold level, (e) any lateral inputs to the basins system are considered negligible in comparison to rainfall inputs, and (f) there is a net water outflow from the system as deep percolation.
The model was implemented using Matlab ® V. R2007a (The MathWorks Inc., USA).

Model structure
The model structure is presented in Fig. 3. The next three sections will describe the model structure according to its three major routines and five layers. The first layer is called "land cover reservoir" and separates the total rainfall into an intercepted loss and a joint throughfall/stemflow component. The second layer or "surface reservoir" regulates the surface runoff. The infiltration process from the second layer is controlled by the joint water content at the third and fourth layers, called "non-saturated root reservoir" and "nonsaturated non-root reservoir", respectively. The evapotranspiration flux is calculated at the "non-saturated root reservoir", while both non-saturated layers control the drainage, the percolation and the non-saturated runoff processes. The fifth and last layer is the "aquifer reservoir", which determines the baseflow and the deep percolation. Finally, we will explain the sum of the total runoff and baseflow components and the routing procedure to generate the modelled streamflow. Let A(t), B(t), C(t), D(t) and E(t) [L] be the water levels at time t in the five reservoirs A, B, C, D and E, respectively (or land cover reservoir, surface reservoir, non saturated root reservoir, non-saturated no-root reservoir and aquifer reservoir). Let A X , B X , C X , D X and E X [L] be the water levels corresponding to the maximum holding capacities for the five reservoirs.

Infiltration and actual evapotranspiration (a) Infiltration
The infiltration process i [LT −1 ] occurs from the second layer (surface reservoir) to the third one (non-saturated root reservoir), and eventually to the fourth one (non-saturated non-root reservoir) when i fills the third one. The infiltration capacity f i (t) [LT −1 ] is a state variable that depends on the water level in the non-saturated root reservoir, given by C(t). In addition to C X (maximum) we define C F as the water level in the root reservoir at field capacity. Then, f i (t) is calculated as (see Fig. 4a): is the infiltration rate at field capacity. The infiltration i both modifies and depends on B (t), which is the water availability in the second reservoir before i is extracted, according to: and and The infiltration module calculates the infiltration i as output variable, using the state variables B(t), C(t) andf i (t). Four parameters (C X , C F , f c and α) are demanded.

(b) Evapotranspiration
The evapotranspiration component ET [LT −1 ] acts directly on the third layer (non-saturated root reservoir) and is the sum of E u [LT −1 ] the understory and soil evaporation, and of T [LT −1 ] the transpirational water uptake by roots.
According to Shuttleworth and Wallace (1985), the fraction of total evapotranspiration originating from the plants is close to 100% of the total evapotranspiration of the ecosystem when LAI > 3 an when the soil is not saturated at its surface, which was always the case in our study. We thus assumed for simplicity that E u , the evaporation from the soil, was nil.
Transpiration T is obtained by solving T from the slightly modified ratio: r = T ET −1 0 [dimensionless] proposed by Granier et al. (1999). We substituted the original Penman potential evapotranspiration PET in that ratio by the Penman-Monteith reference evapotranspiration ET 0 [LT −1 ] (Allen et al., 1998). While ET 0 was calculated at each time step t, we estimated r as a function of the relative extractable water REW(t) [dimensionless], a state variable given by Granier et al. (1999) as: The REW(t) is linked to the soil water content according to: with θ(t): volumetric soil water content [L 3 L −3 ] at time t, θ r : residual soil water content [L 3 L −3 ] and θ f : soil water content at field capacity [L 3 L −3 ]. The parameter REW c [dimensionless] is the critical REW(t) below which the transpiration of the system begins to decrease. Figure 4b shows an example of some r curves as a function of REW(t). Each curve can be defined only by REW c and the rm LAI , a maximum value for the ratio r that depends on the LAI of the system as: where LAI X is the maximum measured LAI during the modelling period and r m is a parameter indicating the maximum ratio T ET −1 0 that can be found in this system. Then: Finally, we find the transpiration as T = rET 0 . The total modelled evapotranspiration including the interception loss, can be calculated as: ETR m = E u + T + R In , with R In [LT −1 ] being the intercepted/evaporated rainfall loss that will be explained in the next section. Hence, ETR m can be directly compared to the evapotranspiration that we measured at the flux tower. This module provides the evapotranspiration ET as a function of the state variable C(t), two input variables (LAI and ET 0 ) and three parameters (C X , REW c and r m ).

Water balance in the model reservoirs (a) Land cover reservoir
The first layer of the model, denoted "land cover reservoir", represents the soil cover in the basin and controls the partition of the total incident rainfall R [LT −1 ] into intercepted (then evaporated) rainfall loss R In [LT −1 ] and the combined troughfall/stemflow R TS [LT −1 ]. A simple water balance of this reservoir is established to calculate a proxy A (t) [L] of the final water level A(t) for each time step t, by adding the incident rainfall R and subtracting the Penman potential evapotranspiration PET [LT −1 ] from the existing land cover humidity level A(t − 1): We calculated the water level A(t) in this reservoir as well as R TS and R In by differentiating three cases: The land cover module calculates at each time t the water level A(t) as a state variable, demanding two input variables (R and PET) and one parameter (A X ). It yields the partition of R into R In and R TS .

(b) Surface reservoir
The second layer is called "surface reservoir" and acts as a sheet top soil with a given roughness and surface runoff delaying properties. The water balance in this surface reservoir for a given interval t is: where R TS [LT −1 ] is the combined throughfall/stemflow component from the previous layer and Q B1 and Q B2 [LT −1 ] are the non-immediate and immediate surface runoffs calculated as: where k B [T −1 ] is a discharge parameter, and: is a function of the water content in the third layer and it is the last component to be evaluated in the surface reservoir. This surface reservoir module calculates the water level B(t) as a state variable and demands one input variable (R TS ), and two parameters for the reservoir (B X and k B ). It produces three output variables: i, Q B1 and Q B2 . The two latter variables constitute the surface runoff in the basin (Fig. 4c).

(c) Non-saturated root reservoir
The "non-saturated root reservoir" is the third layer of the model and it represents a soil layer with presence of root systems from trees and plants. The water balance here is: is the state variable of the water level at a given time t, i is the infiltration from the second layer, ET [LT −1 ] is the evapotranspiration, d 1 and d 2 [LT −1 ] are the non-immediate and immediate drainages to the fourth layer, respectively and Q C [LT −1 ] is the non-saturated runoff from the root reservoir. There will be immediate drainage d 2 if at anytime the R TS component fills the reservoir above C X . Then d 2 = [C (t) − C X ] t −1 goes to the fourth layer and C(t) is reset to C X . Both non-immediate drainage d 1 and non-saturated runoff Q C occur whenever C(t) is higher than the field capacity threshold C F [L] according to: where ρ [LT −1 ] is the total outflow capacity in this reservoir and k C [T −1 ] a discharge parameter. The partition of ρ in d 1 and Q C depends on a parameter β [dimensionless], with 0 < β < 1. Then: The root soil module calculates the water level C(t) as state variable using two input variables (i and ET) and four parameters (C X , C F , k C and β). It provides three outputs (d 1 , d 2 and Q C ).

(d) Non-saturated non-root reservoir
The fourth layer of the model is denoted "non-saturated nonroot reservoir" and represents a soil layer with total absence of root systems and hence, of root water extraction. The water balance here is given by: where D(t) [L] is the state variable of the water level at a given time t, d 1 and d 2 [LT −1 ] are the non-immediate and immediate drainages from the third layer, respectively; Q D [LT −1 ] is the non-saturated runoff from the non-root reservoir and g 2 and g 1 [LT −1 ] are the immediate and nonimmediate percolation to the fifth model layer, respectively. Immediate percolation g 2 will be produced if at anytime the drainage (d 2 and/or d 1 ) fills the reservoir above D X . Then g 2 = [D(t) − D X ] t −1 moves to the aquifer reservoir and D(t) is reset to D X . Both non-immediate percolation g 1 and non-saturated runoff Q D occur whenever D(t) is higher than the field capacity threshold D F [L]: where η [LT −1 ] is the total outflow capacity of this reservoir and k D [T −1 ] a discharge parameter. The partition of η in g 1 and Q D depends on the parameter β [dimensionless]. Then: The non-root soil module calculates the water level D(t) as state variable using two input variables (d 1 and d 2 ) and four parameters (D X , D F , k D and β). It provides three outputs (g 1 , g 2 and Q D ).

(e) Aquifer reservoir
A fifth layer called "aquifer reservoir" represents the groundwater system and controls baseflow and deep percolation. The reservoir is composed by a shallow aquifer that acts whenever the water level in the reservoir is higher than E X , and by a deep aquifer with a permanent contribution. The water balance here is: (25) where E(t) [L] is the state variable of the water level at a given time t, g 1 and g 2 [LT −1 ] are respectively the nonimmediate and immediate percolation from the fourth layer, Q E1 and Q E2 [LT −1 ] are the baseflow from deep and shallow aquifers respectively (Fig. 4d), and DP [LT −1 ] is the deep percolation. If where k E1 , k E2 , and k E3 are discharge parameters controlling deep/shallow aquifers and deep percolation, respectively. This module calculates the water level E(t) as state variable using two input variables (g 2 and g 1 ) and four parameters (E X , k E1 , k E2 and k E3 ), to provide three outputs (Q E1 , Q E2 and DP).

Total runoff, baseflow and streamflow
The components of surface runoff, non-saturated runoff and baseflow are added to obtain the total runoff Q T [LT −1 ]: As explained in Moussa and Chahinian (2009) the streamflow Q [LT −1 ] at the outlet of the basin is obtained by the routing of Q T using a transfer function (to take into account the water travel time). The Hayami (1951) kernel function (an approximation of the diffusive wave equation) is developed to obtain a unit hydrograph linear model for this purpose. That is:

Model parameterization, calibration and validation
Summarizing, this Hydro-SVAT model uses four input variables: rainfall R, Penman-Monteith ET 0 , Penman PET and leaf area index LAI to generate five main output variables: interception R In , infiltration i, evapotranspiration ET, discharge components Q = Q B +Q C +Q D +Q E (from surface, non-saturated and aquifer reservoirs, respectively) and deep percolation DP. Five state variables are calculated for every time step, the water levels in the five reservoirs: A(t), B(t), C(t), D(t) and E(t). A coupled discharge Q CD for the two non-saturated reservoirs is obtained by adding Q C and Q D .
In our experimental basin we applied the model for a oneyear period (2009), a time step t = 30 min (1800 s) and a basin area equal to 0.886 km 2 .
This model contains 20 parameters that are used to calculate infiltration (C X , D X , C F , D F , f c and α), evapotranspiration (REW c and r m ), the exchange between reservoirs (A X , B X , k B , k C , k D , β, E X , k E1 , k E2 and k E3 ) and the basin transfer function (w and z F ).
Four out of these twenty parameters (C X , D X , C F , D F ) were estimated using field data. For instance, two excavation experiments down to 3.5 m showed that very few roots were present below 1.5 m, where the andisol layer turns into a more clayey, compact and stony deposit. Then, the depth of the non-saturated root soil layer was fixed at C H = 1.6 m (for simplicity, equal to the length of our FDR probe tubes). The depth of the non-root layer was estimated in D H = 1.0 m. Following the relationships C X = (θ s − θ r ) C H and D X = (θ s − θ r ) D H , the levels for maximum water holding capacities in the non-saturated reservoirs can be calculated, as well as the levels for field capacities, using the equations C F = (θ f − θ r ) C H and D F = (θ f − θ r ) D H . The volumetric soil water contents θ were estimated as θ r = 0.37: the residual water content equal to the minimum θ observed in the basin during the study period; θ s = 0.63: the θ at saturation for a typical andisol, according to Hodnett and Tomasella (2002); and θ f = 0.43: the θ at field capacity equal to the average van Genuchten value for a matric potential = −10 kPa (Hodnett and Tomasella, 2002).
Three parameters were taken from literature reviews and expert criteria (A X , REW c and r m ). We set the surface reservoir maximum storage capacity A X for our coffee AF system equal to 4 × 10 −4 m using data from Siles et al. .(2010b). The two parameters of the evapotranspiration routine were taken as: REW c = 0.4 from Granier et al. (1999) and r m = 0.8 from field measurements of T ET −1 0 (data not shown). One parameter (B X ) was obtained at the end of the optimization process in order to fit the maximum observed streamflow peak (this parameter is very sensitive as it acts directly on the highest peaks). Two other parameters (w and z F ) were separately estimated by a trial and error procedure, given the low sensitivity of the model to their variation.
The remaining 10 empirical parameters (f c , α, k B , k C , k D , β, E X , k E1 , k E2 and k E3 ) were simultaneously optimized. For this purpose we used the Nelder-Mead (Nelder and Mead, 1965) simplex algorithm included in Matlab (following Lagarias et al., 1998) on 17 520 semi-hourly time steps (one year). Convergence was reached within 1000 runs and the stabilization of all parameter values by the end of the iteration process was checked. A two-step calibration procedure was applied: (a) selection of an initial value for each parameter, falling within the respective range (fourth column, Table 2), and (b) simultaneous estimation of parameter values that maximize an objective function (sixth column of Table 2, identified as M 1 ), in this case the Nash and Sutcliffe (1970) efficiency coefficient. We calibrated the model using one year of streamflow Q, and then validated it using the remaining three measured variables: evapotranspiration ETR, water content in the nonsaturated zone θ and water table level z. For ETR we grouped the measured and modelled values (given in the same units) at the daily time scale, which is the original time scale in Granier's model, and then we excluded the gap-filled values. To calculate the modelled θ from the water level in the root reservoir C(t) we used the relation θ = C (t) C −1 X (θ s − θ r ) + θ r , while the observed values were obtained as the average of the 20 FDR point-measurements throughout the basin. To validate z we proposed an effective porosity of the aquifer n A = 0.39, to be able to directly link piezometric measurements z with the modelled water level in our aquifer reservoir E(t), according to z = E(t)/n A .
We also carried out an independent validation of Q by running the calibrated Hydro-SVAT model over 8 months of data from 2010 and verifying the model performance.
In addition, we performed an analysis of model residuals: zero expectancy, normality, homoscedasticity and standardized residuals.

Uncertainty and sensitivity analysis
In order to investigate the uncertainty in model predictions we performed a Monte-Carlo approach on a restricted subset of parameter combinations, as suggested by Helton (1999). For each of the 10 parameters a range was created with  White et al., 2000;Ines and Droogers, 2002;Lenhart et al., 2002;Zaehle et al., 2005;Droogers et al., 2008), following the "one-at-a-time" method described by Hamby (1994Hamby ( , 1995 and Frey and Patil (2002). Then, we assumed a uniform distribution for all the parameters, and used the Latin Hypercube function of Simlab 2.2 (http://simlab.jrc.ec.europa.eu) in order to produce a sample from a joint probability distribution, of size equal to ten times the number of parameters as recommended by Sim-Lab developers, i.e., 100 parameter combinations. These combinations were introduced into the calibrated Hydro-SVAT model and we retrieved 100 streamflow output series over the 17 520 semi-hourly time steps. Then we generated 17 520 empirical confidence intervals (95% and 99%) from the respective frequency distributions over the 100 parameter combinations.
A sensitivity analysis (SA) was carried out to determine which of the input variables contributed significantly to this modelling uncertainty. Two separate assessments were produced. First, a summary (through the time) index of model performance was studied: the Nash-Sutcliffe coefficient (NS), being evaluated by four sensitivity indexes: Pearson, Spearman, standardized regression and standardized rank regression coefficients (Hamby, 1994(Hamby, , 1995. Then, a second approach of sensitivity analysis was tested to try to follow the behaviour of the Spearman coefficient for each of the 10 parameters included in the SA through all time steps, as in Helton (1999).
The sensitivity indexes are calculated departing from the joint probability sample matrix x ij of size m × n, where m is the sample size and n the number of independent variables (here our 10 parameters) to study. The Monte Carlo evaluation of x ij in the model produces the result vector y i , configuring the matrix system [y i :x ij ]. Then, the Pearson product moment correlation (PEAR) for a given parameter j is the linear correlation coefficient between the variables x ij and y i over the m samples. To account for non-linear relationships that can be hidden by indicators like PEAR, a simple rank transformation is applied, replacing the original x − y series with their corresponding ranks R(x) and R(y). Then, the Spearman coefficient (SPEA) is obtained by calculating the correlation on the transformed data, as SPEA(x, y) = PEAR[R(x), R(y)]. The standardized regression coefficients (SRC) are the result of a linear regression analysis performed on [y i :x ij ] but previously standardizing all the variables. This is useful to evaluate the effect of the independent variables on the dependent one, without regarding their units of measurement, and can be computed by standard statistical methods. Finally, the standardized rank regression coefficients are obtained as SRRC(x, y) = SRC[R(x), R(y)].
In the 0.9 km 2 micro-basin of Mejías creek, within the Cafetalera Aquiares AF coffee farm, we obtained one full year (2009) of comprehensive experimental results of streamflow, evapotranspiration, soil water content and piezometry that we used to calibrate and validate the Hydro-SVAT model. First we present the hydrological behaviour of the basin, then the ecophysiological behaviour, and finally the water balance.

Hydrological behaviour of the basin
The time series of streamflow Q and rainfall R are given at a semi-hourly time-step in Fig. 5a for 2009. Rainfall (Fig. 2) was quite evenly distributed (no marked dry spell), although the period from January to June clearly received less rain (later named the "drier season", in opposition to the "wetter season"). From the total R (3208 mm), the measured Q at the outlet was 2048 mm, yielding an annual streamflow coefficient of 0.64. The Q hydrograph (Fig. 5a) displays a continuous baseflow with episodes of groundwater recharge after rainfall events, followed by marked recessions controlled by the baseflow. Q peaks reached an annual maximum of 0.84 m 3 s −1 , i.e. 28% of the nominal capacity of the flume, indicating that the size chosen for the flume was adequate. It can be observed that similar rainfall events resulted in higher Q peaks during the wetter season. The lower Q peaks in response to rainfall events during the drier season could thus be interpreted as the consequence of higher infiltration rates when belowground was less saturated. Soil water content in the 0-1.6 m layer (Fig. 5d) remained above 37% all-year round and rose up to a maximum of 47% during the transition between the drier and the wetter season. In order to understand the large baseflow shaping the hydrograph of Q, four piezometers were installed in early June 2009 (Fig. 5e) and showed the existence of a permanent aquifer at levels varying between 0.7 and 3.2 m deep, according to their respective distance to water channels. Their behaviour was extremely variable, as expected, from responsive (piezos #1 and #3) to conservative behaviours (piezos #2 and #4). The continuous baseflow observed by the flume originated mainly from an important aquifer, covering a large (although undefined) area within the basin.

Ecophysiological behaviour of the basin
The ecosystem LAI (Fig. 5b) changed seasonally quite severely from 2.8 (in March), as a result of coffee pruning during the drier season and leaf shedding by E. poeppigiana, coffee flowering and new leafing just after the beginning of the wetter season, to a maximum of 4.8 (in September), then coffee leaf shedding during the main coffee-berry harvest (in October).
The actual evapotranspiration (ETR) obtained by eddycovariance (Fig. 5c) accounted for the sum of coffee and shade-tree transpiration, understory evaporation (mainly bare soil), and rainfall interception loss. The total ETR in 2009 amounted to 818 mm (25% of R), fluctuating daily according to atmospheric demand and seasonally according to LAI and canopy conductance. It always stayed below 4.5 mm d −1 and, in average, around 60% of reference ET 0 , clearly invalidating the use of ET 0 as a reliable indicator of ETR in coffee system eco-hydrological models. The crop coefficient (the ETR ET −1 0 ratio) was clearly lower during the drier season, as a consequence mainly of a lower LAI. We did not observe a period during which the relative extractable water (REW) of the soil drop below the critical value (REW c ) of 0.4, confirming that the coffee plants probably encountered no seasonal water stress. Figure 6 shows the water balance partition and closure for 2009, as obtained by the water flows measured by independent experimental methods, rainfall (R), streamflow (Q), and evapotranspiration (ETR) (including the rainfall interception loss). It was observed on cumulative values that the sum of Q + ETR was 11% lower than annual R. However, Q + ETR matched or exceeded R once at the end of April, just after three episodes of lower rainfall, which confirmed the occurrence of an important storage in the basin, namely the soil and aquifer reservoirs, creating a seasonal hysteresis between R and Q. On an annual basis, measured Q represented 64% of R and measured ETR amounted 25%. The remaining 11% was attributed to deep percolation, measurement errors and/or inter-annual changes in soil and aquifer water stocks. Figure 7a shows the result of the streamflow modelling for the year 2009, after calibrating the 10 parameters from Table 2. The model yielded a NS coefficient of 0.89 and R 2 = 0.88 for N = 17 520 semi-hourly time-steps. Baseflow and peakflow events appeared to be satisfactorily represented for the whole time series. The modelled partitioning of streamflow into surface runoff, non-saturated runoff and baseflow is presented in Fig. 7b. It indicated a prominent contribution of baseflow, as already inferred from visual inspection of the Q time series. Hillslope surface runoff and non saturated runoff were a minor part of Q. The water stored in the 4 lowest reservoirs of the model is shown in Fig. 7c to f. The water level at the surface reservoir was rarely greater than zero, displaying some few events of surface storage beyond the modelling time step (30 min). The non-saturated reservoirs were replenished during storm events and depleted by non-saturated runoff and evapotranspitation (this latter acting only in the root reservoir). Finally, the aquifer reservoir displayed the largest magnitude of variation, fluctuating seasonally by a factor of almost 4 (this factor is reduced to 1.3 when the aquifer effective porosity is taken into account).

Discussion
In the next sections we first discuss the validation, uncertainty and sensitivity analysis of the model, second some model simplification attempts, third the main hydrological processes that we observed, and finally we examine the hydrological services in our coffee AF basin.

Validation, uncertainty and sensitivity analysis for
the Hydro-SVAT model

Model validation
Model validation was done by direct comparison of model output variables with three field measurements: actual F. Gómez-Delgado et al.: Modelling hydrological behaviour of coffee agroforestry basin in Costa Rica  evapotranspiration ETR, soil water content θ and aquifer water level z. For ETR, the determination coefficient between the daily sum of observed and modelled ETR was R 2 = 0.79 (Fig. 8a). Concerning θ, the observed and modelled time series appeared to be consistent, but divergent during the driest season (Fig. 8b), reaching a rather low R 2 = 0.35. However, considering that this global θ was obtained as the arithmetic average of only 20 observations for the whole basin, the uncertainty of these measured values is high. The large amount of rocks hindering the tubes might also affect the FDR readings. Finally, we obtained a good approximation for the behaviour of the aquifer (Fig. 8c), with R 2 = 0.84. The use of the two piezometers that displayed the highest stability (probably representing the larger and more relevant aquifer systems) out of the four piezometers installed was crucial at this step. We considered to be recording two different processes, the typical gradual aquifer response (piezos #2 and #4 in Fig. 1c) and the local quick-varying shallow-water accumulations (piezos #1 and #3 in Fig. 1c), which might be associated with rapid changes in soil water contents at smaller and less representative aquifer units. Though the representativeness of these few piezometers of the behaviour of the main basin-aquifer could be questioned, the very high similarity that we found between the modelled values and the average measurements from piezometers #2 and #4 (located in two opposite ends of the basin), supported the idea of a correct performance of both, the field method for aquifer monitoring and the corresponding model routine that is based on a linear reservoir.
Concerning the precision of our model, the NS = 0.89 on streamflow seems to be good enough for a semi-hourly time step, considering the detail with which the hydrological processes need to be described. Some studies that have evaluated hydrological models using this coefficient for different time scales have shown the decline in NS values as the modelling time step is shortened. For instance, at calibration stages, Notter et al. (2007) achieved maximum NS values from 0.8 to 0.69 for decadal to daily time steps (basin area = 87 km 2 ), while Bormann (2006) obtained 0.8, 0.9, 0.85 and 0.73 for annual, monthly, weekly and daily time steps, respectively (basin area = 63 km 2 ). García et al. (2008) reached NS coefficients of 0.93, 0.91 and 0.61 for quarterly, monthly and daily modelling in a 162 km 2 basin.
In addition, the streamflow validation over an independent eight-month test period (January-August 2010) produced a NS = 0.75, which is considered satisfactory (Fig. 8d). During this first eight months of 2010 the total rainfall amounted to 1978 mm (compared to 2222 mm in the first eight months of 2009), the total streamflow summed 1057 mm (against 1445 mm in 2009) and 11 storm events exceeded the threshold of 0.3 m 3 s −1 (compared to 13 events in 2009).
To study possible systematic errors in the Hydro-SVAT model, we examined the distribution of residuals between measured and modelled streamflows (since this was the optimized variable). A t-test with a confidence level of 95% indicated that we cannot conclude that the mean of our residuals (equal to −1 × 10 −4 m 3 s −1 ) is significantly different from zero (p = 0.36). At our short time step it is common to find highly autocorrelated residuals, so we found significant partial autocorrelations up to the seventeen time lag (8.5 h). The Kolmogorov-Smirnov test revealed that the distribution of residuals is not normal, though it is highly symmetrical around zero. Studying residuals as functions of time, rainfall R, streamflow Q, soil water content θ and water table level z (data not shown), we did not find any trends, but noticeable changes in their variability make clear that the homoscedasticity condition was not properly fulfilled. However, it is accepted that these ordinary least squares assumptions are often not satisfied in streamflow modelling (Xu and Singh, 1998).

Uncertainty analysis
If we accept that the model reproduces efficiently the actual streamflow at the outlet of the experimental basin, the uncertainty in the modelled values needs to be known. A Monte-Carlo (MC) uncertainty analysis was produced as detailed in Sect. 3.3 to yield the empirical 95% and 99% confidence limits (CL) (the 95% CL are presented in Fig. 9, along with measured streamflow). 82% of the measured values fell within the 95% CL produced by our model, while 87% of them fell within the 99% CL. The ranges of the 95% confidence intervals (CI) along each of the time steps varied from 20% to 131% of the MC mean value, while for the 99% CI the ranges were from 24% up to 157% of the MC mean. From these analyses it is possible to state that the model is efficient (high NS coefficient), precise (relatively small confidence intervals, Fig. 9) and accurate (given the high percentages of measured Q values falling within the 95% and 99% confidence intervals).  (c) surface reservoir, (d) non-saturated root reservoir (offset by θ r C H to show the absolute water content in the soil layer), (e) non-saturated non-root reservoir (offset by θ r D H ) and (f) aquifer reservoir. Graphs on the right side correspond to a two-month period (June-July 2009) for better illustration.

Sensitivity analysis
The sensitivity analysis was carried out to determine the responsiveness of the model predictions to variations in our main parameters. The first assessment consisted in testing the statistical significance of the relationships between each of our 10 calibration parameters and the NS coefficient, using four dimensionless sensitivity indexes: Pearson (PEAR), Spearman (SPEA), standardized regression and standardized rank regression coefficients (SRC and SRRC, respectively). Table 3 presents the results of these tests revealing that, with the exception of the discharge coefficients from the root nonsaturated and shallow-aquifer reservoirs (k C and k E2 ), and of the shallow-aquifer threshold (E X ), none of the parameters seemed to significantly influence the global model efficiency. However, these results may be misleading given the enormously variable conditions under which the model works through time. Therefore, in a second approach we  selected the Spearman test (a simple rank transformation index that can identify non linear relationships) to assess the influence of each of these 10 parameters on streamflow, at each time step. Figure 10 presents the results, displaying the dimensionless Spearman index in the vertical axis and the time in the horizontal axis. A positive Spearman index reveals a proportional influence of the parameter on the model result, while negative values indicate inverse proportionality. Figure 10a shows the time series for all the parameters, revealing an alternation in their influence on streamflow over time and model state. The two black horizontal lines represent the 95% confidence limits above (or below) which the correlation is significantly different from zero. In Fig. 10b to k the significant values are plotted as black points in contrast to non-significant in grey, for each parameter separately. Figure 10h suggests that the parameter E X is again one of the most influential, because it is permanently displaying high positive or negative correlations, together with the discharge coefficient (DC) for the deep-aquifer k E1 , which most of the time has a proportional influence on model outputs (Fig. 10j). These two parameters reach the highest Spearman indexes during sustained periods (not only during storm events) and are clearly relevant during long streamflow recessions, when the modelled water table level z is around or below E X . The DC for the surface reservoir k B (Fig. 10c) is noticeably influencing model outputs during each individual storm event, while the vertical/lateral split coefficient β (Fig. 10b) is relevant only during the main ones. The DC for the deep percolation k E3 (Fig. 10k) is significant over recession periods, while DC for root reservoir k C is significant when its water content exceeds field capacity (Fig. 10f). The DC for the shallow aquifer k E2 (Fig. 10i) has proportional influence on streamflow when z is above E X , or else, inverse influence during the driest part of the recessions. Finally, the infiltration rate at field capacity (f c ), the coefficient for maximum infiltration rate (α) and the DC for non-root reservoir (k D ) seem to have no relevant effects on the streamflow modelling (Fig. 10d, e and g, respectively).

Model simplification
A model simplification was carried out by assuming that the three less sensitive parameters could be considered as modelling constants. The original 10-parameter model (identified as M 1 ) was converted to a 7-parameter model (called M 2 ) in the search for model parsimony and faster convergence. The parameters f c , α and k D were fixed according to the optimal  values from the original model. Hence, we assumed that f c = 7.45 × 10 −6 m s −1 , α = 102 and k D = 6.65 × 10 −5 s −1 . This simplified model M 2 converged after 1524 runs and reached a NS = 0.88, a performance that is almost as good as that of M 1 . The optimal values for M 2 are presented in the seventh column of Table 2, and are very close to the optimal for M 1 . Given the algorithmic configuration of Nelder-Mead procedure, the convergence of M 2 to the same parametric optimal than M 1 is not assured (since the new model structure defines a completely different simplex problem, on a smaller 7-dimensional space). From this result, an almost identical water balance is obtained from M 2 . In any case, M 2 should be validated under different conditions (other basins and spatio-temporal scales), in order to corroborate that the numerical values that we assigned to the three fixed parameters can be exported to different contexts, and can always be treated like modelling constants. The simplification here undertaken also proves the effectiveness of the time-varying sensitivity analysis here proposed, as a tool to identify the main parametric sources of modelling uncertainty and, therefore, the potentially redundant model parameters. An additional simplification step was carried out by separately removing two lateral flowpaths that could be unnecessary in the model structure (Q B2 and Q E2 ), given that reservoirs B (surface) and E (aquifer) already have primary flowpaths (Q B1 and Q E1 ). However, while suppressing Q B2 affects the simulation of the main peakflow (that decreases from 0.84 m 3 s −1 to 0.52 m 3 s −1 ), the exclusion of Q E2 enormously worsens the model's ability to shape the hydrograph recessions (and therefore the aquifer recharge/discharge processes), yielding a NS = 0.71 for the optimized model. From this exercise we conclude that these two flowpaths (and their respective parameters) play a clearly identifiable and relevant role in the model structure.  . Spearman indexes between Q and: (a) all the parameters: black horizontal lines indicate the 95% confidence interval out of which the Spearman is significant, (b) partition parameter β, c) surface discharge k B , (d) infiltration rate at field capacity f c , (e) maximum infiltration α, (f) root-reservoir discharge k C , (g) non-root discharge k D , (h) threshold for shallow aquifer E X , (i) shallow aquifer discharge k E2 , (j) deep aquifer discharge k E1 and (k) deep percolation discharge k E3 . The black dots indicate the time steps for which the Spearman is significant at the 95% confidence level, while grey dots indicate non-significant correlations.

Hydrological processes in the experimental basin
The main hydrological processes and components observed using this measuring/modelling approach are presented in the following four sub-sections.

Interception, throughfall, stemflow and surface runoff
A review of water balance (WB) partitioning in comparable situations is proposed in Table 4 (interception loss, evapotranspiration, surface/non-saturated runoff, baseflow, change in soil water content and deep percolation). In our basin, the adjusted interception loss (R In ) equalled 4% of input rainfall (R). This value is the same as found by Imbach et al. (1989) in a WB experiment under similar coffee and E. poeppigiana land cover, but is lower than other reports obtained by direct measurements at plot scale in Costa Rica: Jiménez (1986) found 16% under the same AF system; Harmand et al. (2007) found 15% for coffee and Eucalyptus deglupta (but this R In could be lower because stemflow was not separately measured); Siles (2007) found 11-15% under coffee and Inga densiflora. Interception loss can be greatly affected by the local LAI of both layers, the specific architecture of the coffee and trees and the rainfall regime. From the throughfall/stemflow component in our model (equal to 96% of R), 4% of R came out of the basin as surface runoff Q B . This value is not far from other reports at plot scale, like those byÁvila et al. (2004): 1-9% (coffee and E. deglupta); Harmand et al. (2007): 2% and Siles (2007): 3-6% (coffee and I. densiflora). At basin scale, Fujieda et al. (1997) measured 5% on a 0.56 km 2 basin under a three layer forest in Brazil, Lesack (1993) modelled 3% in a rain-forest basin in Brazil (area: 0.23 km 2 ), Kinner et al. (2004) modelled 4% in a Panamanian tropical-forest basin (0.10 km 2 ) and Charlier et al. (2008) modelled 10% in a banana-plantation basin in Guadeloupe (0.18 km 2 ). A major source of surface runoff and discrepancy between plot and basin scale studies could be the presence of roads. In our experimental basin the total length of roads is 10 km and it represents 4.5% of the basin area. Then, assuming a direct road-runoff of 80% (Ziegler et al., 2004), this component could represent up to 95% of the total basin surface runoff.

Infiltration
The non-intercepted and non-runoff fraction of incident rainfall was infiltrated (i = 92% of R). This large i/R ratio and a Q B /i ratio close to 4% give indication of very high infiltration and drainage capacities, which are typical of andic-type volcanic soils (Poulenard et al., 2001;Cattan et al., 2006) and are further enhanced for perennial crops in the absence of tillage and in presence of substantial macroporosity (Dorel et al., 2000). It was indicated by preliminary measurements carried out in our experimental basin with a Cornell infiltrometer (Ogden et al., 1997) that steady state infiltrability values could be as high as 4.7 × 10 −5 m s −1 (168 mm h −1 ) (Kinoshita, personal communication, 2009). Recent experiments found hydraulic conductivity values as high as 3.4 × 10 −5 m s −1 (122 mm h −1 ) and 2.1 × 10 −5 m s −1 (75 mm h −1 ) for andisols in Costa Rica (Cannavo et al., 2010) and Guadeloupe (Charlier et al., 2008), respectively. From our infiltration capacity (that was modelled as a function of soil water content in the nonsaturated reservoirs), we calculated the infiltration/rainfall (i/R) and the runoff/infiltration (Q B /i) ratios at the stormevent scale, analyzing 78 events with cumulative rain higher than 10 mm. No significant changes were found in any of these two ratios as a function of time (or season), and only slight reductions were detected for increasing storm cumulative rainfall or soil water content. This fact, added to the permanently high i/R ratios for each individual event (65-98%), are explained by the constantly high infiltration capacity of andisols and the stable behaviour of soil water content throughout the year. We therefore observed very efficient soil and aquifer recharge mechanisms ( Fig. 7d to f). According to Dorel et al. (2000) the soil properties of non-tilled perennial crops in andisols are mainly determined by wetting-drying cycles and by biological activity in the soil. Those are relatively stable factors in our experimental basin, given the absence of a well-defined dry season on this side of the country (Caribbean influence), and the permanently high organic matter content in these soils.

Evaporation and transpiration
Transpiration of coffee plants and trees T accounts for 20% of R (645 mm y −1 , obtained from modelling). Higher transpiration values were reported by van Kanten and Vaast (2006) in a coffee AF system under E. poeppigiana (29%, 897 mm yr −1 ), and also by Siles (2007), who measured values ranging from 31% (1008 mm yr −1 ) to 34% (905 mm yr −1 ) (see Table 4). In other experimental plots in Costa Rica, estimations of T ranged from 42% (811 mm yr −1 ) to 53% (750 mm yr −1 ) (Imbach et al., 1989;Jiménez, 1986, respectively). T can be highly dependent on local ET 0 , on the effect of drought on stomatal closure, and on LAI. For a better site comparison, we computed a simple "normalized transpiration" index NT = T (ET 0 LAI) −1 in Table 4 and found that our value (0.16) is very close to the respective ratio in the study of Siles (2007).
As mentioned earlier, we measured the actual evapotranspiration (ETR = T + E u + R In ) in our experimental basin by the eddy-covariance method and it represented 25% of R (804 mm yr −1 ). This is the smallest value reported in comparison to the other studies on coffee AF systems (Table 4). The closest value (38%, 956 mm in nine months) was measured by Harmand et al. (2007), but it is about 1.5 times our ETR; while a maximum of 69% (985 mm yr −1 ) was reported by Jiménez (1986). It must be stressed that our study is the only one that brings independent validation of ETR through energy balance closure, whereas many plot studies might carry errors due to the calibration of sapflow, the model of E u or the sampling of R In . At the basin scale, many authors modelled ETR values around 30% and 46% in tropical experimental basins. In absolute terms these values are: 1120, 554-682, 1961-2345 and 1300 mm yr −1 (Lesack, 1993;Fujieda et al., 1997;Genereux et al., 2005;Charlier et al., 2008, respectively).

Drainage, deep percolation and streamflow
Another component of interest at plot scale is drainage (vertical flow beyond root reservoir), which we modelled as 69% of R and seems, on average, slightly higher than the values reported in the literature (Table 4). The modelled deep percolation (or subsurface outflow downstream the basin outlet) was 12%, much lower than the 42% reported by Charlier et al. (2008), but similar to the 6% encountered by Kinner and Stallard (2004). As deep percolation is calculated from water balance closure, its accuracy is enhanced when Q and ETR are precisely measured. The first (or the second) highest WB basin output is usually streamflow Q, which in our experimental basin was recorded and modelled as 64% of R. It seems to be very close to similar measurements in the tropics (Table 4). The baseflow from the aquifer accounted for 87% of total Q, while surface runoff was only 6% and non-saturated runoff 7%.

The coffee agroforestry basin and hydrological services
Modelling the hydrological behaviour of this experimental, coffee AF basin gave some insights on the provision of HES by this system. This is particularly relevant in the Costa Rican context where HES payments have already been implemented as national environmental protection policies (Pagiola, 2008). Two main services related to water quality can be recognized, both linked to the observed high infiltration i (92% of R) and low surface runoff Q B (4% of R). At first, the low Q B in the basin is closely associated to low surface displacement of fertilizers, pesticides and sediments (Cattan et al., 2006Leonard and Andrieux, 1998;Bruijnzeel, 2004). We found very constant Q B /i ratios through the time, or under different rainfall intensities, which may come from the expected stability in soil hydraulic properties (e.g. high infiltration capacity) and the absence of either a marked dry season that controls soil desiccation (Park and Cameron, 2008;Dorel et al., 2000) or mechanized agricultural practices like tillage affecting soil compaction, surface roughness, continuity of pores, macroporosity, soil cover and organic matter content (Le Bissonnais et al., 2005;Chahinian et al., 2006). However, the high drainage capacity of these andisols might be a disadvantage in terms of percolation and groundwater contamination by agrochemicals (Cattan et al., 2007;Saison et al., 2008) given our model estimates for groundwater recharge of around 67% of R. In addition, the modelled surface runoff for the experimental basin includes possible discharges from unpaved roads and ditches, which needs to be controlled to avoid excessive water, sediment and contaminant flux concentrations.
A second HES might be the streamflow regulation function provided by this AF basin through aquifer recharge/discharge mechanisms. With a measured evapotranspiration close to 25% of R (which is presumably much lower than the equivalent for forests), soil and aquifer water depletion seems unlikely under the observed hydrogeological and climatic conditions, favouring water availability during dry seasons (Robinson et al., 2003;Bruijnzeel, 2004). On the other hand, during intense rainfalls and tropical storms the aquifer is efficiently recharged, as we have observed in our piezometric measurements. The result is a homogeneous seasonal distribution of streamflow, with a high rainfall recuperation fraction of 64% (Q/R).

Conclusions
This paper gives some insights into the assessment of Hydrological Environmental Services (HES) by studying the hydrological processes in a particular micro-basin. The water balance partition is proposed as a baseline for analyses and negotiations leading to the payment for HES, for which measuring and modelling approaches are complementary. The understanding of water dynamics supplied a better knowledge about the main services provided by the studied ecosystem, as well as some potential vulnerabilities.
The general behaviour of the coffee AF basin (1 km 2 ) on andisols can be summarized by the fact that 92% of R was infiltrated through the highly permeable andisol, 64% of R was measured as streamflow, 25% of R was measured as evapotranspiration, no major seasonal variation was detected in the soil water stock, and a large aquifer contribution to streamflow was observed in the shape of baseflow (Fig. 7b). These are characteristics of a system prone to generate important HES at basin scale, which is a result infrequently reported for coffee systems.
We proposed an original modelling approach coupling a hydrological and a SVAT model, calibrated using the streamflow at the outlet of the basin, but validated by independent and direct measurements of streamflow (test period), evapotranspiration, soil water content and water table level.
We presented a standard uncertainty analysis in order to built simulation confidence intervals around our modelled streamflow values, as well as a sensitivity analysis to investigate the source of such uncertainty. The first parameterization of the model is considered adequate, though model simplification could be attempted centred on the three less sensitive parameters. Special attention needs to be given to direct measurement of a representative field capacity and the associated probability distribution function.
The conceptual nature of our Hydro-SVAT model allows a wide time/space domain of application, conditional only on knowledge of some general properties of the basin of interest and on the acquisition of basic hydrological data. Different environments can be configured in terms of climate, land cover, soils and hydrogeology, and further applications under different conditions are desired to test the generality of the model. Complementary studies like hillslope and channel surface runoff, basin water losses through roads, temporal variation in soil and ecophysiological properties and ground water dynamics and composition are expected.