Environmental Modelling & Software 83 (2016) 310e325

Contents lists available at ScienceDirect

Environmental Modelling & Software

journal homepage: www.elsevier.com/locate/envsoft

Simulating the spatio-temporal dynamics of soil erosion, deposition,

and yield using a coupled sediment dynamics and 3D distributed
hydrologic model
Tan Zi a, Mukesh Kumar a, b, Gerard Kiely c, Ciaran Lewis c, John Albertson a, d, *
a
Department of Civil and Environmental Engineering, Duke University, NC, USA
b
Nicholas School of the Environment, Duke University, NC, USA
c
Centre for Hydrology, Micrometeorology and Climate Change, Dept. of Civil and Environmental Engineering, University College Cork, Cork, Ireland
d
School of Civil and Environmental Engineering, Cornell University, NY, USA

a r t i c l e i n f o a b s t r a c t

Article history: Since soil erosion is driven by overland ow, it is fair to expect heterogeneity in erosion and deposition in
Received 3 August 2015 both space and time. In this study, we develop and evaluate an open-source, spatially-explicit, sediment
Received in revised form erosion, deposition and transport module for the distributed hydrological model, GEOtop. The model was
11 May 2016
applied in Dripsey catchment in Ireland, where it captured the total discharge volume and suspended
Accepted 3 June 2016
sediment yield (SSY) with a relative bias of
1.2% and
22.4%, respectively. Simulation results suggest
that daily SSY per unit rainfall amount was larger when the top soil was near saturation. Simulated
erosion and deposition areas, which varied markedly between events, were also found to be directly
Keywords:
Distributed hydrologic model
inuenced by spatial patterns of soil saturation. The distinct inuence of soil saturation on erosion,
Soil erosion deposition and SSY underscores the role of coupled surface-subsurface hydrologic interactions and a
Sediment spatio-temporal dynamics need to represent them in models for capturing ne resolution sediment dynamics.
Coupled modeling 2016 Elsevier Ltd. All rights reserved.
Sediment transport model
GEOtop model

Software availability 1. Introduction

Software Name: GEOtopSed Soil erosion by rainfall and overland ow is a widespread threat
Developers: Tan Zi to soil fertility and water quality. Accurate estimation of soil loss
Contact Address: Department of Civil and Environmental and its spatial distribution is often needed for pollutant risk ana-
Engineering, Duke University, Durham, North Carolina, lyses, reservoir management, agriculture productivity forecasts,
27708, US and soil and water conservation. In this regard, several distributed
Email: [email protected] models have been developed to obtain erosion estimates (DeRoo
Year First Available: 2015 et al., 1996; Wicks and Bathurst, 1996; Morgan et al., 1998;
Hardware Required: Desktop/Laptop with 2 GHz CPU, 2 GB RAM or Hessell, 2005; Jain et al., 2005; de Vente et al., 2008). Notably,
more majority of distributed erosion-deposition models e.g., WEPP,
Operating System Required: Macintosh OSX 10.4 or newer; EUROSEM etc., consider simplistic representations of vertical and
Windows XP or newer; Linux lateral subsurface water ow, and often do not account for the
Libraries Required: ASCII, FLUIDTURTLES, GEOMORPHOLOGYLIB, lateral subsurface water movement, or the coupled dynamic
KeyPalette, MATH interactions between vadose zone and the groundwater table, or
Cost: Free the evolution of soil moisture and groundwater with evapotrans-
Source Code: https://sourceforge.net/projects/geotopero/ piration. Given that the detachment, transport, and deposition of
Program Language: C soil are dominantly inuenced by the velocity and volume of
overland ow (Julien and Simons, 1985), which in turn may be
* Corresponding author. Department of Civil and Environmental Engineering, inuenced by antecedent soil moisture conditions (Legates et al.,
Duke University, NC, USA. 2011; Penna et al., 2011; Jost et al., 2012; Chen et al., 2014;
E-mail address: [email protected] (J. Albertson).

http://dx.doi.org/10.1016/j.envsoft.2016.06.004
1364-8152/ 2016 Elsevier Ltd. All rights reserved.
 T. Zi et al. / Environmental Modelling & Software 83 (2016) 310e325 311

Hueso-Gonz
alez et al., 2015), subsurface heterogeneity (Lewis et al., monotonically with precipitation amount and energy? If not, does
2012; Ghimire et al., 2013; Orchard et al., 2013; Zimmermann et al., the hydrologic response of the watershed has a role to play in the
2013; Niu et al., 2014; Tao and Barros, 2014), and groundwater departure from monotonic relation? c) Does the simulated source/
distribution (Kumar et al., 2009; Miguez-Macho and Fan, 2012; sink area of sediments vary spatially from one event to other? If yes,
Rosenberg et al., 2013; Safeeq et al., 2014; von Freyberg et al., 2015), is the variation driven by hydrologic state, specically the surface
it is important to consider the coupled impacts of antecedent hy- soil saturation state? and d) To what extent does the linear relation
drologic states (soil moisture and groundwater distribution) and between erosion and the slope-length factor (product of specic
subsurface hydrogeologic properties on sediment generation and catchment area and slope), which is often used in USLE-based
yield. Failing to do so may limit the applicability of these models to model representations (e.g. USLE (Wischmeier and Smith, 1978),
a few events (Hessel et al., 2006; Mati et al., 2006; Ramsankaran RUSLE (Renard et al., 1991), RUSLE2 (Foster et al., 2005)), hold for
et al., 2013) or to regimes where the dynamic role of antecedent GEOtopSed simulated states and uxes?
conditions and subsurface heterogeneity on erosion are not large
enough. Heppner et al. (2006) made signicant headway in this 2. Process formulation, model implementation, and
direction by coupling sediment processes within an integrated verication
hydrologic model, InHM (VanderKwaak and Loague, 2001). The
study specically evaluated the rainfall splash erosion component 2.1. The GEOtop model: a short review
of the model on a 6 m by 2.4 m plot. Heppner et al. (2007) used the
same model to perform sediment-transport simulations for six The open-source GEOtop model (Rigon et al., 2006) is process
events in a 0.1 km2 rangeland catchment. It is to be noted that InHM based and simulates core hydrological processes such as unsatu-
solves subsurface ow using the variably saturated 3D-Richards rated ow, saturated ow, overland ow, stream ow generation/
equation, while surface ow is simulated using diffusion wave routing, and surface energy balances. Overland ow modeling is
approximation of St. Venant equation. Equations corresponding to performed using the kinematic wave approximation of St. Venant
these coupled processes are spatially discretized using a control equation while subsurface ow and soil moisture simulations are
volume nite element strategy on each unstructured grid. A global performed by solving a variably-saturated representation of 3D
implicit solver is used to perform the simulation. Another notable Richards equation. By solving the Richards equation, GEOtop model
effort in this direction was by Kim et al. (2013), who coupled can simulate the surface runoff generation processes due to both
sediment processes within a hydrologic and hydrodynamic model inltration excess and saturation excess, and can also redistribute
tRIBS-OFM and validated their model against analytical solutions. the sub-surface water both laterally and vertically, as determined
Similar to InHM, tRIBS-OFM is also an unstructured grid based by the head gradient. The model has been extensively tested and
model. The model uses a gravity-dominated formulation (Cabral validated in Bertoldi (2004). The water and energy balance calcu-
et al., 1992) to simulate vadose zone ow and a quasi-3D Boussi- lations in GEOtop were recently rened to account for soil freezing
nesq's equation under the Dupuit-Forchheimer assumptions to and thawing effects (Endrizzi et al., 2014). In summary, with
simulate groundwater ow (Ivanov et al., 2004). The model was detailed water and energy balance modules, GEOtop can provide
used to evaluate sediment yield simulations for 10 events in a accurate simulations of evapotranspiration and soil moisture dy-
0.036 km2 Lucky Hills watershed located in southeastern Arizona, namics (Bertoldi et al., 2014; Della Chiesa et al., 2014), given
USA. Development of these physically-based integrated models of adequate watershed data. By simulating coupled hydrologic states
hydrology and sediment dynamics has opened new opportunities, (e.g. surface ow depth, soil moisture and groundwater) on each
especially in regards to understanding the impact of the hydrologic grid of the model domain, the model is well suited to study the
state on spatio-temporal distribution of erosion, deposition and inuence of watershed properties and subsurface states on
yield. Notably, the aforementioned two models are not open- spatially-distributed runoff, an important control on erosion, at
source. multiple scales. Furthermore, as an open source software (http://
Here, we develop an open-source, spatially-explicit, structured- www.geotop.org/wordpress/), the GEOtop model provides a com-
grid based, sediment erosion/deposition module for a 3D surface- plete hydrological model framework with ease for extensions. One
subsurface hydrologic model, GEOtop (Rigon et al., 2006; Endrizzi such example is the incorporation of landslide occurrence predic-
et al., 2014), and evaluate its applicability in explaining the sedi- tion within the GEOtop framework by Simoni et al. (2008).
ment yield dynamics. Similar to InHM (Heppner et al., 2007), the
GEOtop model also solves subsurface ow using the variably
2.2. Process formulation of the sediment dynamics model
saturated 3D-Richards equation, while surface ow is simulated
using kinematic wave approximation of St. Venant equation. The
The sediment dynamics model developed here takes advantage
sediment dynamics model developed here takes advantage of the
of the GEOtop simulated distributed hydrological states such as
GEOtop simulated distributed hydrological states such as moisture
moisture content, surface ow depth, and ow velocity. Here we
content, surface ow depth, and ow velocity. The model accounts
only highlight the aspects of the model that are most relevant to the
for the inuence of spatial heterogeneities in land surface charac-
sediment erosion, deposition and transport modeling. Readers may
teristics, subsurface hydrogeology, and antecedent conditions in
refer to GEOtop model papers (Rigon et al., 2006; Endrizzi et al.,
the generation of overland ow, and hence on the erosion and
2014) to learn more about the individual process representations.
deposition of sediment in the catchment. The model developed
GEOtop simulates soil moisture in each subsurface layer by
here was applied on a much larger catchment (area 15 km2) and
solving the 3D Richards equation:
for a longer period (simulation duration 2 years) than in Heppner
et al. (2007) and Kim et al. (2013), allowing validation of the vH
coupled model for extended wet and dry periods. The coupled CHf Sw Ss V$
KVH Sw 0 (1)
vt
model is then used synergistically with the observed data to
answer four pointed questions: a) Is the performance of the GEO- where K [m s
1] is the hydraulic conductivity, H [m] is the sum of
topSed model for simulating SSY, dependent on the ow regime pressure and potential head, and Sw is the source/sink mass ux
and the model's ability to capture streamow response? b) Does [s
1], Ss is the specic storage coefcient [m
1], f is porosity [
],
the daily suspended sediment yield (SSY) from the watershed vary and C(H) is the specic moisture capacity function.
312 T. Zi et al. / Environmental Modelling & Software 83 (2016) 310e325

Surface overland ow is routed along the ow direction using calculated using:

kinematic wave approximation of the Saint-Venant equation rep-
resented as: vhc C c vqc C c
In (12)
vt vx
vh vq
qL (2) where qc is the discharge per unit width in channel [m s
1], hc is the
vt vx
water depth of channel [m], Cc is sediment concentration in chan-
where h is depth of overland ow [m], x is a local coordinate system nel [kg m
3], and In is the sediment exchange rate per unit area
oriented along the ow direction [m]; q is ow rate per unit width between adjacent land cells and the channel cell [kg m
2 s
1]. The
[m2 s
1], and qL is vertical inow or outow rate (to the surface suspended sediment yield (SSY) [kg] at the outlet of the catchment
water) per unit area [m s
1]. The ow direction is dened along the was the integration of suspended sediments over the period of
line of steepest downslope head gradient between the grid and its interest:
eight neighbors. The overland ow equation is coupled to the Z
continuity equation of sediment transport at each time step using: SSY Qol Col dt (13)
vhC vqC
Ex (3) where Qol is the discharge rate at outlet [m3 s
1], Col is the sediment
vt vx
concentration at the outlet cell [kg m
3] and dt is an hourly inte-
where C is sediment concentration in the overland ow within each gration time step. The period of integration in the ensuing analyses
cell [kg m
3], and Ex is the exchange rate of sediment per unit varied from hourly to annual scale and has been appropriately
surface area [kg m
2 s
1] at the interface of soil and water. Ex is identied at relevant locations.
composed of three major mechanisms: rainfall splash detachment
(DR), ow detachment (DF) and deposition (DP), as: 2.3. Model implementation
Ex DR DF
DP : (4)
The GEOtop model rst solves the nite difference discretization

2
1
DR [kg m s ] is approximated by DeRoo et al. (1996): of 3D Richards equation (Equation (1)) using Bi-conjugate gradient

 stabilized method (Vorst, 1992). An absolute numerical tolerance of
Ke
1:48h 10
4 mm is used for ensure water balance. This requires performing
DR 0:1033 e 3:58 *I (5)
z simulations at an adaptive time-step with time intervals generally
ranging from minutes to hours. After solving the 3D Richards
where z is soil cohesion [kPa], Ke is rainfall kinetic energy equation, H in the top layer is used to calculate the head gradient for
[J m
2 mm
1], and I is the precipitation intensity [mm h
1]. DF and estimation of inltration/exltration ux in Equations (1) and (2),
DP [kg m
2 s
1] are related to transport capacity (TC) based on the at a user dened interval (hourly in this study). Finite difference
erosion-deposition theory proposed by Smith et al. (1995). discretization of the surface ow, channel routing, and sediment
transport equations (Equations (2), (3) and (12)) are then solved
DF TC
C*y*vs (6) altogether using a forward explicit Euler method. To ensure sta-
bility, Courant-Friedrichs-Lewy (CFL) condition is used. The Cou-
DP TC
C*vs (7) rant number is set to 0.25. If the CFL condition is not satised, the
time step is reduced by 75% adaptively until it gets satised.
where y is an efciency coefcient that is a function of soil cohesion Although the integrated model simulations are performed using an
(DeRoo et al., 1996), and vs is settling velocity of the particles adaptive time-stepping scheme, the model outputs are printed at
[m s
1]. The dependence of soil cohesion on soil moisture (Bullock constant intervals.
et al., 1988) and root tensile strength (Wu et al., 1979) is captured
using:
2.4. Model verication at plot scale

2
q
zs zss (8) GEOtopSed was rst evaluated on a soil ume setup discussed in
qs Ran et al. (2012). The setup included a rainfall simulator, a tilted soil
ume, an overland ow collector, and a set of soil water content
zadd 1:2*RS (9) monitor devices. The rainfall simulator was used to generate rain
events of varying intensity, duration, moving direction, rainfall
z zadd zs (10) position, and no-rainfall interval. The size of the soil ume was 5 m
long by 1 m wide. In the ume, a 30 cm thick soil layer with a xed
where zs, zss, zadd are bare soil cohesion, saturated bare soil cohesion slope of 25 , and containing 13% clay, 58% silt, and 29% sand was set
and additional cohesion by roots respectively, q and qs are the soil on top of a 5 cm thick sand layer. Surface runoff and sediment were
moisture and saturated soil moisture contents respectively, and collected at the end of the soil ume in a metal container placed
RS is the root total tensile strength [kPa]. The transport capacity is slightly below the soil surface. All other boundaries of the soil ume
based on the experiments conducted by Govers (1990): were impermeable. Readers are referred to Ran et al. (2012) to learn
more about the experimental setup. The model simulations were
TC au
ucr b *r (11) performed on a grid discretization with a resolution of
1 m  1 m  10 cm in x, y and z directions. The model results were
where u is the unit stream power [m s
1] (Yang, 1972), ucr is the output at 3 min interval in this experiment.
critical power that initiates ow detachment of soil particles Out of 33 rainfall scenario experiments discussed in Ran et al.
[m s
1], r is the density of soil particles [kg m
3], and a and b are (2012), here we summarize the results for three scenarios (No.10,
empirical parameters related to soil particle size. No.14, and No.18) with markedly different rainfall intensities
The sediment mass balance in each stream channel reach was (Table 1). For scenario 10 (Fig. 1) which included rain events of
 T. Zi et al. / Environmental Modelling & Software 83 (2016) 310e325 313

highest rainfall intensity (>80 mm/h) and shortest duration among antecedent soil moisture have been found to play an important role
the three selected scenarios, the model was able to capture the on runoff generation (Lewis et al., 2013) and nutrient transport
observed timing of both ow and SSC peaks. Observed ow peaks (Warner et al., 2009) in the catchment, the site serves as a good test
were also captured accurately. Further analyses of simulation re- case for validating the integrated model at catchment scale. The
sults suggest that a smaller runoff peak for the rst event, even Dripsey catchment is located approximately 25 km northeast of
though the duration and intensity of all the precipitation events Cork, and has an area of 15 km2. The elevation of this catchment
were almost identical, was because of drier antecedent soil mois- ranges from 60 to 210 m. It is a beef and dairy producing agricul-
ture conditions in the top soil layer. Once the moisture decit of top tural catchment and is almost 100% covered by perennial ryegrass.
layer is fullled, additional precipitation contributes to inltration- The catchment slopes gently, with around 85% of the area having
excess runoff even while the lower soil layers are not yet saturated. less than 3% grade. Gleys and podzols are the two major soil types.
Since the second event occurred only 0.25 h after the rst one, the The experimental eld site is managed by the Hydromet team of
top soil layer was still near saturation resulting in a larger runoff the University College Cork (UCC). There is a meteorological ux
peak than the rst event. For later events, several model cells were tower at the top of the catchment (elevation 192 m) where radia-
fully saturated i.e. all vertical layers were saturated, indicating that tion, wind speed, air temperature, surface temperature, relative
both inltration-excess and saturation-excess processes played a humidity, and soil moisture (at ve depths) have been measured at
role in generation of runoff peaks for these events. The model was 30 min interval since 1998 (Albertson and Kiely, 2001). All the
also able to capture the relatively smaller runoff peak magnitude in meteorological data from the site are archived in the FLUXNET re-
scenarios 14 and 18 due to the reduced intensity of events. Overall, pository (http://uxnet.ornl.gov). At the catchment outlet (eleva-
the timing of simulated SSC peaks also matched the peak in tion 60 m), stream ow was monitored continuously at 30 min
observed records. For scenarios 10, the model was able to capture interval for a period of over two years (2002e2003). Flow-
the magnitude and the decreasing trend in SSC for events 2 to 5. weighted water samples at the catchment outlet were collected
The rst simulated SSC peak was however underestimated. Ana- using an ISCO 6712 auto-sampler with intake set at approximately
lyses of simulated states suggest that the decreasing trend in SSC 0.25 m above the streambed (Lewis, 2003, 2011).
may be explained based on the combined effects of increase in ow The climate in the study region is temperate maritime, and is
volume which reduces SSC, reduction in sediment generation with characterized by high humidity and a lack of temperature extremes
small decrease in peak ow, and an increase in soil cohesion with during the year. The minimum daily temperature for years 2002
increasing soil moisture which impedes sediment generation. For and 2003 was
0.2 C. The mean annual precipitation locally is
scenario 14, the model was again able to capture the decreasing approximately 1400 mm. The annual precipitation in the Dripsey
trend in SSC for the four events. In scenario 18, the SSC was catchment for year 2002 and 2003 was 1823 mm and 1178 mm,
measured only for two events as the rst precipitation pulse did not respectively. Winter and spring were the wetter seasons while
generate any response. The model captured the magnitude of summer was generally dry (Fig. 3). No precipitation was recorded
runoff peak corresponding to the third event, but underestimated during the period with temperature below 0  C during the two
for the second event. It is to be noted that the runoff generated for years. Suspended sediment losses over the catchment were esti-
the second event is negligible and SSC estimates are very sensitive mated from measured data of stream ow volume and suspended
at these magnitudes. sediment concentration. The monthly variations in suspended
In summary, the GEOtopSed model generally captured both the sediment yield are shown in Fig. 3. Notably, the runoff ratio in
trend and the quantity of runoff and suspended sediment con- January 2003 is larger than one. This indicates that antecedent
centration (SSC). It is to be noted that an accurate simulation of groundwater recharge participated in delayed streamow response
runoff and consequently of SSC was possible because of compre- in this month. Marked variations in runoff ratio through the year
hensive representation of surface and subsurface hydrological underscores the need for appropriate partitioning of the water
processes in the model. budget across different stores of the hydrologic continuum.

3. Model application at watershed scale 3.2. Input data

3.1. Site description GEOtop requires a digital elevation model (DEM), land use/land
cover (LULC) map, and soil type map to simulate the hydrological
The coupled model was applied at a small experimental catch- processes for a catchment. The soil type map and soil parameters
ment, in Dripsey, Ireland (Fig. 2). Given that groundwater and were obtained from Irish Forestry soils (IFS) database and in situ soil

Table 1
Summary of plot experiment. The experiment data are from Ran et al. (2012). Numbers in the parenthesis are simulation results from GEOtopSed.

Scenario no. Event Rainfall intensity Rainfall amount Rainfall duration (h) No-rainfall Time of runoff Runoff peak SSC (kg/m3)
(mm/h) (mm) interval (h) peak (h) (ml/s)

10 US001 83.99 21.00 0.25 0 0.12 (0.15) 30.2 (31) 57.08 (25.6)
US002 84.42 21.10 0.25 0.25 0.09 (0.15) 38.8 (41) 30.05 (32.6)
US003 85.05 21.26 0.25 0.5 0.20 (0.25) 40 (40) 19.53 (32.7)
US004 86.08 21.52 0.25 1 0.12 (0.15) 38.9 (35) 14.30 (24.2)
US005 83.80 20.95 0.25 3 0.18 (0.15) 40.8 (36) 10.63 (20.3)

14 UM001 59.15 29.57 0.5 0 0.24 (0.4) 13.6 (17) 69.35 (89.4)
UM002 59.28 29.64 0.5 0.5 0.27 (0.3) 21.95 (21) 36.76 (33.4)
UM003 59.47 29.73 0.5 1 0.21 (0.3) 22.4 (23) 24.20 (18.3)
UM004 61.09 30.54 0.5 3 0.27 (0.5) 21.45 (20) 16.78 (18.7)

18 UW001 29.63 29.63 1 0 0.00 (0) 0 (0) 0.00 (0)

UW002 23.45 23.45 1 1 0.27 (0.65) 0.233 (1) 4.77 (0.001)
UW003 34.35 34.35 1 3 0.33 (0.65) 7.5 (11) 6.71 (6.38)
314 T. Zi et al. / Environmental Modelling & Software 83 (2016) 310e325

Fig. 1. Modeled and observed records for rainfall scenarios 10 (left column), 14. (middle column), and 18 (right column). The rst row are the time series of precipitation and surface
runoff; the second row are time series of precipitation and suspended sediment concentration (SSC); and the third row are time series of plot-average soil moisture of three
subsurface layers. In the top two rows, precipitation (P) is plotted on the secondary axis.

Fig. 2. Location of Dripsey catchment in County Cork, Ireland. The location of the meteorological tower is where Dripsey is marked in yellow. The right top panel is the Google
Earth image of Dripsey catchment. The lower panel shows the conceptualization of rills (discussed in Section 3.3).(For interpretation of the references to colour in this gure
legend,the reader is referred to the web version of this article.).

samples (Lewis, 2011). LULC parameters were derived based on resolution. Meteorological data such as precipitation, temperature,
classications using the Corine land cover 2000 database and land incoming shortwave radiation, air pressure, relative humidity, wind
use data observed in the catchment. The stream channel was speed and direction in half-hourly time steps were collected at the
delineated using DEM processing in GIS. The derived extent of HYDROMET ux tower at Dripsey. Because of the small size of the
stream was validated against the regional channel map. Geomor- catchment (area 15 km2) and absence of any other precipitation
phic properties of the channel were dened based on the DEM data. data ne enough to resolve the heterogeneities within the water-
All thematic maps were resampled at 50 m  50 m spatial shed, the rainfall was assumed to be uniform within the catchment.
 T. Zi et al. / Environmental Modelling & Software 83 (2016) 310e325 315

Fig. 3. Observed monthly variations in precipitation (P), discharge volume (Q) and suspended sediment losses (SSY) in the Dripsey catchment.

The assumption is reasonable given the mild topographic relief and The model accounts for erosion, deposition, and transport of
a uniform land cover within the catchment. The stream ow and suspended sediment on the hillslopes. However, considering that
suspended sediment concentration data, collected by the water the bank and bed erosion and deposition in rst order river chan-
level recorder and water ISCO 6712 auto-sampler at the outlet of nels are not signicant (Golubev, 1982), these processes were not
catchment, were used to both calibrate and validate the model. included in sediment dynamics calculations in the channel. All the
sediment entering the channels was assumed to directly reach the
outlet of the watershed. It is to be acknowledged that this
3.3. Model implementation in dripsey catchment assumption may cause bias in suspended sediment yield estimates,
especially from large basins wherein river beds and banks can be a
The integrated model simulations were performed at signicant source and/or sink of sediment.
50 m  50 m spatial resolution. Simulation results were output at
hourly temporal resolution. While overland ow (from either
saturation excess or inltration excess) is generated over the entire 3.4. Sensitivity analysis
land surface cell in the original GEOtop model, a rill width ratio
(wdx) parameter was introduced in the sediment dynamics model A global multivariate sensitivity analyses was rst conducted to
to account for the ow organization within each cell. The evaluate the role of physics-based parameters that may impact soil
assumption here is that the generated overland ow is transferred erosion and deposition. The methodology is based on a Monte Carlo
into small rills (Fig. 2) in which overland ow gets concentrated and framework, and can be used to analyze the inuence of a parameter
the potential of rill erosion is high. The rill width ratio (wdx) was while also considering the inuence of all other parameters at the
calculated using the following equation: same time (Franks et al., 1997). The parameters chosen for the
sensitivity analysis included hydrologic properties that may inu-
P
i wdi
ence runoff generation (Table 2) and land properties that may in-
wdx (14) uence soil erodibility (Table 3). Ranges of some of the soil
xw
parameter such as saturated conductivity for different layers (K),
where i is the rill index within a cell, wdi is the width of ith rill [m], soil residual water content (qr), wilting point (qw), eld capacity
and xw is the width of the cell [m]. In other words, wdx is the fraction (qfc), saturated water content (qs), and soil median particle size (D50)
of a grid cell covered by overland ow, when overland ow is active. were assigned based on the sand-silt-clay fraction of in situ soil
The soil erosion calculations were accordingly modied to account samples (Lewis, 2011). The LAI range was assigned based on the
for erosion under a redistributed overland ow regime. Rain splash land cover of the catchment. Other parameters (Table 2, Table 3),
detachment occurred in the entire cell while the ow detachment such as Chezy's roughness coefcient (Cm), rill width ratio (wdx),
only occurred within the rills. While rill characteristics can be and root area ratio (Ra) were assigned a conservative range large
observed, the absence of relevant data for the site leads us to enough to encompass the range of parameters used in previous
consider rill characteristics as a calibration parameter. This is an studies. Model simulations were conducted using a single rainfall
important area for future research attention. event for 10,000 random sets of parameters, which were sampled
316 T. Zi et al. / Environmental Modelling & Software 83 (2016) 310e325

Table 2
Ranges and numerical sensitivity index (NSI) of hydrologic parameters.

Variable name Description Range NSI

LI Leaf area index 0.1e4 0.073

CF Canopy fraction 0e1 0.0009
Cm Chezy's roughness coefcient 0.01e5 0.0325
wdx Rill width ratio 0.01e1 0.073
D1 [mm] Soil depth of rst layer 25e300 0.0731
Kh_ L1 [mm s
1] Horizontal saturated conductivity (1st layer) 0.00036e0.12 0.0008
Kv_ L1 [mm s
1] Vertical saturated conductivity (1st layer) 0.00036e0.06 0.0322
qr _L1 Residual water content (1st layer) 0.03e0.06 0.0033
qw_L1 Wilting point (1st layer) 0.06e0.16 0.0040
qfc_L1 Field capacity (1st layer) 0.3e0.5 0.0050
qs_ L1 Saturated water content (1st layer) 0.3e0.6 0.0045
Kh_ L2 [mm s
1] Horizontal saturated conductivity (2nd layer) 0.00036e0.12 0.0012
Kv_ L2 [mm s
1] Vertical saturated conductivity (2nd layer) 0.00036e0.06 0.0001
qr _L2 Residual water content (2nd layer) 0.03e0.06 0.0035
qw_L2 Wilting point (2nd layer) 0.06e0.16 0.0006
qfc_L2 Field capacity (2nd layer) 0.3e0.5 0.0013
qs_ L2 Saturated water content (2nd layer) 0.3e0.6 0.0062
Kh_ L3 [mm s
1] Horizontal saturated conductivity (3rd layer) 0.00036e0.12 0.0076
Kv_ L3 [mm s
1] Vertical saturated conductivity (3rd layer) 0.00036e0.06 0.0119
qr _L3 Residual water content (3rd layer) 0.03e0.06 0.0026
qw_L3 Wilting point (3rd layer) 0.06e0.16 0.0006
qfc_L3 Field capacity (3rd layer) 0.3e0.5 0.0048
qs_ L3 Saturated water content (3rd layer) 0.3e0.6 0.0005

Table 3 The sensitivity of hydrologic and land surface parameters is

Ranges and numerical sensitivity index (NSI) of erosion and sediment dynamics presented individually here. Among the 23 parameters selected for
parameters.
the hydrological sensitivity analysis, leaf area index (LAI), rill width
Variable name Description Range NSI ratio (wdx), depth of the rst soil layer (D1), Chezy's roughness co-
Cm Chezy's roughness coefcient 0.01e5 0.3569 efcient (Cm), and vertical hydraulic conductivity of the rst soil
Rs [MPa] Root tensile strength 0.001e0.1 0.0236 layer (Kv_L1) were found to be among the ve most sensitive pa-
Ra Root area ratio 0.001e0.1 0.0033 rameters for runoff generation (Table 2). Similarly, Cm, wdx and D50
wdx Rill width ratio 0.01e1 0.3569
were observed to be the three most sensitive parameters for sedi-
D50 [mm] Median particle size 5e200 0.3563
z [kPa] Soil cohesion 0e30 0.0237 ment yield simulations (Table 3). Root tensile strength (Rs) and soil
cohesion (z) also were important players in sediment yield simu-
lations, though to a lesser extent than the other three parameters
identied above. Among the three most sensitive parameters for
from uniform distributions across the specied parameter ranges.
sediment yield, D50 can be estimated from soil texture and Cm can
Based on the rank of Nash-Sutcliffe coefcients (NSE) for each
be related to land cover and land use. wdx can be highly varied
parameter set, the 10,000 parameter sets were divided into ten
between different watersheds or land parcels depending on micro-
performance classes, each with 1000 parameter sets. The classes
topographic heterogeneity arising from natural causes or imple-
were ranked by NSE from low to high, i.e. Class 1 consisted of
mentation of agricultural practices (e.g. tillage farming). The other
parameter sets with lowest NSE while Class 10 comprised of
two sensitive parameters for sediment yield simulations: root
parameter sets with highest NSE. A numerical sensitivity index
tensile strength (Rs) and soil cohesion (z), are correlated, as Rs
(NSI) was used to quantify the sensitivity (Montaldo et al., 2003):
proportionally affects z (see Equations (7)e(9)) and both the pa-
rameters have a tendency to reduce soil erosion.
1 X
nc X
N
Having identied the most sensitive parameters, next a uni-
NSI xd i; j (15)
nc Nxrange variate sensitivity experiment was performed to investigate the
i1 j1
role of two sensitive but independent parameters (wdx and z) on
where xd is the difference between the parameter values of two erosion simulations. This involved forcing the model with a 32 mm
classes for the jth cumulative frequency and for ith pair of cumu- precipitation event of 24 h duration. As shown in Fig. 4, increasing
lative frequency curves, nc is the total number of pairs of cumulative wdx caused a reduction in suspended sediment losses. This is
frequency classes (with 10 classes, nc 45), and N is the number of because for a wider rill, overland ow depth and velocity is smaller
cumulative frequency values in each performance class (N 1000 resulting in a smaller shear stress between water and soil surface to
in this case). xrange is the range of parameter values. NSI is an in- detach the soil particles. It is to be noted that the relative changes in
dicator of the distance between each class. A larger value of NSI for a wdx (base value equals to 0.5) were larger than the relative changes
parameter indicates higher sensitivity to it. It is to be noted that in SSY. For example, 20% decrease in wdx increased sediment yield
while generating the parameter sets, extra care was taken to ensure by 15%. This was due to the offset effect of reduction in ow
that the considered values of dependent parameters were physi- detachment area when the rills are narrower.
cally realistic. For example, only the parameter sets with soil re- Five additional sensitivity experiments with varying rainfall
sidual water content (qr), wilting point (qw), eld capacity (qfc) and intensities, but with identical precipitation amount, were con-
saturated water content (qs) in increasing order were considered. It ducted to evaluate potential controls of event characteristics on the
is also possible that only certain ranges of parameters may co-exist. sensitivity of sediment yield to wdx. Different rainfall intensities did
However, because of the absence of any parameter-concurrency not show signicant impacts on sensitivity analysis, indicating that
data for the region, the sensitivity analyses did not consider this. ow detachment rather than splash erosion, was the major driving
 T. Zi et al. / Environmental Modelling & Software 83 (2016) 310e325 317

Fig. 4. Sensitivity of rill width ratio, wdx (a) and soil cohesion, z (b) on suspended sediment yield (SSY) for a range of precipitation intensities.

force for water erosion in Dripsey catchment. This is possibly 0 indicates no agreement at all. Krause et al. (2005) concluded that
caused by the dense grassland vegetation cover within Dripsey IOA of 0.65 could indicate a good model performance, although IOA
catchment which protected the soil surface from detachment by is not sensitive to systematical under- or over-estimation:
raindrops. Similar sensitivity experiment for soil cohesion (z) Pn 2
showed that the sediment yield decreased with increase in z. The i1 Oi
Si
IOA 1
   2 : (18)
inuence of soil cohesion on erosion was the largest when cohesion Pn Si
O Oi
O
value was between 10 and 20 kPa. Notably, as the precipitation i1
intensity increased, the sensitivity to soil cohesion decreased.
RMSE quanties prediction error as:
It is to be noted that sediment yield simulations could be
inuenced both by parameters that control rainfall-runoff pro- s
Pn 2
cesses and parameters that dene the land surface and soil char- i1 OI
Si
RMSE : (19)
acteristics. This makes calibration of a physically based integrated N
erosion-deposition models much more challenging than inte-
A RMSE value of zero indicates perfect t. If the value is less than
grated hydrologic models.
half of the standard deviation of the observations, the model per-
formance may be considered as good (Singh et al., 2005).
3.5. Model performance during calibration and validation periods
3.5.2. Calibration period
3.5.1. Evaluation metrics The calibration was carried out using hourly ow data and
Performance of the model was evaluated using four metrics: suspended sediment concentration data from Jan 1st, 2002 to Mar
Relative bias (RB), NasheSutcliffe efciency (NSE), index of agree- 14th, 2002 (73 days, 10% of the observed data). The GEOtopSed
ment (IOA), and Root mean square error (RMSE). RB evaluates the model spin-up was performed using the meteorological data of Dec
relative distances between observed (O) and simulated (S) data 2001. The calibration process mainly focused on the most sensitive
with respect to the observed mean (O) using: parameters, as identied in the global sensitivity analyses. Cm, D1,
Pn and wdx were the three major parameters calibrated for the ow
i1 Si
Oi simulation. z and D50 were mainly calibrated for erosion simulation.
RB 100 (16)
nO The initial values of the calibration parameters were set to the
median of the ranges in Tables 2 and 3. The model calibration was
where n is the length of data series. NSE (Nash and Sutcliffe, 1970) is
based on both mass balance and goodness of t of hydrograph and
a dimensionless goodness-of-t indicator that ranges from nega-
sedigraph. The goal was to at least limit the relative bias (RB) of total
tive innity to one. A NSE value of 1 indicates a perfect t between
discharge volume and total SSY to less than 15% and to have NSE of
observed and simulated data.
daily average ow rate and daily total SSY to be larger than 0.4.
Pn 2 The observed and simulated stream ow series in calibration
i1 Oi
Si
NSE 1
2 : (17) period are shown in Fig. 5a and b. Observed and simulated total
Pn 
i1 Oi
O
discharge volumes per unit area over the calibration period were
411 mm and 472 mm, respectively (RB ~ 15%). The R2 and Nash-
In the context of watershed hydrologic modeling, Moriasi et al. Sutcliff coefcient (NSE) were 0.76 and 0.48, respectively. The in-
(2007) summarized that the model performance could be consid- dex of agreement (IOA) was 0.92 and the root mean square error
ered as Good if the RB is less than 15% for streamow and 30% for (RMSE) was 3.03 mm/day. The NSE, IOA, and RB values indicate a
sediment, and NSE is larger than 0.65 for monthly time step sim- good model performance based on the threshold proposed by
ulations. For ner time step, the equivalent goodness threshold of Moriasi et al. (2007) and Krause et al. (2005). The RMSE was 71% of
NSE is lower, e.g. NSE value of 0.65 for monthly time step may yield the standard deviation of observed streamow, which means that
a NSE value of 0.4 for daily time step (Fernandez et al., 2005). IOA the model performance was not good (Singh et al., 2005). Large
(Willmott, 1981) is used to detect the differences in the observed RMSE value is because of the model's tendency to overestimate high
and simulated means and variances, especially during the intense ow and underestimate low ow (Fig. 5b).
rainfall events as the metric is sensitive to extreme values. IOA In regards to the estimates of suspended sediment yield (SSY),
ranges from 0 to 1. A value of 1 indicates a perfect match, and the model also showed good agreement with the observations. The
318 T. Zi et al. / Environmental Modelling & Software 83 (2016) 310e325

Fig. 5. Model estimates of daily discharge and daily suspended sediment yield (SSY) during the calibration period. (a) plot of observed (with subscription o) and simulated (with
subscription s) daily average stream ow rate; (b) scatter plot of observed and simulated daily average stream ow rate; (c) plot of observed and simulated daily total SSY; (d)
scatter plot of observed and simulated daily total SSY. Black solid line in panel (b) and (d) is the 1:1 line. Blue dashed line is the best-t line.(For interpretation of the references to
colour in this gure legend,the reader is referred to the web version of this article.).

observed and modeled total suspended sediment losses were 72.76 The model performances of ow and SSY simulation during
and 79.94 tons, respectively, for the calibration period. The RB was calibration and validation periods are summarized in Table 4. The
9.87%. R2 and NSE were 0.79 and 0.65, respectively (Fig. 5c, d). The NSE and IOA values of modeled ow in calibration and validation
RMSE was 0.83 ton/day (56% of the standard deviation of observed periods were similar. The RMSE of simulated ow in validation
SSY) and IOA was 0.93. The model performance based on the NSE, period was smaller than in calibration period. The differences be-
IOA, and RB values can be termed as very good. The SSY during the tween observed and simulated ow rate were smaller during the
peak ow were overestimated, just like the overestimation of high low ows than high ows. This was partly because of the longer
ows by the hydrologic model. duration of low ow during the validation period. The model per-
formance of suspended sediment simulation in validation period
3.5.3. Validation period was not as good as in calibration period, even though the perfor-
The model was validated using data from Mar 15th, 2002 to the mances of ow simulation in both periods were similar. This was
end of 2003. Observed and simulated total discharge per unit area due to the bias in the ow simulation and power-law relationship
were 1303.8 mm and 1287.7 mm, respectively. The time series of between ow rate and water transport capacity, which exacerbated
cumulative observed and simulated discharge volume matched the SSY estimates.
each other for most of the period, and the RB was only
1.2%. In In order to investigate the performance of the integrated model
general, the model reasonably simulated both the timing and during markedly different ow regimes, a 20-day moving average
magnitude of stream ow during the validation period (RMSE: discharge threshold of 0.2 m3/s was used to divide the simulation
1.47 mm/day or 56% of the standard deviation of observed series into wet and dry periods (see Fig. 6a). The total duration of
discharge; NSE: 0.5; and IOA: 0.9). However, the model simulation the dry period was 173 days, almost 27% of the two years. The total
overestimated high ows and underestimated low ows, much like precipitation during the dry period was 326.6 mm, which is almost
it performed during the calibration period. 11% of the total precipitation. Table 5 summaries the performances
The performance of SSY simulation registered similar biases as
that shown in streamow simulation. The observed and simulated
Table 4
total annual suspended sediment losses were 170.66 and 132.48
Performance of stream ow and total suspended sediment yield simulation in
tons, respectively. The total simulated SSY for validation period had calibration and validation periods.
a relative bias of
22.37%. The erosion model did not miss any
Calibration Validation
erosion events during the validation period. Similar to the ow
simulation results, SSY estimates were somewhat larger than Q SSY Q SSY
observed during high ows and smaller during low ows. The RMSE NSE 0.48 0.65 0.5 0.32
was 0.41 ton/day (53% of the observed SSY) and the value of IOA was IOA 0.92 0.93 0.9 0.86
0.86. NSE was 0.32. Relatively poor NSE is because of the metric's RB 15% 9.87%
1.2%
22.37%
RMSE 3.03 0.83 1.47 0.41
sensitivity to overestimation of SSY peaks.
 T. Zi et al. / Environmental Modelling & Software 83 (2016) 310e325 319

Fig. 6. Model estimates of discharge and suspended sediment yield (SSY) during the validation period. (a) plot of observed and simulated daily average stream ow rate. The wet
and dry periods were identied by letter W and D; (b) scatter plot of daily average stream ow rate; (c) plot of observed and simulated daily total SSY; (d) scatter plot of daily total
SSY. Black solid line in panel (b) and (d) is the 1:1 line. Blue dashed line is the best-t line.(For interpretation of the references to colour in this gure legend,the reader is referred to
the web version of this article.).

of the integrated model for both ow and SSY simulations during intensity (Fig. 7a). The R2 between observed daily SSY and daily
wet and dry periods. Both simulations showed better agreement rainfall was only 0.14. The R2 between observed daily SSY and the
with observed data during wet periods than in dry periods. In product of daily total rainfall kinetic energy and maximum 30-min
general, discharge volume was overestimated for dry periods and rainfall intensity (EI30) was relatively better but still modest, and
underestimated in wet periods, while the SSY was underestimated was equal to 0.53. It is to be noted here that the lagged correlation
in both wet and dry periods. It is to be noted however that for large magnitude between observed daily SSY and daily rainfall was
peaks, both ow rate and SSY was overestimated by the model. The maximum for a zero day lag. The large variations in observed daily
negative NSE and RB indicate that the calibrated parameter set did SSY indicate that the daily precipitation magnitude and intensity
not provide accurate enough estimates of ow and SSY during dry only have a marginal potential for estimating daily SSY. The varia-
periods. This is in line with the ndings of Jetten et al. (1999) that tion could be because of the non-linear integrated response of the
showed erosion-deposition models cannot guarantee good per- catchment, which may result in generation of markedly different
formance, especially during dry periods, if the events used for runoff and hence SSY for the same daily precipitation amount and
calibration are not representative of the prevailing ow regime. It is EI30. Since runoff generation is directly inuenced by the soil
to be noted that even though the suspended sediment simulation moisture in the surface soil layer (which in turn is inuenced by
was not good enough during dry period, the bias of suspended meteorological forcings, antecedent hydrologic states, and water-
sediment is only 9.94% of the total bias. shed properties), we explore if the variations in daily SSY with EI30
could be directly related to the prevailing soil moisture. To this end,
4. Results and analyses data in the SSY vs EI30 plot are color coded based on simulated daily
average soil moisture (Fig. 7a). Based on the rst glance, for an
4.1. Temporal variations in modeled estimates identical EI30, SSY value appears to be larger for larger catchment-
average soil moisture.
SSY from the catchment varied markedly from day to day and To explore this further, two daily events of similar magnitude
also non-monotonically with daily precipitation magnitude and but with very different moisture conditions were considered
(Fig. 8). The rst event occurred on 2002/1/27 (identied in Fig. 5)
Table 5 and delivered a total precipitation of 11 mm, a total measured
Performance of stream ow and suspended sediment simulations in wet and dry streamow of 1.3  105 m3 (modelled: 1.8  105 m3) and an event
periods. total suspended sediment yield (SSY) of 2.1 ton (modelled: 3.4 ton).
QW Qd SSYw SSYd The second event occurred on 2002/7/2 (identied in Fig. 6) and
delivered almost identical rainfall (11.6 mm in total) as the rst
NSE 0.33
0.75 0.25
1.52
IOA 0.88 0.59 0.85 0.44 event, but the total streamow response and SSY at the outlet were
RB 4.55% 57.58%
20.65%
95.76% much different. Streamow response for the second event was
Bias
56.32 39.15
34.44
3.80 measured at 1.1  104 m3 (modelled: 1.8  104 m3), while the event
320 T. Zi et al. / Environmental Modelling & Software 83 (2016) 310e325

Fig. 7. (a) Scatter plot of daily observed SSY and EI30, color coded with daily simulated catchment-average soil moisture; the inset gure is a zoom-in of daily observed SSY between
0 and 1 ton/d and EI30 between 0e1.2  104 MJ*mm*ha
1*h
1; (b) Scatter plot of daily observed SSY from 2002 to 2003 and daily simulated catchment average soil moisture.

mean SSY was 0.066 ton (modelled: 0.054 ton). Notably, even responses (both streamow and SSY) were negligible for events
though the rainfall volumes and intensities were similar for these from day-of-year (DOY) 200 to 290 in 2002, as the soil moisture was
two selected events, the rst event lead to much larger runoff and generally below saturation during this period (Fig. 9c). By day 290,
SSY at the outlet (Fig. 8a). The difference was mostly because of the water table was near the soil surface, which resulted in any
strong contrast in the antecedent moisture condition between the additional precipitation to cause saturation excess runoff. To eval-
two cases (Fig. 8c and f). Around 88% of the catchment area was uate the variation of SSY with soil moisture across all events, a
saturated at the start of event 1, with this fraction increasing to scatter plot between SSY and spatially averaged moisture condi-
around 98% during the event. As a result almost the entire catch- tions in the top soil layer (of thickness 10 cm) was drawn (Fig. 7b).
ment participated in runoff generation leading to larger soil loss The gure suggests that runoff generation (from saturation excess
and larger mean SSY. This is not surprising as the catchment is fairly or inltration excess) and hence erosion per unit event magnitude
small. In contrast, only around 18% of the catchment area was is indeed larger when the top soil layer in the catchment is near
saturated before event 2 and the fraction increased to around 36% saturation. These results suggest that accurate prediction of spatial
during the event. Fig. 9b and d further reinforces the narrative that distribution of soil moisture is critical for generating temporally
surface soil moisture inuences SSY as it shows that event ne estimates of SSY. Notably, increase in soil saturation also

Fig. 8. The hourly rainfall, observed discharge and SSY of two selected events (a,d); the simulated spatial distribution of soil loss area for two selected events (b,e), area with soil
losses were black; the changes of saturated ratio of the simulated rst layer soil moisture of the whole catchment (e,f); the plots for rst event were on the top.
 T. Zi et al. / Environmental Modelling & Software 83 (2016) 310e325 321

increases soil cohesion which in turn may reduce erosion (see fraction were larger in wet periods than in dry periods. The varia-
Equation (7)), but its impact is negated by large runoff generation tions of erosion and deposition areal portion were larger in wet
for higher soil saturation cases. It is to be noted that the soil satu- periods (Table 6). Notably, the erosion and deposition areal frac-
ration of the top soil layer is a function of coupled surface water- tions in the catchment were as large as 79% and 36% of the catch-
groundwater-evapotranspiration interactions, and the parameters ment area during the two year simulation periods. These statistics
(such as topography, subsurface soil property, land cover etc.) that highlight that the source and sink areas in a catchment are very
inuence these processes. As such, the expressed role of surface soil dynamic and change at both event and seasonal scales. Also, given
moisture on SSY indirectly highlights the need for modeling of the role of soil saturation area on runoff/erosion generation, dy-
coupled processes in both space and time, much along the lines of namic mapping of source/sink areas can benet from spatially-
process representations implemented in GEOtop. explicit simulations of coupled hydrologic processes.
Suspended sediment yield from the catchment also showed At the annual scale, most deposition occurred near the riparian
ample variations between wet and dry periods. Even though the zone of the stream (dark blue color in Fig. 10a), due to its relatively
wet period (see Fig. 6a) spanned only around 73% of the simulation at topography. Erosive losses were mainly located in the transition
time, it delivered 98.4% (99.9%, modeled) of the total sediment area between the hillslope and the riparian zone (dark red color in
yield. This is in line with previously reported results from both eld Fig. 10a), which had both large ow accumulation area and rela-
(Fu, 1989) and modeling experiments (Baartman et al., 2012; Zhang tively steep slope. The locations of erosion and deposition in terms
et al., 2012) that have shown that a few extreme events may of slope and specic catchment area are shown in Fig. 10c and d.
contribute to a large portion of annual total soil erosion. Notably, Specic catchment area is the ratio of total contributing area of
our research area has a temperate weather and lacks extreme each land pixel and the pixel's width perpendicular to the ow
precipitation events, but nonetheless the dry period (~27% of the direction. The bar plot of average erosion/deposition per unit area
total simulation time) delivered only ~1% of the total SSY. for different specic catchment areas shows that in general, erosion
and deposition rates per unit area were larger with increasing
4.2. Spatially distributed estimates of erosion and deposition specic catchment area. Areas with larger contributing areas
generally had larger overland ow and hence larger transport ca-
Estimates of erosion and deposition simulated by the model pacity thus allowing more entrainment of soil particles (Fig. 10c).
displayed signicant heterogeneity in both space and time. For Additionally, as more suspended sediment was transported to areas
example, for the two selected precipitation events with similar with large contribution area, the possibility that the suspended
amounts of rainfall (Fig. 8a, d), the areal extent of soil loss and its sediment concentration (SSC) exceeded the transport capacity,
spatial distribution locations were very different (Fig. 8b, e). The especially at atter locations, also increased. This explains the in-
percentage area that participated in erosion loss was as high as 72% crease in deposition with larger contributing areas (Fig. 10d).
in event 1, while reaching only around 14% for event 2. Similarly the Notably, total simulated soil erosive losses per unit area (Es) did not
areal fraction of deposition areas were 14% and 7% respectively for show a monotonic trend with slope (Fig. 10c). For atter slopes, Es is
the two events. A larger erosion loss area for event 1 can again be expected to be small because of the lower ow velocities. For larger
attributed to higher antecedent surface soil moisture in the slopes, while the local erosion at the cell under consideration is
catchment, which leads to a larger fraction of catchment area expected to be large, the incoming sediment input may either be
generating runoff. A larger volume of overland ow resulted in, small or large for a given contributing area depending on the
higher shear stress and hence more soil loss from larger area. For spatio-temporal extent of soil saturation area that determines
the two years, the mean areal fraction of erosion and deposition runoff generation. As a result, Es, which is an integrated net-
was 27% and 8% respectively. Both erosion and deposition areal erosional response from the contributing area and the pixel

Fig. 9. (a) the integrated model is able to capture occurrence of runoff events after a dry periods. The red box identies the periods used in panel b, c, and d; (b)observed and
simulated discharge per unit area; (c). observed and simulated soil moisture (top 25 cm) (d)observed and simulated daily SSY.
322 T. Zi et al. / Environmental Modelling & Software 83 (2016) 310e325

Table 6
The mean and standard deviation of erosion and deposition areas.

Mean erosional Mean depositional Standard deviation of erosional Standard deviation of depositional
area fraction area fraction area fraction area fraction

2002e2003 0.27 0.08 0.246 0.0859

Wet periods 0.33 0.11 0.2476 0.0835
Dry periods 0.11 0.02 0.1519 0.0585

Fig. 10. (a) Spatial map of soil erosion and deposition for each land pixel (unit:ton/ha/yr) for the entire two year simulation period; (b) Scatter plot between erosive loss (Es) and
slope-length factor for each land pixel; (c) Color chart of the logarithm of simulated erosive loss (log (ton)) for a range of slope (on x-axis) and specic ow accumulation area (on y-
axis). Also included on the two axes are the bar plots of total erosive losses per unit area for different slope and specic catchment area classes; (d) Color chart of the logarithm of
simulated deposition (log (ton)). Also included on the two axes are the bar plots of total deposition per unit area for different slope and specic catchment area classes. Note: grey
cells in (c) and (d) are used to identify joint classes with no data.

under consideration, may show non-monotonic variation with deposition module for a 3D surface-subsurface hydrologic model,
slope. This non-monotonic variation also has a tangible impact on GEOtop, and evaluated its applicability in explaining the spatio-
the relation between the length-slope (LS) factor and the total temporal distribution of erosion, deposition, and sediment yield
simulated soil erosive losses per unit area (Es). Model results sug- dynamics at both plot and catchment scale. The model uses a
gest that while Es shows an overall increasing trend with LS physically-based representation of coupled surface-subsurface hy-
(Fig. 10b), the relation is not one to one, and shows ample variations drologic processes and a comprehensive land-surface energy and
around the Es-LS line (R2 0.47, s2 5.44). While the result is based water interaction scheme to simulate hydrologic response and
on only two years of simulation, it does highlight that the source consequent sediment dynamics. Because of its fully distributed
efciency or the erosion rates (in homogeneous watersheds/hill- nature, the model can account for spatial heterogeneity in the
slopes with uniform soil property) is not just a simple function of watershed. As the model runs at an adaptive time interval deter-
contribution area and slope, and can be strongly inuenced by the mined by the dynamics of states, it can be used to perform simu-
dynamic spatio-temporal distribution of soil moisture conditions. lations at event to inter-decadal scales. At the plot level, the model
was evaluated at event time scale. Results show that the model was
5. Summary and conclusions able to capture both the trend and the quantity of runoff and sus-
pended sediment concentration (SSC) in response to events of
In this study, we developed an open-source sediment erosion/ varying intensity and duration. At the catchment level, after
 T. Zi et al. / Environmental Modelling & Software 83 (2016) 310e325 323

performing a global multivariate sensitivity analyses to identify parameters (e.g. rill width, soil cohesion etc.), and by implementing
sensitive parameters for hydrologic modeling (such as leaf area the model in varied settings. The presented model version does not
index, rill width ratio, depth of the rst soil layer, Chezy's roughness account for bank-erosion processes, which limits its applicability to
coefcient, and vertical hydraulic conductivity of the rst soil layer) watersheds with small bank erosion w.r.t. total hillslope losses. It is
and sediment dynamics (such as Chezy's roughness coefcient, soil to be noted that 1D representation of ow routing in channel
median particle size, and rill width ratio), the model performance especially limits the calculation of lateral convective gradients and
was evaluated at daily and seasonal time scales. In general, the shear, thus hindering the development of a comprehensive sedi-
model reasonably simulated both the timing and magnitude of ment erosion and transport module. Also, as the GEOtop model
stream ow and suspended sediment yield (SSY) from the catch- uses kinematic wave scheme to solve for overland and channel
ment. Both streamow and suspended sediment yield simulations ow, the model is not well suited for ow and sediment yield
showed better agreement with observed data during wet periods calculations in ow regimes with either Froude number smaller
than in dry periods. The discharge volume was overestimated for than 0.5 or kinematic wave number smaller than 5 (Vieira, 1983).
dry periods and underestimated in wet periods, while the SSY was For using the model as a predictive tool, vegetation dynamics in
underestimated in both wet and dry periods. Notably, the biases in response to changes in meteorological forcings and hydrologic
suspended sediment yield simulation were similar to that shown in states need to be accounted for. In spite of aforementioned limi-
streamow simulation, thus highlighting that accuracy of ow tations, the open-source integrated modeling framework presented
simulations critically inuences the estimation accuracy of SSY. here offers the potential for its use both as an evaluation and
While our research area has a temperate maritime climate and retrospective-prediction tool, and as a virtual laboratory for un-
lacks extreme precipitation events, but still dry period (~27% of the derstanding the role of hydrologic states and parameters on sedi-
total simulation time) delivered only ~1% of the SSY. This un- ment dynamics.
derscores the importance of obtaining robust calibration parameter
sets that at least perform well during wet periods, as is well Acknowledgements
documented in nutrient export studies (Jordan et al., 2005; Nasr
et al., 2007; Ye et al., 2012). It is to be noted that the sensitivity We appreciate the valuable feedbacks from reviewers and the
analysis performed for Dripsey catchment could serve as a refer- editor, which greatly improved this manuscript. This study was
ence for future model applications, and for prioritizing observation supported by the Irish Environmental Protection Agency (EPA)
of parameters to reduce the uncertainties in the erosion model. under the Science Technology Research & Innovation for the
Observed SSY from the catchment showed a non-monotonic Environment (STRIVE) Programme 2007e2013 of Ireland (Soil H:
variation with daily precipitation magnitude. Further examination Interactions of soil hydrology, land use and climate change and
of simulation results revealed that SSY per unit event magnitude their impact on soil quality; 2007S-SL-2-S1). Mukesh Kumar
varied proportionally with the prevailing soil moisture. The result acknowledges the support of National Science Foundation CZO
indicates that accurate prediction of spatial distribution of soil grant (EAR-1331846). We acknowledge the developers of GEOtop
moisture is critical for generating temporally ne estimates of SSY. model, for keeping the software open-source and free.
Larger SSY per unit event magnitude for higher soil saturation also
indicates that large runoff generation for wetter soil moisture References
conditions negated the impact of increase in soil cohesion with
moisture content. The simulation results also showed that the Albertson, J.D., Kiely, G., 2001. On the structure of soil moisture time series in the
context of land surface models. J. Hydrol. 243, 101e119. http://dx.doi.org/
source (erosion) and sink (deposition) areas in a catchment were 10.1016/s0022-1694(00)00405-4.
heterogeneous and dynamic, and could change from one event to Baartman, J.E.M., Jetten, V.G., Ritsema, C.J., de Vente, J., 2012. Exploring effects of
the next. Again, the extent of source/sink areas were found to be rainfall intensity and duration on soil erosion at the catchment scale using
openLISEM: Prado catchment, SE Spain. Hydrol. Process. 26, 1034e1049. http://
inuenced by prevailing moisture conditions, which in turn
dx.doi.org/10.1002/hyp.8196.
determined the quantity of runoff generation. Model results also Bertoldi, G., 2004. The water and energy balance at basin scale: a distributed
suggest that long-term erosion rate from a location was not a modeling approach. Environmental Engineering, University of Trento, p. 202.
Bertoldi, G., Della Chiesa, S., Notarnicola, C., Pasolli, L., Niedrist, G., Tappeiner, U.,
simple function of slope-length. In fact, the relation between
2014. Estimation of soil moisture patterns in mountain grasslands by means of
erosion and slope-length showed ample variations, thus high- SAR RADARSAT2 images and hydrological modeling. J. Hydrol. 516, 245e257.
lighting that the source efciency or the erosion rates (in homo- http://dx.doi.org/10.1016/j.jhydrol.2014.02.018.
geneous catchments/hillslopes with same rainfall forcing and Bullock, M.S., Kemper, W., Nelson, S., 1988. Soil cohesion as affected by freezing,
water content, time and tillage. Soil Sci. Soc. Am. J. 52, 770e776.
uniform soil property) can be severely inuenced by the dynamic Cabral, M.C., Garrote, L., Bras, R.L., Entekhabi, D., 1992. A kinematic model of inl-
spatio-temporal distribution of soil moisture conditions. Since the tration and runoff generation in layered and sloped soils. Adv. Water Resour. 15,
spatio-temporal dynamics of soil moisture is dependent on coupled 311e324. http://dx.doi.org/10.1016/0309-1708(92)90017-V.
Chen, X., Kumar, M., McGlynn, B.L., 2014. Variations in streamow response to large
process interactions such as evapo-transpiration, capillary rise and Hurricane-season storms in a southeastern U.S. Watershed. J. Hydrometeorol.
lateral groundwater ow, aforementioned results make a compel- 16, 55e69. http://dx.doi.org/10.1175/JHM-D-14-0044.1.
ling case for spatially-explicit simulations of coupled hydrologic de Vente, J., Poesen, J., Verstraeten, G., Van Rompaey, A., Govers, G., 2008. Spatially
distributed modelling of soil erosion and sediment yield at regional scales in
processes for estimation of erosion/deposition distribution and Spain. Glob. Planet. Change 60, 393e415. http://dx.doi.org/10.1016/
sediment generation. j.gloplacha.2007.05.002.
The simulation results presented in this study do not account for Della Chiesa, S., Bertoldi, G., Niedrist, G., Obojes, N., Endrizzi, S., Albertson, J.D.,
Wohlfahrt, G., Ho rtnagl, L., Tappeiner, U., 2014. Modelling changes in grassland
uncertainty in parameters. Also, even though the parameters in the hydrological cycling along an elevational gradient in the Alps. Ecohydrol. 7,
presented model are physically-based and can be obtained through 1453e1473. http://dx.doi.org/10.1002/eco.1471.
measurements, because of the sparseness of observed data, DeRoo, A.P.J., Wesseling, C.G., Ritsema, C.J., 1996. LISEM: a single-event physically
based hydrological and soil erosion model for drainage basins .1. Theory, input
disconnect between observation and model scale, and model un-
and output. Hydrol. Process. 10, 1107e1117.
certainty, these parameters still need to be calibrated. This is a big Endrizzi, S., Gruber, S., Dall'Amico, M., Rigon, R., 2014. GEOtop 2.0: simulating the
challenge for spatially-distributed models such as GEOtopSed, as combined energy and water balance at and below the land surface accounting
they are computationally demanding. Further condence in the for soil freezing, snow cover and terrain effects. Geosci. Model Dev. 7,
2831e2857. http://dx.doi.org/10.5194/gmd-7-2831-2014.
modeled estimates could be built by obtaining eld estimates of Fernandez, G.P., Chescheir, G.M., Skaggs, R.W., Amatya, D.M., 2005. Development
states (e.g. ground water, residence time etc.) and related and testing of watershed-scale models for poorly drained soils. Trans. ASAE 48,
324 T. Zi et al. / Environmental Modelling & Software 83 (2016) 310e325

639e652. Montaldo, N., Toninelli, V., Albertson, J.D., Mancini, M., Troch, P.A., 2003. The effect
Foster, G., Yoder, D., Weesies, G., McCool, D., McGregor, K., Bingner, R., 2005. Revised of background hydrometeorological conditions on the sensitivity of evapo-
Universal Soil Loss Equation Version 2. Science Documentation. (Draft). USDA- transpiration to model parameters: analysis with measurements from an Italian
ARS. alpine catchment. Hydrol. Earth Syst. Sci. 7, 848e861. http://dx.doi.org/10.5194/
Franks, S.W., Beven, K.J., Quinn, P.F., Wright, I.R., 1997. On the sensitivity of soil- hess-7-848-2003.
vegetation-atmosphere transfer (SVAT) schemes: equinality and the problem Morgan, R.P.C., Quinton, J.N., Smith, R.E., Govers, G., Poesen, J.W.A., Auerswald, K.,
of robust calibration. Agric. For. Meteorol. 86, 63e75. http://dx.doi.org/10.1016/ Chisci, G., Torri, D., Styczen, M.E., 1998. The European Soil Erosion Model
S0168-1923(96)02421-5. (EUROSEM): a dynamic approach for predicting sediment transport from elds
Fu, B., 1989. Soil erosion and its control in the loess plateau of China. Soil Use and small catchments. Earth Surf. Process. Landforms 23, 527e544.
Manag. 5, 76e82. http://dx.doi.org/10.1111/j.1475-2743.1989.tb00765.x. Moriasi, D., Arnold, J., Van Liew, M., Bingner, R., Harmel, R., Veith, T., 2007. Model
Ghimire, C.P., Bonell, M., Bruijnzeel, L.A., Coles, N.A., Lubczynski, M.W., 2013. evaluation guidelines for systematic quantication of accuracy in watershed
Reforesting severely degraded grassland in the Lesser Himalaya of Nepal: ef- simulations. Trans. ASABE 50, 885e900.
fects on soil hydraulic conductivity and overland ow production. J. Geophys. Nash, J.E., Sutcliffe, J.V., 1970. River ow forecasting through conceptual models part
Res. Earth Surf. 118, 2528e2545. http://dx.doi.org/10.1002/2013JF002888. I d A discussion of principles. J. Hydrol. 10, 282e290. http://dx.doi.org/10.1016/
Golubev, G.N., 1982. Soil erosion and agriculture in the world: an assessment and 0022-1694(70)90255-6.
hydrological implications. Hydrol. Sci. J. e J. Des. Sci. Hydrol. 27, 243e243. Nasr, A., Bruen, M., Jordan, P., Moles, R., Kiely, G., Byrne, P., 2007. A comparison of
Govers, G., 1990. Empirical Relationships for the Transport Capacity of Overland SWAT, HSPF and SHETRAN/GOPC for modelling phosphorus export from three
Flow. IAHS Publication, IAHS Press, Institute of Hydrology, pp. 45e63. catchments in Ireland. Water Res. 41, 1065e1073. http://dx.doi.org/10.1016/
Vorst, HAvd, 1992. BI-CGSTAB: a fast and smoothly converging variant of BI-CG for j.watres.2006.11.026.
the solution of nonsymmetric linear systems. SIAM J. Sci. Stat. Comput. 13, Niu, G.-Y., Pasetto, D., Scudeler, C., Paniconi, C., Putti, M., Troch, P., DeLong, S.,
631e644. http://dx.doi.org/10.1137/0913035. Dontsova, K., Pangle, L., Breshears, D., 2014. Incipient subsurface heterogeneity
Heppner, C.S., Ran, Q., VanderKwaak, J.E., Loague, K., 2006. Adding sediment and its effect on overland ow generationeinsight from a modeling study of the
transport to the integrated hydrology model (InHM): development and testing. rst experiment at the biosphere 2 landscape evolution observatory. Hydrol.
Adv. Water Resour. 29, 930e943. http://dx.doi.org/10.1016/ Earth Syst. Sci. 18, 1873e1883.
j.advwatres.2005.08.003. Orchard, C.M., Lorentz, S.A., Jewitt, G.P.W., Chaplot, V.A.M., 2013. Spatial and tem-
Heppner, C.S., Loague, K., VanderKwaak, J.E., 2007. Long-term InHM simulations of poral variations of overland ow during rainfall events and in relation to
hydrologic response and sediment transport for the R-5 catchment. Earth Surf. catchment conditions. Hydrol. Process. 27, 2325e2338. http://dx.doi.org/
Process. Landforms 32, 1273e1292. http://dx.doi.org/10.1002/esp.1474. 10.1002/hyp.9217.
Hessel, R., van den Bosch, R., Vigiak, O., 2006. Evaluation of the LISEM soil erosion Penna, D., Tromp-van Meerveld, H., Gobbi, A., Borga, M., Dalla Fontana, G., 2011. The
model in two catchments in the East African Highlands. Earth Surf. Process. inuence of soil moisture on threshold runoff generation processes in an alpine
Landforms 31, 469e486. http://dx.doi.org/10.1002/esp.1280. headwater catchment. Hydrol. Earth Syst. Sci. 15, 689e702.
Hessell, R., 2005. Effects of grid cell size and time step length on simulation results Ramsankaran, R., Kothyari, U.C., Ghosh, S.K., Malcherek, A., Murugesan, K., 2013.
of the Limburg soil erosion model (LISEM). Hydrol. Process. 19, 3037e3049. Physically-based distributed soil erosion and sediment yield model (DREAM)
http://dx.doi.org/10.1002/hyp.5815. for simulating individual storm events. Hydrol. Sci. J. 58, 872e891. http://
Hueso-Gonza
lez, P., Ruiz-Sinoga, J., Martnez-Murillo, J., Lavee, H., 2015. Overland dx.doi.org/10.1080/02626667.2013.781606.
ow generation mechanisms affected by topsoil treatment: application to soil Ran, Q., Su, D., Li, P., He, Z., 2012. Experimental study of the impact of rainfall
conservation. Geomorphology 228, 796e804. characteristics on runoff generation and soil erosion. J. Hydrol. 424, 99e111.
Ivanov, V.Y., Vivoni, E.R., Bras, R.L., Entekhabi, D., 2004. Catchment hydrologic Renard, K.G., Foster, G.R., Weesies, G.A., Porter, J.P., 1991. RUSLE: revised universal
response with a fully distributed triangulated irregular network model. Water soil loss equation. J. Soil Water Conserv. 46, 30e33.
Resour. Res. 40 http://dx.doi.org/10.1029/2004WR003218. Rigon, R., Bertoldi, G., Over, T.M., 2006. GEOtop: a distributed hydrological model
Jain, M., Kothyari, U., Raju, K., 2005. Gis based distributed model for soil erosion and with coupled water and energy budgets. J. Hydrometeorol. 7, 371e388.
rate of sediment outow from catchments. J. Hydraulic Eng. 131, 755e769. Rosenberg, E., Clark, E., Steinemann, A., Lettenmaier, D., 2013. On the contribution of
http://dx.doi.org/10.1061/(ASCE)0733-9429(2005)131:9(755). groundwater storage to interannual streamow anomalies in the Colorado
Jetten, V., de Roo, A., Favis-Mortlock, D., 1999. Evaluation of eld-scale and River basin. Hydrol. Earth Syst. Sci. 17, 1475e1491.
catchment-scale soil erosion models. CATENA 37, 521e541. http://dx.doi.org/ Safeeq, M., Mauger, G.S., Grant, G.E., Arismendi, I., Hamlet, A.F., Lee, S.-Y., 2014.
10.1016/S0341-8162(99)00037-5. Comparing large-scale hydrological model predictions with observed stream-
Jordan, P., Menary, W., Daly, K., Kiely, G., Morgan, G., Byrne, P., Moles, R., 2005. ow in the Pacic northwest: effects of climate and groundwater*.
Patterns and processes of phosphorus transfer from Irish grassland soils to J. Hydrometeorol. 15, 2501e2521.
riversdintegration of laboratory and catchment studies. J. Hydrol. 304, 20e34. Simoni, S., Zanotti, F., Bertoldi, G., Rigon, R., 2008. Modelling the probability of
http://dx.doi.org/10.1016/j.jhydrol.2004.07.021. occurrence of shallow landslides and channelized debris ows using GEOtop-
Jost, G., Schume, H., Hager, H., Markart, G., Kohl, B., 2012. A hillslope scale com- FS. Hydrol. Process. 22, 532e545. http://dx.doi.org/10.1002/hyp.6886.
parison of tree species inuence on soil moisture dynamics and runoff pro- Singh, J., Knapp, H.V., Arnold, J., Demissie, M., 2005. Hydrological modeling of the
cesses during intense rainfall. J. Hydrol. 420, 112e124. Iroquois River watershed using HSPF and SWAT. J. Am. Water Resour. Assoc. 41,
Julien, P., Simons, D., 1985. Sediment transport capacity of overland ow. Trans. 343e360.
ASAE 28, 755e762. Smith, R.E., Goodrich, D.C., Quinton, J.N., 1995. Dynamic, distributed simulation of
Kim, J., Ivanov, V.Y., Katopodes, N.D., 2013. Modeling erosion and sedimentation watershed erosion: the KINEROS2 and EUROSEM models. J. Soil Water Conserv.
coupled with hydrological and overland ow processes at the watershed scale. 50, 517e520.
Water Resour. Res. 49, 5134e5154. http://dx.doi.org/10.1002/wrcr.20373. Tao, J., Barros, A.P., 2014. Coupled prediction of ood response and debris ow
Krause, P., Boyle, D., Ba se, F., 2005. Comparison of different efciency criteria for initiation during warm- and cold-season events in the Southern Appalachians,
hydrological model assessment. Adv. Geosciences 5, 89e97. USA. Hydrol. Earth Syst. Sci. 18, 367e388. http://dx.doi.org/10.5194/hess-18-
Kumar, M., Duffy, C.J., Salvage, K.M., 2009. A second-order accurate, nite vol- 367-2014.
umeebased, integrated hydrologic modeling (FIHM) framework for simulation VanderKwaak, J.E., Loague, K., 2001. Hydrologic-Response simulations for the R-5
of surface and subsurface ow. Vadose Zone J. 8, 873e890. http://dx.doi.org/ catchment with a comprehensive physics-based model. Water Resour. Res. 37,
10.2136/vzj2009.0014. 999e1013. http://dx.doi.org/10.1029/2000WR900272.
Legates, D.R., Mahmood, R., Levia, D.F., DeLiberty, T.L., Quiring, S.M., Houser, C., Vieira, J.H.D., 1983. Conditions governing the use of approximations for the Saint-
Nelson, F.E., 2011. Soil moisture: a central and unifying theme in physical ge- Ve
nant equations for shallow surface water ow. J. Hydrol. 60, 43e58. http://
ography. Prog. Phys. Geogr. 35, 65e86. dx.doi.org/10.1016/0022-1694(83)90013-6.
Lewis, C., 2003. Phosphorus, Nitrogen and Suspended Sediment Loss from Soil to von Freyberg, J., Rao, P.S.C., Radny, D., Schirmer, M., 2015. The impact of hillslope
Water from Agricultural Grassland. Department of Civil and Environmental groundwater dynamics and landscape functioning in event-ow generation: a
Engineering, University College Cork. eld study in the Rietholzbach catchment. Switz. Hydrogeol. J. 1e14.
Lewis, C., 2011. Measurement and Modelling of Soil Hydrological Properties for Use Warner, S., Kiely, G., Morgan, G., O'Halloran, J., 2009. Does quantifying antecedent
in the Distributed Rainfall Runoff Model - GEOtop. Department of Civil and ow conditions improve stream phosphorus export estimation? J. Hydrol. 378,
Environmental Engineering, University College Cork. 97e104. http://dx.doi.org/10.1016/j.jhydrol.2009.09.009.
Lewis, C., Albertson, J., Xu, X., Kiely, G., 2012. Spatial variability of hydraulic con- Wicks, J.M., Bathurst, J.C., 1996. SHESED: a physically based, distributed erosion and
ductivity and bulk density along a blanket peatland hillslope. Hydrol. Process. sediment yield component for the SHE hydrological modelling system.
26, 1527e1537. http://dx.doi.org/10.1002/hyp.8252. J. Hydrol. 175, 213e238. http://dx.doi.org/10.1016/s0022-1694(96)80012-6.
Lewis, C., Albertson, J., Zi, T., Xu, X., Kiely, G., 2013. How does afforestation affect the Willmott, C.J., 1981. On the validation of models. Phys. Geogr. 2, 184e194.
hydrology of a blanket peatland? A modelling study. Hydrol. Process. 27, Wischmeier, W.H., Smith, D.D., 1978. Predicting Rainfall Erosion Losses: a Guide to
3577e3588. Conservation Planning. Dept. of Agriculture, Science and Education Adminis-
Mati, B.M., Morgan, R.P.C., Quinton, J.N., 2006. Soil erosion modelling with EURO- tration: for sale by the Supt. of Docs., US Govt. Print. Off.
SEM at Embori and Mukogodo catchments, Kenya. Earth Surf. Process. Land- Wu, T.H., McKinnell Iii, W.P., Swanston, D.N., 1979. Strength of tree roots and
forms 31, 579e588. http://dx.doi.org/10.1002/esp.1347. landslides on Prince of Wales Island, Alaska. Can. Geotechnical J. 16, 19e33.
Miguez-Macho, G., Fan, Y., 2012. The role of groundwater in the Amazon water http://dx.doi.org/10.1139/t79-003.
cycle: 1. Inuence on seasonal streamow, ooding and wetlands. J. Geophys. Yang, C.T., 1972. Unit stream power and sediment transport. J. Hydraulics Div. 98,
Res. Atmos. 117 (1984e2012). 1805e1826.
 T. Zi et al. / Environmental Modelling & Software 83 (2016) 310e325 325

Ye, S., Covino, T.P., Sivapalan, M., Basu, N.B., Li, H.-Y., Wang, S.-W., 2012. Dissolved practices. J. Soil Water Conserv. 67, 390e405. http://dx.doi.org/10.2489/
nutrient retention dynamics in river networks: a modeling investigation of jswc.67.5.390.
transient ows and scale effects. Water Resour. Res. 48 http://dx.doi.org/ Zimmermann, A., Schinn, D.S., Francke, T., Elsenbeer, H., Zimmermann, B., 2013.
10.1029/2011WR010508. Uncovering patterns of near-surface saturated hydraulic conductivity in an
Zhang, Y., Hernandez, M., Anson, E., Nearing, M.A., Wei, H., Stone, J.J., Heilman, P., overland ow-controlled landscape. Geoderma 195e196, 1e11. http://
2012. Modeling climate change effects on runoff and soil erosion in south- dx.doi.org/10.1016/j.geoderma.2012.11.002.
eastern Arizona rangelands and implications for mitigation with conservation