Simulating Soybean–Rice Rotation and Irrigation Strategies in Arkansas, USA Using APEX

With population growth and resource depletion, maximizing the efficiency of soybean (Glycine max [L.] Merr.) and rice (Oryza sativa L.) cropping systems is urgently needed. The goal of this study was to shed light on precise irrigation amounts and optimal agronomic practices via simulating rice–rice and soybean–rice crop rotations in the Agricultural Policy/Environmental eXtender (APEX) model. The APEX model was calibrated using observations from five fields under soybean–rice rotation in Arkansas from 2017 to 2019 and remote sensing leaf area index (LAI) values to assess modeled vegetation growth. Different irrigation practices were assessed, including conventional flooding (CVF), known as cascade, multiple inlet rice irrigation with polypipe (MIRI), and furrow irrigation (FIR). The amount of water used differed between fields, following each field’s measured or estimated input. Moreover, fields were managed with either continuous flooding (CF) or alternate wetting and drying (AWD) irrigation. Two 20-year scenarios were simulated to test yield changes: (1) between rice–rice and soybean–rice rotation and (2) under reduced irrigation amounts. After calibration with crop yield and LAI, the modeled LAI correlated to the observations with R2 values greater than 0.66, and the percent bias (PBIAS) values were within 32%. The PBIAS and percent difference for modeled versus observed yield were within 2.5% for rice and 15% for soybean. Contrary to expectation, the rice–rice and soybean–rice rotation yields were not statistically significant. The results of the reduced irrigation scenario differed by field, but reducing irrigation beyond 20% from the original amount input by the farmers significantly reduced yields in all fields, except for one field that was over-irrigated.


Introduction
Arkansas accounts for 47.9% of the rice grain output and 49% of the rice production area in the United States, making it the largest rice producer in the country [1]. Of its total rice area, 68.5% is planted in rotation with soybean as of 2018 [1]. Arkansas is also the tenth largest producer of soybean in the country, with an estimated planted area of 1,214,057 hectares [2,3]. A two-year soybean-rice rotation is a typical rotation used in Arkansas [4]. Although rice-rice systems have higher gross returns because of the high value of rice, rotations with soybean are generally more profitable due to the lower production costs of soybean [4]. Soybean is known to increase available nitrogen in the soil through nitrogen fixation [5]. Since much of the nitrogen is removed from the field during harvest, the benefits of nitrogen fixation by soybean may be minimal [6]. However, crops in rotation with soybean are still shown to increase yield, perhaps for a variety of reasons not necessarily related to nitrogen [6], such as weed and pest management [5,7]. Thus, studies focusing on comparing the long-term yields in soybean-rice compared to rice-rice rotations are needed.
In addition to crop rotation, the implementation of different irrigation methods to reduce water usage is becoming increasingly necessary. The alluvial aquifer used for agricultural irrigation in Arkansas is decreasing in connection to rice production [8]. Groundwater is used to irrigate approximately 76% of the irrigated rice area in Arkansas, whereas the remaining area uses surface water [1]. The most common rice irrigation method is delayed flooding following drill-seeding, and it generally involves flooding in the early growth stage and maintaining a depth of around 51-102 mm (2-4 inches) [9] or deeper (102-152 mm) to suppress disease [9]. Flooding uses a large amount of water, a significant portion of which may be lost as runoff and evaporation [9]. Flood irrigation typically involves either straight or contour levees, which act as spillways through which water follows the slope of the field to flow into each paddy [10]. Generally, in conventional flooding (CVF) known as cascade or single-inlet irrigation, farmers fill either the highest paddy, allowing water to flow into the lower paddies through the levee gates, or the bottom of the field, serving to "stair-step" the water back up [9]. An alternative is multiple inlet rice irrigation (MIRI), which utilizes polypipe laid along the length of the field to fill each paddy at the same time [10]. The MIRI technique can reduce 22% of the water pumped compared to CVF [11], reduce waste due to runoff, and reduce the wear on levee gates due to over-pumping [9]. The MIRI technique can be used with both straight and contour levees. As of 2018, the MIRI technique is used in around 33% of the total Arkansas rice area, while flooding with levees is used on around 56% [1]. Zero-grade is another practice used in flooded fields, which involves precision-leveling the field and flooding from multiple sides to increase the uniformity of the flood [10]. Zero-grade fields can be flooded quickly but take longer to drain, so this method is mostly used for rice-rice rotation [10]. Zero-grade land management is used on 12.8% of the Arkansas rice area and includes fields that use CVF, MIRI, and alternate wetting and drying (AWD) [1].
Furrow irrigation for rice (FIR), also called row rice, involves pumping water from the top of the field into trenches or furrows dug between beds, allowing the water to flow down the trenches using the slope of the field [12]. Rice seeds are drilled into both the furrows and the ridges [13]. Furrow irrigation is used in part to reduce water use since water is only applied to the furrows [14]. The FIR method has been used increasingly, from 2.7% (2016) to over 10% (2019) of Arkansas' total rice area [1,13]. FIR may be considered for fields that have relatively steep slopes, such that the required number of levees take up too much area or are challenging to flood [9]. While FIR has the capacity to achieve similar yields while using less water compared to flood irrigation, doing so can be difficult. Some previous studies reported a significant reduction from 10 to 40% of rice yield compared to flooded irrigation systems [13]. Studies in this area are inconclusive, and other factors, such as agronomic management, may be involved [9].
Currently, the intermittent flooding or AWD method is used on only 3% of Arkansas's total rice production area [1]. In AWD, a field is flooded and allowed to dry alternately rather than remaining flooded continuously throughout the season. This technique can reduce water use by up to 15-30% while not reducing yield [15][16][17][18][19]. According to Carrijo et al. [20], mild AWD reduced water use by 23.4% without a significant reduction in yield, and even severe AWD resulted in yield losses of 22.6% when compared to continuous flooding. Nalley et al. [18] compared different AWD regimes, the most severe of which resulted in a water use reduction of 44% and a yield reduction of 12.6% compared to continuous flooding [18]. Application of AWD irrigation can be complicated and site-specific, and recommendations on the severity of the dry periods have varied [21]. More studies can help determine optimum irrigation amounts and the risk of yield loss. Modeling increasingly severe irrigation regimes can therefore be useful for finding the limits of AWD irrigation.
Crop models are useful tools in evaluating the environmental and agricultural effects of various management techniques and environmental factors [22][23][24]. In this study, we were interested in the crop yield and how it was affected by different methods of irrigation management and crop rotation of rice-rice and soybean-rice systems under reduced and no-till systems. These questions were investigated by using the Agricultural Policy/Environmental eXtender (APEX) model [25], which is a small to medium field scale watershed model. The scope of this work involves calibration and validation of the APEX model using crop yield at production field sites and vegetation growth using the remote sensing retrieved leaf area index (LAI) from the Moderate Resolution Imaging Spectroradiometer (MODIS) [26]. Only a few studies have used MODIS products for crop model calibration, but none of those studies used APEX [27,28]. Next, we evaluated the effects of increasingly water-stressed scenarios on rice yield, water stress, and nitrogen stress of an annual soybean-rice rotation as compared to a rice-rice rotation. The goals were to (1) determine whether including soybean in the rotation had a significant effect on rice yield or nitrogen stress and (2) find the threshold of water stress before rice yields were significantly affected.

Site Description
The experimental fields in this study consist of five fields in an annual soybean-rice rotation coded as R2, R3, R5, R7, and R8. These fields are located in Faulkner (R7), Lonoke (R2, R5, R8), and Prairie (R3) counties in Central Arkansas ( Figure 1). The R7 field has Perry clay soils, is in hydrologic soil group D, and is located around 80 km from the others. The rest of the fields have Stuttgart (R2, R5), Dewitt (R3), and Calhoun (R8) soil series, are within 13 km of each other, and have silt loam soils in hydrologic soil group C/D. Table 1 summarizes the general characteristics of these fields. a small to medium field scale watershed model. The scope of this work involves calibration and validation of the APEX model using crop yield at production field sites and vegetation growth using the remote sensing retrieved leaf area index (LAI) from the Moderate Resolution Imaging Spectroradiometer (MODIS) [26]. Only a few studies have used MODIS products for crop model calibration, but none of those studies used APEX [27,28]. Next, we evaluated the effects of increasingly water-stressed scenarios on rice yield, water stress, and nitrogen stress of an annual soybean-rice rotation as compared to a rice-rice rotation. The goals were to (1) determine whether including soybean in the rotation had a significant effect on rice yield or nitrogen stress and (2) find the threshold of water stress before rice yields were significantly affected.

Site Description
The experimental fields in this study consist of five fields in an annual soybean-rice rotation coded as R2, R3, R5, R7, and R8. These fields are located in Faulkner (R7), Lonoke (R2, R5, R8), and Prairie (R3) counties in Central Arkansas ( Figure 1). The R7 field has Perry clay soils, is in hydrologic soil group D, and is located around 80 km from the others. The rest of the fields have Stuttgart (R2, R5), Dewitt (R3), and Calhoun (R8) soil series, are within 13 km of each other, and have silt loam soils in hydrologic soil group C/D. Table 1 summarizes the general characteristics of these fields.  Table 1. Soil characteristics of the experimental fields, including soil series, soil texture, hydrologic soil group (HSG), area, and slope. Soil information was retrieved from the Web Soil Survey [29]. Areas and slopes were given by farmers or found using Google Earth [30].

Management
Historically, most of these fields have been planted in soybean-rice rotation for periods ranging from 10 to more than 60 years. Rice management information was collected from the farmers, including tillage, irrigation operations, seeding and fertilizer rates, and planting and harvesting dates  Table 1. Soil characteristics of the experimental fields, including soil series, soil texture, hydrologic soil group (HSG), area, and slope. Soil information was retrieved from the Web Soil Survey [29]. Areas and slopes were given by farmers or found using Google Earth [30].

Management
Historically, most of these fields have been planted in soybean-rice rotation for periods ranging from 10 to more than 60 years. Rice management information was collected from the farmers, including tillage, irrigation operations, seeding and fertilizer rates, and planting and harvesting dates ( Table 2). Two rice varieties, both hybrids, were used: XP753 in R2, R3, and R5 and Gemini 214 CL (G214) in R7 and R8. The rice seeding rates were converted from kg ha −1 to plants ha −1 using bulk seed densities of each cultivar [31]. A germination rate of 80% was used for rice [32]. For soybean, a seeding rate of 345,947 seeds ha −1 was given by one of the farmers. We assumed an 85% germination rate, resulting in 296,526 plants ha −1 , which is reasonable considering seeding rates given by Ashlock et al. [33]. Soybean was generally planted in late May and harvested around the second week of October. No nitrogen fertilizer was applied to soybean. Each soybean field received 92 kg ha −1 of potassium and 24 kg ha −1 of phosphorus, except for R7, which received no fertilizer.

Tillage
Most fields received some sort of reduced tillage, with the exception of R7, which was a no-till field ( Table 2). Reduced tillage is defined as limiting soil disturbance to control plant residue [34]. Reduced tillage is the same as conventional tillage but involves fewer field passes. No-till is defined as no disturbance to the soil to control plant residue [35]. For the farmers in this study, no-till involves no tillage or land labor.

Irrigation
Three types of irrigation, CVF, MIRI, and FIR, were used in this study. The R3, R7, and R8 fields were under CVF, and R2 was managed with MIRI. The only field in FIR was R5. The R2, R3, and R8 fields were considered AWD, meaning the fields were allowed to dry down between irrigations at least once during the growing season rather than staying flooded throughout the season. The R7 field was flooded continuously with more than twice the amount of water than other fields. The total irrigation amounts were given by the farmers (from flowmeters or estimations from their records) or from literature values and divided by the number of irrigation events in order to obtain the weekly amount of water use for our simulation. This method was used to simulate the total amount of water applied in all fields, regardless of water management (CF or AWD). The model offers the option of furrow irrigation; however, there was no option to distinguish between conventional flooding and MIRI. Total irrigation amounts are shown in Table 2.
For soybean, the general start and end times of irrigation were given by farmers, but the amount of water applied was not. Based on information in the Arkansas Soybean Handbook [36], a value of 254 mm (10 inches) per season was used. For the soybean simulation, irrigation was applied around twice per month, with 51 mm of water per application.

APEX and WinAPEX
The APEX model was developed as an improvement on the Environmental Policy Integrated Climate (EPIC) model, which was designed to analyze the relationship between erosion and productivity of soil [37][38][39]. The EPIC model has a wide array of components and simulation capabilities and is able to compare various management systems that are limited to a homogeneous area of up to 100 ha [39]. The APEX model inherited EPIC's capabilities with water routing capacity from one subarea to another to simulate the entire watershed [39]. The APEX model has been continually developed and used worldwide [22,40]. The model has successfully simulated water use and nitrogen management [41,42], as well as a variety of best management practices [43]. It has also been used to model the effects of soil, tillage, irrigation, and cropping systems on productivity [24,[43][44][45][46][47][48]. The APEX model was chosen for this study to simulate soybean-rice rotations for a combination of future utility, ease-of-use, and comprehensiveness in terms of crop and rotation simulating abilities [41,49]. In this study, the WinAPEX interface for APEX version 1501 October 2016 [50] was used to simulate the fields.
The input weather data were downloaded from the Parameter-elevation Regressions on Independent Slopes Model (PRISM) [51] and included daily maximum and minimum temperatures and precipitation from 1981 to 2018. Thus, the Hargreaves equation [52] was used to calculate potential evapotranspiration. The watershed variables include latitude and longitude, weather station, nutrient uptake rates, and others. Within the watershed, there is a subarea component with options for editing county, soil, operation schedule, and slope, among others.
Input for the management, as discussed in Section 2.2, includes field operations, seeding rate, planting date, fertilizer rates, irrigation start and end date, irrigation amount, and harvest date. Each rice and soybean planting was created as a single annual crop in the WinAPEX interface and then combined into a soybean-rice rotation in the Microsoft Access database. The same set of crop parameters was used for rice even though, in reality, there were two different varieties. The differences in plant population were taken into account, and no calibration was performed for parameters related to weed or pest pressures, such as the parameters in the PARM Editor-PARM (9) and PARM (10)-due to the lack of damage in those areas. A 16-year warm-up run was set up, starting in 2001 in the APEX control file, to stabilize the model and to establish the soil parameters from historical land use.

APEX Calibration and Validation
The calibration of the APEX model can be carried out at the scale of a single subarea (field), a landscape, or an entire watershed [40]. In this study, each field was a standalone subarea; thus, there was no water routing between fields. The values for slope and soil series were adjusted for each subarea. The slope and soil series are shown in Table 1, along with the other field characteristics. The LAI from MODIS [26] in 2017-2019 was used to calibrate the crop growth curve and to compare with the APEX output. The MODIS observation is given as a single value covering a 500 m pixel over a 4-day period, representing the best data point between estimates from both of the MODIS satellites, Terra and Aqua. Calibration of the model was also performed by using the observed yields in 2018 for rice and 2017 for soybean. Validation was performed with soybean yields from the available data in 2019. The R7 field was fallow for both 2017 and 2019; thus, observed soybean yields were not available.
An iterative calibration procedure was performed to optimize several crop parameters. The crop parameters, including the harvest index (HI), the fraction of the growing season when leaf area declines (DLAI), the plant population for crops and grass-first point on curve (PPLP1), and plant population for crops and grass-second point on curve (PPLP2) were first adjusted one at a time to determine the effects that they would have on yield and LAI. The curve in the PPLP1 and PPLP2 referred to the plant population curve (plants m −2 versus percent of maximum leaf area index). For LAI, when the simulated LAI curve does not match the observed LAI, the first point on the optimal leaf area development curve (DLAP1) and the second point on the optimal leaf area development curve (DLAP2) can be adjusted. Changing the DLAP1 and DLAP2 changes the LAI during the initial stages of crop growth. These are the optimal (not stressed) leaf area development curves. These parameters each convey two pieces of information following the convention that numbers before the decimal represent the percentage of the growing season, while numbers after the decimal represent the fraction of the maximum potential that LAI reached at that point. In addition, the maximum potential leaf area index (DMLA) and the leaf area index decline rate (RLAD) are the two parameters that can be adjusted for yield and the maximum point on the simulated LAI.

Scenarios
The baseline scenario of this study was the management input from the field for 2017 and 2018. All the crop parameters remained the same as in the calibration process for rice and soybean, regardless of varying irrigation management methods or crop rotation. The rice-rice or soybean-rice management was used to simulate a 20-year period, from 2018 to 2038. The APEX model generates monthly weather data based on the daily historical data (1981-2018), which in turn are used to create APEX weather input data for long-term simulation [53], but it does not anticipate anthropogenic climate change. Two scenarios were performed to evaluate rice yield, focusing on (1) crop rotation and (2) irrigation reduction in rice.
In the crop rotation scenario, rice yield was compared between the soybean-rice and rice-rice rotations. Other than the rotation, the management systems for rice were identical. The hypothesis for this scenario was that the nitrogen fixation associated with soybean would increase the simulated rice yield in the long term. The rotation effects on rice yield were simulated over 20 years, from 2018 to 2038.
In the irrigation reduction scenarios, rice irrigation amounts were reduced in increments of 10% of the baseline irrigation (which varied for each field), up to 60% reduced irrigation for each field in soybean-rice rotation. The goal of this scenario was to find the threshold at which water stress has a significant impact on rice yield. The same period, 2018-2038, was chosen for consistency.

Statistical Methods
Several statistical analyses were used to assess the goodness-of-fit of the model outputs to the observed crop yields. The percent bias (PBIAS) (Equation (1)) is good for continuous, long-term simulations [54]. However, if bias is equal in both directions, it will appear as though the PBIAS is low, which may be deceptive. Therefore, Moriasi et al. [54] recommend using PBIAS, along with other statistical methods. The PBIAS for crop yield and biomass is considered acceptable at ±25%, while, for nutrient loss, it is considered acceptable at ±50% [40]. The PBIAS was also used to analyze the fit of the rice and soybean LAI crop growth curve over the years 2017-2019, but there was no standard mentioned for LAI. The percent difference (Equation (2)) and t-tests were also performed to determine whether there was a significant difference between simulated and observed yields. Since the standard deviation of the observed data was more than twice that of the predicted, a t-test assuming unequal variances was used. These analyses were performed separately for rice and soybean. For the crop rotation and reduced irrigation scenarios, paired t-tests were performed to determine differences in yield and water and nitrogen stress between values in the baseline scenario and the crop rotation or irrigation scenarios. To measure the standard deviation, the error (E) (Equation (3)), which was used for both crop yield and LAI, measures the model's tendency to over or underestimate the observed data [55]. The index of agreement (d) (Equation (4)) is also used as an indication of the degree of model error, with an optimal value of 1 and a minimum of 0 [54]. The coefficient of determination (R 2 ) (Equation (5)) measures the linearity between the observed and simulated data, and R 2 of 0.6 is considered satisfactory for any simulation with the EPIC and APEX models [40]. However, Moriasi et al. [54] suggest the field-scale R 2 values to be satisfactory in the range of (0.70-0.75) and good in the range of (0.75-0.85). The following equations were used for the statistical analyses: where OBS i and SIM i are the ith observed and simulated values, respectively; i = 1,2,3, . . . n; OBS and SIM are the mean observed and simulated values, respectively, and n is the number of samples.

Calibration and Validation
Before the calibration, the simulated yields were very low (3-5 t ha −1 for rice and 2-3 t ha −1 for soybean), so adjustments were made to the crop parameters for rice and soybean ( Table 3). The crop parameters were first adjusted one at a time to determine the effects they would have on yield. The HI and DMLA for rice were changed by less than 25%, as suggested by [56]. The PPLP1, PPLP2 for rice agreed with Le et al. [56]. In addition, the initially simulated rice LAI curve was delayed by around a month compared to MODIS. Thus, two parameters were adjusted: DLAP1 and DLAP2 [48]. The DLAP1 was changed from 30.01 to 15.05, which means that the percentage of the growing season is changed from 30 to 15% at 0.01 to 0.05 of the maximum potential LAI, respectively. For the DLAP2, the percentage of the growing season was changed from 70 to 50% at 0.95 of the maximum potential LAI. The changes made to the crop parameters in Table 3 were the same for each field. The results showed that the LAI values observed by MODIS and simulated in APEX produced consistent trends regarding vegetation growth for all the fields (Figure 2). The LAI values in APEX were zero in the areas between the cropping seasons because no off-season regrowth was simulated between the soybean and rice crops. Meanwhile, some LAI values were around 0.20 in MODIS, which might be due to crop residue or grass. All of the R 2 values met the satisfactory limit [40] and were greater than 0.70, except for the R5 field (R 2 = 0.66) ( Table 4). All d values were greater than 0.75, and E values were all within a magnitude of 0.32. The PBIAS value fluctuated and was quite high in the R7 field. This might be because of some of the outliers that appear in MODIS, although a few points were removed from the MODIS dataset if they suddenly dropped to close to zero in the middle of the growing season. The LAI range of each crop for each field is also presented in Table 4. The average LAI range of all rice fields in APEX is 0 to 4.9, and in MODIS, it is 0.1 to 4.05. Meanwhile, the average LAI range of all soybean fields in APEX is 0 to 3.28, and in MODIS, it is 0.1 to 3.15. For APEX, a daily average was used, while for MODIS, the average was used for every four days. E values were all within a magnitude of 0.32. The PBIAS value fluctuated and was quite high in the R7 field. This might be because of some of the outliers that appear in MODIS, although a few points were removed from the MODIS dataset if they suddenly dropped to close to zero in the middle of the growing season. The LAI range of each crop for each field is also presented in Table 4. The average LAI range of all rice fields in APEX is 0 to 4.9, and in MODIS, it is 0.1 to 4.05. Meanwhile, the average LAI range of all soybean fields in APEX is 0 to 3.28, and in MODIS, it is 0.1 to 3.15. For APEX, a daily average was used, while for MODIS, the average was used for every four days. Note that the APEX model did not simulate the period between the two crops, but the MODIS image provided some values greater than zero that was likely due to grass, weed, or root biomass. The calibration results of the rice yield were very good ( Table 5). The mean simulated yields in all the fields were similar to the observed, and the standard deviation was small (< 1 t ha -1 ). The t-test indicated that the observed and simulated values were not significantly different, with a p-value of 0.65 at α = 0.05. The percent difference between the mean simulated and mean observed yields was -2.47% or 0.75 t ha -1 . The PBIAS between the simulated and observed yields was −2.50, which was considered very good [54], and the E value was −0.22. The model did not consistently over or underpredict the yield.  Note that the APEX model did not simulate the period between the two crops, but the MODIS image provided some values greater than zero that was likely due to grass, weed, or root biomass. The calibration results of the rice yield were very good ( Table 5). The mean simulated yields in all the fields were similar to the observed, and the standard deviation was small (< 1 t ha −1 ). The t-test indicated that the observed and simulated values were not significantly different, with a p-value of 0.65 at α = 0.05. The percent difference between the mean simulated and mean observed yields was −2.47% or 0.75 t ha −1 . The PBIAS between the simulated and observed yields was −2.50, which was considered very good [54], and the E value was −0.22. The model did not consistently over or underpredict the yield. For the soybean calibration, three parameters, namely DLAI, maximum potential leaf area index (DMLA), and leaf area index decline rate (RLAD), were modified during calibration ( Table 3). The DMLA value at 0.8 was close to that obtained by Setiyono et al. [57] and within the maximum range found by Tagliapietra et al. [58]. The soybean calibration results (Table 6)  Both PBIAS values met the 25% satisfactory limit [40] and were considered very good according to the criteria suggested by Moriasi et al. [54]. The t-test, performed for 2017 and 2019 together, indicated that the observed and simulated soybean yields were not significantly different, with a p-value of 0.055 at α = 0.05. Note: OBS-observed yield, SIM-simulated yield, ±-standard deviation, PBIAS-percent bias, % Diff.-percent difference, E --error. "All" refers to all of the fields averaged together. R7 was fallowed in 2017 and 2019.

Crop Rotation Scenarios
According to the paired t-tests, the rice yield, water stress, and nitrogen stress days between the two crop rotation scenarios, soybean-rice and rice-rice, were not significantly different at α = 0.05. The expectation for the rice-rice scenario was that removing soybean from the rotation would be detrimental to the rice yield since soybean is known to increase soil nitrogen [5]. As can be seen in Table 7, there was no significant difference in the rice yields under the two rotation systems. Irrigation was applied on different days for each field, depending on the time of planting and farm management decisions. The AWD in R2 and R8 did not have as much water stress as the others. The AWD in R3 has the highest water stress (Table 7) because it has the least irrigation input ( Table 2). The R5 field has the fewest nitrogen stress days, which could be explained by the plant density and the amount of nitrogen input. The R8 field has the lowest nitrogen input, followed by R3, R2, R7, and R5. The R5 field has the lowest plant density, followed by R3, R2, R7, and R8. Thus, all the AWD fields (R2, R3, and R8) consistently had higher nitrogen stress days because of the high plant density and lower nitrogen inputs. The CVF field with the continuous flood (R7) had the second-highest nitrogen input as well as the second-highest plant density and thus ended up with high nitrogen stress. Note: No significant difference was found between the scenarios within the same field at α = 0.05 using the paired t-test.

Limiting Irrigation Scenario
The goal of simulating increasingly limited irrigation scenarios was to find the threshold at which the effect of water stress on yield became significant and determine whether water was being wasted in the current management strategy. The results varied for each field (Table 8). Among the three AWD fields, R2, R3, and R8 had significant yield decrease at 60, 40, and 30% reduced irrigation, respectively. The only FIR field (R5) had a significant yield decrease at 20% irrigation. The only CVF field (R7) had a significant yield increase starting at 40%. Yield decreases were generally considered significant when they reached a~5% reduction (Table 9). These scenarios provoked only minor changes in nitrogen stress days, with the only statistically significant difference in R7, starting at 40% reduced irrigation (Table A1). The water stress in R3, R5, and R8 was significant for each reduced irrigation (Table A2). The water stress in R2 and R7 became significant at 30 and 60% reduced irrigation, respectively (Table A2). Note: Superscripts of different letters across the same row denote a significant difference at α = 0.05 using the paired t-test. RI-reduced irrigation, AWD-alternate wetting and drying, FIR-furrow irrigation for rice, CVF-conventional flooding. Values of 10, 20, 30, 40, 50, and 60% are the amounts by which irrigation was reduced. "All" refers to all of the fields averaged together.  10,20,30,40,50, and 60% are the amounts by which irrigation was reduced. "All" refers to all of the fields averaged together. Percent change is calculated ((RI-baseline)/baseline) *100%).
The general effect of reduced irrigation was a reduction in yield ( Figure 3). The opposite was true, however, in R7, which consistently had increased yield compared to the baseline. Two fields, R5 and R8, both had consistently decreasing yield with reduced irrigation. The R2 and R3 fields also had general patterns of decreasing yield, with the exception that reductions in the range of 10-30% were relatively minor (<3%). The largest overall effect was in R8, possibly because of the high plant density (Table 2), increasing competition for resources between the plants. Moreover, a significant yield drop was observed in the years 2030 and 2032, which was due to less rainfall and more nitrogen stress in those two years ( Figure A1).
Sustainability 2020, 12, x FOR PEER REVIEW 11 of 17 density (Table 2), increasing competition for resources between the plants. Moreover, a significant yield drop was observed in the years 2030 and 2032, which was due to less rainfall and more nitrogen stress in those two years ( Figure A1).

Calibration and Validation
Upon running a t-test assuming unequal variances, observed yields of both rice and soybean were not significantly different from the simulated yields at α = 0.05. The model did not consistently over or underpredict rice or soybean yields, and the PBIAS was well within the 25% limit deemed satisfactory by Wang et al. [40] and was comparable to those obtained in other studies using APEX [59] and EPIC [56]. The overall E value was smaller in rice (−0.22) compared to soybean (0.38) and was lower than that obtained by Gilardelli et al. using the Water Accounting Rice Model (WARM) [60]. The percent difference between the mean of simulated and observed rice yields was −2.47%, and that of soybean yields was 11.87%, indicating very good calibration. For reference, the yields simulated by Williams et al. [38] using EPIC were all within 7% of observed yields. It is important to note that the calibration was based solely on two years for which data were available, one for rice and one for soybean. Including more years may have made for a more robust calibration. The agreement in the MODIS and APEX LAI was also a great achievement. Ren et al. [27] also used MODIS data to compare with the EPIC model at a global scale and concluded that it is feasible to use this method.

Crop Rotation Scenarios
Given that the rice-rice rotation yields were higher in three of the five simulated fields, the results of the simulation may suggest that a rice-rice rotation leads to higher yields compared to the soybean-rice rotation. According to the paired t-tests, however, none of the differences in yield were significant. Additionally, based on the simulation, there were large differences in nitrogen stress in three of the fields, although it was not consistent whether the rice-rice rotation increased or decreased stress. Again, none of the results on nitrogen stress were statistically significant. Although Anders et al. [61] conclude that including soybean in rotation with rice increases both rice yield and nitrogen uptake, Peoples et al. [6] report that the benefits of nitrogen fixation by soybean may be minimal

Calibration and Validation
Upon running a t-test assuming unequal variances, observed yields of both rice and soybean were not significantly different from the simulated yields at α = 0.05. The model did not consistently over or underpredict rice or soybean yields, and the PBIAS was well within the 25% limit deemed satisfactory by Wang et al. [40] and was comparable to those obtained in other studies using APEX [59] and EPIC [56]. The overall E value was smaller in rice (−0.22) compared to soybean (0.38) and was lower than that obtained by Gilardelli et al. using the Water Accounting Rice Model (WARM) [60]. The percent difference between the mean of simulated and observed rice yields was −2.47%, and that of soybean yields was 11.87%, indicating very good calibration. For reference, the yields simulated by Williams et al. [38] using EPIC were all within 7% of observed yields. It is important to note that the calibration was based solely on two years for which data were available, one for rice and one for soybean. Including more years may have made for a more robust calibration. The agreement in the MODIS and APEX LAI was also a great achievement. Ren et al. [27] also used MODIS data to compare with the EPIC model at a global scale and concluded that it is feasible to use this method.

Crop Rotation Scenarios
Given that the rice-rice rotation yields were higher in three of the five simulated fields, the results of the simulation may suggest that a rice-rice rotation leads to higher yields compared to the soybean-rice rotation. According to the paired t-tests, however, none of the differences in yield were significant. Additionally, based on the simulation, there were large differences in nitrogen stress in three of the fields, although it was not consistent whether the rice-rice rotation increased or decreased stress. Again, none of the results on nitrogen stress were statistically significant. Although Anders et al. [61] conclude that including soybean in rotation with rice increases both rice yield and nitrogen uptake, Peoples et al. [6] report that the benefits of nitrogen fixation by soybean may be minimal because most of the absorbed nitrogen is removed at harvest. Other studies state that including soybean in the rotation does decrease the need for nitrogen fertilizer, if only marginally [62]. Another factor in the nitrogen stress analysis is that the soybean in soybean-rice rotation did not include any nitrogen application. The modeled nitrogen fixation by soybean may have been able to make up for this deficit in some fields and not in others. Weed and pest management should also be considered. Filizadeh et al. [7] and Scherner et al. [5] both attest to the benefits of the soybean-rice rotation with regard to limiting weed growth and pest development.

Limiting Irrigation Scenarios
The results showed a decrease in yield with reduced irrigation in all fields except R7, which is the only field considered to be in continuous flood rather than AWD, as can be seen in the large discrepancy in irrigation amounts applied to those fields ( Table 2). One notable aspect of R7 is that the baseline irrigation amount (1649 mm) is almost twice that of the next highest field (751 mm); this field was over-irrigated. A typical irrigation amount for CVF rice is around 818 mm [10].
The three AWD fields, R2, R3, and R8, responded differently in the percent yield change (Table 9). Starting at RI 30%, the R8 field reduced its yield by 10% for every 10% reduction in irrigation. This might be due to the highest plant density in R8 compared to other fields. Ahmad et al. [63] assert that, indeed, density related competition between plants can affect yield and water use efficiency. The results in R8 seem to be close to those of Yang et al. [64], who found that severe AWD with 38 to 50% reduced irrigation caused a yield reduction of 19 to 35%. However, Yao et al. [65] and Rejesus et al. [16] found that using AWD with 38% reduced irrigation did not significantly reduce yield compared to continuous flood. In this study, reducing up to 20% water use only caused less than 2% yield change in all fields, except R8. Other studies have also been able to reduce irrigation by 15% [66] and 28% [67] without affecting yield. Tabbal et al. [68] even assert that, if soil saturation is maintained and weed control is adequate, irrigation could be reduced by 40% or more, compared to traditional shallow flooding, without affecting yield.
Reducing irrigation seemed to affect nitrogen stress as well, although the pattern was not consistent in each field. In some fields (R3 and R7), nitrogen stress consistently decreased with reduced irrigation, and R7 had the most nitrogen lost to denitrification as simulated by APEX compared to the other fields of the baseline irrigation. Indeed, Yang et al. [69] conclude that controlling irrigation is an effective method of reducing nitrogen loss to the environment. These effects on nitrogen stress could partly explain why the pattern of yield decrease was less clear in R3 and why the yield tended to increase with reduced irrigation in R7. The results in R3 and R7 could indicate that, in some cases, reducing irrigation not only saves water but may also increase yield. The response of rice yield to reduced irrigation is also dependent on environmental factors, such as precipitation, temperature, and soil, so results vary and direct comparisons may be misleading.
The results of R5, the only FIR field, seemed to be affected most by smaller reductions in irrigation. The R5 results could indicate that FIR is more sensitive to irrigation reduction compared to flooded rice. FIR also had lower yields compared to flooded rice in studies from Arkansas [70] and Australia [71]. Henry et al. [9] also stated that maintaining yield can be difficult with furrow irrigated rice. It is also interesting that the FIR field had so much less nitrogen stress than the other fields, probably because it had the highest amount of nitrogen input compared to other fields and less runoff.

Conclusions and Future Studies
Although an independent validation test was not performed for the rice model due to limitations in the observed data, APEX was able to model rice yields based on the calibration results of the yield and LAI. The scenarios in crop rotation showed that there was no statistically significant difference between the soybean-rice and rice-rice systems. Although the results were not significant, they could still have long-term implications concerning profits to the farmers, especially over a longer period, but a detailed economic analysis would be needed to determine whether a long-term rice-rice rotation poses any economic benefits.
The reduced irrigation scenario showed a decrease in yield with decreased irrigation, except for the R7 field, which was over-irrigated. These results raise questions about the detriments of over-irrigation, such as increased nitrogen loss, which could be examined more closely in future studies. The results in R8 also raise questions about whether high plant density can cause a negative effect on yield, which could also be examined more closely in future studies. Furthermore, more detailed future studies could test the effects of different irrigation methods on crop nitrogen stress and consider how these differences can be represented in models.