18th World IMACS / MODSIM Congress, Cairns, Australia 13-17 July 2009http://mssanz.org.au/modsim09

Hillslope and stream connectivity: simulation ofconcentration-discharge patterns using the HYDRUS
model
Smethurst, P.J.1, G. Langergraber2, K. C. Petrone3 and G.K. Holz4
1

Corresponding Author; CSIRO and Landscape Logic CERF, Private Bag 12, Hobart TAS 7001, Australia;
Email: Philip.Smethurst@csiro.au
2
Institute of Sanitary Engineering and Water Pollution Control, University of Natural Resources andApplied Life Sciences (BOKU), Muthgasse 18, A-1190 Vienna, Austria
3
CSIRO and Landscape Logic CERF, Private Bag 5, Wembley WA 6913, Australia
4
Antarctic Climate & Ecosystems CRC, Private Bag 80, Hobart TAS 7001, Australia

Abstract:
Nutrient concentrations and loads in streamflow are sensitive to rapidly changing stream
chemistry and discharge during storms. Mechanistic models that can simulate water and solute movement athillslope scales could be useful for predicting concentration-discharge (C-Q) patterns and thereby improve
our quantitative understanding of terrestrial-aquatic linkages for targeted catchment management. Our
objective was to use the HYDRUS model to represent hydro-biogeochemical processes in soils that drive
seepage of water and solutes from soil profiles into streams. Specifically we compared measurements in theliterature with HYDRUS outputs using two methods for simulating runoff. This model predicts runoff (R) as
rainfall that is instantaneously in excess of infiltration, but it is not designed to route runoff as overland flow.
Post-HYDRUS addition of seepage to runoff was used to simulate the delivery of dissolved or particulate
constituents to a stream (method A). Alternatively, we demonstrated how simulations using HYDRUS could
include a hypothetical layer at the top of the soil profile with extremely high porosity and hydraulicconductivity that enabled overland flow and down-slope infiltration, but in this case only dissolved
constituents could be considered (method B). These methods were evaluated by comparing the simulated
temporal patterns of discharge and concentration with observed patterns. The catchments considered were in
Slovenia (4210 ha) and in Australia (11.9 ha).
Methods A and B were shown to adequately simulate some aspects of published discharge-concentration
patterns, e.g. runoff dilution or concentration effects, but the temporal patterns of discharge for both methods
did not precisely match those measured at small time-steps (e.g. 15 minutes). This limitation was due mainly
to inadequate simulation of the down-slope movement of runoff and down-slope infiltration of a portion ofthis runoff. Method A was generally more useful than method B. Despite this limitation, both methods, if
used carefully, should be adequate for many purposes, especially when simulating longer time-steps.
Additional hypothetical simulations illustrated the significance of soil hydraulic conductivity, soil watercontent, and vertical gradients in solute concentrations in soil. Two temporal types of dischargeconcentration patterns were observed; short-term hysteresis caused by runoff during and shortly after arainfall event, and longer-term trends associated with infiltration and seepage. Clockwise and anti-clockwise
hysteresis was demonstrated to be potentially due to the temporal asynchrony of peak discharge and peakconcentration in runoff. Simulations also demonstrated advantages over using the more common approach ofa 2- or 3-component mixing model.
Our results suggest that the HYDRUS model will be useful for the mechanistic simulation of within-soil
processes that are needed to predict discharge-concentration patterns at hillslope scales.
Keywords: Water quality, nitrate, soil, runoff, overland flow, seepage, drainage

4057

Smethurst et al., Hillslope and stream connectivity: simulations using HYDRUS model1.

INTRODUCTION

Discharge-concentration patterns (C(Q); hysteresis) at a weir can be explained using 2- or 3-component
mixing models that rely on signatures of water chemistry from the various sources e.g. rainfall, soil water andground water, to indicate their proportional contribution to flow (Petrone et al., 2007). Evans and Davies(1998) used this approach to develop nine hypothetical hysteresis patterns from a particular pattern of soil
water, ground water (aquifer) and event water (rain) that each had constant but contrasting concentrations.
These patterns matched some event data from their own sites and others, but the patterns did not encompass
the complete range of hydrologic conditions. Weiler and McDonnell (2006) used hillslope modeling to
reproduce hysteresis patterns of nitrate (NO3) flushing associated with event water that completely infiltrated
highly conductive soil surface horizons. These authors also encouraged the use of mechanistic models in
virtual experiments for simulating hillslope water and solute transport and transformations, including
nitrogen mineralization and nitrification. Such models are needed to fully capture the potential effects of soil
and other characteristics that differ between sites.
Several models attempt to mechanistically represent water and solute transport processes in two or three
dimensions. The HYDRUS model (Šimůnek et al., 2008) is attractive because it incorporates conservativeand non-conservative solutes by using the Richards equation for saturated-unsaturated water flow and the
Fickian-based advective-dispersion equations for solute transport. It also includes options for dual
permeability, 3-dimensional transport, and ground water. HYDRUS also has an attractive and usefulgraphical user interface, and a high degree of spatial and temporal flexibility. We also had an interest in themechanistic simulation of nitrogen dynamics in the soil-plant system, and HYDRUS includes a module withthat potential. These attributes offered the potential that HYDRUS might be useful for mechanistic
simulations of water and solute transport at hillslope and headwater catchment scales, but its simulation of
runoff requires evaluation. Our objective here was to examine how HYDRUS N transport and hydrologic
processes in order to simulate observed discharge-concentration patterns. We compared simulated HYDRUS
patterns of discharge and concentration with those presented in the literature for contrasting rainfall eventsand sites,. We also ran hypothetical simulations to illustrate the importance of some within-soil processes for
producing various short- and long-term patterns.
2.

METHODSHYDRUS (Šimůnek et al. 2008, version 1.05) was used in a 2-dimensional, sloped, rectangular or trapezoidal
configuration. Units used were cm for length and mg/L for concentration. An atmospheric (precipitation)boundary condition was usually specified for the surface, with a vertical seepage face, and no-transfer
boundaries for other faces, unless constant flux and constant concentration groundwater input was simulated
for part of the lower surface. Seepage refers to water movement out of a soil profile at a seepage face with anatmospheric boundary condition (saturation excess), and can include components of interflow soon after
rainfall, stored soil water, and ground water entering the soil profile from an aquifer. We use the term runoff
to specifically mean overland flow in excess of infiltration. Some authors use the term deep seepage to imply
movement of water deep into a soil profile or into an unconfined aquifer. Such a process was not needed inour simulations, but could potentially be simulated in HYDRUS as a drainage boundary condition; we used a
no-flux lower boundary condition except where ground water was included.
Because we wanted to simulate hillslope processes in two dimensions, an average hillslope length was
calculated as catchment area divided by stream length. Total and horizon soil depths, hydraulicconductivities, and rainfall were taken from the literature (Table 1). At least 1,200 spatial nodes were used,
and time-steps started at very low values (10-7 s; to ensure adequate conservation of mass) and increased
during stable periods of a simulation to a maximum of 1-10 days. Simulations were built up by specifying
firstly water only, then by adding transpiration (root water uptake), and followed by solute transport. A low
evaporation rate was included. Before rainfall events were simulated, setting up of a simulation included preruns (up to 200 d), average rainfall and solute inputs that enabled a quasi steady-state to be achieved forseepage rate and concentration. Simulated seepage and runoff fluxes were in two-dimensional units (cm2/d)and converted to three-dimensional output (m3/d) by multiplying by the length of the third dimension(catchment length = catchment area/length of hillslope) and the cm to m conversion factor.

4058

Smethurst et al., Hillslope and stream connectivity: simulations using HYDRUS model
Two methods of simulating runoff were tested as follows. HYDRUS simulates runoff as rainfall that is
instantaneously in excess of infiltration (method A), but it is not designed to route runoff as overland flow.
This method required post-processing of HYDRUS runoff and seepage output using a spreadsheet to add
these fluxes (without runoff delay delay) for each printed time-step, and thereby estimate stream flow.
During post processing, concentrations in runoff needed to be specified in order to be added to seepage andthereby calculate concentrations in stream flow. As an alternative with potentially different outcomes, we
also tested the use in HYDRUS of a hypothetical layer at the top of the soil profile with extremely highporosity (0.98) and hydraulic conductivity (Ksat = 2E7 cm/d, method B), which facilitated simulation of
overland flow and down-slope infiltration. Only dissolved constituents could be considered with method B,
but less post-processing was required because all water and solute output as seepage by HYDRUS included
the runoff component.
Some characteristics of the catchments simulated are summarised in Table 1. Hypothetical variations werealso tested, including the hypothetical scenarios of Evans and Davies (1998).
Table 1. Some characteristics of the two catchments simulated as hillslopes.
Characteristic

Padez

MontaguoCountry, latitude

Slovenia ,45 North

Australia, 41o SouthArea (ha)

4210

11.9

Hillslope length (m)

515

13.2

Rainfall annual (mm/year)

1440

1283

Rainfall events (mm)34.5, 18.2

54.4

Rainfall peak intensity (mm/15 minutes)

1.2, 1.5

6.0

Land-use

forest 92%, pasture 18%

Pasture 100%

Geology

Marine sediments

Colluvium and alluvium

Soil depth (m)12Soil surface texture

silty clay

sandy loam over clay

Soil saturated 3D hydraulic conductivity (cm/d)31

60 over 6

Soil solid-liquid phase partition coefficient (mL/g)

00

References

Rusjan et al (2008)

Holz (2008, pers. comm.)

3.
3.1.

RESULTS AND DISCUSSION
Scenarios based on the Padež catchment

Rusjan et al. (2008) provided two quite contrasting measured hysteresis patterns. The first type was shown inboth their March and April events (Figs. 6a and 6b in Rusjan et al., 2008), in which 2- to 5-fold increases indischarge occurred with little change in concentration except a slight anticlockwise increase. These patterns
occurred with a 34.5 mm rainfall event on relatively dry soil at the end of winter. We simulated a similar
pattern (Fig. 1) using hillslope, soil profile, soil water, and rainfall as described by Rusjan et al (2008), andassuming initial concentrations of NO3 in soil water decreased from 2.5 mg/L at the soil surface to 0 mg/L at
the bottom of the profile. The simulated concentration was 1.4-1.5 mg/L, compared to 1.8-2.4 mg/L,
suggesting that the assumed initial condition for the concentration profile were not quite correct.
Changing to a high soil water content and uniform concentration of 5.8 mg/L NO3 down the profile, with an
18 mm rainfall event (as indicated by Rusjan et al 2008 for their August II-5 event), similar temporal patternsof discharge and concentration were produced to those observed, i.e. sharp increase and return in flowaccompanied by a concurrent decrease and return in concentration (Fig. 2). The resulting C(Q) plot showedno hysteresis while runoff occurred, and was followed by a steady but less substantial increase in discharge
while concentration remained constant at 5.8 mg/L, which reflected high NO3 concentrations observed duringlate spring.

4059

Smethurst et al., Hillslope and stream connectivity: simulations using HYDRUS model
We used this scenario to demonstrate how hypothetical patterns of concentration in runoff could produce
contrasting patterns of hysteresis, i.e. C(Q). By assuming the concentration in runoff peaked at 15 mg/L 1 h
before peak discharge, a clockwise hysteresis was simulated (Fig. 3 left). Conversely, by assuming
concentration in runoff peaked at 15 mg/L 0.5 h after peak discharge, an anticlockwise hysteresis was
simulated (Fig 3 right). Because concentrations in runoff are important for predicting C(Q) patterns, further
advances in these predictions will only be made if we develop a sufficient understanding of the mechanisms
contributing to concentrations of dissolved and particulate constituents in runoff and the resultant temporal
patterns of concentration in discharge. Some progress has been made, but generalised mechanistic predictions
are not yet available (Robertson and Nash, 2008; Vadas et al., 2008).
3

2
2000
1.5
1

1000

0.5
0

0
0

5

10

Concentration (mg/L)

2.5

C(S+R)

Concentration (mg/L)

3000

Discharge Rate (m3/d)

2.0

Q(S+R)

1.8
1.6
1.4
1.2
1.0

15

0

Time (d)

1000

2000

3000

Discharge Rate (m3/d)

Fig. 1. Simulated temporal patterns of discharge and concentration (left) and the resultant dischargeconcentration pattern (hysteresis; right), which are similar to those patterns described by Rusjan et al. (2008)
for their March and April rainfall events. Q = discharge, C = concentration, S+R = combined seepage andrunoff. Runoff method A was used, but insignificant runoff occurred during this event (observed andsimulated). The direction of hysteresis is indicated by arrows.
3000

7

7

6

2000

5

1500

4

1000
0.5

1

6

5

4

3

3
0

Concentration (mg/L)

C(S+R)

2500

Concentration (mg/L)

Discharge Rate (m3/d)

Q(S+R)

1000

1.5

1500

2000

2500

3000

Discharge Rate (m3/d)

Time (d)Fig. 2. Patterns of C(Q) simulated for the Rusjan et al. (2008) August II-5 rainfall event in Slovenia.
Significant runoff occurred during this event (observed and simulated). Simulated flow (left) shows theseparate influence on flow of runoff that lasted a few hours and longer-term seepage at an increased rate due
to infiltration. Runoff method A was used.
We also tested the use of an artificial runoff layer (method B) for the scenario presented in Fig. 2, but
because most ‘runoff’ water reinfiltrated the profile, a very different C(Q) pattern was predicted to thatobserved, i.e. in this case it simulated a slow increase in discharge and only a minor change in concentration(data not shown).

4060

Smethurst et al., Hillslope and stream connectivity: simulations using HYDRUS model

8

Concentration (mg/L)

Concentration (mg/L)

8
7
6
5
4
3
1000

1500

2000

7
6
5
4
3
1000

2500

1500

2000

2500

Discharge Rate (m3/d)

Discharge Rate (m3/d)

Fig. 3. Patterns of C(Q) simulated for the Rusjan et al (2008) August II-5 rainfall event assuming peakconcentrations preceded (left) or followed (right) peak discharge. The direction of hysteresis is indicated byarrows. Runoff method A was used.

3.2.

Inclusion of high-resolution, measured rainfall for predicting NO3 and PO4 patterns

We simulated a transect of one half of a hump-and-hollow pasture system in the Montagu catchment
described by Holz (2008 and pers. comm.). This pasture system is used in flat landscapes by moving soil
from the hollows and depositing to construct the humps, thereby creating mini-catchments. Average annual
rainfall was used to initialise water flow conditions, and solute concentrations in the profile were set atconcentrations measured during baseflow. A particular period was selected for simulation because it included
a complex pattern of rainfall events during a period of a few days, with an associated complex pattern ofdischarge and NO3 and PO4 concentrations. Hence, this period was chosen to test the behaviour of HYDRUS
under complex rainfall patterns. The two methods of simulating runoff evaluated earlier were again tested
here.
Using Method A, Ksat was varied within the measured range to achieve peak flow rates approximately 8.5times higher than baseflow, i.e. similar to the measured ratio of peak flow to base flow. Using these
conditions, the sharp rise and fall in runoff was matched by a decrease in NO3 concentrations during peakrunoff that was similar to those measured, which represented the expected dilution effect (Fig. 4 left).
Conversely, assuming a pattern of PO4 concentration variation with runoff that peaked at a similar value tothat measured, a PO4 concentration effect was simulated during this event (Fig. 4 right). However, simulated
runoff did not occur during other rainfall events during this 3.4 day period, in contrast to that suggested bymeasured patterns of concentrations. This difference between observed and simulated behaviour might reflect
the delays in surface drainage in this low-slope condition, which we did not simulate.
Method B provided two runoff events (not shown). However, using this method, the ratio of peak discharge
to baseflow was much higher than measured and discharge returned to baseflow conditions between thesetwo peaks, which did not match the measured pattern. Again, this difference between observed and predicteddischarge patterns might be a reflection of differences in storage and the delay in drainage of overland flowin this low-slope catchment. Nitrate dilution was simulated using this method, but PO4 patterns were notattempted. This result from the Montagu catchment and those above from the Padež catchment do not render
method B unsuitable for all applications, but indicates that caution at least will be needed when it is used.
3.3.

Hypothetical Scenarios That Included Groundwater

Various boundary conditions are provided as options in HYDRUS (e.g. constant or variable fluxes or heads),
which allows scenarios to include ground water, i.e. water entering the catchment from outside its boundaries
via an aquifer. Evans and Davies (1998) included ground water, soil water and runoff in a hypothetical 3component mixing model to develop nine patterns of hysteresis using one particular temporal pattern offlows for these three components. Using runoff method A, we reproduced their type A3 pattern using
HYDRUS, i.e. anticlockwise, concave and negative for CG>CSO>CSE (Fig. 5; see figure for symbol
definitions), which indicates that groundwater was adequately simulated by HYDRUS, but attempts to
reproduce another pattern (C1) were unsuccessful because, in a hill-slope scenario, we could not exactly

4061Smethurst et al., Hillslope and stream connectivity: simulations using HYDRUS model
reproduce the temporal patterns of discharge and concentration for the three water sources used in the
hypothetical mixing model. Evans and Davies (1998) noted also that their patterns did not completely
describe their own observed system nor those of some other authors. Hence, we abandoned attempts to
reproduce the other seven patterns of these authors.

Q(S+R)

600

Runoff
400

C(S+R)

0.5

200
0Discharge Rate (m3/d)

1.0

Concentration (mg/L)

Discharge Rate (m3/d)

800

1000

1.25

1.3

1.35

Q(S+R)

800

Runoff
600

1.4

1.0

C(S+R)

400

0.5

200
0
1.20

0.0
1.2

1.5

Time (d)

1.25

1.30

1.35

Concentration (mg/L)

1.5

1000

0.0
1.40

Time (d)Fig. 4. Simulated temporal patterns of discharge and concentrations of NO3 (left) and phosphate (right)
during a complex rainfall event in the Montagu catchment that led to runoff.
600

10

8

8

400

6

300

200

4

0.1

1

10

100

1000

Concentration (mg/L)

C

Concentration (mg/L)

Realtive Discharge Rate

Q
500

7

6

5

4
250

350

450

550

Relative Discharge Rate

Time (d)Fig. 5. A rainfall runoff event simulated with three water sources with contrasting concentrations. This scenario isakin to type A3 as described Evans and Davies (1998), i.e. ground-water concentration (CG) > soil water
concentration (CSW)> event concentration (CSE). The direction of hysteresis is indicated by arrows.

3.4.

Difficulties Encountered

Successful HYDRUS simulations require a sound knowledge of soil processes in combination with aworking knowledge of numerical models generally, i.e. use of time steps and spatial nodes, and initial andboundary conditions. Failure of HYDRUS to reach a solution was not uncommon during the development of
this paper, but sometimes failure was due to the user expecting the model to solve unrealistic situations, e.g.
when extreme temporal or spatial gradients in pressure or concentration were unintentionally specified. Trialand error was required to develop solutions under many conditions.
HYDRUS requires boundary conditions to be set prior to a simulation, e.g. an atmospheric boundary
condition is required for rainfall. Although HYDRUS simulates infiltration-excess runoff, simulation ofsaturation-excess runoff required progressive changing of surface nodes from an atmospheric boundary
condition to a seepage face, which greatly complicated maintenance of the hydrological balance, because
converted nodes no longer receive rainfall. Instead, we tested the approach of adding a runoff layer to thelarge catchment in the Rusjan et al. (2008) August II-5 scenario presented in Fig. 2. However, the simulation
did not match observations or the simulation in Fig. 2 due to excessive down-slope infiltration. Althoughsome rapid runoff was observed, it was minor and resulted in only a small dilution of seepage water. In
contrast, when an artificial runoff layer was used to simulate the very small catchment of the hump-and-4062

Smethurst et al., Hillslope and stream connectivity: simulations using HYDRUS model
hollow system in Fig. 4, excessively rapid runoff occurred. Although use of an artificial runoff layer might bewarranted under some conditions (perhaps longer time-scales and with appropriate tuning of surface soil andrunoff layer Ksat values), results of these tests were not encouraging.
4.

CONCLUSIONS

Using simple assumptions and measured values of soil, rainfall and solute characteristics, the HYDRUS
model simulated some generalised aspects of C(Q) patterns, i.e. a large increase in discharge resulting fromrainfall exceeding infiltration, and accompanying dilution or enrichment effects on solute concentrations. It
was useful for evaluating hypotheses about concentration-discharge patterns that offers a means of improving
upon those resulting from mixing model interpretations, because HYDRUS can mechanistically integratemany of the landscape processes that affect these patterns. Overall, HYDRUS was not easy to use in ahillslope context that included runoff. Two methods of simulating runoff and discharge (runoff plus seepage)were evaluated. Both methods had advantages and disadvantages, but method A was more useful for theapplications tested. Applications of HYDRUS to the prediction of stream flow and water quality will bestrongest where temporal patterns are dominated by processes within the soil profile rather than thoseaffecting temporal patterns in the amount or quality of overland flow. HYDRUS simulations should not be
expected to reflect high-resolution patterns of runoff or solute concentrations, but generalised, simulated
longer-term patterns will probably be useful for testing some hypotheses about hillslope processes affecting
stream flow and water quality. Efforts to find alternative approaches should continue.
ACKNOWLEDGEMENTS
Advice and training from Jirka Šimůnek was much appreciated. This research was funded by CSIRO and theLandscape Logic CERF hub.
REFERENCES
Evans, C., and T.D. Davies, (1998) Causes of concentration/discharge hysteresis and its potential as a tool for
analysis of episode hydrochemistry. Water Resources Research 34, 129-137.
Holz, G.K. (2008) Seasonal variation in groundwater levels and quality under intensively drained and grazed
pastures in the Montagu catchment, NW Tasmania. Agricultural Water Management 96, 255-266.
Petrone, K.C., L.D. Hinzman, H. Shibata, J.B. Jones, and R.D. Boone (2007) The influence of fire andpermafrost on sub-arctic stream chemistry during storms. Hydrological Processes 21, 423-434.
Robertson, F.A., and D.M. Nash (2008) Phosphorus and nitrogen in soil, plants, and overland flow fromsheep-grazed pastures fertilized with different rates of superphosphate. Agriculture, Ecosystems andEnvironment 126, 195-208.
Rusjan, S., M. Brilly, and M. Mikos (2008) Flushing of nitrate from a forested watershed: an insight intohydrological nitrate mobilization mechanisms through seasonal high-frequency stream nitrate dynamics.
Journal of Hydrology 354, 187-202.
Šimůnek, J., M.Th. van Genuchten, and M. Šenja (2008) Development and application of the HYDRUS andSTANDMOD software packages and realted codes. Vadose Zone 7:587-600.
Vadas, P.A., L.B. Owens, and A.N. Sharpley (2008) An empirical model for dissolved phosphorus in runoff
from surface-applied fertilizers. Agriculture, Ecosystems and Environment 127, 59-65.
Weiler, M., and J.J. McDonnell (2006) Testing nutrient flushing hypotheses at the hillslope scale: A virtual
experiment approach. Journal of Hydrology 319, 339-356.

4063