A simulation-based interval two-stage stochastic model for agricultural nonpoint source pollution control through land retirement

Science of the Total Environment 361 (2006) 38 – 56 www.elsevier.com/locate/scitotenv

A simulation-based interval two-stage stochastic model for agricultural nonpoint source pollution control through land retirement B. Luo a,*, J.B. Li b, G.H. Huang a, H.L. Li a a b

Environmental Systems Engineering Program, Faculty of Engineering, University of Regina, Regina, Saskatchewan, Canada S4S 0A2 Environmental Engineering Program, University of Northern British Columbia, Prince George, British Columbia, Canada V2N 4Z9 Received 26 May 2005; received in revised form 7 September 2005; accepted 9 September 2005 Available online 20 October 2005

Abstract This study presents a simulation-based interval two-stage stochastic programming (SITSP) model for agricultural nonpoint source (NPS) pollution control through land retirement under uncertain conditions. The modeling framework was established by the development of an interval two-stage stochastic program, with its random parameters being provided by the statistical analysis of the simulation outcomes of a distributed water quality approach. The developed model can deal with the tradeoff between agricultural revenue and boff-siteQ water quality concern under random effluent discharge for a land retirement scheme through minimizing the expected value of long-term total economic and environmental cost. In addition, the uncertainties presented as interval numbers in the agriculture-water system can be effectively quantified with the interval programming. By subdividing the whole agricultural watershed into different zones, the most pollution-related sensitive cropland can be identified and an optimal land retirement scheme can be obtained through the modeling approach. The developed method was applied to the Swift Current Creek watershed in Canada for soil erosion control through land retirement. The Hydrological Simulation Program-FORTRAN (HSPF) was used to simulate the sediment information for this case study. Obtained results indicate that the total economic and environmental cost of the entire agriculture-water system can be limited within an interval value for the optimal land retirement schemes. Meanwhile, a best and worst land retirement scheme was obtained for the study watershed under various uncertainties. D 2005 Elsevier B.V. All rights reserved. Keywords: Land retirement; Nonpoint source pollution; Simulation; Optimization; Two-stage stochastic; Uncertainty

1. Introduction Agriculture is the basic economic component for many countries around the world. However, agricultural land is often associated with nonpoint source (NPS) * Corresponding author. Tel.: +1 306 585 5631; fax: +1 306 525 9989. E-mail addresses: [email protected] (B. Luo), [email protected] (J.B. Li), [email protected] (G.H. Huang), [email protected] (H.L. Li). 0048-9697/$ - see front matter D 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.scitotenv.2005.09.053

pollution which occurs when the surface runoff carries and deposits effluents into downstream water bodies. The agricultural effluents may include sediments from eroded or overgrazed lands, pesticides, nutrients, and other chemicals. The wash-off of such effluents through runoff process may result in serious environmental concerns to the downstream water quality. In Canada, for example, it was found that 48%, 63%, and 100% of the sampled farm dugouts were contaminated with pesticides in the provinces of Alberta, Ontario, and Saskatchewan, respectively (Coote and Gregorich,

B. Luo et al. / Science of the Total Environment 361 (2006) 38–56

2000; Renzetti, 2005). On the other hand, soil erosion resulting from farming practices has also been of increasing concern in several provinces in Canada. It not only represents the loss of a valuable asset to farmers but also may impose costs on the Canadian society by increasing water turbidity and by facilitating the transportation of other contaminants (Renzetti, 2005). As a result, it is of great importance to address the boff-siteQ environmental consequences due to agricultural NPS pollution (Contant et al., 1993; Tonderski, 1996; Horan and Ribaudo, 1999; Trevisan et al., 2000). The mitigation of NPS pollution may need to consider a number of factors, such as the type and intensity of agricultural practices, crop and land management practices, the type and amount of chemicals used, and soil characteristics (Acton and Gregorich, 1995). One effective method to address the adverse NPS pollution impacts is through the introduction of land retirement programs, which have been widely implemented around the world (Barling and Moore, 1994; Fennessy and Cronk, 1997; Wenger, 1999; Dosskey, 2001). Land retirement is the practice of taking highly erodible and/ or sensitive agricultural land out of crop production and/or grazing, and converting it by planting with a permanent vegetative cover (Ribaudo et al., 1994; Uri, 2001). This practice stabilizes the soil and reduces the movement of sediment and other chemicals from the land. Some typical examples of the land retirement program include the National Conservation Buffer Initiative established by the United States Department of Agriculture (USDA) (Dosskey, 2001; Wallender et al., 2002), and the Watershed Evaluation of Beneficial Management Practice (WEBs) and Greencover Canada program established by the Agriculture and Agri-Food Canada (AAFA, 2005). While being able to effectively mitigate the NPS pollution problems, the cost of land retirement could be huge especially when large area of croplands need to be retired. Thus, there is a necessity to systematically evaluate the tradeoff between agricultural revenue and environmental concerns for carrying out any cost-effective land retirement programs. Previously, many research efforts were made to investigate optimal land use planning with the consideration of NPS pollution impacts on the downstream water quality (Das and Haimes, 1979; Chang and Chen, 1997). For example, Hamdar (1999) presented a linear programming model to study the Mississippi’s marginal agricultural land conversion to timber use. Giupponi and Rosato (1999) developed a combined socio-economic and environmental modeling system to examine the effects of alternative agricultural land use scenarios on groundwater environmental consequences in a wa-

39

tershed in northern Italy. Nikkami et al. (2002) integrated a multi-objective linear programming model and the modified universal soil loss equation (USLE) with spatial analysis system to investigate the environmental and economical impacts of soil erosion due to improper management of land-use activities for a watershed in Iran. Furthermore, the modeling system linking optimization model and distributed water quality simulation approaches has been proposed to facilitate cost-effective NPS pollution control efforts (Newham et al., 2004). One advantage of using distributed water quality simulation approaches in this kind of modeling system is to characterize the spatial variations in NPS pollutant contributions so as to identify the most pollution-related sensitive agricultural zones. Many studies were reported to simulate NPS pollution in agricultural watershed with distributed water quality simulation (Lenat and Crawford, 1994; LeBlanc et al., 1997; Fisher et al., 2000; Fitzhugh and Mackay, 2001; Leoaˆn et al., 2001; Wade et al., 2004). In particular, the Hydrological Simulation Program-FORTRAN (HSPF) (Johanson et al., 1984; Donigian et al., 1984) has obtained extensive applications with the capability and robustness in longterm hydrological and contaminant transport simulation (Bicknell et al., 1997; Tong and Chen, 2002; Hayashi et al., 2004). Usually, the agricultural NPS pollution control is associated with various uncertainties. Unlike point source pollution, the effluent discharge from agricultural land has a random nature which is closely attribute to surface runoff (Fujiwara et al., 1988). When surface runoff is low, it may lead to relatively low amount of effluent discharge; conversely, more effluents will be discharged when surface runoff is high. Specifically, this kind of uncertainty can generate a recourse problem associated with land retirement planning. Suppose more land being retired, it will cost more agricultural revenue and consequently, it also means less effluent discharge when the surface runoff is low; however, the effluent discharge can still be high when the surface runoff will be high. Thus, it is difficult to determine the area of cropland that can be retired under such uncertainty. Very few studies were reported to address this inherent uncertainty in NPS pollution control problems (Burn and McBean, 1985; Luo et al., 2005). Two-stage stochastic programming (TSP) is an effective method to deal with the recourse problem in decision-making associated with randomness (Wets, 1966). Under the standard TSP paradigm, the decisions are partitioned into two stages. A first-stage decision is made based on uncertain future events. Subsequently, once the random events have presented themselves, a

40

B. Luo et al. / Science of the Total Environment 361 (2006) 38–56

recursive or corrective action (second-stage decision) is taken against any infeasibilities arising due to a particular realization of uncertainty. The objective of TSP is often to choose the first-stage decision variables in a way that the expected value of the random second-stage function is optimized (Birge and Louveaux, 1988). Studies were reported to apply TSP to NPS pollution control planning (Patry and Marino, 1984; Haith, 1985; Luo et al., 2005). Thus, a simulation-based TSP approach can be developed by linking the TSP with a distributed water quality simulation approach to effectively quantify the randomness and spatial variations in agricultural NPS pollution control through land retirement. In addition, most information in agriculture-water systems is not of sufficient quality to be presented as PDFs. Instead, it is more feasible to describe such uncertainties by intervals (Huang, 1996, 1998; Chang and Chen, 1997). By defining the interval inputs in a linear programming as interval numbers, Tong (1994) introduced the interval number linear program. When uncertain data in the agriculture-water system are specified as interval numbers, the interval programming can be incorporated into the simulation-based TSP approach to quantify such uncertainties. In this respect, this study aims to develop a simulation-based interval two-stage stochastic programming (SITSP) model for optimal agricultural land retirement planning. The modeling framework will be based on the interval TSP model, with its random parameters provided by the modeling outputs from the distributed water quality simulation approach with HSPF. The developed model will be demonstrated in the Swift Current Creek watershed in Canada for soil erosion control through land retirement. The tradeoff between the agricultural revenue and sediment pollution concern will be investigated in a long term point view to obtain the optimal land retirement scheme in the study watershed under various uncertainties.

ment. From the viewpoint of practicability, HSPF version 11.0 is used in this study as an important component of the proposed modeling framework. As a time series management system, HSPF can generate a continuous hydrograph of streamflow at the basin outlet and produce a time series of runoff, sediment load, and nutrient and contaminant concentrations (Donigian and Huber, 1991). HSPF consists of three application modules (PERLND, IMPLND, and RCHRES) and five utility modules (Bicknell et al., 1997). It uses the concept of hydrological response unit (HRU) to divide the watershed into the so-called homogeneous hydrological response units. In each HRU, the soil layer is vertically divided into three storages (upper-zone storage, lower-zone storage, and active groundwater storage). The water flux and evapotranspiration in each HRU are calculated respectively according to the moisture conditions in these three storages. Horizontally, three types of flow components (surface runoff, interflow, and active groundwater) contribute to the streamflow routed by a nonlinear function. The Chezy–Manning equation and an empirical expression are used to route the surface runoff (Bicknell et al., 1997). The hydrological component of HSPF includes a large number of parameters. Table 1 presents the most important and sensitive parameters involved in this component. HSPF has incorporated many NPS models for simulation of pesticides, nutrients, sediment transport, and other contaminant runoff processes on land surfaces and in the soil profiles. It has been used successfully in modeling the stream loadings of sediment, nutrients, Table 1 Important parameters for the hydrologic component of HSPF (Donigian et al., 1984) Parameter

Description

Units

AGWRC INFEXP

Basic groundwater recession rate Exponent parameter in the infiltration greater than one Ratio between the maximum and mean infiltration capacities Index zone nominal storage Interflow inflow parameter Interflow recession parameter Groundwater recession flow Lower zone ET parameter Lower zone storage Lower zone nominal storage Air temperature below which evapotranspiration will arbitrarily be reduced below the value obtained from the input time series Upper zone nominal storage

d 1

2. Methodology INFILD

2.1. The Hydrological Simulation Program-FORTRAN (HSPF) HSPF is a very comprehensive watershed model including functions for simulation of runoff and effluents from various land covers and hydraulic simulation in dams and reservoirs (Donigian et al., 1984; Bicknell et al., 1997). It is usually used as a subbasin-based distributed parameter model, which enables a user to restrict the number of runoff calculation units even in a large basin by lumping data in each subdivided seg-

INFILT INTFW IRC KVARY LZETP LZS LZSN PETMAX

UZSN

mm/h d 1 mm 1 mm mm 8C

mm

B. Luo et al. / Science of the Total Environment 361 (2006) 38–56

and pesticides from agricultural lands (Moore et al., 1988; Laroche et al., 1996). The sediment transport component of HSPF is based on the works of Donigian and Crawford (1976). The simulation of sediment removal by water is through wash-off of detached sediment in storage and scour of matrix soil. The supporting management practice factor which has been added to the soil detachment by rainfall equation is based on the bPQ factor in the Universal Soil Loss Equation (Wischmeier and Smith, 1965). Overland flow plays a major role in sediment transport on watershed surface and is treated as a turbulent flow process. The volume of deposition of sand is calculated as a function of the concentration of SS, the settling velocity, and the shear stress. To calculate the volume of scouring of sand, the power function of the average flow velocity is used in each reach. Table 2 lists the important parameters in the sediment transport component of HSPF. In general, the HSPF parameters fall into two categories including fixed parameters and process-related parameters (Al-Abed and Whiteley, 2002). The fixed parameters represent measurable characteristics that remain constant throughout a simulation period, which can be established from field data and are not involved in the calibration process. Examples of such parameters include the amount of pervious and impervious area, the areal extent of each soil type, and the hydraulic characteristics of the drainage network. The processrelated parameters represent watershed properties that only become apparent during the movement of water through the watershed, such as the water interception amount, soil water amount, infiltration, evapotranspiration rate, and the lag time of different runoff. The process-related parameters do not have directly measurable physical analogues and need to be calibrated. More Table 2 Parameters for the sediment transport component of HSPF (Donigian et al., 1984) Parameter

Description

Units

AFFIX

Fraction by which detached sediment storage decrease each day, as result of soil compaction Fraction of land surface which is shielded from erosion by rainfall Exponent in the matrix soil scour equation Exponent in the soil detachment equation Exponent in the detached sediment wash-off equation Coefficient in the matrix soil scour equation Coefficient in the soil detachment equation Coefficient in the detached sediment wash-off equation Supporting management practice factor

d 1

COVER JGER JRER JSER KGER KRER KSER SMPF

Complex Complex Complex Complex Complex Complex

41

detailed description of HSPF can be found in Bicknell et al. (1997). 2.2. Simulation-based interval two-stage stochastic programming model The simulation-based interval two-stage stochastic programming (SITSP) model is developed here by coupling the HSPF simulation results with the interval two-stage stochastic programming model. It is common that a number of economic, environmental and social factors may affect the planning of a land retirement program. Such factors include the economic costs typically represented by the loss of cropping returns, the erodibility of agricultural land, the environmental objectives such as the abatement of sediment and nutrients and preserving wildlife habitat and aquatic life, the incentive payments for voluntary retirement of land by farmers, and the rural unemployment cost due to the retirement of active farmland (Yang et al., 2003). In this study, we assume that the decision problem in land retirement program from the planners’ perspective is to select a subset of land to be retired from the total eligible land within a watershed, so as to minimize the total system cost resulting from agricultural production loss and environmental abatement of adverse downstream water quality. In this case, only the intrinsic conflict or tradeoff between the goals of maintaining agricultural production and improving the water quality in surface waters is addressed. For example, less retired land may result in less cost due to the loss of agricultural production but more cost of environmental concerns (water quality deterioration in downstream water bodies and high effluent concentration in sensitive sites). On the contrary, more retired land may lead to more cost due to agricultural production loss but less environmental cost. In this respect, the optimal land retirement scheme can be approached by minimizing the mean annual system cost for the entire agriculture-water system (without loss of generality, here we only consider the annual model). However, the random nature of effluent discharge from the watershed can complicate such efforts. When the annual effluent discharge is high, more cropland should be retired; conversely, less cropland needs to be retired when the annual effluent discharge is low. Thus, the decision makers face a problem of determining the area of cropland that needs to be retired in such an uncertain environment. This is a typical recursive problem which can be stated as a two-stage stochastic program (TSP). The first-stage decision is the area of cropland that is expected to be

42

B. Luo et al. / Science of the Total Environment 361 (2006) 38–56

retired, and the second-stage decisions are based on the random effluent discharge. By considering the spatial variations of NPS pollution contributions through subdividing the whole watershed into several agricultural zones (segments), the problem under consideration can be stated as follows: m  X  min f ¼ ai bi xi þ Ey˜ i ½bd y˜ i  þ E t˜½cd t˜ ð1aÞ i¼1

s:t: 0Vxi Vxm i VAi ;

8i

y˜ i zg˜ i ðAi  xi ; e˜ i Þ; T zt˜zh˜ s˜ ;

X

8i

ð1bÞ ð1cÞ

! xu z0

ð1dÞ

uai

where f is total system cost; x i is area of cropland that is to be retired for agricultural zone i (acre); A i is cropped area for zone i before land retirement (acre); xm i is maximum area of cropland that can be retired for zone i; a i is net benefit per tonne crop in zone i ($/ tonne); b i is crop yield per acre of land for zone i (tonne/acre); y˜ i is annual amount of effluents discharged from zone i when its cropped area becomes A i  x i (tonne); e˜ i is annual load of effluents from zone i before land retirement (tonne), which can be obtained through the simulation approach with HSPF; b is environmental costper tonne of effluents discharged into the downstream water bodies ($/tonne); ˜t is number of days per year for a sensitive site within the watershed that has daily effluent concentration above a tolerable level (s) after land retirement (day); s˜ is number of days per year for the sensitive site with daily effluent concentration above a tolerable level (s) before land retirement (day), which can be obtained based on the simulated results e˜ ; c is daily cost of the sensitive site when its daily effluent concentration exceeds the tolerable level ($/day); T is total number of days in a year (i.e., 365 days); g˜ i means that y˜ i is a function of (A i  x i ) and e˜ i , indicating that the annual amount of effluent discharged from zone i is related to that zone’s cropped area and total effluent discharge before land ˜ means that ˜t is a function of s˜ and retirement; hP P uai xu ; uai xu is total area above the sensitive site that can be retired; and E n (n) is expectation value of random variable n, and the parameters with b ˜ Q represent random variables. Constraint (1b) ensures the maximum area of cropland that can be retired in each agricultural zone; Constraint (1c) reflects that the amount of effluent from zone i should be no less than g˜ i since the cropped area A i is usually no larger than the total area of zone i;

Constraint (1d) reflects that ˜t should be no less than h˜ due to the same reason in constraint (1c). Model (1) is a TSP program with the first-stage decision x i and second-stage decisions y˜ i and ˜t . When x i is high, more land need to be retired, leading to low values in y˜ i and ˜t , which means low environmental cost. Conversely, low values in x i means high environmental cost. Therefore, a recursive action needs to be taken in order to get optimal x i . Model (1) has numerical solutions only when the probability density functions (PDFs) of the random variables y˜ i and ˜t are known. To solve this problem with linear programming, such PDFs must be piecewise linearized. According to Loucks et al. (1981),  Pn let y˜ i take values y ij with  discrete ˜ probability pij p ¼ 1; j ¼ 1; 2; . . . ; n j¼1 ij  PK ; and t take values t k  with probability pkV p V ¼ 1; k¼1 k k ¼ 1; 2; . . . ; K : Then, model (1) can be reformulated as follows: " # m n K X X X ai bi xi þ bd pij yij þ cd pkVtk min f ¼ j

i¼1

k¼l

ð2aÞ s:t: 0Vxi Vxm i VAi ;

8i

  yij zgij Ai  xi ; eij ; T ztk zhk sk ;

X

ð2bÞ 8i; j

ð2cÞ

! xu z0;

8k

ð2dÞ

uai

where y ij is the annual amount of effluent from zone i under probability p ij when that zone’s cropped area is A i  x i ; e ij is the annual amount of effluent from zone i under probability p ij before land retirement, and it can be obtained through the fitted PDF of e˜ i ; s k is number of days per year for the sensitive site with daily effluent concentration above a tolerable level before land retirement, which can be obtained through the fitted PDF of s˜; t k is the number of days per year for the sensitive site when its daily effluent concentration exceeds the tolerable level under probability p kV. Both g ij and h k are functions of the deterministic variables that will be defined later. Obviously, model (2) can reflect uncertainty presented as PDFs. However, in most agricultural NPS pollution control problems, uncertainty presented as interval numbers is more straightforward than PDFs due to the poor quality of information that can be obtained (Huang, 1996; Chang and Chen, 1997). For example, some modeling coefficients, such as a i , b i , b, and c in model (2), may not be available as determin-

B. Luo et al. / Science of the Total Environment 361 (2006) 38–56

istic values. Instead, it is more straightforward to specify them as interval numbers in most cases. Thus, by introducing the interval number programming to quantify this kind of uncertainty, model (2) can be transformed into the following simulation-based interval two-stage stochastic programming (SITSP) model: min f F ¼

m X

F F aF i bi x i þ b d

m X n X

i¼1

þ cF d

K X

pij yij

j

i¼1

ð3aÞ

pkVt k

k¼l

s:t: xm i zxi z0;

8i

ð3bÞ

sF k ;

X

Thus, we can have the following expression for g ij in constraint (2c): gij ¼

ðAi  xi Þ eij ; Ai

8i; j:

ð4Þ

Similarly, we assume that the PDF of ˜t has a similar structure as that of s˜ . This implies that a linear relationship can be established between t k and s k under a given probability p kV. Meanwhile, it is also reasonable to assume that t kP has a linear relationship with the total mV cropped area i¼1 ðAu  xu Þ above the sensitive site. We thus can have the following expression: mV P ðAu  xu Þ hk ¼ u¼1 P sk ; uai; 8k ð5Þ mV Au u¼1

  yij zgijF Ai  xi ; eF ij ;

T ztk zhF k

43

8i; j

ð3cÞ

! xu z0;

8k

ð3dÞ

uai F F F F F F F are interval where aF i , b i , c i , b , c , s k , e ij , and f coefficients/objectives that belong to the set of interval numbers (Tong, 1994). For example, letting X i and X i+ be the minimum and maximum bounds of interval number X iF, respectively, we have X i V X i+ and X iF = [X i, X i+]. Although e ijF have random characteristics, they are also presented as intervals. This is because it always has system errors when using the statistical analysis to obtain e ij from simulation outcomes e˜ i . Similarly, sF k are also expressed as intervals with random characteristics. The values of e ijF and sF k can be determined based on those of e ij and t k , respectively.

where mV is number of agricultural zones above the sensitive site (mVV m). By substituting expressions (4) and (5) into model (3), the following interval linear TSP problem can be obtained. m m X n X X F F min f F ¼ aF b x þ b d pij yij i i i i¼1

þ cF d

j

i¼1 K X

ð6aÞ

pkVtk

k¼1

s:t: xm i zxi z0; yij z

ð A i  xi Þ F eij ; Ai mV P

T ztk z

u¼1

8i

ð6bÞ 8i; j

ðAu  xu Þ sF k ; mV P Au

ð6cÞ

uai;

8k:

ð6dÞ

u¼1

2.3. Method of solution Model (3) is a simulation-based interval TSP (SITSP) model, which can only be solved after deterministic model (2) is solved. To solve model (2), constraints (2c) and (2d) should be determined first. Theoretically, both y ij and t k can be determined with a water quality simulation approach. However, the PDFs of both y˜ i and ˜t could shift when values of x i change. This makes it formidable to solve model (2). In this study, we assume that water pollution in the study system is mainly contributed by NPS and y˜ i has a similar format of PDF (e.g. C distribution) as e˜ i has. This implies it has a linear relationship between y ij and e ij . In HSPF, the NPS effluent discharge y ij from zone i is defined to have a linear relationship with cropped area A i  x i .

Model (6) is a simulated-based interval linear TSP problem with minimized objective. According to Tong (1994), the solution of an interval linear problem can be defined to get the most favorable (best) and least favorable (worst) objective-function values. By introducing the maximum and minimum value range inequalities, an interval linear program can be turned into two conventional linear problems so as to obtain an interval number solution. Thus, we can turn model (6) into the following two linear programming models which correspond its best and worst objective function values, respectively. min f  ¼

m X i¼1

  a i b i xi þ b d

m X n X i¼1

j

pij yij þ cd

K X

pkVtk

k¼1

(7a)

44

B. Luo et al. / Science of the Total Environment 361 (2006) 38–56

s:t: xm i zxi z0;

yij z

8i

ðAi  xi Þ  deij ; Ai mV P

T ztk z

u¼1

ð7bÞ

8i; j

ð A u  xu Þ s k ; mV P Au

ð7cÞ

8k

uai;

ð7dÞ

u¼1

and m X

min f þ ¼

þ þ aþ i bi x i þ b d

i¼1

þ cþ d

m X n X

K X

pij yij

j

i¼1

ð8aÞ

pk tVk

k¼1

s:t: xm i zxi z0;

yij z

8i

ðAi  xi Þ þ deij ; Ai mV P

T ztk z

u¼1

ð8bÞ

8i; j

ð A u  xu Þ sþ k ; mV P Au

ð8cÞ

uai;

8k:

ð8dÞ

u¼1

Models (7) and (8) are both conventional linear TSP problems which can be solved by various algorithms (Birge and Louveaux, 1988; Sen, 1993; Edirisinghe and Ziemba, 1994; Ruszczynski and Swietanowski, 1997; Zhao, 2001). In this study, we used the method introduced by Birge and Louveaux (1988) to solve these two conventional TSP problems. As a result, the solutions of models (7) and (8) are presented as  + {xV,i y ijV, t kV, f opt | 8i, j, k}, and {xW, i y ijW, t kW, f opt | 8i, j, k}, respectively. Thus, the solutions of model (6) can be expressed as:

)

h i F  þ fopt ¼ fopt ; fopt xiopt ¼ ðxiV; xiWÞ ; yijopt ¼ yijV; yijW tkopt ¼ ðtkV; tkWÞ

8i; j;

8k

ð9Þ

where solutions expressed in b()Q indicate two solution scenarios of the interval model (6), with the best case solution YiV(xV,i y ijV, t kV | 8i, j, k) and the worst case solution

YiW(x iW, y ijW, t kW | 8i, j, k). In this case, we don’t have the relationship YiV V YiW as interval numbers have as defined above. Fig. 1 illustrates the general framework of the SITSP method. The modeling approach is based on the distributed water quality simulation approach and two system optimization techniques, namely the HSPF simulation, the two-stage stochastic programming (TSP), and interval number programming. Each technique has a unique contribution in enhancing the model’s capability in dealing with complexities and uncertainties in agricultural land retirement planning. For example, the HSPF simulation can reflect the heterogeneous nature of effluent discharge within an agricultural watershed. The understanding of such spatial distribution is of crucial importance for an effective land retirement scheme. In addition, the recursiveness generated by random nature of effluent discharge in analyzing the tradeoff between agricultural revenue and downstream water quality concern were handled through TSP. Furthermore, the uncertainties presented as interval numbers were reflected through interval number programming. The solution algorithm of model (6) can then be summarized by using the following pseudo-code: Step 1: Calibrate and verify the hydrological component of HSPF; Step 2: Use the validated hydrological component of HSPF to estimate time series of streamflow and runoff in sub-watersheds with poor records; Step 3: Use the validated hydrological components of HSPF and surveyed water quality data to validate the water quality simulation components of HSPF in sub-watersheds with water quality records; Step 4: Use the validated HSPF to generate the time series of effluent discharge for each agricultural zone before land retirement (e˜ i ); Step 5: Use the estimated effluent discharge (e˜ i ) to determine the time series of number of days with daily SS at the sensitive site exceeding the tolerable level before land retirement (s˜); Step 6: Fit the PDFs of e˜ i and s˜ to obtain the discretization values e ij and s k , and then calcuF late their interval values eF ij and s k under given probability levels; Step 7: Formulate the SITSP model (6); Step 8: Transform model (6) into two sub-models corresponding to the best and worst objective-function values;

B. Luo et al. / Science of the Total Environment 361 (2006) 38–56

Optimal land retirement planning under uncertainty

Distributed watershed modeling approach with HSPF 11.0

Randomness in effluent discharge

Model calibration and verification

Targeted land retirement scheme ( x ) i

Two-stage stochastic programming (TSP)

Uncertainties in economic and technical data

Interval number programming

Simulation of daily sediment concentration above tolerable level at the sensitive site before

45

Simulation of annual sediment load from each agricultural zone before land retirement ~ ( ei )

land retirement

Simulation-based Interval Two-stage Stochastic Programming (SITSP) Model

SITSP most favorable sub-model (7)

SITSP least favorable sub-model (8)

Time series of number of days with daily sediment concentration at the sensitive site exceeding the tolerable level before land ~ retirement ( τ )

Fitting the probability density functions ~ ~ of e i and τ

± ± Obtain the interval values of e ij and τ k under different discretization of probability levels Solutions of submodel (7)

Solutions of submodel (8)

Optimal land retirement scheme

Fig. 1. Schematic of the SITSP methodology.

Step 9: Solve sub-model (7) and obtain its solutions  {xV,i y ijV, t kV, f opt | 8i, j, k}; Step 10: Solve sub-model (8) and obtain its solutions + {xW, i y ijW, t kW, f opt | 8i, j, k}; Step 11: Obtain the solution of the SITSP model and get the optimal land retirement scheme for each agricultural zone (x iopt); Step 12: STOP. 3. The study system 3.1. The Swift Current Creek watershed The major Canadian agricultural operations are found in the Prairie (Alberta, Saskatchewan, and Manitoba) and the central provinces (Ontario and Quebec). For example, nearly 80% of all Canadian farmland is located in the Prairie (Renzetti, 2005). With the growing recognition of the challenges posed by nonpoint source (NPS) pollution, many provinces in Canada have extended their water quality protection efforts in the direction of regulating land uses. As a result, the investigation of NPS control planning for the Prairie agricultural watersheds is of practical importance. In this study, the developed SITSP modeling approach is applied to the Swift Current Creek watershed

that is located in southern Saskatchewan, Canada (Fig. 2). As a tributary of the South Saskatchewan River, the Swift Current Creek originates from the eastern portion of the Cypress Hills and flows northeasterly through the City of Swift Current (with an area around 16 km2) and finally into Lake Diefenbaker. The entire watershed has a total area of 4934 km2, including the Rock, Jones, Bone, and Lac Pelletier Creeks. The Reid Lake, which is formed by the Duncairn Dam in the middle reach of the Swift Current Creek, is the source of fresh water for the City of Swift Current. This lake is a small reservoir which is not intended to be operated as a flood control facility. It is a semi-arid climate region in the watershed with high temperature and low precipitation in summer. The mean annual snowfall and rainfall in the watershed are around 108 cm and 379 mm, respectively. The main land cover in the watershed is cropland with cultivated crops, including oilseeds (e.g. sunflower and canola), pulses (e.g. bean and pea), cereal (e.g. soft and hard wheat), and forage (e.g. alfalfa and timothy). Along with the water bodies and urban area of the Swift Current Creek city, the watershed is covered with grassland in its upper reach. The grassland consists of spear grass, wheat grass, june grass, blue grama grass, silver sagebrush, and cactus. Effluents from the cropland, such as pesticides, sediment, and total nitrogen,

46

B. Luo et al. / Science of the Total Environment 361 (2006) 38–56

Fig. 2. The Swift Current Creek watershed.

provide major contaminant contributions to the surface waters within and downstream this watershed. 3.2. Problem statement In the Swift Current Creek watershed, soil erosion is a direct consequence of farming and sediment in surface waters is mainly from farmland. The Reid Lake is the source of fresh water of the City of Swift Current and water quality in the Lake is influenced by the sediment. The cost of treating sediment is high especially when it has high suspended solids (SS) concentration. In addition, sediment loads also generate environmental problems in the downstream water body, namely the Lake Diefenbaker (McLaughlin and Yuzyk, 1981; Yuzyk, 1983). To reduce sediment from farmland in the watershed, one effective method is to implement a land retirement program. Considering the huge costs in cutting the cropped area, the tradeoff between agricultural revenue and water quality improvement should be measured together so as to minimize the long-term cost of the entire agriculture-water

system. However, uncertainties presented as randomness and intervals in the system can complicate such efforts. Moreover, the spatial variations of sediment contribution within the watershed should be reflected for an optimal land retirement scheme so as to identify and reduce the most erodible cropped area. In this case, the developed simulation-based interval two-stage programming (SITSP) model can be applied here to deal with this land retirement problem. 3.3. Input data Extensive data are needed for the modeling approach of the SITSP model. The soil data are important for the sediment simulation with HSPF. The study watershed is dominant by three soil types, including chernozemic, solonetzic, and vertisolic. Vertically, most of the soil types consist of three layers. The first layer with a depth of 10 to 15 cm consists of dark topsoil where most of the organic (humus) materials are present. This is the layer where most of the plowing occurs. The second layer is heavier and drier in texture than that of

B. Luo et al. / Science of the Total Environment 361 (2006) 38–56

47

Table 3 Meteorological and hydrometric stations, sub-watersheds, drainage areas, and record periods Hydrometric station

Sub-watershed

Latitude

Longitude

Flow records

Drainage area (km2)

Sediment records

05HD036 05HD031 05HD041 05HD030 05HD007 05HD039 05HD037

1-036 2-031 3-041 4-030 5-007 6-039 7-037

49850V39W N 49857V30W N 5084V51W N 5089V53W N 50816V50W N 50829V38W N 50836V6W N

108828V46W W 108811V0W W 10882V40W W 107853V47W W 107846V40W W 107839V31W W 107843V22W W

1955–2001 1944–1956 1975–1992 1944–1974 1910–1974 1973–2001 1965–1972

1430 1920 2470 2910 3280 3730 3910

1973–1982 1965–1972

Climate station

Latitude

Longitude

Record period

Elevation (m)

4028040 4028060

50818V N 50816V N

107841V W 107844V W

1938–2002 1959–2002

818 825

ric stations. Table 4 lists the information of each zone. The cropped area in each zone was obtained from the Irrigation Crop Diversification Corporation (ICDC) of Saskatchewan. Currently, it is still hard to determine the maximum cropped area that can be retired (xm i ) in each zone due to the complexities in decision-making, and we thus assumed that its value equals to one third of the cropped area in each zone before land retirement. Since topologic information is important for HSPF simulation, slope of the main channel of each agricultural zone was extracted by using ArcGIS with a 100-m resolution digital elevation model (DEM). Table 4 presents the slope, total cropped area, and maximum area that can be retired for each agricultural zone. Data about crops in the study watershed were investigated through communications with local farmers and engineers in the ICDC. It is difficult to determine the cropped area of a particular crop planted in each zone since farmers usually grow different crops in the same land in different years. However, the dominant crop types planted in each zone can be identified, which is listed in Table 5. Meanwhile, the ICDC provides annual reports of statistics of net benefit per land for various crops in the study watershed. For example, Table 5 lists such values of different crop types for year 2004. Based on such statistical values, we estimated the net agricul-

the first layer, and occurs from 10 to 61 cm. The third layer consists mostly of unaltered material that is poor in nutrient value. The daily discharge and sediment data (including daily concentration and load) were obtained from Water Survey of Environment Canada. Fig. 2 presents the study watershed and the seven selected hydrometric stations in it. Table 3 lists the station number and geographic coordinate of each station, and the drainage area of the sub-watershed controlled by each station. The climatic data was obtained from the National Climate Data and Information Archive of Environment Canada. The information of the selected two climate stations is listed in Table 3. The meteorological data used for HSPF simulation include daily maximum and minimum temperature, precipitation, and evapotranspiration. In the study watershed, each of the seven hydrometric stations corresponds to a sub-watershed with its drainage area controlled by the related station. The drainage area of each sub-watershed was obtained from Water Survey of Environment Canada, which is also listed in Table 3. To identify the most erodible cropped areas and locate the cropland to be retired, the whole watershed was then divided into 7 agricultural zones (segments) according to the distribution of dominant crops in the watershed. Except for zone 1, each zone is located between two adjacent hydrometTable 4 Slope and area for each agricultural zone Agricultural zone (i)

Between hydrometric stations

Slope

Total area (km2)

A i (acre)

xmi (acre)

i =1 i =2 i =3 i =4 i =5 i =6 i =7

Above 05HD036 05HD036-05HD031 05HD031-05HD041 05HD041-05HD030 05HD030-05HD007 05HD007-05HD039 05HD039-05HD037

0.00422227 0.00395537 0.00391852 0.00212726 0.00193086 0.00477097 0.00537466

1430 490 550 440 370 450 180

141,344 96,865 108,726 54,363 60,953 66,718 33,359

47,115 32,288 36,242 18,121 20,318 22,239 11,120

48

B. Luo et al. / Science of the Total Environment 361 (2006) 38–56

Table 5 Agricultural and economic data for each agricultural zone Agricultural zone (i)

i =1 i =2 i =3 i =4 i =5

i =6 i =7

Dominant crops

Alfalfa Barley Canola Soft wheat Timothy Hard wheat Lentil Fababean Dry bean Fax Pea Corn silage Durum wheat Mustard Corn silage

F aF i ! b i ($/acre)

Net benefit per land ($/acre) Low value

Average value

37 70 1 6 48 77 12 3 83 7 46 121 31 47 121

112 187 22 26 74 149 30 83 198 423 3 63 84 156 42 87 215 347 68 106 3 39 263 405 12 76 87 127 263 405 c F = [2000, 3000] ($/day)

b F = [250, 350] ($/tonne)

tural benefit per acre land for each of these seven agricultural zones. We also expressed such values F (aF i d b i ) as interval values based on their averaged and extreme (low and high) values as listed in Table 5. The capacity of the water plant in the City of Swift Current is approximately 2  106 gal per day. The intake of the water plant is located above the Duncaim Dam and near the hydrometric station 05HD041. Thus, sediment from sub-watershed 3-041 contributes most of the sediment in this sensitive site. In this study, we set the daily suspended solids (SS) concentration with 10 mg/l as the tolerable level based on the Saskatchewan surface water quality objective (SERM, 1997). Above this level, the water treatment plant will get high cost in treating the sediment. Values of c F in model (6) are estimated based on the operation of water plant in the City of Swift Current and presented in Table 5. However, it is still difficult to measure the environmental cost with economic values for sediment accumulation in the Lake Diefenbaker. Based on previous studies on environmental problems from sedimentation in that Lake (McLaughlin and Yuzyk, 1981; Yuzyk, 1983), values of b F were estimated as shown in Table 5. 4. Modeling approach To get the optimal land retirement scheme for the F study watershed with model (6), values of sF k and e ij should be obtained first. To get these two set interval values, values of s k and e ij should be determined. This can be conducted through statistical analyses with time series data of annual sediment load (e˜ i ) from each agricultural zone and SS concentration in the sensitive

High value [90, 125] [60, 85] [170, 200] [70, 90] [145, 175]

[200, 235] [125, 150]

site (intake of the water plant). In this study, the time periods for statistical analysis were set from year 1972 to 2001 (30 years) which is also the time scale for HSPF simulation. However, except that hydrometric stations 05HD039 and 05HD037 have short period of sediment records (around 7 years), all other stations have no sediment records. The HSPF was then applied to estimate and extend sediment loads and concentrations in stations with no records and stations with short period of records, respectively. It should be noted that in this case study each agricultural zone is viewed as one hydrological response unit (HRU) in the simulation process. Thus, meteorological and land cover data in each zone were lumped in each zone for the simulation process. Both the hydrological component and sediment transport component of HSPF were simulated in a daily time step. To simulate with HSPF, the model was validated first in sub-watersheds with surveyed hydrological and sediment records. The purpose of model validation is to determine the process-related parameters of HSPF in each sub-watershed. The Nash–Sutcliffe coefficient R 2 (Nash and Sutcliffe, 1970) was used as a criteria to measure the simulation performance. Generally, the validation process can be thought as satisfied if the R 2 values from both calibration and verification are close or above 0.8. The streamflow simulation was conducted first since the outcomes of hydrological simulation provide bases for the sediment simulation. During the validation process of hydrological simulation, the study period was divided into a calibration period and a verification period as shown in Table 6. Since stations 05HD036

B. Luo et al. / Science of the Total Environment 361 (2006) 38–56

49

Table 6 Period of calibration, verification, and extended time series in simulation process Sub-watershed

1 2 3 4 5 6 7

Hydrological simulation

Sediment simulation

Calibration period

Verification period

1972–1991 1944–1950 1975–1985 1945–1964 1945–1964 1973–1991 1965–1969

1992–2001 1951–1956 1986–1992 1965–1974 1965–1974 1992–2001 1969–1972

and 05HD039 have long period of daily discharge records (30 years), the hydrological simulation was conducted first in the related sub-watersheds 1-036 and 7-039 to determine the process-based parameters. Then, the calibrated parameters from these two simulation approaches were used as references to calibrate the model in other 5 sub-watersheds. The calibration period represents a combination of dry, average, and wet years. A manual calibration process was conducted to adjust the modeling parameters. The initial parameter values in the calibration process were approximated based on the HSPF User’s Manual (Bicknell et al., 1997). Each of the sensitive parameters identified in previous studies (Donigian et al., 1984; Linsley et al., 1986; Al-Abed and Whiteley, 2002) was changed individually and the corresponding simulated daily discharges were compared with observed data so that the values of R 2 can reach above or close the desired value (0.8). The sensitive parameters for the hydrologic component of HSPF include LZS, LZSN, PETMAX, and INFILT as shown in Table 1. After calibration, the model was verified using the observed data of the verification periods without any change of the calibrated modeling parameters. We then used the validated hydrological component of HSPF to extend the time series of daily streamflow from 1972 to 2001 for sub-watersheds 2031, 3-041, 4-030, 5-007, and 7-037. The simulated streamflow and surface runoff in each sub-watershed were further used in the relevant sediment simulation process. Table 3 shows that only hydrometric stations 05HD039 and 05HD037 have sediment records. The sediment transport component of HSPF was then first validated in sub-watersheds 6-039 and 7-037. Similar to the validation of the hydrologic component of HSPF, the study period of each of these two stations was also divided into a calibration and a verification period (Table 6). The initial parameter values in the calibration process were also approximated on the basis of the HSPF User’s Manual (Bicknell et al., 1997). According to studies by Fontaine and Jacomino (1997), the simulation

Extension period

Calibration period

Verification period

1973–1978 1965–1968

1979–1982 1969–1972

1972–2001 1993–2001 1975–2001 1975–2001 1973–2001

of sediment transport is highly sensitive to the parameters JRER, JSER, and JGER as listed in Table 2, so they were adjusted carefully. The validated model was further applied to extend daily sediment load in these two sub-watersheds from 1972 to 2001, with continuous inputs of meteorological data. The validated model was further applied to estimate time series of daily SS and load (from 1972 to 2001) in the five sub-watersheds without sediment records and SS concentration in the outlet of sub-watershed 3-041 (sensitive site), with no changes in process-related parameters and continuous inputs of meteorological data. As a result, the daily SS concentration and sediment load were obtained in each sub-watershed and the sensitive site. The annual sediment load in each sub-watershed was further calculated based on the relevant simulated daily values. Then, the annual sediment load from each agricultural zone was calculated by the subtraction of annual sediment loads in two adjacent sub-watersheds which encompass the corresponding agricultural zone. For example, the annual sediment load from agricultural zone 2 was obtained through subtracting the annual sediment load in sub-watershed 1-036 to the relevant sediment load in sub-watershed 2-031. For agricultural zone 1, its annual sediment load equals to the relevant value in sub-watershed 1 since it is located in the upper stream of the study watershed. Thus, we obtained the time series of annual sediment load (e˜ i ) for each agricultural zone and the time series of number of days per year (s˜) that daily SS concentration in the sensitive site exceeds the tolerable level. Three candidate distributions, namely the Gamma, the 2-parameter Lognormal and the Log Pearson type III distributions, were used to fit the PDFs of random variables e˜ i and s˜, respectively. Two commonly used parameter estimation techniques for the 2-parameter Lognormal and the Log Pearson type III distributions were applied in this study, including the method of moments and the maximum likelihood estimation (MLE) (Kite, 1988). For the Gamma distribution, the

50

B. Luo et al. / Science of the Total Environment 361 (2006) 38–56

Table 7 Results from streamflow simulation Sub-watershed

1 2 3 4 5 6 7

Hydrometric stations

05HD036 05HD031 05HD041 05HD030 05HD007 05HD039 05HD037

N–S coefficient R 2

Mean annual runoff volume (106 m3)

Calibration

Verification

Observed

Simulated

0.92 0.83 0.85 0.88 0.86 0.93 0.79

0.87 0.80 0.82 0.83 0.83 0.89 0.76

27.054 29.776 31.865 33.865 35.962 38.398 39.824

28.036 29.165 31.683 34.275 36.458 38.862 40.363

maximum likelihood estimation of parameters was utilized (Macuen and Snyder, 1986). After the PDFs of these two random variables were determined, their discretization values (i.e., e ij and s k ) with probability 0.1 were calculated. It was further assumed that values of s k and e ij each have a 15% deviation to their fitted values due to system uncertainties, and then obtained the values of interval number e ijF and s kF. By inputting the values of these two simulated interval numbers and economic data introduced in Section 3.3, the optimal land retirement scheme was obtained with model (6) for the study watershed. 5. Results 5.1. Simulation results with HSPF Since the purpose of this case study is to get the optimal land retirement scheme, here we only present parts of the outcomes from the simulation process with HSPF. Table 7 lists the Nash–Sutcliffe coefficient R 2 values for the hydrologic simulation in each sub-watershed, and it also presents the corresponding observed and simulated mean annual runoff volume. It shows that the longer the hydrological records, the better the validation results. Fig. 3 illustrates the observed and

simulated daily discharge in the verification process for sub-watershed 1-036 in 2001, and Fig. 4 presents the simulated daily discharge for sub-watershed 7-037 in the same year. Table 8 lists the R 2 values for the sediment simulation in the two sub-watersheds with sediment records. It also presents the observed and simulated mean annual sediment load from each sub-watershed. According to the simulated mean annual sediment load from the whole study watershed (i.e., sub-watershed 7-037) and the total cropped area, the mean annual sediment yield per acre cropped land of the study watershed is about 0.06 tonne. Fig. 5 illustrates the observed and simulated daily SS concentration in the verification process for sub-watershed 7-037 in 1972. It shows that there is no sediment load in winter since most precipitation is snowfall during that season. Fig. 6 presents the simulated annual runoff volume and annual sediment load from sub-watershed 7-037. Since the drainage area of sub-watershed 7-037 equals to the drainage area of the whole study watershed, annual sediment load from this sub-watershed represents the relevant annual load into the Lake Diefenbaker from the Swift Current Creek watershed. Based on the simulation results, the annual sediment load from each agricultural zone was calculated. Fig. 7 7 Simulated daily flow

Observed Simulated

2.5 2 1.5 1 0.5

6

Daily discharge (m3/s)

daily discharge (m3/s)

3

5 4 3 2 1

0 1

26 51

76 101 126 151 176 201 226 251 276 301 326 351

Time (year 2001)

0 1

26

51

76 101 126 151 176 201 226 251 276 301 326 351

Time (year 2001)

Fig. 3. Observed and simulated daily discharge for sub-watershed 1-036 in 2001.

Fig. 4. Simulated daily discharge for sub-watershed 7-037 in 2001.

B. Luo et al. / Science of the Total Environment 361 (2006) 38–56

51

Table 8 Results from sediment transport simulation Sub-watershed

N–S coefficient R 2

Hydrometric stations

Calibration 1 2 3 4 5 6 7

05HD036 05HD031 05HD041 05HD030 05HD007 05HD039 05HD037

Mean annual sediment load (tonne) Verification

0.86 0.89

0.83 0.86

Observed

Simulated

24,045 31,712

2092 7877 11,718 15,180 20,034 26,941 33,867

5.2. Fitted probability density functions (PDFs)

5.3. Optimal land retirement scheme

Through fitting with the three candidate distributions, it was found that the 2-parameter Lognormal distribution is the best fit for e˜ i and the Gamma distribution˜ is the best fit for variable s˜. The fitted PDFs of e˜ i and s are expressed as: " # 1 ðln˜e i  lzi Þ2 pffiffiffiffiffiffi exp  pðe˜ i Þ ¼ ð10aÞ 2r2zi e˜ i rzi 2p

Table 10 presents the optimized area of cropland that needs to be retired in each agricultural zone. It also lists the associated mean annual system cost under such optimal land retirement scheme. The results indicate two scenarios of optimal land retirement scheme, the most favorable (best optimum) and least favorable (worst optimum) scheme under uncertain environment. Under the most favorable scheme, the area of cropland that can be retired for zone 1 to 7 is 0, 32,288, 0, 0, 0, 22,239, and 11,120 acre respectively, with total system cost $41.22  106. Under the least favorable scheme, such value for zone 1 to 7 is 0, 32,288, 0, 18,121, 0, 22,239, and 11,120 acre respectively, with mean annual system cost $66.25  106. Such solutions have impor-

s˜ a1 e˜s =b ba CðaÞ

ð10bÞ

Daily SS concentration (mg/l)

where l zi and r zi are the mean and standard deviation of the natural logarithms of e˜ i ; C(a) is the gamma function; and a, b are fitted parameters. 450

Observed Simulated

400 350 300 250 200 150 100 50 0 1

26

51 76 101 126 151 176 201 226 251 276 301 326 351

Time (year 1972)

Fig. 5. Observed and simulated daily SS concentration for sub-watershed 7-037 in 1972.

160000 140000

180000 Annual runoff volume

160000

annual sediment load

140000

120000

120000

100000

100000

80000

80000

60000

60000

40000

40000

20000

20000

0

0

19 6 19 5 68 19 71 19 74 19 77 19 80 19 8 19 3 86 19 89 19 92 19 95 19 98 20 01

pðs˜ Þ ¼

Annual sediment load (tonne)

Table 9 lists the values of fitted parameters l z and r z in formula (10a) for each zone. Fig. 10 shows the cumulative distribution function of random variable s˜ and the fitted parameters a and b in formula (10b). The discretization values of s k and e ij under probability value 0.1 for each discretization interval were then calculated with these two fitted distributions. As a result, we have s kF = [0.85s k , 1.15s k ] and e ijF = [0.85e ij , 1.15e ij ].

Annual runoff volume (103 m3)

presents the annual sediment load from each agricultural zone. It shows that year 1995 and 1998 have high annual loads in every zone, this is due to the heavy floods in these two years. Fig. 8 illustrates the simulated daily discharge and SS concentration in the sensitive site (intake of the water plant) in year 2001. Fig. 9 presents the number of days in every individual year when the daily SS concentration in this sensitive site exceeds the tolerable level.

Year

Fig. 6. Simulated annual runoff volume and annual sediment load from the entire watershed.

B. Luo et al. / Science of the Total Environment 361 (2006) 38–56

Number of days per year (day)

50 80000 70000 60000 50000 40000 30000

Zone 1 Zone 2 Zone 3 Zone 4 Zone 5 Zone 6 Zone 7

20000 10000

45 40 35 30 25 20 15 10 5 0 19 75 19 77 19 79 19 81 19 83 19 85 19 87 19 89 19 91 19 93 19 95 19 97 19 99 20 01

0 19 75 19 77 19 79 19 81 19 83 19 85 19 87 19 89 19 91 19 93 19 95 19 97 19 99 20 01

Annual sediment load in each zone (tonne)

52

Year

Year

Fig. 7. Simulated mean annual sediment load from each agricultural zone.

tant policy meaning for decision-makers. When the most favorable land retirement scenario is implemented, less cropped area would be cut and the total system cost will be low. This policy is more favorable to farmers than to environment groups and local stakeholders (e.g. water plant in City of Swift Current Creek) that undertake cost due to sedimentation and high SS concentration. When the least favorable land retirement scenario is implemented, more cropland would be retired and the total system cost will be high. Such a policy is less favorable to farmers while welcomed by environmental groups and related local stakeholders. Whatever the land retirement schemes will be implemented, the mean annual system cost will be within its best and worst optimum. To ensure the mean annual system cost within the optimized interval, all croplands that can be retired in zone 2 need to be retired since its net agricultural benefit per land is the lowest among all zones, meanwhile, the mean annual sediment load per acre land is also high in

4.5

3.5

160 Daily discharge

140

Daily SS concentration 120

3

100

2.5 80 2 60

1.5

40

1

20

0.5 0

Daily SS concentration (mg/l)

Daily discharge (m3/s)

4

0 1 26 51 76 101126 151 176 201 226 251 276 301 326 351

Year (2001)

Fig. 8. Simulated daily discharge and SS concentration in the Reid Lake for year 2001.

Fig. 9. Number of days per year with high daily SS concentration in the Reid Lake.

this zone. For example, the area that needs to be retired in zone 2 is 32,288 acre, since the multiplication value of a 2F d b 2F has lowest value [60, 85] $/acre in all zones. However, the mean annual sediment load per acre land is also as high as [0.059, 0.80] (tonne/acre) as illustrated in Fig. 11. In Table 10, it is also shown that a large area of cropland needs to be retired in zone 6 even though the value of a 6F d b 6F is high. This is because the mean annual sediment load per acre land is [0.095, 0.128] tonne/acre in that zone, the second highest among all zones as presented in Fig. 11. Fig. 11 also shows that zone 7 has the highest erodible value with mean annual sediment load per land within [0.190, 0.256] tonne/acre. Meanwhile, the multiplication value of a 7F d b 7F is only [125, 150] $/acre in that zone. As a result, all cropped land that can be retired in zone 7 will be retired under optimal retirement scheme. Under the optimal land retirement scheme, the optimized values of y ijopt and t kopt are also obtained under given probabilities. For example, when p ij takes value 0.1 for low level ( j = 1) sediment load in zone 6 (i = 6), y 61opt has best value of 44 tonne and worst value of 59 tonne, indicating the extreme of the total annual sediment load from sub-watershed 6-039 under this probability after land retirement; when p kV takes low level (k = 1) probability value 0.1 for high level SS concentration in the Reid Lake (sensitive site), t 1opt has the best and worst values of 38 and 52 days respectively, indicating the extreme of the number of days per year with

Table 9 Parameters for the 2-parameter lognormal distribution of e˜ i Parameters

i =1

i =2

i =3

i =4

i =5

i =6

i =7

lz rz

5.932 1.815

6.409 2.353

7.021 1.612

6.814 1.850

7.297 1.721

7.626 1.652

7.356 1.71

B. Luo et al. / Science of the Total Environment 361 (2006) 38–56 0.3

Mean annual sediment load per land (tonne/acre)

1 0.9 0.8 0.7 0.6

CDF

53

a = 0.973; b = 14.12 0.5 0.4 0.3 0.2 0.1

low value 0.25

high value

0.2 0.15 0.1 0.05 0

61

56

51

46

41

36

31

26

21

16

6

11

1

0

day

Fig. 10. Cumulative distribution function of random variable s˜.

high SS concentration in the sensitive site under this probability. 6. Discussion In most previous studies, the cropland area within an agricultural watershed that needs to be retired is often determined based on specific return year of the effluent load (Hamdar, 1999), which cannot characterize the long term effects of land retirement on agricultural revenue and environmental consequences. The developed SITSP model can deal with this problem by minimizing the expected value of both economic and environment cost, and then lead to an optimal land retirement scheme by solving the recursiveness generated by land use planning and randomness in effluent discharge. The introduction of interval programming in the modeling framework can effectively treat vague information presented as interval numbers. The interval solution provides two scenarios of land retirement scheme, the most favorable and the least favorable scheme. Model (6) can also be solved through a conventional TSP method by letting all interval inputs be equal to their mid-values for this specific case. The deterministic solution of the Table 10 Optimal land retirement scheme and associated system cost Agricultural zone (i) i =1 i =2 i =3 i =4 i =5 i =6 i =7 F f iopt ($106)

x i opt (acre) b F, c F

Solution with mid-value inputs

0 32,288 0 (0, 18,121) 0 22,239 11,120 [41.22, 66.25]

0 32,288 0 0 0 22,239 11,120 52.63

Zone 1

Zone 2

Zone 3

Zone 4

Zone 5

Zone 6

Zone 7

Agricultural zone

Fig. 11. Mean annual sediment load per acre land in each zone.

conventional TSP problem then lies within the interval solution as presented in Table 10. For example, the total system cost will be a deterministic value of $52.63  106 for the conventional TSP problem which is within the interval value obtained by the SITSP model. However, the optimal land retirement scheme keeps the same as the most favorable scheme obtained by the SITSP model. This is due to the fact that the constraints and coefficients for the conventional optimization model with mid-value inputs and model (7) are different for this specific case study. The policy meaning of such solutions is that the decision-makers need only adjust the cropland area to be retired in zone 4 to get different land retirement schemes based on different balances among various interest groups. The application of the developed SITSP model to the study watershed is heavily based on the simulation parts of HSPF. The simulated outputs from HSPF can be used to identify sediment contribution from each agricultural zone, which can be further used to measure the tradeoff between agricultural revenue and water quality concern so as to locate the cropland that needs to be retired in the whole watershed. According to previous studies (Tong and Chen, 2002; Hayashi et al., 2004), HSPF has advantages of modeling water quantity and quality information in ungauged catchments. Such previous efforts strengthen the performance of simulating sediment information in sub-watersheds without water quality records. In addition, the major human activity in the study watershed is farming with minor changes in cropped area, and the watershed’s landscape has no major changes in the past decades. Thus, the simulation efforts to get the sediment information in ungauged sub-watersheds are reasonable and the simulated results are acceptable. In the developed model, we assume that there is a linear relation between y ij and e ij , and between t k and s k as well, which are expressed in formulas (4) and

54

B. Luo et al. / Science of the Total Environment 361 (2006) 38–56

(5). Since the topography is flat and the land cover is relatively homogeneous in each zone of the study watershed, this assumption is acceptable. However, for watershed which has erratic geomorphology and heterogeneous land cover, such assumptions could become unrealistic and nonlinear relationships should be used. In this case, the contribution of effluent discharge from each parcel of cropland should be simulated to identify such a nonlinear relationship. As a result, a grid-based fully distributed watershed model like SHE model (Abbott et al., 1986) is required to simulate water quality information in each subdivided grid. At the same time, the probability density distributions of y˜ i and ˜t should be automatically fitted in each step of calculating the first-stage decision variable x i . Thus, an integrated computer system with simulation and optimization should be developed to solve these kind of problems. The present study integrates different techniques to investigate the impacts of land retirement on surface water quality. The models and the methodology can provide a simple, but effective management tool for policy makers. In addition, the information derived from this study can have direct application values to decision makers for agricultural NPS pollution control through land retirement. The research findings can also contribute to the existing knowledge of the plausible surface water quality implications of land use. Though the developed approach is only for agricultural land retirement, it is also of value for other land use planning problems such as urban and forestry land use planning under uncertainty. 7. Summary and conclusions A simulation-based interval two-stage stochastic programming (SITSP) model was developed for agricultural NPS pollution control through land retirement. The model is based on the interval two-stage stochastic programming (TSP) and a water quality simulation model. The optimization and simulation models are connected with the random parameters in the optimization model provided by the simulation outcomes, and through the assumed linear relationships in the system. The developed model represents a new effort to obtain optimal land retirement scheme for agricultural NPS pollution control through minimizing the long term total economic and environmental cost. Meanwhile, it deals with uncertainties represented as probability density distributions and interval numbers in the entire agriculture-water system.

The developed SITSP model was demonstrated through its application in the Swift Current Creek Watershed in Canada for soil erosion control with land retirement, with its simulation part provided by the HSPF. Modeling results can provide two extreme scenarios of optimal land retirement scheme, the most favorable and least favorable scenario, according to the optimized mean annual system cost. It can provide as a base for the decision-makers to determine an optimal land retirement scheme based on the interests of different stakeholders, at the same time, to understand the total economic and environmental cost for such a policy. Acknowledgements This work was partly supported by the Climate Change Action Fund of Canada. The authors are very grateful to Len Shumaker in Saskatchewan Watershed Authority and Les Bohron in Irrigation Crop Diversification Corporation (ICDC) of Saskatchewan for help in data collection. References AAFA (Agriculture and Agri-Food Canada). Watershed evaluation of BMPs (WEBs). http://www.agr.gc.ca/env/greencover-verdir/ webs_abstract_e.phtml, accessed on June 8, 2005. Abbott MB, Bathurst JC, Cunge JA, O’Connell PE, Rasmussen J. An introduction to the European hydrological system, systeme hydrologique Europeen, SHE, 1 history and philosophy of a physicallybased, distributed modeling system. J Hydrol 1986;87:45 – 59. Acton DF, Gregorich LJ, editors. The health of our soils: toward sustainable agriculture in Canada. Canada7 Centre for Land and Biological Resources Research, Research Branch, Agriculture and Agri-Food; 1995. Al-Abed NA, Whiteley HR. Calibration of the hydrological simulation program Fortran (HSPF) model using automatic calibration and geographical information systems. Hydrol Process 2002; 16:3169 – 88. Barling RD, Moore ID. Role of buffer strips in management of waterway pollution: a review. Environ Manage 1994;18:543 – 58. Bicknell, B.R., Imhoff, J.C., Kittle, J.L., Donigian, A.S., Johanson, R.C., 1997. Hydrological Simulation Program-Fortran, user’s manual for version 11: Rep. No. EPA/600/R-97/080, EPA, U.S.A. Birge JR, Louveaux FV. A multicut algorithm for two-stage stochastic linear programs. Eur J Oper Res 1988;34:384 – 92. Burn DH, McBean EA. Optimization modeling of water quality in an uncertain environment. Water Resour Res 1985;21:934 – 40. Chang NB, Chen HW. The application of fuzzy interval genetic algorithm for water pollution control in the river basin. Proceedings of ICOQM ’97: 1st annual conference on information system cluster of ICOQM, Jaipur, India, 5–8 Jan; 1997. Contant CK, Duffy MD, Holub MA. Determining tradeoffs between water quality and profitability in agricultural production: implications for nonpoint source pollution policy. Water Sci Technol 1993;28:27 – 34.

B. Luo et al. / Science of the Total Environment 361 (2006) 38–56 Coote DR, Gregorich LJ, editors. The health of our water: toward sustainable agriculture in Canada. Canada, Ottawa, Ontario7 Research Planning and Coordination Directorate, Research Branch, Agriculture and Agri-Food; 2000. Das P, Haimes YY. Multiobjective optimization in water quality and land management. Water Resour Res 1979;15:1313 – 22. Donigian AS, Crawford Jr NH. Modeling nonpoint pollution from the land surfaceRep. No. EPA/600/3/76/083. Athens7 U.S. EPA, Environmental Research Laboratory; 1976. Donigian AS, Huber WC. Modeling of nonpoint source water quality in urban and nonurban areas. Athens, GA7 Environmental Research Laboratory, US Environmental Protection Agency; 1991. Donigian AS, Imhoff JC, Bicknell BR, Kittle JL. Application guide for the hydrological simulation program FORTRAN. Athens, GA7 Environmental Research Laboratory, US Environmental Protection Agency; 1984. Dosskey MG. Toward quantifying water pollution abatement in response to installing buffers on crop land. Environ Manage 2001;28:577 – 98. Edirisinghe NCP, Ziemba WT. Bounds for two-stage stochastic programs with fixed resources. Math Oper Res 1994;19: 292 – 313. Fennessy MS, Cronk JK. The effectiveness and restoration potential of riparian ecotones for the management of nonpoint source pollution, particularly nitrate. Crit Rev Environ Sci Technol 1997;27:285 – 317. Fisher DS, Steiner JL, Wilkinson SR. The relationship of land use practices to surface water quality in the Upper Oconee Watershed of Georgia. For Ecol Manag 2000;128:39 – 48. Fitzhugh TW, Mackay DS. Impact of subwatershed partitioning on modeled source- and transport-limited sediment yields in an agricultural nonpoint source pollution model. J Soil Water Conserv 2001;56:137 – 43. Fontaine TA, Jacomino MF. Sensitivity analysis of simulated contaminated sediment transport. J Am Water Resour Assoc 1997; 33:313 – 25. Fujiwara O, Puangmaha W, Hanaki K. River basin quality management in stochastic environment. J Environ Eng—ASCE 1988;114:864 – 80. Giupponi C, Rosato P. Agricultural land use changes and water quality: a case study in the Watershed of the Lagoon of Venice. Water Sci Technol 1999;39:135 – 48. Haith DA. Event-based procedure for estimating monthly sediment yields. Trans ASAE 1985;28:1916 – 20. Hamdar B. An efficiency approach to managing Mississippi’s marginal land based on the conservation reserve program (CRP). Resour Conserv Recycl 1999;26:15 – 24. Hayashi S, Shogo M, Masataka W, Xu BH. HSPF simulation of runoff and sediment loads in the Upper Changjiang River Basin, China. J Environ Eng—ASCE 2004;130:801 – 15. Horan RD, Ribaudo MO. Policy objectives and economic incentives for controlling agricultural sources of nonpoint pollution. J Am Water Resour Assoc 1999;35:1023 – 35. Huang GH. IPWM: an interval parameter water quality management model. Eng Optim 1996;26:79 – 103. Huang GH. A hybrid inexact-stochastic water management model. Eur J Oper Res 1998;107:137 – 58. Johanson RC, Imhoff JC, Kittle JL, Donigian AS. Hydrological Simulation Program-FORTRAN (HSPF): users manual for release 80. Athens, GA7 Environmental Research Laboratory, US Environmental Protection Agency; 1984.

55

Kite GW. Frequency and risk analyses in hydrology. Littleton, Colo., U.S.A7 Water Resources Publications; 1988. Laroche AM, Gallichand J, Lagace R, Pesant A. Simulating atrazine transport with HSPF in an agricultural watershed. J Environ Eng 1996;122(7):622 – 30. LeBlanc RT, Brown RD, FitzGibbon JE. Modeling the effects of land use change on the water temperature in unregulated urban streams. J Environ Manag 1997;49:445 – 69. Lenat DR, Crawford JK. Effects of land use on water quality and aquatic biota of three North Carolina Piedmont streams. Hydrobiologia 1994;294:185 – 200. Leoaˆn LF, Soulis ED, Kouwen N, Farquhar GJ. Nonpoint source pollution: a distributed water quality modeling approach. Water Res 2001;35:997 – 1007. Linsley RK, Kohler MA, Paulhus JLH, 1986. Hydrology for engineers. New York7 McGraw-Hill; 1986. p. 339 – 56. Loucks DP, Stedinger JR, Haith DA. Water resource systems planning and analysis. Englewood Cliffs, N.J.7 Prentice-Hall; 1981. Luo B, Maqsood I, Huang GH, Yin YY, Han DJ. An inexact fuzzy two-stage stochastic model for quantifying the efficiency of nonpoint source effluent trading under uncertainty. Sci Total Environ 2005;347:21 – 34. Macuen RH, Snyder WM. Hydrologic modeling: statistical methods and applications. Prentice-Hall. McLaughlin F, Yuzyk TR. Lake Diefenbaker, Saskatchewan: a case study of reservoir sedimentation. Can Water Resour J 1981; 7:391 – 2. Moore LW, Matheny H, Tyree T, Sabatini D, Klaine SJ. Agricultural runoff modeling in a small west Tennessee watershed. J Water Pollut Control Fed 1988;60(2):242 – 9. Nash JE, Sutcliffe JV. Riverflow forecasting through conceptual model. J Hydrol 1970;10:282 – 90. Newham LTH, Letchera RA, Jakemana AJ, Kobayashic T. A framework for integrated hydrologic, sediment and nutrient export modelling for catchment-scale management. Environ Model Softw 2004;19:1029 – 38. Nikkami D, Elektorowicz M, Mehuys GR. Optimizing the management of soil erosion. Water Qual Res J Can 2002;37:577 – 86. Patry GG, Marino MA. Two-stage urban runoff forecast model. J Water Resour Pl—ASCE 1984;110:479 – 96. Renzetti S. Canadian agricultural water use and management. Department of Economics Working Paper Series. Canada7 Brock University; 2005. http://139.57.161.145/papers/Ag_Water_Chapter.pdf, accessed on August 10, 2005. Ribaudo MO, Osborn CT, Konyar K. Land retirement as a tool for reducing agricultural nonpoint source pollution. Land Econ 1994;70:77 – 87. Ruszczynski A, Swietanowski A. Accelerating the regularized decomposition method for two-stage stochastic linear problems. Eur J Oper Res 1997;101:328 – 42. Sen S. Subgradient decomposition and differentiability of the recourse function of a two stage stochastic linear program. Oper Res Lett 1993;13:143 – 8. SERM (Saskatchewan Environment and Resource Management). Surface water quality objectives. MB. Regina, Saskatchewan, Canada7 Saskatchewan Environment and Resource Management; 1997. Tonderski A. Landuse-based nonpoint source pollution: a threat to water resources in developing countries. Water Sci Technol 1996;33:53 – 61. Tong SC. Interval number and fuzzy number linear programming. Fuzzy Sets Syst 1994;66:301 – 6.

56

B. Luo et al. / Science of the Total Environment 361 (2006) 38–56

Tong STY, Chen WL. Modeling the relationship between land use and surface water quality. J Environ Manag 2002;66:377 – 93. Trevisan M, Padovani L, Capri E. Nonpoint-source agricultural hazard index: a case study of the province of Cremona, Italy. Environ Manage 2000;26:577 – 84. Uri ND. A note on soil erosion and its environmental consequences in the United States. Water Air Soil Pollut 2001;129:181 – 97. Wade AJ, Whitehead PG, Jarvie HP, Neal C, Prior H, Johnes PJ. Nutrient monitoring, simulation and management within a major lowland UK river system: the Kennet. Math Comput Simul 2004;64:307 – 17. Wallender W, Rhoades J, Weinberg M, Lee S, Uptain C, Purkey D. Irrigated land retirement. Irrig Drain Syst 2002;16:311 – 26. Wenger S. A review of the scientific literature on riparian buffer width, extent and vegetation. Athens, GA7 Univ. of Georgia, Institute of Ecology, Office of Public Service and Outreach; 1999.

Wets RJB. Programming under uncertainty: the solution set. SIAM J Appl Math 1966;14:1143 – 51. Wischmeier WH, Smith DD. Predicting rainfall erosion losses from cropland east of the Rocky Mountains. US Department of Agriculture, Agricultural Handbook. 47 pp. Yang W, Khanna M, Farnsworth R, Onal H. Integrating economic, environmental and GIS modeling to determine cost effective land retirement in multiple watersheds. Ecol Econ 2003;46:249 – 67. Yuzyk, T.R., 1983. Lake Diefenbaker, Saskatchewan: a case study of reservoir sedimentation. Internal Report, Sediment Survey Section, Resources Branch, Inland Waters Directorate, Environment Canada. Zhao G. A log-barrier method with Benders decomposition for solving two-stage stochastic linear programs. Math Program 2001; 90:507 – 36.