Modelling the hydrological behaviour of a coffee agroforestry basin in Costa Rica F. Gomez-Delgado, Olivier Roupsard, G. Le Maire, S. Taugourdeau, A. Pérez, M. van Oijen, P. Vaast, Bruno Rapidel, J.M. Harmand, Marc Voltz, et al. To cite this version: F. Gomez-Delgado, Olivier Roupsard, G. Le Maire, S. Taugourdeau, A. Pérez, et al.. Modelling the hydrological behaviour of a coffee agroforestry basin in Costa Rica. Hydrology and Earth System Sciences, 2011, 15 (1), pp.369-392. ￿10.5194/hess-15-369-2011￿. ￿hal-01189535￿ HAL Id: hal-01189535 https://hal.science/hal-01189535 Submitted on 1 Sep 2015 HAL is a multi-disciplinary open access L’archive ouverte pluridisciplinaire HAL, est archive for the deposit and dissemination of sci- destinée au dépôt et à la diffusion de documents entific research documents, whether they are pub- scientifiques de niveau recherche, publiés ou non, lished or not. The documents may come from émanant des établissements d’enseignement et de teaching and research institutions in France or recherche français ou étrangers, des laboratoires abroad, or from public or private research centers. publics ou privés. Hydrol. EarthSyst. Sci.,15,369–392,2011 Hydrology and www.hydrol-earth-syst-sci.net/15/369/2011/ Earth System doi:10.5194/hess-15-369-2011 ©Author(s)2011. CCAttribution3.0License. Sciences Modelling the hydrological behaviour of a coffee agroforestry basin in Costa Rica F.Go´mez-Delgado1,2,O.Roupsard3,4,G.leMaire1,S.Taugourdeau1,A.Pe´rez4,M.vanOijen5,P.Vaast1, B.Rapidel4,6,J.M.Harmand1,M.Voltz7,J.M.Bonnefond8,P.Imbach4,andR.Moussa7 1CIRAD,UPR80,Avenued’Agropolis,34398Montpellier,France 2ICE,CENPE,10032SanJose´,CostaRica 3CIRAD,UMREco&Sols–EcologieFonctionnelle&Bioge´ochimiedesSols&Agroe´cosyste`mes (SupAgro-CIRAD-INRA-IRD),2PlaceViala,34060Montpellier,France 4CATIE,7170,30501Turrialba,CostaRica 5CEHEdinburgh,BushEstate,Penicuik,EH260QB,UK 6CIRAD,UMRSYSTEM(CIRAD-INRA-SupAgro),2PlaceViala,34060MontpellierCedex1,France 7INRA,UMRLISAH,2PlaceViala,34060Montpellier,France 8INRA,EPHYSE,BP81,33883Villenaved’Ornon,France Received: 30April2010–PublishedinHydrol. EarthSyst. Sci.Discuss.: 17May2010 Revised: 1December2010–Accepted: 20December2010–Published: 28January2011 Abstract. The profitability of hydropower in Costa Rica validation of the water balance partitioning was consistent is affected by soil erosion and sedimentation in dam reser- with the independent measurements of actual evapotranspi- voirs, which are in turn influenced by land use, infiltration ration(R2=0.79,energybalanceclosedindependently),soil andaquiferinteractionswithsurfacewater. Inordertofoster water content (R2=0.35) and water table level (R2=0.84). the provision and payment for Hydrological Environmental Eightmonthsofdatafrom2010wereusedtovalidatemod- Services (HES), a quantitative assessment of the impact of elled streamflow, resulting in a NS=0.75. An uncertainty specific land uses on the functioning of drainage-basins is analysis showed that the streamflow modelling was precise required. The present paper aims to study the water bal- for nearly every time step, while a sensitivity analysis re- ance partitioning in a volcanic coffee agroforestry micro- vealed which parameters mostly affected model precision, basin (1km2, steep slopes) in Costa Rica, as a first step to- depending on the season. It was observed that 64% of the wards evaluating sediment or contaminant loads. The main incident rainfall R flowed out of the basin as streamflow hydrologicalprocessesweremonitoredduringoneyear,us- and 25% as evapotranspiration, while the remaining 11% is ing flume, eddy-covariance flux tower, soil water profiles probablyexplainedbydeeppercolation,measurementerrors and piezometers. A new Hydro-SVAT lumped model is and/orinter-annualchangesinsoilandaquiferwaterstocks. proposed, that balances SVAT (Soil Vegetation Atmosphere Themodelindicatedaninterceptionlossequalto4%ofR,a Transfer)andbasin-reservoirroutines. Thepurposeofsuch surfacerunoffof4%andaninfiltrationcomponentof92%. a coupling was to achieve a trade-off between the expected The modelled streamflow was constituted by 87% of base- performance of ecophysiological and hydrological models, flow originating from the aquifer, 7% of subsurface non- which are often employed separately and at different spa- saturated runoff and 6% of surface runoff. Given the low tial scales, either the plot or the basin. The calibration of surfacerunoffobservedunderthecurrentphysicalconditions the model to perform streamflow yielded a Nash-Sutcliffe (andisol)andmanagementpractices(notillage,plantedtrees, (NS) coefficient equal to 0.89 for the year 2009, while the baresoilkeptbyweeding),thisagroforestrysystemonavol- canic soil demonstrated potential to provide valuable HES, suchasareducedsuperficialdisplacement-capacityforfertil- izers,pesticidesandsediments,aswellasastreamflowregu- Correspondenceto: F.Go´mez-Delgado lationfunctionprovidedbythehighlyefficientmechanisms ([email protected]) PublishedbyCopernicusPublicationsonbehalfoftheEuropeanGeosciencesUnion. 370 F.Go´mez-Delgadoetal.: ModellinghydrologicalbehaviourofcoffeeagroforestrybasininCostaRica of aquifer recharge and discharge. The proposed combina- (Charlier et al., 2008; Fujieda et al., 1997; Genereux et al., tionofexperimentationandmodellingacrossecophysiolog- 2005; Kinner and Stallard, 2004), but no coffee AF basins ical and hydrological approaches proved to be useful to ac- have been equipped so far. Some reports are available for count for the behaviour of a given basin, so that it can be coffee AF systems but at the plot level and for some par- applied to compare HES provision for different regions or ticular fluxes such as throughfall and stemflow (Siles et al., managementalternatives. 2010b), tree and coffee transpiration (Dauzat et al., 2001; vanKantenandVaast,2006),surfacerunoff(Harmandetal., 2007), energy balance and latent heat flux (Gutie´rrez et al., 1994). To our knowledge, there is no comprehensive study 1 Introduction ofthewaterbalancepartitioningofcoffeeAFsystemsatthe basinlevel,includingthebehaviouroftheaquifer. The ability of ecosystems to infiltrate rainfall, sustain Truly balanced combinations of hydrological and eco- aquifers,andavoiderosionisakeydeterminantforthepro- physiological experiments and models remain scarce, al- vision of hydrological environmental services (HES), espe- though they intrinsically carry a more realistic and compre- ciallyinthehumidtropicswheresurfacefluxescanbevery hensiverepresentationofplant,soilandaquifercomponents high (MEA, 2005). Woody plants and, in particular, agro- at plot and basin scales. Most hydrological studies at basin forestry (AF) systems associating shade trees and perennial scale use flumes for monitoring the streamflow and simply cropswithdeeprootsystems,areassumedtoenhancethese estimateevapotranspiration(ET),whichpreventsatrueveri- HESincomparisontotraditionalintensivecroppingsystems ficationofthewaterbalanceclosure. (Siles et al., 2010a; Ataroff and Monasterio, 1997; Vaast et As in the tropics we assumed that ET, including the re- al., 2005). Costa Rica is renowned as a promoter of HES evaporationofinterceptedwater(R ),isanimportantcom- In bychargingwaterusersfortheHEStheyreceivefromland ponent of the water balance, even for precipitations around owners (e.g. forest conservation), focusing on water quality 3000mmyr−1, we decided to measure it directly by eddy- (Pagiola, 2008). Hydropowerproducers, generating78%of covariance (Roupsard et al., 2006; Baldocchi and Meyers, the total electricity consumption in Costa Rica during 2008 1998; Wilson et al., 2001), choosing a 0.9km2 micro-basin (ICE, 2009), are major HES payers. Coffee is one of the embeddedinaveryhomogeneouscoffeeAFplantation. As mosttradedagriculturalcommoditiesintheworldemploying anadditionaladvantage,theeddy-covariancemethodcanbe 100 million people (Vega and Rosenquist, 2001). In Costa validated itself by closing the energy balance (Falge et al., Rica, coffee accounted for 15% of the agricultural exports 2001). in2008andcovered2%oftheterritory(SEPSA,2009). As Lumped,conceptualrainfall-streamflowmodelshavebeen coffeeplantationsarepresentinthemainbasinsusedforhy- usedinhydrologysincethe1960s(e.g.CrawfordandLins- droelectric generation in Costa Rica, the possible trade-offs ley, 1966; Cormary and Guilbot, 1969; Duan et al., 1992; ofthepaymentforHESfromhydropowerproducerstocof- Bergstro¨m, 1995; Donigan et al., 1995; Havnø et al., 1995; feefarmersbecomeevident. Negotiationforthesepayments Chahinianetal.,2005). Thesemodelsconsiderthebasinas isfacilitatedbetweenprovidersandpurchaserswhentheser- anundividedentity,anduselumpedvaluesofinputvariables vice, and/or the impact of a given practice on the provision andparameters. Forthemostpart(forareview,seeFleming, oftheservice,areclearlyevaluated. However,linksbetween 1975; Singh, 1995), they have a conceptual structure based land use, management practices (like tree planting) and hy- ontheinteractionbetweenstoragecompartments,represent- drologyinCostaRicahavenotbeenthoroughlyinvestigated ing the different processes with mathematical functions to by quantitative research (Anderson et al., 2006). There is a describethefluxesbetweenthecompartments. Mosthydro- need of both, experimentation at the basin scale in order to logical models simplify the ET component based on refer- evaluate the main hydrological processes, and of integrated ence ET routines (Allen et al., 1998) or using very empiri- modellingtounderstandthebehaviourofallwatercompart- cal,non-validatedmodelsforactualET.Forinstance,many ments,includinghiddenones(e.g.theaquifer). modelsassumethatactualevapotranspirationequalstheref- The partitioning of the water balance (WB) is a pre- erence ET, or a constant fraction of it, or a variable frac- requisite to evaluate HES such as infiltration, aquifer reg- tion that depends only on the soil water content in the root ulation capacity, erosion control and contaminants reten- zone. However, improper parameterization of the crop co- tion in coffee AF systems. Comprehensive WB stud- efficient may severely affect the parameterization of hydro- ies at basin scale, including closure verification by inde- logicalresistancesandfluxes. Inconstrast,ecophysiological pendent methods, have been carried out in the developed modelsmayoperateefficientlyatplotlevelbutmisstheparti- world and for other land covers, like those reported by tioningbetweenlateralnon-saturated(subsurface)runoffand Roberts and Harding (1996), Dawes et al. (1997), Cebal- verticaldrainage,andthedynamicsofwaterinaquifersand losandSchnabel(1998),Wilsonetal.(2001)andMaedaet rivers. ThisisamajorlimitationfortheassessmentofHES, al.(2006). Someexperimentalbasinsarelocatedinthetrop- whichismainlydesiredatthebasinscale. ics,likethoseinBrazil,CostaRica,GuadeloupeandPanama´ Hydrol. EarthSyst. Sci.,15,369–392,2011 www.hydrol-earth-syst-sci.net/15/369/2011/ F.Go´mez-Delgadoetal.: ModellinghydrologicalbehaviourofcoffeeagroforestrybasininCostaRica 371 Figure 1 Fig.1. (a)LocationofReventazo´nriverbasininCostaRica, CentralAmerica. (b)PositionofexperimentalbasininsideofReventazo´n basin.(c)The“Coffee-Flux”experimentalbasininAquiaresfarmanditsexperimentalsetuptomeasurethewaterbalancecomponents. In the present study we attempted to couple two lumped zoneandwatertablelevel)andthataredescribedinthesub- models into a new integrative one, chosen to be scalable sequentsections. Third,weproposeamulti-variablecalibra- and parsimonious: a basin reservoir model similar to the tion/validationstrategyfortheHydro-SVATmodelsowecal- CREC model (Cormary and Guilbot, 1969) and employing ibrateusingthestreamflowin2009andvalidateusingthere- theDiskinandNazimov(1995)productionfunctionaspro- mainingthreevariablesin2009andthestreamflowin2010. posedbyMoussaetal.(2007a,b),andtheSVATmodelpro- Fourth,wemakeanuncertaintyanalysistoproduceaconfi- posed by Granier et al. (1999). While the basin model was dence interval around our modelled streamflow values, and consideredappropriateforitssimplicityandcapacitytosup- a sensitivity analysis to assess from which parameters this port new routines, the SVAT model was chosen for its par- uncertainty might come. Fifth, we attempt to simplify the simony (three parameters in its basic formulation), its ro- initiallyproposedmodelbasedontheresultsofthesensitiv- bustness (uses simple soil and stand data in order to pro- ityanalysis.Finally,wediscussthemainfindingsconcerning ducemodelrunsformanyyears,avoidinghydraulicparam- thewaterbalanceinourexperimentalbasin. eters that are difficult to measure and scale up), its ability to quantify drought intensity and duration in forest stands, andforitssuccessfulpastvalidationinvariousforeststands 2 Thestudysite andclimaticconditions,includingtropicalbasins(Ruizetal., 2010). 2.1 Location,soilandclimate This paper aims to explain and model the hydrological behaviour of a coffee AF micro-basin in Costa Rica, as- The area of interest is located in Reventazo´n river basin, in sessing its infiltration capacity on andisols. The method- theCentral-CaribbeanregionofCostaRica(Fig.1aandb). ology consists of experimentation to assess the main water It lies on the slope of the Turrialba volcano (central vol- fluxesandmodellingtoreproducethebehaviourofthebasin. canic mountain range of the country) and drains to the First,wepresentthestudysiteandtheexperimentaldesign. Caribbean Sea. The Aquiares coffee farm is one of the Second, we develop a new lumped eco-hydrological model largest in Costa Rica (6.6km2), “Rainforest AllianceTM” with balanced ecophysiological/hydrological modules (that certified, 15km from CATIE (Centro Agrono´mico Tropical we called Hydro-SVAT model). This model was tailored to de Investigacio´n y Ensen˜anza). Within the Aquiares farm, the main hydrological processes that we recorded (stream- we selected the Mej´ıas creek micro-basin (Fig. 1c) for the flow, evapotranspiration, water content in the non-saturated “Coffee-Flux”experiment. Thebasinisplacedbetweenthe www.hydrol-earth-syst-sci.net/15/369/2011/ Hydrol. EarthSyst. Sci.,15,369–392,2011 372 F.Go´mez-Delgadoetal.: ModellinghydrFoliogguircea l2b ehaviourofcoffeeagroforestrybasininCostaRica coordinates −83◦4403900 and −83◦4303500 (West longitude), and between 9◦560800 and 9◦5603500 (North latitude) and is homogeneously planted with coffee (Coffea arabica L., var Caturra)onbaresoil,shadedbyfree-growingtallErythrina poeppigianatrees.Theinitialplantingdensityforcoffeewas 6300plantsha−1,withacurrentage>30years,20%canopy openness and 2.5m canopy height. It is intensively man- aged and selectively pruned (20% per year, around March). Shade trees have a density of 12.8treesha−1, with 12.3% canopy cover and 20m canopy height. The experimental basin has an area of 0.9km2, an elevation range from 1020 Fig. 2. Mean monthly rainfall at Aquiares farm station (period up to 1280ma.s.l. and a mean slope of 20%. Permanent 1973–2009) compared to the rainfall in the experimental basin in streamsextendalong5.6km,implyingadrainagedensityof 2009. 6.2kmkm−2. Theaverageslopeofthemainstreamis11%. AccordingtotheclassificationbyMora-Chinchilla(2000), the experimental basin is located along a 1.3km wide strip Rainfallandclimate: rainfallwasmonitoredat3mabove ofvolcanicavalanchedeposits,characterizedbychaoticde- groundinthemiddleof3transectsofthebasin,usingthree posits of blocks immersed in a matrix of medium-to-coarse lab-intercalibrated ARG100 tipping-bucket (R. M. Young, sand, which is the product of the collapse of the south- MI, USA) connected to CR800 dataloggers (Campbell Sci- easternslopeofTurrialbavolcano’sancientcrater. Thegen- entific, Shepshed, UK), and integrated every 10min. Other eralclassificationgivenbythegeologicalmapofCostaRica climatevariableswereloggedontopoftheeddy-fluxtower (MINAE-RECOPE,1991)describesthegeneralstratigraphy withaCR1000,every30s,integratedhalf-hourlyandusing: asshallowintrusivevolcanicrocks,andtheparticularregion Netradiation: NR-Lite(KippandZonen,Delft,TheNether- as proximal facies of modern volcanic rocks (Quaternary), lands); PPFD: Sunshine sensor BF3 (Delta-T devices Ltd, withpresenceoflavaflows,agglomerates,laharsandashes. UK); temperature and humidity: HMP45C in URS1 shel- SoilsbelongtotheorderofandisolsaccordingtotheUSDA ter (Campbell Scientific); wind-speed and direction: 03001 soil taxonomy (USDA, 1999), which are soils developing WindSentry(R.M.Young,MI,USA).Thetheoreticalevap- from volcanic ejecta, under weathering and mineral trans- otranspiration from a wet grass placed under local climate formation processes, very stable, with high organic matter conditions, ET0, wascomputedinaccordancewithAllenet contentandbiologicalactivityandverylargeinfiltrationca- al.(1998). pacities. Streamflow: a long-throated steel flume (length: 3.9m; According to Ko¨ppen-Geiger classification (Peel et al., width: 2.8m; height: 1.2m)washome-builttomeasurethe 2007), the climate is tropical humid with no dry season streamflowattheoutletoftheexperimentalbasin,torecord and strongly influenced by the climatic conditions in the up to 3m3s−1, the maximum estimated discharge for the Caribbean hillside. The mean annual rainfall in the study study period from an intensity-duration-frequency analysis. regionfortheperiod1973–2009wasestimatedas3014mm TheflumewasequippedwithaPDCR-1830pressuretrans- at the Aquiares farm station (Fig. 2). At the experimental ducer (Campbell Scientific) to record water head at gauge basintherainfallin2009(3208mm)wasclosetotheannual point (30s, 10min integration), while the rating curve was mean, but showed a monthly mean deviation of ±100mm calculated considering the geometric and hydraulic proper- around the historical regime. Mean monthly net radiation tiesoftheflumeusingWinflumesoftware(Wahletal.,2000). ranged in 2009 from 5.7 to 13.0MJm−2d−1, air tempera- Avalidationoftheratingcurvewasmadesuccessfullyusing turefrom17.0to20.8◦C,relativehumidityfrom83to91%, thesaltdilutionmethodaswellasapygmycurrentmeter. windspeed at 2m high from 0.4 to 1.6ms−1 and Penman- Soil water content: a frequency-domain-reflectometry Monteith reference evapotranspiration (Allen et al., 1998) portableprobe(FDRDiviner2000,SentekPtyLtd)wasused from1.7to3.8mmd−1. to survey 20 access tubes distributed in the three study transects to provide the mean volumetric soil water in the 2.2 Experimentalsetup basin. The sensor measures at 10cm intervals, reaching a total depth of 1.6m. A measurement campaign through the The “Coffee-Flux” experimental basin and instrument lay- 20 sites was carried out every week. The sensors were cal- out was designed to trace the main water balance compo- ibrated by digging sampling pits in the vicinity of six test nents employing spatially representative methods (Fig. 1c). tubes,toobtaintheactualvolumetricsoilwatercontentfrom It is part of the FLUXNET network for the monitoring of gravimetriccontentanddrybulkdensity. greenhouse gases of terrestrial ecosystems. The hydrologi- Evapotranspiration:theactualevapotranspirationfromthe calmeasurementswererecordedfromDecember2008upto soil, coffee plants and shade trees was measured at refer- February2010. ence height (26m) on the eddy-covariance tower, similarly Hydrol. EarthSyst. Sci.,15,369–392,2011 www.hydrol-earth-syst-sci.net/15/369/2011/ F.Go´mez-Delgadoetal.: ModellinghydrologicalbehaviourofcoffeeagroforestrybasininCostaRica 373 Table1.Measuredhydrologicvariables,recordandgapsperiods. Variable Frequencyof Measurement Datagaps Gapfillingmethod measurement period (2009) RainfallR 10min 1Jan–31Dec* No – StreamflowQ 10min 1Jan–18Jul 18Jul–23Jul Usingthecurrentmodel.Peakestimationbypeakflow/ 23Jul–31Dec* peakrainfallanalysis Measured 20Hertz 4Mar–17Jul 1Jan–4Mar Penman-Monteithmodeladjustingthecanopy evapotranspirationETR 30Aug–31Dec 17Jul–30Aug conductance Soilwatercontentθ 7to15days 2Apr–7Dec 1Jan–2Apr – Watertablelevelz 30min 2Jun–31Dec 1Jan–2Jun – *AdditionalinformationcorrespondingtotheperiodJanuary–August2010wasrecordedandusedformodelvalidationpurposes. to Roupsard et al. (2006). 3-D wind components and tem- of the actual LAI on the three LAI2000 transects. The LAI peratureweremeasuredwithaWindMastersonicanemome- for shade trees was estimated using their crown cover pro- ter (Gill Instruments, Lymington, UK) at 20Hz. H O fluc- jection(onaverage12.3%overthewholebasin)observedon 2 tuations were measured with a Li-7500 open path (LiCor, averyhighresolutionpanchromaticsatelliteimage(World- Lincoln, NE, USA). Raw data were collected and pre- View image, February 2008, 0.5m resolution). As we did processed by “Tourbillon” software (INRA-EPHYSE, Bor- not have measurements of LAI for shade trees, we consid- deaux, France) for a time-integration period of 300s, then ered this LAI in the order of magnitude of coffee LAI on a post-processed using EdiRe software (University of Edin- crown-projected basis, and therefore we multiplied the ac- burgh, UK) into half-hourly values and quality checked. A tualcoffeeLAImeasuredontransectsby1.123,toestimate validation was made by direct comparison of the measured the ecosystem LAI (tree and coffee). In order to monitor net radiation R with the sum of sensible heat flux (H) the time-course of ecosystem LAI at the basin scale, we n and latent heat flux (λE): at daily time step, this yielded combined these ground measurements with time series of H+λE=0.92 R (R2=0.93) which was considered suffi- remotely-sensedimages. WeusedtimeseriesofNormalized n cientlyaccuratetoassumethatadvectioneffectsonλEcould Difference Vegetation Index (NDVI) from Moderate Reso- be neglected here. Due to lighting and sensor breakdown, lution Imaging Spectroradiometer (MODIS) data products 45 days of data were lost between July and August 2009. MOD13Q1 and MYD13Q1 (16-Day composite data, 250m Togap-fillthemissingperiodweusedthePenman-Monteith resolution). NDVIisknowntobecorrelatedwiththegreen model,whosecanopyconductancewasadjustedusingmea- LAI if it is low, for most ecosystems (Rouse et al., 1974). suredvalues. Twenty-threeMODISpixelscoveringtheexperimentalbasin LeafAreaIndex(LAI):thecoffeelighttransmittancewas wereselected,andtheirNDVItimeseriesweredownloaded. measuredmonthlyindiffuselightconditions,forfiveringsat We filtered the raw NDVI time series according to quality differentzenitalangles(LAI2000,Li-CORCorvallis,USA), criterion given in the MODIS products, and we adjusted a along three 50m-long transects through the flux tower plot, smooth spline function on it as in Marsden et al. (2010). similarly to Roupsard et al. (2008). Effective coffee LAI, Then,alinearregressionbetweenthesmoothedNDVIofthe obtained from this light transmittance, was converted into pixel including the flux tower and ground values of actual actual LAI according to Nilson (1971), using a ratio of ef- ecosystem LAI was calibrated (R2=0.69). This regression fective to actual LAI that was estimated from a dedicated wasusedonotherpixelsofthebasin, andaveragedtohave calibration. The actual coffee LAI was measured directly theannualtime-courseofactualLAI. on a small plot by counting total leaf number of 25 coffee Water table level: four piezometric wells measuring up plants,measuringleaflengthandwidthevery20leavesand to4mdepthwerebuiltinthethreemaintransectsofstudy. usingempiricalrelationshipsbetweenleaflengthandwidth Theywereequippedwithpressuretransducers(Mini-Divers, andleafarea(LI-3100C,Li-COR)(R2>0.95). Onthesame Schlumberger Water Services) that measure and record the small plot, the effective LAI was measured with LAI2000. watertablelevelevery30min. The ratio of effective to actual LAI was then calculated on Period of measurement, data gaps and gap-filling: the thissmallplot(1.75)andwasconsideredtobeconstantwith recordinginformationisgiveninTable1forthefivehydro- time and space in the micro-basin, allowing the estimation logic variables. The frequency of measurement varies, but www.hydrol-earth-syst-sci.net/15/369/2011/ Hydrol. EarthSyst. Sci.,15,369–392,2011 374 F.Go´mez-Delgadoetal.: ModellinghydrologicalbehaviourofcoffeeagroforestrybasininCostaRica is finally calculated at the 30min time step (except for soil 3.1 Modelstructure watercontentthatisanon-continuousmeasurement). When gaps are present in the measurements, a gap filling method The model structure is presented in Fig. 3. The next three wasapplied. sections will describe the model structure according to its threemajorroutinesandfivelayers. Thefirstlayeriscalled “land cover reservoir” and separates the total rainfall into 3 Hydro-SVATlumpedmodel an intercepted loss and a joint throughfall/stemflow com- ponent. The second layer or “surface reservoir” regulates Wedesignedalumped,five-reservoir-layereco-hydrological the surface runoff. The infiltration process from the second model in order to predict the water balance (WB) partition- layeriscontrolledbythejointwatercontentatthethirdand ing (stocks and fluxes) at the scale of the entire basin. It is fourthlayers,called“non-saturatedrootreservoir”and“non- basedonthewaterbalancemodelsdevelopedbyMoussaet saturated non-root reservoir”, respectively. The evapotran- al.(2007a,b)andGranieretal.(1999),andbuilttoreproduce spiration flux is calculated at the “non-saturated root reser- themainhydrologicalprocessesmeasuredattheexperimen- voir”, while both non-saturated layers control the drainage, tal basin, which will be presented in Sect. 4. The model of thepercolationandthenon-saturatedrunoffprocesses. The Moussa et al. (2007a,b) works at the basin scale and simu- fifth and last layer is the “aquifer reservoir”, which deter- latestheecosystemevapotranspirationratherroughly,while mines the baseflow and the deep percolation. Finally, we theoneofGranieretal.(1999)worksattheplotscaleandto- will explain the sum of the total runoff and baseflow com- tallyignoresthelateralwaterfluxesthroughthesoilandthe ponentsand therouting procedure togenerate themodelled roleofthebasinaquifers. ThemainnoveltiesoftheHydro- streamflow. Let A(t), B(t), C(t), D(t) and E(t) [L] be the SVAT model with respect to the model structure of Moussa waterlevelsattimet inthefivereservoirsA,B,C,DandE, etal.(2007a,b)aretheinclusionofalandcoverreservoirto respectively (or land cover reservoir, surface reservoir, non separatetheinterceptedrainfallfromthecombinedthrough- saturated root reservoir, non-saturated no-root reservoir and fall/stemflow component, and the partition of non-saturated aquifer reservoir). Let A , B , C , D and E [L] be the X X X X X soil into two reservoirs, one with and one without roots of waterlevelscorrespondingtothemaximumholdingcapaci- plants and trees. The first of these innovations intends to tiesforthefivereservoirs. takeintoaccountthenon-negligibleinterceptionlossincof- feeAFsystems,asreportedbyJime´nez(1986),Harmandet 3.1.1 Infiltrationandactualevapotranspiration al.(2007)andSiles(2007). Thesecondnewadditiontothe model is to better represent the water dynamics in the non- (a)Infiltration saturatedsoil,giventhatonlyitsupperlayerwilllosehumid- ity by root extraction. The water balance model of Granier The infiltration process i [LT−1] occurs from the second et al. (1999) is incorporated in this superficial reservoir but layer(surfacereservoir)tothethirdone(non-saturatedroot in a simplified form, so that both, the root distribution and reservoir), and eventually to the fourth one (non-saturated thesoilporosity,arehomogeneousthroughthevertical,non- non-root reservoir) when i fills the third one. The infiltra- saturatedprofile. Hence,thewatercontentinthisreservoiris tioncapacityfi(t)[LT−1]isastatevariablethatdependson thestatevariablelinkingourtwoparentmodels. the water level in the non-saturated root reservoir, given by The modelling hypotheses governing the model architec- C(t). InadditiontoCX (maximum)wedefineCF asthewa- turewere: (a)theinterceptionlosscomponentisnotnegligi- terlevelintherootreservoiratfieldcapacity. Then,fi(t)is ble in the WB and is a function of rainfall intensity, (b) in- calculatedas(seeFig.4a): filtration is a function of the soil water content in the non- IfC(t) < C thenf (t) = f + (f − f ) C(t) C−1 (1) saturated reservoirs, (c) evapotranspiration is a significant F i 0 c 0 F component in the WB and is best described using a SVAT IfC(t) ≥ C thenf (t) = f (2) F i c model that couples evapotranspiration to root water extrac- tionfromthesoil,(d)theaquiferhasahigherdischargerate wheref0 [LT−1]isthemaximuminfiltrationcapacity(f0= above a threshold level, (e) any lateral inputs to the basins α fc) and fc [LT−1] is the infiltration rate at field capacity. system are considered negligible in comparison to rainfall Theinfiltrationi bothmodifiesanddependsonB0(t),which inputs,and(f)thereisanetwateroutflowfromthesystemas is the water availability in the second reservoir before i is deeppercolation. extracted,accordingto: ® The model was implemented using Matlab V. R2007a IfB0(t) 1t−1 < f (t) theni = B0(t) 1t−1 (3) i (TheMathWorksInc.,USA). andB(t) = 0 IfB0(t) 1t−1 ≥ f (t) theni = f (t) (4) i i andB(t) = B0(t) − f (t) 1t i Hydrol. EarthSyst. Sci.,15,369–392,2011 www.hydrol-earth-syst-sci.net/15/369/2011/ F.Go´mez-Delgadoetal.: ModellinghydrologicalbehaviourofcoffeeagroforestrybasininCostaRica 375 Figure 3 Fig.3.Thelumpedconceptualhydrologicalmodelproposedfortheexperimentalbasin. Theinfiltrationmodulecalculatestheinfiltrationi asoutput TranspirationT isobtainedbysolvingT fromtheslightly variable,usingthestatevariablesB(t),C(t)andf (t). Four modified ratio: r=T ET−1 [dimensionless] proposed by i 0 parameters(C ,C ,f andα)aredemanded. Granier et al. (1999). We substituted the original Penman X F c potentialevapotranspirationPETinthatratiobythePenman- (b)Evapotranspiration MonteithreferenceevapotranspirationET [LT−1](Allenet 0 al., 1998). While ET was calculated at each time step 1t, The evapotranspiration component ET [LT−1] acts directly 0 weestimatedr asafunctionoftherelativeextractablewater on the third layer (non-saturated root reservoir) and is the REW(t)[dimensionless],astatevariablegivenbyGranieret sum of E [LT−1] the understory and soil evaporation, and u al.(1999)as: ofT [LT−1]thetranspirationalwateruptakebyroots. REW(t) = C(t) C−1 (6) ET = E + T (5) F u TheREW(t)islinkedtothesoilwatercontentaccordingto: According to Shuttleworth and Wallace (1985), the frac- tionoftotalevapotranspirationoriginatingfromtheplantsis REW(t) = θ(t) − θr (7) close to 100% of the total evapotranspiration of the ecosys- θ − θ f r tem when LAI>3 an when the soil is not saturated at its withθ(t):volumetricsoilwatercontent[L3L−3]attimet,θ : surface,whichwasalwaysthecaseinourstudy. Wethusas- r residualsoilwatercontent[L3L−3]andθ :soilwatercontent sumed for simplicity that E , the evaporation from the soil, f u atfieldcapacity[L3L−3]. wasnil. www.hydrol-earth-syst-sci.net/15/369/2011/ Hydrol. EarthSyst. Sci.,15,369–392,2011 376 F.Go´mez-Delgadoetal.: ModellinghydroFloiggicuarleb 4eh aviourofcoffeeagroforestrybasininCostaRica The parameter REW [dimensionless] is the critical c REW(t)belowwhichthetranspirationofthesystembegins to decrease. Figure 4b shows an example of some r curves asafunctionofREW(t). Eachcurvecanbedefinedonlyby REW and the rm , a maximum value for the ratio r that c LAI dependsontheLAIofthesystemas: rm = LAILAI−1r (8) LAI X m whereLAI isthemaximummeasuredLAIduringthemod- X ellingperiodandr isaparameterindicatingthemaximum m ratioT ET−1thatcanbefoundinthissystem. Then: 0 IfREW(t) < REW thenr = rm REW(t) REW−1 (9) c LAI c IfREW (t) ≥ REW thenr = rm (10) c LAI Finally,wefindthetranspirationasT =rET .Thetotalmod- 0 elled evapotranspiration including the interception loss, can becalculatedas: ETR =E +T +R ,withR [LT−1]be- m u In In ing the intercepted/evaporated rainfall loss that will be ex- Fig.4. (a)Infiltrationfromsurfacereservoirasafunctionofsoil plained in the next section. Hence, ETRm can be directly watercontentinthenon-saturatednon-rootreservoir. (b)Theratio comparedtotheevapotranspirationthatwemeasuredatthe r=T ET−1asafunctionofrelativeextractablewaterREWandfor 0 fluxtower. differentvaluesofLAI(here,LAI1>LAI2>...>LAIi). (c)Sur- ThismoduleprovidestheevapotranspirationETasafunc- facerunofffromthesurfacereservoirasafunctionofitswatercon- tionofthestatevariableC(t),twoinputvariables(LAIand tent. (d) Baseflow from the aquifer reservoir as a function of its watercontent. ET )andthreeparameters(C ,REW andr ). 0 X c m 3.1.2 Waterbalanceinthemodelreservoirs (b)Surfacereservoir (a)Landcoverreservoir The second layer is called “surface reservoir” and acts as a sheettopsoilwithagivenroughnessandsurfacerunoffde- Thefirstlayerofthemodel,denoted“landcoverreservoir”, layingproperties. Thewaterbalanceinthissurfacereservoir represents the soil cover in the basin and controls the parti- tion of the total incident rainfall R [LT−1] into intercepted foragiveninterval1t is: (thenevaporated)rainfalllossRIn [LT−1]andthecombined B(t) = B(t − 1) + R − Q − Q − i (15) troughfall/stemflowR [LT−1]. Asimplewaterbalanceof TS B1 B2 TS this reservoir is established to calculate a proxy A0(t) [L] where RTS [LT−1] is the combined throughfall/stemflow of the final water level A(t) for each time step t, by adding componentfromthepreviouslayerandQB1andQB2[LT−1] theincidentrainfallR andsubtractingthePenmanpotential arethenon-immediateandimmediatesurfacerunoffscalcu- evapotranspirationPET[LT−1]fromtheexistinglandcover latedas: humiditylevelA(t−1): Q = k B(t) (16) B1 B A0(t) = A(t − 1) + R − PET (11) wherek [T−1]isadischargeparameter,and: WecalculatedthewaterlevelA(t)inthisreservoiraswellas B RTSandRInbydifferentiatingthreecases: IfB(t) ≤ BXthenQB2 = 0 (17) IfA0(t) ≤ 0thenA(t) = 0and (12) IfB(t) > B thenQ = [B(t) − B ] 1t−1 (18) X B2 X R = A(t − 1) + RandR = 0 In TS If0 < A0(t) < A thenA(t) = A0(t) and (13) If QB2>0 then the water level B(t) is reset to BX. The X infiltrationi [LT−1]isafunctionofthewatercontentinthe RIn = PETandRTS = 0 thirdlayeranditisthelastcomponenttobeevaluatedinthe IfA0(t) ≥ A thenA(t) = A and (14) surfacereservoir. X X This surface reservoir module calculates the water level RIn = PETandRTS = A0(t) − AX B(t) as a state variable and demands one input variable The land cover module calculates at each time t the water (RTS), and two parameters for the reservoir (BX and kB). levelA(t)asastatevariable,demandingtwoinputvariables It produces three output variables: i, QB1 and QB2. The (R andPET)andoneparameter(A ). Ityieldsthepartition twolattervariablesconstitutethesurfacerunoffinthebasin X ofRintoR andR . (Fig.4c). In TS Hydrol. EarthSyst. Sci.,15,369–392,2011 www.hydrol-earth-syst-sci.net/15/369/2011/ F.Go´mez-Delgadoetal.: ModellinghydrologicalbehaviourofcoffeeagroforestrybasininCostaRica 377 (c)Non-saturatedrootreservoir whereη[LT−1]isthetotaloutflowcapacityofthisreservoir andk [T−1]adischargeparameter. Thepartitionofηing D 1 The “non-saturated root reservoir” is the third layer of the andQ dependsontheparameterβ [dimensionless]. Then: D modelanditrepresentsasoillayerwithpresenceofrootsys- g = (1 − β) ηandQ =β η (24) temsfromtreesandplants. Thewaterbalancehereis: 1 D Thenon-rootsoilmodulecalculatesthewaterlevelD(t)as C(t) = C(t − 1) + i − ET − d − d − Q (19) 1 2 C statevariableusingtwoinputvariables(d andd )andfour 1 2 parameters (D , D , k and β). It provides three outputs where C(t) [L] is the state variable of the water level at a X F D (g ,g andQ ). given time t, i is the infiltration from the second layer, ET 1 2 D [LT−1] is the evapotranspiration, d and d [LT−1] are the 1 2 (e)Aquiferreservoir non-immediateandimmediatedrainagestothefourthlayer, respectivelyandQC[LT−1]isthenon-saturatedrunofffrom Afifthlayercalled“aquiferreservoir”representstheground- therootreservoir. water system and controls baseflow and deep percolation. There will be immediate drainage d2 if at anytime the The reservoir is composed by a shallow aquifer that acts RTS component fills the reservoir above CX. Then d2= wheneverthewaterlevelinthereservoirishigherthanEX, [C(t)−CX] 1t−1 goestothefourthlayerandC(t)isreset and by a deep aquifer with a permanent contribution. The to CX. Both non-immediate drainage d1 and non-saturated waterbalancehereis: runoff Q occur whenever C(t) is higher than the field ca- C E(t) = E(t − 1) + g + g − Q − Q − DP (25) pacitythresholdC [L]accordingto: 1 2 E1 E2 F where E(t) [L] is the state variable of the water level at ρ = [C(t) − CF] kC (20) a given time t, g and g [LT−1] are respectively the non- 1 2 whereρ [LT−1]isthetotaloutflowcapacityinthisreservoir immediateandimmediatepercolationfromthefourthlayer, and k [T−1] a discharge parameter. The partition of ρ in QE1andQE2[LT−1]arethebaseflowfromdeepandshallow C aquifers respectively (Fig. 4d), and DP [LT−1] is the deep d andQ dependsonaparameterβ [dimensionless], with 1 C percolation. 0<β<1. Then: IfE(t) ≤ E thenQ = k E(t) andQ = 0 (26) d = (1 − β) ρ andQ = β ρ (21) X E1 E1 E2 1 C IfE(t) > E thenQ = k E and (27) X E1 E1 X TherootsoilmodulecalculatesthewaterlevelC(t)asstate variable using two input variables (i and ET) and four pa- QE2 = kE2 [E(t) − EX] rameters (C , C , k and β). It provides three outputs (d , X F C 1 DP = k E (t) (28) E3 d andQ ). 2 C wherek ,k ,andk aredischargeparameterscontrolling E1 E2 E3 (d)Non-saturatednon-rootreservoir deep/shallowaquifersanddeeppercolation,respectively. This module calculatesthe waterlevel E(t) asstate vari- Thefourthlayerofthemodelisdenoted“non-saturatednon- able using two input variables (g and g ) and four param- 2 1 rootreservoir”andrepresentsasoillayerwithtotalabsence eters(E , k , k andk ), toprovidethreeoutputs(Q , X E1 E2 E3 E1 ofrootsystemsandhence,ofrootwaterextraction. Thewa- Q andDP). E2 terbalancehereisgivenby: 3.1.3 Totalrunoff,baseflowandstreamflow D(t) = D(t − 1) + d + d − Q − g − g (22) 1 2 D 1 2 Thecomponentsofsurfacerunoff,non-saturatedrunoffand where D(t) [L] is the state variable of the water level at baseflowareaddedtoobtainthetotalrunoffQ [LT−1]: a given time t, d and d [LT−1] are the non-immediate T 1 2 and immediate drainages from the third layer, respectively; QT = QB + QC + QD + QE (29) Q [LT−1] is the non-saturated runoff from the non-root D As explained in Moussa and Chahinian (2009) the stream- reservoir and g2 and g1 [LT−1] are the immediate and non- flow Q [LT−1] at the outlet of the basin is obtained by the immediatepercolationtothefifthmodellayer,respectively. routingofQ usingatransferfunction(totakeintoaccount T Immediate percolation g will be produced if at anytime 2 the water travel time). The Hayami (1951) kernel function thedrainage(d and/ord )fillsthereservoiraboveD .Then 2 1 X (an approximation of the diffusive wave equation) is devel- g =[D(t)−D ] 1t−1 moves to the aquifer reservoir and 2 X oped to obtain a unit hydrograph linear model for this pur- D(t)isresettoD . Bothnon-immediatepercolationg and X 1 pose. Thatis: non-saturatedrunoffQ occurwheneverD(t)ishigherthan D t thefieldcapacitythresholdD [L]: Z F Q(t) = Q (τ) H(t − τ) d τ (30) T η = [D(t) − D ] k (23) F D 0 www.hydrol-earth-syst-sci.net/15/369/2011/ Hydrol. EarthSyst. Sci.,15,369–392,2011