Technical note : A hydrological routing scheme for the Ecosystem Demography model ( ED 2 + R ) tested in the Tapajós River basin in the Brazilian Amazon

Land surface models are excellent tools for studying how climate change and land use affect surface hydrology. However, in order to assess the impacts of Earth processes on river flows, simulated changes in runoff need to be routed through the landscape. In this technical note, we describe the integration of the Ecosystem Demography (ED2) model with a hydrological routing scheme. The purpose of the study was to create a tool capable of incorporating to hydrological predictions the terrestrial ecosystem responses to climate, carbon dioxide, and land-use change, as simulated with terrestrial biosphere models. The resulting ED2+R model calculates the lateral routing of surface and subsurface runoff resulting from the terrestrial biosphere models’ vertical water balance in order to determine spatiotemporal patterns of river flows within the simulated region. We evaluated the ED2+R model in the Tapajós, a 476 674 km2 river basin in the southeastern Amazon, Brazil. The results showed that the integration of ED2 with the lateral routing scheme results in an adequate representation (Nash–Sutcliffe efficiency up to 0.76, Kling–Gupta efficiency up to 0.86, Pearson’s R up to 0.88, and volume ratio up to 1.06) of daily to decadal river flow dynamics in the Tapajós. These results are a consistent step forward with respect to the “no river representation” common among terrestrial biosphere models, such as the initial version of ED2.


Introduction
Understanding the impacts of deforestation (e.g., Lejeune et al., 2015;Medvigy et al., 2011;Andréassian, 2004) and climate change (e.g., Jiménez-Cisneros et al., 2014) on the Earth's water cycle has been a topic of substantial interest in recent years given its potential implications for ecosystems and society (e.g., Wohl et al., 2012;Brown et al., 2005).Analyses of climate change impacts on the Earth's water cycle increasingly use terrestrial biosphere models, which are capable of estimating changes in the vertical water balance as a function of climate forcing and/or land-use-induced changes in canopy structure and composition (Zulkafli et al., 2013).Terrestrial biosphere models actively used for hydrological and Earth system sciences include the Joint UK Land Environment Simulator (JULES) (Best et al., 2011;Clark et al., 2011), the Community Land Model (CLM) (Lawrence et al., 2011;Oleson et al., 2010), the Lund-Potsdam-Jena (LPJ) land model (Gerten et al., 2004;Sitch et al., 2003), the Max Published by Copernicus Publications on behalf of the European Geosciences Union.F. F. Pereira et al.: A hydrological routing scheme for the ED2 model Planck Institute MPI-JSBACH model (Vamborg et al., 2011;Raddatz et al., 2007), and the Integrated Biosphere Simulator (IBIS) (Kucharik et al., 2000).
Initial formulations of the hydrological processes within terrestrial biosphere models were based on simple "bucket" model formulations (Cox et al., 1999after Carson, 1982).Moisture within each climatological grid cell of the domain was simulated in a single below-ground pool in which surface temperature and specific soil moisture factors determined evaporation, while runoff was equal to the bucket overflow (Cox et al., 1999;Carson, 1982).Recently, the hydrologic schemes within terrestrial biosphere models have become increasingly sophisticated.In the most recent generation of land surface models, water fluxes in and out of the soil column are vertically resolved and take into account feedback from the different components, for instance, through an explicit formulation of the soil-plant-atmosphere continuum.This enables the models to provide a detailed representation of the interactions between evapotranspiration, soil moisture, and runoff (Clark et al., 2015).
To couple the calculation of the one-dimensional water balance with the estimation of daily river flows, it is necessary to simulate multiple hydrological dynamics involved in the lateral flow propagation through the landscape, ideally including the most complex hydraulic features of floodplains, lakes, and wetlands (Yamazaki et al., 2011).The first step towards representing the finer-scale hydrodynamic processes responsible for patterns in river gauge observations is to consider the topographic and geomorphological features that control water flow (Arora et al., 1999).The coarse spatial resolution of regional land surface models, imposed by computational constraints, does not allow for proper simulation of the complex hydrological dynamics determined by finescale topography in river channels and floodplains (Yamazaki et al., 2011;Kauffeldt et al., 2016).However, the combination of the terrestrial models with routing schemes can be used to simulate the implications of global and regional environmental changes for flood/drought forecasting, water resources planning and management, and infrastructure development (Andersson et al., 2015).Consequently, several terrestrial biosphere models have been integrated with routing schemes.For example, JULES has been integrated with the Total Runoff Integrating Pathways (TRIP) to evaluate the accuracy of its estimates of annual streamflow (Oki et al., 1999).This integrated model was used to investigate the status of the global water budget (Oki et al., 2001).Rost et al. (2008) also used a modeling framework composed of the global dynamic vegetation model, LPJ, and a simple water balance model to quantify the global consumption of water for rain-fed and irrigated agriculture.An offline coupling of the dynamic vegetation model, ISIS, and HYDRA -which simulates the lateral transport of water through rivers, lakes, and wetlands -was proposed in Coe et al. (2008) with the purpose of reproducing linkages between land use, hydrology, and climate.Moreover, Liang et al. (1994) developed and tested the coupling of the well-known VIC model with a general circulation model (GCM) to improve the GCM's ability to capture the interactions between surface hydrology and atmosphere.For the same purpose, the MPI hydrological discharge model was validated with NCEP reanalysis and parameterized for simulating the river routing for climate analysis at global scale (Hagemann and Gates, 2001;Hagemann and Dumenil, 1997).Several routing schemes have been designed to date, including normal depth, modified pulse, simple Muskingum, and Muskingum-Cunge (USACE, 1991).Most notably, the semi-distributed kinematic wave-routing Muskingum-Cunge method has been recognized for its stability over different spatial and temporal modeling resolutions (USACE, 1991;Miller and Cunge, 1975;Cunge, 1969), and has been adopted by the most widely used regional-scale hydrological models, such as VIC, SWAT, and MGB-IPH.
Recent studies have investigated the influence of land use on regional patterns of rainfall and biosphere temperature (Ostberg et al., 2015;Bahn et al., 2014;Pearson et al., 2013).These studies tracked how the occurrence of conversion of land from its natural state over the same time frame as observed fluctuations of rainfall and air temperature occurredaspects fully analyzed by terrestrial biosphere models (Hurtt et al., 2006;Goldewijk, 2001;Ramankutty and Foley, 1999).However, these models assumed that global and regional changes in the biosphere were a result of dynamics of vegetation in a collection of landscapes given by forests, deserts, and farmland only.Inland surface waters (e.g., rivers, lakes, and wetlands) were not considered as an interactive component of the biosphere and hence the climate system (Cole et al., 2007).
The Ecosystem Demography (ED2) is a terrestrial biosphere model that simulates the coupled water, carbon, and energy dynamics of terrestrial land surfaces (Longo, 2014;Medvigy et al., 2009;Moorcroft et al., 2001) to describe the coupled water, carbon, and energy dynamics of heterogeneous landscapes (Hurtt et al., 2013;Medvigy et al., 2009;Moorcroft et al., 2001).ED2's ability to incorporate sub-grid-scale ecosystem heterogeneity arising from landuse change makes the model suited for investigating how the combined impacts of changes in climate, atmospheric carbon dioxide concentrations, and land cover affect terrestrial ecosystems.For example, ED2 was successfully used to simulate the carbon flux dynamics in the North American continent (Hurtt et al., 2002;Albani et al., 2006) and to assess the impacts on Amazonian ecosystems of changes in climate, atmospheric carbon dioxide, and land use (Zhang et al., 2015).Moreover, ED2, coupled with a regional atmospheric circulation component, has also been successfully applied to assess the impacts of deforestation on the Amazonian climate (Knox et al., 2015;Swann et al., 2015).The aforementioned studies were not aimed at assessing hydrological implications of changes in land use and climate.These works demonstrated the validity of ED2 for assessing impacts of global and regional changes on ecosystem function and built the foundations for an integrated tool aimed at analyzing hydrological implications.
In this technical note, we describe the integration of ED2 with a hydrological routing scheme.The hydrological routing scheme chosen was adapted from the MGB-IPH (Collischonn et al., 2007).This exercise aims to calculate the lateral propagation and attenuation of the surface and subsurface runoff resulting from the vertical balance calculations in order to simulate daily river flows through a large river basin.The advantage of the proposed model is its ability to predict the sensitivity of river flows to global and regional environmental changes such as climate and land-use changes.The new product combines the advantages of biosphere and hydrological models, bringing together global-, regional-, and local-scale hydrological dynamics in a single modeling framework.The resulting model is intended to be used in future studies as a computational tool to explore a variety of research questions.In particular, it could be used to analyze how current and future climate and land cover affect water availability in river systems; how land-use-driven changes can influence the water availability for human activities (hydropower, food production, urban supply); and what the implications of those changes are for water and land resources management.
The identified research areas are in line with key problems raised in the literature, focusing on the importance of largescale modeling and remote sensing to fill knowledge gaps in water resources and hydrological dynamics (Alsdorf et al., 2007;Prigent et al., 2007).The product obtained from this exercise was tested in the Tapajós Basin, a large river system in the southeastern Amazon, Brazil.
2 Ecosystem Demography (ED2) model ED2 is a terrestrial biosphere simulation model capable of representing biological and physical processes driving the dynamics of ecosystems as a function of climate and soil properties.Rather than using a conventional "ecosystem as big-leaf" assumption, ED2 is formulated at the scale of functional and age groups of plants.Ecosystem-scale dynamics and fluxes are calculated through a scaling procedure which reproduces the macroscopic behavior of the ecosystem within each climatological grid cell.It simulates ecosystem structure and dynamics as well as the corresponding carbon, energy, and water fluxes (Fig. 1; Hurtt et al., 2013;Med-Figure 2. Schematic representation of the connection between the terrestrial biosphere model and the hydrological routing scheme.Calibrating parameters are circled in red.The reservoirs are used to determine the contribution of streamflow that comes from overland flow, interflow, and groundwater flow.The daily sum of these three reservoirs is then moved from each grid cell into the drainage network.vigy et al., 2009;Moorcroft et al., 2001).ED2 simulates the dynamics of different plant functional types subdivided into tiles with a homogeneous canopy (Swann et al., 2015;Medvigy et al., 2009).The dynamic tiles represent the sub-gridscale heterogeneity in ecosystem composition within each cell.Grid cell size is determined by the resolution of meteorological forcing and soil characteristics data, typically from 1 to 0.001 • (∼ 110 to 1 km).ED2 simulates biosphere dynamics by taking into consideration natural disturbances, such as forest fires and plant mortality due to changing environmental conditions, as well as human-caused disturbances, such as deforestation and forest harvesting (Medvigy et al., 2009;Albani et al., 2006).Disturbances are expressed in the model as annual transitions between primary vegetation, secondary vegetation, and agriculture (cropland and pasture) (Albani et al., 2006).Natural disturbance, such as wildfire, is represented in the model by the transition from primary vegetation (forest in the case of the Amazon) to grasslandshrubland, and subsequently to secondary vegetation (forest regrowth); the abandonment of an agricultural area is represented with the conversion from grassland to secondary vegetation, while forest logging is represented by the transition from primary or secondary vegetation to grassland.The model is composed of several modules operating at multiple temporal and spatial scales, including plant mortality, plant growth, phenology, biodiversity, soil biogeochemistry, disturbance, and hydrology (Longo, 2014;Medvigy et al., 2009).A selection of the main parameters and the input used for this study are presented in Table 1, and for a more complete description of the model, we refer the reader to the literature available (Zhang et al., 2015;Longo 2014;Kim et al., 2012;Medvigy et al., 2009;Moorcroft et al., 2001).

ED2 hydrology module
The hydrological module of the ED2 model is derived from the Land Ecosystem-Atmospheric Feedback model (LEAF-2) (Walko et al., 2000).The model computes the water cycle through vegetation, air canopy space, and soils, yielding daily estimates of subsurface and surface runoff from each grid cell, isolated from the others in the domain.The number of soil layers and their thickness influence the accuracy with which the model is able to represent the gradients near the surface.Soil composition was derived from Quesada et al. (2010) and from the IGBP-DIS global soil data (Global Soil Data Task, 2014).As described in Zhang et al. (2015), the mean fraction values of sand and clay were assigned to each grid cell at 1 km resolution and then aggregated at 1 • resolution.Due to limited data availability, soils were assumed to be homogeneous for a depth of 6 m.Hydraulic conductivity of the soil layers is a function of soil texture and moisture (Longo, 2014).Groundwater exchange is a function of hydraulic conductivity, soil temperature, and terrain topography.Water percolation is limited to the bottom layer by the subsurface drainage, determining the bottom boundary conditions.Vegetation historical records and land-use transitions were derived from the Global Land Use Dataset (Hurtt et al., 2006).A more detailed description of the hydrological subcomponent of the ED2 model is available in Longo (2014).

ED2 runoff routing scheme (ED2+R)
River routing schemes are often used to compute the lateral movement of water over land in hydrological models for large river basins.In this way, the prediction performance of models can be evaluated using river discharge measurements.The use of routing schemes was then extended to Earth system models in order to capture the impacts of man-made structures (e.g., dams and reservoirs) and floodplain wetlands on the climate system (Li et al., 2011;Yamazaki et al., 2011).
Daily runoff estimates from ED2 were computed for specific grid cells independently.A hydrological routing scheme  Zhang et al., 2015;Longo, 2014;Knox, 2012).Additional information about ED2 parameter calibration for the Amazon Basin are available in Zhang et al. (2015) and Longo (2014) to the values for early-, mid-, and late-successional cohorts, respectively.* * * The residence time parameters are dimensionless and used to correct the Kirpich formula for time of concentration as explained in Collischonn et al. (2007).Their magnitude is influenced by the size of the grid cell and its topography.
was then linked to this model in order to estimate flow attenuation and accumulation as water moved through the landscape.The hydrological routing scheme chosen was adapted from the original formulation of the MGB-IPH, a rainfall-runoff model that has been used extensively in large river basins in South America (Collischonn et al., 2007).This model was later developed using hydrodynamic solutions and floodplain coupling (Pontes et al., 2015;Paiva et al., 2013b).Although the later development increased the modeling capabilities of the MGB-IPH in representing finescale dynamics, given the regional application of our tool, for ED2+R we decided to use the typical application of the MGB-IPH characterized by the Muskingum-Cunge approach.The original MGB-IPH model is composed of four different sub-models: soil water balance, evapotranspiration, intra-cell flow propagation, and inter-cell routing through the river network.In the present study, only the catchment and river routing methods were utilized.The resulting ED2+R model computes the daily total volume of water passing through any given grid cell in the resulting drainage network in two separate steps: first, ED2 estimates of daily surface and subsurface runoff from each grid cell are divided into three linear reservoirs with different residence times to represent overland flow, interflow, and subsurface water flow (Fig. 2).The reservoirs are used to determine the contribution and attenuation of river flow by different soil layers, characterized by different routing times.The sum of overland flow, interflow, and subsurface water flow is then moved from each grid cell into the drainage network, designed in the preprocessing phase using data from a digital elevation model (DEM) from the Shuttle Radar Topography Mission (SRTM -USGS, 2016) at a 90 m resolution and the Cell Outlet Tracing with an Area Threshold algorithm (COTAT) (Reed, 2003).Each DEM grid cell therefore becomes part of a flow path, which then accumulates water to a final downstream drainage network outlet.A complete description of the technique for defining drainage networks from DEMs employed in this study can be found in Paz et al. (2006).Once water reaches the drainage network, ED2+R adopts the Muskingum-Cunge numerical scheme for the solution of the kinematic wave equation, which also accounts for flow attenuation, using a finite-difference method as a function of river length, width, depth and roughness, as well as terrain elevation slope (Collischonn et al., 2007;Reed, 2003).Statistical relationships for the river morphology were obtained as a function of the drainage area based on geomorphic data collected by Brazil's National Water Agency (ANA) and the Observation Service for the geodynamical, hydrological, and biogeochemical control of erosion/alteration and material transport in the Amazon Basin (HYBAM) at several gauging stations in the Amazon and Tocantins basins as presented by Coe et al. (2008).Further studies successfully derived geomorphological relations in order to estimate river geometric parameters and carry out hydrodynamic simulations of the Amazon River system using a similar approach (Paiva et al., 2011(Paiva et al., , 2013)).Multiple groups of grid cells with common hydrological features, or hydrological response units, can be created in order to parameterize and calibrate ED2+R.In our approach, hydrological traits associated with soil and land cover are primarily computed in ED2; thus, we calibrated ED2+R at the sub-basin level as delineated based on the DEM.Details about the calibration procedure are provided in the next section.
The model's performance was calculated through the adoption of widely used indicators: (Pearson, 1895), calculated as in Eq. ( 1): where sim and obs are the simulated and observed time series, while n is the number of time steps of the simulation period; volume ratio, calculated as ratio of the simulated (sim) and observed (obs) total water volume in the simulation period without consideration for the seasonal distribution of flow, as in Eq. ( 2): the Nash-Sutcliffe efficiency (NSE) coefficient (Nash and Sutcliffe, 1970), calculated as in Eq. ( 3): where obs i and sim i are the observed and simulated data at time i, obs i is the mean of the observed data, and n is number of time steps of the simulation period; the Kling-Gupta efficiency (KGE) index, both 2009 and 2012 versions, calculated as in Eq. ( 4): where s values are scaling factors (set to 1 in this case), r is Pearson's correlation coefficient, β is the ratio between the mean of the observed values and the mean of the simulated values, and vr is the variability ratio, defined as vr 2009 (simulated vs. observed standard deviation ratio; Eq. 5) for the 2009 method and vr 2012 (ratio of coefficient of variation of simulated and coefficient of variation of observed values; Eq. 6) for the 2012 method (Kling et al., 2012;Gupta et al., 2009).
The optimal value for the Pearson's R, VR, NSE, and KGE indexes is 1; the closer the indexes are to this value, the more accurately the model reproduces the observed values.
Missing observations in the river flow records (HYBAM and ANA) were filled via linear spatial and temporal interpolation between the series in neighboring gauge stations (Eq.7): where z, y, and q are three gauge stations with time series highly correlated (Pearson's R ≥ 0.85), and t expresses time in days.The estimated β coefficients in Eq. ( 7) were used for the estimation of the missing observations in site y (Table 2).The interpolation of the gauge historical records was necessary to have continuous time series with a sufficient number of observations to calibrate and validate the ED2+R application in the basin.
For the presentation of the results, in order to compare the simulated and observed values, we also used flow duration curves (FDCs).FDCs are cumulative frequency plots that show the percentage of simulations steps (days in the case presented in this study) in which the discharge is likely to equal or exceed a specific value, without taking into consideration the sequence of the occurrence.

Case study: Tapajós River basin
The ED2+R formulation was parameterized and evaluated for the Tapajós River basin, the fifth largest tributary of the Amazon.This basin drains an area of 476 674 km 2 in the southeastern Amazon, within the Brazilian states of Mato Table 2. Statistics about the gauge information filling procedure (correlation with the station to be filled, number of original observations, and filled number of observations).

Sub-basin name
Main river gauge station -z in Eq. ( 7) Original number of daily gauge records (number of daily observations) Gap-filling station 1 -q in Eq. ( 7) (correlation with z) Gap-filling station 2 -y in Eq. ( 7  Grosso, Pará, and Amazonas.The main rivers in the basin are the Tapajós (with a length greater than 1800 km and average discharge of 11 800 m 3 −1 ), Juruena (length of approximately 1000 km and discharge of 4700 m 3 s −1 ), and Teles Pires (also known by the name São Manoel, about 1600 km long and average discharge of 3700 m 3 s −1 ).The river system flows northwards, with terrain elevation ranging from about 800 m a.s.l.(above sea level) in the southern region, to a few meters above sea level at its confluence with the Amazon River (ANA, 2011).The basin ecosystems are mainly represented by tropical evergreen rainforests in the northern region (in the states of Amazonas and Pará), and Cerrado dry vegetation in the south (Mato Grosso).Precipitation ranges from about 1500 mm year −1 in the headwaters (southern region) to about 2900 mm year −1 towards the basin's outlet (Fig. 3a, b).Rainfall temporal distribution is characterized by a clear seasonal distinction; total precipitation in the wet season (September to May) could be as high as 400 mm month −1 in the most tropical areas, whereas in the dry season (June to August), precipitation is close to zero in the Cerrado and as low as 50 mm month −1 in the wetter areas (Mohor et al., 2015).As a result of the large rainfall seasonal variability, river flows are also extremely variable: the mean monthly flow of the Tapajós River ranges between about 2300 and 28 600 m 3 s −1 according to the historical records used for the calibration of the ED2+R model.Soils vary from those typically seen in the Brazilian shield in the south of the basin to alluvial sediments in the north.Land use, almost completely represented by primary forest until the 1970s, was radically changed in recent decades.As estimated from the land-use/land-cover dataset used in this study (Hurtt et al., 2006), in the late 2000s only about 56 % of the basin (270 000 km 2 ) was covered by the original vegetation cover.Large parts of the basin lying in the territory of Mato Grosso were cleared to make room for agricultural and livestock production, while vast areas around the border between the state of Pará and Mato Grosso were cleared for cattle production.The northern portion of the basin is largely protected by natural parks or indigenous lands, but significant deforestation hotspots could be identified around the cities of Santarém and Itaituba and along the main transportation routes (Fig. 3c).For a more detailed description of the basin's physical characteristics and historical analysis of trends in deforestation, precipitation, and discharge, we refer the reader to Arias et al. (2017) and Farinosi et al. (2017).
For calibration purposes, the basin was divided into seven sub-basins, each of them with a corresponding gauge for which historical daily river flow observations were available (Fig. 4a).The domain was gridded with a spatial resolution of 0.5 • by 0.5 • , roughly corresponding to 55 km by 55 km.Simulations were carried out for the period 1970-2008.The ED2 model was forced using reconstructed climate (Sheffield et al., 2006) and land-use/land-cover data (Hurtt et al., 2006;Soares-Filho et al., 2006) at 1 • spatial resolution.The original meteorological dataset has a 3 h temporal resolution, which was downscaled to an hourly resolution, as described in Zhang et al. (2015).In this technical note, we describe the calibration of the flow-routing component of ED2+R.The parameterization of the ED2 terrestrial biosphere model was developed and evaluated independently using eddy-flux tower observations of carbon, water, and energy fluxes and forest inventory observations of above-ground biomass dynamics.Further details are available in Zhang et al. (2015) and Longo (2014).

ED2+R model calibration
The ED2+R model was manually calibrated using gauge observations (ANA, 2016; Observation Service SO HYBAM, 2016) spanning a period of 17 years from 1976 to 1992 (the period 1970-1975 was not considered in order to avoid simulation initiation effects) through a two-step procedure, as highlighted in Fig. 2. The first step is partitioning the flows from the two reservoirs (surface and subsurface) of the ED2 biosphere model into the three reservoirs (surface, intermediate, base) of the ED2+R routed biosphere model (parameters α and β in Fig. 2).In particular, α (ranging from 0 to 1, or from 0 to 100 %) represents the share of ED2 surface runoff allocated to the ED2+R surface reservoir.The remaining part (1 − α) is allocated to the ED2+R intermediate reservoir.β represents a similar partitioning coefficient for the ED2 subsurface reservoir to the ED2+R intermediate and base reservoirs.The second step relates to the adjustment of the residence times of the water flows in the three reservoirs for each of the grid cells in each of the sub-basins (overland, intermediate, and subsurface water flows -represented by the adjustment parameters CS, CI, and CB in Fig. 2).
In the first step, following the methodology described by Anderson (2002), the sensitivity of the α and β parameters was tested by running the model multiple times (∼ 35).For each run, the NSE (Nash and Sutcliffe, 1970) was quantified by comparing the results of the simulation to historical flow observations.The combinations of the α and β parameters characterized by the largest NSE were selected.Parameters α and β were assumed to be uniform for the whole basin.Figure 5 shows the different combinations of the α and β parameters introduced in Fig. 2. The color bar indicates the NSE resulting from the comparison between the simulated and observed river flow values obtained using different combinations of the parameters α (x axis) and β (y axis).The chosen combination (indicated by an x in Fig. 5) lies in one of the optimal combination areas (NSE ∼ 0.8).
In the second step, the residence times (τ ) of flow within the ED2+R reservoirs of each grid cell in the domain were calibrated through the adjustment of the non-dimensional parameters (CS, CI, and CB in Fig. 2) used to correct the Kirpich formula for time of concentration (as explained in Col-      , 2007).The calibration procedure characterizing the second step is similar to the previous one but in this case the calibration is repeated for each sub-basin sequentially.The calibration process was conducted from the furthest upstream sub-basins -headwaters -to the final outlet of the basin (Anderson, 2002).The model was run multiple times (between 30 and 50 per sub-basin) with different combinations of the three parameters (CS, CI, and CB in Fig. 2); for each run, the goodness of fit was quantified.This allowed us to design a sensitivity curve of the model to different combinations of the three parameters for each of the seven subbasins and to select the combination that best approaches the historical observations.Figure 6 shows how the model is sensitive to marginal variation in initial conditions of baseflow, particularly in the upstream section (i.e., UTP -Upper Teles Pires, UJ -Upper Juruena, and LTP -Lower Teles Pires).Changes in initial subsurface water were controlled by the 5-year initialization period; thus, contributions the downstream part of the basin had minimal impact (i.e., UT and LT -Upper and Lower Tapajós).Figure 7 describes the calibration of the residence time adjustment parameters for each of the sub-basins, as well as an approximate calculation of the corresponding time of concentration for each of the reservoirs in the cell.The different combinations of the values assigned to the parameters CS, CI, and CB significantly affect the overall goodness of fit of the river flow simulations (NSE indicator).The calibration process was conducted from the furthest upstream sub-basins -headwaters -(UTP -Upper Teles Pires, UJ -Upper Juruena, and JA -Jamanxim) to the final outlet of the basin (LT -Lower Tapajós).The different combinations are marked with the corresponding NSE value; the optimal combination is marked in red (Fig. 7).

Results
The integration of the routing scheme with ED2 increases the ability of the model to reproduce the observed temporal variations in river flows at the basin outlet (Fig. 8).This statement applies to all of the sub-basins, as the application of the routing scheme improved the model's performance between simulated and observed values with respect to all four measures selected (NSE, KGE, Pearson's R correlation, and volume ratio) (Table 3).Both routed (ED2+R) and non-routed (ED2) simulation results manage to reproduce the observed water availability (quantity of water available) in the basin in terms of volume.The volume ratio at the furthest downstream sub-basin (Lower Tapajós), in fact, ranges around the optimal value for both validation and calibration periods (ED2: 1.11-1.13;ED2+R: 1.06-1.13).The routing scheme improves the ability of the model to reproduce the spatiotemporal distribution of water flows across the basin: both the NSE and the KGE indexes reached values ranging between 0.76 and 0.86 in the calibration, and between 0.68 and 0.80 in the validation period (Table 3).Also, the correlation values confirm the results of the other indexes, reaching 0.88 for the calibration and 0.86 for the validation period.The performance of the presented tool is evident also when analyzing FDCs (Fig. 9a-g).The adoption of the river routing scheme allows a more realistic representations of the high discharge values (flow equaled or exceeded 0 to 20-30 % of the time), and low discharge values (flow equaled or exceeded 60 to 100 % of the time) in all the sections of the basin (Fig. 9).The model's performance in simulating river flows is generally Hydrol.Earth Syst.Sci., 21, 4629-4648, 2017 www.hydrol-earth-syst-sci.net/21/4629/2017/  3).FDCs, representing the probability of the flow values to exceed a specific discharge, highlight the positive effect of the application of the routing scheme in ED2+R across the entire range of flow variability (Fig. 9).The simulated FDCs follow the same shape of the observed ones in the furthest upstream sub-basins, especially in the cases of the Upper Juruena and Upper Teles Pires, implying that the routing scheme is effective in maintaining the simulated discharge range (Upper Juruena: 1200-2480 m 3 s −1 , Upper Teles Pires: 393-4130 m 3 s −1 ) in line with the observations (1030-2400 and 302-2767 m 3 s −1 , respectively).This is especially true for the lowest flows, where the error between simulated and observed curves is lower than 15 % (Figs.9a, b, A1).Regarding the intermediate sub-basins, Lower Juruena and Lower Teles Pires, flood duration curves show that the model overestimates the lowest values of the distribution by approximately 30 % of the observed values (flow equaled or exceeded 60 to 100 % of the time in Fig. 9c, d).Similar overestimation of the model could be noticed in the furthest downstream subbasins, Upper and Lower Tapajós (Fig. 9e-g).The overestimation of the lower discharge values highlighted in Fig. 9g is also evident in the multiyear hydrograph (Fig. 8), which shows that the ED2+R simulation results overestimate (by about 40 % on average in the discharge values included in the range of 60 to 100 % in Fig. 9g) the observations during the dry seasons of the period under consideration.

Discussion
As the results in Table 3 and Figs.8-9 show, the one-way integration of ED2 with a routing scheme improves the performance of simulated daily discharges.Although this could appear obvious from a hydrological modeling perspective, the significance of this study lies in the fact that terrestrial biosphere models, which are widely applied to examine the impacts of climate and land use on the hydrology of the land surface, are typically "no-river representation" models.bon dioxide, and land-use changes simulated by terrestrial biosphere models with hydrological modeling improves the representation of the hydrological characteristics of basins characterized by large forest cover and/or large deforestation rates.In applications in the tropics, the one-way integration of the terrestrial biosphere model and the routing scheme (i.e., the two tools are not fully coupled) could lead to a partially inaccurate representation of the seasonally flooded ecosystems, a relevant aspect as documented in the literature (e.g., Cole et al., 2007).
As seen in Fig. 9, the performance of the model in simulating river flows in the basin is generally higher in the downstream sub-basins and poorer in the headwaters.Several factors are likely to cause this issue, both from the simulation of the hydrological dynamics in ED2, the flow partitioning (α and β parameters), and the basin hydraulic characteristics in ED2+R.The accurate calibration of the biosphere model with flux tower observations (Zhang et al., 2015;Longo, 2014) and the optimization of the flow partitioning make us believe that this variation in performance is due to the relatively coarse spatial resolution of the model in combination with the limitations typical of most land surface models in capturing the interactions with deep groundwater (Lobligeois et al., 2014;Zulkafli et al., 2013;Smith et Hydrol. Earth Syst. Sci., 21, 4629-4648, 2017 www.hydrol-earth-syst-sci.net/21/4629/2017/ al., 2004).We believe that the error is arising from the complexities associated with deep soils present in the headwaters of the Tapajós Basin.In particular, in the model application developed, soil layers are represented to a depth of 6 m (Table 1), which might be too shallow to realistically represent the conditions in the headwaters of the basin.The importance of groundwater is also evident from the calibration of the residence time parameter of the subsurface water flow: as shown in Fig. 7, in fact, especially in the headwaters, even small variations in the CB parameter greatly affect the model performance (specifically quantified with NSE in Fig. 7).The combined effect of groundwater interactions and spatial resolution is more evident in the upstream sub-basins because of the greater marginal contribution of baseflow in these areas.Surface flow accumulation, in fact, is lower in the headwaters.Therefore, in relative terms, the role of baseflow is more relevant in this portion of any basin.Further downstream, the effect of groundwater interactions and spatial resolution is, at least in part, masked by the larger rainfall-runoff contribution and the overall flow accumulation from the upstream sub-basins.Other recent hydrological simulations of the Tapajós have obtained higher accuracy (e.g., Mohor et al., 2015;Collischonn et al., 2008;Coe et al., 2008); however, these simulations were set up discretizing the basin into a finer spatial resolution grid (9 to 20 km vs. ∼ 55 km grid cells) and using hydrological tools able to reproduce highly detailed hydrodynamic characteristics of complex river systems (i.e., floodplain, lakes, wetlands, backwater effects) that are out of the scope of the tool presented in this study.The advantage of the ED2+R model is the ability to study the sensitivity of the river flows to global and regional changes as computed by traditional terrestrial biosphere models, but adding a more detailed hydrological feature with respect to a very simplistic-or no-river representation.The coarse spatial resolution of the global datasets used as input for ED2+R is, however, a limiting factor.Higher-resolution climatological data, vegetation, and land-use datasets, which would allow a finer resolution of the hydrological grid, are expected to improve the performance of the model by providing more detailed hydrological processes.On the other hand, a finer spatial resolution of the hydrological grid would also require a more detailed representation of the subsurface water in the model.In general, the tool can be used to study how different hydrological systems are being affected by changes in climate forcing and changes in ecosystem composition and structure arising from the combination of changing climate, rising atmospheric carbon dioxide, and land-use transformation.Additionally, ED2+R could potentially bridge one of the missing gaps for diagnosing and assessing feedback between atmosphere and biosphere with inland surface waters being represented as a dynamic system.

Conclusion
In this technical note, we present the integration of the terrestrial biosphere model Ecosystem Demography 2 (ED2) with the Muskingum-Cunge routing scheme.We tested the integrated model (ED2+R) in the Tapajós River basin, a large tributary of the Amazon in Brazil, for the period 1970-2008.
The results showed that the integration of a biosphere model with a routing scheme improves the ability of the land surface simulation to reproduce the hydrological and river flow dynamics at the basin scale.The main limitations highlighted in this case study were linked to the relatively coarse spatial resolution of the model and the rough representation of subsurface water flow typical of these kinds of models.Moreover, the terrestrial biosphere model ED2 and the routing scheme are presented here in a one-way integration.The full coupling of the routing scheme and ED2 could further improve the tool's ability to reproduce the water balance considering flooded ecosystems, a relevant feature in the simulation of environments like the tropical forest, where local evapotranspiration plays a primary role in the specific ecosystem's dynamics.In this first integration, our goal was to give the terrestrial biosphere model the ability to reproduce river flows through a routing scheme.With a fully coupled (i.e., twoway) integration, the model would be able to determine the grid cells that are likely to be saturated and use this information for the modeling of the ecosystem's dynamics.For instance, this could determine the increase of the mortality rate of plants that are sensitive to inundation.An additional limitation of the model could be identified in its inability to reproduce highly detailed hydrological dynamics of complex river systems (as, for instance, floodplain hydraulic features or backwater effects).However, such a detailed hydrological complexity was out of the scope of this study.Future efforts will address the highlighted limitations, while upcoming studies will use ED2+R to understand historical changes and future projections of the impacts of climate change and deforestation on the Amazon's water resources.

Figure 1 .
Figure 1.Schematic of the enthalpy fluxes (all arrows) and water fluxes (all but solid black arrows) that are solved in ED2.The schematic is based on Walko et al. (2000) and Medvigy et al. (2009).(Figure courtesy of Marcos Longo.)

Figure 5 .
Figure 5. Calibration of flow partitioning (α and β parameters in Fig. 2) between the ED2 and the ED2+R reservoirs.The color bar indicates the NSE values of the simulated vs. the observed river flow values (0: very different; 1: very similar).

Figure 7 .
Figure 7. Calibration of the residence times (τ ) of the flow within the ED2+R reservoirs of different grid cells in the domain through the adjustment of the non-dimensional C parameters.Overland, intermediate, and subsurface water flows are calibrated, respectively, through the adjustment parameters CS, CI, and CB (Fig. 2).The color bars refer to the model performance (NSE) of the specific parameter combination in the specific sub-basin.In red are the chosen combinations.At the bottom of each graph, there is a table providing the corresponding approximate average time of concentration (in days) for the cells in the sub-basin.(a) Upper Juruena (UJ); (b) Upper Teles Pires (UTP); (c) Lower Juruena (LJ); (d) Lower Teles Pires (LTP); (e) Upper Tapajós (UT); (f) Jamanxim (JA); and (g) Lower Tapajós (LT).

Figure 8 .
Figure 8. Calibration and validation of the river flow (m 3 s −1 ) at Itaituba (farthest downstream river gauge -Lower Tapajós sub-basin).ED2 output (green line), ED2+R (red line), and observations (blue dotted line).The dotted black line splits the calibration and validation periods.Similar comparison for each of the seven sub-basins is available in Appendix A.
Based on Poorter et al. (2006) Allometry for above-ground biomass Based on Eq. (2) of Baker et al. (2004) Allometry for leaf biomass Based on Cole and Ewel (2006) and Calvo-Alvarado et al. (2008) * Radiation-dependent parameters are given in the format xPAR, xNIR, and xTIR corresponding to values for photosynthetically active, near infrared, and thermal infrared, respectively.* * PFT-dependent parameters are given in the format xETR, xMTR, and xLTR corresponding