Comparison of performance of tile drainage routines in SWAT 2009 and 2012 in an extensively tile-drained watershed in the Midwest

Subsurface tile drainage systems are widely used in agricultural watersheds in the Midwestern US and enable the Midwest area to become highly productive agricultural lands, but can also create environmental problems, for example nitrate-N contamination associated with drainage waters. The Soil and Water Assessment Tool (SWAT) has been used to model watersheds with tile drainage. SWAT2012 revisions 615 and 645 provide new tile drainage routines. However, few studies have used these revisions to study tile drainage impacts at both field and watershed scales. Moreover, SWAT2012 revision 645 improved the soil moisture based curve number calculation method, which has not been fully tested. This study used long-term (1991–2003) field site and river station data from the Little Vermilion River (LVR) watershed to evaluate performance of tile drainage routines in SWAT2009 revision 528 (the old routine) and SWAT2012 revisions 615 and 645 (the new routine). Both the old and new routines provided reasonable but unsatisfactory (NSE < 0.5) uncalibrated flow and nitrate loss results for a mildly sloped watershed with low runoff. The calibrated monthly tile flow, surface flow, nitrate-N in tile and surface flow, sediment and annual corn and soybean yield results from SWAT with the old and new tile drainage routines were compared with observed values. Generally, the new routine provided acceptable simulated tile flow (NSE= 0.48–0.65) and nitrate in tile flow (NSE= 0.48–0.68) for field sites with random pattern tile and constant tile spacing, while the old routine simulated tile flow and nitrate in tile flow results for the field site with constant tile spacing were unacceptable (NSE= 0.00–0.32 and −0.29–0.06, respectively). The new modified curve number calculation method in revision 645 (NSE= 0.50–0.81) better simulated surface runoff than revision 615 (NSE=−0.11–0.49). The calibration provided reasonable parameter sets for the old and new routines in the LVR watershed, and the validation results showed that the new routine has the potential to accurately simulate hydrologic processes in mildly sloped watersheds.


Introduction
Subsurface drainage systems have been built up since 1870 and become common practices in agricultural watersheds in the Midwest to alleviate the damage caused by uneven drainage (Jaynes and James, 2007).Subsurface drainage allows excess water to leave the soil profile through perforated tubes installed below the soil surface.Water flows into the perforated tubes through cracks between adjacent tiles or holes in the tube and drains away when the water table is higher than the tile (Sugg, 2007).The soil horizontal hydraulic conductivity is increased and makes water drainage from soils to ditches or subsurface drains effective; the soil Published by Copernicus Publications on behalf of the European Geosciences Union.vertical hydraulic conductivity is large enough to prevent crop damage from flooding (Mitchell et al., 2003;Guo et al., 2012a, b).In this way, subsurface drainage systems enable large regions of the Midwestern US to become some of the most productive agricultural lands.
Subsurface drainage plays an important role in water balance in the poorly drained soils of Midwestern agricultural lands.For instance, the Little Vermilion River (LVR) watershed, an extensively tile-drained watershed in Illinois, has altered hydrology from an extensive subsurface drainage system network (Zanardo et al., 2012).Lal et al. (1989) studied tillage-caused alterations in water balance and sediment transport for a corn-soybean rotation in Ohio, and the results demonstrated that the percentages of annual precipitation drained by tiles in plowed conditions and on no-till plots are 33 to 58 % and 28 to 59 %, respectively.
Intensive tile drainage systems also create environmental problems, due to contaminants such as nitrate-N and pesticides in the water they transport.Subsurface drainage systems of agricultural fields in the Midwest have provided the majority of the nitrate that enters the Mississippi River and contributes to hypoxia in the Gulf of Mexico (Guo et al., 2018;Jaynes and James, 2007;Kalita et al., 2006).Subsurface tile drainage systems can increase nitrate and pesticide transport, because they remove excess water and convey soluble nitrate-N from the crop root zone.Nitrate coming from tile drains, especially for storm nitrate losses through tile drainage, has been considered the main source of nitrate in rivers and streams in the Mississippi River basin and nitrate export rate has been shown to be positively correlated with precipitation amount (P < 0.01) (Cuadra and Vidon, 2011).Additionally, 89-95 % of nitrate losses in a ditch catchment were transported by the tile drainage system of the catchment (Tiemeyer et al., 2008).Water discharge and nutrient loads of the Mississippi River reduce light penetration, and increase aquatic habitat loss and hypoxia in the northern Gulf of Mexico, the largest zone of oxygen-depleted coastal waters in the US (Diaz and Solow, 1999;Rabalais et al., 1999).A decrease in nutrient loading from Mississippi River discharge could alleviate Gulf of Mexico hypoxia (Rabalais et al., 1999).
Field-and watershed-scale models have been used to evaluate nutrient reduction strategies, and it is important to accurately simulate hydrological processes of tile drainage systems for evaluation of hydrological and water quality impacts of conservation practices in watersheds in the Midwest (Guo et al., 2018).The Soil and Water Assessment Tool (SWAT), a physically based and watershed-scale hydrological model, has been widely used to simulate land use change impacts on water quantity and quality (Basheer et al., 2016;Guo et al., 2015;Luo et al., 2012;Shope et al., 2014;Teshager et al., 2016;Wang et al., 2016;Yin et al., 2016), but studies on simulation of tile drainage impact at field and watershed scales using the new tile drainage routine from SWAT2012 are few (Boles et al., 2015;Du et al., 2005Du et al., , 2006;;Moriasi et al., 2005Moriasi et al., , 2012)).For instance, Sui and Frankenberger (2008) quantified the impact of tile drains on nitrate loss in an extensively tiledrained watershed, and showed that simulated nitrate loss results by SWAT2005 could be used for simulation of nitrate reductions at the watershed scale.Moriasi et al. (2012) used the new tile drain equations in SWAT to evaluate hydrology of a watershed in Iowa and determined value ranges for the new tile drain parameters, finding that Hooghoudt steadystate and Kirkham tile drain equations could be alternative tile drain simulation methods in SWAT.Boles et al. (2015) tested a new tile drainage routine in a watershed in Indiana using SWAT and found that the new tile drainage routine in SWAT2012 has the potential to predict tile flow and nitrate transported by tiles.
It is important to accurately simulate tile drains in hydrological models to correctly predict hydrologic processes and model nutrient reduction strategies in the Mississippi River system to alleviate hypoxia in the Gulf of Mexico.More information about application of realistic parameters for SWAT2012 tile drainage is needed.It is necessary to test and calibrate the new drainage routines in SWAT2012 and compare the modeled results with those by the old routines in SWAT2009, especially for dynamics between longterm monthly observed water quantity and quality data at both field and watershed scales.Therefore, the main objective of this study is to compare simulated flow, tile flow, runoff, nitrate in tile flow, and sediment load results for the new tile drainage routines in SWAT2012 and the old one in SWAT2009 at subsurface, surface, and river stations in an extensively tile-drained watershed and determine which routine provides a better model fit with observed values, to help improve understanding of tile drainage systems and evaluate the impacts of conservation practices on nitrate load reductions at both field and watershed scales in the Mississippi River system in further research.We calibrated and validated SWAT models with the new and old tile drainage routines to simulate tile flow and nitrate in tile flow at subsurface stations, surface runoff, and sediment and nitrate in surface runoff at surface stations, and streamflow, and sediment and nitrate in streamflow at the river station, and compared their performance.We also considered the new tile drainage routine with the improved curve number calculation method.We then determined which tile drainage routine can provide a better model fit.The research results have been used for simulation of the impacts of bioenergy crop scenarios on streamflow, tile flow, and sediment and nitrate losses in the LVR (Guo et al., 2018), and has also provided guidance to the selection of parameter sets for phosphorus reduction simulation at agricultural fields in northwestern Ohio using the Agricultural Policy/Environment eXtender (APEX), to achieve the nutrient reduction goal in Lake Erie.Moreover, the research results can allow selection of the most appropriate tile drainage routine and reasonable parameter sets, to guide evaluation of performance of conservation practices in reducing nutrient loads at both field and watershed scales in mildly Hydrol.Earth Syst.Sci., 22, 89-110, 2018 www.hydrol-earth-syst-sci.net/22/89/2018/ sloped watersheds in the Midwest with subsurface drainage systems.
2 Materials and methods

Tile drainage routines in SWAT
The SWAT model has simulated subsurface drainage since early versions (Boles et al., 2015).Arnold et al. (1999) enhanced SWAT2000 with a subsurface tile flow component and tested the enhanced model (SWAT2002) at a field scale with satisfactory results.However, because pothole impacts had not been included in SWAT2002 and the tile drainage routines were old, the SWAT2002 tile drainage method was not adequate for simulating tile flow and streamflow at a watershed scale (Arnold et al., 1999;Du et al., 2005).The equation used for tile drainage simulation in SWAT2002 (Neitsch et al., 2011) is where tile wtr is the amount of water removed from the layer on a given day by tile drainage (mm H 2 O), SW ly is the water content of the layer on a given day (mm H 2 O), FC ly is the field capacity water content of the layer (mm H 2 O), and t drain is the time required to drain the soil to field capacity (hs) (Neitsch et al., 2011). Du et al. (2005) created an impervious layer and improved the simulation of water table dynamics, and monthly flow and subsurface tile drainage simulated by SWAT2005 are much better than those simulated by SWAT2002.The time to drain soils to field capacity (TDRAIN) was used to determine the flow rate.Additionally, a new coefficient GDRAIN, the drain tile lag time, was introduced and used as the portion of the flow from tile drains into the streams on a daily basis (Du et al., 2006).Some studies have shown that the tile drainage routine in SWAT2005 could simulate the influence of subsurface drainage on hydrology at a watershed scale (Koch et al., 2013;Sui and Frankenberger, 2008).However, using a drawdown time (TDRAIN) method to simulate tile drains is simplified and limited.Equation (2) (Neitsch et al., 2011) is used for tile drainage simulation in SWAT2005: where tile wtr is the amount of water removed from the layer on a given day by tile drainage (mm H 2 O), h wtbl is the height of the water table above the impervious zone (mm), h drain is the height of the tile drain above the impervious zone (mm), SW is the water content of the profile on a given day (mm H 2 O), FC is the field capacity water content of the profile (mm H 2 O), and t drain is the time required to drain the soil to field capacity (hs) (Neitsch et al., 2011).
A new drainage routine which includes the use of the Hooghoudt and Kirkham drainage equations was used to simulate real-world drainage systems more accurately (Moriasi et al., 2005(Moriasi et al., , 2012)).Based on measured streamflow data from a watershed in Iowa, SWAT with the new tile drain equations was evaluated.The water balance components were simulated, and the results showed that the modified SWAT with the Hooghoudt steady-state and Kirkham tile drain equations simulated flow well (Moriasi et al., 2012).The new tile drainage routines (Eqs.3, 4, and 5) added to SWAT2005 are shown below.
When the water table is below the surface and ponded depressional depths are below a threshold, the Hooghoudt steady-state equation is used to compute drainage flux: where q is the drainage flux (mm h −1 ), m is the midpoint water table height above the drain (mm), K e is the effective lateral saturated hydraulic conductivity (mm h −1 ), L is the distance between drains (mm), and d e is the equivalent depth of the impermeable layer below the tile drains.When the water table completely fills the surface and ponded water remains at the surface for long periods of time, drainage flux is computed using the Kirkham equation (Moriasi et al., 2012(Moriasi et al., , 2013)): where t is the average depressional storage depth (mm), b is the depth of the tile drain from the soil surface (mm), r is the radius of the tile drain (mm), and δ is a dimensionless factor, determined by an equation developed by Kirkham (1957).
When predicted drainage flux is greater than the drainage coefficient, then the drainage flux is set equal to the drainage coefficient where q is the drainage flux (mm h −1 ) and DRAIN_CO is the drainage coefficient (mm day −1 ) (Moriasi et al., 2012(Moriasi et al., , 2013)).These new tile drainage routines have been used to model tile flow at watershed scale since SWAT2009 (Boles et al., 2015).Additionally, the drainage coefficient (DRAIN_CO) was included in the new tile drainage routine in SWAT2012 to control peak drain flow.

Study area
The LVR watershed (Fig. 1) is located in eastern-central Illinois and drains approximately 518 km 2 .Eighty-five percent of the watershed area is in eastern Vermilion County, 13 % of the watershed is in Champaign County, and 2 % of the watershed is in Edgar County.The LVR watershed consists of flat topography, with elevations ranging from 235 m in the headwaters to 174 m at the watershed outlet and with an average slope reaching at most 1 % (Zanardo et al., 2012).The  long-term (1991-2000) average annual precipitation for the watershed is 990 mm yr −1 (Kalita et al., 2006).
The watershed was subdivided into two subwatersheds, the upstream contributing areas of Georgetown Lake and the LVR.Ninety percent of the LVR watershed is agricultural land used for corn and soybean production, and the remainder consists of grassland, forest land, roadways, and farmsteads (Kalita et al., 2006).Annual area planted with soy- beans is equal to the area for corn planting (Algoazany et al., 2007).The dominant soil associations in the LVR watershed are Drummer silty clay loam and Flanagan silt loam (Zanardo et al., 2012;Keefer, 2003), and the dominant hydrologic soil groups are B and C. The LVR watershed is a typical tile-drained watershed in Illinois.Water quantity and quality data for this watershed are available from a long-term (1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003) monitoring project through which data were collected from several subsurface stations, surface stations, river stations and wetland sites in the watershed (Mitchell et al., 2003;Kalita et al., 2006).Based on long-term field observation data (1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000) from the watershed, Mitchell et al. (2003) and Kalita et al. (2006) studied the hydrology of flat upland watersheds in Illinois and demonstrated that the water could remain ponded on the soil surface until it would evaporate, seep or flow to the subsurface when the precipitation rate exceeds the infiltration rate of rainfall events, and surface runoff could flow into the streams directly during extremely large rainfall events.

Sites and data for model setup
Two subsurface stations (B and E), two surface runoff stations (Bs and Es), and one river station (R5), with drainage areas of 0.03, 0.076, 0.03, 0.023, and 69 km 2 , were selected for this study (Fig. 1 and Table 1).The drainage areas were determined from hand drawn tile layout with locations of tile lines and monitoring stations at the subsurface and surface stations (Singh et al., 2001a).Subsurface sites B and E were close to surface stations Bs and Es, respectively.B and E had similar land use, cropping systems and tile drainage systems with Bs and Es, respectively (Table 1).Elevation, soil, land use and weather data were used for SWAT model setup (Table 1).Daily water discharge data are available at subsurface, surface runoff, and river stations monitored by the Illinois Agricultural Experiment Station, University of Illinois at Urbana-Champaign.Water samples were obtained biweekly, while additional samples were taken by pump samplers during increased flow (Kalita et al., 2006).Daily nitrate and sediment load was computed by multiplying water dis- Negative value for CN II, value for USLE_C{19}, and USLE_C{ 56} is relative change to the default value.{19} and {56} represent corn and soybean, respectively.
charges with nitrate and sediment concentration, respectively (Yuan et al., 2000).Nitrate and sediment concentrations were not available for every day that water discharge occurred, and available data contained more water discharge measurements than nitrate and sediment concentration measurements.Nitrate and sediment loads were computed by multiplying the concentration at a specific time by half the flow volume since the last concentration measurement plus half the flow volume from the concentration measurement to the next concentration measurement (Kalita et al., 2006).Daily tile flow, surface runoff, nitrate load in tile flow, surface runoff, streamflow, and sediment load in surface runoff and streamflow were aggregated into monthly data and adopted in this study for model calibration and validation (Table 1).Other stations were not considered due to the quality of their data (Zanardo et al., 2012).Corn and soy-bean planting, harvest, and tillage practice data were collected from landowners (Table 1).

Modification to the soil moisture retention parameter calculation method
The tile drainage routine based on drawdown time in SWAT2009 Revision 528 (Rev.528) was called the "old routine" in this study.The tile drainage routine based on the Hooghoudt and Kirkham equations with a DRAIN_CO in SWAT2012 Revision 615 (Rev.615) was called the "new routine" in this study.SWAT Revision 645 (Rev.645) added a retention parameter adjustment factor (R2ADJ) to Rev. 615 to modify the soil moisture retention parameter calculation method (Eqs.5 and 6) (Neitsch et al., 2011).
where S is the retention parameter for a given day (mm), CN is the curve number for the day, S max is the maximum value the retention parameter can achieve on any given day (mm), SW is the soil water content of the entire profile excluding the amount of water held in the profile at wilting point (mm H 2 O), and w 1 and w 2 are shape coefficients.
rto3 is the fraction difference between CN III and CN I retention parameters, rtos is the fraction difference between CN = 99 (CN max ) and CN I retention parameters, SW fc is the amount of water held in the soil profile at field capacity, and SW sa is the amount of water held in the soil profile at saturation.In Rev. 645, R2ADJ was used to modify shape coefficients, w1 and w2, to increase S and thus decrease CN.R2ADJ ranges from 0 to 1 (Eqs.10, 11 and 12).When R2ADJ is 0, CN II is calculated when soil water content is at field capacity.When R2ADJ is 1, CN II is calculated when soil water content is at saturation.In this case, CN is decreased gradually based on soil from capacity to saturation, which is more reasonable than decreasing CN directly.In reality, CN II could be calculated when soil water content is near saturation (CN II < 100) rather than exactly at satura-tion (CN II = 100) (Neitsch et al., 2011).
MSW fc is the modified amount of water held in the soil profile at field capacity, and R2ADJ is the newly added retention parameter adjustment factor.

Model setup
SWAT2012 in conjunction with ArcGIS10.1 was used to simulate the LVR watershed.The 30 m National Hydrography Dataset (NHD) was used to generate a clipped stream layer for the LVR watershed into the simulation, and subbasins in the LVR watershed were delineated.Land-use data (NLCD 2006) for the study area were obtained from USGS.
The National Map Viewer and SSURGO from USDA Web Soil Survey were added to ArcSWAT (Table 1).HRUs were defined using the following thresholds: 0 % land use, 10 % soil, and 0 % slope.Daily precipitation data from rain gauge stations at sites B and E and 6 km southeast of site R5 were added in ArcSWAT and used for simulation at sites B and Bs, sites E and Es, and site R5, respectively (Table 1).Daily temperature, solar radiation, wind speed, and relative humidity data from an Illinois State Water Survey (ISWS) station (Champion Station, latitude: 40.08 • , longitude: −88.24 • , elevation: 219 m) closest to the LVR watershed were used (Table 1).
Management operation data for corn and soybean growth at sites B and E were collected (Table 1).Fertilizer was applied 10 days before planting at the rates of 218 kg ha −1 for anhydrous ammonia and 67 kg ha −1 for P 2 O 5 .Atrazine was applied at 2.2 kg ha −1 3 days before planting during corn growing years.P 2 O 5 fertilizer was applied at 56 kg ha −1 14 days before planting during soybean production years.
Tile drainage area was determined in HRUs where corn or soybeans were the current land use, slope was lower than 5 %, and soil drainage was somewhat poorly drained, poorly drained, or very poorly drained (Boles et al., 2015;Sugg, 2007;Sui and Frankenberger, 2008), and the tile drained area of the LVR watershed is about 75 %.

Parameter adjustments before model calibration
Plant growth parameters for corn and soybean growth simulation at sites B and E were adjusted.Radiation-use efficiency (BIO_E) and harvest index for optimal growing conditions (HVSTI) values for corn growth ranged from 32 to 39, and from 0.41 to 0.54, respectively, based on various studies (Edwards et al., 2005;Kiniry et al., 1998 (Edwards and Purcell, 2005;Mastrodomenico and Purcell, 2012;Sinclair and Muchow, 1999).The plant growth parameters for corn and soybean growth simulation of sites B and E were adjusted (Table 2).Cibin et al. (2016) adjusted BIO_E and potential heat units (PHU) for corn growth, and PHU, minimum temperature for plant growth (T_BASE), HVSTI, and normal fraction of phosphorus in yield (CPYLD) for soybean growth (Table 2) to reasonably simulate corn and soybean yields for two watersheds in the Midwestern US.This study adopted the same adjustment for corn and soybean growth simulation.
For surface runoff simulation, curve number calculation based on the soil moisture (ICN = 0) method was included in model calibration.Tile drainage simulation parameters were adjusted for the new routine.For Revs. 615 and 645, tile depth ranged from 1.05 to 1.1 m at various sites (Drablos et al., 1988;Singh et al., 2001a), and tile depth (DDRAIN) was set as 1.075 m in the model.The maximum depressional storage selection flag/code (ISAMX) was used to control the method used to calculate the static maximum depressional storage parameter (SSTMAXD), representing the surface storage.When ISMAX is 0, SSTMAXD is allowed to be defined by the user, while when ISMAX is 1, SSTAMXD is dynamically calculated based on rainfall and tillage practices (Moriasi et al., 2005;Moriasi et al., 2012).In this study, ISMAX was set as 0 and SSTMAXD was set as 12 mm, based on the previous DRAINMOD (Skaggs et al., 2012) and SWAT studies (Boles et al., 2015).DRAIN_CO, the amount of water drains in 24 h, was set as 20 mm day −1 , describing the size of the main collector drain pipes and the outlet (Sui and Frankenberger, 2008).Tile spacing (SDRAIN) for site B and site R5 was set as 28 000 mm, the same as the observed tile spacing for site E. Effective radius (RE) was used to simulate the entrance resistance into the perforations of tile drains, and was set as 15 mm (Boles et al., 2015;Moriasi et al., 2012).

Model calibration and validation
Revisions 528, 615, and 645 simulated tile flow at sites B and E were compared with the observed values to evaluate tile drainage simulation performance of the old and new routines and the new routine with a modified curve number calculation method.Revisions 528 and 615 simulated nitrate in tile flow at sites B and E was compared with the observed values to evaluate nitrate in tile flow simulation performance of the old and the new routines.Revisions 615 and 645 simulated surface runoff at sites Bs and Es was compared with the observed values to evaluate surface runoff simulation performance of the default soil moisture based curve number calculation method and modified curve number calculation method.Revisions 528 and 645 simulated flow at site R5 were compared with the observed values to evaluate flow simulation performance of the old and new routines.Revision 645 was not used for flow simulation at river station R5, because Rev. 645 could not run successfully for the mainly tile drained river station R5.This was thought to be because depth to impervious layer (DEP_IMP) values were too low and the impervious layer was too close to the soil profile, which may have affected the functionality of Rev. 645 in simulating groundwater and tile flow on a watershed level.
For Rev. 528, the calibrated values for tile flow simulation parameters at site B, time to drain soil to field capacity (TDRIAN), drain tile lag time (GDRIAN), and DEP_IMP were used for flow simulation at site R5.For Rev. 615, the calibrated values for tile flow simulation parameters at site B, DEP_IMP and multiplication factor to determine lateral saturated hydraulic conductivity (LATKSATF), were modified at site R5, to accurately simulate flow and obtain reasonable water budget results.

Model performance evaluation
Model outputs, annual corn and soybean yield, monthly tile flow and nitrate in tile flow at sites B and E, monthly surface runoff, sediment and nitrate in surface runoff at sites Bs and Es, and monthly flow, sediment and nitrate in flow at site R5 from the old and new routines were compared with observed values for model calibration and validation.Com-T.Guo et al.: Comparison of performance parison between simulated results from the old and new routines and observed values were plotted.The statistical methods used for verifying model performance included Percent bias/Percent error (P BIAS (%)), the coefficient of determination (R 2 ), the Nash-Sutcliffe model efficiency coefficient (NSE), the modified NSE (MSE) and the Kling-Gupta efficiency (KGE) (Eqs.13-17).
where Obs and Sim represent the ith observed and simulated monthly data, respectively.n is the total number of months.Obs and Sim represent the average values of the observed and simulated monthly data, respectively.α = σ Sim /σ Obs , β = µ Sim /µ Obs , and r is the linear regression coefficient between simulated and observed data (Eq.15).Percent bias (Gupta et al., 1999) can measure the average tendency of the simulated data to deviate from the observed data.A value of 0.0 is optimal for P BIAS , representing accurate model simulation.Negative values represent model overestimation bias, and positive values indicate model underestimation bias.If P BIAS ± 25 % for streamflow, ±55 % for sediment, and ±70 % for N and P, model simulation results can be considered satisfactory (Moriasi et al., 2007).The R 2 value indicates the strength of the linear relationship between the simulated and observed data.A R 2 value of greater than 0.5 is considered reasonable model performance (Moriasi et al., 2007).The NSE (Nash and Sutcliffe, 1970) can represent how well the plot of observed versus simulated data fits the 1 : 1 line.The NSE value ranges from −∞ to 1, and the optimal value is 1.A NSE value of greater than 0.5 is considered satisfactory model performance (Moriasi et al., 2007).A NSE value of 0 means that the simulated values are as accurate as the mean of the observed data, and a negative NSE value shows that the mean value of observed data is a better predictor than the simulated data, meaning unacceptable performance (Moriasi et al., 2007); 0.36 ≤ NSE ≤ 0.72 and NSE ≥ 0.75 have also been considered satisfactory and good simulated data, respectively (Larose et al., 2007;Van Liew et al., 2003).A modified form of the NSE (Eq.12) could decrease the oversensitivity of the NSE to extreme values (Krause et al., 2005), and is sensitive to chronic overor under-predictions.The KGE computes the Euclidian distance of the correlation, the bias, and a measure of variability.The use of KGE (Eq.13) improves the bias and the variability measure considerably and decreases the correlation slightly compared to the NSE (Gupta et al., 1999).The KGE value ranges from −∞ to 1.The closer to 1, the more accurate the model is.A KGE value of greater than 0.5 is considered a satisfactory simulated result (Gupta et al., 1999).

Calibration and validation results for corn and soybean yields
Simulated annual corn and soybean yields were compared with observed values during the calibration and validation periods at sites B and E (Fig. 2), and model performance in simulating crop yields was evaluated (Table 4).Performance of the simulated corn and soybean yields from Rev. 615 at sites B (Fig. 2a, b) and E (Fig. 2c, d) during the calibration and validation was satisfactory (Table 4).Simulated annual corn and soybean yields fit observed values well (Fig. 2).P BIAS values of corn and soybean yields during the calibration and validation at sites B and E ranged from −2 to 13 %, indicating accurate model simulation.During the calibration, R 2 , NSE, MSE, and KGE values for corn and soybean yields at both sites ranged from 0.75 to 0.99.During the validation, R 2 , NSE, MSE, and KGE values for corn and soybean yields at both sites ranged from 0.71 to 0.92 (Table 4).Adjusted crop growth parameters (Table 2) in Rev. 615 provided good predictions of corn and soybean yields.Generally, simulated corn and soybean yield results were improved compared to the simulated results from Root Zone Water Quality Model (RZWQM) at sites B and E (Singh et al., 2001b), since SWAT incorporated more details of crop management practices, such as pre-plant and post-harvest fertilizer application (Neitsch et al., 2011).

Calibration and validation results for tile flow, surface runoff and flow
This section outlines calibration and validation performance for monthly tile flow at subsurface stations B and E, surface runoff at surface stations Bs and Es, and flow at river station R5.

Calibration and validation results for monthly tile flow at sites B and E
Simulated monthly tile flows were compared with observed values during the calibration and validation at sites B (Fig. 3a, b) and E (Fig. 3c, d).Model performance in sim- ulating monthly tile flow and nitrate in tile flow at sites B and E was evaluated (Table 4).Performance of the simulated monthly tile flow from Revs.528, 615 and 645 at site B, and Revs.615 and 645 at site E during the calibration and validation was satisfactory (NSE > 0.5), except that NSE (0.48) from Rev. 615 during the validation at site E was slightly under the acceptable limit (Table 4).Generally, simulated tile flow results for the old routine from Rev. 528 were better than those for the new routine from Revs.615 and 645 at site B. However, simulated tile flow results from Revs.615 and 645 were better than those from Rev. 528 at site E. The modified curve number calculation method in Rev. 645 improved surface runoff simulation and then improved tile flow simulation compared to the default curve number calculation method based on soil moisture in Rev. 615 (Fig. 3a-d, and Table 4).Generally, P BIAS values of tile flow results from the three versions ranged from −16 to 24 % during the calibration and validation at two sites, indicating accurate model simulation, except that P BIAS (−28 %) from Rev. 645 during the calibration at site B, and P BIAS (49 %) from Rev. 528 during the validation at site E represented slightly overestimated and underestimated results, respectively.Generally, R 2 , NSE and KGE values for tile flow from the three versions were satisfactory (> 0.5), except that NSE (< 0.5) from Rev. 528 during the calibration and validation and NSE (0.48) from Rev. 615 during validation at site E were unacceptable, and slightly under the acceptable limit, respectively (Table 4).Revision 528 simulated tile flow fit observed data very well at site B (Fig. 3a, b, and Table 4).However, Rev. 528 simulated tile flows were overestimated at tile flow peaks in May 1996, June 1998, and May 2002 (Fig. 3c, d).Revision 528 simulated tile flows were underestimated from May to October in 1992, from June to November in 1994, from July in 1995 to March in 1996, from May in 1999 to February in 2000, from May to August in 2001, and from July to December in 2002 (Fig. 3c, d).The old routine in Rev. 528 had different performance at sites B and E, which was mainly caused by different soil and weather characteristics, tile pattern and cropping systems of the two sites (Table 1), and the physical process of simulating tile flow in the old routine.For instance, site B had clay silt loam soil, random tile pattern and reduced-tillage practices, while site E had silt loam soil, constant tile spacing and no-tillage practices (Table 1).The old routine in Rev. 528 has the potential to overestimate tile flow peaks, since simulated tile flow by the old routine was controlled by a simple drawdown time parameter (TDRAIN), and tiles were allowed to carry an unlimited maximum of water no matter how intense the rainfall.The calibrated TDRAIN values were 26 and 25 h for sites B and E, representing that it would take 26 and 25 hs to drain soils from saturation to field capacity at sites B and E, respectively (Table 4).The calibrated drain lag time (GDRAIN) values were 25 and 26 h for sites B and E, representing that there were 25 and 26 hrs lag time between water enters the tiles from soil and water enters the main channel from the tiles at sites B and E, respectively, which was used to smooth the tile flow hydrograph (Table 4).However, using a draindown time (TDRAIN) to determine tile flow rate was too simplified, since TDRIAN was a static value for tiles no matter whether there was a large storm or not.Thus, the old routine overestimated tile flow peaks for site E (Fig. 3c, d), which was consistent with tile flow simulation using the old routine in the Matson Ditch watershed in Indiana (Boles et al., 2015).Moreover, the old routine was used to simulate tile flow on days when the simulated height of the water table exceeded the height of the tile drain (Neitsch et al., 2011).Tile drainage systems can cause water table recession in tile-drained soil.Water table was lower when respiratory activity was highest in summer (Muhr et al., 2011), which may be lower than the depth of subsurface tiles during long dry summer periods.Water table depth calculation based on change in the soil water for the whole soil profile tended to overestimate the distance between water table and the soil surface when long-term simulations were performed, most commonly in cases where days without rainfall dominated (Moriasi et al., 2013).Thus, Rev. 528 simulated tile flow was zero during long dry summer periods.
The calibrated DEM_IMP for all three revisions at sites B and E ranged from 2765 to 3000 mm, representing the depth to the impervious layer, which also could determine the percent of potential seepage flows through this layer (0.0-1.0).Compared to the old routine in Rev. 528, the new routine in Revs.615 and 645 incorporates the DRAIN_CO, and tile flow peaks can be limited by the radius of the tile.In reality, subsurface drainage systems are designed with a drainage coefficient (DRAIN_CO), which is the amount of water that can be drained in 24 h.In this case, the tiles could flow for a slightly longer period, and simulated tile flow matched well with observed values at site E (Fig. 3c, d).In this study, the more physically based equations and the DRAIN_CO (20 mm day −1 ) in the new routine in Revs.615 and 645 can reduce the flashiness of the tile flow simulation and result in lower tile flow peak and longer recession.Moreover, the Kirkham equation was used in the new routine to calculate drainage flux when water table was lower than tiles, which improved tile drainage calculation during dry periods to the old routine.Simulated monthly tile flow Rev. 615 was similar to observed values at two sites, except that Rev. 615 sim-ulated tile flow could not capture tile flow peaks well in May of 1996 and1998 (Fig. 3a) and May of 2002 (Fig. 3b) at site B. Soil moisture was reduced during long dry periods from June of 1995 to April of 1996.Subsurface tile drains can lower the water table (Sui and Frankenberger, 2008), and long-term water depletion may drop the water table lower than the depth of tiles (DDRAIN, 1075 mm).For long-term water table depth simulation (19 years in this study), the computed water table depth may gradually drop as profile soil water decreases due to periods of higher ET, which makes it harder for the water table to rise to the surface after rain events (Moriasi et al., 2013).When water storage is higher than the height of the surface storage threshold (20 % of the static maximum depressional storage (SSTMAXD), 0.24 mm in this study) and water table is near the bottom of the soil surface, the Kirkham equation is used to calculate drainage flux (Boles et al., 2015).In this study, overestimation of water table depth might have caused the new routine not to trigger the Kirkham equation to calculate tile flow drainage even though 1996 was a wet year (annual precipitation was 1008 mm).The new routine in Rev. 615 resulted in decreased tile flow peaks and longer storage time (Boles et al., 2015).The new routine in Rev. 645 captured tile flow peaks well at two sites, although the differences between simulated and observed tile flow values were large in July 2003 at site B (Fig. 3a).For Rev. 645, the calibrated values of a newly added curve number calculation retention parameter adjustment factor (R2ADJ) were 0.80 to 0.85 at sites B and E, respectively.CN II value was calculated when soil water content was near saturation (Eq.10), which was realistic for a mildly sloped watershed with extensive tile drainage systems.With R2ADJ, CN II was decreased gradually based on soil from capacity to saturation, which was more reasonable than decreasing CN directly (Moriasi et al., 2013).The newly added curve number calculation retention parameter adjustment factor in Rev. 645 calculates curve numbers reasonably well based on the soil moisture retention curve, and can partition surface runoff and tile flow well.Thus, simulated tile flow results from Rev. 645 captured peaks well, and the differences between simulated and observed tile flow values were small after long dry periods (Fig. 3a-d).
Besides DRAIN_CO, DDRAIN, and SSTMAXD, the new routines in Revs.615 and 645 incorporate more tile drainage simulation parameters, SDRAIN (28 000 mm), RE (15 mm) and multiplication factor to determine lateral saturated hydraulic conductivity (LATKSATF), to represent tile drainage systems more realistically than the old routine in Rev. 528.LATKSATF was used to determined lateral hydraulic conductivity using the saturated hydraulic conductivity for each soil layer and soil type.The calibrated LATKSATF values for Revs.615 and 645 at subsurface stations B and E, surface station Bs and Es and river station R5 ranged from 1.0 to 1.4, which were reasonable based on the previous tile drainage studies in Iowa and the recommendations value (1.4) by the Iowa Drainage Guide (Cooperative Extension Service, 1987;Singh et al., 2007;Singh et al., 2006;Singh and Helmers, 2008).Simulated monthly tile flow results for Rev. 615 were better than the previous DRAINMOD and RZWQM simulated results at sites B and E (Singh et al., 2001a), since both DRAINMOD and RZWQM models overestimated daily tile flow at these sites to obtain an acceptable R 2 value (> 0.5), but they did not match well with the observed values generally from 1993 to 1998.Simulated monthly tile flow results for Rev. 615 at sites B and E were similar to the observed values, and obtained acceptable P BIAS , R 2 , NSE and KGE generally from 1991 to 2003.

Calibration and validation results for monthly surface runoff at site Bs and Es
This section describes calibration and validation performance for monthly surface runoff at surface sites Bs and Es (Fig. 4).The LVR watershed is dominated by agricultural land with an extensive tile drainage system.Direct surface runoff was a small percentage (≤ 15 %) of the streamflow in the LVR watershed, and was nearly zero for years 1995 and 1997, even though there was sufficient precipitation (Mitchell et al., 2003).Thus, it is challenging to simulate surface runoff in the LVR watershed.Performance of the modeled monthly surface runoff at sites Bs and Es during the calibration and validation was satisfactory from Rev. 645, and was unsatisfactory from Rev. 615.Generally, simulated surface runoff results from Rev. 645 with the improved curve number calculation method were better than those from Rev. 615 with the default soil moisture based curve number calculation method.Simulated surface runoff results from Rev. 645 were better than those from Rev. 615 for two sites (Fig. 4a-d and Table 4).Generally, simulated monthly surface runoff from Rev. 645 was similar to observed values.Revision 615 simulated surface runoff results were higher than observed values (Fig. 4a-d).For Rev. 615, the calibrated curve number (CN II) value (60) at sites Bs and Es was realistic for a watershed dominated by agricultural land based on the previous studies (Boles et al., 2015;Moriasi et al., 2012;Neitsch et al., 2011), and simulated surface runoff was overestimated (Fig. 4a-b).For Rev. 645, the calibrated values of the newly added curve number calculation retention parameter adjustment factor (R2ADJ) were 0.81 to 0.83 at sites Bs and Es, respectively.In this case the CN II value was calculated when soil water content was near saturation (Eq.10), which was reasonable for a mildly sloped watershed with low runoff (Neitsch et al., 2011).P BIAS values of surface runoff results from Rev. 615 during the calibration and validation at two sites ranged from −143 to −82 %, representing overestimated simulation results.P BIAS values of surface runoff results from Rev. 645 during the calibration and validation at field site Bs were 13 and 12 %, indicating accurate simulation results.P BIAS values of surface runoff results from Rev. 645 during the calibration and validation at site Es were  −25 and 28 %, indicating slightly overestimated and overestimated simulation results, respectively.Generally, NSE, MSE, and KGE values for simulated surface runoff results from Rev. 615 at the two sites were unacceptable (< 0.5) (Table 4).R 2 , NSE, MSE, and KGE values for simulated surface runoff results from Rev. 645 at the two sites were acceptable (> 0.5) (Table 4), except that MSE during the validation at site Bs (0.49) and MSE during the calibration at site Es (0.42) were slightly under the acceptable limit (Table 4).In this watershed with flat topography and dominated by tile drainage, surface runoff was small and nearly zero from 1994 May to 1996 March and from 1999 March to 2002 April for surface station Bs (Fig. 4a, b), and nearly zero from 1994 June to 1995 April and from 1998 July to 2002 March for surface station Es (Fig. 4c, d).

Calibration and validation results for monthly flow at site R5
Performance of the modeled monthly flow from Revs.528 and 615 at site R5 during the calibration and validation was satisfactory.Simulated monthly flow results from Rev. 528 were slightly better than those from Rev. 615 at site R5 (Fig. 3e, f, and Table 4).Generally, simulated monthly flow was similar to observed values (Fig. 3e, f).However, Rev. 528 simulated flow values were higher than observed values in May 1996 and December 1997 (Fig. 3e), which was mainly caused by overestimation of tile flow during these periods.Simulated tile flow by the old routine in Rev. 528 was controlled by a simple drawdown time parameter (TDRIAN), no matter how intense the rainfall.Thus, Rev. 528 has the potential to overestimate tile flow peaks.Revisions 528 and 615 simulated flow values were slightly higher than observed values from June to November of 1994, 1996, 1998 (Fig. 3e), and 2002 (Fig. 3f), which was mainly because of the overestimation of surface runoff during these periods.The calibrated CN II values were reduced by 20 % to accurately simulate streamflow at the river station (Table 4).Moreover, the calibrated surface runoff lag coefficient (SURLAG) ranged from 0.2 to 0.3 for three revisions at all sites, showing that the model allowed a small portion of surface runoff to reach the main channel when the time of concentration is greater than 1 day, and could smooth the simulated flow hydrograph at site R5 (Neitsch et al., 2011).The calibrated soil evaporation compensation factor (ESCO) values at five sites ranged from 0.88 to 0.91 at all sites, which meant that the reduction of ESCO would allow lower soil layers to compensate for a water deficit in upper layers and increase ET and reduce surface runoff (Jha, 2011).These values were reasonable for a watershed dominated by agricultural land based on the previous studies (Boles et al., 2015;Moriasi et al., 2012;Neitsch et al., 2011) and have provided reasonable water balance component proportions when used for modeling the impacts of various bioenergy crop scenarios on hydrology and water quality in the LVR (Guo et al., 2018).Revisions 528 and 615 simulated flow values were lower than observed values from January 2000 to February 2001 (Fig. 3e), which was mainly caused by underestimation of tile flow.Since the water table was lower than the tiles after the long dry periods in 1999, the old routine in Rev. 528 could not simulate tile flow, and the new routine in Rev. 615 could not use the Kirkham equation to calculate tile drainage flux.P BIAS values of flow results from Rev. 528 during the calibration and validation were −16 and −20 %, respectively, representing fairly accurate simulation results.P BIAS values of flow results from Rev. 615 during the calibration and validation were −26 and −37 %, respectively, indicating overestimated simulation results.Generally, R 2 , NSE, and KGE values for simulated flow results from Revs.528 and 615 were satisfactory (> 0.5), except that MSE (0.30) from Rev. 615 during the calibration was under the acceptable limit (Table 4).Simulated average annual tile flow values from Rev. 528 (128 mm) and Rev. 615 (129 mm) were 14 and 15 % of total precipitation, respectively, over the period from 1992 to 2003.Simulated average annual ET values from Rev. 528 (585 mm) and Rev. 615 (571 mm) were 71 and 69 % of total precipitation, respectively.Simulated average annual water yield values from Rev. 528 (248 mm) and Rev. 615 (265 mm) were 27 and 29 % of total precipitation, respectively.Flow partitioning appeared reasonable for simulated results from Revs.528 and 615 based on the previous watershed-scale tile drainage simulation studies (Boles et al., 2015;Moriasi et al., 2012;Moriasi et al., 2013).Major flow paths are important in determining sediment and nitrate loads.
The new tile drainage routine in Rev. 615 was improved compared to the old routine in Rev. 528 (Fig. 3e, f, and Table 4).Revision 528 could not simulate tile flow once the water table was lower than tile depth, while Rev. 615 could simulate tile flow by the Hooghoudt equation once the water table dropped after long dry periods during the validation (Fig. 3e).Revision 615 incorporated tile parameters, such as DRAIN_CO, DDRAIN, LATKSATF, RE, and SDRAIN, to represent characteristics of the tile drainage system which simulate tile flow more realistically.Some processes in Rev. 615 could be improved.For instance, DEP_IMP can represent depth to impervious layer and soil permeability and can be separated in the model.Water table depth calculation can determine which equation will be used for tile flow simulation, and water table depth calculation during long dry periods can be improved to better simulate tile flow.

Calibration and validation results for sediment losses in surface runoff and flow
This section outlines calibration and validation performance for monthly sediment losses in surface runoff at surface stations Bs and Es, and monthly sediment losses in flow at river station R5.The calibrated minimum value of USLE C factor for water erosion (USLE_C) was increased from 12 to 15 % for corn and from 6 to 7 % for soybean at sites Bs, Es, and R5 (Table 3), to increase the generation of sediment (Qiu et al., 2012).The calibrated peak rate adjustment factor for sed-iment routing in the subbasin (ADJ_PKR) ranged from 1.1 to 1.2 at sites Bs, Es, and R5 (Table 3), which was used to adjust the amount of erosion generated in the HRUs (Neitsch et al., 2011).The exponent parameter for calculating sediment re-entrained in channel sediment routing (SPEXP) was calibrated as 1.5 for Revs.528 and 615 at site R5 (Table 3), which could be used to determine the maximum amount of sediment that can be re-entrained during channel sediment routing (Sexton et al., 2011).The calibrated value of channel erodibility factor (CH_COV1) was 0.3 for Revs.528 and 615 at site R5 (Table 3), which could alter channel erosion and sediment re-entrainment (Qiu et al., 2012).Performance of the modeled monthly sediment load from Rev. 645 was reasonable during the calibration and satisfactory during the validation for site Bs (Fig. 5a, b, Table 4), and was reasonable during the calibration and validation for site Es (Fig. 5c, d, Table 4).Simulated monthly sediment load from Rev. 645 was similar to observed values at two sites (Fig. 5a-d), except that simulated sediment load was lower than the observed value for March 1999 at site Bs (Fig. 5a), and for May 1996, and April and May of 2002 at site Es (Fig. 5c, d).Revision 645 overestimated sediment load for June 1998 and May 2002 at site Bs (Fig. 5a, b), and for June 1998 at site Es (Fig. 5c).P BIAS values of Rev. 645 simulated sediment load results during the calibration were 20 and −8 % at sites Bs and Es, respectively, indicating accurate simulation results (Table 4).P BIAS values of Rev. 645 simulated sediment load results during the validation were −77 and 86 % at sites Bs and Es, indicating overestimated and underestimated simulation results, respectively (Table 4).R 2 , NSE, MSE, and KGE values for simulated sediment during the calibration and validation were satisfactory at site Bs (> 0.5), except that R 2 (0.38) and NSE (0.27) during the calibration were under the acceptable limit (Table 4).R 2 , NSE, MSE, and KGE values for simulated sediment during the calibration and validation were unsatisfactory at site Es (< 0.5) (Table 4), except that KGE (0.78) during the calibration and R 2 (0.56) and MSE (0.50) during the validation were acceptable.Simulated sediment could not capture the sediment peak well at the two sites (Fig. 5a-d), and performance evaluation methods are sensitive to high values.The magnitude of sediment load for the mildly sloped sites Bs and Es was small; thus, Rev. 645 simulated results were reasonable even though it had difficulty in capturing sediment load peaks well (Fig. 5a-d).
Performance of the modeled monthly sediment load in flow from Revs.528 and 615 at site R5 during the calibration and validation was reasonable (Fig. 5e, f, and Table 4).Simulated monthly sediment loads in flow results from Rev. 528 were better than those from Rev. 615 during the calibration at site R5 (Fig. 5e and Table 4).Simulated monthly sediment load from Revs.528 and 615 matched observed values fairly well, except that both the old and new routines could not capture sediment load peaks well (Fig. 5e, f).This was caused by the failure to predict surface runoff well.4).Generally, R 2 , NSE, MSE, and KGE values for simulated sediment were unsatisfactory (< 0.5), except that R 2 (0.96) from Rev. 528 and R 2 (0.95) from Rev. 615 during the validation were acceptable (Table 4).However, the LVR watershed is a mildly sloped watershed with extensive tile drainage systems, which was dominated by tile flow, and surface runoff and sediment in surface runoff were low, and it was challenging to simulate sediment load accurately.Revisions 528 and 615 simulated sediment load had difficulty in matching sediment load peaks (Fig. 5e, f), and performance evaluation results were unacceptable generally (Table 4).However, simulated sediment load can still be considered reasonable, since the magnitude of sediment load in this mildly sloped watershed was small (Fig. 5e, f).

Calibration and validation results for
nitrate-nitrogen losses in tile flow, surface runoff and flow This section outlines calibration and validation performance for monthly nitrate-nitrogen losses in tile flow at subsurface stations B and E, monthly nitrate-nitrogen losses in surface runoff at surface stations Bs and Es, and monthly nitratenitrogen losses in flow at river station R5.The amount of nitrate removed from surface runoff relative to the amount removed through percolation was controlled by the nitrate percolation coefficient (NPERCO), the calibrated value of which ranged from 0.12 to 0.15 at all sites (Table 3).As NPERCO decreased from the default value (0.20) to the calibrated values, nitrate concentration in runoff would be decreased (Neitsch et al., 2011).The calibrated denitrification threshold water content (SDNCO) at all sites ranged from 0.9 to 1.1, which was in a reasonable range based on the previous studies (Boles et al., 2015;Moriasi et al., 2013) (Table 3).
Hydrol.Earth Syst.Sci., 22, 89-110, 2018 www.hydrol-earth-syst-sci.net/22/89/2018/ SDNCO was used to determine denitrification level and the calibrated SDNCO could allow reasonable amounts of denitrification to occur and then realistic amounts of nitrate to exit tiles (Boles et al., 2015;Neitsch et al., 2011).Then the denitrification level could be fine-tuned by the denitrification exponential rate coefficient (CDN), the calibrated values of which ranged from 0.05 to 0.06 (Table 3).Generally, calibrated nitrate loads in tile flow results were improved compared to nitrate concentration simulation using RZWQM and the nutrient component of DRAINMOD (DRAINMOD-N) at sites B and E (Singh et al., 2001b).Performance of the simulated monthly nitrate in tile flow from Rev. 528 at site B and Rev. 615 at sites B and E during the calibration and validation was satisfactory, and from Rev. 528 at site E during the calibration and validation it was unsatisfactory (Fig. 6a-d, and Table 4).Generally, simulated nitrate in tile flow results for the old routine from Rev. 528 were slightly better than those for the new routine from Rev. 615 at site B (Fig. 6a, b, and Table 4).However, simulated nitrate in tile flow results for the new routine were better than those for the old routine at site E (Fig. 6c, d, and  Table 4).Generally, simulated monthly nitrate in tile flow matched observed values well, except that Revs.528 and 615 simulated nitrate in tile flow at site B could not capture peaks well in May 1996, February 1997and May 2002 (Fig. 6a, b), and Rev. 615 simulated nitrate in tile flow at site E was overestimated in November and December 1994 and April 2001, and underestimated in in July 1992, May 1995, and from May to October 2002 (Fig. 6c, d), which was caused by the failure to predict tile flow correctly during these periods (Fig. 3a, b).Generally, Rev. 528 simulated nitrate in tile flow at site E during the calibration and validation was underestimated (Fig. 6c, d), which was likely caused by the failure to predict accurate tile flow (Fig. 3c, d).P BIAS values of nitrate in tile flow results from Revs.528 and 615 at site B and Rev. 615 at site E during the calibration and validation ranged from 24 to 39 %, indicating accurate model simulation (Table 4).Generally, R 2 , NSE, MSE and KGE values for simulated nitrate in tile flow from Revs.528 and 615 at site B and Rev. 615 at site E were satisfactory (> 0.5).However, NSE (0.48) from Rev. 615 during the validation at site E was slightly under the acceptable limit, and MSE from Rev. 615 during the calibration (0.32) and validation (0.31) at site E were unacceptable (Table 4), which was due to underestimated nitrate values in tile flow after long dry periods (Fig. 6c, d).R 2 , NSE, MSE and KGE values for simulated nitrate in tile flow from Rev. 528 at site E were unsatisfactory (< 0.5), except that R 2 (0.5) during the validation was acceptable.
Performance of the modeled monthly nitrate load in surface runoff from Rev. 645 at site Bs during the calibration and validation and at site Es during the calibration was reasonable, and at Es during the validation it was satisfactory (Fig. 7a-d, and Table 4).Simulated monthly nitrate load was similar to observed values, except that simulated nitrate load values were lower than the observed values in May of 1996and 1998, January 1999, and May 2002at site Bs (Fig. 7a, b), in May 1996, April 1999, June 2000, October 2001, and April and May 2002 at site Es; P BIAS values of nitrate load results ranged from 59 to 99 % during the calibration and validation at sites Bs and Es, indicating underestimated model simulation.Generally, R 2 , NSE, MSE, and KGE values for simulated nitrate load during the calibration and validation at site Bs and during the calibration at site Es were unsatisfactory (< 0.5), except that MSE (0.54) and KGE (0.55) during the validation at site Bs were acceptable (Table 4).R 2 , NSE, and MSE values for simulated nitrate load during the validation at site Es were satisfactory (> 0.5) (Table 4).Revision 645 simulated nitrate in surface flows at sites Bs and Es was reasonable, as nitrate in surface runoff was low given the watershed was dominated by tile flow and surface runoff rarely occurred.
Performance of the modeled monthly nitrate load in flow from Revs.528 and 615 at site R5 during the calibration and validation was satisfactory (Fig. 6e, f, and Table 4).Simulated monthly nitrate loads in flow results from Rev. 615 were better than those from Rev. 528 during the calibration, and Rev. 528 simulated results were better than those from Rev. 615 during the validation at site R5 (Fig. 6e, f, and Table 4).Simulated monthly nitrate load was similar to observed values, except that Rev. 528 simulated nitrate load values were higher than observed values in May 1996, December 1997, and May 2002 (Fig. 6e, f), which was mainly caused by overestimation of tile flow during these periods.Revisions 528 and 615 simulated nitrate load values were lower than observed values during June 1997, January and February of 1999, May and June of 2000, and June 2002 (Fig. 6e, f), which was mainly caused by underestimation of tile flow during these periods.P BIAS values of nitrate load results were 27 and 6 % from Rev. 528 during the calibration and validation, and 31 and 23 % from Rev. 615 during the calibration and validation, indicating fairly accurate model simulation.Generally, R 2 , NSE, MSE, and KGE values for simulated nitrate load were satisfactory (> 0.5).However, NSE (0.43) and MSE (0.41) from Rev. 528 during the calibration were unsatisfactory, and MSE (0.49) from Rev. 615 during the calibration was slightly under the acceptable limit (Table 4).
Limitations of this work include limited observed rainfall data for site R5, water table depth calculation after long dry periods, and difficulty in simulating surface runoff, sediment, and nitrate in surface runoff from this extensively tile drained, mildly sloped watershed.Observed rainfall data for site R5 were from the closest rain gauge station located 6 km southeast of site R5, which may impact the accuracy of flow simulation.There is an opportunity to improve the representation of tile drainage systems in SWAT, especially for individual tiles, and improve Rev. 645 functionality at watershed scales.

Conclusions
In this study, the old tile drainage routine in SWAT2009 (Rev.528) and the new tile drainage routine in SWAT2012 (Revs.615 and 645) were used to simulate monthly tile flow, nitrate in tile flow, surface runoff, and sediment and nitrate in surface runoff at field sites, and monthly flow, sediment and nitrate in flow at a river station.Performance of both the old and new routines was evaluated and compared with observed values.
The results showed that Rev. 615 satisfactorily simulated corn and soybean yields at field sites, and both the old and new routines provided satisfactory tile flow and nitrate in tile flow results at subsurface sites, satisfactory flow and nitrate load in flow, and reasonable sediment load in flow results at the river station after model calibration.Revision 645 with an improved curve number calculation method provided satisfactory surface runoff, and reasonable sediment and nitrate load in surface runoff results at surface stations.
Generally, simulated tile flow results for the old routine from Rev. 528 were better than those for the new routine from Rev. 615 at the subsurface station with random pattern tiles, while simulated tile flow results from the new routine were better than those from the old routine at the subsurface station with constant tile spacing.Nitrate in tile flow results from the new routine from Rev. 615 were better than those from the old routine from Rev. 528 at both sites.Simulated flow and nitrate in flow results from the new routine were better than those from the old routine at the river station on the county line.The new routine provided more realistic and accurate simulation of tile drainage, and the new curve number retention parameter adjustment factor in Rev. 645 improved surface runoff simulation, and is suitable for surface runoff simulation in mildly sloped watersheds.
The results determined which tile drainage routine can provide a better model fit, and provided representative parameter sets in SWAT for simulation of tile flow, nitrate in tile flow, surface runoff, sediment and nitrate in surface runoff at field scale, simulation of streamflow, and sediment and nitrate in streamflow at watershed scale in tile-drained watersheds.The results provide guidance for selection of tile drainage routines and related parameter sets for accurate simulation of tile drainage systems for hydrologic processes at both field and watershed scales, and can be used for tile flow, runoff, and sediment and nitrate loss simulation of mildly sloped watersheds in the Midwestern US.It is necessary and important to test tile drainage routines and related parameter sets before their applications in hydrological and water quality modeling.To improve representation of the tile drainage system in the model, DEP_IMP can be separated into two parameters, depth to impervious layer and soil permeability factor.Water table depth calculation during long dry periods can also be improved.The new routine and the improved curve number calculation method can be tested for more individual tiles and watersheds.
Data availability.The observed data used in this study are not publicly accessible.These data have been collected by personnel from the Department of Agricultural Engineering, University of Illinois, at Urbana-Champaign, with funds from various cooperative sources.If anyone would like to use these data in any form in papers and/or publications, they should contact J. Kent Mitchell, Prasanta Kalita, Michael C. Hirschi, or Richard A. C. Cooke to obtain permission.They should also include one or more of the aforementioned authors responsible for these data as co-author(s).Moreover, they should also incorporate the same acknowledgments as those listed in this study, and they should also use the specific names for the stations in the publication as required by the data usage requirements.
Competing interests.The authors declare that they have no conflict of interest.
Special issue statement.This article is part of the special issue "Coupled terrestrial-aquatic approaches to watershed-scale water resource sustainability".It is not connected with a conference.

Figure 1 .
Figure 1.Cali and Vali represent calibration and validation, respectively.
were for model warm-up.Model outputs and annual corn and soybean yield from 1991 to 1997, and from 1998 to 2003, at sites B and E were compared with the observed values for model calibration and validation, respectively.Monthly tile flow and nitrate in tile flow from 1992 to 2000 and from 2002 to 2003 at site B were compared with the observed values for model calibration and validation, respectively.Monthly tile flow and nitrate in tile flow from 1991 to 2000 and from 2001 to 2002 at site E were compared with the observed values for model calibration and validation, respectively.Monthly surface runoff, sediment and nitrate in surface runoff at sites Bs and Es, and monthly flow, sediment, and nitrate in flow at site R5 from 1993 to 2001 and from 2002 to 2003 were compared with the observed values for model calibration and validation, respectively.Multi-objective autocalibration of the model for all sites was performed in Rstudio 1.0.136 with Monte Carlo simulation

Figure 2 .
Figure 2. Calibration and validation results for annual crop yields at sites B (a, b) and E (c, d).Obs and Rev. 615 represent Observed and Rev. 615, respectively.

Figure 4 .
Figure 4. Calibration and validation results for monthly surface runoff at sites Bs (a, b) and Es (c, d).Obs and Revs.615 and 645 represent Observed and revisions 615 and 645, respectively.

Figure 5 .
Figure 5. Calibration and validation results for monthly sediment losses in surface runoff at sites Bs (a, b) and Es (c, d), and monthly sediment losses in flow at site R5 (e, f).Obs and Revs.528, 615, and 645 represent Observed and Revision 528, 615, and 645, respectively.

Figure 6 .
Figure 6.Calibration and validation results for monthly nitrate-nitrogen losses in tile flow at sites B (a, b) and E (c, d), and monthly nitrate-nitrogen losses in flow at site R5 (e, f).Obs, Revs.528, and 615 represent Observed, Revision 528, and 615, respectively.

Figure 7 .
Figure 7. Calibration and validation results for monthly nitrate-nitrogen losses in surface runoff at sites Bs (a, b) and Es (c, d).Obs and Rev. 645 represent Observed and Revision 645, respectively.

Table 1 .
Monitored subsurface, surface and river stations (a), data for simulation in the LVR watershed (b), and cropping and tillage practices for sites B and E (c).

Table 2 .
Adjusted parameter values for plant growth simulation (a), and parameters used for model calibration (b).

Table 3 .
Calibrated values of adjusted parameters for tile flow and nitrate-N calibration of SWAT at sites B, E, Bs, Es, and R5.

Table 4 .
Performance evaluation of the calibrated and validated results at sites B, E, Bs, Es, and R5.