simulation model, drainmod-n ii - PubAg - USDA

7 downloads 0 Views 92KB Size Report
research over many years (Gilliam et al., 1999). ...... Brevé, M. A., R. W. Skaggs, J. E. Parsons, and J. W. Gilliam. 1997. ... Chichester, U.K.: John Wiley and Sons.
SENSITIVITY ANALYSES OF THE NITROGEN SIMULATION MODEL, DRAINMOD-N II X. Wang, M. A. Youssef, R. W. Skaggs, J. D. Atwood, J. R. Frankenberger

ABSTRACT. A two-step global sensitivity analysis was conducted for the nitrogen simulation model DRAINMOD-N II to assess the sensitivity of model predictions of NO3 -N losses with drainage water to various model inputs. Factors screening using the LH-OAT (Latin hypercube sampling - one at a time) sensitivity analysis method was performed as a first step considering 48 model parameters; then a variance-based sensitivity analysis was conducted for 20 model parameters, which were the parameters ranked 1 to 14 by the LH-OAT method, five organic carbon (OC) decomposition parameters, and the empirical shape factor of the temperature response function for the nitrification process. DRAINMOD-N II simulated a continuous corn production on a subsurface drained sandy loam soil using a 40-year (1951-1990) eastern North Carolina climatological record. Results from the first 20-year period of the simulations were used to initialize the soil organic matter pools, and results from the last 20-year period of the simulations were considered for the sensitivity analyses. Both yearly and 20-year average model predictions of NO3 -N losses through drainage flow were used in the analyses. Both sensitivity analysis methods indicated that DRAINMOD-N II is most sensitive to denitrification parameters, especially those controlling temperature effect on process rate. Results also indicated that the model is mildly sensitive to the parameters controlling OC decomposition and associated N mineralization/immobilization. The use of different sensitivity analysis methods with dissimilar theoretical foundations increases the confidence in key parameters identification. More efforts should be focused on quantifying key parameters for more accurate model predictions. Keywords. Extended FAST, Drainage, DRAINMOD-N II, Nitrate-N, Sensitivity analysis.

D

eveloping agricultural management practices that maintain high crop productivity without deteriorating environmental quality has been the focus of research over many years (Gilliam et al., 1999). This complex task requires understanding of nitrogen (N) dynamics in the agro-ecosystem, which is regulated by a large number of interacting and sometimes highly dynamic physical, chemical, and biological processes. Computer models have been developed and used to simulate N dynamics in the soil-water-plant systems at different spatial and temporal scales (Singh and Kanwar, 1995; Cavero et al., 1999). DRAINMOD-N II (Youssef, 2003; Youssef et al., 2005a) is a field-scale, process-based model that was recently developed to simulate carbon (C) and N dynamics and turnover in drained agricultural lands. The model has been

Article was submitted for review in June 2005; approved for publication by the Soil & Water Division of ASABE in October 2005. The authors are Xiuying Wang, Postdoctoral Research Associate, Texas Agricultural Experiment Station, Blackland Research and Extension Center, Temple, Texas; Mohamed A. Youssef, ASABE Member Engineer, Research Assistant Professor, and Richard W. Skaggs, ASABE Fellow Engineer, WNR and Distinguished University Professor, Department of Biological and Agricultural Engineering, North Carolina State University, Raleigh, North Carolina; Jay D. Atwood, Agricultural Economist, USDA Natural Resources Conservation Service, Resource Inventory and Assessment Division, Temple, Texas; and Jane R. Frankenberger, ASABE Member Engineer, Associate Professor, Department of Agricultural and Biological Engineering, Purdue University, West Lafayette, Indiana. Corresponding author: X. Wang, Texas Agricultural Experiment Station, Blackland Research and Extension Center, 720 East Blackland Rd., Temple, Texas 76502; phone: 254-774-6105; fax: 254-774-6145; e-mail: [email protected].

successfully tested using a six-year data set from an artificially drained agricultural field research site in southeastern Indiana (Youssef et al., 2005b). The model can be used to study the effects of drainage design and management, cropping systems, N fertilization, and tillage on C and N transformations in the soil-water-plant system and subsequently on N losses from drained croplands under a wide range of soil types and climatic conditions. Thus, it can be a useful tool for the development and evaluation of methods and management practices for sustainable agricultural systems on drained lands. The model, however, requires a relatively large number of input parameters. Like any other model, proper parameterization of DRAINMOD-N II is essential for making use of its modeling capabilities. In order to facilitate model parameterization, a sensitivity analysis is necessary to rank various model inputs according to their influence on model predictions. In practical model applications, users should make measurements, if feasible, to accurately estimate the most influential model inputs, while other input parameters can be estimated from the literature. Sensitivity analysis for models having large numbers of input parameters usually starts with parameter screening to exclude less influential inputs from further analysis. The LH-OAT sensitivity analysis method (Griensven and Meixner, 2003), which combines the Latin hypercube (LH) sampling technique (McKay et al., 1979) with the traditional one-factor-at-a-time (OAT) methods (Morris, 1991), is an efficient tool for parameter screening. The LH sampling technique is commonly used in stochastic water quality models, which usually require large numbers of input parameters (Weijiers and Vanrolleghem, 1997; Vandenberghe et al., 2001), due to its efficiency and robustness. The LH-OAT method provides a simple and

Transactions of the ASAE Vol. 48(6): 2205−2212

2005 American Society of Agricultural Engineers ISSN 0001−2351

2205

effective way to assess the influence of a parameter across its full range of values while varying the values of the other parameters. The method, however, does not address the variance in the output due to the variance in model inputs. Variance-based sensitivity analysis methods, e.g., the Fourier amplitude sensitivity test (FAST), are more reliable for analyzing non-linear and non-monotonic water quality models (Melching and Yoon, 1996). The application of variance-based methods allows the determination of the individual effect of input factors as well as the interactive effect among factors (Ratto et al., 2001). The main objective of this study was to conduct sensitivity analyses for DRAINMOD-N II using the LH-OAT and the extended FAST techniques to identify model inputs that have the greatest impact on model predictions of N losses from drained croplands. The results will help potential users of the model make informed decisions about model parameterization.

MATERIALS AND METHODS DESCRIPTION OF DRAINMOD-N II DRAINMOD-N II (Youssef, 2003) is a field-scale, process-based model that simulates C and N dynamics and turnover in drained agricultural lands under a wide range of soil types, management practices, and environmental conditions. The model is an enhanced version of the nitrogen model DRAINMOD-N (Brevé et al., 1997). A brief description of the model is given in this section. The model is described in detail by Youssef et al. (2005a). DRAINMOD-N II simulates C and N cycles in the soil-water-plant system and operates at different levels of complexity. The C cycle considered in the model is similar to that is used in the CENTURY model (Parton et al., 1987; Parton et al., 1993) and consists of three soil organic matter (SOM) pools (active, slow, and passive), two above- and below-ground litter pools (metabolic and structural), and a surface microbial pool. Each pool is characterized by its OC content, potential rate of decomposition, and carbon-to-nitrogen (C:N) ratio. The N cycle considered in the model includes three N pools: nitrate-nitrogen (NO3-N), ammoniacal-nitrogen (NHx-N), and organic-nitrogen (ON). The model simulates processes of atmospheric deposition, application of mineral N fertilizers, soil amendment with ON sources (plant residues/animal manure), N plant uptake, OC decomposition and associated N mineralization/immobilization, nitrification, denitrification, ammonia (NH3) volatilization, and NO3-N and NHx-N losses with subsurface drainage and surface runoff (Youssef et al., 2005a). The model simulates N reactive transport using a finite difference solution to a multi-phase form of the one-dimensional advection-dispersion-reaction equation. Driving hydrologic parameters of DRAINMOD-N II are based on the output of the water management model DRAINMOD 5.1 (Skaggs, 1978, 1999). DRAINMOD-N II simulates N fertilizer application and associated short-term processes including fertilizer dissolution, urea hydrolysis, and soil pH change caused by the application of NH4-forming fertilizers such as urea and anhydrous NH3 (Youssef et al., 2005a). It simulates the effects of plant residue management and tillage practices on C and N dynamics using a tillage intensity factor that reflects the level of soil disturbance caused by tillage operations. It

2206

also simulates the application of animal waste to the soil-plant system. Mineral and organic N portions of the animal waste are handled by the mineral N fertilizer and the organic matter components of the model, respectively. In DRAINMOD-N II, the potential N uptake for each crop is computed (Youssef et al., 2005a) and empirically distributed over the days of the growing season (Brevé et al., 1997). Both NO3 and NH4 in the root zone are assumed to be equally available to plants. If one form of N is consumed before satisfying plant needs, the rest is taken from the other form. If total mineral N does not meet the uptake demand and the plant is not a legume, actual N uptake is set equal to the available amount of mineral N (Youssef et al., 2005a). Nitrogen mineralization/immobilization is modeled in DRAINMOD-N II as a consequence of C cycling in the system during OM decomposition. Carbon transformations among individual pools are assumed to follow first-order kinetics (Parton et al., 1987; Parton et al., 1993). The model simulates the nitrification and denitrification processes using Michaelis-Menten kinetics. The effect of nitrification inhibitors on the nitrification process is simulated using a response function that modifies the process rate according to inhibitor concentration. The influence of OC decline with depth on the denitrification process is implicitly simulated using an exponential depth function that reduces potential process rate with soil depth (Youssef et al., 2005a). The effect of environmental factors, including soil temperature, soil moisture, and soil pH, on C and N biological transformations is modeled by defining, for each factor, a dimensionless response function. An overall environmental response function, composed of a linear combination of the individual response functions, is used to modify optimum process rates for non-optimum environmental conditions (Youssef et al., 2005a). DRAINMOD-N II simulates temporal changes in soil pH using a simplified approach in order to determine the form of NHx-N in the soil profile. It calculates the amounts of OH-ion added to or taken from the system due to fertilizer application, nitrification, and N plant uptake and estimates Table 1. Summary of simulated site conditions. Climatic Parameters Average annual precipitation 1300 mm Average annual potential ET 900 mm Average annual temperature 15.5°C Properties of top soil pH Bulk density Organic carbon content Water content at welting point Water content at saturation Lateral saturated hydraulic conductivity Drainage system Drain depth Drain spacing Surface depressional storage

7 1.4 g cm−3 4% 0.12 cm3 cm−3 0.41 cm3 cm−3 3.8 cm h−1 115 cm 20 m 0.5 cm

Agronomic practices Tillage and seedbed preparation Corn planting N fertilizer application of 50 kg N ha−1 N fertilizer application of 100 kg N ha−1 Corn harvesting

10 April 15 April 15 April 1 June 25 August

TRANSACTIONS OF THE ASAE

the resulting pH change using a soil buffer factor (Youssef et al., 2005a). Model inputs include soil properties, crop parameters, management parameters, and C and N processes and transformations parameters. Model output includes NO3-N and NHx-N concentrations in soil solution and drain flow, OC content of the top 20 cm soil layer, and cumulative rates of simulated N processes.

SENSITIVITY ANALYSIS The sensitivity analysis of DRAINMOD-N II was conducted for a continuous corn production on an artificially drained Portsmouth sandy loam soil using a 40-year period (1951-1990) of climatological record at the Tidewater Experiment Station near Plymouth, North Carolina (table 1). Preliminary simulations were conducted to study the effects of changing model inputs during the sensitivity analysis on

Table 2. DRAINMOD-N II input parameters and their ranges considered in the LH-OAT sensitivity analysis. Uniform Distribution Minimum Value

Maximum Value

Crop (corn)

Harvest index, HI Root-to-shoot ratio, RSR Grain N content, Ngrain (%) Shoot N content, Nshoot (%) Root N content, Nroot (%) Shoot lignin content, Lshoot (%) Root lignin content, Lroot (%)

0.44 0.10 1.2 0.4 0.4 4 6

0.56 0.22 1.8 1.0 1.0 12 15

Parameter

Management

Tillage intensity factor, TI

0.1

1.0

Transport

Longitudinal dispersivity, λ (cm) Tortuosity, τ

5 0.3

30 0.7

pH and NH3 volatilization

Maximum soil buffering capacity, BCmax (pH g mole−1) Empirical factor for resistance to NH3 transport from soil surface to atmosphere, γ (s cm−1)

1 × 104 3

20 × 104 1000

Nitrification

Maximum reaction rate, Vmax-nit (µg g−1 d−1) Half saturation constant, Km-nit (µg g−1) Optimum temperature, Topt-nit (°C) Empirical shape factor, βnit Value of soil water function at wilting point, fwp-nit Value of soil water function at saturation, fsat-nit Empirical exponent for soil water function, enit

4.0 1.0 19.0 0.272 0.1 0.0 0.5

38.0 180.0 33.0 0.676 0.4 0.2 2.0

Denitrification

Maximum reaction rate, Vmax-den (µg g−1 d−1) Half saturation constant, Km-den (mg L−1) Optimum temperature, Topt-den (°C) Empirical shape factor, βden Empirical exponent of SOM depth function, α Relative saturation below which denitrification ceases, sden Empirical exponent for soil water function, eden

2.0 1.0 18.5 0.093 0.2 0.5 1.5

60.0 60.0 74.0 0.372 0.5 0.9 2.5

Fertilizer dissolution

Threshold soil water content below which dissolution ceases, θdis (cm3 cm−3) Zero-order dissolution rate coefficient, Kdis (d−1)

0.120 0.5

0.276 2.0

Maximum reaction rate, Vmax-hyd (µg g-1 d-1) Half saturation constant, Km-hyd (mg L-1) Optimum temperature, Topt-hyd (°C) Empirical shape factor for temperature function, βhyd Value of soil water function at wilting point, fwp-hyd Value of soil water function at saturation, fsat-hyd Empirical exponent for soil water function, ehyd

80.0 30.0 48.0 0.061 0.3 0.6 0.5

4500.0 1750.0 54.0 0.204 1.0 1.0 2.0

5.3 × 10−3 2.0 × 10−2 8.2 × 10−3 6.7 × 10−3 2.5 × 10−2 1.0 × 10−2 2.7 × 10−4 6.2 × 10−6 18.5 0.093 0.1 0.4 0.5

21.4 × 10−3 8.1 × 10−2 32.9 × 10−3 26.9 × 10−3 10.1 × 10−2 4.0 × 10−2 11.0 × 10−4 24.7 × 10−6 74.0 0.372 0.4 0.7 2.0

Urea hydrolysis

OC decomposition

Surface structural litter pool, Kdec1 Surface metabolic litter pool, Kdec2 Surface microbial pool, Kdec3 Potential rate of Below−ground structural litter pool, Kdec4 decomposition (d−1) for: Below−ground metabolic litter pool, Kdec5 Active SOM pool, Kdec6 Slow SOM pool, Kdec7 Passive SOM pool, Kdec8 Optimum temperature, Topt-dec (°C) Empirical shape factor for temperature function, βdec Value of soil water function at wilting point, fwp-dec Value of soil water function at saturation, fsat-dec Empirical exponent for soil water function, edec

Vol. 48(6): 2205−2212

2207

soil organic matter. A set of long-term (50 years) simulations was conducted using minimum and maximum values of parameters that most likely affect soil organic matter. Simulation results indicated that the shortest period needed for the system to reach a new “quasi” state of equilibrium is 20 years. Therefore, it was decided to conduct the sensitivity analysis by running 40−year simulations and only using model output during the last 20 years of each simulation for the analysis. A total of 48 input parameters (table 2) were selected for a first step screening using the LH-OAT sensitivity analysis method (Griensven and Meixner, 2003). The effects of these parameters on model prediction of NO3-N losses via subsurface drainage were investigated using both yearly outputs and long-term (20 years) average outputs. LH-OAT SENSITIVITY ANALYSIS A modified version of the public domain FORTRAN code of the LH-OAT sensitivity analysis method developed by Griensven and Meixner (2003) was used for DRAINMOD-N II parameter screening. The code was modified to calculate relative sensitivity instead of absolute sensitivity in order to provide a basis for comparing various parameters. The code was also customized to: (1) update the DRAINMOD-N II input file for each set of inputs, (2) execute DRAINMOD-N II simulations, and (3) read yearly outputs, calculate average values, and write the output into a file formatted for the LH-OAT code. The LH-OAT method uses the LH samples as initial points for the one-factor-at-a-time (OAT) method. The LH is a stratified sampling technique, in which the range of each input parameter is divided into N intervals of equal probability, 1/N. Then, one sample is randomly generated within each interval, resulting in a total of N non-overlapping samples for each input parameter. In this study, each of the input parameters was divided into N = 10 intervals of equal probability. The parameters were assumed to be triangularly distributed in order to make the intervals around the means narrower. Therefore, a total of 10 non-overlapping samples for each of the 48 input parameters were generated. Each LH sample set, j, was then used as an initial point for the OAT design. In the OAT method, only one factor, xi , varies at a time while other factors are fixed. The change in model output can then be unambiguously attributed to such a change in the factor xi . A relative sensitivity index, defined as the ratio between the relative normalized change in output to the normalized change in related input, was calculated to facilitate a direct comparison and to avoid difficulties concerning the different orders of magnitude in input parameters (Brunner et al., 2004). The relative sensitivity, Rs , (eq. 1) developed by McCuen (1973) was slightly modified (eq. 2) to consider the absolute change in model output and related input: ∆O x Rs = ⋅ (1) ∆x O Sij =

O( x1 ,..., xi + ∆xi ,..., x p ) − O( x1 ,..., xi ,..., x p ) O( x1 ,..., xi + ∆xi ,..., x p ) + O( x1 ,..., xi ,..., x p )



xi (2) ∆xi

where Sij is the relative partial effect of parameter xi around the LH point j, p is the total number of parameters (equals 48 in this study),and O is the model output. Parameter xi is ran-

2208

domly increased or decreased by  xi = fxi , where f is a dimensionless fraction. In this study, f was set equal to 0.05 for all parameters. The final effect (sensitivity index) of xi is calculated by averaging these partial effects: N

∑ Sij

S xi =

j =1

(3)

N

VARIANCE-BASED SENSITIVITY ANALYSIS A detailed description of variance-based sensitivity analysis can be found in Saltelli et al. (2000a). Variancebased sensitivity analysis methods are based on dividing the variance of model output into partial variances depending on single factors and factor interactions of increasing order (Ratto et al., 2001). Variance-based sensitivity indices are estimated as ratios between conditional and unconditional variances for model output (Schwieger, 2004; Ratto et al., 2001): σ 2E (O| X ) i (4) Si = 2 σO where Si is called the correlation ratio, importance measure, 2 is the total variance of or first-order sensitivity index, σO 2 model output O, and σ E (O| X ) is the variance of the expected i

value of model output O by holding Xi fixed. The expected value, E (O | X i ) , has to be evaluated over the whole interval of Xi in order to obtain a global sensitivity measure (Schwieger, 2004). The second-order sensitivity index, Si,j , represents the sensitivity of O to the interaction between parameters Xi and Xj (Ratto et al., 2001; Schwieger, 2004): S i, j = where σ 2E (O| X

σ 2E (O| X

i,X j)

i,X j)

− σ 2E (O| X ) − σ 2E (O| X i

j)

(5)

2 σO

is the variance of the expected value of

model output O by holding Xi and Xj fixed, and σ 2E (O| X

j)

is

the variance of the expected value of model output O by holding Xj fixed. Higher-order terms can be defined in a similar way to equation 5. Obtaining all higher-order terms, however, is computationally expensive and might be impracticable, especially for complex models (Schwieger, 2004; Ratto et al., 2001). Alternatively, the total sensitivity index, STi , can be used to indicate the overall impact of parameter Xi on model output Y through all the interactions involving Xi . The value of STi can be expressed as (Saltelli et al., 2000a; Schwieger, 2004): σ 2E (O| X ) ~i (6) STi = 1 − 2 σO where X~i indicates all the parameters, but Xi , σ 2

E (O| X ~i )

is the

variance of the expected value of model output O by holding all input parameters other than Xi fixed. Different methods have been used for estimating Si and STi with special sampling schemes, such as the extended Fourier

TRANSACTIONS OF THE ASAE

amplitude sensitivity test (FAST) (Cukier et al., 1973; Saltelli et al., 2000a; Saltelli et al., 1999) and Sobol’s method (Sobol, 1993). The extended FAST technique has the advantage of a small sample size in comparison to Sobol’s method (Schwieger, 2004). In this study, the extended FAST method was used to determine the first-order and total-order sensitivity indices. The FAST method is based on performing numerical calculations to obtain the expected value and variance of a model output. It uses a Fourier transformation to convert a multidimensional integral over all the uncertain model inputs into a one-dimensional integral (Saltelli and Bolado, 1998; Saltelli et al., 1999). A detailed description of the extended FAST method can be found in Saltelli et al. (1999). The SIMLAB software, developed by the Joint Research Centre of the European Commission (SIMLAB, 2003), provides an interface for the extended FAST sampling. A total of 2,000 parameter sets were generated using SIMLAB from the given ranges and distribution specified in table 2. A FORTRAN program was written to: (1) update parameter values using the generated FAST samples, (2) execute DRAINMOD-N II simulations, and (3) provide output files formatted for FAST sensitivity analysis. The SIMLAB Statistical Post Processor was used to calculate FAST sensitivity indices.

RESULTS AND DISCUSSION The 48 model input parameters considered in the LH-OAT screening step were ranked according to their influence on model output (table 3). Results indicated that DRAINMODN II is most sensitive to inputs that control the denitrification process. The optimum temperature for denitrification, Topt-den , and the empirical shape factor, den , of the temperature response function for denitrification ranked first and second, respectively. Michaelis-Menten parameters for denitrification had high rankings as well. The maximum reaction rate, Vmax-den , ranked fourth, and the half saturation constant, Km-den , for denitrification ranked seventh. The empirical exponent of the depth function, , that characterizes the decline in the optimum denitrification rate with depth came fifth, and the threshold relative saturation, sden , at which denitrification starts ranked sixth. Model inputs that parameterize the temperature response function for OC decomposition were also among the high-ranking parameters. The optimum temperature ranked third, and the empirical shape factor had the eighth rank. Then followed three nitrification parameters: Km-nit ranked ninth, Vmax-nit ranked tenth, and Topt-nit ranked eleventh. The least influential inputs were urea hydrolysis parameters (ranks 39, 41-44, 46, 47), OC decomposition parameters for surface litter pools, and the empirical resistance factor,  (rank 48). The yearly sensitivity indices for the denitrification parameters and OC decomposition parameters are plotted in figures 1 and 2, respectively, which illustrate that sensitivities are dynamic in the temporal dimension, as indicated by Ratto et al. (2001). In the group of the denitrification parameters, Topt-den and den were the dominant parameters for the whole time period (fig. 1), while Topt-dec and dec were the dominant parameters in the group of OC decomposition parameters (fig. 2). Changes in sensitivity from year to year occurred because climate variables, temperature and precipitation, are model inputs that interact with other model parameters considered in the sensitivity analysis.

Vol. 48(6): 2205−2212

Table 3. Sensitivities of average annual nitrate-N losses to DRAINMOD-N II parameters using LH-OAT. S xi Parameter Rank Topt-den βden Topt-dec Vmax-den α Sden Km-den βdec Km-nit Vmax-nit Topt_nit eden Kdec7 λ Ngrain Lshoot HI Nshoot Nroot RSR enit Kdec6 Kdec8 fsat-dec edec Kdec4 Kdec5 θdis BCmax Kdis fsat-nit TI Kdec2 Lroot fwp-dec fwp-nit τ fwp-hyd Vmax-hyd Kdec3 fsat-hyd βhyd βnit ehyd Kdec1 Topt-hyd Km-hyd γ

1.365 0.677 0.505 0.486 0.472 0.426 0.408 0.233 0.185 0.146 0.113 0.112 0.105 0.100 0.076 0.064 0.062 0.056 0.052 0.046 0.045 0.039 0.035 0.033 0.029 0.028 0.028 0.023 0.022 0.022 0.019 0.018 0.018 0.015 0.014 0.010 0.008 0.007 0.006 0.006 0.006 0.004 0.004 0.003 0.003 0.002 0.001 0.000

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48

A total of 20 parameters were selected for further analysis using the extended FAST technique. They were the parame− ters ranked 1 to 14 by the LH-OAT method, five OC decomposition parameters, and the empirical shape factor of the temperature response function for the nitrification process. Table 4 lists the ranking of the parameters based on the FAST technique, for both first order and total order. The FAST first order and the LH-OAT method calculate the main influence (singular influence). Parameter rankings according to the FAST first-order sensitivity indices were similar to those obtained by the LH-OAT method (tables 3 and 4). The 14 most influential parameters identified by both methods differed only in two parameters. Michaelis-Menten parame−

2209

5 Topt−den

b den Vmax −den a Sden Km−den eden

LH−OAT sensitivity index

4.5 4 3.5 3 2.5 2 1.5 1 0.5 0

1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990

Year Figure 1. Sensitivities of DRAINMOD-N II predicted yearly nitrate-N losses to denitrification parameters using LH-OAT. 5

LH−OAT sensitivity index

4.5

Topt−dec Topt−dec βdec bdec Kdec7 K dec7 Kdec6 K dec6 Kdec8 K dec8 Kdec4 Kdec4 K Kdec5 dec5 K dec2 Kdec2

4 3.5 3 2.5 2 1.5 1 0.5

0 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987 19881989 1990

Year Figure 2. Sensitivities of DRAINMOD-N II predicted yearly nitrate-N losses to OC decomposition parameters using LH-OAT. Table 4. Sensitivities of average annual nitrate-N losses to DRAINMOD-N II parameters using extended FAST technique. FAST First Order FAST Total Order Parameter Topt-den βden Topt-dec α Vmax-den Sden eden βdec Kdec7 Km-den Kdec4 Kdec6 λ Topt_nit Km-nit Vmax-nit Kdec8 Kdec5 Kdec2 βnit

2210

Si

Rank

STi

Rank

0.624 0.090 0.087 0.044 0.042 0.036 0.020 0.003 0.003 0.003 0.002 0.002 0.001 0.001 0.001 0.001 0.000 0.000 0.000 0.000

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

0.702 0.141 0.154 0.075 0.064 0.061 0.091 0.022 0.02 0.016 0.051 0.076 0.031 0.019 0.013 0.041 0.016 0.015 0.016 0.014

1 3 2 6 7 8 4 12 13 15 9 5 11 14 20 10 16 18 17 19

ters for the nitrification process, Vmax-nit holding the 9th rank and Km-nit holding the 10th rank according to LH-OAT, were moved back to the 15th and 16th ranks, respectively. FAST advanced potential rates of decomposition for the below-ground structural litter pool (Kdec4 , ranking 26 according to LH-OAT) and the active SOM pool (Kdec6 , ranking 22 according to LHOAT) to ranks 11 and 12, respectively. Almost half of the 14 most influential parameters identified by both methods hold the same rank (tables 3 and 4). This similarity between the LH-OAT and the FAST first-order rankings indicates that both methods are robust in identifying the key parameters with main effect. Both methods indicated that DRAINMOD-N II is most sensitive to the parameters controlling the denitrification process. The parameters controlling OC decomposition and associated N mineralization/immobilization came after the denitrification parameter. It should be mentioned that the relative importance of the parameters controlling the two processes depends on the conditions of the system being simulated. The model is expected to be far more sensitive to the denitrification parameters for agricultural systems, like the simulated corn production, in which N fertilizers are the dominant source of mineral N. On the other hand, model sensitivity to OC decomposition parameters is expected to increase for systems that rely more on OC decomposition as

TRANSACTIONS OF THE ASAE

FAST total order sensitivity index

0.8 0.7 Topt−den T opt−den Topt−dec T opt−dec bden b den eeden den Kdec6 K dec6 a Vmax−den V max−den SSden den

0.6 0.5 0.4 0.3 0.2 0.1

0 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990

Year Figure 3. FAST total-order sensitivities of yearly nitrate-N losses to the 20 DRAINMOD-N II parameters.

a source of mineral N (e.g., soils rich in organic matter, and crop rotations that include legumes). The results obtained in this study are consistent with results of sensitivity and uncertainty analyses previously conducted for the original version of the model (Breve et al., 1997; Haan and Skaggs, 2003). Using a simple one-at-a time approach, Brevé et al. (1997) showed that model prediction of NO3-N loss is most sensitive to the denitrification rate coefficient, followed by the net N mineralization rate coefficient. Haans and Skaggs (2003) estimated that the uncertainty in the denitrification rate coefficient accounts for more than 70% of the uncertainty in model predictions of NO3-N drainage losses. The extended FAST total-order sensitivity indices reveal the overall influence, which includes both additive and interaction influence. Model parameters interact within the linked system of equations of the model to produce simulation outputs. Interaction occurs where the parameter effect on output depends on the values that the other parameters may take. The interaction terms usually grow with the number of parameters and with the range of variation of the parameters (Saltelli et al., 2000b). The increase in FAST total-order sensitivity indices over first-order sensitivity indices reveals the influence of the interaction among parameters. The parameter with small main effect (first-order sensitivity index) but high total effect, for example, the potential rate of decomposition for the active SOM pool (Kdec 6), influence the model output mainly through interaction. Although the variance-based global sensitivity analysis provides quantitative criteria for choosing “leading” parameters based on main and total effect, such criteria do not necessarily provide a direct, complete representation of the interaction (Ratto et al., 2001). The first eight ranked parameters according to the FAST total-order sensitivity indices are shown in figure 3 for their yearly sensitivity indices. In general, the coefficients of variation of the extended FAST total-order sensitivity indices were smaller than those of the FAST first-order sensitivity indices and those of LH-OAT sensitivity indices (table 5). This indicates that the extended FAST total-order sensitivity measure, which takes into account the total effect of a parameter with respect to all other parameters, is more stable.

Vol. 48(6): 2205−2212

Parameter Topt-den Topt-dec βden eden Kdec6 α Vmax-den Sden

Table 5. CV (%) of sensitivities of yearly nitrate-N losses to important parameters. LH-OAT FAST First Order FAST Total Order 19.0 26.3 15.2 22.3 42.0 21.7 18.0 19.0

8.2 26.4 15.0 31.8 38.2 20.8 14.8 42.7

5.2 24.3 9.4 16.1 8.4 17.6 9.3 36.0

Although DRAINMOD-N II has a fairly large number of inputs, the model is most sensitive to only a few of them (tables 3 and 4). This means that modelers can exclude a large number of model inputs from the calibration procedure and focus only on the few important ones.

CONCLUSIONS A two-step sensitivity analysis was conducted for the nitrogen simulation model DRAINMOD-N II. Forty-eight model input parameters were selected for a first-step screening using the LH-OAT method, from which 20 parameters were identified for further analysis using a variancebased sensitivity analysis technique, the extended FAST. DRAINMOD-N II simulated a continuous corn production on a subsurface drained sandy loam using a 40-year eastern North Carolina climatological record. The sensitivity indices for the yearly and the 20-year average annual nitrate-N losses through drain flow were calculated using both sensitivity analysis methods. The use of different sensitivity analysis methods with different theoretical foundations increases the confidence in the results. Both sensitivity analysis methods indicated in these conditions that DRAINMOD-N II was most sensitive to denitrification parameters, especially those controlling temperature effect on process rate. Results also indicated that the model was mildly sensitive to the parameters controlling OC decomposition and associated N mineralization/immobilization. Although not proven by this study, the relative importance between the two influential sets of inputs is expected to vary according to the local conditions of the agricultural system being simulated.

2211

This study demonstrated a widely applicable procedure for conducting sensitivity analysis for hydrologic/water quality models with large input requirement that uses the efficient LH-OAT method and a variance-based sensitivity analysis technique. This procedure can effectively identify the significant and less significant parameters during the screening step. It can also reveal parameter influence resulting from interaction with other parameters even though the parameter itself might not be influential individually. This is particularly important since the selection of input parameters for model calibration should be based on the total effect of the parameter, which accounts for both the individual parameter effect as well as the effect due to parameter interaction.

REFERENCES Brevé, M. A., R. W. Skaggs, J. E. Parsons, and J. W. Gilliam. 1997. DRAINMOD-N, A nitrogen model for artificially drained soils. Trans. ASAE 40(4): 1067-1075. Brunner, A. C., S. J. Park, G. R. Ruecker, R. Dikau, and P. L. G. Vlek. 2004. Catenary soil development influencing erosion susceptibility along a hillslope in Uganda. CATENA 58(1): 1-22. Cavero J., R. E. Plant, C. Shennan, D. B. Friedman, J. R. Williams, J. R. Kiniry, and V. W. Benson. 1999. Modeling nitrogen cycling in tomato-safflower and tomato-wheat rotations. Agric. Systems 60(2): 123-135. Cukier, R. I., C. M. Fortuin, K. E. Shuler, A. G. Petschek, and J. H. Schaibly. 1973. Study of the sensitivity of coupled reaction systems to uncertainties in rate coefficients: Part I. Theory. J. Chem. Physics 59(8): 3873-3878. Gilliam, J. W., J. L. Baker, and K. R. Reddy. 1999. Water quality effects of drainage in humid regions. In Agricultural Drainage, 801-830. Agronomy Monograph 38. R. W. Skaggs and J. van Schilfgaarde, eds. Madison, Wisc.: ASA, CSSA, and SSSA. Griensven, A., and T. Meixner. 2003. LH-OAT sensitivity analysis tool. Tucson, Ariz.: University of Arizona, Department of Hydrology and Water Resources. Available at: www.sahra.arizona.edu/software/index_main.html. Accessed 26 November 2004. Haan, P. K., and R. W. Skaggs. 2003. Effect of parameter uncertainty on DRAINMOD predictions: II. Nitrogen loss. Trans. ASAE 46(4): 1069-1075. McCuen, R. H. 1973. The role of sensitivity analysis in hydrologic modeling. J. Hydrol. 18(1): 37-53. McKay, M. D., R. J. Beckman, and W. J. Conover. 1979. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics 21(2): 239-245. Melching, C. S., and C. G. Yoon. 1996. Key sources of uncertainty in QUAL2E model of Passaic river. ASCE J. Water Res. Planning and Mgmt. 122(2): 105-113. Morris, M. D. 1991. Factorial sampling plans for preliminary computational experiments. Technometrics 33(2): 161-174. Parton, W. J., D. S. Schimel, C. V. Cole, and D. S. Ojima. 1987. Analysis of factors controlling soil organic matter levels in Great Plains grasslands. SSSA J. 51(5): 1173-1179.

2212

Parton, W. J., M. O. Scurlock, D. S. Ojima, T. G. Gilmanov, R. J. Scholes, D. S. Schimel, T. Kirchner, J.-C. Menaut, T. Seastedt, E. G. Moya, A. Kamnalrut, and J. I. Kinyamario. 1993. Observations and modeling of biomass and soil organic matter dynamics for the grassland biome worldwide. Global Biogeochem. Cycles 7(4): 785-809. Ratto, M, S. Tarantola, and A. Saltelli. 2001. Sensitivity analysis in model calibration: GSA-GLUE approach. Computer Physics Comm. 136(3): 212-224. Saltelli A., and R. Bolado. 1998. An alternative way to compute Fourier Amplitude Sensitivity Test (FAST). Computational Stat. and Data Analysis 26(4): 445-460. Saltelli, A., S. Tarantola, and K. Chan. 1999. A quantitative, model independent method for global sensitivity analysis of model output. Technometrics 41(1): 39-56. Saltelli, A., K. Chan, and E. M. Scott, ed. 2000a. Sensitivity Analysis. Chichester, U.K.: John Wiley and Sons. Saltelli, A., S. Tarantola, and F. Campolongo. 2000b. Sensitivity analysis as an ingredient of modeling. Stat. Sci. 15(4): 377-395. Schwieger, V. 2004. Variance-based sensitivity analysis for model evaluation in engineering surveys. Presented at INGEO 2004 and FIG Regional Central and Eastern European Conference on Engineering Surveying. Available at: www.fig.net/pub/bratislava/papers/ts_01/ts_01_schwieger.pdf. Accessed 20 January 2005. SIMLAB. 2003. Manual of SIMLAB Ver 2.2. European Commission, Joint Research Centre. Available at: http:/webfarm.jrc.cec.eu.int/uasa. Accessed 20 March 2004. Singh, P., and R. S. Kanwar. 1995. Simulating NO3-N transport to subsurface drain flows as affected by tillage under continuous corn using modified RZWQM. Trans. ASAE 38(2): 499-506. Skaggs, R. W. 1978. A water management model for shallow water table soils. Report 134. Raleigh, N.C.: North Carolina State University, North Carolina Water Resources Research Institute. Skaggs, R. W. 1999. Drainage simulation models. In Agricultural Drainage, 469-500. Agronomy Monograph 38. R. W. Skaggs and J. van Schilfgaarde, eds. Madison, Wisc.: ASA, CSSA, and SSSA. Sobol, I. M. 1993. Sensitivity estimates for nonlinear mathematical models. Math. Modeling and Computational Experiments 1(4): 407-414. Vandenberghe, V., A. Van Griensven, and W. Bauwens. 2001. Sensitivity analysis and calibration of the parameters of ESWAT: Application to the River Dender. Water Sci. and Tech. 43(7): 295-301. Weijiers, S. R., and P. A. Vanrolleghem. 1997. A procedure for selecting parameters in calibrating the activated sludge model No. 1 with full-scale plant data. Water Sci. and Tech. 36(5): 69-79. Youssef, M. A. 2003. Modeling nitrogen transport and transformations in high water table soils. Unpublished PhD diss. Raleigh, N.C.: North Carolina State University, Department of Biological and Agricultural Engineering. Youssef, M. A., R. W. Skaggs, G. M. Chescheir, and J. W. Gilliam. 2005a. The nitrogen simulation model, DRAINMOD-N II. Trans. ASAE 48(2): 611-626. Youssef, M. A., E. J. Kladivko, R. W. Skaggs, X. Wang, and J. R. Frankenberger. 2005b. Field-testing of DRAINMOD-N II for an Indiana silt loam. ASAE Paper No. 052023. St. Joseph, Mich.: ASAE.

TRANSACTIONS OF THE ASAE