Rotating maize reduces the risk and rate of nitrate leaching

There is a strong link between nitrate (NO3-N) leaching from fertilized annual crops and the rate of nitrogen (N) fertilizer input. However, this leaching-fertilizer relationship is poorly understood and the degree to which soil type, weather, and cropping system influence it is largely unknown. We calibrated the Agricultural Production Systems sIMulator process-based cropping system model using 56 site-years of data sourced from eight field studies across six states in the U.S. Midwest that monitored NO3-N leaching from artificial subsurface drainage in two cropping systems: continuous maize and two-year rotation of maize followed by unfertilized soybean (maize-soybean rotation). We then ran a factorial simulation experiment and fit statistical models to the leaching-fertilizer response. A bi-linear model provided the best fit to the relationship between N fertilizer rate (kg ha−1) and NO3-N leaching load (kg ha−1) (from one year of continuous maize or summed over the two-year maize-soybean rotation). We found that the cropping system dictated the slopes and breakpoint (the point at which the leaching rate changes) of the model, but the site and year determined the intercept i.e. the magnitude of the leaching. In both cropping systems, the rate of NO3-N leaching increased at an N fertilizer rate higher than the N rate needed to optimize the leaching load per kg grain produced. Above the model breakpoint, the rate of NO3-N leaching per kg N fertilizer input was 300% greater than the rate below the breakpoint in the two-year maize-soybean rotation and 650% greater in continuous maize. Moreover, the model breakpoint occurred at only 16% above the average agronomic optimum N rate (AONR) in continuous maize, but 66% above the AONR in the maize-soybean rotation. Rotating maize with soybean, therefore, allows for a greater environmental buffer than continuous maize with regard to the impact of overfertilization on NO3-N leaching.


Introduction
Globally, nitrogen (N) fertilization has been a crucial variable in increasing crop yields. This increase, however, comes at a high environmental cost, considering that around 15% of the N fertilizer applied to maize (Zea mays L.) leaches into the groundwater as nitrate (NO 3 -N) [1][2][3][4][5].
There is inherent risk to the environment in applying N fertilizer at any rate as the ecosystem processes that determine the fate of N fertilizer are both complex and dynamic. Within a given field and season, however, there is a threshold N rate or 'breakpoint' above which the risk increases substantially [1,[6][7][8][9][10][11][12][13]. The factors that determine the N rate at which that breakpoint occurs are not well understood as the NO 3 -N leached each year is not only derived from that year's N inputs, but also from soil organic matter as well as residual N fertilizer from previous applications and returned crop residue from previous seasons [14,15].
Nevertheless, it has been widely conjectured that the breakpoint occurs at or around the yieldoptimizing N rate (known as the agronomic optimal N rate or AONR) [1,6,12,16]. The ecosystem dynamics that determine the AONR are similar to those of the breakpoint in that they are the result of complex interactions with the weather, soil, and residue from the previously planted crop [17,18]. Therefore, as both the breakpoint and the AONR are products of complex systems, their relationship may not be as simple as previously conceived.
The aim of this study was to determine how the response of NO 3 -N leaching to increases in the N fertilizer rate relates to the AONR. This relationship can both define the environmental costs of over-fertilization and provide targeted guidance to improve management strategies.
The soil N storage capacity provides a proportional buffer between the AONR and the leaching breakpoint [19,20]. Any analysis of the leaching breakpoint-AONR relationship, therefore, needs to be conducted across a wide range of environments and management variables as differences in soil characteristics, climate, tillage, and cropping system may be influencing both the AONR and the leaching breakpoint and, thus, the size of the corresponding buffer [6]. Understanding how and why the breakpoint may differ across sites is also important when considering adjustments to N fertilizer recommendations to mitigate NO 3 -N leaching at a regional scale [21,22].
As NO 3 -N leaching originates from residual soil N, mineralized organic matter, and additional applications of N fertilizer, it is likely that the leaching breakpoint in a given season is influenced by the previously planted crop and management, that is the cropping system. The two cropping systems that account for about 72% of cultivated land in the top 12 contiguous maize-producing states in the U.S. (comprising the U.S. Midwest) are continuous maize and a two-year rotation of maize with soybean in which, typically, no N fertilizer is applied to the soybean phase [23,24]. Maize rotated with soybean yields 15% higher than continuous maize despite receiving 30% less N fertilizer input, raising questions about how the fate of N inputs differs between the two systems [25]. Moreover, the organic N inputs, both from previous crop residue and soybean N fixation may have an additive effect on how much N is lost from N fertilizer. Previous studies that have looked at system effects on leaching are often short-term (3-4 years) and have noted that extreme variability in seasonal precipitation appears to influence the fate of NO 3 -N in the soil profile, particularly during a soybean year, potentially masking any differences in the relationship between fertilizer N rate and NO 3 -N leaching load in these two cropping systems [11,13,[26][27][28][29][30]. There is need to expand these findings to understand if and how cropping system selection alters the impact N fertilizer rate on NO 3 -N leaching and if this relationship is consistent across multiple environments (climate and soil).
Cropping system models can successfully simulate NO 3 -N leaching [31]. These process-based models can also investigate the factors driving NO 3 -N leaching independent of each other without losing the complexity of the full cropping system. Specifically, Agricultural Production Systems sIMulator (APSIM) ( [32]; www.apsim.info) has been found to accurately predict both AONR [18] and NO 3 -N leaching [15] in the U.S. Midwest. In a 2014 analysis, Laan et al [31] called upon researchers to use such processbased models in conjunction with data from multiyear field studies to improve our understanding of the processes driving NO 3 -N leaching. Roberts et al [33] found that using a statistical model in combination with a process-based model significantly improved the accuracy and scope of the process-based model's predictions. Therefore, in this study, by fitting a statistical model to the outputs from APSIM, we distilled APSIM's complexity into select parameters that we then used to answer two targeted questions: (a) is the leaching breakpoint related to the AONR? (b) does this relationship differ with site location and/or crop  [32]. The platform is customizable to different environments and management strategies and can be run subsequentially to capture multi-year legacy effects [34]. It comprises of integrated crop and soil modules that accurately simulate water, temperature, carbon (C), and N dynamics [35][36][37]. In this project, we used a modified APSIM version 7.9 with waterlogging capabilities [38,39].

Model setup and calibration
We calibrated the APSIM model using 56 siteyears of NO 3 -N leaching data sourced from eight artificially subsurface-drained field experiments located in the U.S. Midwest (figure 1; study details can be found in supplementary table 1 (available online at stacks.iop.org/ERL/16/064063/mmedia)). These experiments measured leaching as well as yield and/or drainage responses to changes in cropping system, tillage, fertilizer N rate, and subsurface drain spacing. As there is a wide variety of approaches to measuring drainage flow and NO 3 -N leaching, we selected only studies that used artificial subsurface drain outlets to gather their data rather than lysimeters or suction cups in order to maintain consistency in the data across sites. The other criteria in study selection were that they were multi-year experiments and applied inorganic (not manure) N fertilizer with a single late winter/spring application. All experiments were rainfed systems. The research design details and data used to model these experiments were primarily sourced from their corresponding publications [11,[40][41][42][43][44][45][46][47].
Input to the model included field-specific soil and weather data (daily air temperature, precipitation, and solar radiation) extracted from publications or public soil-weather sources ( [48][49][50], https:// power.larc.nasa.gov/). The soil profiles and weather summaries are provided in the supplementary materials (supplementary tables 2-3). The model was calibrated using the management scheme of the original study (i.e. cropping system, N source/timing/rate, tillage, subsurface drain depth/spacing etc.).

Model simulation experiment
Two cropping systems were considered: continuous maize and maize rotated with unfertilized soybean. For the two year rotation, the model was run twice: one set of simulations with maize followed by soybean and the other with soybean followed by maize such that we could compare the outputs from rotated maize (maize in the two-year maize-soybean rotation) with continuous maize each year.
For each site and cropping system combination, the calibrated model was run for five years (1994)(1995)(1996)(1997)(1998)(1999) with standard management per location in the corresponding rotation to stabilize fast decomposing N and C pools and estimate initial conditions [51][52][53]. The model was then run for 20 years (2000-2019) at seven different N fertilizer rates (0, 56, 112, 168, 224, 250, and 300 kg N ha −1 ). Maize and soybean cultivars were sourced from the APSIM database with different maturity groups per location (supplementary tables 4-5 [54]). The USDA-NASS 50% planting dates for each respective state were used for each season and simulation [55]. Tillage, N fertilizer source/timing, subsurface drain depth/spacing and soils parameters were kept consistent with the original calibrated study.

Statistical analysis 2.2.1. Metrics for APSIM model calibration
In the calibration process, to test the model's accuracy in capturing grain yield, drainage, flow-weighted NO 3 -N, and NO 3 -N load (kg ha −1 ), we compared the simulated data from the calibrated models with the corresponding observed data. The statistical analysis of the model's performance was conducted using R version 3.6.2 [56]. Goodness of fit analysis was measured by calculating r 2 , root mean square error (RMSE) and modelling efficiency (ME) [57].

Analysis of APSIM model outputs
All processing of the APSIM simulation outputs and statistical analyses were done using the tidyverse meta package [58] in R version 3.6.2 [56]. All other packages are cited below.
APSIM runs on a calendar year, but to facilitate comparisons between the rotations we wanted drainage and NO 3 -N leaching to be summed from crop sowing to crop sowing. In order to accommodate this timeline, drainage and NO 3 -N leaching were categorized into two groups: season/postharvest (sowing to the end of year, 31 December) and pre-season (1 January to the day before sowing). The season/post-harvest data of one calendar year were then combined with the pre-season data of the next year to capture the full effect of the sown crop on drainage and leaching both during the season and following harvest.
For all systems, the leaching from maize sowing to maize sowing was compared. For the cropping system consisting of maize rotated with soybean, the data from the following soybean year were added onto that of the maize year to include in the analysis any residual N from the fertilizer applied in the maize year that leached out during the soybean year. A diagrammatic representation of this method is presented in supplementary figure 1.

Leaching
Data were grouped into experimental units, which were defined as a site-year-rotation combination, resulting in a total of 294 units. We fit three candidate non-linear models to N leaching as a function of N fertilizer rate using the nlraa package [59]. The three models investigated were (a) bi-linear, (b) exponential, and (c) exponential-linear [60]. We chose these three models because they offer meaningful parameters and/or are commonly found in literature. Final model selection was based on Akaike's Information Criteria [61], examination of residual plots, and knowledge of biophysical constraints to these systems.
To assess the impact of rotation on the statistical model parameters, we fit a non-linear mixed effect model using the nlme package [62]. We used a fixed effect for rotation (continuous, rotated) and random effects for site and experimental unit nested within site. Estimates and contrasts for the effect of rotation on the parameters were assessed using the emmeans package [63]. Contributions of random effects were assessed using interclass correlation.

Yields
Several models have been proposed to determine the AONR for maize. We chose to use a bi-linear model for its simplicity and applicability to a range of scenarios [64]. While this model may not capture the nuances of real-world relationships, in our simulations it provided a more robust estimation of the yield breakpoint, which is representative of the AONR. We fit a non-linear mixed effect model as described above.
A bi-linear model outputs lower AONR values than the more conventionally used quadratic plateau model [64]. As such, we have included quadratic plateau model-derived AONR values for all site-years in the supplementary materials for comparative purposes (supplementary table 6).

Relationship between leaching breakpoint and AONR
The difference between the leaching and yield breakpoints was defined as the 'buffer,' with positive values indicating the leaching breakpoint occurred at a higher N fertilizer rate than the yield breakpoint. The conditional breakpoint parameter estimates for the leaching and yield experimental units were calculated on a per-experimental-unit basis. The buffer for each rotation was compared using a linear model with rotation as a fixed effect.

Model calibration and simulation
The APSIM model simulated yield, drainage, flowweighted NO 3 -N, and NO 3 -N leaching load well with ME values falling primarily between 0.7 and 0.95. In two studies, however, model accuracy evaluation metrics reflected modelling limitations caused by a lack of reported data: leaching and flowweighted NO 3 -N data from MI [43] were reported as three year averages, resulting in an ME of 0.46 and yield data from IN2 [47] were reported as six year averages, resulting in a high RMSE Yield (11.22 kg ha −1 ) (supplementary figures 2-5; supplementary table 7). Most importantly, the model captured treatment differences such as till vs no-till, narrow vs wide subsurface drain spacing, and cropping systems (supplementary figures 6-8). There was a wide range of AONR values found for both continuous (66-170 kg N ha −1 ) and rotated maize (24-119 kg N ha −1 ); the average AONR, however, was 111 kg N ha −1 and 70 kg N ha −1 for continuous and rotated maize, respectively (figure 2(b); supplementary figure 9; supplementary table 6).

Leaching model
Of the 294 simulated site-years, 277 had all three tested non-linear models (bi-linear, linear-exponential, and exponential) converge for the NO 3 -N leaching load's response to fertilizer N rate. For 92% of those site-years, the bi-linear model fit was best and for 8%, the linear-exponential fit was best. The exponential fit was not found to be the best fit for any of the site-years. The rest of our analysis, therefore, will be based on the bi-linear model fit (figure 2(c)).
When no fertilizer was applied, significantly more NO 3 -N was leached from the maize-soybean rotation than from continuous maize (13 kg N ha −1 and 6 kg N ha −1 , respectively). Both systems' leaching loads varied significantly by site (28% of variation) and, to a lesser extent, by year (16%) ( figure 2(c)).
The parameters defining the slope below and above the breakpoint as well as the breakpoint itself differed significantly with cropping system. Location and year, however, did not have a significant influence on these parameters. Below the breakpoint, continuous maize lost 0.08 kg NO 3 -N kg −1 N applied while the maize-soybean rotation lost 0.1 kg NO 3 -N. There was then a breakpoint at the fertilizer N rate 129 kg N ha −1 (SE: 0.6) for continuous maize and at 116 kg N ha −1 (SE: 1.9) for the maize-soybean rotation at which the rate of leaching changed. Above the breakpoint, continuous maize lost an average of 0.6 kg NO 3 -N kg −1 N applied (95% CI: 0.54-0.63), but the maize-soybean rotation lost only 0.4 kg NO 3 -N kg −1 N applied (95% CI: 0.37-0.43).

Relationship between AONR and leaching breakpoint
There was a greater margin for error in rotated maize than in continuous maize around overestimating a given field's AONR without drastically increasing the rate of NO 3 -N leaching. The leaching breakpoint in the model occurred at an N rate 46 kg N ha −1 (95% CI: 44-49 kg N ha −1 ) above the AONR in the maizesoybean rotation, but only 17 kg N ha −1 (95% CI: 14-20 kg N ha −1 ) above the AONR in continuous maize ( figure 2(a)). The N rate that minimized the yieldscaled leaching load (i.e. the NO 3 -N leaching load per unit maize grain yield) was lower than the leaching breakpoint in both cropping systems (supplementary figure 10).
Within the maize-soybean rotation, at N rates above the breakpoint, leaching loads during soybean following maize were higher than those in maize following soybean, pointing to the residual effects of applying excessive amounts of N during maize ( figure 3). However, the leaching per year from the maize-soybean rotation system (be it the maize or soybean year) was always lower than that of continuous maize.

N fertilizer-leaching relationship
Research suggests the N leaching-fertilizer rate breakpoint occurs around the AONR, however, we found that the leaching breakpoint occurs at N rates 16% above the AONR in maize following maize and 66% above the AONR in the maize-soybean rotation. There is, therefore, less risk in a maize-soybean rotation that farmers will have to under-fertilize while trying to minimize environmental N losses than in a continuous maize system. A meta-analysis also found that the leaching breakpoint occurs at least 15% above the AONR [1]. Our results expand upon this previous finding as it made no distinction between continuous and rotated maize: we found that system selection impacts where that breakpoint occurs. Moreover, in our analysis, we found that the breakpoint was not a function of site or year, only cropping system, whereas AONR was strongly influenced by both site and year.
Our results suggest that reductions in N fertilizer inputs have little risk of reducing productivity per unit N loss. We found that the N rate that optimized the leaching load per unit grain produced was lower than the leaching breakpoint in both systems. Because the N leaching breakpoint occurred well above the AONR and yield-scale N leaching was minimized at a N rate below the N leaching breakpoint, there should be opportunities to reduce N leaching without losing productivity. Basing our analysis on long-term simulated experiments rather than on data from shorterterm experiments strengthens our confidence in these conclusions.

Rotation effect
Above the breakpoint, leaching per unit N applied in a single year of continuous maize increased at a rate 1.5 times that of a full two-year cycle of a maize-soybean rotation. There was also a larger buffer between the AONR and breakpoint in the maize-soybean rotation than in continuous maize.
In APSIM, the size of the leachable N pool in each soil layer is calculated as the balance of mineralized soil N (either sourcing from soil organic matter or N fertilizer) that is not immobilized, denitrified, or taken up by the plant. The actual fate of this N is driven by the drainage of water through the soil profile. The model showed that at low N rates, the mineralized soil N pool was smaller under continuous maize than under the rotated system. When the fertilizer input exceeded the breakpoint, the mineralized N pool under continuous maize increased to the same size as that under rotation. This increase has also been found in field experiments where increasing the applied N rate above a certain threshold N rate (between 100 and 200 kg N ha −1 ) increased the residual soil N pool under continuous maize systems more than under rotated systems [65,66]. This systemdependent effect of fertilizer on the mineralized N pool size is further evidence of a greater leaching risk under continuous maize in response to increases in N fertilizer rates.
In the model, rotated maize took up more N than the continuous maize. This gap in uptake efficiency increased as the N rate increased. Varvel and Peterson [65,67] found that in a maize-soybean rotation, maize took up around 50% of the N applied to the crops regardless of the N rate. Meanwhile, when grown continuously, maize took up 50% of N applied at low N rates but only 30% at higher N rates [65]. The smaller buffer we found between the AONR and the leaching breakpoint under continuous maize, therefore, stems from both N saturated soil and low N uptake efficiency.
When evaluating the fate of N in a two-crop rotation, it is important to consider how much leaching occurs throughout the two-year cycle, as much of the NO 3 -N that was not leached during the maize year may be lost during the soybean year [30]. Our simulations found that within the two-crop rotation, soybeans experienced more drainage and leaching than rotated maize (figure 3). Higher drainage levels correspond to more leaching, but that is likely not the only reason for the greater leaching loads during the soybean phase. While soybean took up the same amount of N as rotated maize, much of that N sources from N fixation, leaving the residual N fertilizer from the maize phase of the rotation to leach from the soil profile [68]. Carryover N from the previous maize year and high drainage levels during soybean year can increase the leaching load in soybean years almost to the level of continuous maize [4,7].
At high N rates, like those above the breakpoint, continuous maize accumulated a greater pool of mineralized soil N (leachable N) relative to the rotated cropping system. High drainage loads in the soybean years of the maize-soybean rotation resulted in the loss of more N than in rotated maize years, but still less than in the continuous maize system.

Model robustness
We found that the cropping system dictated the shape of the leaching response and breakpoint while the site and year (which encompasses annual differences in plant growth, temperature, precipitation, and soil N pool size/composition) explained the variability in the magnitude of the leaching. Therefore, while our leaching model can be applied across multiple sites and years to determine the degree to which increasing the N fertilizer rate impacts the leaching load, use of an average site-year prediction should be avoided in calculations of total leaching load from multiple sites/years at a given N rate.
The large range of baseline yields predicted may be indicative of minor limitations in the accuracy of APSIM model. Puntel et al [18] found that while, in general, APSIM simulates yields well across the range of the N rates used in this analysis (average Relative RMSE Yield = 14.9% and 11.6% in continuous and rotated maize, respectively), its accuracy is lower when 0 kg N ha −1 is applied to continuous maize (Relative RMSE Yield = 30.8%) than at non-zero N rates (Relative RMSE Yield <16% at N rates ranging from 67 to 268 kg N ha −1 ).
The significant effect of site and year on the baseline leaching load is further evidence of the important role soil and weather have in determining the fate of residual N over multiple seasons. Residual N can be leached from the soil profile or taken up by the plant, resulting in high leaching loads even in seasons when no N fertilizer is applied. These legacy effects may explain why previous studies have found inconsistent and/or inconclusive differences in leaching response to N fertilizer between cropping systems depending on the site, system, and year [13,[26][27][28][29][30].
Establishing and maintaining long-term field drainage experiments with as many different N-rates and cropping systems as we simulated (7 N rates × 3 phases × 4 replication = 84 plots × 20 years) is practically impossible as it is expensive and laborious. Moreover, artificially drained experiments tend to be constrained by field size more than other field experiments. For instance, the drainage experimental site with the most plots in the US Midwest has 72 plots, but the plots are small (∼0.05 ha) [11]. Using cropping system and statistical models to expand upon field experiments is the only way to delineate the complex relationships between N rate, AONR, and leaching.

Implications of findings
The smaller buffer between the AONR and the breakpoint in continuous maize suggests that the risk of negatively impacting water quality via overfertilization in continuous maize is much greater than in a maize-soybean rotation.
A larger buffer between the AONR and breakpoint can reduce the risk that fertilizer overapplication leads to groundwater contamination. After conducting a state-wide survey, Sellars et al [69] reported that 67% of fields in Illinois received N fertilizer rates that exceeded the AONR. Combining their findings with our own demonstrates the potential impact replacing continuous maize with a maize-soybean rotation can have on NO 3 -N leaching in the U.S. Midwest. If the 67% statistic observed in Illinois is assumed consistent across continuous maize cropping systems in the U.S. Midwest, changing from continuous maize to rotated maize could greatly reduce the leaching from around 5 million hectares of cropland across the region.
This cropping system approach to reducing N leaching is an inexpensive but effective N leachingmitigating strategy [70]. Nevertheless, there is room for improvement as there has been very little research on what might be a more desirable rotation crop in simple U.S. Midwestern two-crop rotations from a leaching perspective.

Conclusion
The cropping system selection plays a significant role in defining the degree to which fertilizer N rate impacts NO 3 -N leaching as well as the rate above which leaching increases significantly. Our findings were robust across both environment and climate variability and can provide sound guidance for how farmers can improve downstream water quality.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: https:// github.com/vanichols/ERL2021_data-analysis.