Comparison of runoff generation methods for land use impact assessment using the SWAT model in humid tropics

Hydrological responses due to deforestation in a humid tropical catchment were analyzed using two runoff genera‐ tion methods available in the Soil Water Assessment Tool (SWAT) model: the Curve Number (CN) and the GreenAmpt (GA) methods. The calibrated model, which per‐ formed well in simulating runoff under present land use condition in the Batanghari River Basin, Indonesia (42,960 km2), was then used to simulate runoff using past and future land use scenarios. Simulations showed similar changes in the annual water budget: decreasing evaporation and increasing total discharge. However, the two methods showed opposite changes in flow regimes: high flow increased (13%) under the CN while low flow increased (27%) under the GA. These results are associated with dif‐ ferences in runoff generation mechanisms, where surface runoff contributes to total discharge to a much larger extent under the CN (43%) than the GA (4%). Land use changes caused a reduction in infiltration rate, leading to higher high flow under the CN, while high flow did not change under the GA. Instead, lower evapotranspiration increased groundwater flow under the GA, and thus the steady low flow increased. This study suggests that the runoff genera‐ tion method should be selected carefully based on the dom‐ inant flow pathway of a catchment, particularly for land use impact studies in the humid tropics.


INTRODUCTION
Conversion of forests to other land uses due to economic activities is increasing more rapidly in developing countries, particularly in the humid tropics. These issues are often identified as the cause of the water disasters such as floods and droughts. The raised concerns surrounding the impact of rapid deforestation has urged policy makers to enact regional planning which allows economic growth as well as forest protection.
To date, most hydrological models were made based on experimental studies carried out in temperate regions. Most of the experimental studies in temperate regions show increases in water yield after forest conversion to pasture lands due to lower evapotranspiration of the new vegetation compared with forest trees. Further, most of these studies show reduced infiltration rates due to soil compaction by agricultural machinery. This suggest less infiltration which decreases low flow (Bruijnzeel, 2004). In the recent decades, a number of experiments have been conducted in the humid tropics showing different hydrological characteristics to temperate regions (Bonell et al., 1981;Nortcliff and Thornes, 1981;Noguchi et al., 1997a;1997b;Chappell et al., 1998;Elsenbeer and Vertessy, 2000). According to Bonell (2005), as the principal driver of runoff generation processes, "the short term, maximum rainfall intensities by storm event and daily rain totals are commonly an order of magnitude higher in the humid tropics than those which are experienced in humid temperate areas." Nortcliff and Thornes (1981) pointed out that such high rainfall rates are linked with high surface permeability of rainforests in the region. From a pedological perspective, Tomasella and Hodnett (2004) emphasized the 'unusual' properties of kaolinitic tropical soils compared to temperate clayey soils: low bulk density, high permeability and low available water capacity. Vieira and Santos (1987) showed that high rainfall and continuous leaching changes the chemical properties of kaolinitic tropical soils to form stable micro-aggregates from the size of silt to fine sand and thus the clayey soil in humid tropics can be very permeable. In addition, the soil in these regions is typically deeper due to fast weathering under constant high temperature and rainfall. Elsenbeer and Vertessy (2000) and Bonell (2005)  At a river basin scale, a typical approach to study the impact of deforestation is through the application of a hydrological model (Leemhuis et al., 2007;Olang and Fürst, 2011;Memarian et al., 2014;Marhaento et al., 2017). These studies include the following two main steps: 1) calibration and validation using present land use and 2) assessment of the impact of land use by changing parameters. Many of these studies do not discuss the representation of runoff generation before changing the parameters.
One of the widely used models to assess land use changes in tropical regions is the Soil Water Assessment Tool (SWAT) (Arnold et al., 1998;2012). Almost all the land use impact studies using the SWAT model in the tropics use the Curve Number (CN) as a runoff generation method, which is largely derived from empirical data sets in the US. Considering the difference in soil properties and dominant flow pathways in humid tropics, there is still limited knowledge on the applicability of a conceptual, empirically-based model such as CN, particularly towards land use change impact studies. In order to highlight this matter, this study uses one hydrological model (SWAT model) and compares an empirically-based (CN) and a process-based (Green-Ampt (GA)) runoff generation method to address the following questions: 1) What are the similarities and differences in hydrological responses to deforestation between the two methods?
2) Why do these runoff generation mechanisms generate similarities/differences in land use impact studies?

STUDY AREA
The Batanghari River Basin (42,960 km 2 ) is located in the western part of Indonesia in Jambi and West Sumatra Provinces ( Figure 1). Originating from the Bukit Barisan range at Kerinci and Solok counties, the river flows eastward and drains on the east coast of Sumatra into the South China Sea (Tan and Lim, 1998) With increasing economic activities, the annual total water demand in the river basin is projected to increase from 12 billion m 3 in 2008 to 13 billion m 3 in 2028 (Ministry of Public Works, 2012). As it runs through four protected national parks, which harbors some of the last viable popu-Muara Tembesi Station Figure 1. Map of the study area lations of endangered mammals and vegetation, the river also plays a crucial role for the ecosystem.
The topography of the basin varies from the mountainous area in the western part to flat floodplains in the eastern part. The mountainous area, the source of the Batanghari River tributaries, is part of Bukit Barisan range with its highest peak of 3,805 m a.s.l. at Mount Kerinci. The largest part of the basin consists of undulating hills covering 60% of the area with elevation between 10-100 m a.s.l. (Ministry of Public Works, 2003). The tributaries' confluence in the middle part of the basin close to Muara Tembesi station flows to the east through flat swamp land up to 200 km from the coast (NASA and METI (National Aeronautics and Space Administration and Ministry of Economy Trade and Industry of Japan), 2017). Parallel with the coastline, before the river delta, the river basin is characterized by peatlands of up to 6 m depth (Supiandi, 1988).
The climate in the basin varies following its geographical conditions (Ministry of Public Works, 2012), but the whole river basin can be categorized as a humid tropical basin based on Chang and Lau (1993). The temperature in the basin ranges from 22.2 ± 0.2°C (upstream) to 26.8 ± 0.2°C (downstream). Data provided by Indonesian Meteorology, Climatology and Geophysical Agency (BMKG) at Depati Parbo and Sultan Thaha stations covers the period 2001 to 2013 (BMKG, 2019). The basin average annual rainfall is 2021 ± 247 mm (BMKG and Ministry of Public Works' data of 23 stations, 2019) which is characterized into two rainy seasons (Eguchi, 1983;Hamada et al., 2002) with monthly average above 100 mm . The rainfall of most areas in the basin, based on Aldrian and Susanto (2003), is strongly influenced by two monsoons namely the wet northwest monsoon from November to March and the dry southeast monsoon from May to September while a small part (the northeast and northwest tips of the basin) is influenced by the southward and northward movement of the inter-tropical convergence zone. The FAO (Food and Agriculture Organization)/UNESCO (The United Nations Educational, Scientific and Cultural Organization) soil map shows that the upper and middle river basin areas are dominated by Acrisols and Ferrasols and the downstream part of the basin, particularly along the river banks, is dominated by Fluvisols. Prior studies on soil properties in the upstream area carried out by Araki (2019) at a hillslope scale and Susiwidiyaliza (2015) at a subcatchment scale show thick soil layer ranging from 1,100 mm to 4,500 mm depth with high hydraulic conductivity up to 320 mm/h. Supiandi (1988) carried out soil sampling in the downstream area (east of Jambi to the coastline) and showed soil depths ranging from about 1.4 m to more than 3 m.

METHODS
This study uses the SWAT model to simulate runoff in the Batanghari River Basin. We compare two runoff generation methods, CN and GA. These methods separate rainfall to surface runoff and soil water. Water entering the soil layers is redistributed throughout the soil profile using a storage routing technique until the whole soil profile has E.M.S. YAMAMOTO ET AL.
uniform water content. The redistribution is computed simultaneously with soil evaporation, plant uptake and transpiration, lateral flow and percolation to shallow and deep aquifers. The lateral subsurface flow is calculated and subtracted from each soil layer using a kinematic storage model (Sloan and Moore, 1984). The potential evapotranspiration was calculated collectively using Penman-Monteith equations (Monteith, 1965), and the evapotranspiration from plants and soil layers were calculated separately as described by Ritchie (1972). Percolation from the bottom soil layer enters the vadose zone before it is partitioned between shallow and deep groundwater. Water in shallow groundwater contributes to river discharge as return flow or is taken by deep rooting plants. All runoff components are summed and routed through the channel network using Muskingum methods. Further detailed description on the SWAT model and its runoff generation methods can be viewed in the supplement (Text S1).

Model design and setup
This study divided the Batanghari River Basin into 359 sub-basins/HRUs (Hydrologic Response Units) in the SWAT model using DEM (Digital Elevation Model), land use and soil maps as explained in the data description section. Furthermore, we modified the crop and management database to fit the plantation growth in tropical regions based on Kilonzo (2013) and Alemayehu et al. (2017). The values were adjusted based on LAI (Leaf Area Index) from ECOCLIMAP (Eco-climatic Map) (Masson et al., 2003;Marhaento et al., 2017) (Table SII).
Two SWAT models were set up with GA and CN based on present land use conditions in 2015. A four-year simulation (2012-2015) was carried out with a two-year warm up period (2010)(2011). This study uses a three step calibration, which is further explained in the supplement (Text S2). After the calibration (2012-2013) and validation process (2014)(2015), the optimized parameter sets for each model were applied using past (1990, 1997, and 2005) and future (2030, 2035, and 2040) land use maps with the same rainfall period (2010)(2011)(2012)(2013)(2014)(2015).

DATA DESCRIPTION
Due to unavailability of sub-daily data from rain gauges in the basin, satellite based gridded hourly rainfall was used. This study uses GSMaP (Global Satellite Mapping of Precipitation) Reanalysis version 6 (January 2010-February 2014) and GSMaP MVK (Motion Vector Kalman filter) version 6 (March 2014-December 2015) with a 0.1-degree resolution. The GSMaP Reanalysis v.6 was analyzed and used by Yamamoto et al. (2019) in the same river basin. In this study, GSMaP MVK v.6 was also compared with available rain gauge data. The comparison demonstrated a similar spatial pattern with the gauged data. The monthly average of GSMaP MVK followed similar seasonal pattern as the gauged data as well as the GSMaP Reanalysis (see Supplement Text S3 for comparison of GSMaP datasets and available rain gauge data). The daily climatic data wind speed, temperature, surface pressure, specific humidity, downward long-and shortwave radiation fluxes were obtained from WFDEI (WATCH-Forcing-Data-ERA-Interim) (Weedon et al., 2014) with 0.5-degree resolution. For model calibration and validation, water level data from Muara Tembesi Station was converted to daily discharge using rating curves provided by the Ministry of Public Works of Indonesia.
The simulation used the DEM data from HydroSHEDS (Lehner et al., 2008) with 30-arc second resolution. It also used the Digital Soil Map of the World (DSMW) and derived soil properties (FAO/UNESCO, 2003). In the STEP 3 calibration, this study used the DSMW soil distribution map and other soil properties and optimized the soil depth and hydraulic conductivity based on field observations (see Supplement Text S2). Past and present land use maps were classified from four Landsat images from 1990,1997,2005 and 2015 with spatial resolutions of 30 m . Future land use maps were classified using the dynamic model CLUE (Conversion of Land Use and its Effect) (Utami, 2017) (Figure S1 and Table SI). Figure 2 shows the simulated discharge resulting from the calibration at the Muara Tembesi station. These hydrographs show that, with calibration of parameters, the CN and GA methods can results in similar patterns. Based on Nash-Sutcliffe Efficiency (NSE), both methods showed satisfactory performance (GA: 0.59 and CN: 0.64) in the The calibrated models were then used to simulate past and future land uses to assess the hydrological response to the land use impact. Figure 3a shows the trend in annual water budgets using land use maps from 1990 to 2040. The simulations from both methods gave similar trends in the water budget, i.e. evapotranspiration (ET) decreasing and discharge (Q) increasing. The ET decreased by 15 mm (1.5%) using the GA method and 23 mm (2.7%) using the CN method. The Q increased by 47 mm (3.9%) using the GA method and 71 mm (5.5%) using the CN method. The decrease in evapotranspiration and increase in water yield as an impact of forest conversion to agricultural area is in agreement with the general consensus paired catchment studies in humid tropics (Kuraji and Paul, 1994;Malmer, 1996;Costa and Foley, 1997;Kleinhans and Gerold, 2004;Panday et al., 2015). In addition, the results are also in agreement with the discharge observed data (Figure 4) where a Mann-Kendall test shows that there is no trend within a significant level of 95% with a Sen's slope of 0.017.

RESULTS
Although the annual water budgets of both methods showed similar changes, the flow duration curves (FDCs) showed opposite changes ( Figure 5). The CN method showed that high flow (Q 1 ) increased by 13% or 472 mm while the low flow (Q 99 ) decreased by 4% or 28 mm with the conversion of forest to agriculture. Q 1 and Q 99 are discharges at 1 and 99 percentiles. On the other hand, the GA method showed Q 1 shows no significant changes (3% or 111 mm) and Q 99 increased by 27% or 123 mm.
The different changes in FDCs between the two methods is likely associated with the difference in the dominant runoff components (Figure 3b). The CN method has a much higher fraction of surface runoff in discharge compared to the GA method. Before deforestation, the fraction is 35% in CN while in GA it is 2%. After deforestation, it increases to 45% in CN and to 5% in GA. These differ-  1989 1991 1993 1995 1997 1999 2001 2003 2005 2007 2009   ences were pointed out in earlier studies (Kannan et al., 2007;Bauwe et al., 2016;Tasdighi et al., 2018). For the CN method, where we use the default initial CN 2 values and assume that the soil is dry, the surface runoff starts after rainfall exceeds 21.7 mm/day in a forest sub-basin. In the same sub-basin, the surface runoff starts only when the rainfall exceeds 65.6 mm/h (effective hydraulic conductivity (K e )) using the GA method. The frequency of rainfall exceeding 21.7 mm/day is 0.6% (including dry days) in our study area between 2012 and 2015, while the frequency exceeding 65.6 mm/h is only 0.002%. Furthermore, the higher CN value under agriculture allows surface runoff to start with lower rainfall amount (21.7 mm/day); thus, the fraction of surface runoff becomes much higher with deforestation. In contrast with the CN method, the largest fraction of discharge in the GA method is subsurface flow (Figure 3b). The calculation of infiltration in the GA method is directly affected by the infiltration rate based on observation data (see Table SIII) in the Batanghari River Basin, which is not the case with the CN method. We found that for the GA method, the fraction of subsurface runoff components changed after the land use change. The groundwater flow increased by 23% and lateral flow from the soil layer decreased by 18%. The opposing trends between groundwater and lateral flow are due to lower evapotranspiration from shallow aquifers in agricultural fields compared to forests. The increase in groundwater flow in the GA method causes an increase in mean annual low flow due to lower evapotranspiration from shallow aquifers. Bruijnzeel (2004) proposed two 'infiltration trade-off' hypotheses for land use change impact in hydrology: 1) the surficial and soil properties remains resilient following for-est conversion so as to keep the dominant pathways virtually unchanged and vertical percolation to groundwater remains uncompromised, 2) the surficial and soil structure are collapsed following forest conversion, corresponding to a significant reduction in vertical percolation flux in the uppermost soil layers.

DISCUSSION
The humid tropics seem to have a wider spectrum of dominant flow pathways (Bonell, 1998). However, the occurrence of infiltration exceeding overland flow are generally restricted by occasionally high magnitude rainfall (Wierda et al., 1989;Malmer, 1993). Our recent field observations  show that even with one order of magnitude lower infiltration rate after forest conversion, surface runoff remains negligible -suggesting that the surficial and soil properties remain resilient and maintain subsurface flow as the dominant pathway. Observations using isotope tracers in a tropical montane cloud forest region by Muñoz-Villers and McDonnell (2013) also showed that even with one or two order of magnitude lower surface hydraulic conductivity in pasture, storm runoff is still dominated by groundwater sources. Bonell (1998) summarized evidence from earlier works in the humid tropics that showed that delayed flow has a more significant role than the stormflow component in increasing water yield due to forest conversion. Kleinhans and Gerold (2004) presented similar hydrographs as a result of field observations in Central Sulawesi, Indonesia where low flow in the natural forest was lower than in mixed forest. Based on these studies, the land use changes in humid tropics, especially in the Batanghari River Basin, seem to follow the first hypothesis of Bruijnzeel, which is in line with simulation results using the GA model.
Our results show that both CN and GA can be calibrated to fit the observed hydrograph and give similarly satisfactory NSE results, even though the dominant runoff components are different. Beven (1991) stated that it is not so difficult to shape a catchment hydrograph by modifying

RUNOFF GENERATION IN HUMID TROPICS
parameters of a model in simulating the 'loss' and 'delay' processes. However, the results also show that when these calibrated models are applied to different land use maps for land use impact assessment, the hydrographs change in different ways ( Figure 6). The way the hydrograph changes with land use change depends on the calculation of runoff generation for each type of land use. A process-based model such as GA attempts to represent physical properties of each land use in the calculation of runoff generation. If the model parameters from soil properties are set properly, runoff generation could be represented more reasonably.
On the other hand, the CN model uses an empirical parameter to directly calculate its runoff components based on a combination of land cover, soil group and land management; the original parameters for which have been derived mostly in temperate regions. Particularly in the humid tropics where conditions are different, the use of these parameters in land use change impact study requires further confirmation to provide an accurate representation of flow pathway in each land use case.

CONCLUSIONS
Changes to flow regimes due to deforestation is one of the main concerns in land use impact studies. It is particularly useful knowledge in policy making for regional planning. This study highlights that different runoff generation methods with the same SWAT model can give opposing change in flow regimes, even though they are calibrated to have similar hydrographs in the calibration period. The differences in flow regime changes are related to the dominant flow pathways resulting from different runoff generation mechanisms. The CN method shows a higher proportion of surface flow, while in the GA method the subsurface flow is higher. A field study conducted by Araki et al. (2019) in the Batanghari River Basin indicated that the dominant flow pathway is subsurface flow. Hence, in the case of our study site, the use of empirical models may give misleading results if they are not calibrated to represent the correct dominant flow pathway. Compared to CN, the calibrated GA method agrees with our perceptual model in this particular basin.
Overall, this study suggests that model selection, particularly for land use impact studies, should not be solely based on goodness of fit to observed hydrographs but should also accurately represent the dominant flow pathway in the basin.

SUPPLEMENTS
Text S1. Detailed description of hydrologic model and runoff generation methods (CN and GA) Text S2. Three steps calibration and results Text S3. Comparison of GSMaP datasets and available rain gauge data Figure S1. Present and future land use classification figures Figure S2. Spatial and temporal distribution of GSMaP Reanalysis v.6 and gauged data Figure S3. Spatial distribution of GSMaP MVK v.6 and gauged data