Simulating the Effects of Agricultural Management on Water Quality Dynamics in Rice Paddies for Sustainable Rice Production—Model Development and Validation

: The Agricultural Policy/Environmental eXtender (APEX) model is widely used for evaluating agricultural conservation efforts and their effects on soil and water. A key component of APEX application in Korea is simulating the water quality impacts of rice paddies because rice agriculture claims the largest cropland area in the country. In this study, a computational module called APEX-Paddy (National Academy of Agricultural Sciences, Wanju, Korea) is developed to simulate water quality with considering pertinent paddy management practices, such as puddling and ﬂood irrigation management. Data collected at two experimental paddy sites in Korea were used to calibrate and validate the model. Results indicate that APEX-Paddy performs well in predicting runoff discharge rate and nitrogen yield while the original APEX highly overestimates runoff rates and nitrogen yields on large storm events. With APEX-Paddy, simulated and observed ﬂow and mineral nitrogen yield (QN) are found to be highly correlated after calibration (Nash & Sutcliffe Efﬁciency ( NSE ) = 0.87 and Percent Bias ( PBIAS ) = − 14.6% for ﬂow; NSE = 0.68 and PBIAS = 2.1% for QN). Consequently, the APEX-Paddy showed a greater accuracy in ﬂow and QN prediction than the original APEX modeling practice using the SCS-CN (Soil Conservation Service-Curve Number) method. potential evapotranspiration (PET). provides a set of new subroutines simulate ponding, control, AET possibly than under managements. a transpiration while Daily AET values of paddy


Introduction
Rice is traditionally the most preferred grain crop in Korea. Even though the preference for rice has noticeably diminished recently with industrialization and the introduction of the western diet, rice remains as the predominant grain crop of the country. The total area of paddy fields in Korea is estimated to be 55% of the cultivated lands, claiming over one million hectares [1,2]. Rice cultivation starts in the spring season with the transplantation of rice seedlings and continues through the fall season. Paddies are kept flooded during most of the growing season (spring through summer) with intensive irrigation, fertilization, and weeding. Seasonally, about 25 to 50% of the annual rainfall occurs during the monsoon season (June and July) in Korea. Because the growing season for rice includes the monsoon period, the paddy fields are vulnerable to discharging sediment, nutrients, and other chemicals to downstream water bodies.
Rice production practices (e.g., improper application of chemicals, over-pumping of groundwater, and loss of biodiversity in rice paddies) are linked to water contamination, excessive salinity, and other environmental problems [3]. Intense fertilization during the rice cultivation period, especially in early

Description of Study Site
Two experimental paddy fields operated by the National Institute of Agricultural Science were selected for case study areas as these paddy fields are rich in monitoring data, such as water balance, management records, weather data, and yield amounts [18,19]. The Icheon site is located in Gyeonggi-do province and the Gimje site is in Jeollabuk-do province. The Icheon (37 • 18 20.34 N, 127 • 30 40.46 E) and Gimje (35 • 44 56.13 N, 126 • 51 50.21 E) sites are both irrigated paddy fields representing interior and costal characteristics of the paddy environment for rice production in Korea (Oryza sativa L.) with cultivation areas of 15 and 0.5 ha, respectively ( Figure 1). The majority of excessive water irrigated to paddy fields discharges into drainage canals. The two sites' landscapes are commonly defined by alluvial plains with poor soil drainage. The predominant soils are Seogcheon series at the Icheon site and Jeonbug series at the Gimje site. The topsoil of Seogcheon series is coarse loamy, which is somewhat poorly drained. The Jeonbug series is poorly drained silty soil built on the fluvio-marine plain. The average annual temperature is 12.8 • C and the average annual rainfall is 1024 mm. Seasonally, rainfall at the study sites is highly influenced by monsoon, a rainy season during the months of July and August in Korea. Approximately 45% of the annual rainfall occurs in the study area during this time. Because the growing season for rice in part includes the monsoon period, paddy fields are vulnerable to discharging sediment and nutrients to downstream water bodies.
Water 2017, 9, 869 3 of 16 Seogcheon series at the Icheon site and Jeonbug series at the Gimje site. The topsoil of Seogcheon series is coarse loamy, which is somewhat poorly drained. The Jeonbug series is poorly drained silty soil built on the fluvio-marine plain. The average annual temperature is 12.8 °C and the average annual rainfall is 1024 mm. Seasonally, rainfall at the study sites is highly influenced by monsoon, a rainy season during the months of July and August in Korea. Approximately 45% of the annual rainfall occurs in the study area during this time. Because the growing season for rice in part includes the monsoon period, paddy fields are vulnerable to discharging sediment and nutrients to downstream water bodies. Rice is seeded and grown in nursery beds for 25-30 days and then transplanted after 1.5-2 leaves emerge. Transplanting is the most popular method of rice cultivation, and is undertaken in the paddy fields between mid-May and mid-June. The period of land preparation (i.e., puddling) and transplanting occur from May to August, during which water demand reaches the peak rate. In paddy cultivation, paddy fields are irrigated sufficiently to maintain water ponding. The depth of ponding water is managed at various levels between 5 and 10 cm during the growing season to control weeds and to maximize rice yield. Typically, the rice is harvested in October. Continuous monitoring of irrigation and surface discharge water at the paddy field were conducted from May to September in 2002 and 2003 in the Icheon area, and from May to September in 2014, at the Gimje site, for the purpose of quantifying and qualifying water and nutrients.

APEX-Paddy Model
APEX is a watershed model for simulating the impact of agricultural managements to water balance, nutrient cycling, carbon dynamics in soil-plant systems that runs on a daily time step [16,20]. Rice is seeded and grown in nursery beds for 25-30 days and then transplanted after 1.5-2 leaves emerge. Transplanting is the most popular method of rice cultivation, and is undertaken in the paddy fields between mid-May and mid-June. The period of land preparation (i.e., puddling) and transplanting occur from May to August, during which water demand reaches the peak rate. In paddy cultivation, paddy fields are irrigated sufficiently to maintain water ponding. The depth of ponding water is managed at various levels between 5 and 10 cm during the growing season to control weeds and to maximize rice yield. Typically, the rice is harvested in October. Continuous monitoring of irrigation and surface discharge water at the paddy field were conducted from May to September in 2002 and 2003 in the Icheon area, and from May to September in 2014, at the Gimje site, for the purpose of quantifying and qualifying water and nutrients.

APEX-Paddy Model
APEX is a watershed model for simulating the impact of agricultural managements to water balance, nutrient cycling, carbon dynamics in soil-plant systems that runs on a daily time step [16,20].
APEX defines land areas with unique soil property, vegetation or land use, and the land slope as subareas. Therefore, a watershed is delineated into a system of subareas that are hydrologically connected by a stream network.
In this work a set of subroutines was developed and coded into the current APEX (version 1501) (Texas A&M AgriLife Research, Temple, TX, USA, https://epicapex.tamu.edu/apex/) to assess the impact of paddy management practices on rice growth, water balance, and water quality. This paddy sub-model is named APEX-Paddy. The greatest challenge in implementing APEX for paddy simulation lies in the fact that APEX Subareas are not coded to permit ponding conditions by flooding water. Therefore, the subarea module in APEX was enhanced to accommodate water ponding conditions with diking and outlet controls. Puddling is a tillage operation with shallow ponding water (saturating the soil is) in a paddy field. The top soil layer in the paddy becomes soft and suitable for transplanting rice seedlings after the puddling.
The standard spacing of rice in Korea is 14 cm × 30 cm. Thus, 23.8 plants per square meters are planted with three to four seedlings per plant. As rice grows after transplanting, the paddy rice goes into tillering as the fourth leaf is fully emerged. Tillers emerging from nodes in the main culm have the potential to produce panicles. APEX's crop growth module does not simulate tillering for rice. In crops such as rice, wheat and sugarcane that produce higher numbers of yielding tillers compared to the number of seeds or shoots planted, the plant population must be estimated based on the number of tillers producing the final yield. We used the following steps to calculate the rice population in the APEX model. According to the Korea Crop Test Report [21], the average effective tillering is 18.4.  Table 1. APEX-Paddy contains physically based algorithms to mechanically simulate these paddy operations and their effects on rice growth and nutrient cycling. The schematic of the APEX processes with the paddy module is illustrated in Figure 2. With the paddy enhancement, a subarea is flexible for setting land as dry or wet (flooded) as prescribed in the operation schedule file that is attached to the subarea. APEX estimates surface runoff discharge and infiltration using the SCS-CN method [22], and soil erosion is estimated by the USLE (Universal Soil Loss Equation) equation [23] or other variations of the USLE. Daily actual evapotranspiration (AET) is calculated as the sum of plant transpiration and soil evaporation, which is equal or less than the estimated daily potential evapotranspiration (PET). APEX-paddy provides a set of new subroutines to simulate water ponding, discharge control, AET possibly greater than PET under wet condition, and paddy managements. During the rice period, the development of leaf areas promotes a greater transpiration while reducing the exposure of the soil surface to the incoming solar radiation.
Daily AET values of rice paddy can often exceed the estimated daily PET during the mid-season after the canopy is established due to high solar radiation, wind speed, and root water uptake. Vu et al. [24] found AET values up to 50% greater than PET values during mid-seasons depending on varieties and regional conditions. When a subarea is set to simulate paddy rice management, APEX-Paddy switches from the SCS-CN method to a weir discharge function to control water balance in the paddy. Sediment yield is estimated based on the settling rate after puddling and then residual sediment concentration in the ponding water. Daily AET may exceed PET during the summer season after rice establishes a full canopy. Infiltration of ponding water occurs continuously based on the saturated hydraulic conductivity of the top soil layer. During off-seasons or when paddy field management does not implement water ponding, APEX-Paddy switches back to default subarea modules to simulate upland non-ponding land processes, such as the SCS-CN method for runoff estimation. A puddling operation results in a rapid resuspension of sediment and nutrient, thus making sediment and nutrient concentration high in the ponding water before the suspended solids resettle. Growth characteristics of paddy rice were developed and calibrated based on field measurement data at the Suwon site in Gyeonggi-do. Leaf area index, dry weight and root weight ratio were measured.
Water 2017, 9,869 5 of 16 estimated daily potential evapotranspiration (PET). APEX-paddy provides a set of new subroutines to simulate water ponding, discharge control, AET possibly greater than PET under wet condition, and paddy managements. During the rice period, the development of leaf areas promotes a greater transpiration while reducing the exposure of the soil surface to the incoming solar radiation. Daily AET values of rice paddy can often exceed the estimated daily PET during the mid-season after the canopy is established due to high solar radiation, wind speed, and root water uptake. Vu et al. [24] found AET values up to 50% greater than PET values during mid-seasons depending on varieties and regional conditions. When a subarea is set to simulate paddy rice management, APEX-Paddy switches from the SCS-CN method to a weir discharge function to control water balance in the paddy. Sediment yield is estimated based on the settling rate after puddling and then residual sediment concentration in the ponding water. Daily AET may exceed PET during the summer season after rice establishes a full canopy. Infiltration of ponding water occurs continuously based on the saturated hydraulic conductivity of the top soil layer. During off-seasons or when paddy field management does not implement water ponding, APEX-Paddy switches back to default subarea modules to simulate upland non-ponding land processes, such as the SCS-CN method for runoff estimation. A puddling operation results in a rapid resuspension of sediment and nutrient, thus making sediment and nutrient concentration high in the ponding water before the suspended solids resettle. Growth characteristics of paddy rice were developed and calibrated based on field measurement data at the Suwon site in Gyeonggi-do. Leaf area index, dry weight and root weight ratio were measured.

Evapotranspiration
The APEX-Paddy model uses the equations of Penman-Monteith [25] and Stockle et al. [26] for calculating plant transpiration. For paddy simulation, the actual daily evapotranspiration method as suggested by Sakaguchi's [27], is incorporated to represent crop coverage and flood condition.

Evapotranspiration
The APEX-Paddy model uses the equations of Penman-Monteith [25] and Stockle et al. [26] for calculating plant transpiration. For paddy simulation, the actual daily evapotranspiration method as suggested by Sakaguchi's [27], is incorporated to represent crop coverage and flood condition.
The amount of evaporation and plant transpiration from the soil and the leaf surface are calculated using following equations: where, λ is the latent heat of vaporization (MJ(Mega Joule)/kg), E 0 is daily potential evapotranspiration (mm/d) for the reference crop, alfalfa with a 40 cm height, EP is daily crop transpiration (mm/d), ∆ represents the slope of the saturation vapor pressure temperature relationship (kPa/ • C), and R n is net solar radiation (MJ/m 2 /d). ρ a is the mean air density at constant pressure (kg/m 3 ), VPD represents the vapor pressure deficit of the air (kPa), u 10 is daily mean wind speed from a point 10 m above the ground, is the psychrometric constant (kPa/ • C), r c is the plant canopy resistance (s/m), and r a is the diffusion resistance of the air layer (s/m). Canopy resistance is estimated by dividing the minimum surface resistance for a single leaf by one-half of the canopy leaf area index [28]. The aerodynamic resistance to sensible heat and vapor transfer, r a , is calculated based on wind velocity and surface roughness factors following Williams et al. [29]. The amount of evaporation from the water surface when a rice paddy is flooded is calculated using Sakaguchi's equations [27]: where, V evap is daily evaporation (mm/d), EVAP is daily evapotranspiration (mm/d), LAI is leaf area index, and LAI evap is leaf area index at which evaporation from water surface does not occur. η refers to the coefficient of evaporation from the water surface. The default coefficient value was set at 0.6, which is applied to the reservoir module of the Soil and Water Assessment Tool (SWAT) model [30]. LAI evap was set at 4.0, based on the work of Miyazaki et al. [31], which was modified by Sakaguchi et al. [27].

Puddling Simulation
In East Asian countries including Korea, rice cultivation often starts with transplantation. Most other crop models can only simulate the growth of rice in flooded paddy fields, but are incapable of simulating the environmental effects of paddy managements such as puddling operations on discharge water quality. Puddling is performed prior to transplanting as a part of field preparation, to break clods and to flatten the rice paddy for transplantation. Somura et al. [32] reported that the discharge of pollutants (e.g., soil, nitrogen, and phosphorous) could significantly increase during (or shortly after) a puddling practice in Japan. APEX-Paddy was designed to simulate sediment resuspension with a puddling operation. Sediment concentration in effluent water is calculated using a modified Stokes equation and residual concentration of sediment [33]: where, C sed,f is final soil concentration in the water body (mg/L), C sed,i is soil concentration in the water body during day one of puddling (mg/L), t is time of occurrence (day), and C sed,rsd is concentration of remaining soil in water body (mg/L). d 50 refers to the diameter of the soil particles (µm), and F clay , F silt , F sand , respectively, refer to the ratios of clay, silt, and sand grain sizes in the surface soil. In conjunction with paddy water balance and nutrient transport modules, discharge of nitrogen and phosphorous can be calculated based on fertilizer management, irrigation schedule, and daily weather conditions.

Transplanting Simulation
While other upland crops are directly seeded on cultivated land, paddy rice is seeded in seedbeds and grows for 25-30 days in a nursery before the seedlings are transplanted into a puddled field. The rice transplanting promotes higher yields and less weeding. The Leaf Area Index (LAI), during the emergence and falling of leaves, was calculated using the following equation after Williams et al. [29]: In Equation (8), LAI 0,i and LAI i are the daily initial LAI, and final LAI, of a crop (i), respectively, XLAI i is the maximum LAI, TLAI is the total LAI during a growth period, REG i is the stress factor of the crop, and ∆HUF i is the daily amount of change of the heat unit. In Equation (9), HUI i refers to the heat unit index of the crop, HUF i refers to the heat unit factor of the crop, and α and β refer to coefficients related to the growth characteristics of the crop. A transplantation operation added to APEX-Paddy was designed to simulate the start of the crop growth with non-zero LAI, biomass amount, and controlled the plant population, which otherwise are forced to start from zero at seed-planting in APEX version 1501(APEX1501).

Model Calibration and Validation
For model calibration, the APEX-auto Calibration and UncerTainty Estimator (APEX-CUTE) programs (3.0 version, Texas A&M AgriLife Research, Temple, TX, USA, https://epicapex.tamu.edu/ apex/) were used. APEX-CUTE is a calibration method, which uses a dynamically dimensioned search (DDS) algorithm [34]. The DDS algorithm drives the model by sequentially applying candidate values as input data until the maximum value of the objective function is derived, and then calculates the actual statistics for the model. Random changes are applied to the optimized value to produce another candidate value, which is applied to the model. APEX-CUTE is programmed to repeat the above process for as many times as the testing of the predetermined objective function.
The Nash-Sutcliffe Efficiency (NSE) [30], Percent Bias (PBIAS), and coefficient of determination (R 2 ) were used as the objective functions for model calibration. The NSE indicates that the means of field data, rather than the simulated values, should be used for more accurate results, when the derived value ranges from 0 to 1.0. The closer the simulated values and field data are, the closer the derived value is to 1.0 [35,36]. The NSE was calculated using Equation (10): where, Y obs i and Y sim i are observed and simulated data at time i, respectively, Y mean is the mean of observed data for the constituent being evaluated, and n is the total number of observations [37]. PBIAS evaluates central trends of simulated output by estimating the sum of residuals, the difference in data points between the observed and predicted normalized by the sum of observed data [38]. The optimal value for PBIAS is 0, with low-magnitude values indicating accurate model simulation [37]. Positive and negative values indicate model underestimation and overestimation bias, respectively.
Moriasi et al. [37] suggested the criteria for evaluation (Table 2) based on the performance value of the watershed model. The simulation result of the model in this study was evaluated based on these criteria.
Note: a Includes stream flow, surface runoff, and base flow for watershed-and field-scale models.

Characterization of Paddy Rice for APEX Simulation
In the APEX-Paddy model, the LAI value of rice at transplanting was estimated based on the growing degree units of seedlings. The estimated LAI of 0.1 was then used as the initial value in the APEX simulation as a part of the transplanting operation. Additionally, the heat unit factor of crop (HUF) during the cultivation period was reduced by the amount of heat units accumulated during the nursery period. The growth characteristic curve according to Equations (8) and (9) for transplanted paddy rice was calibrated to field data, as the default coefficient values reflect seed-planted upland rice ( Figure 3). When compared to upland rice (default in APEX), paddy rice showed a faster growth rate in early growing season. As depicted in Figure 3, the LAI curve for paddy rice diverts from upland rice at 30% of the growing season and establishes the full canopy before reaching 80% of the maturity.
Note: a Includes stream flow, surface runoff, and base flow for watershed-and field-scale models.

Characterization of Paddy Rice for APEX Simulation
In the APEX-Paddy model, the LAI value of rice at transplanting was estimated based on the growing degree units of seedlings. The estimated LAI of 0.1 was then used as the initial value in the APEX simulation as a part of the transplanting operation. Additionally, the heat unit factor of crop (HUF) during the cultivation period was reduced by the amount of heat units accumulated during the nursery period. The growth characteristic curve according to Equations (8) and (9) for transplanted paddy rice was calibrated to field data, as the default coefficient values reflect seedplanted upland rice (Figure 3). When compared to upland rice (default in APEX), paddy rice showed a faster growth rate in early growing season. As depicted in Figure 3, the LAI curve for paddy rice diverts from upland rice at 30% of the growing season and establishes the full canopy before reaching 80% of the maturity. The weight of the rice root during the initial cultivating period after transplantation was 35% of the total biomass and the growth of roots almost ceased after the heading season. Consequently, the root weight proportion of the whole plant decreased to 16% through the heading season, and to 7% by the harvesting season. Using these values, the root growth characteristics of paddy rice were refined for the model. The root-shoot weight ratio is calculated using Equation (12): where, RWi is the total root weight of the crop (i) (tons/ha), DMi is the total biomass of the crop (i) (tons/ha), and ε and γ refer to the coefficient of the crop. Important crop parameters for simulating paddy rice compared with upland rice are summarized in Table 3. Detailed descriptions of parameters are in the work of Stenglich and Williams [39]. The weight of the rice root during the initial cultivating period after transplantation was 35% of the total biomass and the growth of roots almost ceased after the heading season. Consequently, the root weight proportion of the whole plant decreased to 16% through the heading season, and to 7% by the harvesting season. Using these values, the root growth characteristics of paddy rice were refined for the model. The root-shoot weight ratio is calculated using Equation (12): where, RW i is the total root weight of the crop (i) (tons/ha), DM i is the total biomass of the crop (i) (tons/ha), and ε and γ refer to the coefficient of the crop. Important crop parameters for simulating paddy rice compared with upland rice are summarized in Table 3. Detailed descriptions of parameters are in the work of Stenglich and Williams [39].

Effects of Paddy Management on Water and Nitrogen Balance
The Icheon paddy model was calibrated based on daily discharge rates that were measured at the outlet of the paddy field from May to September of 2002. R 2 was used as the main objective function for calibrating the model to extract the optimum parameter value. As depicted in Figure 4, the performance of the paddy model is demonstrated to be excellent in predicting runoff discharge rate during the calibration period (R 2 = 0.88, NSE = 0.87, PBIAS = −14.57%), which proves that the simulation effectively reflects the field data. The calibrated Icheon site model was evaluated for runoff estimation during the cropping season (June to September) of 2003. The performance statistics indicated good performance (R 2 = 0.80, NSE = 0.65, PBIAS = 9.6%), which also supports the accurate simulation of field data.

Effects of Paddy Management on Water and Nitrogen Balance
The Icheon paddy model was calibrated based on daily discharge rates that were measured at the outlet of the paddy field from May to September of 2002. R 2 was used as the main objective function for calibrating the model to extract the optimum parameter value. As depicted in Figure 4, the performance of the paddy model is demonstrated to be excellent in predicting runoff discharge rate during the calibration period (R 2 = 0.88, NSE = 0.87, PBIAS = −14.57%), which proves that the simulation effectively reflects the field data. The calibrated Icheon site model was evaluated for runoff estimation during the cropping season (June to September) of 2003. The performance statistics indicated good performance (R 2 = 0.80, NSE = 0.65, PBIAS = 9.6%), which also supports the accurate simulation of field data.     Nitrogen (N) fertilizer is the main source of nitrogen yield in paddy discharge water. As summarized in Table 1, typically a total of 100-120 kg N/ha of fertilizer is applied during a growing season. As the N fertilizer is applied in mineral form, poorly managed paddy irrigation and discharge can make significant contribution to QN loads in the downstream of paddy fields. Only the Icheon site was monitored for QN during 2002 and 2003 and data was available for model evaluation. At the monitoring site, the two important factors that affect QN were the runoff discharge rate and the timing of fertilizer application. As depicted in Figure 6, the largest discharge event in 2002 resulted in the greatest daily QN for the year. Furthermore, the hikes of QN were often initiated by fertilizer applications. In general, the model calibrated for runoff discharge successfully simulated the timing and peak rates of QN when compared to observed values (R 2 = 0.66, NSE = 0.63, PBIAS = 2.1%). In the validation year (2003), the performance statistics for QN were measured well beyond acceptable ranges (R 2 = 0.64, NSE = 0.43, PBIAS = 4.5%). The predicted values satisfy the criteria in evaluating the model's performance recommended by Moriasi et al. [37].  Nitrogen (N) fertilizer is the main source of nitrogen yield in paddy discharge water. As summarized in Table 1, typically a total of 100-120 kg N/ha of fertilizer is applied during a growing season. As the N fertilizer is applied in mineral form, poorly managed paddy irrigation and discharge can make significant contribution to QN loads in the downstream of paddy fields. Only the Icheon site was monitored for QN during 2002 and 2003 and data was available for model evaluation. At the monitoring site, the two important factors that affect QN were the runoff discharge rate and the timing of fertilizer application. As depicted in Figure 6, the largest discharge event in 2002 resulted in the greatest daily QN for the year. Furthermore, the hikes of QN were often initiated by fertilizer applications. In general, the model calibrated for runoff discharge successfully simulated the timing and peak rates of QN when compared to observed values (R 2 = 0.66, NSE = 0.63, PBIAS = 2.1%). In the validation year (2003), the performance statistics for QN were measured well beyond acceptable ranges (R 2 = 0.64, NSE = 0.43, PBIAS = 4.5%). The predicted values satisfy the criteria in evaluating the model's performance recommended by Moriasi et al. [37]. Nitrogen (N) fertilizer is the main source of nitrogen yield in paddy discharge water. As summarized in Table 1, typically a total of 100-120 kg N/ha of fertilizer is applied during a growing season. As the N fertilizer is applied in mineral form, poorly managed paddy irrigation and discharge can make significant contribution to QN loads in the downstream of paddy fields. Only the Icheon site was monitored for QN during 2002 and 2003 and data was available for model evaluation. At the monitoring site, the two important factors that affect QN were the runoff discharge rate and the timing of fertilizer application. As depicted in Figure 6, the largest discharge event in 2002 resulted in the greatest daily QN for the year. Furthermore, the hikes of QN were often initiated by fertilizer applications. In general, the model calibrated for runoff discharge successfully simulated the timing and peak rates of QN when compared to observed values (R 2 = 0.66, NSE = 0.63, PBIAS = 2.1%). In the validation year (2003), the performance statistics for QN were measured well beyond acceptable ranges (R 2 = 0.64, NSE = 0.43, PBIAS = 4.5%). The predicted values satisfy the criteria in evaluating the model's performance recommended by Moriasi et al. [37].

Performance of APEX-Paddy over APEX1501
A prevailing advantage of the proposed APEX-Paddy algorithm for simulating hydrological processes over the conventional SCS-CN method [22] is the improved simulation of water management and of related hydrologic processes calculation specific to paddy environment. As accurate prediction of water balance is essential to estimating water quality variables such as sediment yield and nutrient yield, the improved water balance estimation provides a synergetic effect in predicting water quality output.
For comparison, a baseline simulation was conducted for the Icheon field, in which runoff was estimated using the SCS-CN method (APEX1501) and was compared with the APEX-Paddy. The uncalibrated Icheon dataset was rerun by APEX1501 model and calibrated for runoff and QN using the APEX-CUTE program [40]. A total of 10,000 iterations were made using the DDS algorithm to calibrate 30 APEX parameters. Outputs from both APEX-Paddy and APEX1501 were collected and compared for relative performance estimation. Results indicate that APEX1501 significantly overestimates peak rates of runoff and QN.
The best result of the SCS-CN method predicted the biggest runoff event in July 2002, which was greater than observed runoff rate by twice the magnitude. Table 4 presents the model calibration and validation results related to daily runoff. Overall the simulated runoff by APEX1501 was unsatisfactory as evidenced by performance statistics: R 2 = 0.57, NSE = −1.91, PBIAS = −80.9%. The negative NSE value indicates that there exists a notable difference between simulated runoff and the mean observed value while the predicted runoff has meaningful correlation with observed data. In particular, the largest storm event in August 2002 was highly overestimated for peak discharge estimation with the APEX1501 estimation because the paddy defined with the model had no outlet control and therefore had zero surface storage capacity to retain the stormwater (see Figure 7). The high peak discharge estimation is well reflected in the APEX1501 output on other small to intermediate storm events. The negative PBIAS value implies that predicted runoff was highly overestimated. According to Moriasi et al. [37], the performance of the APEX1501 is unsatisfactory.

Performance of APEX-Paddy over APEX1501
A prevailing advantage of the proposed APEX-Paddy algorithm for simulating hydrological processes over the conventional SCS-CN method [22] is the improved simulation of water management and of related hydrologic processes calculation specific to paddy environment. As accurate prediction of water balance is essential to estimating water quality variables such as sediment yield and nutrient yield, the improved water balance estimation provides a synergetic effect in predicting water quality output.
For comparison, a baseline simulation was conducted for the Icheon field, in which runoff was estimated using the SCS-CN method (APEX1501) and was compared with the APEX-Paddy. The uncalibrated Icheon dataset was rerun by APEX1501 model and calibrated for runoff and QN using the APEX-CUTE program [40]. A total of 10,000 iterations were made using the DDS algorithm to calibrate 30 APEX parameters. Outputs from both APEX-Paddy and APEX1501 were collected and compared for relative performance estimation. Results indicate that APEX1501 significantly overestimates peak rates of runoff and QN.
The best result of the SCS-CN method predicted the biggest runoff event in July 2002, which was greater than observed runoff rate by twice the magnitude. Table 4 presents the model calibration and validation results related to daily runoff. Overall the simulated runoff by APEX1501 was unsatisfactory as evidenced by performance statistics: R 2 = 0.57, NSE = −1.91, PBIAS = −80.9%. The negative NSE value indicates that there exists a notable difference between simulated runoff and the mean observed value while the predicted runoff has meaningful correlation with observed data. In particular, the largest storm event in August 2002 was highly overestimated for peak discharge estimation with the APEX1501 estimation because the paddy defined with the model had no outlet control and therefore had zero surface storage capacity to retain the stormwater (see Figure 7). The high peak discharge estimation is well reflected in the APEX1501 output on other small to intermediate storm events. The negative PBIAS value implies that predicted runoff was highly overestimated. According to Moriasi et al. [37], the performance of the APEX1501 is unsatisfactory.  Predicted daily discharge hydrographs indicate that APEX-Paddy successfully simulates high peaks in daily discharge at the Icheon site by allowing surface retention of stormwater. The APEX1501 model with the SCS-CN method over-predicted peak daily discharges up to 250% of the measured value.   In reality, even a best management practice (BMP) for controlling the Non-Point Sources (NPS) pollution may not be effective if the BMP results in meaningful a decrease in crop yield. This means that recommendations on controlling NPS pollution must be developed in a way that minimizes influencing crop production and promotes economic sustainability for farmers. Table 6 presents simulated and statistical data (1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014) of rice production at the Icheon and the Gimje site [41]. The predicted rice yield by APEX-Paddy compares well with observed values (PBIAS = 4.0%) for the Icheon field. Similarly, APEX-Paddy performed well for the Gimje field in predicting rice yield (PBIAS = 0.75%), which is in part attributable to the successful improvement in crop parameterization for root growth and LAI development conducted as part of the modeling study.

Conclusions
With agricultural management including fertilizer application and irrigation, paddy fields are potentially a source of non-point pollution. This study found that agricultural management significantly influences the discharge water quality of rice paddy fields. The ability to simulate paddy management processes using physically based algorithms had to be incorporated to the APEX model to successfully calibrate flow and nitrogen yield. The Agricultural Policy/Environmental eXtender model was modified for simulating biophysical processes in paddy fields under flooded (and dry) conditions. Agricultural management for growing rice was prescribed in scheduled operations, and relevant biophysical processes were constructed in a paddy algorithm.
The ability to assess the effects of paddy management on water quantity and quality in simulation models allows for a better understanding of the interactive processes of management, water, and pollutant loads. Fertilizer management can also improve discharge water quality while maintaining rice production. Two experimental paddy fields, Icheon and Gimje in Korea, were selected for evaluating the paddy module in the APEX model. As demonstrated in the case study at the two paddy fields in Korea, incorporating physical processes of paddy management into APEX resulted in substantial improvement in predicting water balance and nitrogen yield. Refined characterization of paddy rice growth might have improved water quality output as well as rice yield prediction. However, nutrient dynamics may need further improvement to better represent microbial effects in standing water and in the anoxic zone of the subsoil layers.
In Korea, paddies are often located in the vicinity of large waterbodies such as rivers or lakes, in order to meet the substantial irrigation demand. Intense agricultural management of paddy rice including fertilizer and other chemical applications has been identified as possible sources of water pollution. The generally short travel time for paddy discharge to receiving a waterbody could cause significant water quality impacts during cultivation because of puddling, transplanting, and other scheduled management practices. These environmental effects can be amplified in places where paddies are predominant land uses because it is common that such scheduled paddy managements are near-concurrently operated. The proposed APEX-Paddy model showed a significant improvement in predicting the environmental processes over the SCS-CN method. Though beyond the scope of this paper, the ultimate goal of the model development is to build a reliable simulation tool for developing innovative conservation practices for agricultural NPS pollution that can account for local climate condition, soil properties, and agricultural management specifics. As a follow-up, APEX-Paddy will continue to be applied to demonstrate the effectiveness of BMPs implementation for paddy discharge water in Korea at the national scale, and the findings will be incorporated in developing agricultural policies.