The Soil and Water Assessment Tool (SWAT) model has been successfully used to predict streamflow, evapotranspiration and soil water. The crop growth model in SWAT was adapted from the EPIC model and is similar in concept to the crop growth models in Agricultural Policy/Environmental Extender Model (APEX, ), ALMANAC, and WEPP, which have undergone significant crop yield validation. SWAT crop yields have been validated for several grain crops .
Preliminary data suggest that while the hydrologic balance in each watershed may be accurately simulated with SWAT, the SWAT model tends to over or under predict wheat yield responses to N-fertilizer application. For example, Haney found that simulated wheat yield increased strongly with N-fertilizer additions (r2 = 0.80) and that yields at higher N fertilizer rates were overestimated and at lower N fertilizer rates were underestimated. In addition, when N fertilizer was not applied during simulation, SWAT predicted yields were close to 0 Mg/ha. These results indicate that SWAT may not be properly accounting for soil N cycling processes. Yield underestimates at low fertilization rates could occur if modelled N―mineralization rates are underestimated, causing underprediction of plant N availability and over-estimates of N limitation. Underestimates of plant N availability can compound yield errors by suppressing yield responses to simulated soil water variation. Furthermore, SWAT neglects the contribution of the soil microbial population to the plant-available N pool, resulting in an underestimation of yield and possible over or underestimation of N runoff from natural and agricultural landscapes.
The SWAT model now has three different N simulation options, SWAT-flush , N routines derived from the CENTURY model , and a one-pool C and N model option . SWAT-flush cycles N through three organic N pools (fresh residue, stable and active organic) and two inorganic N ( and ) pools with an added flush of N after significant rainfall events (greater than 26 mm). The variation of the CENTURY model is more complex than SWAT-flush, simulating microbial, slow and passive soil organic N, surface microbial N, above and below ground structural and metabolic N, and mineral N . On the other hand, the one-pool C model merges C, N, and P soil organic matter (SOM) pools within each soil layer, as well as separate residue and manure pools in the topsoil and subsoil.
In this study, wheat yield values from a long-term wheat yield data from fertilizer study research plot in north-central OK, were compared to simulated wheat yield values from SWAT-flush, SWAT-C, and SWAT-One. The objective of this study is to assess the ability of various N cycling sub-routines within SWAT to predict yield at a long-term fertilizer study in Oklahoma.
2. Materials and Methods
SWAT-flush utilizes three organic N pools (fresh, stable and active organic) and two inorganic N pools ( and ) and an added flush of after rainfall events greater than 26 mm (Figure 1). The SWAT-flush model algorithms were derived from the PAPRAN (Production of Arid Pastures limited by Rainfall and Nitrogen) model . Mineralization, decay, and immobilization equations are first order kinetics, which are based on the substrate amount, determined by a model “warm up” period of several years prior to the years of interest.
The sizes of the organic N pools are assigned assuming that the C:N ratio for humic materials is 14:1. The concentration of humic organic nitrogen is determined based on the soil organic C (SOC) values from soil data contained in SWAT input files. The soil data must be entered by the user and can either be obtained from soil sampling or publicly available data sets such as the Soil Survey Geographic Data Base (SSURGO). SWAT-flush then assigns 20% of the organic N to the active pool and 80% of the organic N to the stable pool . The initial residue (fresh) pool is assigned to the top 10 mm of the soil profile and is set to 15% of the initial amount of residue on the soil surface, and does not include root biomass. After initialization, the fresh pool is determined based on simulated management practices. The simulated N resulting from decomposition and mineralization of the fresh pool is partitioned as 20% to the active organic and 80% to the pool. Decomposition and mineralization in SWAT-flush depend on the residue decomposition rate, the C:N and C:P ratios of the residue in the soil layer, and soil temperature and water content. N cycling processes are calculated for each soil layer within the profile.
Initial concentration is an exponential function of soil depth. The pool is initially set to zero and only contributes to the NO3 pool when urea fertilizer is added to the soil. Nitrification and volatilization describe the conversion of to either or NH3, respectively. SWAT-flush simulates both processes simultaneously then partitions the calculated values between the two processes . The nitrification process in SWAT-flush depends solely on the soil water and temperature. While temperature and soil moisture are critical forcing factors on the nitrification process, SWAT-flush does not specifically account for soil microbial activity, soil pH, or the water-extractable soil C or N content, which form the C and N source for the microbial population. Volatilization simulation in SWAT-flush depends on soil temperature and depth and includes a default cation exchange factor. Volatilization is also strongly affected by soil pH, wind conditions, and soil clay content and type .
SWAT-flush incorporates an addition of to the pool after a rainfall event was based on the water soluble organic C and N (WSOC and WSON) and microbial activity determined using 1 − d CO2 evolution (Equation (1)).
The flush of N is added to the top 10 mm of soil to simulate rapid changes in soil moisture, temperature and N cycling at the soil surface. After a significant rainfall event (greater than 26 mm) occurs on sufficiently dry soil (based on soil matric potential), a flush of N is added to the pool.
Figure 1. The N cycle as defined in the SWAT-flush model .
The CENTURY-based N simulation option was first incorporated into the Environmental Policy Integrated Climate (EPIC) model , and was then incorporated into SWAT for testing at the watershed scale and referred to hereafter as SWAT-C . The CENTURY option is a multi-pool model whose strength lies in the linkage between organic C and N dynamics. The CENTURY option in SWAT (SWAT-C) includes a residue pool consisting of lignin, non-lignin, and metabolic residue, each having its own decomposition rates (Figure 2). Residue dynamics occur at the surface of the soil and in the top 10-mm layer of soil. SOM is simulated as microbial, slow, and passive pools, each with their own turnover rates. The microbial pool occurs in all soil layers, while the slow and passive pools exist in all soil layers except the top 10 mm. Decomposition of residue and mineralization of SOM depends upon lignin content of the residue, soil temperature, texture and moisture, tillage effects, and O2 content. Depth profiles of O2 in SWAT-C differ from both those of CENTURY model and SWAT-flush. Residue composition and lignin content are calculated based on plant age. All mineralization and decomposition processes result in the simultaneous transformation of C and N and ultimate release of CO2 .
As with other models where N dynamics are based on first order kinetics
Figure 2. Carbon and N cycling in the SWAT-C subroutine in SWAT (Recreated from Zhang, X., R. C. Izaurralde, J. G. Arnold, J. R. Williams, and R. Srinivasan, Modifying the Soil and Water Assessment Tool to simulate cropland carbon flux: Model development and initial evaluation, Page 812, 2013, with permission of Elsevier).
(basic SWAT), C and N flows in the SWAT-C are controlled by the size of the pools. It is therefore critical that the various organic pools are initialized and tracked correctly. It has been reported the CENTURY model successfully simulates daily CO2 fluxes except during rewetting periods, which is when important N mineralization fluxes occur . SWAT-C was tested by Zhang et al. by comparing simulated results to corn and soybean crop yields on lands across the U.S. Midwest. They found that SWAT-C performed well in its simulations of annual crop yield for sites where detailed management data was known. On the site where management data was not available or sparse, model performance was reduced. In general, Zhang et al. found that long-term average crop production (corn and soybean) was predicted well using SWAT-C.
The third N cycling option is SWAT-One, a one-pool C, N, and P model . SWAT-One simulates decomposition of a lumped C, N, and P soil organic matter (SOM) pool within each soil layer, as well as residue and manure pools in the topsoil and subsoil (Figure 3). Decomposition of residue and manure follows first order kinetics and results in either mineralization or immobilization depending upon the humification rate and C:N and C:P ratios of the residue, the manure and the SOM. Manure and residue C is either incorporated into the soil C pool or respired as CO2, and their decomposition rates are functions of soil temperature and moisture. Maximum formation of humus from residue is 0.18 g g−1, and the manure maximum humification rate is 1.6 times higher. N mineralized is transferred to the soil pool. The C:N and C:P ratios of newly formed SOM vary throughout simulation depending upon available mineral N and residue or manure C:N ratios. If there is not enough organic N to supply the microbial N needed for decomposition with a continuously changing soil C:N ratio, mineral N is immobilized. SOM decomposes depending upon a tillage factor and soil moisture. Mineralized N from the SOM is transferred to the pool and is always positive. Testing of the SWAT-One option has been minimal to date.
The various SWAT simulations were compared to data obtained from Oklahoma State University’s long-term wheat yield study (Experiment 502) in Lahoma, OK. The Experiment 502 plot research is conducted
(http://nue.okstate.edu/Long_Term_Experiments/E502.htm) at the North Central Agricultural Research Station near Lahoma, OK (36.42˚N, 97.87˚W) in Garfield county . The OSU study was established in 1970 to study the response of wheat grain yield to varying rates of long-term N, P, and K fertilizer application. The randomized complete block (4 replications) designed experiment is conducted on continuous winter wheat grown under conventional tillage on a Grant silt loam (fine-silty, mixed, thermic Udic Argiustoll). The soil has an average pH of 5.7 in the top 30 cm . Soil depth, texture, slope, albedo and SOC content were obtained from SSURGO . Mean average temperature at the research site is 15.6˚C with an average annual rainfall of approximately 800 mm . Nitrogen was applied as Urea (46-0-0) at pre-plant rates of 0, 22.4, 44.8, 67.3, 89.7, 112.1 kg N/ha annually. Phosphorus was applied as triple superphosphate
Figure 3. The one-pool C, N, and P submodel (SWAT-One) within SWAT (Reprinted from Kemanian, A. R., S. Julich, V. S. Manoranjan, and J. R. Arnold, Integrating soil carbon cycling with that of nitrogen and phosphorus in the watershed model SWAT: Theory and model testing, Page 1915, 2011, with permission of Elsevier).
(0-46-0) at the rates of 9.9, 19.7, 29.6, 39.5 kg P/ha annually. Fertilizer application occurred between early August and early October and planting followed between late September to late October. The wheat seeding rate varied between 0.07 and 0.08 Mg/ha. Grain was harvested from early June to early July, depending upon weather conditions.
Yield simulations were performed by constructing a set of SWAT input files using local weather and soils data in the Texas Best Management Practice Evaluation Tool (TBET, ). Weather data were obtained from the National Oceanic and Atmospheric Administration (NOAA) National Weather Service Cooperative Observer Program, Lahoma Research Station (USD00344950) weather station (Latitude: 36.3894, Longitude: −98.1061, Elevation: 388.6 m). Simulations were run for 28 years from 1981 to 2012, for which yield data and most actual planting and harvest dates were available. 1985 and 1986 served as warm up years, to allow initial fractions of SOC and other variables to stabilize prior to simulation of the period of study. When dates were unavailable, an average October 21 planting date and June 13 harvest date was used (6 cases where one was missing). Simulations included 12 combinations of N and P rates and forms that correspond with the fertilizer rates used in Experiment 502.
Simulations were performed with uncalibrated SWAT models. Previous research has indicated that the SWAT model can successfully predict crop yield without calibration . In addition, we were interested in seeing the raw results from an uncalibrated model for comparison to actual field data. Yield data obtained from each N modeling option in SWAT (SWAT-flush, SWAT-C, and SWAT-One) were compared to historical yield data using linear regression analysis, descriptive statistical analyses, percent bias (PBIAS), Nash-Sutcliffe efficiency (NSE) and nitrogen use efficiency (NUE) analysis, Pearson Correlation Coefficients and Analysis of Variance (ANOVA, ).
Nitrogen Use Efficiency was calculated by taking the average yield at the 22, 45, 67, 90 and 112 kg N/ha fertilizer application rates, subtracting the control (0 kg N/ha) yield and dividing by the fertilizer application rate. NUE is chiefly regulated in all SWAT model N variations using attributes listed in the plant growth database (crop.dat). This database includes plant classification (i.e. warm season annual), radiation-use efficiency, harvest index, maximum potential leaf area index (LAI), optimal and base temperature for plant growth, maximum rooting depth and canopy height, the fraction of N in the harvested portion of the biomass, and potential heat unit information at various stages of growth . The optimal N that should be in plant biomass on a given day is calculated by first determining the fraction of N in the plant as a function of growth stage under optimal growing conditions. Specifically, the fraction of N in a plant on a given day is determined based on the fraction of heat units accumulated on that day and the fraction of N at emergence, maturity, and midseason which were determined experimentally for winter wheat by The University of Saskatchewan . Optimal biomass N for the day is the product of the fraction of N in the plant on a given day and the biomass on the same day:
where bioN,opt is the optimal mass of nitrogen stored in plant material for the current growth stage (kg N/ha), frN is the optimal fraction of nitrogen in the plant biomass for the current growth stage, and bio is the total plant biomass on a given day (kg・ha−1). Potential N uptake is subsequently determined using the following equation:
where Nup is the potential nitrogen uptake (kg N/ha), bioN,opt is the optimal mass of nitrogen stored in plant material for the current growth stage (kg N/ha), bioN is the actual mass of nitrogen stored in plant material (kg N/ha), frN,3 is the normal fraction of nitrogen in the plant biomass at maturity, and Dbio is the potential increase in total plant biomass on a given day . Daily control of N uptake depends upon biomass growth each day and the amount of available N in the soil. The amount of available N is determined by initial N in the soil, fertilizer applications, leaching and surface N runoff.
The Nash-Sutcliffe efficiency (NSE) was calculated for simulated versus actual yields for each year, averaged over all fertilizer treatments, to determine how well the observed values versus simulated values fit the 1:1 regression line. NSE is calculated as follows:
where is the simulated yield value at time t, yt is the actual yield at time t, and is the mean of the observed data values for the entire evaluation time period. NSE values range from −∞ to 1 and the larger the NSE values, the better the model performance . Percent bias (PBIAS) was used to statistically measure the average propensity of simulated data to be larger or smaller than observed values . Percent bias is calculated as:
where ft is the simulated yield value at time t, and yt is the actual yield at time t. Smaller PBIAS values are desired. Negative PBIAS values indicate model under-estimation, while positive values indicate model over-estimation bias . PBIAS values less than 15% are considered acceptable.
3. Results and Discussion
Actual and simulated yields averaged over 28 years were positively correlated with N fertilizer additions (r2 > 0.96) except SWAT-One (r2 = 0.33, Figure 4). When no fertilizer was applied, actual wheat yields averaged 1.71 Mg/ha. SWAT-C average simulated yield was 1.97 Mg/ha at 0 kg N/ha applied and was closest among the submodels to simulating actual yield at this fertilizer rate (PBIAS, 13%). SWAT-C most closely simulated the effect of fertilizer on average actual wheat yield over a 28-year period according to regression and PBIAS analysis (PBIAS, 2%; Table 1). SWAT flush had an improved average NSE value (NSE, −0.05) over SWAT-C (NSE, −0.52). SWAT-One overestimated yield across simulated fertilizer application (PBIAS, 61%). SWAT-flush underestimated yield below the 67.3 kg N/ha fertilizer treatment, but over estimated at higher fertilization levels (PBIAS −19% at 0 kg N/ha and 9% at 112.1 kg N/ha). Srinivasan et al. found that PBIAS values of simulated yield varied from region to region depending upon the soil data used by SWAT to simulate soil processes. The data in this study; however, indicate that the PBIAS of simulated yield can also vary drastically depending upon the way that N cycling is treated in the model. The yield underestimates from SWAT-flush at lower levels of N fertilization suggest that his submodel underestimates plant N availability (or over-estimates N losses) and would lead to overestimates of the N-fertilizer inputs needed to achieve a given yield.
Yield under- or over-estimates result in erroneous estimates of percent Nitrogen Use Efficiency (NUE). NUE ranges from 23% - 50% in winter wheat cropping systems . Oklahoma State University reports that NUE for Experiment 502 averages 32% , although we calculated the average actual NUE values at less than 25%, decreasing with increasing fertilizer applications to 14% at 112.1 kg N/ha treatment. Simulated NUE values for the submodels also decreased with increasing amounts of N fertilizer (25% to 20% for SWAT-flush, 15% to 14% for SWAT C, and all negative NUE values for SWAT One). The percent of N removed in grain relative to N uptake was 66% for SWAT-C and 70% and SWAT-flush. This value was 37% for SWAT-One, partially explaining the negative NUE values for SWAT-One. SWAT-flush most accurately represented field NUE values. Nitrogen use efficiency can be an important indicator of N dynamics in the soil and is reflective of nitrification, management, weather and plant growth . Factors affecting low NUE in the field include losses of N from volatilization, which can be as great as 50% when urea or urea-containing products are applied . In addition, N runoff losses range between 1% and 13% . Certainly, biomass growth, pest, weed, temperature, and moisture stress also affects NUE in the field.
Deviations between simulated and actual yields and NUE values can be partly explained by variations in the way each model handled specific nitrogen pools and transformations (Figure 5). Volatilization varied among the sub-models. On average 81% of fertilizer applied was lost to volatilization with SWAT-One versus a 37% loss with SWAT-flush and SWAT-C. All sub-routines utilize the same volatilization/nitrification subroutine; however, SWAT-One may have simulated higher volatilization compared to SWAT-flush because all N mineralized is added to the pool, instead of the pool. SWAT-flush does not simulate volatilization unless an based fertilizer is used. We expected SWAT-C to have greater volatilization, but the average values were similar to those from SWAT-flush.
Table 1. Percent bias (PBIAS) Nash-Sutcliffe efficiency (NSE) and values of model simulated values at 0 kg N/ha, 112 kg N/ha, and the average of all fertilizer treatments.
Figure 4. Relationship between fertilizer additions and simulated and actual yield averaged over 28 years.
Figure 5. Simulated N cycling values for each of the three N-cycling subroutines.
Simulated nitrification also varied among the N sub-routines. SWAT-One simulated the greatest amount of nitrification, followed by SWAT-C and SWAT-flush. The fact that SWAT does not account for pH of the soil (average pH 5.7) in any N subroutine complicates replication of the N cycle and plant growth operations, which could have an effect on yield prediction capabilities. Nitrification at pH 5.7 should be low as rates fall distinctly below pH 6 and nitrification should be almost non-existent below pH 5.0 . SWAT-flush nitrification values were equal to that of the fertilizer added minus the volatilized. SWAT-one and SWAT-C will continually have significantly higher volatilization and nitrification values than SWAT-flush because they simulate the transformation of organic N to . It may be beneficial to add a pH control to these processes because the volatilization and nitrification values were unrealistically high given the actual pH of the soil simulated. Furthermore, SWAT-flush may benefit from converting mineralized N into versus to more realistically simulate field N transformations.
Denitrification did not occur in any of the simulations, and therefore did not contribute to simulated yield estimates. Denitrification occurs in the absence of O2, and varies depending upon soil moisture, temperature, organic matter content, C and concentration . SWAT usually only simulates denitrification under flooded conditions, although it is well documented that this process occurs in small pockets of the soil profile where anaerobic conditions can take place, regardless of the level of the soil water table or complete saturation of the soil.
Simulation of NH4-N or NO3-N pools is also critical to accurate yield estimates. N fertilizer applied in excess of plant uptake should increase soil N pools . During simulations, total soil N (organic N + -N) was only increased when yield was simulated using SWAT-One. SWAT-One simulated a marked increase in organic N and -N in the soil. SWAT-C simulated an overall decrease in organic N. -N values decreased when using SWAT-flush, but increased when using SWAT-One, and SWAT-C. Because this is a conventional till, wheat-fallow cropping system, we would expect that organic N values would decrease over time as a result of long-term losses in SOC . Because of the soil texture, we would expect values in the soil to decrease on average. Surprisingly, excess N was not lost to leaching using any of the sub-models even though the soil is a silt-loam and should drain well. Based on these data, it appears that none of the N subroutines are adequately simulating N cycling processes in the soil.
There was a significant correlation between SWAT-flush (r2 = 0.34, p < 0.001) and SWAT-C (r2 = 0.20, p < 0.001) predicted and actual annual yields (Figure 6). Although the trend is significant, the variability around the regression line indicates that neither model is precise in its predication of annual wheat yield. SWAT-One annual predicted yields were not correlated with actual annual yield. These data suggest that the N cycling models may be ineffective at simulating mineralization, decomposition or the conversion of urea.
As it is in nature, plant growth is moderated in the SWAT model due to water, nutrient, and temperature stress. SWAT calculates the amount of stress for water, temperature, N and P stress on a daily basis and reduces plant growth as a percentage of optimal growth when the plant is not dormant. Potential biomass production for each day is calculated as the potential increase in total plant biomass on a given day multiplied by the plant growth factor :
where greg is the plant growth factor (0.0 - 1.0), wstrs is the water stress for a given day, tstrs is the temperature stress for a given day expressed as a fraction of optimal plant growth, nstrs is the nitrogen stress for a given day, and pstrs is the phosphorus stress for a given day. Potential leaf area added on a given day is also adjusted daily for plant stress in the same manner.
All three N cycling options utilize the one maximum stressor for each day. For example, if temperature stress is at 30% and water stress is at 20% for the day, the 30% temperature stress is used to regulate plant growth on that day. Annually, the wheat plants in all simulations were under some kind of stress (below optimal conditions for plant growth) between 126 and 177 days per year (Table 2). Phosphorus stress (not shown) was negligible and therefore not reported. Overall, SWAT-One had the least amount of stress days (especially N), which corresponds to its consistent over-prediction of yield. SWAT-flush had the highest N stress days, most likely due to under-prediction of N mineralization at
Figure 6. Regression between yearly SWAT-flush and SWAT-C predicted and actual yield for all fertilizer treatments.
Table 2. Average annual water, temperature, N, and P stress days for each N cycling routine for the 28-year run.
low N fertilization rates. Temperature stress was the same throughout the three N submodels. It appears that water stress was low for all N subroutines and was overshadowed by temperature or N stress on most days.
These results indicate that, except under extreme wet or dry conditions, temperature and N have a stronger influence over simulated yield on a yearly basis than soil moisture. Simulated yields did not differ significantly regardless of yearly precipitation, which indicates that rainfall is not the significant controlling factor in predicted yield. In fact, research has shown that annual precipitation, spring precipitation, and growing season precipitation were not directly correlated to predicted or actual crop behavior (unpublished data). Lobell et al. found that in the field soil variability is greater than variability in weather when water availability is not a limiting factor. Based on the results shown in Table 2, it appears that the SWAT model may be overly sensitive to temperature stress, thereby reducing the importance of both N and water stress.
We found that although multi-year average simulated crop yields were well correlated with actual average yields, SWAT-flush underestimates yield at low N fertilizer levels then over-estimates at higher N fertilization. SWAT-flush most accurately represented field NUE values. On average, SWAT-One was unsuccessful at predicting yields. SWAT-C most closely estimates average yield according to calculated PBIAS values, while NSE calculations indicate that SWAT-flush is more capable of predicting average yield. The N removed in yield relative to N uptake and N volatilization were surprisingly similar between SWAT-C and SWAT-flush; however, nitrification, final in soil, and the amount of water and N stress varied between the two models. Annually, SWAT-flush and SWAT-C yields were correlated with actual yield, but showed a high degree of variability indicating that these submodels may not be reliable for predicting annual yield.
Overall, this research indicates that SWAT-C or SWAT-flush provides the most accurate prediction of average wheat yield and can be used for wheat cropland yield assessment. However, none of the N-cycling routines included in the SWAT model predict annual variations in wheat yield with great certainty. Further research is needed to determine the effectiveness of SWAT-C and SWAT-flush in determining average and annual yield in various farming regions and with numerous agronomic crops.
This research was funded with a cooperative agreement funded by USDA-ARS with NRCS agreement number 60-3098-5-006.
 Gassman, P.W., Williams, J.R., Wang, X., Saleh, A., Osei, E., Hauch, L.M., Izaurralde, R.C. and Flowers, J.D. (2010) The Agricultural Policy/Environmental eXtender (APEX) Model: An Emerging Tool for Landscape and Watershed Environmental Analyses. Transactions of the ASABE, 53, 711-740. https://doi.org/10.13031/2013.30078
 Haney, E.B., Haney, R.L., Arnold, J.G., White, M.J., Srinivasan, R. and Senseman, S.A. (2016) Spatial Analysis and Modeling the Nitrogen Flush after Rainfall Events at the Field Scale in SWAT. American Journal of Environmental Sciences, 12, 102-121.
 Zhang, X., Izaurralde, R.C., Arnold, J.G., Williams, J.R. and Srinivasan, R. (2013) Modifying the Soil and Water Assessment Tool to Simulate Cropland Carbon Flux: Model Development and Initial Evaluation. Science of the Total Environment, 463-464, 810-822.
 Kemanian, A.R., Julich, S., Manoranjan, V.S. and Arnold, J.R. (2011) Integrating Soil Carbon Cycling with that of Nitrogen and Phosphorus in the Watershed Model SWAT: Theory and Model Testing. Ecological Modelling, 222, 1913-1921.
 Seligman, N.G. and van Keulen, H. (1981) PAPRAN: A Simulation Model of Annual Pasture Production Limited by Rainfall and Nitrogen. In: Frissel, M.J. and Van Veen, J.A., Eds., Simulation of Nitrogen Behaviour of Soil-Plant Systems, Pudoc, Wageningen, 192-220.
 Neitsch, S.L., Arnold, J.G. and Kiniry, J.R. (2009) Williams. Soil and Water Assessment Tool: Theoretical Documentation. Version 2009. http://swat.tamu.edu/documentation/
 Raun, W.R., Johnson, G.V., Phillips, S.B. and Westerman, R.L. (1998) Effect of Long-Term N Fertilization on Soil Organic C and Total N in Continuous Wheat Under Conventional Tillage in Oklahoma. Soil Tillage Research, 47, 323-330.
 White, M.J., Harmel, R.D. and Haney, R.L. (2012) Development and Validation of the Texas Best Management Practice Evaluation Tool (TBET). Journal of Soil and Water Conservation, 67, 525-535. https://doi.org/10.2489/jswc.67.6.525
 Srinivasan, R., Zhang, X. and Arnold, J. (2010) SWAT Unguaged: Hydrological Budget and Crop Yield Predictions in the Upper Mississippi River Basin. Transaction of the ASABE, 53, 1533-1546.
 Kiniry, J.R., Arnold, J.G. and Xie, Y. (2002) Application of Models with Different Spatial Scales. In: Lajpat, R.A., Ma, L. and Howell, T.A., Eds., Agricultural System Models in Field Research and Technology Transfer, CRC Press, Boca Raton.
 Gupta, H.V., Sorooshian, S. and Yapo, P.O. (1999) Status of Automatic Calibration for Hydrologic Models: Comparison with Multilevel Expert Calibration. Journal of Hydrologic Engineering, 4, 135-143. https://doi.org/10.1061/(ASCE)1084-0699(1999)4:2(135)
 Thomason, W.E., Raun, W.R. and Johnson, G.V. (2000) Winter Wheat Fertilizer Nitrogen Use Efficiency in Grain and Forage Production Systems. Journal of Plant Nutrition, 23, 1505-1516.
 Macnack, N., Khim, B.C., Mullock, J. and Raun, W. (2014) In-Season Prediction of Nitrogen Use Efficiency and Grain Protein in Winter Wheat (Triticum aestivum L.). Communications in Soil Science and Plant Analysis, 45, 2480-2494.
 Johnson, G.V. and Raun, W.R. (1995) Nitrate Leaching in Continuous Winter Wheat: Use of Soil-Plant Buffering Concept to Account for Fertilizer Nitrogen. Journal of Production Agriculture, 8, 486-491. https://doi.org/10.2134/jpa1995.0486
 Lobell, D.B, Ortiz-Monasterio, J.I. and Anser, G.P. (2004) Relative Importance of Soil and Climate Variability for Nitrogen Management in Irrigated Wheat. Field Crops Research, 87, 155-165. https://doi.org/10.1016/j.fcr.2003.10.004