Using Response Surface Methodology for Economic and Environmental Trade-offs at the Farm Level

Abstract

United States farmers typically spend over $10 billion annually on commercial fertilizer. Chemical inputs such as nitrogen (N) are essential for maintaining crop yields; however, farmers often apply excessive N inputs as an insurance policy. Nitrogen fertilizer consumption in the U.S. quadrupled from 3 million metric tons in 1961 to over 12 million metric tons in 2004, and per ha N fertilizer use quadrupled. Increase in N use has been associated with the impairment of U.S. streams, lakes, and aquifers. The objective of this research study was to develop an integrated farm-level economic/environmental risk framework for trade-off analysis between farm profitability and environmental externalities (impacts). Results indicated that there was no single point of optimal trade-off between farm profitability and the environment. Additionally, trade-offs between farm profit and environmental impacts varied significantly depending on the choice of cropping or tillage system.

Keywords

response surface method trade-off analysis environmental impacts nitrogen optimization

Introduction

Federal policies rely primarily on voluntary approaches by farmers and ranchers to minimize environmental externalities of agricultural production. Research has shown that in many cases the loss of agricultural chemicals from farm fields can be substantially reduced by conservation or best management practices (BMPs), and that significant improvements in water quality would result from widespread BMP adoption.¹ However, farmers face the challenge of balancing environmental considerations with farm profitability. Effective quantification of environmental versus economic trade-offs should result in tools that are based on the best available science and can be easily utilized by practitioners who may face time constraints. Comprehensive process-based agro-ecosystem models, such as the Nitrogen Loss and Environmental Assessment Package (NLEAP, Shaffer et al²) model and the Root Zone Water Quality Model 2 (RZWQM2, Ahuja et al³), have been used to accurately predict nitrogen (N) fate and crop yield in agricultural soils. However, direct use of these models is too difficult and time-consuming for conservation planners and land managers. Black box models have been developed based on statistical or empirical relations expressing the load of nitrate-N (NO₃-N) to groundwater or surface waters as a function of relatively few parameters such as land use, rainfall and soil type. The empirical relations are often developed from linear or non-linear regression analysis (eg, Simmelsgaard and Djurhuus,⁴ Hansen et al,⁵ Haberle et al,⁶ Pedersen et al⁷). The black box models are only applicable for simulating and predicting NO₃-N loads under similar conditions to those that they were calibrated under. Furthermore, they cannot be extrapolated to simulation of scenarios where considerable changes in climate, soils, management practices, are considered because the calibrated model relies fully on the observed data (Shepherd et al⁸). In summary, an opportunity exists for development of user-friendly tools that may lead to increased adoption of conservation and BMPs that both enhance water quality and maintain farm profitability.

Malone et al⁹ used 10 years of data (1994-2003) from a long-term study near Nashua, Iowa to develop multivariate polynomial regression equations predicting crop yield, NO₃-N concentration, drainage volume, and NO₃-N loss in subsurface drainage from corn (Zea mays L.) and soybean (Glycine max (L.) Merr.) crop rotations. The regression equations described over 87%, 85%, 94%, 76%, and 95% of variation in soybean yield, corn yield, subsurface drainage, NO₃-N concentration, and NO₃-N loss in subsurface tile drainage, respectively. While the above approaches concentrate strictly on environmental concerns, another approach to this problem has been to develop Pareto-efficient frontier tools (eg, Tóth et al¹⁰) that map trade-offs of economic gains (or losses) to environmental losses (or gains) for a range of conservation or BMPs. In this context, a management practice is Pareto-efficient if none of the associated environmental outputs or the associated cost can be improved (ie, increased for environmental outputs or decreased for costs) without compromising another output. Typically, these frontiers do not consider more than 1 environmental variable at a time because they are 2-dimensional (2-D) and all factors except economic ($) and environment (E) are fixed relative to $ and E. That is, input factors such as tillage type, soil type, and fertilizer application rate/timing are not considered. The notion of Pareto-optimality is critical because it helps towards finding management options that lead to BMPs in the most cost-efficient way possible.

This study utilizes response surface methodology (RSM)¹¹ to examine trade-offs when there are multiple E variables (ie, N measured in tile drainage flow and total N measured in the soil profile) and by allowing for factor inputs to remain variables. This research builds on previous work (eg, Downing et al,¹² Kelly et al¹³) by expanding trade-off curves to include 2 E variables instead of 1 and by improving the accuracy of (and information about) the Pareto-optimal frontier by simultaneously estimating the impacts of all input factors for each point on the frontier. That is, we estimate a 3-dimensional (3-D) frontier surface for $ and 2 E variables where the variables not seen in the frontier, like the type of tillage system or the N application rate, are simultaneously optimized with the surface. Burdick and Naylor¹⁴ argued that RSM is a useful alternative to classical optimization and mathematical programming techniques in the presence of experimental design data and that there is a need for an econometric model describing the economic system. They presented a mathematical framework for an economic constrained RSM model. Fitzgerald et al¹⁵ used a response surface approach to create ground sampling designs from input imagery in order to develop regression equations for predicting crop height and width attributes in a 3.4-ha cotton field. Predictions of height and width from regressions between the imagery and ground sampling at the calibration locations gave coefficients of determination (R²) for height ranging from 0.34 to 0.90 and width ranging from 0.30 to 0.94. All regression models but 1 were statistically significant at the α = 0.01 level. Lesch¹⁶ presented an objective sampling and simplified modeling strategy for predicting soil property information from spatially referenced sensor data. A model-based sampling strategy for estimating an ordinary linear regression model in the spatial setting was developed incorporating a traditional response surface design into an iterative, space-filling type algorithm for the purpose of selecting sample site locations. A case study confirmed the effectiveness of the model-based strategy over a more traditional simple random sampling strategy, and Lesch¹⁶ concluded that the strategy was highly effective at ensuring efficient regression model parameter estimates and reliable prediction maps. Frey and Patil¹⁷ summarized previous studies on using RSM and concluded it could be successfully used to represent the relationship between a response variable (output) and 1 or more explanatory inputs. They also showed that RSM could be used in a probabilistic analysis for identifying curvatures in the response surface and to account for higher-order effects. Dhungana et al¹⁸ used RSM in connection with the CERES-Wheat model and the HADCM2 climate simulation model to identify optimal configurations of plant traits and management practices that maximize yield of winter wheat in high CO₂ environments. They found that RSM used in conjunction with crop and climate simulation models was useful in understanding the complex relationships between wheat genotypes, climate, and management practices. Isukapalli et al¹⁹ compared a stochastic RSM with traditional Monte Carlo and Latin Hypercube Sampling (LHS) methods in a study of uncertainty propagation for environmental and biological systems. The traditional methods required a prohibitive number of model simulations compared to the stochastic response surface for a wide range of problems; however, the results obtained in all 3 methods agreed closely.

Nitrogen is essential to agricultural sustainability but poor N management may lead to adverse environmental impacts (eg, excessive leaching to groundwater).⁹ Clearly N application rate, type of tillage system, and soil physical properties significantly affect environmental (eg, N transport in the surface and subsurface) and economic (eg, crop yield and subsequently gross margin) trade-offs, but the combined or interactive effects of these variables is not clear. Accurate quantification of environmental and economic tradeoffs, as affected by multiple variables, is a first step toward developing relatively simple predictive management tools. Therefore, the primary objective of this research was to utilize RSM to develop an integrated farm-level economic/environmental risk framework for trade-off analysis between farm profitability and environmental externalities (impacts). We used 14 years (1990-2003) worth of economic budget and environmental data collected from 36 0.4 ha experimental plots at the Iowa State University Northeast Research and Demonstration farm near Nashua, Iowa. We consider maximization of economic profit and minimization of environmental impairment simultaneously by constraining the optimization problem at specific empirical limits for NO₃-N in tile drainage flow and total soil profile residual N. The Nashua plot data suggests NO₃-N leaching and/or corn and soybean yield was significantly affected by N application amounts, cropping system, tillage system, and precipitation.⁹ Other desirable features of the Nashua data set are that it has been thoroughly investigated with multiple peer-reviewed manuscripts, includes numerous N management treatments, is relatively long-term (1990-2003) with a range of yearly precipitation, and the soil associations present on the 0.4-ha experimental plots represent approximately 575,000 ha where corn or soybean is commonly grown in Iowa, USA.

Materials and Methods

Site Description and Management

Data for our study were obtained from 35 treatments on 36 plots (0.4 ha each) located at the Iowa State University Northeast Research Station near Nashua, IA, USA (43.0° N, 92.5° W). The experimental plots were established to quantify the impact of management practices on crop production and water quality.^20,21 The soils are predominantly Floyd loam (fine-loamy, mixed, mesic Aquic Hapludolls), Kenyon silty-clay loam (fine-loamy, mixed, mesic Typic Hapludolls) and Readlyn loam (fine-loamy, mixed, mesic Aquic Hapludolls) with 30-40 g kg^-1 (3%-4%) organic matter.²² These soils vary from moderately well to poorly drained, lie over loamy glacial till and belong to the Kenyon-Clyde-Floyd soil association. Soil slopes varied from 1% to 3% among the various plots. The field experiments were established on a 15 ha research site in 1977 using a randomized complete block design with 3 replications. The seasonal water table at the site fluctuates from 20-160 cm and subsurface drainage tubes/pipes (10 cm diameter) were installed in the fall of 1979 at 120 cm depth and 29 m apart. 3 experimental phases were conducted from 1978-1992, 1993-1998, and 1999-2003. From 1978-1992, there were 4 tillage treatments (chisel plow, moldboard plow, no-till, ridge-till) under 2 different cropping sequences (continuous corn and both phases of a corn-soybean rotation). Crop yield was the primary measurement from 1978-1989. Experimental data collected from 1990 included tile drain flow, NO₃-N concentration in tile drain flow, residual N in soil, crop yield, biomass and plant N uptake. From 1993-1998, there were 2 tillage treatments (chisel plow and no-till), with 8 N management treatments (eg, different rates, times of application, fertilizer type and/or swine manure) for chisel plow and 4 N treatments for no-till with no change in the number of crop sequences. The experimental data collected remained essentially the same as from 1990-1992 with the addition of runoff. Continuous corn was replaced with both phases of the corn-soybean rotation in 1999, and the experiments were continued along with 10 fertilizer and swine manure treatments in the chisel plow system and 2 swine manure treatments in the no-till system. All plots received swine manure and/or urea-ammonium-nitrate (UAN) fertilizer each cropping season, with the swine manure applied in either fall or spring using application rates based on N or phosphorus (P) needs for the corn-soybean/soybean-corn rotations. Experimental measurements from 1999-2003 again focused on tile drain flow, NO₃-N concentration in drain flow, soil N, and crop yield, biomass and N uptake. Similar to Malone et al,⁹ tillage treatments were combined in this study to simplify the analysis and reduce the number of variables in the developed response surface equations. However, predominate tillage indirectly affects N in the soil profile and subsurface drainage and therefore the crop yield.²⁰ Table 1A lists the major management practices for each of the 35 treatments (ie, tillage and cropping systems) from 1990-2003 for the Nashua, IA, experiment. The major management practices for each of the 36 plots from 1990-2003 are listed in Table 1B.

Table 1A.

Major management practices by treatment at the Northeastern Research and Demonstration Farm, Nashua, IA from 1990-2003.^#

Treatment no.	Treatment period	Cropping system	Tillage system	Treatment no.	Treatment period	Cropping system	Tillage system
1	1990-1992	CC	NT	19	1993-1998	CC	CP
2	1990-1993	CS	NT	20	1994-2003	CS	CP
3	1990-1992	SC	NT	21	1993-2003	SC	CP
4	1990-1992	CC	CP	22	2000-2003	CS	CP
5	1990-1993	CS	CP	23	2001-2003	SC	CP
6	1990-1992	SC	CP	24	1993-1998	CC	CP
7	1990-1992	CC	MP	25	1994-2003	CS	CP
8	1990-1992	CS	MP	26	1993-2003	SC	CP
9	1990-1992	SC	MP	27	1999	CC	CP
10	1990-1992	CC	RT	28	2000-2003	CS	CP
11	1990-1992	CS	RT	29	2000-2003	SC	CP
12	1990-1992	SC	RT	30	2000	CC	CP
13	1994-1998	CS	NT	31	2001-2003	CS	CP
14	1993-2000	SC	NT	32	2001-2003	SC	CP
15	1994-1999	CS	CP	33	2000-2003	CS	NT
16	1993-2000	SC	CP	34	2001-2003	SC	NT
17	1994-1999	CS	NT	35	1999-2000	SC	CP
18	1993-1998	SC	NT

Abbreviations: ^#CS, corn-soybean rotation with corn during even years; SC, soybean-corn rotation with corn during odd years; CC, continuous corn; CP, chisel plow; RT, ridge-till; MP, moldboard plow; NT, no-till.

Table 1B.

Major management practices by plot at the Northeastern Research and Demonstration Farm, Nashua, IA from 1990-2003.^#

Plot no.	Dominant soil type	Cropping system			Tillage system
		1978-1992	1993-1998	1998-2003	1978-1992	1993-1998	1998-2003
1, 7, 30	Readlyn, Kenyon	CS	CS	CS	CP	CP	CP
2, 16, 20	Readlyn, Kenyon	CS	CS	CS	MP	NT	NT
3, 24, 28	Readlyn, Kenyon	SC	SC	SC	NT	NT	CP
4, 18, 33	Kenyon	SC	CS	CS	CP	CP	CP
5, 21, 26	Readlyn, Kenyon	CC	CC	SC	CP	CP	CP
6, 32, 36	Readlyn, Kenyon	CC	SC	SC	RT	CP	CP
8, 9, 19	Readlyn, Floyd	CS	CS	CS	RT	CP	CP
10, 15, 29	Kenyon	CS	CS	CS	NT	NT	CP
11, 23,27	Kenyon	SC	SC	SC	RT	CP	CP
12, 17, 34	Kenyon, Floyd	SC	SC	SC	MP	CP	CP
13, 22, 35	Readlyn, Floyd	CC	CC	CS	MP	CP	CP
14, 25, 31	Readlyn, Kenyon	CC	SC	SC	NT	NT	NT

Economic Budgets

Economic budgets by plot for 1990-2003 (36 plots x 14 years = 504 plot-years) were developed as part of the web-based USDA Natural Resources Conservation Service (NRCS)—EconDoc exchange tool. Primary data sources for the study included both Nashua experimental records and USDA National Agricultural Statistical Services (NASS) published data. Variable and farm machinery inputs for each tillage system were based on actual management practices and equipment used at the Nashua agricultural experiment station. Similar to Williams et al,²³ equipment ages were assumed to be half of their depreciable life (the half-life ranged from 5 to 7 years old depending on equipment type), with equipment prices deflated to the appropriate year the machine was purchased to calculate the original value for depreciation. Variable costs included farm input costs and variable machinery costs. Farm input costs for tillage and variable inputs, such as seed, fertilizer, or chemicals, were assumed similar for each tillage system and estimated based on actual management practices obtained from the experiment stations. Variable machinery costs were based on equipment usage and included oil and lubrication (estimated at 20% of total fuel costs), fuel consumption, and repair and maintenance. Following Yiridoe et al,²⁴ machinery fixed costs including annual depreciation, interest on investment, and equipment storage, insurance, and repairs were allocated to each crop in the rotation based on usage (ie, if the same equipment was used for both corn and soybean production, the associated annual fixed cost was split between the crops). Machinery insurance and storage were based on 2% of the purchase price. Labor requirements were calculated using the time required for field operations with the equipment used for the operation. Additional details on the total cost of production for each tillage system for the period prior to the experimental phase analyzed in this study are described in Chase and Duffy.²⁵ To determine gross returns, we used average annual prices for corn and soybeans from NASS county data records and annual yields reported by the Nashua experiment station. Gross margins for each of the 4 tillage systems were then calculated by subtracting the operational costs from the corresponding gross returns. Gross margin represents the enterprise's contribution towards covering the fixed costs and generation of profit after operational costs have been covered.²⁶ The gross margin data were then averaged across the experimental replications and discounted to reflect the net present values.

Response Surface Problem Formulation

RSM was first introduced by Box and Wilson²⁷ and explores the relationships between several explanatory variables and one or more response variables. Response surface models are general linear models where attention is focused on characteristics of the fitted response function, particularly where optimum estimated response values occur. The main idea of RSM is to use carefully designed empirical data sets to obtain an optimal response surface. The major advantage of the RSM approach, with the availability of today's high-speed computational resources, is that a potentially computationally intensive model can be reduced to a simplified form that enables much faster model run times. The canonical analysis and eigenvectors analysis of the RSM model and the values of its coefficients may provide a useful indication of explanatory variable sensitivities.

A constrained RSM approach, including selected optimization algorithms (ie, the steepest descent or ascent) was applied to find the optimum surface regions of corn and soybean profitability from the Nashua, IA experimental plots subject to two constraints representing environmental externalities: NO₃-N measured in tile drainage and total N measured in the soil profile from each experimental plot. The RSM approach in this study uses a surface regression least squares method to fit linear, quadratic, and cross product response combined surfaces. 5 primary steps were performed for the analysis and were taken from multiple sources in the literature (eg, Box and Draper,²⁸ Khuri and Cornell,²⁹ SAS Institute Inc.³⁰):

Create a simple linear model and explore the nature of the response, β coefficients (ie, the slope of each explanatory variable), magnitude, and direction and level of significance.

Continue the exploration by moving towards the optimum zone response surface by following the favorable optimization (ie, the steepest descent or ascent) direction. For example, add more data, investigate relationships between the response and explanatory variables, and add scientifically sound (ie, statistically significant) variables.

If there are a large number of “hills” and “valleys” in the response surface graph, then a first-order model is inappropriate and a more complex model (ie, a second- or third-order model) will be needed to obtain a meaningful RSM approximation.

Perform a (lack of) fitness test using the F statistical test to measure overall performance of the model.

Use the final response surface model to perform the analysis and then summarize implications of the results.

For this study, 4 optimization models were constructed. 2 targets were maximized (yield and gross margin) on 2 crops (corn and soybeans), subject to limits of observed NO₃-N in tile drainage flow and soil profile N. The overall model problem was estimated using the RSREG regression procedure in the Statistical Analysis System (SAS Institute²⁸) software program, which uses the method of least squares to fit quadratic response surface regression models. For this study, the RSREG procedure was formulated as follows:

Maximize RSREG = Y_ij subject to Y_ik with limits L_L to L_U

where:

RSREG = Calculated response surface regression model.

L_L = Lower empirical level of NO₃-N in tile drainage flow or soil profile N.

L_U = Higher empirical level of NO₃-N in tile drainage flow or soil profile N.

The general format of the production or gross margin function is:

Y_{ij} = β_{0} + β_{1} x_{1} + β_{2} x_{2} + β_{3} x_{1}^{2} + β_{4} x_{2}^{2} + β_{5} x_{1} * x_{2} + ε

(1)

where ∊ is the error term and Y_ij are the dependent variables [i is the crop (corn or soybean) and j is the yield or gross margin] subject to Y_ik = NO₃-N in tile drainage or soil profile N (for corn or soybean). The model explanatory variables for corn were:

x₁ = crop rotation (classification variable for rotation: continuous corn, corn-soybean, or soybean-corn).

x₂ = tillage system (classification variable for tillage system type: chisel plow, reduced-till, moldboard plow, or no-till).

x₃ = planting rate (# seeds/m of row).

x₄ = total profile soil water content (calculated using the RZWQM2 model for each experimental plot during the 1990-2003 experimental analysis period).

x₅ = total applied N content (by plot) calculated as: N content = weight %^*quantity applied or application rate.

where: weight% = atomic weight of N/molecular weight of the fertilizer)^*100].

x₆ = year (classification variable for the years during the 1990-2003 experimental analysis period).

x₇ = effective soil porosity, ie, fraction of total volume where water flow effectively occurs (calculated using RZWQM2 for each plot during the 1990-2003 experimental analysis period).

Profile soil water content and effective soil porosity values were taken from simulations of the Nashua, IA experimental data plots in which the RZWQM2 model with 26 years of data (1978-2003) was used to evaluate year-to-year crop yield, water, and N balances.³¹ RZWQM2 was calibrated using data from 2 0.4 ha plots and evaluated by comparing simulated values for crop yield, soil water content, tile drainage flow, NO₃-N loading in tile drainage flow, and soil profile N with the corresponding measured data from the experimental plots. Root mean square errors (RMSE) for simulated soil water storage, water table, and annual tile flow were 3.0, 22.1, and 5.6 cm, respectively. These values were close to the average RMSE for the measured data between replicates (3.2, 22.4, and 5.7 cm, respectively). RMSE values for simulated annual N loading and soil profile N were 16.8 and 47.0 kg N ha^-1, respectively, which were somewhat higher than the average RMSE for measurements among replicates (12.8 and 38.8 kg N ha^-1, respectively). The model explanatory variables for soybeans were:

x₁ = crop rotation (classification variable for rotations: corn-soybean or soybean-corn).

x₂ = tillage system (classification variable for tillage system type: chisel plow, reduced-till, moldboard plow, or no-till).

x₃ = total profile soil water content (calculated using RZWQM2 for each plot during the 1990-2003 experimental analysis period).

x₄ = year (classification variable for the years during the 1990-2003 experimental analysis period).

x₅ = effective soil porosity, ie, fraction of total volume where water flow effectively occurs (calculated using RZWQM2 for each plot during the 1990-2003 experimental analysis period).

The above models were optimized using the RSREG command in SAS, which automates most of these steps, ie, a model is fitted to a surface for Y₁, subject to the limits on Y₂ and Y₃, while also solving for the independent variables of interest. The RSREG procedure uses the lack-of-fit concept to confirm the parameters of interest, a canonical analysis to optimize the shape of the function (accounting for local optima), and then ridge analysis to search for regional optima.³⁰ In general, a lack-of-fit test determines if the response surface model needs to be improved through the addition or subtraction of independent or explanatory variables. The test compares the variation around the model with “pure” variation within replicated observations. This measures the adequacy of the response surface model. In particular, if there are n replicated observations Y_il, …, Y_in of the response all at the same values x_k of the factors, then the true response can be predicted at x_i either by using the predicted Ŷ value based on the model or by using the mean Ῡ of the replicated values. Thus, the test for lack-of-fit²⁸ decomposes the residual error into a component due to the variation of the replications around their mean value (the “pure” error), and a component due to the variation of the mean values around the model prediction (ie, the “bias” error):

\sum_{i = 1}^{n} \sum_{j = 1}^{k} {(Y_{ij} - {\hat{Y}}_{i})}^{2} = \sum_{i = 1}^{n} \sum_{j = 1}^{k} (Y_{ij} - {\bar{Y}}_{i}) + \sum_{i} k({\bar{Y}}_{i} - {\hat{Y}}_{i})^{2}

(2)

If the model is adequate, then both components estimate the nominal level of error; however, if the bias component of error is much larger than the pure error, this constitutes evidence that there is significant lack-of-fit.

Final Response Surface Estimation

The response surface model steepest ascent or steepest descent process combines canonical and ridge analysis to develop a comprehensive surface over the range of values of interest for the 3 dependent variables. It is an iterative process. The search for an optimum response surface begins with simple linear models and simple designs to explore the nature of the response function in the vicinity/region of where we think the final factor combination settings would be located.²⁸ The search for the optimum response starts with testing whether a first-order linear model provides adequate fitness (using the lack-of-fitness process explained above) for the response surface at the initial response region. The next step is usually to select local optimum points along the direction of steepest ascent (or descent, if response is to be decreased) and continue the search as long as the response is behaving as expected, ie, the process is continuing to move in the current direction. When the response begins to decrease rather than increase (or vice versa), then new factor points are identified (using the SAS RIDGE procedure) and a new direction of steepest ascent or descent is determined (provided that a first-order model is satisfactory). The RIDGE procedure computes the ridge of optimum response starting at a given point x_o; the point on the ridge at radius r from x_o is the collection of factor settings that optimizes the predicted response at this radius. The ridge is analogous to climbing or falling as fast as possible on the surface of predicted response. Thus, the ridge analysis is useful as a tool to help interpret an existing response surface or to indicate the direction in which further experimentation should be performed.

If a first-order model is inadequate, then the search for the optimum points is enlarged by conducting additional runs at appropriately selected experimental points such that a second-order model may be fitted and its coefficients estimated.³² The process continues to a third-order model (ie, cross product of the explanatory variables) until the best fit model is obtained. The RSM process is interactive in the sense that the response surface modeler should test for first-order model fitness, add (scientifically sound) or eliminate (scientifically unsound) explanatory variables, modify input data (eg, biophysical model data), and re-test for the fitness again. A modeler would not typically move to a second-order model unless the first-order model showed a lack-of-fitness (and so forth to the third-order model). This process is applicable only in a situation such as this study where there is lack of a predefined functional form to reflect the complexity of the research problem. That is, the complexity in this study is primarily due to the lack of a predefined single functional form for biophysical variables, environmental variables, and economic variables that can be combined together in a single model.

Example of using Rsm with Agro-Ecosystem Modeling

To illustrate the above discussion of RSM, the following is a simple example of how the approach could be used in conjunction with an agroecosystem simulation model (such as RZWQM2) to approximate the process of identifying new plant cultivar technology as a result of potential change in irrigation management. Assume that a scientist would like to find levels of 2 inputs, D (a cultivar genetic trait such as days to maturity) and P (a management variable such as planting rate), which maximize yield (Y) of a particular crop for a different irrigation management scenarios. Furthermore, let yield be only a function of inputs D and P and a vector of irrigation management-related variables (eg, irrigation rate and timing) I or more specifically, Y = f(D,P,I) where a simulation model such as RZWQM2 accurately quantifies this functional relationship. A predicted yield value is then generated for a particular set of inputs D and P by simulating for a particular irrigation management scenario I and using D, P and I in the RZWQM2 model to generate crop yields over time. In a broad sense, RSM proceeds by successively adjusting D and P until maximum yield is achieved which approximates the scientist's search for improved cultivars and irrigation management practices.

In this simple example, the RSM procedure begins with an initial factorial first-order experiment, centered at the “current” levels of D and P (point a in Fig. 1) where each of the 4 design points are identified by the levels of D and P in Figure 1 and yield is measured on a third axis that is perpendicular to the D, P plane. Point a represents the current cultivar and irrigation management practice. At each of the four design points about a, yield values are simulated over a number of years and the mean yield (Y) is obtained for each point. Based on these 4 Y values, a simple first-order statistical model: Y = b0 + b₁X₁ + b₂X₂ is fit to give a planar yield response surface as a function of D and P. Using this surface, the path of steepest ascent is determined as the line that predicts the steepest increase in yield (line ab in Fig. 1). Yield values are simulated sequentially at various D and P values along the path of steepest ascent until yield decreases substantially. Another factorial experiment is then conducted near the point of highest yield on the path of steepest ascent (point c, Fig. 1), another first-order statistical model is fit and a second path of steepest ascent is identified (line cd, Fig. 1). The process is continued until there is little increase in yield at which point a final second-order experiment is conducted to identify the values of D and P that maximize yield (D_opt.|P_opt. in Fig. 1). This final experiment is normally a central composite design and data from this experiment are analyzed with a second-order model y = b0 + b₁X₁ + b₂X₂ + b₁₁X₁² + b₂₂X₂² + b₁₂X₁X₂. Differentiating the estimated equation with respect to D, setting the result equal to zero and solving gives the optimal inputs (D_opt.|P_opt.). Evaluation of the fitted response surface then determines the nature of the surface and the nature of the D_opt.|P_opt. values, ie, if they are maximum, minimum, or saddle-points. Ideally, the final values of D_opt. (days to maturity) and P_opt. (planting rate) will maximize yield under future irrigation management scenarios at this location. Even though the values of D_opt. and P_opt. are based on agro-ecosystem model simulations, they can be useful in understanding the types of cultivars and irrigation management practices that may be efficacious in the future. Although this example is somewhat unsophisticated, the method is quite general since any cropping system model with any number of input variables and any climate scenario can be used as long as the output (Y), the inputs (D, P) and the irrigation management variables (C) are clearly identified.

Figure 1.

Example of RSM (path of steepest ascent) to identify maximum yield and optimal value for two explanatory variables based on hypothetical agroecosystem model simulation experiments.

Results and Discussion

Table 2 summarizes the results for the 4 response surface models (eg, corn yield, corn gross margin, soybean yield, and soybean gross margin) used in this study for the first local optima reached in the RSM (ie, the point closest to the origin). This solution was found to be a local, but not necessarily global optimum using the steepest ascent method; ridge analysis produced other local optimum through additional searching of the factor space. In Table 2, the local optimum yield for corn is 7.6 t/ha, which would generate an optimum gross margin of $308/ha. The environmental externality of the optimum corn yield and gross margin values is 20 kg/ha tile drainage NO₃-N and 54 kg/ha total soil profile N. The optimum yield for soybean is 2.9 t/ha which would generate an optimum gross margin of $380/ha. The environmental externality of the optimum soybean yield and gross margin values is 14 kg/ha tile drainage NO₃-N and 41 kg/ha total soil profile N (Table 2). In this case, the economic versus N environmental impact trade-off clearly varies by crop grown.

Table 2.

Optimum values for the dependent variables (corn/soybean yield and gross margin) versus tile drainage NO₃-N and soil profile N.^#

	Yield (t/ha)	Gross margin ($/ha)	Tile drainage NO₃-N (kg/ha)	Soil profile N (kg/ha)
Corn
Mean	7.6	308	20	55
RMSD	0.7	69	11	28
R²	0.8	0.73	0.81	0.64
CV	9.2	22	53	51
Soybean
Mean	2.9	380	14	41
RMSD	0.3	85	8	23
R²	0.7	0.75	0.70	0.67
CV	9.6	22	56	56

Abbreviations: ^*RMSD, root mean square deviation; R², coefficient of determination; CV, coefficient of variation.

Table 3 summarizes the statistical analysis of the 4 response surface models as they are expanded from the simple linear model to the full quadratic model with cross products. The R² of the 4 models ranges between 0.73 and 0.84, with the corn yield model having the highest R² and the corn gross margin and soybean yield models having the lowest. The linear model in all 4 cases contributed the most to overall model significance between the linear, quadratic, and cross product components of the response surface model. However, the response surface models were all improved when made more complex, ie, through the addition of soil water content (data not shown). Table 4 presents the lack-of-fitness measures for the 4 final response surface models. The results indicate the absence of lack-of-fitness (ie, the models are statistically sound based on the F and P statistics for all of the response surface models with the soybean gross margin model performing the best (ie, it has the highest F and lowest P values). This result indicates that using multiple response method regression analysis demands detailed data on both the dependent variable and the explanatory variables in order to accurately simulate response surface models. The differences in the results between yield and gross margin in Table 4 also indicate the influence of input and output price in the trade-off analysis. Farmers face variation in their crop and input prices at the marketplace. Therefore, market prices play an important role in the farmer's decision to accept the tradeoff between using higher levels of N fertilization and environmental conservation efforts. Unfortunately, a full investigation of price variability is beyond the scope of this investigation.

Table 3.

Response surface statistical analysis.

Regression model	Degrees of freedom (DF)	Type I sum of squares	Coefficient of determination (R²)	F value
Corn yield
Covariates (treatments)	1	0.01	0.00	0.02
Linear	7	513.70	0.65	148.25
Quadratic	7	58.48	0.07	16.88
Cross product	21	92.49	0.12	8.90
Total model	36	664.67	0.84	37.30
Corn gross margin
Covariates (treatments)	1	64,556	0.01	13.49
Linear	7	1,451,045	0.32	43.32
Quadratic	7	653,825	0.14	19.52
Cross product	21	1,206,973	0.26	12.01
Total model	36	3,376,399	0.73	19.60
Soybean yield
Covariates (treatments)	1	0.18	0.00	2.27
Linear	5	24.44	0.44	61.26
Quadratic	4	9.42	0.17	29.52
Cross product	10	6.81	0.12	8.54
Total model	20	40.86	0.73	25.60
Soybean gross margin
Covariates (treatments)	1	109,860	0.02	15.07
Linear	5	1,613,175	0.29	44.27
Quadratic	4	1,213,498	0.22	41.62
Cross product	10	1,269,317	0.23	17.42
Total model	20	4,205,851	0.75	28.85

Table 4.

Response surface fitness analysis.

Response surface model	Degrees of freedom (DF)	Sum of squares	Mean square	F value	P value
Corn yield	253	126	0.50	2.25	0.36
Corn gross margin	253	1,217,254	4811.28	3.18	0.27
Soybean yield	187	15	0.08	2.77	0.16
Soybean gross margin	187	1,386,897	7416.56	5.74	0.05

Table 5 shows that all selected explanatory variables were significant at the P < 0.10 level with the exception of effective porosity. However, the inclusion of effective porosity in the response surface models plays a critical role in changing the lack-of-fitness measure from positive to negative (data not shown). In other words, the success of including an explanatory soil property variable (ie, effective porosity) in the response surface models indicates the importance of using biophysical measures when conducting a trade-off analysis. Therefore, this study shows both detailed economic and biophysical variables may be required to robustly analyze trade-offs.

At this point, various economic and environmental trade-offs can be explored. The complete surfaces for yield and gross margins are shown for corn and soybeans in Figures 2–5. As would be expected, the surfaces are complex and exemplify many potential trade-offs. It is difficult to discern any clear patterns; however, local optima for yield or gross margin can be easily observed where the surface points upward. Fortunately, ridge analysis revealed that there were only 4 local optima for both yield and for gross margin, which greatly reduces the dimensionality of the problem without any loss of relevant information. As presented in Table 6, local optima for corn yield and gross margin ranges from 10.4 t/ha to 7.6 t/ha and $485/ha to $308/ha, respectively. Local optima for soybean yield and gross margin ranges from 4.5 t/ha to 2.9 t/ha and $546/ha to $380/ha, respectively. The level of environmental impact is also provided for each local optimum listed in Table 6, which allows for easy comparison of the trade-offs between the economic (ie, yield and gross margin) and environmental (ie, tile drainage NO₃-N and total soil profile N) objectives. As an example for soybeans, tile drainage NO₃-N could be reduced from 53 kg/ha to 40 kg/ha and total soil profile N from 218 kg/ha to 62 kg/ha if producers were willing to accept a reduction in gross margin from $546/ha to $421/ha. For the trade-off analysis, all inefficient choices were removed from consideration. The remaining local optima (Table 6) require a choice between economic gain or loss and corresponding environmental impacts. Furthermore, in each case for corn or soybean yields/gross margins, the local optima show a reduction in returns required to achieve a marginal reduction in nitrogen.

Table 6.

Trade-offs for the dependent variables (corn/soybean yield and gross margin) and tile drainage NO₃-N and soil profile N.

Yield vs. N trade-off	Yield (t/ha)	Yield reduction (%)	Tile drainage NO₃-N (kg/ha)	Tile drainage NO₃-N reduction (%)	Soil profile N (kg/ha)	Soil profile N reduction (%)
Corn	10.4^#		55		114
	8.7	16	49	11	93	18
	7.8	25	44	20	72	37
	7.6	27	20	64	55	52
Soybean	4.5^#		40		160
	3.8	16	34	15	135	16
	3.2	29	13	41	85	47
	2.9	36	14	59	41	74
Gross margin vs. N trade-off	($/ha)
Corn	485^*		77		156
	451	7	55	29	93	40
	376	23	44	43	62	60
	308	37	20	74	55	65
Soybean	546^*		53		218
	450	18	47	11	181	17
	421	23	40	25	62	72
	380	30	14	74	41	81

Maximum possible corn/soybean yield within the RSM constrained region (only local optima are shown).

Maximum possible corn/soybean gross margin within the RSM constrained region (only local optima are shown).

Figure 2.

3-D trade-offs between corn yield, tile drainage N, and soil profile N.

Figure 3.

3-D trade-offs between corn gross margin, tile drainage N, and soil profile N.

Figure 4.

3-D trade-offs between soybean yield, tile drainage N, and soil profile N.

Figure 5.

3-D trade-offs between soybean gross margin, tile drainage N, and soil profile N.

Table 6 also shows that the soybean trade-off results are much different than the corn trade-off results. The optimum response surface soybean yield is 2.9 t/ha (Table 2). At this level, the expected tile drainage NO₃-N is 14 kg/ha while the total soil profile N is 41 kg/ha (Table 6). Increasing soybean productivity to 4.5 t/ha would add 26 kg/ha NO₃-N to tile drainage and 119 kg/ha N to the soil profile. This indicates that very high levels of productivity in soybean production can cause proportionally significantly high levels of N in the soil profile and therefore potentially higher levels of environmental damage (eg, excessive N leaching to groundwater) as compared to corn. Table 6 also shows the percent reduction in yield, tile drainage NO₃-N, and soil profile N for each of the optima compared to the maximum possible corn/soybean yield and gross margin within the RSM optimum constrained region.

It is not easy to fully visualize the 3-D economic and environmental trade-offs in Figures 2–5. Therefore, the information in Table 6 was used to create visual 2-D representations of the economic and environmental trade-off frontiers. The 2-D trade-offs for corn yield and corn gross margin are presented in Figures 6–9, showing the 4 local optima for each trade-off case. For example, a farm manager must give up a substantial amount of corn yield (10.4 t/ha -8.7 t/ha = 1.7 t/ha) to reduce total soil profile N from 114 to 93 kg/ha (or 21 kg/ha annually) while reducing tile drainage NO₃-N by 6 kg/ha. Note that the optimum response surface level of corn yield is 7.6 t/ha (Fig. 6). At this level, both tile drainage NO₃-N and total soil profile N reach their minimum possible levels (20 kg/ha and 55 kg/ha, respectively, as shown in Table 6). As expected, when prices and production costs enter the equation in the form of gross margin as the trade-off factor, the levels of N reductions in the economic case are not proportional to the yield reduction case to the fact that yield and gross margin have a nonlinear relationship. As an example for corn, the farm manager must sacrifice approximately $177/ha (from $485/ha to $308/ha) in gross margin in order to reduce tile drainage NO₃-N by 57 kg/ha and total soil profile N by 101 kg/ha annually, or 65% (Table 6 and Fig. 7). This result is consistent with the study of Koikkalainen et al³³ that showed a trade-off of about € 52/ha when farmers reduced N use by over 50%.

Figure 6.

2-D trade-off frontier between corn yield, tile drainage N, and soil profile N.

Figure 7.

2-D trade-off frontier between corn gross margin, tile drainage N, and soil profile N.

Figure 8.

2-D trade-off frontier between soybean yield, tile drainage N, and soil profile N.

Figure 9.

2-D trade-off frontier between soybean gross margin, tile drainage N, and soil profile N.

Conclusions and Limitations

Our results indicate that trade-offs between farm profitability and environmental externalities are complex and vary significantly depending on the crop. In this study, we present detailed tabular and graphical economic and environmental trade-offs between corn and soybean yield/gross margin and N in tile drainage flow and the soil profile. An important finding of this research is that there was no one single point of trade-off between economics and the environment. Furthermore, the trade-offs reflect economic and environmental returns to scale. That is, if a farm manager reduces the amount of N applied by 1 half, then the environmental benefits would more than double. However, if the same farm manager doubled the amount of N applied, the economic benefits would be more than double. We believe that this concept could have direct and significant policy consequences. For instance, a policy that targets support of farmers who are willing to trade-off economic gains for environmental benefits should consider both economic returns to scale and the (multiplicative) environmental benefits due to a reduction in resource use (the implications of which are shown herein).

Response surface methodology (RSM) provides a useful mechanism to quantitatively evaluate tillage and cropping system treatments under northeast Iowa climatic conditions. However, the underlying regression equations are limited in that extrapolation beyond the experimental data may be unwarranted, ie, applying the RSM regression equations beyond the experimental conditions would require comprehensive re-testing and modification of the regression models for those conditions. Other important weather variability that drives NO₃-N loading such as rainfall distribution within the season was not considered because additional variables would increase the complexity and the limited years of observations poses difficulty. The use of process-based models was recommended by Malone et al⁹ to overcome some of the limitations of regression-based modeling and to more comprehensively evaluate and quantify yield and NO₃-N loading. The reasoning behind their recommendation was: (1) process-based models allow extrapolation of management and climate effects to conditions (climate, soil, management) where observed data is sparse or nonexistent, and (2) process-based models also allow cause and effect analysis because observed data is necessarily limited.⁹ However, direct use of process-based models to predict NO₃-N transport in artificially drained soil is too time-consuming for conservation planners and land managers. Prior research has suggested using a meta-modeling (eg, polynomial regression, splines, and neural networks) approach to upscale field scale modeling results of N leaching to regions (eg, Wu and Babcock,³⁴ Borgesen et al,³⁵ Haberlandt et al³⁶). Development of the response surfaces for the Nashua experimental plot data suggests meta-models can be developed to quantify NO₃-N leaching and crop yield under a variety of climate and management conditions in artificially drained soil.

In summary, farm managers face challenging decisions when considering trade-offs between economics and the environment. Trade-off analysis is a complex issue and involves a large numbers of factors that influence the decision maker. The use of RSM was found to be appropriate in this study for analyzing this type of complex problem in that quantifiable independent variables were identified (eg, crop rotation, tillage system, planting rate, profile soil water content, amount of N applied) that significantly affected the dependent variables corn/ soybean yield and gross margin. These variables were predicted with reasonable accuracy; therefore, the developed RSM regression equations are a step toward development of a simple, accurate, and objective method to quantify management and climate effects on NO₃-N loading and crop yield for a region.

Funding

Author(s) disclose no funding sources.

Author Contributions

Conceived and designed the experiments: JCA, EMF, DLH. Analyzed the data: EMF, JCA. Wrote the first draft of the manuscript: EMF, JCA. Contributed to the writing of the manuscript: JCA, EMF, DLH. Agree with manuscript results and conclusions: JCA, EMF, DLH. Jointly developed the structure and arguments for the paper: JCA, EMF, DLH. Made critical revisions and approved final version: JCA, EMF, DLH. All authors reviewed and approved of the final manuscript.

Competing Interests

Author(s) disclose no potential conflicts of interest.

Disclosures and Ethics

As a requirement of publication the authors have provided signed confirmation of their compliance with ethical and legal obligations including but not limited to compliance with ICMJE authorship and competing interests guidelines, that the article is neither under consideration for publication nor published elsewhere, of their compliance with legal and ethical guidelines concerning human and animal research participants (if applicable), and that permission has been obtained for reproduction of any copyrighted material. This article was subject to blind, independent, expert peer review. The reviewers reported no competing interests.

Footnotes

Acknowledgements

We acknowledge the Department of Agricultural and Biosystems Engineering at Iowa State University for initiating and maintaining the experimental study at the Nashua Northeast Research Station, collecting the data, and graciously sharing the data with us. In addition, we acknowledge Kevin Boyle with the USDA NRCS for his role in initially summarizing the Nashua economic data under the NRCS EconDocs exchange tool.

References

D'Arcy

Frost

The role of best management practices in alleviating water quality problems associated with diffuse pollution.

Science of the Total Environment. 2001; 265(1-3): 359–67.

Shaffer

M.J.

Halvorson

A.D.

Pierce

F.J.

Nitrate leaching and economic analysis package (NLEAP): Model description and application. In: Folett

R.F.

Keeney

D.R.

Cruse

R.M.

, editors. Managing Nitrogen for Groundwater Quality and Farm Profitability. 1991; Madison: SSSA; 1991: 285–322.

Ahuja

L.R.

Rojas

K.W.

Hanson

J.D.

Shaffer

M.J.

, editors. Root Zone Water Quality Model. 2000; Englewood: Water Resources Publications.

Simmelsgaard

S.E.

Djurhuus

An empirical model for estimating nitrate leaching as affected by crop type and the long-term N fertilizer rate.

Soil Use and Management. 1998; 14: 37–43.

Hansen

Kristensen

E.S.

Grant

Høgh-Jensen

Simmelsgaard

S.E.

Olesen

J.E.

Nitrogen leaching from conventional versus organic farming systems—a systems modelling approach.

European Journal of Agronomy. 2000; 13: 65–82.

Haberle

Kusá

Svoboda

Klír

The changes of soil mineral nitrogen observed on farms between autumn and spring and modelled with a simple leaching equation.

Soil and Water Research. 2009; 4: 159–67.

Pedersen

J.K.

Bjerg

P.L.

Christensen

T.H.

Correlation of nitrate profiles with groundwater and sediment characteristics in a shallow sandy aquifer.

Journal of Hydrology. 1991; 124: 263–77.

Shepherd

Harper

Millington

Modelling catchment-scale nutrient transport to watercourses in the UK.

Hydrobiologia. 1999; 396: 227–37.

Malone

R.W.

Karlen

D.L.

Empirical analysis and prediction of nitrate loading and crop yield for corn-soybean rotations.

Geoderma. 2007; 140: 223–34.

10.

Tóth

S.F.

Rabotyagov

S.S.

Ettl

G.J.

Experimental testbeds for ECOSEL: A market framework for private provision of forest ecosystem services.

Agricultural and Applied Economics Association. 2009; 49565.

11.

Box

G.E.P.

Draper

N.R.

Empirical Model-Building and Response Surfaces. Wiley Series in Probability and Mathematical Statistics. 1987; New York: Wiley & Sons, Inc.

12.

Downing

D.J.

Gardner

R.H.

Hoffman

F.O.

An examination of response-surface methodologies for uncertainty analysis in assessment models.

Technometrics. 1985; 27(2): 151–63.

13.

Kelly

T.C.

Teasdale

Economic-environmental trade-offs among alternative crop rotations.

Agriculture, Ecosystems, and Environment. 1996; 60: 17–28.

14.

Burdick

D.S.

Naylor

T.H.

Response surface methods in economics.

Review of the International Statistical Institute. 1969; 37(1): 18–35.

15.

Fitzgerald

G.J.

Lesch

S.M.

Barnes

E.M.

Luckett

W.E.

Directed sampling using remote sensing with a response surface sampling design for site-specific agriculture.

Computers and Electronics in Agriculture. 2006; 53: 98–112.

16.

Lesch

S.M.

Sensor-directed response surface sampling designs for characterizing spatial variation in soil properties.

Computers and Electronics in Agriculture. 2005; 46: 153–79.

17.

Frey

H.C.

Patil

S.R.

Identification and review of sensitivity analysis methods.

Risk Analysis. 2002; 22: 553–78.

18.

Dhungana

Eskridge

K.M.

Weiss

Baenziger

P.S.

Designing crop technology for a future climate: An example using response surface methodology and the CERES-Wheat model.

Agricultural Systems. 2006; 87: 63–79.

19.

Isukapalli

S.S.

Roy

Georgopoulos

P.G.

Stochastic response surface methods (SRSMs) for uncertainty propagation: Application to environmental and biological systems.

Risk Analysis. 1998; 18(3): 351–63.

20.

Karlen

D.L.

Berry

E.C.

Colvin

T.S.

Kanwar

R.S.

Twelve year tillage and crop rotation effects on yields and soil chemical properties in northeast Iowa.

Communications in Soil Science and Plant Analysis. 1991; 22: 1985–2003.

21.

Bakhsh

Kanwar

R.S.

Karlen

D.L.

Tillage and nitrogen management effects on crop yield and residual soil nitrate.

Transactions of the ASAE. 2000; 44: 1589–95.

22.

Voy

K.D.

Soil Survey of Floyd County, Iowa. USDA-SCS in cooperation with Iowa Agriculture and Home Economics Experimental Station Cooperative Extension Service. 1995; Ames: Iowa State University.

23.

Williams

J.R.

A stochastic dominance analysis of tillage and crop insurance practices in a semiarid region.

American Journal of Agricultural Economics. 1988; 70: 112–20.

24.

Yiridoe

E.K.

Weerksink

Hooker

D.C.

Vyn

T.J.

Swanton

Income risk analysis of alternative tillage systems for corn and soybean production on clay soils.

Canadian Journal of Agricultural Economics. 2000; 48: 161–74.

25.

Chase

C.A.

Duffy

M.D.

An economic analysis of the Nashua tillage study: 1978-1987.

Journal of Production Agriculture. 1991; 4: 91–8.

26.

Kay

R.D.

Farm Management: Planning, Control, and Implementation. 2nd ed. 1986; New York: McGraw-Hill Book Company.

27.

Box

G.E.P.

Wilson

K.B.

On the experimental attainment of optimum conditions.

Journal of the Royal Statistical Society. Series B (Methodological). 1951; 13: 1–45.

28.

Box

G.E.P.

Draper

N.R.

Response Surfaces, Mixture and Ridge Analysis. 2007; New York: Wiley-Interscience.

29.

Khuri

A.I.

Cornell

J.A.

Response Surface: Designs and Analyses, Second Edition. 1996; New York: CRC Press.

30.

SAS Institute Inc. SAS 9.2 Help and Documentation Manual. 2008; Cary: SAS Institute.

31.

Malone

R.W.

Heilman

RZWQM simulation of long-term crop production, water and nitrogen balances in Northeast Iowa.

Geoderma. 2007; 140: 247–59.

32.

Littell

R.C.

Milliken

G.A.

Stroup

W.W.

Wolfinger

R.D.

Schabenberber

SAS for Mixed Models. 2nd ed. 2006; SAS Institute Publishing.

33.

Helin

Laukkanen

Koikkalainen

Abatement cost for agricultural nitrogen and phosphorus loads: a case study of southwestern Finland.

MTT Discussion Papers. 2006: 8.

34.

Babcock

B.A.

Metamodeling potential nitrate water pollution in the central United States.

Journal of Environmental Quality. 1999; 28: 1916–28.

35.

Borgesen

C.D.

Djurhuus

Kyllingsbaek

Estimating the effect of legislation on nitrogen leaching by upscaling field simulations.

Ecological Modelling. 2001; 136: 31–48.

36.

Haberlandt

Krysanova

Bardossy

Assessment of nitrogen leaching from arable land in large river basins Part II: regionalisation using fuzzy rule based modeling.

Ecological Modelling. 2002; 150: 277–94.