Author's personal copy - Sevilleta LTER - LTER Network
This article appeared in a journal published by Elsevier. The attached copy is furnished to the author for internal non-commercial research and education use, including for instruction at the authors institution and sharing with colleagues. Other uses, including reproduction and distribution, or selling or licensing copies, or posting to personal, institutional or third party websites are prohibited. In most cases authors are permitted to post their version of the article (e.g. in Word or Tex form) to their personal website or institutional repository. Authors requiring further information regarding Elsevier’s archiving and manuscript policies are encouraged to visit: http://www.elsevier.com/copyright
Author's personal copy
Journal of Hydrology 388 (2010) 258–272
Contents lists available at ScienceDirect
Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol
Hydrology, erosion and nutrient transfers over a transition from semi-arid grassland to shrubland in the South-Western USA: A modelling assessment Laura Turnbull a,*, John Wainwright b, Richard E. Brazier c a
Global Institute of Sustainability, Arizona State University, 800 South Cady Mall, Tempe, AZ 85287-5402, USA Sheffield Centre for International Drylands Research, Department of Geography, University of Sheffield, Winter Street, Sheffield S10 2TN, UK c Department of Geography, University of Exeter, Amory Building, Rennes Drive, Exeter EX4 4RJ, UK b
a r t i c l e
i n f o
Article history: Received 26 August 2009 Received in revised form 27 April 2010 Accepted 3 May 2010
This manuscript was handled by L. Charlet, Editor-in-Chief, with the assistance of Bernhard Wehrli, Associate Editor Keywords: Runoff Erosion Land degradation Nitrogen Phosphorus Ecogeomorphology
s u m m a r y Land degradation in arid and semi-arid areas, as a consequence of the invasion of grasslands by shrubs, is often associated with an increase in runoff and erosion and a change in nutrient transport. Modelling of nutrient transport during runoff events (in particular particulate-bound nutrients), is especially important, since the spatial redistribution of nutrients (in addition to water and sediment) can have significant implications for vegetation dynamics in these ecosystems. In this study, Mahleran (Model for Assessing Hillslope to Landscape Erosion Runoff, And Nutrients) is extensively evaluated against runoff and erosion data from four plots (representative of different stages of land degradation) over a transition from grassland to shrubland, at the Sevilleta National Wildlife Refuge in New Mexico, USA. A new particulate-bound nutrient module was developed to include a representation of particulate-bound nutrient dynamics, which is an important form of nutrient transport in these ecosystems. Understanding dynamics of both dissolved and particulate-bound nutrient dynamics during runoff events is imperative, because of their differing roles in terms of nutrient bioavailability and potential implications for plant dynamics. Results of the model evaluation show that the runoff and erosion components of Mahleran perform reasonably well, as does the new particulate-bound nutrient sub-model, though not consistently. Performance of the particulate-bound nutrient model was better for the end-member plots, because of better parameterization data available for end-member vegetation types. Since the particulate-bound nutrient sub-model is by necessity strongly dependent on the simulated erosion rate, the performance of the particulate-bound nutrient model is dependent on the performance of the erosion component of Mahleran, so that when erosion is well represented by the model, so typically are particulate nutrient transfers. The performance of the dissolved nutrient component of Mahleran was poor in this application, which indicates that the process representation for this semi-arid environment and the parameterisation of the dissolved nutrient component were inadequate. Results from the model evaluation suggest that an improved understanding of dissolved nutrient dynamics during runoff events and simulation if interevent nutrient dynamics is required, in order to improve the level of process representation within modelling approaches and thus the ability to simulate dissolved nutrient dynamics and their subsequent effects on other ecosystem processes. Ó 2010 Elsevier B.V. All rights reserved.
1. Introduction Recent research on changes in runoff and erosion over a grassland to shrubland transition (Turnbull et al., 2010, in preparation-a,-b) has provided new insights into changes in ecosystem structure and function. The observed changes in ecosystem structure and function over a semi-arid, grass–shrub ecotone demonstrate that changes in ecosystem structure and function do not necessarily follow a linear continuum across a grassland to shrubland transition. As the distribution and percentage cover of * Corresponding author. E-mail address: [email protected] (L. Turnbull). 0022-1694/$ - see front matter Ó 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.jhydrol.2010.05.005
vegetation decreased over the grassland to shrubland transition, concurrent changes in the connectivity of these surface properties have profound effects on hydrological function (Turnbull et al., 2010). Determining more fully the interactions between ecosystem structure and function, and interactions/feedbacks between biotic and abiotic components of the ecosystem is imperative in order to ascertain the effects that climate and land-use change will have on these ecosystems in the future (Turnbull et al., 2008). Modelling-based studies provide a potential means of expanding further our understanding of interactions and feedbacks between different components of an ecosystem, under a multitude of climatic and land-use conditions, that would not be feasible to investigate using field-based experimentation alone. However,
Author's personal copy
L. Turnbull et al. / Journal of Hydrology 388 (2010) 258–272
the limitation at present of carrying out such modelling-based studies is that the transfers of matter and energy among units are rarely included in models of vegetation change, but are critically important in arid ecosystems in which resource redistribution plays a significant role (Peters and Herrick, 2001). Therefore, the development of existing modelling approaches is imperative in order to understand further the dynamics of land degradation in semi-arid areas. The aim of this study is to develop a process-based model to understand better the dynamics of dissolved and particulate-bound nutrient transport processes in runoff, with particular emphasis on the testing of the model algorithms. The specific model to be developed is Mahleran (Model for Assessing Hillslope to Landscape Erosion, Runoff And Nutrients; Müller et al., 2007a,b; Wainwright et al., 2008a,b). 2. Modelling approach In the following sections the key components of Mahleran are outlined: (i) the hydrology component (Wainwright and Parsons, 2002), (ii) the erosion component (Wainwright et al., 2008b), (iii) the dissolved nutrient component (Müller et al., 2007a,b), and (iv) a new particulate-bound nutrient component of Mahleran which is developed here, based on the improved understanding of particulate-bound nutrient fluxes outlined in Turnbull et al. (in preparation-b). 2.1. Hydrology The hydrological component of Mahleran (Wainwright and Parsons, 2002) uses a simple infiltration model to generate infiltration- and saturation-excess runoff. Infiltration-excess runoff is generated by determining the difference between the rainfall rate and the infiltration rate. The infiltration rate is simulated using the Smith-Parlange approach (Smith and Parlange, 1978). Runoff is routed over the hillslope using a kinematic wave approximation of the St Venant equations (Wainwright and Parsons, 2002), with flow routing in the direction of steepest descent from cell to cell (in cardinal directions) over a finite-difference grid (Scoging et al., 1992). The flow is routed over the grid using a finite-difference solution (Euler backward difference form: Scoging et al., 1992). Flow velocity is determined dynamically using the DarcyWeisbach flow equation. This modelling approach differs from most kinematic wave approaches, which assume a fixed rating function, and thus poorly represent feedbacks of flow hydraulics with runoff and infiltration. Full details of the hydrology and hydraulics components are found in Wainwright and Parsons (2002); Parsons et al. (1997), and Wainwright et al., 2008b. 2.2. Erosion Mahleran considers the interaction of both concentrated and unconcentrated erosion as bedload and in suspension, which is particularly important since the relative balance of the different processes is the most critical control on the resulting pattern of erosion (Wainwright et al., 2008a). Mahleran employs an approach to erosion modelling that is based upon the concept of entrainment and travel distances of particles from six different size classes (u = 1 is 12 mm). Sediment transport is modelled by a form of the continuity equation, in which sediment belonging to these different size classes are simulated separately (Wainwright et al., 2008a).
@hs;u @qs;u þ eu þ du ¼ 0 @t @x
ð1Þ
259
where hs is the equivalent depth of sediment transport (m), qs is the unit discharge of sediment (m2 s-1), e is the rate of erosion of the surface (m s1), d is the rate of deposition (m s1), u is an index relating to a specific size class of sediment, t is time (s) and x is location in space (m). To account for the temporal scaling characteristics of this model, qs,u is defined as
qs;u ¼ hs;u v p;u
ð2Þ
where vp,u is the virtual velocity of transport of sediment of particlesize class u (m s1), which can be derived from empirical observations of sediment transport distance per unit time or by theoretical calculation by differentiation of transport distance per unit of time (Wainwright et al., 2008a). Mahleran simulates three sediment detachment and transport conditions (see Wainwright et al., 2008a, for more detail) (i) erosion as a function of raindrop detachment and transport occurs by splash when no flow is present; (ii) unconcentrated overland flow erosion simulated using raindrop detachment rates that are modified to account for the protective effects of the surface-water layer; (iii) transport by mixed unconcentrated and concentrated flow (500 6 Re 6 2000) (where Re is the flow Reynolds number); or (iv) transport by concentrated flow (Re > 2000). Deposition is modelled similarly, using the transport-distance approach, since the distribution function of travel differences of particles under specific transport mechanisms and flow conditions enables the deposition rate to be determined directly at each point along the transport pathway (Wainwright et al., 2008a). Rates of deposition depend on whether material being moved is in splash, as bedload in unconcentrated or concentrated overland flows or in suspension (Wainwright et al., 2008a).
2.3. Nutrients Dissolved nutrient concentrations in runoff are affected by the nutrient content in the upper soil layers, which depends upon the nutrient content of the soil and the soil-exchange characteristics (Barger et al., 2006). A common approach employed to model the transfer of solutes from soil solution to runoff is based on the principle that nutrients in the soil solution will diffuse into runoff water. The simplest assumption in determining the mass transfer of nutrients between the soil surface and overland flow is that the rate of exchange is directly proportional to the difference in concentration between the soil solution and the surface water (Walton et al., 2000b). Wallach and van Genuchten (1990) developed a modelling approach in which the diffusion of solutes from the soil solution to runoff water was coupled with a rate-limited chemical transfer process across a laminar boundary at the soil surface/runoff water interface. Solute release into runoff is limited by the initial mass of nutrient near the soil surface and the depth at which surface water interacts with the soil (Havis et al., 1992). The depth of soil that interacts with overland flow through turbulent mixing, convection and diffusion, that contributes chemicals to the surface water during a rainfall event, is defined as the mixing zone (Havis et al., 1992). Solute exchange through this boundary layer was modelled using a mass-transfer coefficient, K (Wallach and van Genuchten, 1990). Walton et al. (2000b) modelled the mass transfer of nutrients between the soil surface and overland flow, assuming that the rate of change of mass transfer is directly proportional to the difference between the concentration in the soil solution and the concentration in the surface water.
Author's personal copy
260
L. Turnbull et al. / Journal of Hydrology 388 (2010) 258–272
Müller et al. (2007a,b) developed the dissolved nutrient-transport component of Mahleran, in which NO3, NH4, and PO4 are modelled by an advection–dispersion model (Havis et al., 1992; Walton et al., 2000b). In the advection–dispersion model the mass transfer of nutrients from the soil surface to overland flow is driven by: (i) diffusion of dissolved nutrients from the soil interstices by movement of soil water into the overland flow, (ii) desorption of the nutrients from soil particles into the overland flow, (iii) dissolution of solid phase nutrients into the soil water or overland flow, and (iv) scouring of solid phase nutrients by hydraulic forces and subsequent transport and moving dissolution (Müller et al., 2007a). An advection–dispersion model is used to model nutrient transport, in which a mass-transfer coefficient is used to lump together the aforementioned mechanisms of mass transfer (Wallach and van Genuchten, 1990). The nutrient mass balance is given by.
@ðChÞ @ðqCÞ þ ¼ KðC s CÞ fC @t @x
ð3Þ
where C is the nutrient concentration in overland flow (mg l1), Cs is the nutrient concentration in the soil-solution zone (mg l1), h is the depth of overland flow (mm), q is the unit discharge (mm2 s1) and f is the infiltration rate (mm h1) (Müller et al., 2007a). However, in this modelling approach, the effect of the nutrient content of the incoming rainfall is neglected. Turnbull et al. (in preparation-b) demonstrate that the concentrations of dissolved nutrients in rainfall are often of the same order of magnitude as those measured in runoff. Therefore, when modelling diffusion of nutrients from the soil-solution zone into runoff, the consideration of the concentration of nutrients in rainfall is important. Otherwise, it is possible that net losses of nutrients from the soil-solution zone will be over-estimated. To account for the effects of rainfall-nutrient concentration on the diffusion of nutrients from the soil-solution zone into runoff, Eq. (3) can be re-written as
@ðChÞ @ðqCÞ þ ¼ KðC s CÞ fC þ RC r @t @x
ð4Þ
where R is the rainfall rate (mm h1) and Cr is the concentration of the nutrient in the rainfall (mg l1) (Havis et al., 1992; Walton et al., 2000a). Havis et al. (1992) related the potential mass (M) of a nutrient over an area (A).Cs is related to the nutrient concentration of the soil by
M ¼nA
Z
x
C s ðzÞdz
ð5Þ
0
where M is the mass of available nutrient in the soil (g), A is the model cell size (mm2), n is porosity (%), x is the depth of the mixing zone (mm) and z is the depth through the soil (mm). Once soluble nutrients are entrained by runoff, they are routed conservatively (e.g. Viney et al., 2000). Although data from field studies, for example data presented in Turnbull et al. (in preparation-b) and Lister (2007) indicate that NO3–N, NH4–N and PO4–P in runoff are not conservative, at present there is insufficient information to represent this non-conservative behaviour in the modelling approach, which would possibly induce greater errors than employing the simple conservative approach. The extent to which this assumption may be the case will be evaluated in the context of the results of this modelling exercise Turnbull et al. (in preparation-b) demonstrate the significance of particulate-bound nutrient transport, which must thus be incorporated into the modelling approach. Very few semi-arid, field-based studies have investigated particulate-bound nutrient transfers in runoff. Since the significance of particulate-bound nutrient transfers in runoff has not been well recognised, very few models have represented this particulate phase of nutrient transfer. A key problem with employing detailed, process-based approaches to the modelling of nutrient transfers in semi-arid
environments is that understanding of the processes of nutrient transfers between states remains limited. In addition, the use of a modelling approach that has high parameterization demands is impractical, and may induce increased uncertainty into the processes being modelled (where data are not sufficient to constrain the additional uncertainty introduced by extra parameters, or more complex model structures – see Brazier et al., 2000).Fierer and Gabet (2002) calculated sediment-bound nutrient loss as the product of sediment detachment rate and the percentage C and percentage N associated with the sediment. They determined the percentage of C and N associated with the sediment by carrying out rainfall-simulation experiments at their semi-arid field site in California, where they sampled the sediment transported by runoff, which was then weighed to determine the sediment concentration in runoff and analysed for C and N content on a Fisons NA1500 C/N analyzer. The approach used to model particulate-bound nutrient transport in Mahleran elaborates upon that of Fierer and Gabet (2002) by determining the particulate-bound nutrient transport as a function of the nutrient concentration of sediment, s, in each particle-size class and the amount of sediment transported in each particle-size class. The flux of a nutrient of type n, adsorbed onto sediment in each particle-size class per time-step is related to the discharge of sediment (qs,u) and the nutrient concentration of sediment in each particle-size class (Cn,u) by:
qn;u ¼ qs;u C n;u
ð6Þ
The modelling of particulate-bound nutrient fluxes is conservative, since it is assumed that there is no cycling of adsorption/ desorption of particulate-bound nutrients when they are being transported (Viney et al., 2000). 3. Parameterizing Mahleran for Grassland–Shrubland Plots at the Sevilleta LTER Site and methods of model testing In order to evaluate the different model components in the context of changing vegetation structures, Mahleran is parameterized for four sites over a transition from black grama (Bouteloua eriopoda) grassland to creosotebush (Larrea tridentata) shrubland at the Sevilleta Long-Term Ecological Research site in central New Mexico, USA, which is lies on the northern margin of the Chihuahuan desert. Sites were located within approximately 1 km of each other and have similar slopes and aspects. The Sevilleta receives an average of 256 mm of rainfall annually, of which 55% typically falls during the summer monsoon rainfall season (July–September) (Dahm and Moore, 1994). At each site, a runoff plot was constructed measuring 10 30 m, from which flow was measured at 1-min intervals. The outflow from each plot was instrumented with an ISCO 6700 autosampler which sampled flow leaving plots at 1-min intervals. In addition to this instrumentation, a collection tank at the base of each plot trapped all runoff and sediment exported from the plots, to provide a secondary measure of total plot runoff and a measure of the total amount of sediment eroded from each plot. A more detailed description of the field site is provided in Turnbull et al. (2010, in preparation-a,-b) and Sevilleta LTER (2009). 3.1. Hydrology Parameterizing final infiltration rate (Ks) is notoriously difficult since Ks is variable over space and time (Müller et al., 2008). Attempts have been made to determine Ks in the field, using infiltrometers, but this technique has been found to overestimate Ks (Müller, 2007). An alternative and preferable method of determining Ks is to carry out rainfall-simulation experiments. However, this approach is not straightforward and was not possible within
Author's personal copy
261
L. Turnbull et al. / Journal of Hydrology 388 (2010) 258–272
this study. Thus, Ks could not be parameterized directly at the field site and therefore was estimated by using the best available information for similar sites where desert pavement progressively occupies the surface as vegetation changes from grass to shrubs. Previous studies have shown that well-developed desert pavements exert a dominant control on infiltration rates (Abrahams et al., 1989). Following Wainwright et al. (2008a) the rainfall-simulation experiments of Parsons et al. (1996) were used to parameterize Ks (mm min1) by
K s ¼ 0:351 þ 0:010 rain 0:006P%
ð7Þ
where rain is rainfall intensity (mm h1) and P% is pavement cover (%). There are limitations inherent in using this relationship to determine Ks for the field site in question, since this relationship was derived for a different field site, albeit one with similar soils. However, the benefit of employing this relationship to determine Ks is that it accounts for the effect of varying stone-pavement cover and rainfall intensities on Ks. Although Ks is likely to vary over the grass to shrub transition, appropriate data are not available to constrain the likely differences, and therefore the above equation is used to parameterize Ks as a first approximation. The quality of this approximation is evaluated below Numerous approaches have been described previously for the determination of the friction factor, ff (see Wainwright et al., 2008b), which differ in their levels of complexity and parameterization requirements. The Darcy-Weisbach friction factor (dimensionless) can be parameterized from field experiments (e.g. Parsons et al., 1994), although these are very time consuming and labour intensive. As a first order approximation, the simple depth–feedback parameterization of Scoging et al. (1992) was used
ff ¼ 14 0:8h
ð8Þ
where h is the depth of runoff (mm). The depth to wetting front must be known a priori, and is a difficult parameter to measure (Wainwright et al., 2008c). Soil thickness represents effective depth to wetting front, which was thus estimated to be 50 cm, based on measurements of soil depth. The wetting front suction parameter was estimated to be 5 cm. The saturated soil-moisture content was set to 39%. Again, while saturated soil-moisture content is likely to be spatially variable, no data were collected to constrain this spatial variability, thus a uniform value was used. The kinetic energy (KE in J m2 mm1) of the rainfall arriving at the ground surface is related to vegetation characteristics. A function for the relationship between KE and vegetation has been derived from Wainwright et al. (1999)
KEv ¼ KEð1 8:1 103 V%Þ
ð9Þ
where KEv (J m2 mm1) is the KE of rainfall arriving at the ground surface that is adjusted to account for vegetation cover and V% is the vegetation canopy cover (%). This function was derived from experiments carried out for creosotebush. No data has been collected to determine the relationship between KE and black grama canopy cover. Since this function was derived for creosotebush, it accounts for effect of leaf drip energy, which will be negligible for grass-covered areas. However, during preliminary testing of the model the sensitivity of erosion to leaf drip energy over grass-covered areas was found to be very small. Therefore, the function for the relationship between KE and vegetation derived from Wainwright et al. (1999) is employed for all vegetation types. Analysis of the spatio-temporal volumetric soil-moisture content dataset (Turnbull et al., 2010) showed that there are often significant differences in soil-moisture content between bare and vegetated soils. In view of these findings, Mahleran was developed to improve the process representation of soil-moisture dynamics, to enable spatially distributed parameterization of initial soil-
moisture content. Thus, rather than a spatially uniform parameterization of soil-moisture content, as had been carried out previously (Wainwright et al., 2008b), in this modified version of Mahleran, soil-moisture content at the start of the event is parameterized using a spatio-temporal dataset of volumetric soil-moisture content for different surface-cover types (Turnbull et al., 2010). 3.2. Erosion Parameterization of the soil particle-size distribution in Mahleran is made spatially distributed, based on the stone-pavement cover (made up of u = 5and u = 6), which was mapped digitally for the field site from aerial photographs (Turnbull et al., in preparation-a). The proportion of sediment in the remaining particle-size classes is also spatially distributed, based on the particle-size characteristics of the soil measured over each plot (see Turnbull et al., in preparation-a). The parameterization of sediment detachment and transport parameters for each u class were taken from Wainwright et al. (2008b). 3.3. Nutrients In order to parameterize the nutrient content of sediment from each particle-size class, the data of Lister (2007) are used, obtained from rainfall-simulation experiments carried out over 1–m2 plots over grass and shrub end-members at the Jornada Experimental Range, New Mexico. Therefore, where grass formed the dominant vegetation cover (plots 1 and 2), the grass end-member parameterization was used, and where shrubs formed the dominant vegetation cover (plots 3 and 4), the shrub end-member parameterization was used. The nutrient content associated with each particle-size class of eroded sediment was determined. The flux (g min1) of particulate-bound nutrients is modelled here as a function of the total amount of sediment flux, the size composition of that sediment and the nutrient content of sediment in each particle-size class. To parameterize the mass-transfer coefficient (K), the values used by Müller (2007) were used (Table 1), which were found to perform well for a similar setting (Müller et al., 2007b). The nutrient concentration in the soil-solution zone (Cs) was determined using Eq. (5). The resultant values used to parameterize Cs for each site for each dissolved nutrient are shown in Table 1. 3.4. Analysis techniques used to test model performance The Nash and Sutcliffe (NS) coefficient of efficiency (Nash and Sutcliffe, 1970) is used to measure the performance of the model compared to the observed runoff, erosion and nutrient fluxes. Rigorous quality control of observed fluxes was carried out (see Turnbull et al., 2010) in order to ensure that model evaluations were carried out against data in which there was a high degree of confidence. The performance of the hydrology model is evaluated against the goodness of fit of the modelled hydrograph with the monitored hydrograph, and against the total volume of modelled runoff with that monitored. Nash–Sutcliffe model efficiencies were
Table 1 Parameterisation of Cs and K for NH4–N, NO3–N and PO4–P. The parameterisation of K is the same for all plots. Plot 1 is the grass end-member, and plot 4 the shrub endmember.
Cs (mg l1)
K (mm s1)
Plot Plot Plot Plot
1 2 3 4
NH4–N
NO3–N
PO4–P
11 12 9 11
11 15 10 11
61 61 61 61
0.00003
0.00007
0.00004
Author's personal copy
262
L. Turnbull et al. / Journal of Hydrology 388 (2010) 258–272
also calculated for the upper and lower uncertainty margins of the total runoff monitored from the plots (see Turnbull et al., 2010) in order that uncertainties in the monitored data are taken into account in the model evaluation. To test the spatial predictions of flow versus no flow conditions by the model, for events over each plot, contingency tables were used. In order to approximate the modelled ‘no flow’ conditions, it is assumed that a total event discharge of 0.005 m3 or less, over a cell of 0.5 m2 is analogous to no flow conditions, since total event
discharges below this amount will not have been sampled by minisamplers (Parsons et al., 2003) which were designed to monitor flow depth at fifteen points over each plot. Kappa coefficients were calculated to evaluate the agreement of monitored and modelled results against those that might be expected by chance (Landis and Koch, 1977; Monserud and Leemans, 1992). Values of j < 0 indicate a poor strength of agreement, 0–0.2 indicates a slight chance of agreement, 0.21–0.4 indicates a fair chance of agreement, 0.41–0.6 indicates a moderate chance of agreement,
Table 2 Summary of model simulation results for events monitored over plots 1–4. The Nash–Sutcliffe efficiency could not be calculated for events with an asterisk in the Nash–Sutcliffe column because equipment failures prevented fitting against the hydrograph. Nash–Sutcliffe model efficiencies were also calculated for the upper and lower bounds of uncertainty of monitored runoff. Spearman’s rank-correlation coefficients of monitored versus modelled runoff and Nash–Sutcliffe efficiency statistic for event totals of monitored and modelled runoff are given. Nash–Sutcliffe efficiency statistics are given for upper and lower uncertainty bounds of monitored data that were derived from analysis of hydrograph uncertainty. Plot/event Plot 1 07/09/05 29/09/05 05/07/06 28/07/06 01/08/06 11/08/06 29/08/06 07/09/06 All events Upper uncertainty Lower uncertainty Plot 2 07/09/05 31/07/06 01/08/06 04/08/06 11/08/06 15/08/06 29/08/06 07/09/06 All events Upper uncertainty Lower uncertainty Plot 3 28/09/05 29/09/05 05/07/06 28/07/06 31/07/06 01/08/06 11/08/06 29/08/06 07/09/06 All events Upper uncertainty Lower uncertainty Plot 4 07/09/05 28/09/05 29/09/05 05/07/06 28/07/06a 28/07/06b 29/07/06 31/07/06 01/08/06 04/08/06 11/08/06 23/08/06 29/08/06 All events Upper uncertainty Lower uncertainty
Total runoff (l) Monitored
Total runoff (l) Modelled
5230 (+2250, 1950) 124 (+18, 18) 277 (+18, 18) 177 (+18, 18) 1011 (+18, 18) 213 (+18, 18) 2820 (+1249, 639) 1426 (+18, 18)
9848 1590 372 572 1110 445 1818 1586
2294 (+540, 261) 83 (+17, 17) 464 (+17, 17) 38 (+16, 16) 149 (+17, 17) 391 (+17, 17) 995 (+17, 17) 1260 (+17, 17)
363 (+17, 17) 122 (+17, 17) 66 (+17, 17) 380 (+17, 17) 1799 (+17, 17) 2030 (+881, 0) 637 (+17, 17) 551 (+17, 17) 6550 (+766, 1623)
2174 (+635, 0.0) 142 (+18, 18) 443 (+18, 18) 2174 (+785, 0.0) 1504 (+907, 468) 452 (+378, 169) 199 (+237, 95) 887 (+178, 18) 2268 (+833, 94) 127 (+157, 62) 450 (+18, 18) 1387 (+18, 18) 6366 (+3435, 1875)
3296 719 2736 56 151 341 816 608
697 1388 679 1162 2075 3503 922 742 5084
1091 576 1039 1954 1955 888 338 1239 4594 105 786 2413 5528
Nash–Sutcliffe
0.74 94.14 0.45 * * * 0.70 0.91 0.09 0.73 3.86
0.83 41.32 17.82 * 0.28 0.90 0.76 0.01 0.49 0.12 1.54
0.32 43.04 142.02 1.51 0.74 0.73 0.82 0.83 0.72 0.81 0.81 0.73
0.09 * 5.19 0.64 0.10 0.65 0.45 0.11 0.16 0.28 0.73 0.39 0.20 0.75 0.70 0.49
Spearman’s rank (sig.)
0.548 (0.160)
0.714 (0.047)
0.783 (0.013)
0.908 (0.000)
Author's personal copy
263
L. Turnbull et al. / Journal of Hydrology 388 (2010) 258–272
0.61–0.8 indicates a substantial chance of agreement and 0.81–1 indicates an almost perfect chance of agreement (Landis and Koch, 1977). The erosion component of Mahleran was evaluated against data collected on sediment loss from the four plots, for runoff events of varying magnitude and duration (Turnbull et al., 2010). The dissolved and particulate-bound nutrient transport model was evaluated against measured nutrient losses at time increments during
the runoff events, again, for runoff events of varying magnitude and duration Turnbull et al. (in preparation-b).
4. Results and discussion In the following sections the hydrological, sediment transport, dissolved nutrient and particulate-nutrient-transport components
Fig. 1. Examples of monitored and modelled hydrographs for plots 1–4 (grass to shrub) that had a Nash and Sutcliffe coefficient of efficiency of >0.65. The black line shows the modelled hydrograph. The grey line is the monitored hydrograph and the circles show the precise points at which measurements were taken. The dashed grey lines show upper and lower uncertainty bounds of the monitored hydrograph.
10000
Modelled runoff (l)
10000
1:1
Plot 1 8000
8000
6000
6000
4000
4000
2000
2000
0
0 0
2000
4000
6000
8000
10000
10000
0
2000
4000
6000
8000
10000
1:1
Plot 3 Modelled runoff (l)
1:1
Plot 2
1:1
Plot 4
8000
8000
6000
6000
4000
4000
2000
2000
0
10000
0 0
2000
4000
6000
Monitored runoff (l)
8000
10000
0
2000
4000
6000
8000
10000
Monitored runoff (l)
Fig. 2. Monitored and modelled runoff (l), with bars showing potential uncertainty in the volume of total runoff that was monitored derived using the approach outlined in Turnbull et al. (2010).
Author's personal copy
264
L. Turnbull et al. / Journal of Hydrology 388 (2010) 258–272
of the model are evaluated against data collected on water, sediment and nutrient fluxes over the grassland to shrubland transition at the Sevilleta (Turnbull et al., 2010, in preparation-b). Suitable field-based observations could not be collected to evaluate quantitatively the within-plot, spatial predictions of erosion and particulate-bound nutrients, although the spatial predictions of erosion were evaluated in Wainwright et al. (2008b). 4.1. Performance of the hydrological component of Mahleran The hydrological characteristics of the events used to test the performance of the hydrological component of Mahleran are shown in Table 2. The Nash–Sutcliffe (NS) model efficiency values are also shown in Table 2. Model performance was considered to be good with a NS model efficiency of NS P 0.65 (see Fig. 1). The performance of the hydrological components of Mahleran over the range of events tested is markedly varied (see Supplementary material). While the model performed reasonably well for some of the smaller runoff events on most of the plots (e.g. plot 2 on 15/08/2006, plot 3 on 29/08/2006 and plot 4 29/07/2006), it appears that modelling the smaller runoff events is more prone to error than that of larger magnitude events (see also Wainwright et al., 2008c). One explanation for the over-estimation of the modelled hydrograph for small events is the over-representation of the connectivity of runoff-generating areas, possibly due to the resolution of the digital elevation model, which, although it is relatively fine-scale (0.5 m2) is still perhaps too coarse to represent the finescale microtopography that is a feature of semi-arid grasslands (Parsons et al., 1997). Furthermore, pit (topographic lows or sinks) removal of the DEM is carried out in the pre-processing of the elevation data, which potentially results in the over-estimation of the connectivity of runoff-generating areas, since typically, during rainfall events pits store water resulting in localised disconnectivity of runoff-generating areas. This problem arises from not knowing a priori which pits actually occur in the landscape (particularly as their occurrence may be dynamic) and which are simply artefacts of the production of the DEM. Since reasonable confidence can be placed in the parameterization of soil moisture in the testing of the model here (as demonstrated by the good predictions of runoff generation), some other aspect of model parameterization or process representation under certain conditions must be inadequate or poorly represented in the model structure. It is also possible that small-scale variability of parameters at a resolution less than 0.5 m2 have an effect on the hydrological modelling (Müller et al., 2007a), although the existence of such variability could not be quantified. For plot 4, the model tends to underestimate the steepness of the rising limb of the hydrograph, which has the effect of causing the timing of the modelled peak discharge to be later than the monitored value. The late modelling of the peak discharge has the effect of a lagged recessional limb, which is further emphasized by the kinematic wave model (Singh, 1996). Thus, it is apparent that some factor that causes the rapid onset of runoff generation over shrub-covered surfaces is either not adequately represented in the model or an aspect of the model parameterization is incorrect. The poor simulation of the rising limb of the hydrograph for plot 4 is indicative of either poor prediction of infiltration or ff being set too high. The latter seems less likely, as Eq. (8) was derived for a plot with conditions most similar to plot 4. It is possible that within-storm rainfall characteristics, at a temporal resolution finer the 1-min resolution at which rainfall was monitored, will affect the quality of model output (Wainwright and Parsons, 2002), since intense bursts of rainfall experienced at the sub 1-min resolution may be significant in terms of runoff generation. Other processes, such as the development and breakdown of surface seals throughout runoff events may also affect infiltra-
tion rates and hence runoff generation (Luk et al., 1993). At present, Mahleran does not account for processes that might lead to the development of surface seals within events. However, before a process-based representation of surface sealing can be incorporated into Mahleran, field-based understanding of these processes and ability to describe these processes numerically must be improved. A comparison of monitored versus modelled total discharge for the events tested over each plot (Fig. 2) shows that although there are discrepancies between the monitored and modelled hydrographs, overall the simulation of total discharge is very good over plots 3 and 4, but variable over plots 1 and 2, with NS model efficiencies of 0.09, 0.49, 0.81, and 0.73 for plots 1 to 4 respectively (see Table 2). In order to take into account the potential uncertainty in the monitored data, NS model-efficiency statistics were determined for the upper and lower-limit Q margins. For the upper uncertainty limits, the NS model efficiencies were 0.73, 0.12, 0.81 and 0.70 for plots 1–4 respectively, and for the lower limit uncertainty the NS model efficiencies were 3.86, 1.34, 0.73 and 0.49 for plots 1–4 respectively. The Spearman’s rank-correlation coefficient for each plot (Table 2) shows that the modelled discharges are significantly correlated with monitored discharges, for plots 2, 3 and 4 (p < 0.05), and that overall the model simulates the correct ranking of events in terms of amount of total discharge compared to that monitored. Most difficulty was encountered in modelling the smaller runoff events and events over plots 1 and 2. It is likely that the model structure is largely responsible for the poorer simulation of smaller runoff events, due to the difficulty of simulating correctly the infiltration rates. Wainwright et al., (2002) demonstrated that because of error propagation in the underlying kinematic wave model structure, runoff estimates tend to be poorest at the start of the event, and the quality of runoff estimates is highly sensitive to uncertainty in infiltration estimates. Since the amount of runoff generated during smaller events is highly sensitive to infiltration rates, particularly as the soil is wetting-up, it is likely that the aforementioned difficulties are responsible for these poorer simulations. Furthermore, empirical relationships used in the model (such as that used to determine saturated hydraulic conductivity) were derived from experiments carried out over shrubland (see Parsons et al., 1996), and may thus account in part for the poorer simulation of hydrology over the grass-dominated plots. All too frequently, the testing of the performance of the withinplot modelling of hydrology against observations is ignored, since monitoring within-plot spatial variations in hydrology is notoriously difficult to achieve (Wainwright et al., 2008b). However, it is important to determine if the spatial dynamics of the hydrolog-
Table 3 v2 and j statistics showing the testing of modelled predictions of flow (F) and no flow (NF) against observed conditions at 15 locations within each plot for event for which data are available. Plot/event
Plot Plot Plot Plot Plot Plot Plot Plot Plot Plot Plot Plot
1 1 2 2 2 3 3 3 3 4 4 4
29/08/06 07/09/06 07/09/05 15/08/06 29/08/06 31/07/06 11/08/06 29/08/06 07/09/06 05/07/06 23/08/06 29/08/06
Monitored:Modelled NF:NF
NF:F
F:NF
F:F
2 2 3 5 5 4 3 6 1 1 1 2
3 0 2 2 0 0 1 0 0 2 1 1
1 3 1 5 4 5 5 3 3 1 2 1
9 10 9 3 6 6 6 6 11 11 11 11
v2
p
j
1.875 4.615 4.261 0.134 5.000 6.636 1.029 6.667 2.946 1.298 1.298 5.104
0.171 0.032 0.039 0.714 0.025 0.057 0.310 0.010 0.060 0.255 0.255 0.024
0.33 0.47 0.33 0.09 0.50 0.39 0.22 0.62 0.33 0.29 0.23 0.58
Author's personal copy
265
L. Turnbull et al. / Journal of Hydrology 388 (2010) 258–272
ical modelling are correct, particularly given the importance of the localised redistribution of resources that govern the finer-scale biotic–abiotic interactions during runoff events. An attempt is made here to test the modelled spatial variations in hydrology within-plot using a qualitative approach. In the following analysis, only the runoff events with the highest NS efficiency values for each plot hydrograph are tested, since we are seeking to determine how successful spatial predictions of flow within the plot are for events that are simulated well at the plot outlet. The qualitative testing of the modelled spatial hydrology uses point based observations of hydrology, describing whether or not flow was monitored at 15 points within each plot using mini-samplers (Parsons et al., 2003). Since the other components of Mahleran are driven by hydrology, it is assumed that if the spatial simulation of hydrology is reasonable, then so too will the spatial
simulation of erosion, particulate-bound and dissolved nutrient fluxes. However, the need remains for the successful collection of spatial data on water, nutrient and sediment fluxes in order to test quantitatively how well the current level of process representation in Mahleran is able to simulate correctly the spatial fluxes. The modelled spatial dynamics of runoff over all the plots vary considerably. Contingency tables of spatial flow results are shown in Table 3. The j coefficients for the events tested show that the model produces varied results in predicting the spatial distribution in runoff. Overall, there was a tendency for the model to predict no flow at locations where flow was observed over the transition plots (plots 2 and 3), which indicates the need for improved data on model parameterisation for mixed grassland and shrubland hillslopes. The lowest j coefficient (0.09) was derived for Plot 2 on 15/08/2006, which, of all the events evaluated, had the lowest total
Table 4 Comparison of monitored and modelled erosion and particle-size distribution of eroded sediment. Spearman’s rank-correlation coefficients and Nash–Sutcliffe model-efficiency statistics of monitored versus modelled sediment yield are given for each plot. Plot/event date
Calibrated plot runoff (l)
Plot 1 07/09/05 29/09/05 05/07/06 28/07/06 01/08/06 11/08/06 29/08/06 07/09/06
Observed sediment yield (kg)
4.10 0.10 0.90 0.20 1.50 0.40 2.00 0.90
Proportions by model sediment size class Observed u=1
u=2
u=3
u=4
u=5
u=6
0.53 – 0.43 0.29 0.41 0.29 0.60 0.50
0.40 – 0.44 0.60 0.46 0.57 0.34 0.41
0.04 – 0.05 0.08 0.07 0.09 0.03 0.05
0.03 – 0.04 0.02 0.06 0.03 0.03 0.03
0.00 – 0.04 0.00 0.01 0.01 0.01 0.01
0.00 – 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.02 0.00 0.01
0.00 0.00 0.00 0.00 0.00 0.00
Modelled sediment yield (kg)
Proportions by model sediment size class Modelled u=1
u=2
u=3
u=4
u=5
u=6
0.29 0.38 0.42 0.69 0.41 0.47 0.26 0.31
0.68 0.62 0.58 0.31 0.58 0.53 0.73 0.68
0.02 0.00 0.00 0.00 0.00 0.00 0.01 0.01
0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
3.41 0.01 0.09 0.23 1.06 0.99
0.39 0.95 0.65 0.53 0.42 0.40
0.60 0.05 0.35 0.47 0.58 0.59
0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00
27.1 0.2 0.6 0.1 0.3 0.4 2.6 1.7
Summary statistics for all events: Spearman’s rank 0.826 sig. 0.011 NS efficiency 5.78 Plot 2 07/09/05 04/08/06 11/08/06 15/08/06 29/08/06 07/09/06
6.25 0.02 0.53 0.54 1.26 1.12
0.46 0.73 0.32 0.38 0.46 0.39
0.42 0.25 0.49 0.47 0.42 0.49
0.07 0.01 0.11 0.07 0.06 0.07
0.05 0.00 0.07 0.06 0.06 0.05
Summary statistics for all events: Spearman’s rank 1.000 sig.
Author's personal copy
Journal of Hydrology 388 (2010) 258–272
Contents lists available at ScienceDirect
Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol
Hydrology, erosion and nutrient transfers over a transition from semi-arid grassland to shrubland in the South-Western USA: A modelling assessment Laura Turnbull a,*, John Wainwright b, Richard E. Brazier c a
Global Institute of Sustainability, Arizona State University, 800 South Cady Mall, Tempe, AZ 85287-5402, USA Sheffield Centre for International Drylands Research, Department of Geography, University of Sheffield, Winter Street, Sheffield S10 2TN, UK c Department of Geography, University of Exeter, Amory Building, Rennes Drive, Exeter EX4 4RJ, UK b
a r t i c l e
i n f o
Article history: Received 26 August 2009 Received in revised form 27 April 2010 Accepted 3 May 2010
This manuscript was handled by L. Charlet, Editor-in-Chief, with the assistance of Bernhard Wehrli, Associate Editor Keywords: Runoff Erosion Land degradation Nitrogen Phosphorus Ecogeomorphology
s u m m a r y Land degradation in arid and semi-arid areas, as a consequence of the invasion of grasslands by shrubs, is often associated with an increase in runoff and erosion and a change in nutrient transport. Modelling of nutrient transport during runoff events (in particular particulate-bound nutrients), is especially important, since the spatial redistribution of nutrients (in addition to water and sediment) can have significant implications for vegetation dynamics in these ecosystems. In this study, Mahleran (Model for Assessing Hillslope to Landscape Erosion Runoff, And Nutrients) is extensively evaluated against runoff and erosion data from four plots (representative of different stages of land degradation) over a transition from grassland to shrubland, at the Sevilleta National Wildlife Refuge in New Mexico, USA. A new particulate-bound nutrient module was developed to include a representation of particulate-bound nutrient dynamics, which is an important form of nutrient transport in these ecosystems. Understanding dynamics of both dissolved and particulate-bound nutrient dynamics during runoff events is imperative, because of their differing roles in terms of nutrient bioavailability and potential implications for plant dynamics. Results of the model evaluation show that the runoff and erosion components of Mahleran perform reasonably well, as does the new particulate-bound nutrient sub-model, though not consistently. Performance of the particulate-bound nutrient model was better for the end-member plots, because of better parameterization data available for end-member vegetation types. Since the particulate-bound nutrient sub-model is by necessity strongly dependent on the simulated erosion rate, the performance of the particulate-bound nutrient model is dependent on the performance of the erosion component of Mahleran, so that when erosion is well represented by the model, so typically are particulate nutrient transfers. The performance of the dissolved nutrient component of Mahleran was poor in this application, which indicates that the process representation for this semi-arid environment and the parameterisation of the dissolved nutrient component were inadequate. Results from the model evaluation suggest that an improved understanding of dissolved nutrient dynamics during runoff events and simulation if interevent nutrient dynamics is required, in order to improve the level of process representation within modelling approaches and thus the ability to simulate dissolved nutrient dynamics and their subsequent effects on other ecosystem processes. Ó 2010 Elsevier B.V. All rights reserved.
1. Introduction Recent research on changes in runoff and erosion over a grassland to shrubland transition (Turnbull et al., 2010, in preparation-a,-b) has provided new insights into changes in ecosystem structure and function. The observed changes in ecosystem structure and function over a semi-arid, grass–shrub ecotone demonstrate that changes in ecosystem structure and function do not necessarily follow a linear continuum across a grassland to shrubland transition. As the distribution and percentage cover of * Corresponding author. E-mail address: [email protected] (L. Turnbull). 0022-1694/$ - see front matter Ó 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.jhydrol.2010.05.005
vegetation decreased over the grassland to shrubland transition, concurrent changes in the connectivity of these surface properties have profound effects on hydrological function (Turnbull et al., 2010). Determining more fully the interactions between ecosystem structure and function, and interactions/feedbacks between biotic and abiotic components of the ecosystem is imperative in order to ascertain the effects that climate and land-use change will have on these ecosystems in the future (Turnbull et al., 2008). Modelling-based studies provide a potential means of expanding further our understanding of interactions and feedbacks between different components of an ecosystem, under a multitude of climatic and land-use conditions, that would not be feasible to investigate using field-based experimentation alone. However,
Author's personal copy
L. Turnbull et al. / Journal of Hydrology 388 (2010) 258–272
the limitation at present of carrying out such modelling-based studies is that the transfers of matter and energy among units are rarely included in models of vegetation change, but are critically important in arid ecosystems in which resource redistribution plays a significant role (Peters and Herrick, 2001). Therefore, the development of existing modelling approaches is imperative in order to understand further the dynamics of land degradation in semi-arid areas. The aim of this study is to develop a process-based model to understand better the dynamics of dissolved and particulate-bound nutrient transport processes in runoff, with particular emphasis on the testing of the model algorithms. The specific model to be developed is Mahleran (Model for Assessing Hillslope to Landscape Erosion, Runoff And Nutrients; Müller et al., 2007a,b; Wainwright et al., 2008a,b). 2. Modelling approach In the following sections the key components of Mahleran are outlined: (i) the hydrology component (Wainwright and Parsons, 2002), (ii) the erosion component (Wainwright et al., 2008b), (iii) the dissolved nutrient component (Müller et al., 2007a,b), and (iv) a new particulate-bound nutrient component of Mahleran which is developed here, based on the improved understanding of particulate-bound nutrient fluxes outlined in Turnbull et al. (in preparation-b). 2.1. Hydrology The hydrological component of Mahleran (Wainwright and Parsons, 2002) uses a simple infiltration model to generate infiltration- and saturation-excess runoff. Infiltration-excess runoff is generated by determining the difference between the rainfall rate and the infiltration rate. The infiltration rate is simulated using the Smith-Parlange approach (Smith and Parlange, 1978). Runoff is routed over the hillslope using a kinematic wave approximation of the St Venant equations (Wainwright and Parsons, 2002), with flow routing in the direction of steepest descent from cell to cell (in cardinal directions) over a finite-difference grid (Scoging et al., 1992). The flow is routed over the grid using a finite-difference solution (Euler backward difference form: Scoging et al., 1992). Flow velocity is determined dynamically using the DarcyWeisbach flow equation. This modelling approach differs from most kinematic wave approaches, which assume a fixed rating function, and thus poorly represent feedbacks of flow hydraulics with runoff and infiltration. Full details of the hydrology and hydraulics components are found in Wainwright and Parsons (2002); Parsons et al. (1997), and Wainwright et al., 2008b. 2.2. Erosion Mahleran considers the interaction of both concentrated and unconcentrated erosion as bedload and in suspension, which is particularly important since the relative balance of the different processes is the most critical control on the resulting pattern of erosion (Wainwright et al., 2008a). Mahleran employs an approach to erosion modelling that is based upon the concept of entrainment and travel distances of particles from six different size classes (u = 1 is 12 mm). Sediment transport is modelled by a form of the continuity equation, in which sediment belonging to these different size classes are simulated separately (Wainwright et al., 2008a).
@hs;u @qs;u þ eu þ du ¼ 0 @t @x
ð1Þ
259
where hs is the equivalent depth of sediment transport (m), qs is the unit discharge of sediment (m2 s-1), e is the rate of erosion of the surface (m s1), d is the rate of deposition (m s1), u is an index relating to a specific size class of sediment, t is time (s) and x is location in space (m). To account for the temporal scaling characteristics of this model, qs,u is defined as
qs;u ¼ hs;u v p;u
ð2Þ
where vp,u is the virtual velocity of transport of sediment of particlesize class u (m s1), which can be derived from empirical observations of sediment transport distance per unit time or by theoretical calculation by differentiation of transport distance per unit of time (Wainwright et al., 2008a). Mahleran simulates three sediment detachment and transport conditions (see Wainwright et al., 2008a, for more detail) (i) erosion as a function of raindrop detachment and transport occurs by splash when no flow is present; (ii) unconcentrated overland flow erosion simulated using raindrop detachment rates that are modified to account for the protective effects of the surface-water layer; (iii) transport by mixed unconcentrated and concentrated flow (500 6 Re 6 2000) (where Re is the flow Reynolds number); or (iv) transport by concentrated flow (Re > 2000). Deposition is modelled similarly, using the transport-distance approach, since the distribution function of travel differences of particles under specific transport mechanisms and flow conditions enables the deposition rate to be determined directly at each point along the transport pathway (Wainwright et al., 2008a). Rates of deposition depend on whether material being moved is in splash, as bedload in unconcentrated or concentrated overland flows or in suspension (Wainwright et al., 2008a).
2.3. Nutrients Dissolved nutrient concentrations in runoff are affected by the nutrient content in the upper soil layers, which depends upon the nutrient content of the soil and the soil-exchange characteristics (Barger et al., 2006). A common approach employed to model the transfer of solutes from soil solution to runoff is based on the principle that nutrients in the soil solution will diffuse into runoff water. The simplest assumption in determining the mass transfer of nutrients between the soil surface and overland flow is that the rate of exchange is directly proportional to the difference in concentration between the soil solution and the surface water (Walton et al., 2000b). Wallach and van Genuchten (1990) developed a modelling approach in which the diffusion of solutes from the soil solution to runoff water was coupled with a rate-limited chemical transfer process across a laminar boundary at the soil surface/runoff water interface. Solute release into runoff is limited by the initial mass of nutrient near the soil surface and the depth at which surface water interacts with the soil (Havis et al., 1992). The depth of soil that interacts with overland flow through turbulent mixing, convection and diffusion, that contributes chemicals to the surface water during a rainfall event, is defined as the mixing zone (Havis et al., 1992). Solute exchange through this boundary layer was modelled using a mass-transfer coefficient, K (Wallach and van Genuchten, 1990). Walton et al. (2000b) modelled the mass transfer of nutrients between the soil surface and overland flow, assuming that the rate of change of mass transfer is directly proportional to the difference between the concentration in the soil solution and the concentration in the surface water.
Author's personal copy
260
L. Turnbull et al. / Journal of Hydrology 388 (2010) 258–272
Müller et al. (2007a,b) developed the dissolved nutrient-transport component of Mahleran, in which NO3, NH4, and PO4 are modelled by an advection–dispersion model (Havis et al., 1992; Walton et al., 2000b). In the advection–dispersion model the mass transfer of nutrients from the soil surface to overland flow is driven by: (i) diffusion of dissolved nutrients from the soil interstices by movement of soil water into the overland flow, (ii) desorption of the nutrients from soil particles into the overland flow, (iii) dissolution of solid phase nutrients into the soil water or overland flow, and (iv) scouring of solid phase nutrients by hydraulic forces and subsequent transport and moving dissolution (Müller et al., 2007a). An advection–dispersion model is used to model nutrient transport, in which a mass-transfer coefficient is used to lump together the aforementioned mechanisms of mass transfer (Wallach and van Genuchten, 1990). The nutrient mass balance is given by.
@ðChÞ @ðqCÞ þ ¼ KðC s CÞ fC @t @x
ð3Þ
where C is the nutrient concentration in overland flow (mg l1), Cs is the nutrient concentration in the soil-solution zone (mg l1), h is the depth of overland flow (mm), q is the unit discharge (mm2 s1) and f is the infiltration rate (mm h1) (Müller et al., 2007a). However, in this modelling approach, the effect of the nutrient content of the incoming rainfall is neglected. Turnbull et al. (in preparation-b) demonstrate that the concentrations of dissolved nutrients in rainfall are often of the same order of magnitude as those measured in runoff. Therefore, when modelling diffusion of nutrients from the soil-solution zone into runoff, the consideration of the concentration of nutrients in rainfall is important. Otherwise, it is possible that net losses of nutrients from the soil-solution zone will be over-estimated. To account for the effects of rainfall-nutrient concentration on the diffusion of nutrients from the soil-solution zone into runoff, Eq. (3) can be re-written as
@ðChÞ @ðqCÞ þ ¼ KðC s CÞ fC þ RC r @t @x
ð4Þ
where R is the rainfall rate (mm h1) and Cr is the concentration of the nutrient in the rainfall (mg l1) (Havis et al., 1992; Walton et al., 2000a). Havis et al. (1992) related the potential mass (M) of a nutrient over an area (A).Cs is related to the nutrient concentration of the soil by
M ¼nA
Z
x
C s ðzÞdz
ð5Þ
0
where M is the mass of available nutrient in the soil (g), A is the model cell size (mm2), n is porosity (%), x is the depth of the mixing zone (mm) and z is the depth through the soil (mm). Once soluble nutrients are entrained by runoff, they are routed conservatively (e.g. Viney et al., 2000). Although data from field studies, for example data presented in Turnbull et al. (in preparation-b) and Lister (2007) indicate that NO3–N, NH4–N and PO4–P in runoff are not conservative, at present there is insufficient information to represent this non-conservative behaviour in the modelling approach, which would possibly induce greater errors than employing the simple conservative approach. The extent to which this assumption may be the case will be evaluated in the context of the results of this modelling exercise Turnbull et al. (in preparation-b) demonstrate the significance of particulate-bound nutrient transport, which must thus be incorporated into the modelling approach. Very few semi-arid, field-based studies have investigated particulate-bound nutrient transfers in runoff. Since the significance of particulate-bound nutrient transfers in runoff has not been well recognised, very few models have represented this particulate phase of nutrient transfer. A key problem with employing detailed, process-based approaches to the modelling of nutrient transfers in semi-arid
environments is that understanding of the processes of nutrient transfers between states remains limited. In addition, the use of a modelling approach that has high parameterization demands is impractical, and may induce increased uncertainty into the processes being modelled (where data are not sufficient to constrain the additional uncertainty introduced by extra parameters, or more complex model structures – see Brazier et al., 2000).Fierer and Gabet (2002) calculated sediment-bound nutrient loss as the product of sediment detachment rate and the percentage C and percentage N associated with the sediment. They determined the percentage of C and N associated with the sediment by carrying out rainfall-simulation experiments at their semi-arid field site in California, where they sampled the sediment transported by runoff, which was then weighed to determine the sediment concentration in runoff and analysed for C and N content on a Fisons NA1500 C/N analyzer. The approach used to model particulate-bound nutrient transport in Mahleran elaborates upon that of Fierer and Gabet (2002) by determining the particulate-bound nutrient transport as a function of the nutrient concentration of sediment, s, in each particle-size class and the amount of sediment transported in each particle-size class. The flux of a nutrient of type n, adsorbed onto sediment in each particle-size class per time-step is related to the discharge of sediment (qs,u) and the nutrient concentration of sediment in each particle-size class (Cn,u) by:
qn;u ¼ qs;u C n;u
ð6Þ
The modelling of particulate-bound nutrient fluxes is conservative, since it is assumed that there is no cycling of adsorption/ desorption of particulate-bound nutrients when they are being transported (Viney et al., 2000). 3. Parameterizing Mahleran for Grassland–Shrubland Plots at the Sevilleta LTER Site and methods of model testing In order to evaluate the different model components in the context of changing vegetation structures, Mahleran is parameterized for four sites over a transition from black grama (Bouteloua eriopoda) grassland to creosotebush (Larrea tridentata) shrubland at the Sevilleta Long-Term Ecological Research site in central New Mexico, USA, which is lies on the northern margin of the Chihuahuan desert. Sites were located within approximately 1 km of each other and have similar slopes and aspects. The Sevilleta receives an average of 256 mm of rainfall annually, of which 55% typically falls during the summer monsoon rainfall season (July–September) (Dahm and Moore, 1994). At each site, a runoff plot was constructed measuring 10 30 m, from which flow was measured at 1-min intervals. The outflow from each plot was instrumented with an ISCO 6700 autosampler which sampled flow leaving plots at 1-min intervals. In addition to this instrumentation, a collection tank at the base of each plot trapped all runoff and sediment exported from the plots, to provide a secondary measure of total plot runoff and a measure of the total amount of sediment eroded from each plot. A more detailed description of the field site is provided in Turnbull et al. (2010, in preparation-a,-b) and Sevilleta LTER (2009). 3.1. Hydrology Parameterizing final infiltration rate (Ks) is notoriously difficult since Ks is variable over space and time (Müller et al., 2008). Attempts have been made to determine Ks in the field, using infiltrometers, but this technique has been found to overestimate Ks (Müller, 2007). An alternative and preferable method of determining Ks is to carry out rainfall-simulation experiments. However, this approach is not straightforward and was not possible within
Author's personal copy
261
L. Turnbull et al. / Journal of Hydrology 388 (2010) 258–272
this study. Thus, Ks could not be parameterized directly at the field site and therefore was estimated by using the best available information for similar sites where desert pavement progressively occupies the surface as vegetation changes from grass to shrubs. Previous studies have shown that well-developed desert pavements exert a dominant control on infiltration rates (Abrahams et al., 1989). Following Wainwright et al. (2008a) the rainfall-simulation experiments of Parsons et al. (1996) were used to parameterize Ks (mm min1) by
K s ¼ 0:351 þ 0:010 rain 0:006P%
ð7Þ
where rain is rainfall intensity (mm h1) and P% is pavement cover (%). There are limitations inherent in using this relationship to determine Ks for the field site in question, since this relationship was derived for a different field site, albeit one with similar soils. However, the benefit of employing this relationship to determine Ks is that it accounts for the effect of varying stone-pavement cover and rainfall intensities on Ks. Although Ks is likely to vary over the grass to shrub transition, appropriate data are not available to constrain the likely differences, and therefore the above equation is used to parameterize Ks as a first approximation. The quality of this approximation is evaluated below Numerous approaches have been described previously for the determination of the friction factor, ff (see Wainwright et al., 2008b), which differ in their levels of complexity and parameterization requirements. The Darcy-Weisbach friction factor (dimensionless) can be parameterized from field experiments (e.g. Parsons et al., 1994), although these are very time consuming and labour intensive. As a first order approximation, the simple depth–feedback parameterization of Scoging et al. (1992) was used
ff ¼ 14 0:8h
ð8Þ
where h is the depth of runoff (mm). The depth to wetting front must be known a priori, and is a difficult parameter to measure (Wainwright et al., 2008c). Soil thickness represents effective depth to wetting front, which was thus estimated to be 50 cm, based on measurements of soil depth. The wetting front suction parameter was estimated to be 5 cm. The saturated soil-moisture content was set to 39%. Again, while saturated soil-moisture content is likely to be spatially variable, no data were collected to constrain this spatial variability, thus a uniform value was used. The kinetic energy (KE in J m2 mm1) of the rainfall arriving at the ground surface is related to vegetation characteristics. A function for the relationship between KE and vegetation has been derived from Wainwright et al. (1999)
KEv ¼ KEð1 8:1 103 V%Þ
ð9Þ
where KEv (J m2 mm1) is the KE of rainfall arriving at the ground surface that is adjusted to account for vegetation cover and V% is the vegetation canopy cover (%). This function was derived from experiments carried out for creosotebush. No data has been collected to determine the relationship between KE and black grama canopy cover. Since this function was derived for creosotebush, it accounts for effect of leaf drip energy, which will be negligible for grass-covered areas. However, during preliminary testing of the model the sensitivity of erosion to leaf drip energy over grass-covered areas was found to be very small. Therefore, the function for the relationship between KE and vegetation derived from Wainwright et al. (1999) is employed for all vegetation types. Analysis of the spatio-temporal volumetric soil-moisture content dataset (Turnbull et al., 2010) showed that there are often significant differences in soil-moisture content between bare and vegetated soils. In view of these findings, Mahleran was developed to improve the process representation of soil-moisture dynamics, to enable spatially distributed parameterization of initial soil-
moisture content. Thus, rather than a spatially uniform parameterization of soil-moisture content, as had been carried out previously (Wainwright et al., 2008b), in this modified version of Mahleran, soil-moisture content at the start of the event is parameterized using a spatio-temporal dataset of volumetric soil-moisture content for different surface-cover types (Turnbull et al., 2010). 3.2. Erosion Parameterization of the soil particle-size distribution in Mahleran is made spatially distributed, based on the stone-pavement cover (made up of u = 5and u = 6), which was mapped digitally for the field site from aerial photographs (Turnbull et al., in preparation-a). The proportion of sediment in the remaining particle-size classes is also spatially distributed, based on the particle-size characteristics of the soil measured over each plot (see Turnbull et al., in preparation-a). The parameterization of sediment detachment and transport parameters for each u class were taken from Wainwright et al. (2008b). 3.3. Nutrients In order to parameterize the nutrient content of sediment from each particle-size class, the data of Lister (2007) are used, obtained from rainfall-simulation experiments carried out over 1–m2 plots over grass and shrub end-members at the Jornada Experimental Range, New Mexico. Therefore, where grass formed the dominant vegetation cover (plots 1 and 2), the grass end-member parameterization was used, and where shrubs formed the dominant vegetation cover (plots 3 and 4), the shrub end-member parameterization was used. The nutrient content associated with each particle-size class of eroded sediment was determined. The flux (g min1) of particulate-bound nutrients is modelled here as a function of the total amount of sediment flux, the size composition of that sediment and the nutrient content of sediment in each particle-size class. To parameterize the mass-transfer coefficient (K), the values used by Müller (2007) were used (Table 1), which were found to perform well for a similar setting (Müller et al., 2007b). The nutrient concentration in the soil-solution zone (Cs) was determined using Eq. (5). The resultant values used to parameterize Cs for each site for each dissolved nutrient are shown in Table 1. 3.4. Analysis techniques used to test model performance The Nash and Sutcliffe (NS) coefficient of efficiency (Nash and Sutcliffe, 1970) is used to measure the performance of the model compared to the observed runoff, erosion and nutrient fluxes. Rigorous quality control of observed fluxes was carried out (see Turnbull et al., 2010) in order to ensure that model evaluations were carried out against data in which there was a high degree of confidence. The performance of the hydrology model is evaluated against the goodness of fit of the modelled hydrograph with the monitored hydrograph, and against the total volume of modelled runoff with that monitored. Nash–Sutcliffe model efficiencies were
Table 1 Parameterisation of Cs and K for NH4–N, NO3–N and PO4–P. The parameterisation of K is the same for all plots. Plot 1 is the grass end-member, and plot 4 the shrub endmember.
Cs (mg l1)
K (mm s1)
Plot Plot Plot Plot
1 2 3 4
NH4–N
NO3–N
PO4–P
11 12 9 11
11 15 10 11
61 61 61 61
0.00003
0.00007
0.00004
Author's personal copy
262
L. Turnbull et al. / Journal of Hydrology 388 (2010) 258–272
also calculated for the upper and lower uncertainty margins of the total runoff monitored from the plots (see Turnbull et al., 2010) in order that uncertainties in the monitored data are taken into account in the model evaluation. To test the spatial predictions of flow versus no flow conditions by the model, for events over each plot, contingency tables were used. In order to approximate the modelled ‘no flow’ conditions, it is assumed that a total event discharge of 0.005 m3 or less, over a cell of 0.5 m2 is analogous to no flow conditions, since total event
discharges below this amount will not have been sampled by minisamplers (Parsons et al., 2003) which were designed to monitor flow depth at fifteen points over each plot. Kappa coefficients were calculated to evaluate the agreement of monitored and modelled results against those that might be expected by chance (Landis and Koch, 1977; Monserud and Leemans, 1992). Values of j < 0 indicate a poor strength of agreement, 0–0.2 indicates a slight chance of agreement, 0.21–0.4 indicates a fair chance of agreement, 0.41–0.6 indicates a moderate chance of agreement,
Table 2 Summary of model simulation results for events monitored over plots 1–4. The Nash–Sutcliffe efficiency could not be calculated for events with an asterisk in the Nash–Sutcliffe column because equipment failures prevented fitting against the hydrograph. Nash–Sutcliffe model efficiencies were also calculated for the upper and lower bounds of uncertainty of monitored runoff. Spearman’s rank-correlation coefficients of monitored versus modelled runoff and Nash–Sutcliffe efficiency statistic for event totals of monitored and modelled runoff are given. Nash–Sutcliffe efficiency statistics are given for upper and lower uncertainty bounds of monitored data that were derived from analysis of hydrograph uncertainty. Plot/event Plot 1 07/09/05 29/09/05 05/07/06 28/07/06 01/08/06 11/08/06 29/08/06 07/09/06 All events Upper uncertainty Lower uncertainty Plot 2 07/09/05 31/07/06 01/08/06 04/08/06 11/08/06 15/08/06 29/08/06 07/09/06 All events Upper uncertainty Lower uncertainty Plot 3 28/09/05 29/09/05 05/07/06 28/07/06 31/07/06 01/08/06 11/08/06 29/08/06 07/09/06 All events Upper uncertainty Lower uncertainty Plot 4 07/09/05 28/09/05 29/09/05 05/07/06 28/07/06a 28/07/06b 29/07/06 31/07/06 01/08/06 04/08/06 11/08/06 23/08/06 29/08/06 All events Upper uncertainty Lower uncertainty
Total runoff (l) Monitored
Total runoff (l) Modelled
5230 (+2250, 1950) 124 (+18, 18) 277 (+18, 18) 177 (+18, 18) 1011 (+18, 18) 213 (+18, 18) 2820 (+1249, 639) 1426 (+18, 18)
9848 1590 372 572 1110 445 1818 1586
2294 (+540, 261) 83 (+17, 17) 464 (+17, 17) 38 (+16, 16) 149 (+17, 17) 391 (+17, 17) 995 (+17, 17) 1260 (+17, 17)
363 (+17, 17) 122 (+17, 17) 66 (+17, 17) 380 (+17, 17) 1799 (+17, 17) 2030 (+881, 0) 637 (+17, 17) 551 (+17, 17) 6550 (+766, 1623)
2174 (+635, 0.0) 142 (+18, 18) 443 (+18, 18) 2174 (+785, 0.0) 1504 (+907, 468) 452 (+378, 169) 199 (+237, 95) 887 (+178, 18) 2268 (+833, 94) 127 (+157, 62) 450 (+18, 18) 1387 (+18, 18) 6366 (+3435, 1875)
3296 719 2736 56 151 341 816 608
697 1388 679 1162 2075 3503 922 742 5084
1091 576 1039 1954 1955 888 338 1239 4594 105 786 2413 5528
Nash–Sutcliffe
0.74 94.14 0.45 * * * 0.70 0.91 0.09 0.73 3.86
0.83 41.32 17.82 * 0.28 0.90 0.76 0.01 0.49 0.12 1.54
0.32 43.04 142.02 1.51 0.74 0.73 0.82 0.83 0.72 0.81 0.81 0.73
0.09 * 5.19 0.64 0.10 0.65 0.45 0.11 0.16 0.28 0.73 0.39 0.20 0.75 0.70 0.49
Spearman’s rank (sig.)
0.548 (0.160)
0.714 (0.047)
0.783 (0.013)
0.908 (0.000)
Author's personal copy
263
L. Turnbull et al. / Journal of Hydrology 388 (2010) 258–272
0.61–0.8 indicates a substantial chance of agreement and 0.81–1 indicates an almost perfect chance of agreement (Landis and Koch, 1977). The erosion component of Mahleran was evaluated against data collected on sediment loss from the four plots, for runoff events of varying magnitude and duration (Turnbull et al., 2010). The dissolved and particulate-bound nutrient transport model was evaluated against measured nutrient losses at time increments during
the runoff events, again, for runoff events of varying magnitude and duration Turnbull et al. (in preparation-b).
4. Results and discussion In the following sections the hydrological, sediment transport, dissolved nutrient and particulate-nutrient-transport components
Fig. 1. Examples of monitored and modelled hydrographs for plots 1–4 (grass to shrub) that had a Nash and Sutcliffe coefficient of efficiency of >0.65. The black line shows the modelled hydrograph. The grey line is the monitored hydrograph and the circles show the precise points at which measurements were taken. The dashed grey lines show upper and lower uncertainty bounds of the monitored hydrograph.
10000
Modelled runoff (l)
10000
1:1
Plot 1 8000
8000
6000
6000
4000
4000
2000
2000
0
0 0
2000
4000
6000
8000
10000
10000
0
2000
4000
6000
8000
10000
1:1
Plot 3 Modelled runoff (l)
1:1
Plot 2
1:1
Plot 4
8000
8000
6000
6000
4000
4000
2000
2000
0
10000
0 0
2000
4000
6000
Monitored runoff (l)
8000
10000
0
2000
4000
6000
8000
10000
Monitored runoff (l)
Fig. 2. Monitored and modelled runoff (l), with bars showing potential uncertainty in the volume of total runoff that was monitored derived using the approach outlined in Turnbull et al. (2010).
Author's personal copy
264
L. Turnbull et al. / Journal of Hydrology 388 (2010) 258–272
of the model are evaluated against data collected on water, sediment and nutrient fluxes over the grassland to shrubland transition at the Sevilleta (Turnbull et al., 2010, in preparation-b). Suitable field-based observations could not be collected to evaluate quantitatively the within-plot, spatial predictions of erosion and particulate-bound nutrients, although the spatial predictions of erosion were evaluated in Wainwright et al. (2008b). 4.1. Performance of the hydrological component of Mahleran The hydrological characteristics of the events used to test the performance of the hydrological component of Mahleran are shown in Table 2. The Nash–Sutcliffe (NS) model efficiency values are also shown in Table 2. Model performance was considered to be good with a NS model efficiency of NS P 0.65 (see Fig. 1). The performance of the hydrological components of Mahleran over the range of events tested is markedly varied (see Supplementary material). While the model performed reasonably well for some of the smaller runoff events on most of the plots (e.g. plot 2 on 15/08/2006, plot 3 on 29/08/2006 and plot 4 29/07/2006), it appears that modelling the smaller runoff events is more prone to error than that of larger magnitude events (see also Wainwright et al., 2008c). One explanation for the over-estimation of the modelled hydrograph for small events is the over-representation of the connectivity of runoff-generating areas, possibly due to the resolution of the digital elevation model, which, although it is relatively fine-scale (0.5 m2) is still perhaps too coarse to represent the finescale microtopography that is a feature of semi-arid grasslands (Parsons et al., 1997). Furthermore, pit (topographic lows or sinks) removal of the DEM is carried out in the pre-processing of the elevation data, which potentially results in the over-estimation of the connectivity of runoff-generating areas, since typically, during rainfall events pits store water resulting in localised disconnectivity of runoff-generating areas. This problem arises from not knowing a priori which pits actually occur in the landscape (particularly as their occurrence may be dynamic) and which are simply artefacts of the production of the DEM. Since reasonable confidence can be placed in the parameterization of soil moisture in the testing of the model here (as demonstrated by the good predictions of runoff generation), some other aspect of model parameterization or process representation under certain conditions must be inadequate or poorly represented in the model structure. It is also possible that small-scale variability of parameters at a resolution less than 0.5 m2 have an effect on the hydrological modelling (Müller et al., 2007a), although the existence of such variability could not be quantified. For plot 4, the model tends to underestimate the steepness of the rising limb of the hydrograph, which has the effect of causing the timing of the modelled peak discharge to be later than the monitored value. The late modelling of the peak discharge has the effect of a lagged recessional limb, which is further emphasized by the kinematic wave model (Singh, 1996). Thus, it is apparent that some factor that causes the rapid onset of runoff generation over shrub-covered surfaces is either not adequately represented in the model or an aspect of the model parameterization is incorrect. The poor simulation of the rising limb of the hydrograph for plot 4 is indicative of either poor prediction of infiltration or ff being set too high. The latter seems less likely, as Eq. (8) was derived for a plot with conditions most similar to plot 4. It is possible that within-storm rainfall characteristics, at a temporal resolution finer the 1-min resolution at which rainfall was monitored, will affect the quality of model output (Wainwright and Parsons, 2002), since intense bursts of rainfall experienced at the sub 1-min resolution may be significant in terms of runoff generation. Other processes, such as the development and breakdown of surface seals throughout runoff events may also affect infiltra-
tion rates and hence runoff generation (Luk et al., 1993). At present, Mahleran does not account for processes that might lead to the development of surface seals within events. However, before a process-based representation of surface sealing can be incorporated into Mahleran, field-based understanding of these processes and ability to describe these processes numerically must be improved. A comparison of monitored versus modelled total discharge for the events tested over each plot (Fig. 2) shows that although there are discrepancies between the monitored and modelled hydrographs, overall the simulation of total discharge is very good over plots 3 and 4, but variable over plots 1 and 2, with NS model efficiencies of 0.09, 0.49, 0.81, and 0.73 for plots 1 to 4 respectively (see Table 2). In order to take into account the potential uncertainty in the monitored data, NS model-efficiency statistics were determined for the upper and lower-limit Q margins. For the upper uncertainty limits, the NS model efficiencies were 0.73, 0.12, 0.81 and 0.70 for plots 1–4 respectively, and for the lower limit uncertainty the NS model efficiencies were 3.86, 1.34, 0.73 and 0.49 for plots 1–4 respectively. The Spearman’s rank-correlation coefficient for each plot (Table 2) shows that the modelled discharges are significantly correlated with monitored discharges, for plots 2, 3 and 4 (p < 0.05), and that overall the model simulates the correct ranking of events in terms of amount of total discharge compared to that monitored. Most difficulty was encountered in modelling the smaller runoff events and events over plots 1 and 2. It is likely that the model structure is largely responsible for the poorer simulation of smaller runoff events, due to the difficulty of simulating correctly the infiltration rates. Wainwright et al., (2002) demonstrated that because of error propagation in the underlying kinematic wave model structure, runoff estimates tend to be poorest at the start of the event, and the quality of runoff estimates is highly sensitive to uncertainty in infiltration estimates. Since the amount of runoff generated during smaller events is highly sensitive to infiltration rates, particularly as the soil is wetting-up, it is likely that the aforementioned difficulties are responsible for these poorer simulations. Furthermore, empirical relationships used in the model (such as that used to determine saturated hydraulic conductivity) were derived from experiments carried out over shrubland (see Parsons et al., 1996), and may thus account in part for the poorer simulation of hydrology over the grass-dominated plots. All too frequently, the testing of the performance of the withinplot modelling of hydrology against observations is ignored, since monitoring within-plot spatial variations in hydrology is notoriously difficult to achieve (Wainwright et al., 2008b). However, it is important to determine if the spatial dynamics of the hydrolog-
Table 3 v2 and j statistics showing the testing of modelled predictions of flow (F) and no flow (NF) against observed conditions at 15 locations within each plot for event for which data are available. Plot/event
Plot Plot Plot Plot Plot Plot Plot Plot Plot Plot Plot Plot
1 1 2 2 2 3 3 3 3 4 4 4
29/08/06 07/09/06 07/09/05 15/08/06 29/08/06 31/07/06 11/08/06 29/08/06 07/09/06 05/07/06 23/08/06 29/08/06
Monitored:Modelled NF:NF
NF:F
F:NF
F:F
2 2 3 5 5 4 3 6 1 1 1 2
3 0 2 2 0 0 1 0 0 2 1 1
1 3 1 5 4 5 5 3 3 1 2 1
9 10 9 3 6 6 6 6 11 11 11 11
v2
p
j
1.875 4.615 4.261 0.134 5.000 6.636 1.029 6.667 2.946 1.298 1.298 5.104
0.171 0.032 0.039 0.714 0.025 0.057 0.310 0.010 0.060 0.255 0.255 0.024
0.33 0.47 0.33 0.09 0.50 0.39 0.22 0.62 0.33 0.29 0.23 0.58
Author's personal copy
265
L. Turnbull et al. / Journal of Hydrology 388 (2010) 258–272
ical modelling are correct, particularly given the importance of the localised redistribution of resources that govern the finer-scale biotic–abiotic interactions during runoff events. An attempt is made here to test the modelled spatial variations in hydrology within-plot using a qualitative approach. In the following analysis, only the runoff events with the highest NS efficiency values for each plot hydrograph are tested, since we are seeking to determine how successful spatial predictions of flow within the plot are for events that are simulated well at the plot outlet. The qualitative testing of the modelled spatial hydrology uses point based observations of hydrology, describing whether or not flow was monitored at 15 points within each plot using mini-samplers (Parsons et al., 2003). Since the other components of Mahleran are driven by hydrology, it is assumed that if the spatial simulation of hydrology is reasonable, then so too will the spatial
simulation of erosion, particulate-bound and dissolved nutrient fluxes. However, the need remains for the successful collection of spatial data on water, nutrient and sediment fluxes in order to test quantitatively how well the current level of process representation in Mahleran is able to simulate correctly the spatial fluxes. The modelled spatial dynamics of runoff over all the plots vary considerably. Contingency tables of spatial flow results are shown in Table 3. The j coefficients for the events tested show that the model produces varied results in predicting the spatial distribution in runoff. Overall, there was a tendency for the model to predict no flow at locations where flow was observed over the transition plots (plots 2 and 3), which indicates the need for improved data on model parameterisation for mixed grassland and shrubland hillslopes. The lowest j coefficient (0.09) was derived for Plot 2 on 15/08/2006, which, of all the events evaluated, had the lowest total
Table 4 Comparison of monitored and modelled erosion and particle-size distribution of eroded sediment. Spearman’s rank-correlation coefficients and Nash–Sutcliffe model-efficiency statistics of monitored versus modelled sediment yield are given for each plot. Plot/event date
Calibrated plot runoff (l)
Plot 1 07/09/05 29/09/05 05/07/06 28/07/06 01/08/06 11/08/06 29/08/06 07/09/06
Observed sediment yield (kg)
4.10 0.10 0.90 0.20 1.50 0.40 2.00 0.90
Proportions by model sediment size class Observed u=1
u=2
u=3
u=4
u=5
u=6
0.53 – 0.43 0.29 0.41 0.29 0.60 0.50
0.40 – 0.44 0.60 0.46 0.57 0.34 0.41
0.04 – 0.05 0.08 0.07 0.09 0.03 0.05
0.03 – 0.04 0.02 0.06 0.03 0.03 0.03
0.00 – 0.04 0.00 0.01 0.01 0.01 0.01
0.00 – 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.02 0.00 0.01
0.00 0.00 0.00 0.00 0.00 0.00
Modelled sediment yield (kg)
Proportions by model sediment size class Modelled u=1
u=2
u=3
u=4
u=5
u=6
0.29 0.38 0.42 0.69 0.41 0.47 0.26 0.31
0.68 0.62 0.58 0.31 0.58 0.53 0.73 0.68
0.02 0.00 0.00 0.00 0.00 0.00 0.01 0.01
0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
3.41 0.01 0.09 0.23 1.06 0.99
0.39 0.95 0.65 0.53 0.42 0.40
0.60 0.05 0.35 0.47 0.58 0.59
0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00
27.1 0.2 0.6 0.1 0.3 0.4 2.6 1.7
Summary statistics for all events: Spearman’s rank 0.826 sig. 0.011 NS efficiency 5.78 Plot 2 07/09/05 04/08/06 11/08/06 15/08/06 29/08/06 07/09/06
6.25 0.02 0.53 0.54 1.26 1.12
0.46 0.73 0.32 0.38 0.46 0.39
0.42 0.25 0.49 0.47 0.42 0.49
0.07 0.01 0.11 0.07 0.06 0.07
0.05 0.00 0.07 0.06 0.06 0.05
Summary statistics for all events: Spearman’s rank 1.000 sig.
