Modeling the effects of agricultural BMPs on sediments, nutrients, and water quality of the Beaurivage River watershed (Quebec, Canada)

Abstract

Agriculture has evolved into the largest non-point source of surface water pollution in Canada as a result of intensification over the past forty years. The Canadian WEBs project (Watershed Evaluation of Beneficial Management Practices, BMPs) was launched to evaluate the environmental and economic performance of BMPs as a means to mitigate agricultural sediment and nutrient issues. In this paper, the Gestion Intégrée des Bassins versant à l’aide d’un Système Informatisé (GIBSI) (or Integrated Watershed Management using a Computer System) integrated modeling framework was used to evaluate the effects of different BMPs on sediment and nutrient yields, as well as water quality in the Beaurivage River watershed in the province of Quebec. A reference scenario was developed that describes the current situation (i.e., base case scenario) of the watershed by calibrating the models used within GIBSI, namely HYDROTEL for hydrology, the Revised Universal Soil Loss Equation (RUSLE) for soil erosion, the Erosion-the Productivity Impact Calculator (EPIC) of the Soil and Water Assessment Tool (SWAT) for contaminant transport and fate, and QUAL2E for stream water quality. The effects of four BMPs were studied: (1) vegetated riparian buffer strips, (2) precision slurry application, (3) grassland conversion of cereal and corn fields, and (4) no-till (on corn fields). Simulation results indicate that BMPs can be effective in reducing nutrient and suspended sediment loads in both surface runoff and stream flow. More specifically, buffer strip and crop rotation showed better efficiency than hog-slurry management and no-till on corn BMPs. Moreover, results highlight the need for further investigation of sediment dynamics in the stream network as well as in the riparian buffer strips.

Au cours des dernières décennies, l’intensification de l’agriculture a été telle que cette activité est devenue la source de pollution diffuse des eaux de surface la plus importante au Canada. Le programme canadien d’Évaluation des pratiques de gestion bénéfiques (PGB) à l’échelle des bassins hydrographiques (EPBH) a été lancé pour analyser les performances environnementale et économique des PGB liées aux nutriments et sédiments agricoles. Dans cet article, nous utilisons le système de modélisation intégrée Gestion Intégrée des Bassins versant à l’aide d’un Système Informatisé (GIBSI) pour évaluer les effets de différentes PGB sur les charges de sédiments et nutriments aux champs et la qualité de l’eau du bassin de la rivière Beaurivage au Québec. Nous avons développé un scénario de base qui décrit la situation présente du bassin en calant les modèles hydrologique (HYDROTEL), d’érosion (Revised Universal Soil Loss Equation (RUSLE)), de transport des contaminants (Soil and Water Assessment Tool (SWAT)/Erosion-Productivity Impact Calculator (EPIC)), et de la qualité de l’eau (QUAL2E) du système GIBSI. Nous avons étudié les effets liés à quatre PGB: (1) bandes riveraines végétalisées, (2) application précise de lisier, (3) conversion en pâturage et foin de champs de maïs et céréales, et (4) semi-direct de maïs. Les résultats de simulation indiquent que les PGB peuvent être efficaces pour réduire les charges de sédiments et nutriments aux champs et en rivière. Plus précisément, la simulation des bandes riveraines ainsi que les rotations de cultures semblent plus efficaces que la gestion de l’épandage du lisier de porcs à l’aide de rampes à pendillards ou le semis-direct sur le maïs. De plus, ces PGB soulignent le besoin de se pencher davantage sur la dynamique des sédiments tant dans le réseau hydrographique que dans les bandes riveraines.

Introduction

Human activity such as intensive agriculture (crop and animal production), timber harvesting, production of hydroelectric power and wastewater treatment can affect flow paths (Bracken and Croke 2007) and water quality at the watershed scale. The many uses of water, including recreational activities such as swimming and water sports, drinking water, and water resources management for agriculture, can lead to conflicts between different users. To render water use more sustainable, it is important to analyze all factors that can potentially influence water quality in the watershed, and evaluate the environmental benefits and economic feasibility of alternative management practices (Yang et al. 2007). This is of great importance for agricultural activities in particular, in view of the need to feed an increasing world population with crops and livestock produced on a relatively small area, and extract drinking water from scarce water resources.

To improve the understanding of the effects of agricultural Beneficial Management Practices (BMPs) on water quality, the Watershed Evaluation of BMPs (WEBs) project was launched in 2004 (Agriculture and Agri-Food Canada 2007, 2011). This project revolves around the in-depth process analysis of nine watersheds in nine different Canadian provinces, and aims to evaluate the environmental effectiveness and economic feasibility of small-scale BMPs (watershed drainage area less than 10 km²) as recommended by the National Farm Stewardship Program (NFSP). As such, WEBs has objectives similar to those of the Conservation Effects Assessment Project (CEAP) in the United States (Duriancik et al. 2008). One of the ultimate goals of WEBs is to upscale, on a regional basis, findings at the micro-watershed scale to the regional watershed scale using an integrated hydrologic/economic modeling framework. The emerging framework will support the development of BMP adoption scenarios that can be compared with reference scenarios without BMPs.

This paper explores the regional effects of BMPs, and in particular their capacity to mitigate water quality issues related to sediment and nutrient loads in the Beaurivage River watershed in the province of Quebec, Canada (drainage area 718 km²), with respect to a reference scenario that represents the current situation of this watershed. The approach consists of: (1) developing a reference scenario for the Beaurivage River watershed and calibrating an integrated modeling system, and (2) distributed modeling of the effects of BMPs on nutrient and suspended sediment loads and concentrations in surface water using the Gestion Intégrée des Bassins versant à l’aide d’un Système Informatisé (GIBSI) (or Integrated Watershed Management using a Computer System) integrated modeling framework. The growing season is the primary concern, because it is assumed that this is the period of maximum water use and agricultural activity.

Material and methods

The GIBSI integrated modeling framework for simulating watershed processes

The GIBSI integrated modeling framework was used to simulate watershed processes and to evaluate the impact of BMPs on surface water quality (Villeneuve et al. 1998; Rousseau et al. 2000, 2005; Quilbé and Rousseau 2007). It was designed as a modeling and decision support system (DSS) for water management, and can also be used as an impact assessment tool for studying the effects of management scenarios on surface water quality (Quilbé and Rousseau 2007). Within the WEBs hydrologic and integrated modeling component, GIBSI was selected to characterize watershed processes and to examine water quality benefits of BMPs within the Beaurivage River watershed. This selection is based on numerous existing studies using GIBSI on the Chaudière River watershed which includes the Beaurivage River. Note that other pilot studies within the WEBs framework (e.g. Ahmad et al., 2011) have made use of Soil and Water Assessment Tool (SWAT), which is a modeling platform similar to GIBSI. It is noteworthy that GIBSI uses the pollutant transport and fate modules of Soil and Water Assessment Tool (SWAT) / Erosion-Productivity Impact Calculator (EPIC) and the sediment routing algorithm of a previous version of SWAT (Arnold and Williams 1995). Both watershed models (GIBSI and SWAT) simulate soil erosion using similar, revised versions of the Universal Soil Loss Equation (USLE) (Wischmeier and Smith 1978), and in-stream nutrient transport and water quality are based on the QUAL2E (Brown and Barnwell 1987) model. Finally, while SWAT simulates hydrology using the Soil Conservation Service curve number approach on hydrological response units (spatial units with homogeneous slope, land cover and soil characteristics), GIBSI is based on a distributed hydrological model, HYDROTEL (Fortin et al. 2001a), that simulates hydrological processes on relatively homogeneous hydrological units (elementary sub-watertersheds with a dominant soil texture, slope and relative percentages of land cover; see below for more details).

GIBSI model structure

A flow chart of the GIBSI integrated modeling system is provided in Figure 1. GIBSI includes (1) a MySQL database server (Dubois et al. 2005) and geographic information systems (GIS) algorithms based on Geo-graphic Resources Analysis Support System (GRASS) (GRASS Development Team 2008) for managing geographical data and corresponding attributes, as well as spatial input and output; (2) model components that allow for simulating watershed processes (including runoff, erosion and contaminant transport, water quality and biological integrity); (3) management tools (for providing information on reservoirs, land cover, agricultural management, livestock management, point sources of pollution, valuation of environmental benefits and costs of implementing BMPs), and (4) a graphical user interface.

Figure 1 GIBSI (Gestion Intégrée des Bassins versant à l’aide d’un Système Informatisé or Integrated Watershed Management using a Computer System) model structure (adapted from Quilbé and Rousseau, 2007). P is daily precipitation; T_min and T_max and are daily minimum and maximum temperatures, respectively; Q is daily stream flow; WQ is daily water quality; BI is the biological integrity using the Eastern Canada Diatom Index (ECDI); WQGs are water quality guidelines.

There are two types of spatial units for each of which state variables are calculated: (1) river segments, for which GIBSI simulates in-stream processes, and (2) relatively homogeneous hydrological units (RHHUs), which are the elementary land surface units made up of three soil layers for which GIBSI simulates soil water and contaminant processes such as mass balance, overland flow and erosion. RHHUs correspond to sub-watersheds that are each connected to one downstream river segment. As presented in Figure 1, GIBSI includes separate modeling components: (1) hydrology (HYDROTEL), (2) soil erosion (RUSLE), (3) nutrient cycling and transportation by overland flow (SWAT/EPIC), (4) stream water quality (QUAL2E) and sediment routing (SWAT), (5) pathogen transport, and (6) biological integrity.

Modeling components (5) and (6) account for the simulation of transport and fate of pathogens following manure application (e.g., E. coli) (Rogel et al. in preparation), and the evaluation of biological integrity, based on calculation of the Eastern Canadian Diatom Index (ECDI, also known as Indice Diatomées de l’Est du Canada) (Grenier et al. 2010a, 2010b, 2010c). These two components were not used in this study.

Simulations are performed at a daily time step and the meteorological input includes time series for daily precipitation, and minimum and maximum temperatures. GIBSI can generate output in the form of daily stream flow, water quality and maps that illustrate the probability of exceeding water quality guidelines, as well as comparative analyses, which are useful for locating stream segments with poor water quality and biological integrity.

The remainder of this section provides examples of studies for which GIBSI was used, and descriptions of the hydrological model, the phosphorus (P) and nitrogen (N) conversion and transportation model, the stream water quality model and the sediment routing model. For more detailed descriptions of the model components, refer to Villeneuve et al. (1998); more details on the soil database are provided in Duchemin et al. (2001).

Applied studies with GIBSI

GIBSI was developed and first applied to the Chaudière River watershed (6682 km²) in southern Quebec (Canada), and has been used to: (1) quantify the effects of a municipal clean water program on the water quality in the Chaudière River watershed (Mailhot et al. 2002); (2) predict the impact of deforestation on the hydrologic response of the 728-km² Famine River watershed (Lavigne et al. 2004); (3) develop a framework for estimating risk-based total maximum daily load (TMDL) (Rousseau et al. 2002); (4) illustrate the importance of evaluating costs and environmental benefits associated with nutrient management in the Beaurivage River watershed (Salvano et al. 2004, 2006); (5) evaluate the response of a watershed to climate change and possible future land covers (Quilbé et al. 2008a); (6) estimate the effect of climate change on stream flow (Quilbé et al. 2008b); (7) investigate the effects of historical land cover changes on runoff and low flows (Savary et al. 2009); (8) define for six Canadian watersheds the achievable performance standards for pesticides in response to different BMPs (Rousseau et al. 2012), and (9) predict the effect of BMPs on biological integrity in a small agricultural watershed (Grenier et al. 2010c). In the context of the WEBs project, GIBSI has also been used to evaluate economic models associated with the implementation of BMPs (Ghazalian et al. 2010).

Hydrological model (HYDROTEL)

HYDROTEL is a semi-distributed physically-based hydrological model compatible with GIS and remote sensing data (Fortin et al. 2001a, 2001b; Turcotte et al. 2003, 2007; Bouda et al. 2012). The model is currently used for inflow and hydrological forecasting by Hydro-Quebec, Quebec’s largest producer of electricity, and the Quebec Hydrological Expertise Centre which is in charge of the management and safety of publicly owned dams (Turcotte et al. 2004). HYDROTEL includes six computational modules that are run one after another: (1) a module for weather data interpolation based on either the weighted mean of the nearest three stations or Thiessen polygons; (2) a snow accumulation and snow melt module based on a mixed (degree-day) energy-budget method; (3) a potential evapotranspiration module based on the Thornthwaite equation, Linacre equation, Penman-Monteith equation, Priestley-Taylor equation, or Hydro-Québec equation; (4) a module for vertical water budget based on BV3C; (5) a module for overland flow in sub-watersheds towards river networks based on the kinematic wave equation, and (6) a channel or stream flow routing module based on either the kinematic wave equation or the diffusion wave equation.

Input data for simulations are maximum and minimum air temperatures and precipitation expressed in water equivalent, and the calculation interval can be set to three hours or one day. It is also noteworthy that the present study uses an empirical equation for potential evapotranspiration formulated by Hydro-Québec that requires only air temperature, and which has previously been applied to the Chaudière River watershed (Fortin et al. 2001a).

Soil erosion model (RUSLE)

Soil erosion (sediment yield) on cropland is calculated with the Revised Universal Soil Loss Equation (RUSLE) (Wischmeier and Smith 1978; Renard et al. 1997), as inherent sediment transport is simulated with the sediment transport equations of Yalin (1963). For this study RUSLE, an empirical equation, is applied to estimate average soil loss based on a multiplicative chain of factors including: (1) R (daily rainfall-runoff erosivity factor), (2) K (soil erodibility factor), (3) L (slope length factor), (4) S (slope steepness factor), (5) C (cover-management factor) and (6) P (support practice factor). Yalin’s equation evaluates the sediment transport capacity of the surface runoff shear stress applied on eroded sediments. The ensuing transport capacity is then compared to the estimated soil loss to determine the sediment yield to the stream network.

Model for P and N cycling and transport by overland flow (SWAT/EPIC)

Models for P and N cycling, pesticide fate, and transport and transformations from cropland are simulated with the contaminant transport algorithms of SWAT (Neitsch et al. 2005) and EPIC (Williams 1995; Arnold and Williams 1995; Izaurralde et al. 2006). These models account for the predominant processes of P and N cycling in the soil, while keeping the number of parameters to a justifiable minimum. A more detailed description of the P and N cycles as used in GIBSI can be found in Lasbleis et al. (2008).

Phosphorus cycling and transport

Five different P pools were simulated (Figure 2), two for organic P and three for inorganic P. Fresh organic P is associated with crop residues and microbial biomass while active organic P is associated with soil humus. Soil inorganic P is divided into a dissolved P pool, an active P pool, and a stable P pool. Dissolved P is in rapid equilibrium with active P, while active P is in slow equilibrium with stable P.

Figure 2 Phosphorus (P) pools and processes simulated in GIBSI (Gestion Intégrée des Bassins versant à l’aide d’un Système Informatisé or Integrated Watershed Management using a Computer System) (adapted from Neitsch et al. 2005).

Variables simulated for the P cycle include: (1) the initial distribution between the different chemical forms; (2) estimation of crop P demand, which is supplied by fertilizer; (3) decomposition and mineralization of fresh organic P; (4) mineralization of active organic P; (5) immobilization of active mineral P; (6) conversion between mineral pools; (7) leaching of dissolved mineral P, and (8) sorption to sediment particles and subsequent transport of P. P-related processes are most active in the topsoil, and therefore are simulated exclusively for the first soil layer in the model (0.1 m).

Nitrogen cycling and transport

GIBSI also simulates four N pools (Figure 3). Three pools represent organic N while one pool accounts for inorganic N (nitrate, NO₃ ^–). Fresh organic N is associated with crop residue and microbial biomass, while active and stable organic N are associated with humus in the topsoil. Organic N associated with humus is either active or stable, and both are available for mineralization.

Figure 3 Nitrogen- (N-) pools and processes simulated in GIBSI (Gestion Intégrée des Bassins versant à l’aide d’un Système Informatisé or Integrated Watershed Management using a Computer System) (adapted from Neitsch et al. 2005).

Variables simulated for the N cycle are: (1) initial distribution between the different chemical forms; (2) estimation of crop N demand, which is supplied by fertilizer; (3) N input from rainfall; (4) decomposition and mineralization of fresh organic N; (5) conversion between active and stable organic N; (6) mineralization of organic N; (7) immobilization of N; (8) denitrification and volatilization of N; (9) NO₃ losses by dissolution in surface runoff and lateral flow in the soil; and, (10) sorption to sediment particles and subsequent transport of N. Since N is more soluble and mobile than P, N cycling and transport are simulated for all three soil layers represented in GIBSI.

Processes that involve nitrates, nitrites and ammonium, such as volatilization, immobilization and denitrification, are all represented by the nitrate pool, hereby simplifying the simulation of conversion between different forms of N [NO₃ ^–, nitrite (NO₂ ^–), and ammonium (NH₄ ⁺)]. Because this simplification can lead to excessive losses by means of denitrification, volatilization and dissolution in surface runoff and lateral flow in the soil, a maximum allowed mineralization rate was defined to reduce this effect (Lasbleis et al. 2008).

Initial P and N levels are based on crop and soil characteristics observed in the Beaurivage River watershed (Beaudin 2006) and on typical values for crops that are produced in this watershed (Brady 1990). The P and N transport and conversion model also accounts for plant growth and nutrient uptake, in particular: (1) plant growth, (2) energy balance (heat units), (3) root development, (4) uptake of mineral P and N, (5) temperature stress, (6) water stress and (7) crop yield (Arnold and Williams 1995; Williams 1995).

Stream water quality (QUAL2E) and sediment routing model (SWAT)

Stream water quality is simulated with QUAL2E (Brown and Barnwell 1987) a one-dimensional steady-state, non-uniform flow water quality model which simulates the following variables: dissolved oxygen (DO) and re-aeration, biological oxygen demand (BOD), water temperature, P and N cycling, growth and fate of coliforms and chlorophyll-a in algae, dispersion and advection of up to three non-conservative constituents and up to two conservative constituents. Also, the sediment routing algorithm of SWAT (Arnold et al. 1995) was added to QUAL2E to account for stream deposition and bed degradation. QUAL2E governing equations are advection-dispersion-reaction equations with external sources and sinks. A flow balance is assumed, and steady non-uniform flow is used to solve the advection equation.

To discuss the model equations here would go beyond the scope of this paper. For more information on the models used in GIBSI refer to cited papers.

The Beaurivage River watershed

The Beaurivage River is an 87-km long tributary of the Chaudière River, draining an area of 718 km² (Figure 4) in southern Quebec (Canada). The climate is humid continental (cool summers and cold snowy winters) with an average annual precipitation of 1100 mm and a growing season of 120–150 days. Soil texture ranges from loamy sand to clay loam and the terrain is flat to gently sloping. Agricultural activity (accounting for 36% of the watershed area) is characterized by intensive livestock operations and dairy farming. In 1994, there were a total of 30,782 pigs or 42.9 pigs km^–2 in the watershed (Bédard et al. 1998). Agricultural land is composed (Rousseau et al. 2008) of pastures/hay (60%), cereal (23%) and cornfields (17%).

Figure 4 Location of the Beaurivage River watershed in the province of Quebec (Canada), the stream flow monitoring stations and water quality sampling sites (Quilbé et al. 2006).

Simoneau et al. (1998) observed P-limited eutrophication of the Beaurivage River, which they attributed for 86% to agricultural sources. Repeated applications of large quantities of organic fertilizer (manure) on cropland have caused P enrichment of the soil, thereby significantly increasing the risk of water contamination (Simard et al. 1995). In 2001, annual P and N application rates were estimated at 128 kg ha^–1 and 24.6 kg ha^–1, respectively, in excess of crop demand. The degradation of water quality can to some extent also be ascribed to point sources of pollution, notwithstanding the efforts made during recent decades to treat industrial and municipal wastewater (Mailhot et al. 2002). However, as Patoine et al. (2012) illustrated recently, seven of eight Quebec rivers in livestock areas have phosphorus levels that are at least twice as large as than the guideline for protection of rivers against eutrophication.

Reference scenario for the Beaurivage river watershed

Use of GIBSI to evaluate the effect of BMPs on nutrient and suspended sediment loads in the Beaurivage River watershed requires the definition of a reference scenario that describes the current situation of this watershed with respect to runoff, sediment and nutrient yield, water quality, reservoirs, land cover, agricultural management (tillage, pest control, nutrient management) and point sources of pollution. The following approach was adopted:

1.	Construction of a GIS database of the Beaurivage River watershed (PHYSITEL/GIBSI);
2.	Calibration of GIBSI’s hydrological model (HYDROTEL);
3.	Calibration of the model for nutrient cycling and transport by overland flow (SWAT/EPIC); and,
4.	Calibration of the stream water quality (QUAL2E) and sediment routing model (SWAT).

Building a spatial database of the watershed (PHYSITEL/GIBSI)

Database construction involved the development of relevant map layers of the watershed using PHYSITEL, a GIS specifically designed for this purpose (Turcotte et al. 2001; Rousseau et al. 2011), and subsequent integration of these map layers and attribute data into GIBSI.

Spatial subdivision of the Beaurivage River watershed (PHYSITEL)

Drainage area was determined from a 20-m digital elevation model (DEM) using the D8-LTD algorithm (Orlandini et al. 2003) in PHYSITEL (Turcotte et al. 2001; Rousseau et al. 2011). The algorithm first identified the grid cell with the greatest upstream drainage area (i.e. the watershed outlet), and subsequently the neighboring cells contributing to it. The watershed was found by repeating this procedure for every grid cell until the drainage divide was reached, and the stream network was generated based on all grid cells with an upstream contributing area greater than a given threshold. Depending on the threshold value selected by the user, the stream network may be more (small threshold area) or less (large threshold area) ramified.

A threshold area of 60 ha was selected for the reference scenario, resulting in a stream network of 609 river segments (average length 1.14 km) and 675 RHHUs (with an area of 105 ha on average). Hydrological processes were simulated for three soil layers: a first soil layer with a thickness of 0.1 m, a second layer with a thickness of 0.1 m and a third layer with a thickness of 1.3 m. Groundwater flow and transport processes below these layers (deeper than 1.5 m) were not accounted for.

Preparation of land cover and soil property layers (PHYSITEL)

Land cover was determined for each RHHU using a series of Landsat images (1976, 1981, 1987, 1990, 1995, 1999, 2003) (Savary et al. 2009) complemented with data from 2007 (provided by Financière Agricole du Québec). Hydrodynamic and chemical properties of the soil were also determined for each RHHU. An overlay of the digitized soil map (Duchemin et al. 2001) and RHHU map was used to calculate the average hydrodynamic and chemical properties of the soil, weighted according to the portion of the RHHU area intersecting with a given soil map unit. Note that RHHUs are elementary units: the approach used does not account explicitly for the presence of cropland, or for the spatial distribution of land cover and soil type within these elementary units.

Construction of the spatial database (GIBSI)

After construction of the geographic layers with the watershed boundary, DEM, slope, stream network, RHHU, land cover, hydrodynamic and chemical properties of the soil, and point data, they were subsequently imported in GIBSI. The land cover map, DEM and slope map were formatted as raster maps, the other layers as vector maps. Point data included the locations of weather stations, streamflow stations and water quality sampling sites, water purification and treatment plants, point sources of pollution, and water recreation sites divided into primary contact sites (swimming and kayaking), secondary contact sites (canoeing and fishing), and visual contact sites (hiking trails and parks along the riverside).

The spatial map data were complemented with attribute data on agricultural management and costs, livestock data, input weather data (daily precipitation and daily minimum and maximum temperatures) and, after completing a simulation, the model output (daily stream flow and water quality). An overview of data sources is given in Table 1, and more detailed information about the construction of the spatial database is provided by Villeneuve et al. (1998) and Rousseau et al. (2008).

Table 1. Sources of spatial data used for the reference scenario of the Beaurivage River watershed.

Data	Source
DEM (Digital Elevation Model)	CCOG GeoBase (www.geobase.ca)
Stream network	CCOG GeoBase (www.geobase.ca)
Land cover	Savary et al. (2009) (Landsat Image and Financière Agricole du Quebec)
Soil type	Duchemin et al. (2001a)
Weather	Environment Canada
Stream flow	Centre d’expertise hydrique du Québec, Quebec Government
Water quality	Ministère du Développement Durable de l’Environnement et des Parcs (MDDEP), Quebec Government
Point sources of pollution	Ministère des Affaires Municipales des Régions et de l’Occupation du Territoire (MAMROT), Quebec Government
Agricultural management practice and livestock counts (animal units per municipality)	Ministère de l’Agriculture des Pêcheries et de l’Alimentation, MAMROT, Quebec Government, Financière Agricole du Québec, Guide de Référence en Fertilisation (CRAAQ 2003).

Results and discussion

Calibration of GIBSI for the reference scenario

Calibration and validation of the hydrological model (HYDROTEL)

HYDROTEL was calibrated manually according to the approach developed by Fortin et al. (2001b) and Turcotte et al. (2003) for the Chaudière River and Beaurivage River watersheds. Daily simulated stream flow was calibrated using daily stream flow measured at the Saint-Étienne-de-Lauzon automated gauging station 5 km upstream of the watershed outlet during 1984–1989. Subsequent validation involved the 1989–1994 period.

Evaluation of simulation performance was based on: (1) the Nash-Sutcliffe (NS) efficiency coefficient (Nash and Sutcliffe 1970), (2) visual comparison of simulated and observed stream flow and (3) the coefficient of determination (R²), which describes the fit between simulated and observed stream flow. NS ranges from –∞ to 1, with NS = 1 representing an optimal fit. It is generally assumed that NS > 0.65 represents a satisfactory model performance for daily stream flow (Moriasi et al. 2007). The determination coefficient (R²) ranges between 0–1, with 1 representing an optimal fit.

Figure 5 shows the calibration results for a part of the calibration period (hydrological years 1984–1985 and 1988–1989). HYDROTEL performed relatively well with an NS efficiency coefficient between 0.70 and 0.87 (Table 2), although high flows during springtime were underestimated. This can be explained by an underestimation of frozen precipitation leading to underestimation of snow melt. Some storms are not simulated at all due to the fact that input data originates from weather stations outside the watershed, and may therefore not capture all storms. Similar behavior was observed for the period of validation, of which a part of the results is shown in Figure 6 (hydrological years 1991–1992 and 1992–1993). The model performance for the entire period of validation is given in Table 3.

Figure 5 Simulated and measured stream flow for the period of calibration: (a) hydrological years 1984–1985 and (b) 1988–1989.

Table 2. HYDROTEL performance for the period of calibration (1984–1989).

Calibration period (Hydrological year running from 1 October – 30 September)
Nash-Sutcliffe coefficient	1984–1985	1985–1986	1986–1987	1987–1988	1988–1989	Entire period
	0.87	0.77	0.70	0.76	0.82	0.79

Figure 6 Simulated and measured stream flow for the period of validation: (a) hydrological years 1991–1992 and (b) 1992–1993.

Table 3. HYDROTEL performance for the period of validation.

Period of validation (Hydrological year running from 1 October – 30 September)
Nash-Sutcliffe coefficient	1989–1990	1990–1991	1991–1992	1992–1993	1993–1994	Entire period
	0.51	0.75	0.70	0.80	0.84	0.75

Daily model performance was also evaluated for a longer period (1979–2004), in particular for the low-flow period between 1 June and 30 September (Table 4), the effect of BMPs are discussed in a later section. Although the simulated mean annual stream flow was 100 mm more than that measured, model performance was comparable to the calibration and validation results as discussed in the previous paragraph. In order to optimize the Nash-Sutcliffe coefficient, the underestimated stream flow for the spring season was compensated by an overestimation for low-flow periods. A statistical analysis showed that over 50% of the prediction errors for simulated stream flow on an annual basis were within a range of +/– 35% of the measured stream flow and centered on zero (non-biased).

Table 4. Long-term performance of HYDROTEL (1979–2004).

Performance criterion	Annual	Low-flow periods (1 June – 30 September)
Coefficient of determination (R^2)	0.77	0.77
Nash-Sutcliffe coefficient	0.77	0.76
Measured annual runoff (mm)	688	123
Simulated annual runoff (mm)	588	142

Calibration of the soil erosion model (RUSLE)

Because distributed data for suspended sediments are very scarce for the Beaurivage River watershed, simulated annual sediment yield on cropland was based on historical weather data and land cover images (1979 to 2004). The mean annual simulated soil loss for the Beaurivage River watershed during this period was 1.5 t ha^–1 yr^–1. This sediment yield reflects the relatively flat topography and relatively low percentage of agricultural land (36% of the total watershed area), and is less than the acceptable limit for agricultural soil loss in North America which is set at 6 t ha^–1 yr^–1 in (Wall et al. 2002), and between the 2.5 t ha^–1 yr^–1 calculated by (Bédard et al. 1998) using the USLE and 0.3 t ha^–1 yr^–1 estimated by Quilbé et al. (2006) for the same watershed. The order of magnitude is equal to that of studies on similar watersheds in the province of Quebec (Gangbazo and Babin 2000).

Calibration of the nutrient cycling and transport model (SWAT/EPIC)

The nutrient cycling and transport model was calibrated by simulating P and N pools for each agricultural RHHU, using weather data from the period 2000–2004. The amount of applied nutrients was based on crop N and P demand: 40 kg N ha^–1 and 22 kg P ha^–1 for grassland, 50 kg N ha^–1 and 30 kg P ha^–1 for cereal fields, and 150 kg N ha^–1 and 45 kg P ha^–1 for corn fields (Centre de Référence en Agriculture et Agroalimentaire du Québec [CRAAQ] 2003). Crop rotation was based on typical rotation schemes for the Beaurivage River watershed, such as hay, alternated with either cereals or corn. N input from rainfall was based on reference values for temperate climates: between 5 and 10 kg ha^–1 yr^–1 (Brady 1990).

Phosphorus cycle

Figure 7 shows the simulated dynamics of P pools and P uptake for all agricultural RHHUs combined over the period 2000–2004 (5 years). Total P was at all times greater than 275 kg ha^–1 (ha of agricultural land), which is less than values commonly observed for pasture, cereal and corn (500 to 2000 kg ha^–1, Brady 1990). Variations in stable mineral P are related to the balance between P removal and conversion between stable and active mineral P. P uptake, which is assumed equal to crop P demand, showed annual variations corresponding to the growing season, starting with a steady increase during the initial growth stage in early May until reaching a maximum uptake rate in the mature stage. Overall, the simulation of P cycling produced coherent results with values within plausible ranges.

Figure 7 SWAT (Soil and Water Assessment Tool) /EPIC (Erosion-Productivity Impact Calculator) simulation of P pools for 2000–2004 (5 years).

Phosporus cycling was simulated for an extended period, again using land cover images for 1979–2004 as was done for the calibration of the soil erosion model. The mean annual total P loss was 3.2 kg ha^–1 yr^–1, of which 93% was mineral P. Losses with the same order of magnitude were found by Gangbazo and Babin (2000) (0.4 kg ha^–1 yr^–1 to 2.2 kg ha^–1 yr^–1) and Quilbé et al. (2006) (1.0 kg ha^–1 yr^–1).

Nitrogen cycle

The simulated evolution of N pools and N uptake for all agricultural RHHUs combined over the period 2000–2004 (5 years) is shown in Figure 8. Total N was over 1450 kg ha^–1, which is within the range of values (1000 to 5000 kg ha^–1, Brady 1990) observed for pasture, cereal and corn, and stable organic N was between 1050 and 1150 kg ha^–1.

Figure 8 Soil and Water Assessment Tool (SWAT)/Erosion-Productivity Impact Calculator (EPIC) simulation of N pools for 2000–2004 (5 years).

As was the case for P, the growing season can be recognized in the evolution of simulated N uptake (calculated as crop N demand) and fresh organic N, as well as in the sudden decrease of nitrates due to uptake during the initial growth stage. N uptake reached a maximum rate in the mature stage, similarly to P uptake. Small variations in N uptake are explained by changes in temperature and water stress. The amount of mineral N remained constant at around 80 kg ha^–1 yr^–1, which represents 6% of the amount of total organic N, a value that compares well to values found by Brady (1990). Overall, the simulation of N cycling produced coherent results with values within plausible ranges.

The N cycling simulation performed for the extended period 1979–2004, again using land cover images for that period, resulted in a mean annual simulated total N loss of 27 kg ha^–1 yr^–1, which represents nearly 2% of total N in the soil, a value comparable to those found by Brady (1990). Total N remained within the same order of magnitude as values presented by Gangbazo and Babin (2000) (4.7kg ha^–1 yr^–1 to 23.9 kg ha^–1 yr^–1) and Quilbé et al. (2006) (8.5 kg ha^–1 yr^–1). The simulated mean annual loss of nitrates (NO₃ ^–) was equal to 17 kg ha^–1 yr^–1, which represents over 63% of total N loss.

Calibration of the stream water quality (QUAL2E) and sediment routing models (SWAT)

QUAL2E was calibrated for the period 1988–2004 for concentrations and loads of in-stream P, in-stream N, and suspended sediments. Observed loads were computed with measured stream flows at the Saint-Étienne-de-Lauzon automated gauging station (located 5 km upstream of the watershed outlet) and water samples were taken at a station located 1.1 km upstream of the stream flow gauging station. Water quality was sampled once per week from 1989–1990, twice per month from 1991–1995, and once or twice per month from 1996–2004, and involved total suspended sediment, total P, nitrate-nitrogen (NO₃ ^–-N) and nitrite-nitrogen (NO₂ ^–-N); the latter two were analyzed jointly.

Irregular sampling rates present difficulties in the estimation of annual contaminant loads at the watershed scale, (e.g., Moatar et al. 2006; Quilbé et al. 2006; Mailhot et al. 2008), which is why daily constituent concentrations and loads were estimated with LOADEST (Crawford 1991, 1996; Runkel et al. 2004). This computer program fits a regression model to simulated concentrations based on stream flow data and observed constituent concentrations. Note that the daily comparisons remain limited because observed concentrations reflect only instant conditions and simulated concentrations with the QUAL2E model representing daily mean concentrations. With such scarce measurements, daily concentrations and loads were statistically grouped using monthly boxplots. This approach is consistent with the observations made by Borah and Bera (2004) and Borah et al. (2006) that daily concentrations simulated by watershed models might be more difficult to reproduce than monthly observations. Simulations of stream water quality are based on the assumption that pollution by sewer overflows is negligible and wastewater treatment plants produce a constant output rate of contaminants.

Suspended sediments

Figure 9 shows the monthly distributions of simulated (QUAL2E), estimated (LOADEST) and measured daily suspended sediment concentrations and loads. Suspended sediment concentrations and loads were at their lowest during winter when snowfall prevailed, and increased rapidly with the snowmelt in spring.

Figure 9 Monthly distributions of simulated (n = 6940), estimated (n = 6118) and measured (n = 345) daily suspended sediment concentrations (a) and loads (b). The boxes correspond with the 2^nd and 3^rd quartiles and the whiskers represent the 10^th and 90^th percentiles.

A Wilcoxon signed-rank test was performed to provide a measure of agreement between monthly distributions of simulated and measured sediment concentrations. This test showed that the results are satisfactory: median simulated suspended sediment concentrations are not significantly different from median observed suspended sediment concentrations (p < 0.05), except for periods with enhanced flow in spring (April, May) and fall (October, November). Similar results were found for the comparison of median simulated loads with median observed loads, although model performance for sediment concentration or loads is similar. Because stream flow is one of the input variables of QUAL2E, errors in simulated runoff partly explain the errors produced by the water quality model. It is difficult to correct for these errors due to the small number of sampling sites for water quality and their low sampling frequencies. The mean in-stream suspended sediment simulated for the Saint-Étienne-de-Lauzon water sampling site was equal to 0.9 t ha^–1 yr^–1 (per ha of agricultural land), which represents 60% of the simulated overland sediment yield for the watershed (i.e., a sediment delivery ratio of 60%).

In-stream phosphorus

Monthly distributions of simulated (QUAL2E), estimated (LOADEST) and measured daily total phosphorus (TP) concentrations and loads are plotted in Figure 10. Again, concentrations and loads were at their lowest during winter, reaching a peak simultaneously with the runoff peak in the spring.

Figure 10 Monthly distributions of simulated, estimated and measured daily total phosphorus (TP) concentrations (a) and loads (b). The boxes correspond with the 2^nd and 3^rd quartiles and the whiskers represent the 10^th and 90^th percentiles.

The Wilcoxon signed-rank test showed that median simulated TP concentrations are not significantly different from median observed TP concentrations for the period March to September (p < 0.05), meaning that the model performs quite well for most of the year. The largest model errors occurred for the period December–February, which can partly be ascribed to (1) the one-layer representation of P-related processes accounting for only the first 0.1 m below the soil surface, and (2) the combination of very low TP concentrations (compared to in-stream nitrogen and suspended sediments) and the presence of an ice cover that affects stream flow measurements. Model errors during summer were mostly the result of incomplete knowledge of fertilizer application schedules. Overall, model performance for TP concentration and loads remain similar. Mean in-stream P simulated for the Saint-Étienne-de-Lauzon water sampling site was equal to 2.7 kg ha^–1 yr^–1 (per ha of agricultural land), which represents 80% of the simulated overland TP yield for the watershed (i.e., a P delivery ratio of 80%).

In-stream nitrogen

Figure 11 gives the monthly distributions of simulated (QUAL2E), estimated (LOADEST) and measured daily concentrations and loads of NO₂ ^–-N and NO₃ ^–-N combined. Overall, measurement trends or annual tendencies are generally reproduced by simulated values with lower values throughout the summer and larger values for the remaining years.

Figure 11 Monthly distributions of simulated, estimated and measured daily concentrations (a) and loads (b) of NO₂ ^–-N and NO₃ ^–-N combined. The boxes correspond with the 2^nd and 3^rd quartiles and the whiskers represent the 10^th and 90^th percentiles.

The Wilcoxon signed-rank test showed that median simulated concentrations and loads are not significantly different from median observed concentrations and loads for the period March to September (p < 0.05). NO₂ ^–-N and NO₃ ^–-N peaks in fall (October–November) were slightly overestimated, whereas winter peaks (January–February) were slightly underestimated. These model errors are believed to be mostly related to the fact that only NO₂ ^–-N and NO₃ ^–-N were simulated by the model used for nutrient cycling and transportation by overland flow (SWAT/EPIC), whereas other chemical forms such as NH₄ ⁺ are also present in the soil.

Mean in-stream NO₂ ^–-N and NO₃ ^–-N simulated for the Saint-Étienne-de-Lauzon water sampling site was equal to 16 kg ha^–1 yr^–1 (per ha of agricultural land). This value represents approximately 60% of the simulated mean annual overland load of total N (27 kg ha^–1 yr^–1).

Developing scenarios of beneficial management practices

Four scenarios with four different BMPs were developed, and GIBSI was subsequently used to evaluate their influence on loads of nutrients and suspended sediments in the Beaurivage River watershed.

Description of BMP scenarios

The following BMP scenarios were generated (the area affected by these BMP scenarios is given in Table 5):

1.

Riparian buffer strips – Transport of suspended load and contaminants from farmland to rivers is enhanced by steep ditch banks, annual row cropping, and lack of erosion control measures. Vegetated riparian buffer strips (also called vegetated filter strips) tend to reduce the amount of sediments and nutrients transported to streams by surface runoff (Lowrance et al. 1984; Correll 2005; Lowrance and Sheridan 2005; Cho et al. 2010). Quebec legislation presently requires a 3-m buffer strip for farmland (Ministère du Développement Durable, de l’Environnement et des Parcs [MDDEP] 2007). In this scenario 1-, 3-, and 5-m wide buffer strips on all agricultural RHHUs were accounted for and the trapping efficiency was evaluated as (Neitsch et al. 2005):

(1)

where T_eff is the fraction of the constituent loading trapped by the buffer strip and W_BS is the width of the buffer strip (m). This equation is used in SWAT and has previously been applied by Cho et al. (2010) and Volk et al. (2009). The trapping efficiency values calculated for 1, 3, and 5 m wide buffer strips equal 0.37, 0.51, and 0.59, respectively.

Table 5. Absolute (ha) and relative area (% of agricultural land) affected by the implementation of Beneficial Management Practices (BMP) scenarios in the Beaurivage River watershed.

BMP scenario	Area affected by implementation
	ha	Percent of agricultural land^∗
Riparian buffer strips (1, 3, 5 m)	25.77	100%
Precision slurry application	25.77	100%
Grassland conversion (of cereal and corn fields)	10.20	40%
No-till (corn fields)	4.32	17%
^∗Agricultural land represents 36% of the total watershed area.

2.	Precision slurry application – Slurry applicators with a trailing shoe mounted on a slurry tanker can apply pig slurry (a mixture of manure and water) directly onto the ground surface through a series of hoses (also known as low-ramp application). Precision slurry application leaves the slurry in rows spaced 15–20 cm apart and tends to reduce N losses through volatilization with respect to conventional surface spreading with, for example, a splash-plate system (Smith et al. 2000). For this BMP, GIBSI’s livestock management tool is used to specify the manure application practice and switch from conventional to surface application (the latter corresponding to a low-ramp with trail hoses application). Such change has one specific impact on the volatilization rate of N from 1.40 to 1.30 for conventional to surface (low-ramp) application (CRAAQ 2003). This difference will induce a reduced loss of N to the atmosphere and might increase N yields to the stream network. In both the slurry management scenario and the reference scenario, manure application corresponded with the crop required P for optimal maturation, and is representative of present-day crop production in the Beaurivage River watershed.
3.	Grassland conversion (of cereal and corn fields) – Planting perennial crops by alternating corn and cereals can help to reduce total P and total N, whereas grassland reduces overland flow and erosion altogether (Sood and Ritter 2010). A grassland scenario was simulated by converting all corn and cereal fields to hay and pasture (grassland conversion) with the land cover management tool. This GIBSI tool allows the user to first choose any land cover (corn and cereal) within a specific area (entire Beaurivage watershed) and then convert it to the desired land cover. This drastic BMP scenario reduces the application of fertilizer because grassland requires a smaller amount than cereal and corn. This BMP represents a radical change in economic and environmental conditions.
4.	No-till (corn fields) – To protect soil and water resources and save time and energy costs, many farmers have already adopted reduced tillage or no-till systems (Burgess et al. 1999). In a typical no-tillage system, corn is planted directly into the soil with no tillage between harvests, and crop residues are left on the ground as a protection against soil erosion. This practice allows the soil structure to re-establish itself and reduces sediment yield (Altieri 1995). No-till is implemented by reducing the C factor in the RUSLE (the soil erosion model), which is the cropping factor that describes the growth period. This BMP was applied to corn fields only.

Approach for the evaluation of BMP scenarios

Rousseau et al. (2002) proposed a methodology to evaluate the effect of BMPs on loads of nutrients and suspended sediments. Accordingly, each BMP scenario was assessed using a long-term simulation period; in this study, the period 1970–2006. Results were evaluated for 1 May to 31 October of each simulation year because this is the period of maximum water use and agricultural activity. The regional effect of the BMPs on water quality was evaluated at the outlet of the Beaurivage River watershed and compared with the reference scenario. The in-stream effects of each BMP scenario were expressed in terms of the probability that the nutrient and suspended sediment concentrations exceed the water quality guidelines (WQGs) listed hereafter (MDDEP 2008):

•	Suspended sediment concentration (SSC): The guideline for the protection of aquatic life allows an average of 5 mg/L above the natural background concentration in a given stretch of a river, which was estimated as the mean of all available measurements near the outlet of the Beaurivage River watershed. The SSC guideline was subsequently established at 70 mg/L.
•	Total phosphorus (TP): The guideline for the protection of aquatic life, recreational activities and aesthetics limits TP to 0.03 mg/L (MDDEP 2008) in order to prevent excessive growth of algae and aquatic plants in rivers and streams.
•	Nitrogen (NO3 –-N and NO2 –-N): The guideline for preventing the contamination of water and aquatic organisms is 10 mg/L (MDDEP 2008).
•	Nitrites (NO2 –-N): The guideline for the protection of aquatic life is 0.02 mg/L (MDDEP 2008).

Although no guidelines exist for water quality of surface runoff, evaluation of overland yields may provide to be very useful for a better understanding of the effects of BMPs. Therefore, the loads of total P, N and suspended sediment of surface runoff generated on RHHUs (the elementary land surface units in GIBSI) were also analyzed.

Simulating the effect of BMPs in the Beaurivage River watershed

Effect of BMPs on overland loads

Table 6 shows for different BMP scenarios the reduction in overland loads with respect to the watershed reference scenario during the period between 1 May and 31 October of the years 1970–2006. Simulations of scenarios with riparian buffer strips show a reduction in overland loads of TP, TN and suspended sediments. Reductions in TN and TP are smaller than those for sediments because these nutrients are not solely adsorbed to sediments but are also transported by soil water and not just surface runoff. The precision slurry application scenario resulted in an increase of TN load, whereas the grassland conversion scenario led to a considerable decrease in sediment yield but not quite as successful in trapping TP and TN compared to a 3-m wide riparian buffer strip. Finally, the no-till scenario for cereal fields produced the lowest reduction of overland loads because only 17% of the agricultural land was affected. Note the two dashes represent no impact of the BMP.

Table 6. Effects of Beneficial Management Practices (BMP) scenarios on surface runoff loads of nutrients and suspended sediment in the Beaurivage River watershed, simulated for the period 1970–2006 (1 May to 31 October of each year). TP, total phosphorus; TN, total nitrogen.

BMP scenario	Per cent reduction of surface runoff load
	Suspended sediment	TP	TN
Riparian buffer strip (1 m)	37%	28%	30%
Riparian buffer strip (3 m)	51%	39%	37%
Riparian buffer strip (5 m)	59%	46%	40%
Precision slurry application	-	-	+3.3% (increase)
Grassland conversion (of cereal and corn fields)	58%	34%	13%
No-till (corn fields)	8%	3%	1%

Effect of BMPs on in-stream loads

Based on available measurements, the effects of BMPs on the probability of exceeding water quality guidelines was evaluated on the basis of daily average in-stream loads of TP, NO₂ ^–-N and NO₃ ^–-N combined, NO₂ ^–-N only, and suspended sediments. The three scenarios with riparian buffer strips reduced in-stream loads and probabilities of exceeding water quality guidelines (Table 7). Reduction of in-stream load was greatest for TP and combined NO₂ ^–-N and NO₃ ^–-N, while the probability of exceeding water quality guidelines for combined NO₂ ^–-N and NO₃ ^–-N concentrations remained zero (i.e., NO₂ ^–-N and NO₃ ^–-N concentrations were always lower than the guideline). Reduction of loads and probabilities of exceedance were directly related to the width of riparian buffer strips and in agreement with the effect of riparian buffer strips on overland loads (Table 6). Conversely, no-till and precision slurry application had little effect on the probability of exceeding water quality guidelines at the watershed outlet. Precision slurry application led to an increase of in-stream combined NO₂ ^–-N and NO₃ ^–-N loads, while grassland conversion resulted in the greatest reduction of in-stream sediment load with a reduction of in-stream TP comparable to the TP reduction by a 3-m riparian buffer strip. The effect of grassland conversion on N loads was minimal.

Table 7. Effects of Beneficial Management Practices (BMP) scenarios on the probability of exceeding water quality guidelines (WQGs) and percent reduction of in-stream load in the Beaurivage River watershed, simulated for the period 1970–2006 (1 May to 31 October of each year). A probability of 1 means that WQGs are exceeded during all 184 simulated days of the year.

Probability of exceeding of WQGs and percent reduction of in-stream load
Scenario	Suspended sediments	TP	NO2 ^–-N and NO3 ^–-N	NO2 ^–-N
	> 70mg/L	> 0.03mg/L	> 10mg/L	> 0.02mg/L
Reference	0.37	0.91	0	0.47
Riparian buffer strip (1 m)	0.31 (10%)	0.87 (28%)	0 (34%)	0.46 (4%)
Riparian buffer strip (3 m)	0.30 (13%)	0.85 (39%)	0 (41%)	0.45 (6%)
Riparian buffer strip (5 m)	0.28 (15%)	0.83 (46%)	0 (43%)	0.44 (7%)
Precision slurry application	-	-	0 (+3.9%)	0.47 (+0.04%)
Grassland conversion (of cereal and corn fields)	0.27 (16%)	0.88 (32%)	0 (11%)	0.46 (2%)
No-till (corn fields)	0.37 (0%)	0.91 (3%)	0 (0.02%)	0.47 (0.5%)
TP, total phosphorus; NO2 ^–-N, nitrite-nitrogen; NO3 ^–-N, nitrate-nitrogen.

General outcome and future challenges

Given the satisfactory calibration results, namely daily stream flow by the hydrological model, annual yield of pollutant at the watershed scale compared to similar studies on this watershed, and monthly in-stream pollutants loads and concentrations for the crop season, the distributed aspect of the model can be used to investigate the impacts of BMPs on sediments, P and N. Indeed, distributed models allow the user to obtain results on diverse spatial scales instead of analyzing simulation results solely at the watershed outlet.

In all scenarios, the reduction of in-stream loads by BMPs were significant for suspended sediments (up to 16%), TP (up to 46%), and combined NO₂ ^–-N and NO₃ ^–-N (up to 43%). Although reductions of overland loads were even greater (up to 59%), none of the BMP scenarios was able to significantly reduce the probability of exceeding water quality guidelines in the Beaurivage River and tributaries (see Tables 6 and 7). Note that the utilization of the empirical equation of SWAT for buffer strip efficiency concurs with other studies (e.g., Cho et al. 2010). Even though such an equation does not bear any physical meaning, the relationship between this representation of the buffer strips’ efficiency based on width and their impact on pollutant loading may provide useful guidance for designing optimum riparian buffers within the watershed (Cho et al. 2010).

In a broader context, there are still challenges that need to be addressed by the modeler, as the representation of BMPs in the model should reflect field measurements where such practices are being implemented on the studied watersheds.

Effect of BMPs on suspended sediment loads

Even with a significant reduction of suspended sediment load, water quality guidelines could not be met. This is illustrated by simulations with 1-, 3-, and 5-m riparian buffer strips, which reduce the number of days during which water quality guidelines are exceeded by 11, 13 and 17 days, respectively, where the reference scenario that represents the current situation of the Beaurivage River watershed counted 68 days of exceedance (note that the analysis includes simulations of the period May 1 to October 31 – 184 days – of the years 1970–2006).

The grassland conversion scenario reduced the number of days when water quality guidelines were exceeded by 18 days. Again the reduction of in-stream loads (16%) was less than the reduction of surface runoff loads (58%), which can be explained by the fact that the elementary land surface units (RHHUs) for which overland loads were simulated are also the source areas of these loads: after erosion and transport, the heaviest load is rapidly deposited when flowing over grassland due to the drag force of the grass blades that slows down surface runoff. In the Ouachita River Basin, in the northern part of Louisiana, the 1480-km² Bayou LaFourche watershed recorded much better results with a 97% sediment load reduction (Yuan et al. 2006). This other grassland scenario modeled with the United States Department of Agriculture (USDA) Annualized Agricultural Non-Point Source Pollution model (AnnAGNPS) is taking place in different conditions, despite a similar part of the watershed used for agriculture (41% of cotton and soybeans mostly for the Bayou LaFourche watershed versus 36% of pasture/hay cereals and corn for the Beaurivage River watershed) and comparable slope conditions. This may suggest an underestimation of the grassland scenario results simulated with GIBSI. Moreover, instead of converting all cropland to grass, the modeler could select the highest eroded croplands as carried out in Jha et al. (2010).

Finally, the no-till scenario has the lowest effect on in-stream sediment, despite a slight decrease (4%) in overland loads from corn fields. A decrease of this order is insufficient to affect the in-stream sediment load at the outlet of the Beaurivage River watershed and matches, in that manner, the conclusion of Moore et al. (1992).

Differences in the reduction of overland loads and in-stream loads at the watershed outlet demonstrate that the GIBSI sediment routing model (SWAT) tends to simulate stream deposition of suspended sediment during most of the year and re-suspension mostly during heavy rainfall episodes with overland flow. For this reason, 95% of the suspended sediments simulated in the river segments near the watershed outlet were re-suspended. The stream network acts as a sediment buffer that is able to store and release sediments, but the rates at which storage and release occur are not necessarily proportional to the mitigating effects of upstream BMPs on in-stream loads.

Effect of BMPs on phosphorus loads

The presence of 1-, 3- and 5-m riparian buffer strips reduced the number of days during which the TP guideline was exceeded in the Beaurivage River watershed by 7, 11 and 15 days, respectively (note again that the analysis includes simulations of the period 1 May to 31 October – 184 days – of the years 1970–2006), compared to the reference scenario with 167 days of exceedance. Contrary to suspended sediments, the corresponding reductions of TP load were the same for surface runoff and stream flow (28%, 39% and 46% for buffer widths of 1, 3 and 5 m, respectively). Similar behaviors between overland loads and in-stream loads were also observed for the scenarios of grassland conversion and no-till, and can be explained by the relatively conservative nature of P.

The grassland conversion scenario was able to reduce the number of days during which the TP guideline was exceeded by 6 days. Here too, the corresponding simulated reduction of TP load was similar for surface runoff and stream flow (34% and 32%, respectively), and this reduction can be explained by the smaller amount of TP applied to grassland as fertilizer (22 kg ha^–1) compared to what was applied to cereal fields (30 kg ha^–1) and corn fields (45 kg ha^–1). No-till had hardly any effect on in-stream and overland loads of TP (3% decrease), again because only 17% of the agricultural land was affected. Sood and Ritter (2010) showed greater impact of grassland conversion on P than no-till on all farm land.

In addition to the agricultural non-point sources of P, the influence of point sources of P associated with wastewater treatment plants in the Beaurivage River watershed were estimated by developing an additional scenario for GIBSI with point sources only (agricultural non-point sources of P were not accounted for). A comparison with the reference scenario demonstrated that on average 25% of TP concentrations and 2.3% of TP loads simulated for the watershed outlet could be traced back to wastewater treatment plants. Point sources of P alone accounted for a probability of 0.12 of exceeding water quality guidelines with respect to TP.

Effect of BMPs on nitrogen loads

The probability of exceeding the guideline for combined NO₂ ^–-N and NO₃ ^–-N load was not affected by BMPs because loads were already below the critical value in the reference scenario. A minor improvement was simulated for NO₂ ^–-N-only loads: the presence of riparian buffer strips of 1, 3 and 5 m reduced the number of days during which the nitrite guideline was exceeded (86 days in the reference scenario) by 2, 4 and 6 days, respectively. The reductions of in-stream combined NO₂ ^–-N and NO₃ ^–-N loads (34%, 41% and 43%, respectively) were larger than the reduction of overland loads of TN (30%, 37% and 40%, respectively), while the reductions of in-stream nitrite loads (4%, 6% and 7%) were smaller. This pattern is explained by two factors: firstly, TN reduction already includes combined nitrate-nitrite reduction, which in turn includes nitrite reduction; secondly, the fate of N (all forms) in the stream network is simulated by QUAL2E using nonlinear relations for transport and conversion.

The grassland conversion scenario did not affect in-stream nitrite load, while the reduction of in-stream combined nitrate-nitrite load (11%) was similar to that of overland load of TN (13%) because the amount of manure applied to grassland (40 kg ha^–1 for TN) is less than for cereal fields (50 kg ha^–1 for TN) and corn fields (150 kg ha^–1 for TN).

The no-till scenario had almost no effect on the probability of exceeding the water quality guideline for nitrites. The reduction of in-stream loads (0.5%) was similar to that of overland load of TN on corn fields (1%). The precision slurry application scenario resulted in a slight increase of in-stream N loads at the outlet of the Beaurivage River watershed, which can be ascribed to reduced volatilization [ammonia (NH₃) losses] during and after low-ramp slurry application and the subsequent increase of available N in the soil and in surface runoff (Smith et al. 2000). Also Burgess et al. (1999) suggested that tillage system generally had relatively little effect on soil NO₃ ^– under till and no-till studied condition.

Figure 12 shows the probability maps of exceeding the TP guideline of 0.03 mg/L for the reference scenario and for the BMP scenario with 5-m riparian buffer strips. This type of map, produced by a built-in feature of GIBSI, represents a powerful tool in the evaluation of the level of agricultural stress in a particular region and the effect of BMPs on water quality. This type of map can be ultimately shared with stakeholders in order to evaluate the advantages/drawbacks of an agricultural practice.

Figure 12 Maps illustrating the probability of exceeding the TP guideline of 0.03 mg/L for the reference scenario (a) and for the Beneficial Management Practices (BMP) scenario with 5-m riparian buffer strips (b).

Conclusion

The objectives of the integrated modeling component of the WEBs project and of this paper are to evaluate the regional effects of BMPs, and their ability to mitigate water quality issues related to suspended sediments and nutrient loads. The GIBSI integrated modeling framework was used in a case study on the Beaurivage River watershed in Quebec, Canada, in which a reference scenario was developed and subsequently used to evaluate water quality changes during the growing season following a 36-year simulation of four scenarios with different BMPs: (1) riparian buffer strips (1, 3 and 5 m wide), (2) precision slurry application, (3) grassland conversion of cereal and corn fields, and (4) no-till on corn fields. Simulations based on input related to land cover, agricultural management, livestock management, soil distribution, elevation (DEM) and point sources of pollution allowed us to assess the influence of these BMPs on overland loads and in-stream loads, and their effect on the probability of exceeding water quality guidelines.

This study suggests that riparian buffer strips and grassland conversion are the most effective of the four evaluated BMPs in terms of reduction of P, N and suspended sediment loads at the watershed outlet, while precision slurry application has a tendency to increase in-stream N loads. It must be noted, however, that even with a significant reduction of loads, none of the BMPs were able to significantly reduce the probability of exceeding water quality guidelines in the Beaurivage River. This is not a problem for NO₂ ^–-N and NO₃ ^–-N, since water quality guidelines for these nutrients were never exceeded in simulations.

Because reduction rates of riparian buffer strips are very promising and their implementation as BMP is more feasible compared to the (highly efficient) grassland conversion BMP, further research on the use of riparian buffer strips is recommended. The nonlinear empirical formula used for calculating the trapping efficiency of buffer strips as a function of their width is too much of a simplification of reality and, therefore, future work should focus on alternative approaches in which the optimal width for riparian buffer strips is based on the desired trapping efficiency and physical parameters such as slope, vegetation characteristics, sediment characteristics and runoff (Gumiere et al. 2011; Gumiere and Rousseau 2011).

Further research is also required to characterize the dynamics of stream deposition and re-suspension of sediments. In its current form, the QUAL2E water quality model used in GIBSI simulates stream deposition and subsequent re-entrainment of sediments throughout the year, whereas it is more plausible that stream deposition occurs during high flow, shortly after the initiation of soil erosion in upstream areas. Also, stream deposition and stream bed erosion are more likely to occur in the bends of a stream (erosion on the outside of a bend and deposition on the inside) rather than being distributed uniformly. Future model components may be able to simulate time-variant, distributed stream deposition-erosion.

This case study of the Beaurivage River watershed demonstrates that the evaluation of the effects of BMPs on water quality cannot be based on the analysis of overland loads alone but must also focus on the in-stream concentrations and loads as non-linear behaviors characterize the outcomes. This can only be addressed using an integrated modeling approach that accounts for a wide range of distributed variables (e.g., land cover, agricultural management, point sources of pollution, nutrient cycling and water quality of surface runoff) and provides a framework to assist policy makers and farmers in improving water quality in watersheds with similar characteristics.

Acknowledgements

This research was supported by the Watershed Evaluation of Beneficial Management Practices project (WEBs), initiated in 2004 as a partnership between Agriculture and Agri-Food Canada (AAFC) and Ducks Unlimited Canada. The authors wish to thank Shane Gabor (Ducks Unlimited Canada), Brook Harker, David Kiely, Éric van Bochove, and Terrie Scott (all Agriculture and Agri-Food Canada), Wanhong Yang and Yongbo Liu (both University of Guelph), for their collaboration.