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

25 Pages·2017·4.07 MB·English

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