Uncertainty in data for hydrological ecosystem services modelling: Potential implications for estimating services and bene ﬁ ciaries for the CAZ Madagascar

This study assesses the di ﬀ erences in modelled total runo ﬀ and modelled runo ﬀ delivered to people between di ﬀ erent rainfall and population datasets in the Ankeniheny Zhamena Corridor (CAZ) of Eastern Madagascar. Runo ﬀ is estimated using the WaterWorld hydrological model driven by six rainfall datasets, and population is derived from ﬁ ve population datasets. Model results for runo ﬀ under di ﬀ erent rainfall datasets lead to variability in runo ﬀ (coe ﬃ cient of variation) up to 99% for single months and 60% in the dry season. These di ﬀ erences are much larger than di ﬀ erences in estimated runo ﬀ between baseline and complete deforestation scenario for each rainfall dataset. Population estimates for the CAZ range from 1.2 to 2 million between the population datasets. Di ﬀ erences in runo ﬀ under di ﬀ erent rainfall datasets lead to an average of 356,000 people estimated to receive 90% more runo ﬀ and nearly 750,000 people estimated to receive 50% more or less runo ﬀ relative to a baseline rainfall dataset. Therefore, the choice of rainfall data in hydrological ecosystem services modelling has a large in ﬂ uence on estimates of ecosystem service ﬂ ows highlighting the need for modellers to justify their data choices and report on uncertainties in results, particularly in light of potential policy decisions based on modelled outcomes.


Introduction
Ecosystem services (ES), or the benefits that people obtain from ecosystems, is a concept that has been widely used in science and policy making (Chaudhary et al., 2015) and is increasingly used to support decisions on natural resource management (Vorstius and Spray, 2015).To support policy applications of ecosystem services based approaches, an improved understanding of the biophysical underpinning of ecosystem processes that drive service delivery is necessary.Such an understanding requires the description and quantification of interactions of ecosystem components and their effects on a service or services (Martin-Ortega et al., 2015).
Hydrological processes have been identified as delivering ecosystem services that are fundamental to both human well-being and the maintenance of biodiversity.Freshwater ecosystems provide society with essential services of water supply for extractive uses (e.g.municipal, agricultural and industrial) as well as in situ uses such as hydropower generation and recreation.The recognition of these services and their value to society has led to considerable growth in payment for ecosystem services schemes in the past two decades (Martin-Ortega et al., 2015).
Assessments of hydrological ecosystem services (defined as the benefits to people produced by terrestrial ecosystem effects on freshwater; Brauman et al., 2007) are ideally based on long-term observational data on e.g.land-use, runoff, storage and groundwater levels as well as detailed data on water use.However, such approaches are typically time consuming and costly.Therefore, hydrological ecosystem services assessments are often carried out using modelling tools.Such tools can be used to assess current situations but can also help to better understand potential future changes due to e.g.climate or land cover change over time and space (Nelson et al., 2009;Burkhard et al., 2013) and thus make water resource decisions more effective.However, while methodological uncertainties are often addressed in these studies, very few specifically focus on uncertainties surrounding model input data that are fundamental to their application (see Schulp et al., 2014;Pagella and Sinclair, 2014 supporting this statement).Including an uncertainty assessment in ecosystem services analysis gives greater validity to findings and benefits decisions they are intended to inform (Hamel and Bryant, 2017;Wright et al., 2017).
Hydrological ecosystem services are strongly driven by climate (Mulligan et al., 2015) and hence to reliably quantify these ecosystem services, good quality climate data is necessary, particularly for rainfall inputs, the most fundamental of hydrological inputs.However, in many cases, and in particular for developing countries, good quality, long term precipitation data is not readily available from local sources (e.g.Maidment et al., 2017).This holds true for Madagascar as well (see e.g.Villa et al., 2014;Neugarten et al., 2016).While observational data from a limited number of monitoring stations may be available, there is often not enough data to create high resolution gridded rainfall fields as required by many hydrological models.This then necessitates the use of regional or global scale gridded rainfall datasets, based on collations of observed data (e.g.CRU; New et al., 2002, WorldClim;Hijmans et al., 2005) or remote sensing based products such as those based on the Tropical Rainfall Measurement Mission (TRMM; Huffman et al., 2007) or its successor, the Global Precipitation Measurement Mission (GPM, Hou et al., 2014).
The choice of rainfall data to use in hydrological modelling is often driven by considerations of availability, representativeness for the study region or model requirements in terms of temporal or spatial resolution.For this reason, many model assessments of hydrological ecosystem services only use a single rainfall dataset (e.g.Portela et al., 2012;Neugarten et al., 2016) and thus do not take into account any uncertainty in this critical input dataset, even though uncertainties in rainfall datasets have been assessed for many regions around the world (Tian and Peters-Lidard, 2010;Levy et al., 2017), including Madagascar (Awange et al., 2016).These studies revealed large relative differences between datasets at global scale but also found that rain gauge observation-based datasets generally outperform those based on remote sensing.
In this paper we assess the potential differences in the potential delivery of hydrological ecosystem services, expressed as water provision (runoff as the downstream accumulation of positive water balances on a monthly and annual basis) between six different global gridded rainfall datasets using the WaterWorld hydrological ecosystem services assessment model.While runoff is not a service per se, estimation of benefits for services like sediment retention and flood control will be affected by the same spatial variability, though the magnitude may be different due to localized dynamics specific to each service.To provide context to the differences between rainfall datasets, we also assess potential changes in runoff under a scenario of complete deforestation for each rainfall dataset.In addition, using different spatially explicit datasets on population densities, we assess the number of beneficiaries for modelled runoff under all six rainfall datasets.We apply this to the hydrological influence area (catchment) of the Ankeniheny Zahamena Corridor (CAZ), which is a newly established protected area and REDD + pilot project area in Eastern Madagascar which was the focus of a large multidisciplinary study into the feasibility of a payments for ecosystem services scheme: Can Paying for Global Ecosystem Services reduce poverty (P4GES).

Study area
Our study area fully encompasses the watersheds covering the CAZ protected area, located in eastern Madagascar (Fig. 1).The CAZ was formally granted protected area status (IUCN category VI) in 2015, linking up a number of previously existing protected areas (Poudyal et al., 2016).The CAZ measures around 381,000 ha and contains remnants of Madagascar's humid rainforest providing important ecosystem services for the surrounding area, in particular water regulation and provision for the agricultural plains on both the east and west sides of the corridor as well as sediment control for two hydroelectric plants that supply electricity for Madagascar's two largest cities.The rivers in the CAZ provide water directly to an estimated 350,000 residents located within the CAZ protected area as well as to residents in the surrounding watersheds and in the provincial capital Toamasina via a network of dams and aquifers.No official estimates for domestic, agricultural and industrial consumption exist for the CAZ but for Madagascar as a whole, 99.98% of surface water withdrawals are for agriculture with nearly all of the irrigated land used for rice production.Irrigation infrastructure is highly susceptible to flooding and siltation (Portela et al., 2012).
Remaining forest in the CAZ is mainly threatened by slash-and-burn agriculture and mining activities.Increasing population pressure has contributed to reduced fallow periods leading to reduced crop yields and further pressure on the forests (Ghimire et al., 2017) threatening the essential ecosystem services it provides.
The area has a hot and humid tropical climate with average annual rainfall of 1625 mm at the Andasibe station (mid elevation, 990 m.a.s.l.) for the period 1983-2013(Météo Madagascar, 2013).Rainfall is highly variable throughout the area, ranging between around 1000 and 3500 mm a year for Andasibe, Brick-ville and Ambatondrazaka stations.Roughly 75% of total annual rainfall is delivered in the rainy season lasting from November to April with the remainder delivered in the dry season.Cyclones bringing high winds and heavy rainfall occur occasionally in the November to May cyclone season (Tadross et al., 2008).
An estimated 1.5-2 million people live in the study area (Table 3).Around 350,000 people are estimated to live in the CAZ protected area, based on a 2010 household survey (World Bank, 2012;Portela et al., 2012).A large proportion of people in the wider study area live in the provincial capital Toamasina.
WorldClim V1 is a global 1-km resolution long-term monthly mean climate surface based on interpolation of some 50,000 observational stations with records ranging from 1950-2000 using a thin plate spline with elevation as a co-variable.Recently, this dataset has been updated using an increased number of observation stations (up to 60,000) for a more limited temporal range from 1970-2000 and including satellite derived variables in the interpolation (WorldClim V2).The TRMM rainfall climatologies processed by Mulligan (2006) are monthly mean climatologies based on analysis of the full dataset for TRMM 2B31 (combined PR, TMI profile) from 1997-2006 representing atmospheric rainfall with observations for each pixel about every ten days.A total of 50,000 satellite swaths (the area imaged on the surface of the earth at each overpass of the satellite) were used and resampled from the native resolution of 5-km to a 1-km spatial resolution.B-TRMM extends this to a 12 year series as described in Bookhagen (2018).CHPclim is a 0.05 deg (∼10 km) monthly rainfall climatology based on satellite fields, gridded physiographic indicators and two sets of long-term rainfall datasets from the FAO and the Global Historical Climate Network (GHCN) resulting in a global monthly rainfall dataset for roughly the period 1980-2009.Finally, the CHELSA dataset is a 30-arc second (∼1-km) resolution monthly global rainfall dataset for the time period 1979-2013 based on a quasi-mechanistic statistical downscaling of the ERA interim global circulation model with a Global Precipitation Climatology Centre bias correction.A summary of these datasets and the mean annual rainfall over the CAZ study area are provided in Table 1.

Runoff data
Runoff data from the Global Runoff Data Centre (GRDC) was available for four stations within the CAZ: Bac Ampitabe on the Laroka river, Rogez and Andekaleka on the Vohitra river and Ringaringa on the Ivondro river (Fig. 1).Table 2 provides more details on the length of records and mean annual flows of these rivers at these locations.Data records range between 14 to 30 years, with no data available after 1983.

Population density data
We used five different spatial population datasets.The default population dataset used by the WaterWorld model is the Landscan 2007   (Landscan 2007, Oak Ridge National Laboratory) population density dataset.These data represent spatially modelled distributions of population based on land use/land cover, high resolution imagery analysis, transportation networks, elevation, slope, etc. Landscan has the finest temporal resolution data available producing global population datasets on an annual basis.However, this dataset is proprietary and for this analysis it was not possible to obtain a more recent dataset.To be able to analyse more recent population datasets we used the 2010 and 2015 population distribution datasets developed by the Worldpop project adjusted to match UN national estimates (Worldpop 2013;Linard et al., 2012;Stevens et al., 2015).These are also modelled population distributions based on combining land cover data with official census data for Madagascar.In addition, we used two more datasets for the year 2015: the Gridded Population of the World (GPW) v4 global gridded population dataset at 1-km resolution (CIESIN, 2016) which is based on an area-weighting approach and the Global Human Settlement Layer (GHSL) population dataset 2015 (EC-JRC/CIESIN, 2015) which uses the GPW4 population estimates but disaggregates these data based on the distribution and density of built-up areas as mapped in the Global Human Settlement Layer (GHSL).A summary of these datasets and the estimated number of people in the study area according to each is provided in Table 3.

Methods
We used the WaterWorld (V2) hydrological model at 1-km spatial resolution to assess differences in modelled runoff when driven with the different rainfall datasets.WaterWorld (Mulligan, 2013a) is a fully distributed, process-based hydrological model that utilises remotely sensed and globally available datasets to support hydrological analysis and decision-making at national and local scales across the globe, with a particular focus on un-gauged and/or data-poor environments.The model includes modules for rainfall distribution based on wind interaction, fog inputs based on cloud cover and potential and actual evapotranspiration based on climate and vegetation cover (MODIS Vegetation Continuous Fields; Hansen et al., 2006;Sexton et al., 2013).Runoff is calculated as the downstream accumulation of water balance along a drainage network calculated from a digital elevation model.Simulations are carried out for four diurnal time-steps representing a daily cycle (morning, afternoon, evening and night) for one day each month, totaling 48 timesteps.The model's equations and processes are described in more detail in Mulligan and Burke (2005) and Mulligan (2013b).The model parameters are not routinely calibrated to observed flows as it is designed for hydrological scenario analysis and parameter values calibrated to existing conditions may be inappropriate for scenario situations (Mulligan, 2013a).
In order to assess which input rainfall dataset provides the closest agreement with observed runoff, modelled runoff for all six rainfall datasets was compared with long term mean observed runoff at the GRDC locations for which data were available (Table 2).Modelled runoff at the nearest location on the hydrological flow network to the reported GRDC location was compared with the GRDC data for the location, noting that they cover different time periods.
To assess the number of people affected by changes in runoff, we ran a baseline simulation using the WorldClim V1 rainfall dataset which is the default rainfall dataset in the model, keeping all other input parameters constant.We then uploaded the five other rainfall datasets to the model as scenario data to understand the change in runoff with these datasets relative to WorldClim V1.The different population density datasets were then overlaid with the runoff results by uploading them to WaterWorld as a Metric of Interest (MOI).This model feature allows for the calculation of statistics for areas identified by the MOI of any output variable, i.e. in this case the total number of people receiving more or less runoff.For pixels with scenario runoff greater or less than the baseline rainfall runoff, the population was summed to determine the number of people with water yields less or more than the baseline.Results are presented as the relative difference in runoff relative to the baseline and the number of people receiving more or less runoff.Whilst this approach does not specifically address beneficiaries in the study area, it does provide an indication of the number of people potentially receiving more or less runoff under different rainfall datasets.The differences in water received for each population data pixel depends on the hydrological influence area upstream of that pixel (cf.Servicesheds; Tallis et al., 2012;Mandle et al., 2015) which is small for most pixels but can be large for pixels located on the main streamflow network (Fig. 2).
In most hydrological ecosystem services studies, the focus is not on analysing differences between input rainfall datasets but rather on scenario simulation of changes in land use, such as deforestation.In order to provide context to the extent of the uncertainty using different rainfall datasets in hydrological ecosystem services modelling, we also ran simulations removing all tree cover in the area, a scenario which can lead to increased runoff in all seasons due to the reduced vegetative water use.The built-in scenario tool (Mulligan and Clifford, 2015;Mulligan et al., 2015) in WaterWorld was used to set up these scenarios, using a rule based approach to remove all tree cover based on MODIS 2010 fractional tree cover data (Sexton et al., 2013) and replace with bare land.Modelled runoff under scenario conditions at the GRDC gauge locations was then analysed and compared with model results for the different rainfall datasets.

Differences in rainfall for the CAZ study area
The six rainfall datasets show relatively large monthly differences on average for the CAZ study region (Fig. 3), particularly in the wet season months of January to April with a range of rainfall amounts for the month of March between 170 and 350 mm for TRMM and CHELSA respectively.As can be expected, both WorldClim datasets are very similar, as are the TRMM and B-TRMM datasets.The CHELSA dataset  consistently shows the highest rainfall while TRMM is lowest for nine out of twelve months.Spatially, the differences are more clearly visible (Fig. 4a).Broadly, all datasets are able to distinguish the relatively dry north west of the study area from the much wetter coastal region, although the CHELSA dataset shows a much wetter north east than the other datasets.Both TRMM datasets have greater ranges within the study area ranging between 241-4759 mm/yr and 399-4765 mm/yr for TRMM and B-TRMM respectively.To analyse some of the spatial variability between the various rainfall datasets, the standard deviation (SD) and coefficient of variation (CV) of annual rainfall among the six datasets are provided in Fig. 4b, showing a variability for some pixels up to 1107 mm/yr and a variability relative to the mean of up to 76%.

Impacts on runoff
Impacts on runoff were assessed for all four sub catchments for which GRDC observational runoff data was available.For all four GRDC locations, modelled monthly runoff in the dry season tends to be lower than observed (Fig. 5), while wet season runoff (Nov-March) tends to be overestimated by the model for some but not all rainfall datasets.This is likely due to the lack of a soil storage component in the model which leads to higher peak flows and reduced dry season flows supported by storage.In addition, fog inputs through cloud forest contributes considerably to the water balance in this region, in particular in the Ringaringa subcatchment.Therefore, loss of forest which has occurred over the last three decades is a plausible explanation for the overall lower modelled runoff found for this subcatchment.A validation of modelled monthly runoff using WorldClim V1 precipitation data for Bac Ampitabe shows a good fit with the observed runoff (R 2 = 0.83, Nash-Sutcliffe efficiency = 0.67, bias = −2.5%).Variability in runoff for the different rainfall datasets for all four locations (Fig. 6) is very high with coefficient of variation ranging between 9% for the month of February in the Rogez catchment and 99% for the month of June for the Ringaringa catchment.In general, variability in modelled runoff using different input datasets is higher in the dry season and lower in the wet season.However, even in the wet season, it can be as much as 60% e.g. in the month of November for Ringaringa.

Impact of deforestation on runoff
A complete deforestation scenario has the largest impact on the Ringaringa subcatchment, losing 68% of tree cover relative to the baseline while the Bac Ampitabe subcatchment has the smallest relative loss of forest cover (46%) in this scenario (Table 4).Relative differences between the baseline modelled annual runoff and the complete deforestation scenario for each rainfall dataset at the four GRDC locations are very small (Table 4).These relatively small changes are partly the result of the loss of fog inputs and increased water availability through reduced loss of water by evapotranspiration cancelling each other out.For the four subcatchments, the average loss of fog inputs amounts to 47 mm/yr (2.4% of baseline rainfall) while evapotranspiration decreases with 44 mm/yr (2% of baseline rainfall).The largest differences for the deforestation scenario are found for the Ringaringa subcatchment under WorldClim V2 and TRMM rainfall data (−3%).In the Bac Ampitabe subcatchment differences between baseline and deforested runoff are a maximum of -0.8% across all rainfall datasets.Compared by rainfall dataset, these differences are smallest for the CHELSA rainfall data (−1.5%) and greatest for the TRMM data (−2.3%)across the four subcatchments.In contrast, relative differences in modelled annual runoff between WorldClim V1 and the five other rainfall datasets for each subcatchment are much larger (Table 5) with up to 94% difference in annual runoff for the Ringaringa subcatchment between WorldClim V1 and CHELSA rainfall data.The sign of the change is also variable between different subcatchments.Relative differences in modelled runoff between WorldClim V1 and the WorldClim V2 rainfall datasets are also greater than the deforestation scenario runs for the Bac Ampitabe and RingaRinga locations, indicating greater uncertainty in modelled water yield based on different input rainfall datasets than between significant land use change scenarios in this rainfall dominated environment in which rainfall is significantly greater than evapotranspiration differences between different land covers.

Distribution of population
The different spatial distribution of rainfall in these six datasets leads to different model estimates of runoff in each subcatchment.These different runoff volumes can then have strongly varying effects on the local populations contained in these subcatchments.The spatial distribution of population in the study area differs significantly between the different population datasets available (Fig. 7).In particular the GHSL 2015 population density data which is modelled using GPW4 2015 data but limited to built-up areas has large areas without any population and some very high densities in populated places such as the city of Toamasina on the East coast which can clearly be seen on the GPW4 2015 and GHSL 2015 datasets.Whilst the other population datasets also have higher densities in city areas, their estimates are more homogeneously distributed and therefore do not stand out as clearly.The Landscan 2007 data shows higher densities around road networks as a result of the modelling approach which uses road data.

Variability in runoff to beneficiaries
Model estimates of runoff under the different rainfall datasets are highly variable in magnitude compared to a WorldClim V1 baseline runoff.The large differences in modelled runoff between rainfall datasets combined with the large variability in population distribution data mean there is high uncertainty in the estimated runoff delivered to beneficiaries (Fig. 8).As a result of this variability, populations could experience more runoff relative to the modelled baseline while others could experience less runoff.Whilst runoff is different everywhere in the study region for each precipitation dataset and thus affect all people in the population datasets to some extent, the magnitude of differences in runoff varies greatly between rainfall datasets (distribution of colours in Fig. 8).Model results for the CHELSA rainfall data show particularly large differences in estimates of water yield compared to WorldClim V1, where an average of 356,000 people across the different population datasets are estimated to receive 90% more runoff and nearly 750,000 people are estimated to receive 50% more runoff.Differences between WorldClim V1 and the TRMM rainfall data runoff estimates are also large.On average, across the population datasets 250,000 people are estimated to receive 50% less runoff under the TRMM rainfall data.The CHPCLIM rainfall dataset has very similar positive numbers as the CHELSA data whereas the two TRMM datasets show similar negative patterns (i.e. less runoff compared with the WorldClim V1 baseline).As expected the WorldClim V2 dataset shows the smallest changes compared to its previous version in terms of quantity of runoff changes.More people are estimated to receive less runoff, corresponding with the lower total annual rainfall in this dataset.
The distribution of the changes across the different population datasets is fairly uniform in the CHELSA and CHPCLIM data for the larger changes (> 10% more or less runoff) as these affect larger areas due to the larger differences in rainfall with WorldClim V1 for these datasets.The TRMM datasets have larger spatial differences with WorldClim V1 rainfall.The estimated number of beneficiaries receiving 10% more runoff under TRMM rainfall compared to WorldClim V1 varies between 2% of total population for GHSL 2015 and 8% of total population for Landscan 2007.Worldclim V2 is spatially and in magnitude very similar to WorldClim V1.Under this rainfall dataset, the estimated number of people receiving between 10% and 50% less runoff annually varies between 4% and 16% of total population across the population datasets highlighting the differences in spatial distribution between these datasets.
The number of beneficiaries estimated to receive more or less runoff does not vary much between dry and wet months for the CHELSA rainfall and Landscan 2007 population data (70% and 75% of total population respectively receive > 10% more runoff) as the changes between the baseline and this rainfall dataset are relatively large and runoff changes are assessed in relative terms.For WorldClim V2, the differences in terms of estimated number of beneficiaries receiving more or less runoff between wet and dry season months are more pronounced.In the driest month, 6% of total Landscan 2007 population are estimated to be affected by large negative changes in runoff (> 10% less runoff) while in the wettest month, this figure is much higher (17%).For the GHSL 2015 dataset these differences are 17% and 24% respectively.

Discussion
While there are several other gridded rainfall datasets available globally (e.g.CRU; New et al., 2002), and regionally (e.g.TAMSAT;Thorne et al., 2001), our study focuses on a number of high resolution (1 km), widely used or recently published datasets and is thus not fully exhaustive.However, given the spread in rainfall estimates between the datasets used it is likely that such differences also exist between other datasets.Mean annual rainfall for the gridded datasets used in our analysis ranges from 1229 mm/yr to 2528 mm/yr.Observed rainfall at the Andasibe, Brick-ville and Ambatondrazaka stations within the CAZ shows an even larger range though, between 1000 mm/yr and 3500 mm/yr which is due to the highly variable terrain.In such complex environments, estimating rainfall is challenging.Rain gauge products need a dense network of observational stations to accurately interpolate rainfall while satellite based products often underestimate rainfall in these environments (Awange et al., 2016;Sun et al., 2018) which is the case for the TRMM and B-TRMM data used in our study.
Our results show that the variability of rainfall between datasets is considerable which can lead to different absolute outcomes when assessed in an ecosystem service delivery context.For example, modelled mean annual runoff at the location of the Andekaleka dam on the Vohitra river varies between 64 m 3 /s and 127 m 3 /s for TRMM and CHELSA rainfall data respectively.Such large differences clearly have implications for absolute estimates of hydropower production by the different datasets or can have implications for dam management.Our model results for a complete deforestation scenario within the CAZ study area show relatively small differences in modelled runoff relative to those shown by the rainfall data, up to a maximum of around 3.0% whereas differences in input rainfall datasets can lead to modelled runoff differences of up to 94%.The small relative differences in runoff under a deforestation scenario are partly due to the loss of fog inputs cancelling out the increase in water availability due to reduced evapotranspiration.In other contexts or areas where fog capture is not relevant, large scale deforestation tends to lead to increased water yield.In most cases with increases more than 10% of annual rainfall (Brown et al., 2005).In our case study however, annual evapotranspiration changes under a scenario of complete deforestation are still only around 3% of annual rainfall.In such rainfall dominated environments, the absolute modelled runoff can therefore be highly sensitive to the input rainfall dataset used.Uncertainties in modelled runoff due to the selection of rainfall data have also been observed in other studies (e.g.Tuo et al., 2016;Ren et al., 2018) WaterWorld is designed to examine relative changes from baseline to scenario (for example land use change) expressed as change in runoff or percent change in runoff from baseline.This strategy (see Mulligan et al., 2015) effectively reduces the dependence on accurate absolute values of rainfall for the baseline and scenario, though there will be some effects No recent runoff data was available for the analysis, therefore the modelled runoff based on the long term means of different rainfall datasets may not be coincident with the period of observed runoff, meaning that model results and data may not be comparable.This is particularly the case for the TRMM datasets which have a much shorter record starting in the late nineties.The mean total rainfall for the more recent observational data for the three rainfall stations in the study area for the period 1983-2013 is 6% lower than the mean total for WorldClim V1.In addition, land use has been changing in these basins and the WaterWorld model uses land use and land cover data for a more recent time period (2010).Therefore, some differences in runoff response between the time period for which observational runoff data is available and the modelled period are likely as a result of differences in vegetation and climatic changes in the sub basins.Furthermore, Fig. 6.Monthly coefficient of variation (CV) in runoff between model runs with six different rainfall input datasets for four subcatchments within the CAZ study region.

Table 4
Loss of tree cover for complete deforestation scenario, and relative differences in modelled annual runoff between baseline and deforestation scenario for six rainfall datasets for four subcatchments in the CAZ study area.differences between the GRDC runoff and modelled runoff are in part likely due to the inaccuracy of the location of the GRDC monitoring stations which required manual adjustment to fit to the stream network.
The availability of more recent runoff data could potentially allow for a better assessment of model performance using the different rainfall datasets.However, the main focus of this analysis is to highlight how the different rainfall datasets and their spatial distribution can lead to differences in modelled runoff.To understand the most appropriate rainfall dataset, measured runoff from all basins in the CAZ and for the same time period as the different rainfall datasets would be required.Total population counts between years also show large differences which is not unsurprising as Madagascar has a high annual population growth of around 3% (UN/DESA, 2017).Using this growth rate to extrapolate the Landscan 2007 population for the CAZ to 2015 would result in around 1.5 million people, which is comparable to the Worldpop population dataset.The GPW4 and GHSL estimates for 2015 are much higher however (around 2 million).This is mainly the result of the underlying methods of deriving pixel-based population densities between the different datasets.Since GPW4 and GHSL both attribute people to built-up areas, adjusted to national census data, their estimates may be more correct at national scale but not for specifically defined areas such as the CAZ study area as there may not be many built-up areas.The GPW4 and GHSL datasets do have very similar total population as the latter uses GPW4 data but utilises a different method of disaggregating this to pixel level.The differences in spatial distribution of the population densities between datasets can also lead to considerable uncertainty in the volume of annual runoff received under the different rainfall datasets.Larger differences in rainfall with the WorldClim V1 baseline (e.g.CHELSA) lead to similar (large) relative numbers of people (of total population estimates) receiving more or less runoff due to these changes covering wider areas.Smaller differences in runoff (e.g.WorldClim V2) affect different (but smaller) relative numbers of people across the population datasets as these differences are more variable in the landscape.Seasonally, these differences vary as well, e.g.under WorldClim V2 between 6% and 17% receiving less runoff in the driest month for LandScan 2007 and GHSL 2015 population data respectively.These differences highlight the complexity of assessing ecosystem service delivery to beneficiaries in a spatial context.Realised hydrological ecosystem services are determined by the distribution of rainfall of population and of land cover, use and management.Whilst most studies are interested in the impacts of land cover, use and management, these impacts are mostly dependent on population and rainfall if expressed in absolute terms.
Our analysis used a relatively simple approach for estimating the number of beneficiaries affected by simply overlaying the differences in modelled runoff with the population distribution data to understand how many people are affected by increases or decreases in runoff.Therefore it only provides an indication of the number of people affected, in particular rural populations who use runoff directly for agricultural or domestic purposes which is the case for most people living in the CAZ.Impacts due to more or less water availability are hard to assess as little data exist on water use in the area.A study by Harvey et al. (2014) in three regions in Madagascar including the CAZ found that around 68% of farmers have experienced severe droughts versus 44% experiencing severe floods.Therefore, both increased and decreased flows are likely to impact on agriculture and livelihoods although flood events are linked to increased peak flows which are not captured in our monthly analysis.Even though the CAZ is generally a wet area, water shortages are common during the main rice growing season, in particular during the months of Dec and Feb-March (McHugh et al., 20002).Therefore, reduced runoff during this period is likely to have a larger impact on beneficiaries.Of course, not all beneficiaries of hydrological ecosystem services are located directly near rivers or even within the CAZ study area, for example the Andekaleka dam which also delivers energy to the capital Antananarivo.Specifically for this beneficiary, continuous water supply is critical as it determines the amount of electricity produced (Onofri et al., 2017).This is particularly important in the dry season, when low water supply frequently affects electricity production in the country (World Bank, 2016).Obtaining actual water use data by beneficiaries was beyond the scope of this study as such data is typically difficult to obtain.However, incorporating such data would provide a much better picture of the uncertainty in the estimated delivery of hydrological ecosystem services to beneficiaries under different rainfall datasets.

Conclusions
The main aim of this study is to highlight the variability in rainfall and population datasets and the potential effect the interaction of their uncertainties can have on the estimation of hydrological ecosystem services.Our study has shown that for our study region, model results for runoff under different rainfall datasets lead to differences in runoff up to 99% for single months in the dry season and up to 60% in the dry season as a whole.These differences are much greater than differences in modelled runoff between a baseline and a scenario of complete deforestation of the catchment for each single rainfall dataset.Since many hydrological ecosystem service modelling assessment studies only use a single rainfall dataset as input, any conclusions on the absolute volume of water produced and delivery to beneficiaries are thus highly dependent on the selection of rainfall data, particularly in tropical regions such as Madagascar.Model studies on hydrological ecosystem services should therefore take into account the uncertainty in this critical input parameter or focus on relative changes between baseline and scenario rather than absolute magnitude of differences.Absolute magnitudes are highly dependent on input data values whereas relative changes between baseline and scenario are much more dependent upon model logic and scenario characteristics and thus less sensitive to differences between input data values (see Mulligan, 2013a).Given that WaterWorld includes several input rainfall datasets, exploring such uncertainties can easily be done which is important for practitioners as they can then convey this uncertainty to stakeholders to support decision-making (Willcock et al., 2016).
Our study also highlights the difficulties in analysing the beneficiaries of ecosystem services, particularly in a spatial context.In order to connect ecosystem service flows to beneficiaries, information on the distribution of these beneficiaries is essential, particularly in areas where beneficiaries are directly dependent on surface water such as in our study region.Some of the best available global datasets on population distribution were analysed and these show very different total populations as well as spatial distributions.These differences between population datasets have important implications for the assessment of hydrological ecosystem services.Assessments of ecosystem service impacts of land use change such as deforestation and conversion to agriculture on water provision will therefore also affect different numbers of people depending on the population distribution dataset selected.In addition, the value (economic or otherwise) of decisions on land use management or water infrastructure will be highly dependent on the number of people affected by these decisions and will therefore need to take into account the uncertainty in these datasets.
In many cases, local or regional data may be available which could help in mitigating some of these data uncertainties.However, these datasets may not be in readily available formats or will also need spatial disaggregation to make them usable for spatial modelling which introduces uncertainties.Therefore, assessments based on readily available (global) datasets will remain necessary but will need to take into account the uncertainties in critical input datasets.It is therefore advisable for practitioners to carefully review input datasets or consider using multiple input datasets and report on the uncertainty.

Fig. 1 .
Fig. 1.CAZ study area consisting of the watersheds of the CAZ protected area (dashed green) in Eastern Madagascar showing elevation, main rivers, CAZ protected area, rainfall gauge locations and GRDC runoff gauge locations and catchments.

Fig. 2 .
Fig.2.Schematic representing potential impact of runoff on beneficiaries in population datasets (for 2 example pixels).Runoff is only shown for main streams.Pixel 1 (green serviceshed) receives runoff from a small number of pixels upstream while pixel 2 (sand serviceshed) receives runoff from a much larger area.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 4 .
Fig. 4. a) Annual rainfall (mm) and b) standard deviation (SD) and coefficient of variation (CV) across datasets for six global gridded rainfall datasets for the CAZ Eastern Madagascar.

Fig. 5 .
Fig. 5. Monthly observed and modelled runoff at four locations within the CAZ study area for six global gridded rainfall datasets and a deforestation scenario using WorldCclim V1 rainfall data.

Fig. 7 .
Fig. 7. Population distribution in the CAZ study area for five different datasets (number of people per square kilometre).White indicates no data (population zero).

Fig. 8 .
Fig. 8. Number, direction and extent of runoff beneficiaries affected under different population and rainfall datasets expressed as annual differences from WorldClim V1 rainfall baseline.Positive values represent more runoff available, negative values represent less runoff.

Table 1
Global gridded rainfall datasets used in analysis.

Table 2
Details of GRDC runoff gauges within the CAZ study area.

Table 3
Details of population distribution datasets used in analysis.

Table 5
Relative differences in modelled annual runoff between WorldClim V1 modelled runoff and five other rainfall datasets for four subcatchments in the CAZ study area.