Calibration of a crop model to irrigated water use using a genetic algorithm

Near-term consumption of groundwater for irrigated agriculture in the High Plains Aquifer supports a dynamic bio-socio-economic system, all parts of which will be impacted by a future transition to sustainable usage that matches natural recharge rates. Plants are the foundation of this system and so generic plant models suitable for coupling to representations of other component processes (hydrologic, economic, etc.) are key elements of needed stakeholder decision support systems. This study explores utilization of the Environmental Policy Integrated Climate (EPIC) model to serve in this role. Calibration required many facilities of a fully deployed decision support system: geo-referenced databases of crop (corn, sorghum, alfalfa, and soybean), soil, weather, and water-use data (4931 well-years), interfacing heterogeneous software components, and massively parallel processing (3.8 ×109 model runs). Bootstrap probability distributions for ten model parameters were obtained for each crop by entropy maximization via the genetic algorithm. The relative errors in yield and water estimates based on the parameters are analyzed by crop, the level of aggregation (countyor well-level), and the degree of independence between the data set used for estimation and the data being predicted.


Introduction
Regionally, short-term consumption of groundwater in the High Plains Aquifer provides for a dynamic bio-socioeconomic system through irrigated agriculture.In the long Correspondence to: S. M. Welch (welchsm@ksu.edu)term, transition to sustainable usage that matches natural recharge rates will impact ecologies, economies, demographics and the landscape.Recharge is that portion of precipitation not lost as evaporation from foliage, run off (affected by ground cover), or root uptake.This problem is of global significance as National Geographic (Montaigne, 2002) declared the High Plains Aquifer to be one of 22 worldwide "critical areas" for "annual renewable water".
This problem has been studied extensively from various disciplinary perspectives with disparate unaligned concepts, viewpoints, vocabulary, models and data.Stakeholder decision makers in this system are equally distributed across a mix of governmental agencies, administrative units, private sector enterprises, and farmers.Disjoint disciplinary science leaves these decision makers ill equipped to understand how consequences of management actions and policies impact and cascade through the integrated system.Thus, all stakeholders share a common need for integrated, science-based, quantitative informational tools that, collectively, target their differing individual responsibilities.Toward this end, researchers at Kansas State University, in conjunction with stakeholder groups, have begun to integrate economic, agronomic, and hydrologic models, supported by geodatabases, to aid in these diverse decision making processes (Steward et al., 2005(Steward et al., , 2009a;;Steward and Bernard, 2006a,b;Bernard et al., 2004Bernard et al., , 2005;;Yang et al., 2009).
Plants (Fig. 1) form the foundation of the human-natural system in the High Plains.Irrigation to meet transpiration needs comprises over 95% of groundwater use in portions of the High Plains Aquifer (Wilson et al., 2000).Statistical crop yield estimators used as economic production functions or in data summarization (Berck and Helfand, 1990;Frank et al., 1990;Paris, 1992) often do not explicitly represent physical water fluxes and may not partition landscapes in ways directly related to hydrological features or patterns of diversion.In contrast, physiological, parcel-based crop simulators, including the Environmental Policy Integrated Climate (EPIC) model 1 (Sharpley and Williams, 1990;Williams et al., 1990), are well suited for linkage to other models and to Geographic Information Systems (Lal et al., 1993;Engel et al., 1997;Yang et al., 2003).Our specific objective is to estimate irrigation needs for large numbers of representative, geo-referenced parcels.This flux couples directly to hydrological models.We have used Sheridan County, Kansas, as a study area to prototype this estimation process.This study evaluates the suitability of the EPIC model for providing crop simulation to decision support systems at the regional scale.In this work we chose a county size as this represents a standard land unit size for aggregation of information reported about crop production (yields, etc.), information that was needed for this study.We are also evaluating extending these methods to larger scales such as the Ogallala Aquifer portion of Kansas.
The first step is to calibrate the EPIC model.Although basically a parameter estimation process, heavy computational requirements mandate development of much of the infrastructure required by a fully deployed decision support system.This includes geo-referenced databases of crop, soil, weather, and water-use data, interfacing heterogeneous software components (i.e.model, optimizer, data retrieval), and, most importantly, distributed parallel processing.The 1 Originally named the Erosion Productivity Impact Calculator.
latter is required because benchmark runs indicated that calibration would require 4 to 14 d on a 40-CPU computing cluster (a significant underestimate as events later proved).On a computations-per-minute basis, this appeared representative of the computing intensity required for multi-year, spatially distributed, water policy analyses.The following sections present the elements of our calibration approach and the results obtained.

The EPIC model
Plant processes have been extensively modeled (Bowen, 1992;Bouman et al., 1996).The EPIC model simulates the physiology of all major forages and crops in the study area.Using a daily time step, three major processes are represented: (1) phenological development; (2) dry matter production and partitioning to plant tissues, resulting in growth; and (3) economic yield.Outputs that are relevant to this study are crop yield and water use reported in t/ha and mm-ha, respectively.The model reproduces the results of irrigation, fertilization, tillage, variety selection, alternative production calendars, etc. EPIC also includes an economic component for evaluating and optimizing management outcomes.(In our research, however, a more robust economic forecasting submodel is being used (Peterson and Steward, 2006) to suggest crop management choices for EPIC to simulate.)Because its original focus was erosion-related, EPIC can simulate decade-scale or longer intervals.These features have suited EPIC to a broad range of applications, including plant nutrition studies (Cole et al., 1987;Dautrebande et al., 1999); national and international assessments of agroecological change impact (Brown and Rosenberg, 1999;Brown et al., 2000;Bernardos et al., 2001), including the High Plains Aquifer (Easterling III et al., 1993); irrigation planning and water use (Bryant et al., 1992;Ellis et al., 1993;Evers et al., 1998;Guerra et al., 2005); and regional studies (Geleta et al., 1994;Cabelguenne et al., 1995;Fortin and Moon, 2000).
EPIC can be divided into nine subroutines of which hydrology and plant growth were of interest for our simulations (Williams et al., 1990).The hydrology subroutine is composed of surface runoff, percolation, lateral sub-surface movement and evaporation.These processes in our simulations were controlled by the parameters used to describe the soil groups.Slope, and NRCS runoff curves, soil water content and rainfall amounts determine runoff.Percolation and lateral sub-surface movement is controlled by the soil layer data.Potential evaporation was estimated using the Penman-Monteith method.
Plant growth is determined on a daily time step based on intercepted solar radiation.Daily plant growth is estimated as a function of intercepted solar energy and plant leaf area.Daily dry matter is accumulated for the growing season that is controlled by heat units or environmental conditions (typically freeze events for summer crops) and yield is estimated using a total biomass to grain ratio, which is referred to as a harvest index.Species-specific parameters distinguish between crops.
EPIC is able to simulate multiple crops because it embodies a generic plant model that can be re-parameterized to represent different species.Table 1 describes the subset of these parameters that were estimated for four major regional crops: corn, grain sorghum, alfalfa, and soybean 2 .A priori parameter value ranges are given in Table 2 -the first column contains values suggested in the EPIC Users Manual (Williams et al., 1990).The other columns are based on authors' experience within the study area.The former are wider because of the crop and geographic diversity of EPIC utilization.Although some of the chosen parameter ranges exceed the typical ranges stated in the EPIC documentation, these typical ranges do not reflect the limits that the model is capable of simulating.
Variables that control plant growth and canopy development (and subsequent water use) were selected for optimizations.These were WA, TB, TG, DLAI, RLAD, and RBMD (Table 2).Variables that affect irrigation timing (IRI, BIR, ARMN, and ARMX) were also optimized.Other inputs that might affect hydrology (soil runoff curve and slope) were based on the soil groups.Growth-related parameters such as fertilizer uptake were not altered as plant nutrition was not of interest and simulations were managed in a manner such that nutrient stress did not affect plant growth and canopy development.
The irrigation-related parameters deserve special mention.The dates and amounts of individual irrigation applications are rarely available to modelers.Thus, most crop simulators include "automatic irrigation" options under which irrigations are simulated on dates when preset soil moisture or water stress levels are reached.We sought values for the parameters defining this option that reproduce annual water usage for wells in the study area -in effect attempting to solve for indexes of irrigator behavior.

Calibration data sets
Data on all irrigated land parcels in Kansas are available in a unique database maintained by the Kansas Division of Water Resources (KDWR).Irrigators in Kansas are required to report their yearly water use by parcel to the KDWR.The annual report for each parcel also includes the type of irrigation system, the crop(s) grown, the number of acres irrigated, and the yearly irrigated water volume.These reports are compiled and distributed via a publicly available GIS data product, the Water Information Management and Analysis System (WIMAS, Fig. 2).Although WIMAS data are spatially comprehensive and detailed, they have some shortcomings.
2 Estimation of wheat parameters was deferred because the extra programming complexity required to split activities across two calendar years was not seen as necessary to evaluate the approach.First, several variables needed to identify production relationships (crop yields and the levels of other inputs such as fertilizers and pesticides) are not included.Second, if multiple crops are grown on a given parcel, the irrigator is not required to report the subdivision of acreage or water allocation.Third, although points of diversion are increasingly metered, water use reports from un-metered sites may have significant error.These limitations mean that estimation must be robust in the face of uncertain data (see Sect. 4).
Obtaining good parameter estimates requires multiple years' of data to overcome annual weather variation.WIMAS data are sparse before 1990 so the 11-year period 1990-2000 was used.Historic weather data collected from five weather stations in and around the county were used in the simulations.The simulation of each well used the weather data from the station nearest the well.Sheridan  County soils were combined into two groups using data from the Soil Survey Geographic (SSURGO, http://soildatamart. nrcs.usda.gov)database maintained by the Natural Resources Conservation Service of the United States Department of Agriculture.Group I contains silt loam soils with slopes between 1 and 3%.This group contains ca.76% of the land area and 663 of 779 total wells in the county.The second group includes all soils with slopes greater than 5%.This group consists of loam, silt loam and silty clay loam soils and accounts for ca.23% of the land and 93 wells.The remaining ca.1% of the land and 23 wells were discarded.Thus, each simulation of the crop production associated with a well utilizes one of the five sources of weather data and one of the two soil groups.Of the possible 8316 well-years in the simulation period (756 wells over 11 years), 4931 of them were used in the calibration.The WIMAS data reported that the other 3385 well-years either did not grow a crop or grew a crop that was not one of the four included in the this study (or multiple crops were grown).We ensured that each crop was represented by at least one well in each year.Table 3 shows the distribution of all 4931 well-year combinations and the corresponding breakdown of total irrigation water usage.Irrigation accounted for 98.8% of all water pumped in Sheridan County during the study period and the well-years in this study totaled 81.3% of all irrigation usage.
Because policy analysis will ultimately entail large amounts of computing power (Steward et al., 2009a), it is desirable to understand the relationship between sample size and estimation outcomes.Thus, a ca.ten percent sub-sample of corn well-years was randomly selected from five clusters identified in each soil group via the Getis and Ord (1992) G * i statistic.The well clusters were based on (1) maximum reported volumes of water pumped in any of the 11 years and (2) spatial propinquity (Fig. 3).
For a coupled hydrology-crop-economic model to be useful, it is clearly important to mimic crop yields as well as water use (Steward et al., 2009b).County average annual yield data by crop were obtained from the National Agricultural Statistics Service Quick Stats database (http://www.nass.usda.gov)that contains reports from the Kansas Agricultural Statistics Service.The different aggregation levels for the yield and water use data reflect a frequent occurrence in interdisciplinary regional studies.Aggregation effects on model accuracy have been studied both theoretically and empirically (Theil, 1954;Grunefeld and Griliches, 1960;Zellner, 1969;Aigner and Goldfeld, 1974;Sasaki, 1978;Pesaran et al., 1989).While any level of aggregation is possible, extreme modeling approaches are (1) to aggregate all the data in the study region and perform the analysis at macro-level or (2) to downscale variables available only in aggregate form into many small regions and conduct a micro-level analysis with a unique sub-model for each decision maker.However, the most accurate aggregation level for a real problem must be found empirically as it depends on unobservables like data measurement errors.While cognizant of these issues, we have elected not to investigate them at this time.Instead, we are utilizing an estimation method that (1) does not require all data to be at the same scale and (2) yields unambiguous indicators of parameter uncertainty.

Maximum entropy estimation
Maximum entropy (ME) (Golan et al., 1996) estimation entails maximizing an information theoretic measure of uncertainty (entropy, H ) subject to constraints imposed by data.The results are probabilistic estimates of parameter values that are as certain as the data allow, but no more so.ME equations remain solvable even in cases where sparse data render the corresponding Least Squares (LS) and Maximum Likelihood (ML) equations indeterminate.ME estimation has become increasingly popular in many situations, particularly in models where the data are incomplete because the variables of interest are measured at high levels of aggregation.In the field of production economics, several researchers have invoked ME estimation to recover disaggregated production relationships (e.g., crop yield as a function of field-level inputs) from aggregate data (Howitt and Msangi, 2002;Lence and Miller, 1998;Lansink et al., 2001;Golan et al., 1994Golan et al., , 1996)).The parameter values obtained through the estimation procedure were not estimated for individual fields, as the values were assumed to be representative across the study region, which is relatively homogeneous in terms of soils, water availability, farming practices/technology, etc.
Analytically, EPIC can be represented as a mapping from the combined spaces of model inputs and parameters to a space of outputs.More formally, if there are a total of J input variables (soil characteristics, weather conditions, irrigation amounts, etc.), and all inputs are real numbers, then EPIC inputs can be represented as a vector x∈ J .Similarly, if there are K parameters, each of which is known to lie in an interval with finite bounds, then the parameters are a vector β∈B, a hypervolume in K .For simplicity, assume initially that the only model output of interest is crop yield, y, which is always non-negative.EPIC is then a mapping F y : J ×B→ + , and a yield prediction for a given situation can be written as y=F y (x; β).
The ME procedure estimates the probability distributions of the unknown parameters β.Let z k be an M-dimensional vector of support points along the kth dimension of B, and let p k be the corresponding vector of probabilities; i.e., p mk =Prob[β k =z mk ], m=1, . . ., M. For a given p k , the estimated value of β k can be written as z T k p k , where z T k is the transpose of z k .The simplest specification is where there are two support points for each parameter, corresponding to upper and lower bounds of the known range.In this case, the estimate of the kth parameter is In the general case, the entire parameter vector can be written compactly as β=Zp, where p=(p 1 , . . ., p K ) and Z=diag(z T 1 , . . ., z T K ).If all input data, x ij , and yield data y ij =F y (x ij ; β) were available for each of i=1, . . ., n years and j i =1, . . ., m i wells where K i m i , then β could be estimated via ME, LS or ML (see Welch et al. (2002) or Steward et al. (2008) for an LS example).However, in the current situation, only where H is entropy and ω j is the area irrigated by well j , and 1 M is an M-dimensional vector of ones.This is solvable no matter how sparse the yield data because the constraint that probabilities sum to one is, alone, sufficient to determine a p (namely, the uniform distribution where p mk =(MK) −1 ).The foundations of entropy estimation are described in detail in Golan et al. (1996).The ability to handle limited data differs drastically from LS and ML where the rank of the governing equations is determined by the number of observations and must at least equal the dimensionality of the vector of unknown parameters.
The description just given readily generalizes to the case of multiple dependent variables (here, crop yield, y (t/ha), and water use, w (mm-ha)) or data statistics (e.g., one might also seek to mimic higher moments of inter-annual water use).Additional information is simply added as further optimization constraints.This flexibility to tailor the problem statement to the information available is another reason for ME's increasing utilization.

Optimization methodology
The genetic algorithm (GA) (Goldberg, 1989), an optimization procedure based on Darwinian evolution, was applied to maximize H . GA possess a number of desirable features that make them an attractive optimization method for this study.They are quite robust in the presence of local optima and both theoretical studies as well as simulations with real-world problems suggest that they are quite effective in obtaining very good solutions.In addition, the vast existing literature on the topic allows one to choose from a variety of operators to suit the needs of a particular optimization problem.
For these reasons, GA have been widely employed in the parameter estimation of models such as EPIC.Zhang et al. (2009) found GA to perform well compared to other optimization algorithms in the calibration of the SWAT model.Multi-objective GA have been used in the calibration of this model (Bekele and Nicklow, 2007;Whittaker et al., 2007) as well.GA have also been used to calibrate runoff models such as HBV (Seibert, 2000) and TOPSIS (Cheng et al., 2006) as well as crop (Dai et al., 2009) and crop-related models such as SWAP (He et al., 2007).
GA operates on a population of trial solutions (100 in this study), each of which is a vector of probabilities, p.The population was initially seeded with random vectors that were normalized to become probabilities summing to unity (Fig. 4).In each of 200 iterations (called generations) per estimation run, the solutions were updated by means of two operators, mutation and recombination.Mutation was accomplished by adding a small, random perturbation to each probability in the solution.The sum of all perturbations added to any vector was zero so mutated vectors still summed to one; appropriate limits kept probabilities between zero and one.During recombination, existing pairs of solutions (called parents) were used to create new ones (offspring) according to In the above equation, 1 p and 2 p are the two parent vectors and 1 p , 2 p the offspring.The quantity λ was randomly distributed in the interval [0, 1].
The selection of parents during crossover was done in a manner motivated by Darwinian survival of the fittest.To obtain each parent, two candidate solutions were drawn randomly from the population and the fitter one chosen.The fitness of any solution, p, was where the first term is the entropy and c(p) was a function that penalized solutions which violated any data constraint.Although not uncommon in mathematical optimization, penalty functions must be carefully designed.Substantial constraint violations may result when penalties are too mild.Excessively harsh penalties can prevent the computer from finding any solution, even when one exists.Additionally, (1) variation in the numerical magnitudes among model outputs will differentially affect the penalty and (2) the investigators may prefer to predict some variables more accurately than others.To ameliorate differences in numerical magnitudes, the penalty function was defined as the weighted sum of the relative absolute errors where (η w /η y ) is the ratio of investigator preference for errors in predicted water use over errors in yield and λ H scales the penalty relative to the entropy.Because c(p) penalizes the closure error of what are intended to be equality constraints, λ H must be as large as possible while still allowing the entropy to influence the optimizer.Thus λ H was chosen so that entropy accounted for 5% of the total fitness.Although our procedure imposes only two constraints (zero yield and water-use error) on 10 free parameters, an exact solution is not possible because the constraints pull the free Hydrol.Earth Syst.Sci., 13,[1467][1468][1469][1470][1471][1472][1473][1474][1475][1476][1477][1478][1479][1480][1481][1482][1483]2009 www.hydrol-earth-syst-sci.net/13/1467/2009/  variables in opposite directions: solutions that satisfy the yield constraint result in high water error (and vice-versa).For this reason, we sought to minimize and balance the error of both constraints.An arbitrary, but not unreasonable value for (η w /η y ) is one for which the relative error in water use, when aggregated from the well to the county level, equals that of yield.To identify a suitable value, a series of 20 estimation runs with randomized initial conditions was performed for each of 14 different values of (η w /η y ) using the ten percent sample corn data.Figure 5 shows that (η w /η y )=14 is an appropriate weight ratio; it was used in all subsequent runs.
The entire investigation entailed 3.841 billion executions of the EPIC odel; the pilot study to set (η w /η y ), alone, required 616 million simulations.Such numbers are ca.three orders of magnitude greater than those in studies reported just a few years ago (e.g., Irmak et al., 2000;Welch et al., 2002).Computation of this scale demands the use of highperformance computing.The GA was designed in a masterslave parallel fashion (Cantu-Paz, 2000) and implemented as a scalable system that hybridized several software components.The interdisciplinary discussion and design was facilitated by writing the GA in a high-level mathematical scripting language (Scilab, http://www.scilab.org/).On the other hand, parallel execution of the model was coordinated by a client written in C to achieve high-performance.The model itself is a legacy Fortran code.The system executes on both dedicated clusters via MPI (Graham et al., 2006) and in a loosely-coupled, distributed fashion via Condor (Thain et al., 2005).The simulations were performed on a 200 CPU Beowulf cluster at Kansas State University and a 200 node Condor pool at the University of Oklahoma.Software performance measures and scalability were reported in (Bulatewicz et al., 2007).There are several general questions to ask in a parameter estimation study of this type.First, what are the resulting estimates and what can be said about their uncertainty?Second, how reasonable are the results in terms of both the individual estimates and their interrelationships?Third, ME integrates all sources of information in determining its results, including both the data as well as the prior information available to the investigators and expressed in the initial ranges set for the parameters.What has been the relative influence of these two factors on the outcome?Fourth, what level of predictability has been achieved?

Estimates and reasonableness
To address these issues in an integrated way, a bootstrap procedure was used (Efron and Tibshirani, 1993).For each of 250 replicates, a random sample of 11 years was selected with replacement from the data and a set of parameter estimates obtained.On average, each replicate contained ca.seven unique years.To ensure that the variation within the resample structurally reflected that within the original data, the wells in each year were further re-sampled by soil type and weather station.This was done for all crops, including separate runs for the ten percent sample and complete corn data sets.Averaging the best estimates of the 250 replications gives the final estimates (Tables 4 to 7).The standard deviations of the 250 estimates are the bootstrap standard errors of the parameters.The shape of the probability distribution for each parameter can be approximately visualized by plotting a histogram of the estimates from the individual bootstrap replications.The histograms indicate the number of times that a parameter value was estimated to be in a given range out of 250 bootstrap replications for each crop .
An immediate result is that the parameter values estimated using the complete corn data were found to closely match the values estimated using only the ten percent sample of corn data.The difference in each estimated value was less than 1% for most parameters, with the largest difference being 5% for the maximum volume per irrigation (ARMX) parameter.
Optimization of the three variables that are associated with irrigation system management (IRI, ARMN, and ARMX, Table 1) resulted in similar values across the four crops (Tables 4 to 7).The minimum application interval (IRI) for all crops was approximately 10.3 d (Tables 4 to 7) and is longer than typically experienced under current production and environmental constraints.Given sufficient data, EPIC applies reasonable total amounts of water on a countywide basis (e.g. corn in Table 8), but appears to do so in fewer, larger, applications.The average well capacity for all wells 45 ha.Using these values, it is possible to apply an irrigation depth of 37 mm every 5.4 d.A common practice is for constant irrigation during critical growth stages so as to maximize yields.Additional model runs indicate that EPIC outputs are not par-ticularly sensitive to IRI in the range of 4 to 10 d, but sensitivity increases above this point (data not shown).A shallow minimum in f (p) slightly above 10 d could therefore account for our results.To the best of our knowledge, no prior studies have calibrated IRI so this model feature has not previously been detected.It is potentially important in that it suggests that irrigation scheduling could possibly reduce water use while maintaining yields.The minimum and maximum volumes for automatic irrigation (ARMN and ARMX, respectively) resulted in mean values across all four crops of 10.6 and 37.6 mm per application, respectively.While realistic in terms of typical irrigation practices, these limits are broad enough to encompass a model shift from more numerous smaller applications to fewer, larger ones.
The physiologically-based variables followed trends previously reported from field and greenhouse experiments as well as other crop modeling efforts.The water stress level to trigger irrigation (BIR) is specified in terms of biomass production: irrigation occurs on days where the ratio of biomass produced to potential production (given adequate water) falls below BIR.The BIR values for corn, grain sorghum, alfalfa, and soybean are 0.86, 0.87, 0.87, and 0.85, respectively -not unreasonable for irrigated cropping conditions where water stress is less likely to occur.
Optimized parameters that not unexpectedly had speciesspecific ranges are biomass to energy conversion ratio (WA), optimum temperature for growth (TB), and minimum temperature of growth (TG).Optimized values for WA were 47.0 for corn (complete data), 33.4 for grain sorghum, 29.2 for alfalfa, and 31.2 for soybean.Reported corn WA values ranged from 14.5 to 43.3 t ha −1 MJ −1 m −2 (Cantarero et al., 1999;Muchow, 1990;Hammer et al., 1998;Kiniry et al., 2004;Idinoba et al., 2002;Sinclair and Muchow, 1999) with the range being attributed to differences in environmental conditions and calculation method.Reported WA specifically for use in crop models for corn have ranged from 43.3 in CERES-Maize (Yang et al., 2004), 39.8 for ALMANAC (Kiniry et al., 2004) and 35.0 to 40.0 t ha −1 MJ −1 m −2 for CropSyst (Stockle et al., 2003).
Alfalfa WA measurements and its inclusion in simulation suites are limited compared with corn, grain sorghum or soybean.Collino et al. (2005) reported WA values between 12.0 and 15.0 t ha −1 MJ −1 m −2 for field grown alfalfa in Argentina and Whitfield et al. (1986) reported a WA value of 17.2 t ha −1 MJ −1 m −2 .Confalonieri and Bechini (2004) used 30.0 for a WA value after calibrating CropSyst for use in northern Italy.Although higher than reported field values, this WA value is nearly identical to the 29.2 value that resulted from our optimization process.Soybean WA values from field experiments have been reported to be 20.0 (Sinclair and Muchow, 1999).CropSyst initial WA values are reported as 20.0 to 25.0 t ha −1 MJ −1 m −2 which is lower than the WA of 31.2 that resulted from the optimization process we used.The potential radiation use efficiency (WA) in EPIC is assumed to be for unstressed plants and includes root growth, which are often cited as reasons for field measured values being lower than those finally published as being used in most crop simulation models.Soybean WA was the only optimized value that was higher than reported for use by simulation models.The estimated optimum temperatures for growth (TB) were 27.2 • C for corn (complete data), 32.5 • C for grain sorghum, 30.4 • C for alfalfa and 29.1 • C for soybean.Several researchers have reported optimum temperatures for growth in corn to be 22.5 • C (Wilhelm et al., 1999) and 25 • C ( Grzesiak et al., 1981).CERES-Maize uses 26 • C as the optimum temperature for growth (Jones et al., 1986) while Yang et al. (2004) currently use 30 • C in the Hybrid Maize model for maximum growth and assimilation.Our optimized corn temperature for optimum growth is within the range of the values used by other simulation models and slightly higher than those reported from research trials.Our estimated TB for grain sorghum (32.5 • C ) agrees with the results of Prasad et al. (2006) and Chowdhury and Wardlaw (1978) who reported optimum temperatures for growth in grain sorghum to be 32 and 30 • C, respectively.Our TB value is lower than the 44 • C that is currently used in SORKAM (Rosenthal et al., 1989), a sorghum simulation model.
The TB estimate of 30.4 • C for alfalfa is higher than the reported values of 20 to 27 • C cited by several others (Arbi et al., 1979;Fick, 1984;Bula, 1972;Ueno and Smith, 1970;Pearson and Hunt, 1972).However, our TB value is similar to the 30 • C used by Confalonieri and Bechini (2004) in their optimization of CropSyst.Our estimate of 29.1 • C for soybean is higher than the 24 to 27.5 • C reported by Seddigh and Jolliff (1984) and Ghazali and Cox (1981), but is in agreement with Gibson and Mullen (1996) and Grimm et al. (1994) who reported optimum temperatures for photosynthesis and growth to range from 29 to 35 • C. Our TB is less than that used in CROPGRO-Soybean, which uses 40 • C as an optimum temperature for photosynthesis (Pedersen et al., 2004).
The estimated minimum temperatures for growth (TG) were 8.2 • C for corn (complete data), 6.2 • C for grain sorghum, 0.5 • C for alfalfa, and 10.7 • C for soybean.These estimates are similar to those reported by others from field or growth chamber research or currently being used in other crop simulation models.Reported corn TG values range from 7.2 to 8 • C (Hesketh and Warrington, 1989;Yang et al., 2004).For grain sorghum, a TG of 8.5 • C is reported by both Craufurd et al. (1998) and Hammer et al. (1989) while SORKAM (Rosenthal et al., 1989) uses 7 • C as a base temperature, all of which are marginally higher than the 6.2 • C estimated here.
Estimates of TG for alfalfa of 0.5 • C are lower than the 5 • C used by Confalonieri and Bechini (2004) in calibrating CropSyst for simulating alfalfa in Italy.Our TG of 10.7 • C for soybean is greater than that used in CROPGRO-Soybean, which uses 8 • C as base temperature for photosynthesis (Pedersen et al., 2004).
Estimates of when leaf area begins to decline (DLAI), the rate of decline (RLAD) and the rate at which WA declines (RBMD) could be deemed reasonable based on typical production scenarios.For all four crops, the values for DLAI indicate that leaf area begins to decline after around 86% of the growing season has occurred.Leaf area typically begins a gradual decline in corn and grain sorghum within approximately 1 week of anthesis, but this largely occurs in the lower canopy and the loss represents older leaves that will not contribute to final yields.However, more rapid leaf loss near physiological maturity typically happens in all three grain crops and appears to be captured appropriately in our results for these three variables.In alfalfa, it is difficult to determine how these values compare to reality since this perennial crop does not normally senesce as it is typically harvested at 10% bloom stage.The final harvest timing recommendation is after cold temperatures induce dormancy in the crop.

Influence of prior information
Maximum entropy estimation makes use of the complete corpus of available information.Thus, the standard errors reported in Tables 4 to 7 reflect notonly the data but also the prior information implicit in the selected probability support points.As just documented, the optimizations reported herein did produce reasonable output values.Even so, it is possible that the data do not constrain the estimates either because (1) they are too fragmentary or (2) the parameter's influence on actual outcomes is too weak or indirect.It is therefore useful to ask (Q1) do the data detectably influence the outputs and (Q2) how strongly does prior information affect the estimates?
If a parameter has no influence on the model outputs then it cannot affect the penalty function values.In this situation f (p) is optimized when H is maximized.In a two-point distribution this happens when the parameter estimate is the midpoint of the support interval no matter where the endpoints are located.Based on this fact, metrics for Q1 and Q2 were developed and applied to the ten percent corn sample.The Q1 metric calculates a two-tailed, binomial distribution p-value for the null hypothesis that the median of the 250 re-sampled parameter estimates is the midpoint of the support range.A p-value of less than 0.05 is interpreted as a significant data influence.The easiest way to measure dependence on the support point choice is to make a different choice.Thus, the Q2 metric is a finite difference estimate of the derivative of the parameter estimate with respect to the midpoint of the parameter range.Values close to unity are consistent with the estimate being completely dependent on prior information.Ten sets of 20 estimations were run with the range of a single parameter changed in each set.The ranges were shifted up or down depending on whether the majority of the original estimates did or did not exceed the midpoint.Range widths were preserved unless doing so resulted in an endpoint that was (1) outside the range suggested by EPIC, or (2) conflicted with the range of another parameter.
The results are in Table 9.All but three parameters (DLAI, RLAD and RBMD) are influenced by the data with more significant median departures from the support interval midpoints being generally associated with lower sensitivities to prior information.It is clear, however, that the Q1 and Q2 metrics measure different properties of the estimation process, as illustrated by ARMN which responds both to prior information and to data.TG shows a similar pattern but with a lower dependence on prior information.The unit sensitivity to prior information displayed by RLAD and RBMD coupled with their close adherence to the support interval midpoint, suggests that these parameters are poorly constrained by the data.

Model verification
To what extent can estimates obtained by the methods reported herein be used to predict outcomes in other situations?This was analyzed in several ways.First, simulations were performed using each of the 250 bootstrap estimates to see how well they could reproduce the complete data and how Hydrol.Earth Syst.Sci., 13,[1467][1468][1469][1470][1471][1472][1473][1474][1475][1476][1477][1478][1479][1480][1481][1482][1483]2009 www.hydrol-earth-syst-sci.net/13/1467/2009/ well they could predict the years that had been omitted in each replicate ("Fitted" and "Predictive", respectively, in Table 8).Depending on the crop, the variable (yield or water use) and, for water, the level of aggregation (county-vs.well-level) attempting to predict independent, out-year data increased the error from 3 to 9% of full scale.The predictive water use errors were quite high at the well level ranging from 35% for the complete and ten percent corn sample to 49% for sorghum.However, for aquifer modeling, accuracy at this fine a geographic scale may not be needed.Predictive errors at the county level were smaller.Second, we used the mean parameter estimates obtained from the ten percent corn sample to estimate the behavior of the remaining 90% of the data.The performance was essentially identical, confirming the ability of a small sample to enable accurate estimation of a much larger set from the same years.The selection procedure for the ten percent sample was designed to achieve good representation.There are two ways this could have failed: the estimates might have deviated from those of the larger sample or the full data set might have had much greater variability.The latter would have resulted in larger relative errors even if the mean predictions were accurate.Clearly neither of these happened.However, the bootstrap relative errors show that the inter-annual variability is such that samples of 7 out of 11 years can only capture it to the degree shown in the Predictive columns, at least at the single county scale.Of course, incorporating data from more than one county may serve to offset temporal variation.Finally, the relative error in total county irrigation use for these four crops is 13%, which is the same as that for corn, due to its heavy preponderance in the county.
Third, the county-level yield and water use errors are within 5% of each other in all cases but sorghum, in which the error difference was 14%.This suggests that a penalty ratio of 14 was effective in balancing thecounty-level yield and water use error, but that crop-specific penalty ratios may give better results.Finally, whereas relative errors are dimensionless fractions, Table 10 reports county-level yield and water use errors in physical units.The two left columns are the root mean square errors (RMSE) resulting from running simulations using each of the 250 estimates against the complete data.The two right columns are the root mean square errors of prediction (RMSEP) and are calculated in the same way except that the 250 estimates were used to predict the years that had been omitted in each replicate.These latter numbers may be taken as an indication of expected model performance using data collected from a single county.These RMSE of water use and the well-level relative error (Table 8) both reflect the accuracy of the simulations at the field scale.When aggregated to the county scale, there is a significant reduction in the error (compare well-and county-level relative error of water use in Table 8).

Conclusions
Informatic technology now provides the means to interweave information sources (data and models) in ways conducive to integrative hydrologic investigations.At one level, the work reported here relates to fitting an existing model to a set of data.But, from another perspective, it prototypes www.hydrol-earth-syst-sci.net/13/1467/2009/ Hydrol.Earth Syst.Sci., 13, 1467-1483, 2009 T. Bulatewicz et al.: Calibration of a crop model to irrigated water use the integration of a range of methodologies whose crossdisciplinary fusion enables a much broader range of activities.This integration reaches across diverse disciplines: agronomy, agricultural economics, civil engineering, computer science, and electrical engineering.The specific technologies are correspondingly diverse: crop physiological simulation, maximum entropy estimation, genetic algorithmbased optimization, alternative parallel processing methodologies, GIS databases, and cluster analysis, to name a few.Through the integration of these technologies, a genetic algorithm that employed ME was able to identify realistic values for ten crop model parameters including both hydrological and physiological inputs.The estimated parameter values were not only realistic, but they were able to reproduce observed irrigation water use to within 13% with an estimated predictive accuracy of 19% (complete corn data) given sufficient data.In addition, several of these parameters represent irrigation practice, effectively recovering detailed irrigator behavior from annual water use reports.Knowledge of this kind is otherwise unavailable yet is essential for the prediction of water use in alternative scenarios as part of an investigation into effective water policy for sustainable usage.Such a policy analysis would entail considerable computational requirements.These demands could be mitigated through the use of carefully selected samples of data with only slight increases in error, as demonstrated in this work through the evaluation of the ten percent sample of corn data.The implementation of these techniques has provided a framework within which additional models and data will be integrated.The particular team assembled for this paper crosses the spectrum of hydrologists, agronomists, economists, and computer scientists, and the results are being shared and translated by disciplinary specialists to interested collaborators, stakeholders and agencies with which the team is working.Examples include applications of the calibrated model in both irrigation studies and in the assessment of the economic impacts of water policy on farmers in western Kansas as part of an integrated model.Bordogna (2003) wrote: "The trick in evolving the capability for providing emergent infrastructure is the selective use of established models and the rapid generation and testing of new models.This is a process of institutional and organizational learning, and the ability to learn rapidly is itself a kind of social infrastructure that is required to pursue the cyberinfrastructure vision."We believe that studies such as this one are important steps along that learning path.

Fig. 1 .
Fig. 1.Plants produce carbohydrates from carbon dioxide, water, and light energy.They grow and develop at rates that are nonlinearly dependent on resources and temperature.All but ca.one percent of water use is for transport or cooling and is transpired.

Fig. 2 .
Fig. 2. Well locations in Sheridan Co. KS.Wells are visually coded to show the crop grown in 2000.The most common crop is corn.

( a )
Maximum water pumped by each well.(b) Clustered wells.A banded northwest-to-southeast (greater to lesser) trend in pumping is evident.

Fig. 3 .
Fig. 3. Wells in Sheridan Co. KS on a map of Group I (largely silt loam) soils.

( a )
Ratios ranging form 0.01 to 100.(b) A close up of the crossover range.

Fig. 5 .
Fig. 5. Mean relative error for the ten percent corn sample as a function of the preference ratio between county-level water (dashed line) and yield (solid line) residuals.

Table 1 .
EPIC parameters to be estimated.

Table 2 .
Parameter ranges used in the estimations.

Table 3 .
Distribution by crop of well-years by soil group and water use.

Table 4 .
Parameter estimates for corn.

Table 5 .
Parameter estimates for alfalfa.

Table 6 .
Parameter estimates for sorghum.

Table 7 .
Parameter estimates for soybean.

Table 8 .
Relative errors.Not separately bootstrapped.b Different grains are not commensurate.c Single wells are single crops.d Analysis would require averaging (number of bootstrap reps) (number of crops) scenarios (ca.3.9×10 9 ).

Table 9 .
Analysis of the influence of prior information on corn parameters.