| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
|
|
|
|||||||
|
|||||||||
a USDA-ARS-JPCNRCC, Watkinsville, GA 30677
b Dep. of Crop and Soil Sci., Univ. of Georgia, Athens, GA 30602-7274
c Usda-Ars-Grl, El Reno, Ok 73036
d Usda-Ars-Ngprl, Mandan, Nd 58554
e USDA-NRCS-ITC, Fort Collins, CO 80526-8121
f Dupont Research, Wilmington, DE 19898
g Biol. and Agric. Eng. Dep., Univ. of Georgia, Griffin, GA 30223-1797
* Corresponding author (dstark{at}uga.edu)
Received for publication June 15, 2004.
 |
ABSTRACT |
|---|
 TOP ABSTRACT INTRODUCTIONMATERIALS AND METHODSRESULTS AND DISCUSSIONCONCLUSIONSAPPENDIXREFERENCES |
|---|
soil water in a 60-cm profile during the critical peak bloom period. A detailed analysis of the calibration procedure and parameters used in this study will aid subsequent users of the model as well as aid in a subsequent evaluation of the model's performance for simulations of tile drainage and nitrate leaching in Georgia Piedmont cotton production systems.
Abbreviations: CI, confidence interval • CT, conventional tillage • EF, tile drain express fraction • LAB, sorptivity factor for lateral infiltration • MSEA, management system evaluation areas • NT, no tillage • PTFs, pedotransfer functions • RRMSE, relative root mean square error • TDR, time domain reflectometry • WT, wetting thickness
|
INTRODUCTION |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSRESULTS AND DISCUSSIONCONCLUSIONSAPPENDIXREFERENCES |
|---|
Calibration of a model includes parameterization based on direct measurements, pedotransfer functions, or direct or indirect fitting of the model to measured data. Pedotransfer functions (PTFs) provide estimates of parameters such as soil hydraulic properties that are often difficult and time-consuming to measure but accurate enough for many applications (Pachepsky and Rawls, 2005). However, Vereecken et al. (1992) showed that >90% of the variability in simulations of a soil map unit was due to the variability in the estimated hydraulic parameters using PTFs.
Recently, some modeling studies have begun to provide more details on the calibration approach or procedure that was used (Abrahamson et al., 1999; Cornelis et al., 2004; Zhang, 2004; FAO, 2004). This is due in part to the recognition of the need for standardization of calibration procedures, and subsequent guidelines that have been developed (Clarke et al., 1994; Hanson et al., 1999; Dubus et al., 2002; Saseendran et al., 2003; Bouman and van Laar, 2004). Though the modeling process can be defined procedurally, processes such as calibration and validation are completely subjective and open to best professional judgment, and modelers have no obligation to meet a standardized set of criteria (Corwin et al., 1999). A lack of emphasis on the process used for calibration may have resulted in assumptions or conclusions by readers and subsequent users of a model that may or may not be accurate. It may not be clear in the reported documentation if parameters were based on measured data, if parameters were adjusted during calibration, if all major processes in the model were parameterized, or if sensitivity analyses were performed. The lack of reporting of the calibration process may leave a reader with limited information to discern the model's ability to comprehensively address the system tested. Adjustments made to parameters during calibration may impact other processes in the model that do not concern the current modeler, but may not be suitable under different conditions that would be of interest to another modeler.
This modeling study is based on a current water quality field experiment initiated in 1991 at the USDA-ARS J. Phil Campbell, Sr. Natural Resources Conservation Center in Watkinsville, GA. Objectives of the study included the water quality impacts of maize production based on the effects of conventional tillage (CT) or no tillage (NT), cover crop, and nutrient source. A model that could accurately simulate the sensitivity of drainage and nitrate leaching to these management practices would provide a valuable tool for testing and evaluating different agricultural production scenarios in Cecil and associated series soils which occupy approximately two-thirds of the cultivated land in the southern Piedmont region (Hendrickson et al., 1963).
Using this same study, Johnson et al. (1999) tested the LEACHN model (Hutson and Wagenet, 1992) for maize production for its ability to simulate soil NO3–N and NH4–N content, and drainage and leached nitrate under CT or NT management with and without a winter rye cover crop. Using modifications based on laboratory estimates for input parameters, LEACHN generally underestimated soil NH4–N and NO3–N during the winter and overestimated soil NH4–N during the summer. The model also overestimated cumulative drainage and leached nitrate during both seasons (Johnson et al., 1999). The overestimation of leached nitrate in a wetter than normal year was attributed to the absence of a soil macropore–matrix exchange component in the model. We chose to evaluate The Root Zone Water Quality Model (RZWQM) because it includes a macropore component as well as an exchange component between the soil matrix and macropore walls. Visible macropores and preferential flow patterns are found in Cecil soils (Gupte et al., 1996), and we expected that the RZWQM might be able to better simulate drainage and leached nitrate based on the results of the Johnson et al. (1999) study. In addition, the RZWQM includes an option for tile drainage, an important consideration when tile drains have been used in the field due to changes in the soil water dynamics caused by artificial drainage systems (Skaggs, 1978; Ritzema, 1994).
The hydrology, pesticide, and nitrate movement, crop growth, and several agricultural management practices in the original version of the RZWQM model published in 1992 have been tested nationally and internationally with data collected from 1972 to 1996 (Ahuja et al., 2000). Tillage effects on hydraulic properties, manure management, crop yield response to water stress, and tile drainage are just some of the refinements present in the version of the model used in our study (USDA-ARS-GPSR RZWQM development team, personal communication, 2004). Conclusions drawn from some of the early applications in the literature may not be strictly valid, and may not represent typical behavior of the current model (Ma et al., 2001). In addition, soils and climate in the Southeast are very different from the Great Plains and Midwest regions of the USA. This paper reports results of the calibration of the most recent version of the RZWQM (v. 1.3.2004.213) for simulations of tile drainage and nitrate leaching in maize and cotton production with a winter rye cover crop as well an analysis of the effect of macroporosity on tile drainage and nitrate leaching under conventional tillage management practices in the southeastern USA.
The main objective of this study was to calibrate the RZWQM for its ability to simulate tile drainage and nitrate leaching in a Cecil soil in maize and cotton production with a winter rye cover under conventional tillage management practices in the Georgia Piedmont region. A second objective was to evaluate the model's sensitivity to soil macroporosity in relation to tile drainage since regions of preferential flow are found in Cecil soils of the Piedmont region (Gupte et al., 1996). Finally, we aimed to provide a detailed explanation of our calibration procedures for other modelers and user groups who are interested in the process of calibration that might be useful before model evaluation. Clarification of calibration procedures provides a better understanding of the parameterization process, and the sensitive parameter adjustments that are discovered during the process. It may have implications for potential users of the model if any specific parameters or parameter adjustments have significantly influenced test results. In addition, this study contributes toward the standardization of the calibration phase of modeling. A standard calibration protocol supplements the current protocol of parameterization, calibration, and testing with an independent data set, with guidelines that for now are left somewhat arbitrarily up to the modeler (Dubus et al., 2002).
|
MATERIALS AND METHODS |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSRESULTS AND DISCUSSIONCONCLUSIONSAPPENDIXREFERENCES |
|---|
The treatment plots were 10 by 30 m and instrumented with 10 cm i.d. PVC drain tiles installed 2.5 m apart at 75- to 100-cm depths on a 1% slope. The plots were hydrologically isolated from each other with polyethylene sheets extending from the soil surface to a depth of 1 m and with plastic borders 10 cm deep. The volume of water drained from a plot was measured by tipping bucket gauges and digitally recorded by automated dataloggers. The tipping bucket gauges had a sampling slot that subsampled drainage and routed it to a beaker. For every 2 mm of cumulative drainage, a sample was pumped from the beaker into a polyethylene bottle inside a refrigerated sequential waste water sampler (Isco Model 3700 FR, Lincoln, NE). An aliquot of this effluent was stored frozen in polyethylene vials and later analyzed for nitrate using the Griess-Ilosvay method (Keeney and Nelson, 1982). The samples were filtered through a 0.45-µm filter before analysis (McCracken et al., 1995; Johnson et al., 1999).
The soil was a Cecil sandy loam. The pH normally ranged from 5.5 to 5.8 as measured at the study site; therefore, lime was applied approximately every 3 yr to maintain a pH of 6.0 to 6.3 in the surface horizon to avail plant nutrients and prevent aluminum toxicity. Since these soils are variably charged, positively charged soil particles can attract anions such as nitrate that can be weakly held in the soil matrix. Nitrate may bypass the soil matrix via soil macropores. However, Gupte et al. (1996) found regions of preferential flow in dye-stained soil cores from the study site that were not necessarily associated with distinct open macropores observed from the mean cross-sectional areas of the soil columns.
In April 1991, the plots were plowed, disked, and planted to maize. In October 1991, maize was harvested, and six plots were no-till planted to rye and six plots left fallow through the winter. In April 1992, three plots from each of the rye cover and fallow treatments were placed under either CT or NT management. The CT plots were mowed, moldboard plowed, and disked. On 24 Apr. 1992, plots were planted to maize in 76-cm rows at the rate of 60000 seeds ha–1. Ammonium nitrate fertilizer was applied at a rate of 168 kg N ha–1 on 26 Apr. 1992. Maize was harvested on 7 Oct. 1992 and rye was planted on 30 Oct. 1992. Rye was sampled and killed with paraquat (1,1'-dimethyl-4,4'-bipyridinium ion) on 12 Apr. 1993, CT plots were plowed and disked on 13 April, and maize was again planted on 14 Apr. 1993. Maize was harvested on 14 Sep. 1993 and rye was planted on 29 Sep. 1993. Maize and rye yields and N uptake were measured from biomass samples before each field harvest (McCracken et al., 1995). The same procedure of planting maize followed by winter rye was used until Nov. 1994 when winter wheat was planted as the cover crop followed by the first cotton crop in May 1995.
To calibrate the RZWQM for cotton growth, and its ability to simulate tile drainage and leached nitrate from cotton production during the period when the water quality study was planted to cotton, we used parameters from a field experiment in cotton production in 1997 adjacent to the water quality study. The calibration study site was planted on 16 May 1997 on a 1.3-ha watershed using a no-till drill. A winter rye cover crop was planted in late October following cotton harvest. Soil moisture was measured in 15-cm increments to a soil depth of 90 cm using time domain reflectometry (TDR) (Moisture Point, ESI, Victoria, BC, Canada). Ammonium nitrate fertilizer was applied after cotton planting at a rate of 67 kg N ha–1, and winter rye was fertilized after planting with 54 kg N ha–1. Cotton biomass and leaf area samples were collected seven times during the growing season beginning on 16 July 1997 through 23 Sep. 1997. Plant height and populations were also estimated at each sampling date (Schomberg and Endale, 2004).
Model Input and Parameters
The RZWQM model uses a Windows interface and can initially be set up with a minimum dataset using readily available data. The required soil properties are texture and bulk density. Parameters for soil crusting, macroporosity, tile drainage, and various soil hydraulic properties can be supplied by the user or, where data are limited or unknown, the model will use default values based on known research documented in an extensive user help utility. The model has been applied to simulate best management practices for the Management Systems Evaluation Areas (MSEA) research project for maize, soybean, and wheat (Ahuja et al., 2000). The calibrated maize, soybean, and wheat crop parameters in the model can be adjusted during the calibration procedure to simulate crop growth for the area of interest to the modeler. Other crops may be added to the generic plant growth submodel and parameterized by the user. Daily weather data can be generated with the CLIGEN stochastic model (USDA-ARS, 2003) based on nearby historic weather station parameters when measured data is not available. However, we used measured rainfall and weather data from the Georgia Environmental Monitoring Network for Watkinsville located approximately 15 m from the study site (Hoogenboom, 2003).
We parameterized the physical properties of the soil in the RZWQM model from measurements made near the study site by Bruce et al. (1983) and Gupte et al. (1996). Seven distinct layers to a depth of 1.25 m were parameterized based on measured properties of each layer. The initial soil water content at the beginning of the simulation period on 1 Jan. 1991 was set to the measured field capacity for each layer (Table 1). The van Genuchten (1980) equation parameters, 
and n were fitted using PROC NLIN (SAS Institute, 2000) based on measured soil water content and pressure head for each depth where residual 
was estimated as that of the wilting point at h = –15000 cm. The parameters were then converted to the Brooks-Corey parameters, S2 and A2, the bubbling pressure and pore size distribution index, respectively, based on the RZWQM documentation (Ahuja et al., 2000). We included a soil crusting option with a crust hydraulic conductivity rate set to 0.68 cm h–1 based on measurements of a Cecil sandy loam crust under simulated rainfall conditions (Chiang et al., 1993). The initial soil NO3–N and NH4–N concentrations used are described in Johnson et al. (1999) from soil data collected from the study site in November 1991. We used 1 Mg ha–1 as the amount of initial surface residue based on fallow conditions and on one season of maize production before the first winter rye cover crop in October 1991. The fraction of surface residue mass that would be incorporated by natural means was set to 2% based on model references. The field area used was 0.03 ha based on the size of a plot, and the slope was 2%. Input data and initial parameter values used are listed in the Appendix.
|
Soil Hydraulics
The soil profile can have up to 12 distinct horizons. The profile can be parameterized based on distinct horizons or as distinct layers within horizons. Three grids are then created—one for defining hydraulic properties, a second nonuniform layering system for redistribution of water and nutrients, and a third 1-cm grid that only functions during infiltration. Hydraulic properties in the model are defined by the soil water content–matric suction relationship and the unsaturated hydraulic conductivity–matric suction relationship described by Brooks and Corey (1964) with slight modifications. The model estimates soil hydraulic properties from soil texture, bulk density, and soil water content at a suction of 33 or 10 kPa when measured data for Brooks-Corey parameters are not available. If soil water content at a suction of 33 or 10 kPa is unknown, the parameters for the hydraulic function properties are taken from Rawls et al. (1982) based on the soil texture class and then adjusted based on bulk density. The user has the option of using a minimum description of these properties or a full Brooks-Corey description to account for the effects of trapped air in the soil, which can reduce Ks by as much as 50% during infiltration (Bouwer, 1969). The field saturated Ks is divided by a viscous resistance correction factor of 2.0 so that the infiltration rate at any given time is a function of this reduced Ks in the Green-Ampt infiltration equation (Green and Ampt, 1911). Between rainfall events, soil water is redistributed using the Richards equation minus a sink term for root water uptake and tile drainage flux. These terms are described in more detail in other sections of the paper.
The model includes an option to define soil macroporosity in terms of size and number of macropores present in the soil. The user supplies the macropore number and size (radius) for each soil layer. If data on macroporosity are unavailable, it is best to run RZWQM assuming no macropores (Ahuja et al., 2000). If the user does not select the macroporosity option, then drainage occurs through the soil matrix only.
Water can only enter the macropores at the surface, and the model allows preferential flow through macropores to go directly to the tile drain when the water table resides above the tile drains. Macropore flow may also exchange the soil solution with the soil matrix by miscible displacement through macropore walls. The water solution in the macropore is subject to lateral absorption into the drier soil matrix, and the chemicals in solution are also subject to adsorption or desorption from the macropore walls. Maximum flow-rate capacity of macropores is calculated using Poiseuille's law assuming gravity flow. Lateral absorption into the macropore walls is simulated using Green-Ampt equations (Green and Ampt, 1911; Childs and Bybordi, 1969; Hachum and Alfaro, 1980). The user may also adjust the fraction of microporosity in each soil layer though to not less than 1% of total porosity.
Other than measured values of macroporosity including macropore size and number, the adjustable parameters in the model that can affect macropore flow are the sorptivity factor for lateral infiltration, the effective lateral infiltration wetting thickness, and the tile drain express fraction. To account for the effect that compaction or lining of macropore walls may have in reducing the ability of a soil to absorb water and chemicals, the calculated Green-Ampt radial (lateral) infiltration rate or sorptivity rate will be multiplied by a user-specified sorptivity factor ranging from 0 to 1. The lateral infiltration wetting thickness into a macropore wall can be adjusted to a value between 0 and 2 cm. The tile drain express fraction can be adjusted to a value between 0 and 0.1 to vary the percentage of macropore flow that follows the path into the tile drains and is not subject to absorption into the soil matrix.
Tile Drainage
If the user chooses to simulate tile drainage in the model, flux out of the drains will occur when the water table in the soil profile is above the depth of the drains. The depth of the water table is defined as the depth at which the pressure head first becomes negative, and all heads below that depth are non-negative. When tile drainage is selected, the system will automatically set the bottom boundary condition for the Richards equation to a constant flux condition described by the Buckingham-Darcy equation (Buckingham, 1907) where the total head is the sum of the matric potential and gravitational heads, h + z, in the form:
 |
Lateral flow to tile drains can introduce error in the measurement of unsaturated zone parameters. However, Radcliffe et al. (1996) found that tile drain breakthrough curves can be used in Cecil soils to determine field-scale unsaturated zone transport parameters if a model accounts for two-dimensional flow in the saturated zone. The drainage rate to the tile drains in the RZWQM is calculated according to the Hooghoudt equation as applied by Skaggs (1978) to correct for the two-dimensional flux to the drains. The RZWQM adds the flux to root uptake to become a sink term at the equivalent depth of the drains.
There are two restrictive layers in the Cecil soils at the study site beginning at depths of 35 to 40 cm and at depths of 85 to 90 cm (Radcliffe et al., 1996). We set the tile drain depth in the model to 80 cm, which places them in the middle of the 30-cm soil layer that resides directly above the second restrictive layer. The model calculates the effective depth of the tile drains by calculating effective lateral hydraulic conductivity using the Ks of the soil layer where the drain resides as well as the layer beneath the drain layer to represent the transmissivity of both layers. Figure 1 depicts how we implemented the tile drainage system at the study site in the model to best represent the soil profile and tile drainage system for our simulations.
|
Soil Nutrient Cycling
The submodel for Organic Matter and Nitrogen cycling (OMNI) is linked to other related submodels in the RZWQM such as soil chemistry, solute transport, and plant growth. Significant use of concepts and principles found in nutrient models such as NTRM (Shaffer and Larson, 1987), Phoenix (Juma and McGill, 1986), CENTURY (Parton et al., 1983), and Frissel's N model (Frissel and van Veen, 1981) were also used (Shaffer et al., 2000). Organic Matter and Nitrogen cycling (OMNI) accounts for all N and C processes and pools, with a subset of these processes modeled independently by rate equations. The remaining processes are modeled as functions of specified zero-order and first-order rate equations. The user may adjust many of these rates; however, the model documentation recommends against adjustments of these rates without carefully considering the complexity of the process as implemented in the RZWQM (Shaffer et al., 2000).
The initial dry mass of surface crop residue is user-specified. The model determines the mass incorporated into the surface soil residue pools for initializing the nutrient chemistry model. Initialization of microbial and humus pools will determine how most C and N cycling processes function during the first several years of a simulation. During the simulation, flat surface residue is made available for decomposition after incorporation by the specified tillage operation in CT systems. Standing dead residue becomes flat residue using an exponential decay function after the previous harvest. Nitrifying bacteria are assumed to have full access to NH4 ions (adsorbed + solution). The concentration of NO3– increases at the rate of nitrification minus the assimilation rate of NH4–N for microbial biomass production. The model does not contain a soil anion exchange process, and transport of chemicals under saturated conditions is simulated as piston flow in the mesopore regions of the soil matrix.
Crop Growth
The RZWQM has a single generic plant growth submodel that can be parameterized to simulate different crops. One can choose any of the crops that have already been parameterized for their simulation. Maize, soybean, and winter wheat crops have already been parameterized for the Management Systems Evaluation Areas (MSEA) sites in the midwestern USA (Hanson, 2000). The RZWQM also provides a second option submodel for simulation of crop growth referred to as the Quikplant model. It is a simple grass model that requires inputs such as maximum leaf area index and rooting depth of the crop, total seasonal N uptake, and harvest date. The plant reaches peak LAI, height, and maximum N use in the middle of the growing season. The Quikplant model includes the root input distribution supplied by the user for extraction of water and N from the soil. However, Quikplant is not a detailed growth model and should only be used to simulate water and soil N extraction, and when simulating crop production is not the primary aim of the modeler (Ahuja et al., 2000).
The RZWQM model calculates potential transpiration and soil evaporation using the extended Shuttleworth and Wallace (S-W) model (Farahani and Ahuja, 1996). The extended S-W model includes the effect of surface residue on soil evaporation and partitions evaporation into the bare soil and residue-covered fractions. Actual rates of soil evaporation and canopy transpiration are controlled by the soil water transport and crop growth components of the model (Farahani and DeCoursey, 2000). Water uptake by the roots is evaluated using the approach of Nimah and Hanks (1973), and the equation is solved iteratively by varying the effective root water pressure head until the potential transpiration demand is met based on the ability of the soil to supply the demand. The sum of the sink term cannot exceed the potential transpiration demand. The pressure head reaches a minimum value where it is held steady, and the sum of the sink term for root water uptake from all soil layers then resides below potential demand. The sink term for the Richards equation consists of both the distributed sink due to root uptake, and a point sink arising from tile drainage.
Nitrogen is passively taken up by the plant in proportion to plant transpiration and in quantities necessary to satisfy the present N demand. The amount of N that passively enters the plant is determined by the concentration of N in soil water extracted by the root system from each soil layer. If inadequate N is brought into the plant via transpiration, active N uptake occurs in a manner similar to the Michaelis-Menten substrate model. The total amount of additional N available to the plant through uptake is the sum of passive and active uptake. Available N is hierarchically allocated to roots and then to the other plant organs. Any N remaining after plant demands are met is placed into a storage pool and subtracted from plant N demand the following day (Hanson, 2000).
Model Calibration
General Procedure
After entering all of the model inputs and parameters, we ran the model for a period of 12 yr (3 yr of climate and rainfall data repeated four times) to initialize the organic N pools (rapid, medium, and slow decomposition pools) as suggested in the model documentation. This step is performed before actual calibration begins. The only parameters that we adjusted after the initialization procedure were the soil nitrate and soil ammonium nitrate values for each soil layer beginning on 1 Jan. 1991, the first day of the period for simulating tile drainage and leached nitrate for the analysis period (Table 2). The reason for this was that after the initialization procedure, the model greatly over- or underestimated these values although we had used a value of 1 Mg ha–1 and a wheat cover factor type based on model references and conventional till management practices during parameterization before running the initialization procedure. The measured mineral soil N data had been collected immediately after winter rye was planted for the first time as a cover crop at the study site in the fall of 1991 (Johnson et al., 1999; McCracken et al., 1995); therefore, the measurements reflected the previous 2 yr of winter and spring fallow conditions followed by a maize crop in the summer of 1991. Including a winter rye cover crop as part of the management practices during the 12-yr initialization procedure created more residue for simulated decomposition and, therefore, more mineralized soil N than that measured in the Fall of 1991. However, the simulation period for calibration that began after initialization of the model (1 Jan. 1991–14 Apr. 1993) included conventional till and winter rye cover crop management practices. Resetting the initial values of soil mineral N to their measured values after the initialization procedure before we began the simulations of tile drainage and leached nitrate on 1 Jan. 1991 reflected the soil N conditions at the study site just before the introduction of the winter rye cover crop in the fall of 1991.
|
|
 |
We first calibrated the model without the macroporosity option, and then with macroporosity because measurements of macroporosity were available from the study site (Gupte et al., 1996). We ran the model with the macroporosity option to determine whether or not the model could simulate tile drainage more accurately with macropores since work by Gupte et al. (1996) showed preferential flow in Cecil and related soils of the Piedmont that was not necessarily associated with distinct open macropores but also occurred by way of infilled macropores. We followed the same general procedure for calibration with and without the macroporosity option, and compared the results of simulated tile drainage and leached nitrate with and without macroporosity.
Water Balance Calibration for Tile Drainage
No Macroporosity
To calibrate the water balance, we chose to adjust the ground water leakage rate (water table leakage rate), vw, or the water that will flow out of the bottom of the user-designated soil profile. We used this parameter for calibration because there were no available measurements of it from the study area. We chose a parameter that had not been measured because one of the strengths of our study was the number of field measurements that were available. According to Corwin et al. (1999), the definition of calibration is a test of a model with known input and output information that is used to adjust or estimate parameters for which there is no measured data. We increased the ground water leakage rate beginning with a value of 0 cm h–1 until total simulated tile drainage was within the prescribed 15% range of total observed tile drainage. During this step, we also observed the effect this adjustment had on leached nitrate since chemicals in the soil move with the soil solution. In addition, this assured that simulation of leached nitrate stayed within our target error range of ±15% of mean measured leached nitrate. The measurement period used for comparison after this adjustment was November 1992 through April 1993, after all conventionally tilled plots were in winter rye cover for 1 yr and all of the drain tiles had been installed for at least 2 yr.
With Macroporosity
After we calibrated the model without macroporosity, we calibrated the model with macroporosity by first testing three adjustable macroporosity parameters in the model. The parameters for adjusting the amount of macropore flow that occurs in the soil include the wetting thickness or effective lateral infiltration into the macropore wall (WT), the tile drain express fraction (EF), or the proportion of macropore water that flows to the tile drains, and the sorptivity factor for lateral infiltration (LAB), an adjustment to the calculated Green-Ampt lateral infiltration rate. These parameters were chosen because there was no measured data available for predetermination of possible values, and preliminary runs of the model that showed tile drainage was sensitive to them.
One of the most common forms of sensitivity analysis is to vary model parameters around their base values by some fixed percentage (Silberbush and Barber, 1983; Ma et al., 2000). We chose values of each of the three macroporosity parameters based on the range of values allowed by the model and created a matrix parameter set varying each parameter by approximately 50%. In the case of EF and LAB, initial and final values were increased or decreased from the 50% target value to avoid unreasonable combinations of parameter values. For example, a wetting thickness of zero and an absorption rate of zero with an express fraction of 0.09 would result in 9% of macropore water flowing into the tile drains. However, there would be no absorption into the macropore wall. The parameter set that we used consisted of values of WT ranging from 0.5 cm to 2.0 cm by 0.5 cm; EF values of 0.01, 0.05, and 0.09; and LAB values of 0.1, 0.5, and 1.0, which would reduce lateral absorption calculated with Green-Ampt by either 10, 50, or 0%, respectively, for a total of 36 simulations. The results of each parameter set on total simulated tile drainage and leached nitrate were compared to find the best combination for reducing errors between simulated and measured tile drainage and leached nitrate. Our target error rate of ±15% or less was used for differences between total simulated and total measured tile drainage and total simulated and total measured leached nitrate. We tested each macroporosity parameter and parameter combination set for its sum of squares contribution to the model sum of squares, described below in the model evaluation section, to determine if parameter values needed to be adjusted to a more narrow range of values. Final adjustments of these parameters to best simulate tile drainage for our study in conjunction with crop development provided a better understanding of how macroporosity functions and influences drainage in Cecil soils under conditions modeled, e.g., conventional tillage in maize or cotton production.
Leached Nitrate Calibration
After total simulated and measured drainage were within the 15% error range, we adjusted the sensitive plant parameters in an attempt to bring the simulated aboveground biomass N of the maize crop within, or as close as possible to, 15% of the measured value. We then evaluated the simulated N balance relative to N mineralization to begin refining the calibration for leached nitrate in drainage if needed. In plots with tile drains, Groffman (1984) found that tile drainage in Cecil soils increased aeration and thereby increased mineralization while decreasing gaseous N losses, resulting in a greater supply of nitrate in the drains. Based on available measured data, we evaluated simulated net mineralization for the period from 6 Nov. 1991 through 13 Apr. 1993 as follows:
 |
Crop Growth Calibration
Since plant production was part of the N balance and tightly coupled to the other processes, we followed the procedure for calibrating plant growth recommended for the model by Hanson (2000) when using the generic plant growth submodel. This procedure is based on adjustments to five sensitive plant parameters including active N uptake rate (µ1), daily respiration as a function of photosynthesis (
), the biomass to leaf area conversion coefficient (CLA), and the age effect for plants during the propagule stage and the seed development stage (Ap and As). We used the generic plant growth submodel for both the maize and cotton calibrations, and based adjustments of these parameters for maize within the range of values used for calibration of the MSEA sites (Hanson, 2000). The calibration for cotton development included adjustments of these same sensitive parameters as well as changes to some of the physiological and phenological parameters described below and used in the plant production input file. The calibration for each crop then proceeded by varying each of the sensitive parameters until total biomass and yield were within the 15% range of measured values. During adjustment of these parameters to improve yield simulations to reflect the observed values, we also checked the effect on simulated tile drainage and leached nitrate. This process was used iteratively as depicted in Fig. 2 until simulated tile drainage, leached nitrate, and maize yield were within, or as close as possible, to the desired 15% error range of observed values.
The parameters for the Quikplant model to simulate the winter rye cover crop were obtained from local crop measurements or estimates based on measurements of rye crops (Univ. of Georgia College of Agric. and Environ. Sciences, 1998; Blount et al., 2000). The parameters included total seasonal N uptake, length of growing season (days), maximum crop height, leaf area index, rooting depth, stover after harvest, the C/N ratio of fodder material, and winter dormancy recovery day of year (Appendix A).
After the model was calibrated for tile drainage and leached nitrate in maize and winter rye production, we held all parameters constant and added cotton to the generic plant growth submodel. Parameters were obtained from the cotton field study conducted adjacent to the water quality site (Schomberg and Endale, 2004), and from literature values (Carns and Mauney, 1968; Miley and Oosterhuis, 1990; Univ. of Georgia College of Agric. and Environ. Sciences, 2000; Nyakatawa and Reddy, 2000; Nyakatawa et al., 2000; Reddy et al., 2004). The cotton calibration simulation period was 1 Jan. 1997 through 31 Dec. 1997. The parameters adjusted for cotton in the generic plant growth input file included the physical dimensions of the plant; the maximum, minimum and optimum temperature for growth; maximum leaf area index; and the minimum number of days the plant required to transect each physiological growth stage (Appendix A7). Through iterative adjustments of these parameters, we compared simulated and observed cotton total biomass until simulated values were within 15% of observed values. Since we did not have measures of tile drainage or leached nitrate from the study used to calibrate for cotton, we compared simulated and calculated PET and water use to ensure that the model could simulate both reasonably well as a basis for later evaluating simulated tile drainage in the water quality study. Calculated PET was based on the Priestley-Taylor equation (Priestley and Taylor, 1972) from the weather station near the water quality study site, approximately 100 m from the cotton calibration study site (Hoogenboom, 2003). We calculated observed and simulated water use for cotton as:
 |
Soil Water = Soil wateri+1 – Soil wateri, where soil water was measured or simulated in a 60-cm profile, and i = day of year. Observed soil moisture in cotton showed little or no change below 60 cm in the field study used for calibration (Schomberg and Endale, 2004). The rooting depth for cotton is shallow in acid soils because it is one of the most sensitive crops to aluminum toxicity, which frequently occurs in acid subsoils such as those in Georgia (Mitchell et al., 1991; Sumner, 1994; Gascho and Parker, 2001).
Evaluation of Simulation Results
For the analysis of the macroporosity parameters, we tested for the main effects and interactions of the three parameters for macroporosity and selected the most significant effects based on the Type I sum of squares each contributed to the model sum of squares (SAS Institute, 2000). Based on this information, we identified the parameter or combination of parameters with the highest correlations and highest probabilities for simulated and observed tile drainage. We determined at this point whether further testing was needed within a more narrow range of parameters. Since one of our objectives was to try to simulate how macropore flow may affect drainage in Cecil soils, we chose to refine the range of the parameters as much as possible to improve our understanding by way of the simulation process.
For the analysis of simulated and observed values of tile drainage and of nitrate leaching, we regressed the final observed values on simulated values using linear regression analysis (SAS Institute, 2000), to compare r2 values, slopes, and intercepts. We calculated the relative root mean square error (RRMSE), standard error of the mean difference (Addiscott and Whitmore, 1987), maximum error, average and standard deviation of measured and simulated values to characterize systematic over- or under-prediction, and used graphical displays to show trends and distribution patterns (Loague and Green, 1991). The RRMSE, which is the RMSE relative to the mean of the observed values, is calculated as follows:
 |
For the cotton water use analyses, we compared observed and simulated water use for the period of peak water use in cotton from first square to first bloom and from first bloom to peak bloom. Our criteria for the acceptable differences between daily simulated and observed water use was ±15% or less based on a minimum daily value of 6.4 mm of water use for cotton during these periods, a total of approximately 55 d in July and August (NCSU-CES, 2004).
|
RESULTS AND DISCUSSION |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSRESULTS AND DISCUSSIONCONCLUSIONSAPPENDIXREFERENCES |
|---|
|
The analysis of simulated and observed soil mineral N for each day of 12 field-measured values from November 1991 through April 1993 revealed that 3 of the 12 simulated soil mineral N predictions were outside the 95% confidence interval (CI) of observed soil mineral N. Total simulated tile drainage and leached nitrate for the final analysis period of November 1992 through April 1993 were 6 and 5% of total observed values, respectively. Since we met our objective of obtaining a difference between simulated and observed values for tile drainage, leached nitrate and maize yield of 15% or less for the final analysis period, we considered the calibration acceptable as the final calibrated scenario for maize production with a winter rye cover crop without macroporosity.
The analysis of simulated tile drainage and leached nitrate for the calibrated scenario revealed that cumulative simulated tile drainage followed the pattern of cumulative observed tile drainage (Fig. 4) . However, 7 of 12 simulated drainage events were outside of the 95% CI of observed tile drainage events (Fig. 5) . Simulated leached nitrate increased at the same rate as observed leached nitrate during the first five drainage events and then leveled out at or near zero for the remaining seven events while observed leached nitrate continued to increase slightly (Fig. 4). Six out of 12 simulated leached nitrate events were outside the 95% CI of observed leached nitrate (Fig. 5). The RRMSE showed a large percentage deviation from the mean observed values, reflecting the fact that the majority of simulated events for tile drainage and half of the events for leached nitrate were outside of the 95% CI (Table 3). Linear regression analysis of total observed tile drainage on total simulated tile drainage revealed that simulated tile drainage explained 37% of the variation in observed tile drainage. Analysis of total observed leached nitrate on total simulated leached nitrate revealed that simulated leached nitrate explained 90% of the variation in observed leached nitrate. The slopes of the regressions for both response variables were not significantly different from one, and the intercepts were not significantly different from zero at the 0.05 probability level (Table 4).
|
|
|
|
With these three parameters selected for macroporosity, the model simulated a very large amount of nitrate with large increases in leached nitrate and net mineralization and smaller increases in the other N balance components for the N balance analysis period (Table 5). We tried six other combinations of the macroporosity parameters that also showed high correlations between simulated and observed tile drainage and leached nitrate for the final analysis period (November 1992–April 1993). In each case, simulated net mineralization increased, and the system produced too much nitrate and increased one or more of the N components by large amounts. The N balance became very volatile with the inclusion of the macroporosity option, and we were not able to simulate the N balance components, including net mineralization, to within ±15% of observed values during the N balance analysis period. Adjustments to one or more of the other sensitive plant parameters such as N uptake (µ1) or the proportion of photosynthate to respire () could cause the model to suddenly generate unreasonably high amounts of nitrate in one or more N components such as leached nitrate. We also found that more than one combination of values for the sensitive plant parameters would simulate yield and possibly simulate one other N component such as leached nitrate within 15% of observed values, but again would create large changes in other components of the N balance such as crop N uptake. This would then create a situation that required an endless number of iterative adjustments to bring simulate leached nitrate, tile drainage, and yield back to within a reasonable range of observed values. After several attempts to adjust the macroporosity components and the sensitive plant parameters to simulate leached nitrate and net mineralization within our 15% target error range without success, we set both the nitrification and denitrification rates to zero to allow the model to produce mineral N by way of organic matter decay and microbial biomass N mineralization and decay (Shaffer et al., 2000). Under these conditions, the OMNI submodel will test for sufficient NH4+ and NO3– in the system and shut down the decay process if net immobilization is occurring, limiting the amount of NH4+ that can be released by the microbial biomass decay process. In contrast, nitrifying autotrophic bacteria have full access to NH4+ in the model in both adsorbed and solution phases so that as long as mineralization is occurring, NH4+ will be nitrified. Setting the nitrification and denitrification rates to zero decreased soil nitrate N and increased soil ammonium N. This also reduced leached nitrate, although it was still 48 kg ha–1 greater than observed leached nitrate, and N uptake by the second winter rye crop increased 17 kg ha–1. Finally, we set the nitrification and denitrification values back to the model default values, and increased the denitrification rate incrementally to decrease the amount of nitrate in the system (Table 5). Using a Latin Hypercube Sampling technique to determine the sensitivity of crop N uptake, silage yield, and nitrate leaching below the root zone in the RZWQM, Ma et al. (2000) found that all of the responses were negatively related to the denitrification constant. In addition, the authors found that a combination of mean irrigation and manure application rates simulated leached nitrate concentrations from 0 to 755 kg N ha–1. They described the outcome of combining irrigation and manure rates as the worst scenario for simulating their response variables. By using the model default nitrification rate and increasing the denitrification rate, we were able to stabilize the N balance components, and to simulate leached nitrate, tile drainage, and maize yield more accurately for the final analysis period of November 1992 through April 1993. However, we were not able to simulate any of the response variables to within ±15% of observed values during the N balance analysis period (Table 5).
|
The analysis of simulated tile drainage and leached nitrate with macroporosity for the calibrated scenario revealed that cumulative simulated tile drainage followed the pattern of cumulative observed tile drainage (Fig. 4) with 8 of 12 simulated drainage events outside the 95% CI of observed tile drainage events (Fig. 5). There were no significant differences in the means or the variances between tile drainage simulated with or without macroporosity. Simulated leached nitrate increased at the same rate as observed leached nitrate during the first 5 of 12 drainage events following the same pattern as simulated leached nitrate without macroporosity (Fig. 4). Six of the 12 simulated leached nitrate events were outside the 95% CI of observed leached nitrate, as was the case with no macroporosity (Fig. 5). The RRMSE showed a large percentage deviation from the mean observed values, reflecting the fact that the majority of simulated events for tile drainage and half of the simulated leached nitrate events were outside of the 95% CI of measured events (Table 3). Linear regression analysis of total observed tile drainage on total simulated tile drainage, and total observed leached nitrate on total simulated leached nitrate revealed nearly the same relationship as the regressions without macroporosity. The slopes were not significantly different from one, and the intercepts were not significantly different from zero at the 0.05 probability level (Table 4). There were no significant differences between the means or the variances with and without macroporosity for simulated leached nitrate.
Though it was more difficult to calibrate the model with macroporosity than without macroporosity due to the volatile nature of the N balance with macroporosity, the differences between simulated tile drainage and leached nitrate relative to macroporosity indicated that macroporosity did not have a significant influence on the amount of tile drainage that occurred in these soils. In a study of intact dye-stained soil cores from the study area in conventional tillage, Gupte et al. (1996) found little evidence of preferential flow in the upper 45 cm of the cores. Preferential flow often occurred in regions of soil and in-filled macropores at depths between 45 and 60 cm rather than through open macropores. In addition, the presence of tile drains below 60 cm in our study would influence the way drainage occurs both in the field and in model simulations due to the difference in the flow patterns created when tile drains are present (Skaggs, 1980; Ritzema, 1994). Any preferential flow that occurs due to the presence of macropores near the depth of the tile drains would be difficult to quantify separately from the impact of tile drains on soil water flow.
The contribution of macropore flow to simulations of nitrate leaching was also difficult to quantify because the amount of nitrate leached was greatly affected by changes to other parameters such as the plant parameters, the nitrification and denitrification rates, and the macroporosity parameters. This was in spite of the fact that we narrowed the combination of adjustable parameters for macroporosity to those that best simulated our response variables before adjustments to any of the plant parameters. A sensitivity analysis using all of the combined parameters that appeared to affect nitrate leaching with the macroporosity option might be effective in the case of calibration for nitrate leaching for one scenario or one study. However, based on our experience in this study, it is likely the model will not perform consistently if the conditions tested are different than those under which the model was calibrated due to the volatile nature of the N balance once macroporosity is introduced. Ma et al. (2000) concluded that the interdependency of various parameters can introduce high variability in response variables that are tested with the RZWQM, but that model output responses can be much less sensitive to variations in one parameter than in the other. A closer examination of this variability is needed where the model produces large amounts of nitrate with minor changes to crop parameters or N rates before the model can be expected to perform in a reliable manner in subsequent simulations of nitrate leaching with the macroporosity option.
Cotton Calibration
After cotton was included in the generic plant growth submodel, the sensitive parameters, µ1, Ap, As, , and CLA, and leaf stomatal resistance were iteratively adjusted as well as the minimum number of days required for the vegetative and reproductive growth stages to simulate cotton biomass development as closely as possible to observed development (Table 6). We also adjusted the albedo of a mature plant to 0.2 based on model references to bring total simulated PET at the end of the cotton growth period as close as possible to total calculated PET from the weather station near the study site (Hoogenboom, 2003). The result of this adjustment was a difference of <3 mm between total simulated and total calculated PET for the cotton growth period. Adjustments to the minimum number of days for each of the vegetative and reproductive growth phases were particularly sensitive in our efforts to achieve a growth pattern and values for simulated cotton biomass and leaf area index that matched observed values on measurement days. It was not possible to simulate total biomass to within 15% of total observed biomass despite numerous iterative adjustments and combinations of the phenology parameters. This was because the model could not produce the large increase in observed biomass between Day 245 and Day 266 without adjusting the plant parameters to rapidly increase total biomass early in the season (before the first bloom period for cotton) (Fig. 6)
. When we adjusted the parameters to rapidly accumulate biomass early in the season, after simulated vegetative growth peaked, biomass accumulation would begin to decline as the simulated plant entered the reproductive stage followed by leaf senescence. The optimum balance for the number of days in each of the vegetative and reproductive growth stages to achieve a simulated pattern of development that matched observed development for cotton resulted in a period of 115 d for the vegetative stage and 40 d for the reproductive stage. This allowed the model to simulate cotton biomass accumulation and leaf area similarly to observed biomass and leaf area during the majority of the growing season by slowing biomass accumulation until the last 21 d of observed cotton boll development when simulated biomass began to decline (Fig. 6).
0.2 cm diameter for soil column depths from 
= distance from the water table to the impermeable layer, m = water table height above the drains, d = distance from the drain to the impermeable layer, and L = distance between drains. Design is based on the Hooghoudt steady state equation to estimate the flux at the center of the drains and correct for two-dimensional flow.
