Simulated nitrate leaching in annually cover cropped and perennial living mulch corn production systems

Materials and Methods

Site Description. The study was performed from October of 2014 through February of 2017 at the J. Phil Campbell Research and Education Center in Watkinsville, Georgia (33°52′09.5″ N 83°26′59.8″ W). The soil is classified as a Cecil sandy loam (fine, kaolinitic, thermic typic Kanhapludults), which was confirmed by pedological analysis of a soil profile adjacent to the research plots. A soil moisture release curve based on the van Genuchten (1980) equation was created using the evaporation method (Arya 2002) with a Decagon HYPROP device (Decagon Devices, Pullman, Washington) from several soil cores collected from the plot area. The soil moisture release curve indicated that field capacity (−0.03 MPa) was 0.24 cm³ water cm⁻³ soil and wilting point (−1.5 MPa) was 0.11 cm³ water cm⁻³ soil. Weather data were collected by The Georgia Environmental Monitoring Network station adjacent to the plot area (http://weather.uga.edu/).

Experimental Design. The cereal rye (CR), crimson clover (CC), and LM treatments were assigned to a randomized complete block design with six replications. Soils in the experimental area were sampled and pH and mineral deficiencies amended according to soil test. The CR treatment received a total of 280 kg N ha⁻¹ supplemental fertilization that was split into two applications based on a recommendation from Lee et al. (2015). The CC and LM received lower supplemental N rates based upon N credits from the cover crops established in a previous study to provide a similar amount of available N as the CR treatment (Andrews et al. 2018; Sanders et al. 2017). CR plots were fertilized with 56 kg N ha⁻¹ at planting using granular urea. At the V6 stage of corn development, the CR, CC, and LM treatments received 224, 168, and 56 kg N ha⁻¹, respectively, also as granular urea.

Agronomic Management. Prior to cover crop planting, soils were tested for pH and extractable phosphorus (P) and potassium (K) by the University of Georgia Soil Testing Laboratory using the Mehlich 1 method (Beck et al. 2004). Soils were limed and fertilized so that pH was above 6.2 and extractable P and K were above 45 and 140 ppm, respectively. On October 17, 2014, land was disked and firmed with a cultipacker before CR, CC, and LM cover crop plots were seeded by hand at rates of 100, 28, and 13 kg ha⁻¹, respectively. The CC and CR plots were killed with a broadcast application of dicamba (3,6-Dichloro-2-methoxybenzoic acid) and glyphosate (N-(phosphonomethyl) glycine) at respective rates of 1.20 and 1.12 kg a.i. ha⁻¹ on March 16 and April 8, 2015. All herbicides were applied by a state licensed pesticide applicator, and all applications were made during environmental conditions when risk for off-target movement of pesticides was low. Additionally, the pesticide formulations, application equipment, and application methods that were used ensured that the risk of off-target movement of pesticides was further reduced (Carlsen et al. 2006). A 20 cm banded application of the herbicides at the same rates was applied to the LM plots on April 8, 2015, using a hooded sprayer. The banded application was centered on 90 cm rows. A John Deere 7300 MaxEmerge no-till planter (John Deere, Moline, Illinois) was used to plant corn (DeKalb DKC64-69, GENVT3P) at a population density of 90,000 plants ha⁻¹ on April 21, 2015. Plot size allowed for eight rows of corn per plot. Corn was harvested on August 16, 2015.

Annual cover crops were reestablished in the same plot areas in October of 2015 using the same seeding rates as 2014. On March 23 and April 4, 2016, herbicide applications were made to the CR and CC plots, while herbicide applications were made to the LM plots on April 14. All herbicides were applied as previously described. Corn was planted at the same seeding rate and methods as the previous year on April 28, 2016. Nitrogen fertilization also occurred in the same rate and method as 2015. Corn from all plots was harvested on August 8, 2016, and cover crops were reseeded on the same plots in November of 2016.

Field Instrumentation, Irrigation Scheduling, and Soil Sampling. Campbell Scientific Model CS625 reflectometers (Logan, Utah) were installed under the corn row to monitor soil water content in the center of each plot. The rods were 30 cm in length and installed at 30 degree angles to measure water content at the 0 to 15 and 15 to 30 cm depths. Water content data were measured and recorded on 10 minute intervals and stored on data loggers within three replications of the experiment. Irrigation was provided on an as-needed basis during the cropping season, which was determined by the soil moisture sensors. Overhead irrigation was used to maintain the water content of the surface 0 to 15 cm between 40% to 90% available water content in 2015, but demand for water from other researchers exceeded the ability of the water resources during 2016, a drought year. Thus, there were occasions when the soil moisture content fell below the desired range.

Two porous cup suction samplers (lysimeters) were placed 1 m apart in the center of each plot alongside the corn row in the other three replications. A Giddings soil core sampler was used to remove soil to a restrictive Bt2 layer, approximately 75 cm deep. A 1 m long lysimeter was installed with a slurry of Bt horizon soil to cover the sampler's ceramic cup, and a slurry of A horizon soil and kaolinite filled the remainder of the augered hole to prevent preferential flow down the side of the lysimeters. Polyvinyl chloride (PVC) pipes with caps were placed over the tops of the lysimeters to prevent damage to the lysimeters and minimize the opportunity for preferential flow.

Lysimeters were purged weekly and a 7 kPa vacuum was applied using a handheld vacuum pump and kept under vacuum for 18 to 36 hours. Water was aspirated from the lysimeters using a handheld vacuum pump, and a 25 mL aliquot of water from each lysimeter was combined for NO₃-N analysis. Occasionally only one lysimeter per plot contained water from which a 50 mL aliquot was sampled for analysis. Eight soil cores were randomly sampled to a depth of 15 cm weekly from each plot using a 1.5 cm diameter handheld soil probe. Soil cores were taken from the center two rows of the plot, combined, air dried, and stored at 4°C.

Laboratory Analysis. Five grams of soil were extracted with 40 mL of 1 M potassium chloride (KCl) for NO₃-N concentration. Water samples and soil extracts were analyzed for NO₃-N concentration using a Timberline TL-2800 Ammonia Analyzer (Boulder, Colorado).

Weekly changes in clover mass in the LM plots were determined using a rising plate meter (RPM) as described by Sanders et al. (2017). Ten to twelve RPM height measurements were made within the center two rows in each LM plot weekly to determine clover mass. Hand harvested clover samples were subjected to Kjeldahl digestion and N content determined using a Timberline TL-2800 Ammonia Analyzer. Total N in the white clover was calculated weekly as the product of white clover mass and the total N content of the clover. Decline in clover mass from one week to the next was used to estimate N release from the clover.

Model Selection and Description. Porouscup suction lysimeters are a cost-efficient means of obtaining leachate from a large number of treatments (Andraski et al. 2000). However, they are not capable of quantifying water drainage. Therefore, a modeling approach was used with solute concentration data to evaluate NO₃-N movement. HYDRUS-1D was selected for simulating water and NO₃-N loss in the different cover crop corn production systems in our study due to the relative ease of use and use in similar research (Tafteh and Sepaskhah 2012; Wang et al. 2010). The model was used to simulate water and NO₃-N movement from planting in April of 2015 and 2016 through February of the following year based on mean data from each treatment. One model was used to simulate each corn-cover crop system per period. The one-dimensional model was developed to simulate the vertical movement of soil water, heat, and solutes in variably saturated-unsaturated media (Šimůnek et al. 2016) using a numerical solution to the Richards (1931) equation 1: 1 where h (cm) is the pressure head, z (cm) is the gravitational potential head, t (d) is time, and S (cm d⁻¹) is root water uptake rate.

The van Genuchten (1980) soil water retention equation 2 was used to describe water retention in the experimental soil: 2 where α (cm⁻¹), m (dimensionless), and n (dimensionless) are fitted parameters of the van Genuchten equation; θ(h) is the volumetric water content (cm³ cm⁻³); θ_s is the saturated volumetric water content (cm³ cm⁻³); and θ_r is the residual volumetric water content (cm³ cm⁻³). The van Genuchten (1980) unsaturated hydraulic conductivity equation 3 was used to describe unsaturated flow in the experimental soil: 3 where K_s is the saturated hydraulic conductivity (cm d⁻¹), m is the fitted parameter from equation 2, and it is assumed that m = 1 − 1/n.

Initial estimates of the water retention parameters for equations 2 and 3 were determined by the evaporation method (Arya 2002) using a Decagon HYPROP device (Pullman, Washington) from soil samples taken from each soil layer in a pit adjacent to the experiment. Soil water retention parameters were optimized through auto-calibration in the top two soil horizons across all models using collected soil moisture data and the inverse solution function in HYDRUS-1D. Since there were no soil moisture data for the lower two soil layers and the inverse solution function in HYDRUS-1D required data from those depths, the soil water retention parameters in the two remaining soil horizons were not changed from the initial parameters. Soil physical properties and final soil water retention parameters are in table 1. Soil samples for bulk density (ρ_b) (g cm⁻³) were also collected and measured following standard procedures (USDA 1996).

Model Domain and Boundary Conditions. The HYDRUS-1D model used a 100 cm profile with four materials to simulate four soil horizons and their water retention parameters (table 1). A total of 101 nodes were used in the model space with an equal density of one node per centimeter of model space. The upper boundary condition was an atmospheric boundary with surface runoff and the lower boundary condition was free drainage. Observation nodes were placed at 8 and 23 cm below the soil surface to represent the center of time domain reflectometers installed within the corn rows. Observed data were input at each observation node as the average daily volumetric water content (cm³ cm⁻³) of the three replications.

Solute transport had a concentration flux boundary condition at the top of the model space with no concentration gradient as the lower boundary condition. Observation nodes were placed at 8 and 75 cm below the soil surface to represent points where soil water NO₃-N concentration was sampled either by soil sampling or suction samplers. Observed data were input in each observation node as the daily average NO₃-N concentration of the three replications. Granular urea fertilizer N was used as an input time variable for the upper boundary condition and the soil water concentration calculated based upon the amount of irrigation water applied. Nitrogen release from the clover in LM plots was added when clover mass declined from one week to the next.

Water Uptake and Evapotranspiration. The Feddes et al. (1978) model (equation 4) was used to describe root water uptake in HYDRUS-1D: 4 where α(h) is the coefficient of root water uptake and S_max is the potential water uptake (d⁻¹). Root water uptake is zero at saturation (h₀) due to lack of oxygen (O) and increases linearly as pressure heads decrease to h₁. Root uptake was optimum (no water stress) for pressure heads within the range of h₁ to h₂. From h₂ to h₃, root uptake decreased linearly and stopped at pressure heads below h₃ due to insufficient soil water. Based on R² provided in model output, the root water uptake parameters that provided the best fit to observed data and subsequently used in this model for all crops were h₀ = −15 cm, h₁ = −30 cm, h₂ = −400 cm, h₃ = −15,330 cm, and S_max = 0.5 cm d⁻¹. For reference, the parameters used for corn in Wesseling et al. (1991) are h₀ = −15 cm, h₁ = −30 cm, h₂ = −600 cm, h₃ = −8,000 cm, and S_max = 0.5 cm d⁻¹.

Input rooting depth for corn increased linearly from 0 cm at planting to a final depth of 72 cm after 82 days of growth and remained at that depth until harvest (Mengel and Barber 1974). An additional 10 cm of rooting depth was added to the corn in the LM treatment to simulate competitive water uptake by the intercropped clover (Kurtz et al. 1952). Rooting depth increased linearly from 0 to 50 cm for the CC and CR cover crops during both model periods, and from 0 to 40 cm in the LM cover crop during the 2015/2016 model period (Caradus 1990; Sainju et al. 1998). Initial to final rooting depth for all cover crops occurred from planting to 60 days after planting (DAP). In the 2016/2017 model period, rooting depth of the LM cover crop was input as 10 cm at corn harvest and then increased linearly to 40 cm after 45 days of growth. As HYDRUS-1D only allows for one root distribution and two crops are simulated in each model, a root distribution to represent corn and all three cover crops was used. Seventy percent of roots were placed in the top 20 cm of soil, which generally occurs in all four crops (Caradus 1990; Mengel and Barber 1974; Sainju et al. 1998). Rooting distribution decreased linearly from 20 to 72 cm (the top of a restrictive clay layer) for the remaining 30% of roots. While the same rooting distribution was used for multiple cover crop systems, it should be noted that each of the three corn-cover crop simulations was performed independently and not incorporated into one model.

Potential evapotranspiration (PET) was separated into evaporation and transpiration in the HYDRUS-1D model by the input of leaf area index (LAI). In the CC and CR models, LAI of corn increased linearly from 0 to 3.0 from 7 DAP to 67 DAP, and then from 3.0 to 3.5 linearly over the next 30 days and remained until harvest (Timlin et al. 2014). For the LM models, an additional LAI of 1.5 was added to account for competitive water uptake by the white clover (Kurtz et al. 1952). Leaf area index of the CC, CR, and LM plots increased linearly from 0 to 3.0 from 7 days after cover crop seeding (DAS) to 69, 37, and 80 DAS, respectively. Leaf area index increased linearly from 3.0 to 5.1 and 3.0 to 4.6 until 82 and 89 DAS for the CC and LM treatments, respectively (den Hollander et al. 2007; Ellen 1993). The LAI remained constant thereafter until the end of February in the following year (den Hollander et al. 2007). In the next model period, white clover in the LM treatment was not reseeded, so LAI followed the same development as the previous LM model, but began with an LAI of 1.0. Therefore, LAI of the white clover increased linearly from 1.0 to 3.0 after 52 DAS (though no white clover was planted), and then from 3.0 to 4.6 during the 56 to 65 DAS time period.

Soil Nitrogen Transport and Transformation. Linear adsorption and chemical equilibrium in the advection–dispersion equation (ADE) was used to estimate solute transport in HYDRUS (5): 5 where c is the dissolved concentration of the solute (g cm⁻³), t is time (d), K_d is an adsorption coefficient (cm³ mg⁻¹), D_e is the effective dispersion coefficient (cm² d⁻¹), z is the vertical dimension (cm), J_w is the vertical Darcy water flux (cm d⁻¹), and μ is the first-order rate constant for solute transformation processes (d⁻¹).

Mineralization and nitrification can occur rapidly, and concentrations of potentially mineralizable N (PMN) and ammonium (NH₄) are low in Cecil soils (Dubey 1968; Anderson 1960; Cabrera 1993). Therefore, PMN and NH₄-N in the soil profile were considered together as a bulk solute. Other relevant HYDRUS-1D NO₃-N leaching research used a similar strategy involving transformations of PMN and NH₄-N into NO₃-N (Tafteh and Sepaskhah 2012; Wang et al. 2010).

An N chain model was used to simulate N transformations in the CC, LM, and CR models. The model included the nitrification of NH₄-N and PMN to NO₃-N as well as denitrification from NO₃-N to dinitrogen/nitrous oxide (N₂/N₂O). The intermediate product of nitrification, NO₂-N, was ignored to simplify the chain model (Bradshaw et al. 2013). The first-order reaction rates for the change in NH₄-N + PMN (equation 6) and NO₃-N (equation 7) concentrations that result from the N chain models were 6 and 7 where k and μ are the nitrification and denitrification rate coefficients (d⁻¹), respectively, and t is time (d). It was assumed that these coefficients only applied to N in solution. Based on percentage bias (PBIAS), R², and d, the nitrification rate constant of 0.2 d⁻¹ was found to provide an acceptable fit to observed data and was within ranges reported in the literature (table 2) (Hanson et al. 2006; Iskandar and Selim 1981; Ling and El-Kadi 1998; Lotse et al. 1992). Complete denitrification of NO₃-N was assumed to occur under saturated conditions, and rate constants were set the same as for nitrification (table 2). Adsorption coefficients for NH₄-N/PMN and NO₃-N were 10 and 0.78 cm³ g⁻¹ for all soil horizons (table 2), and were based on values from Bradshaw et al. (2013). Based on R² in model output, the adsorption coefficient was optimized across models using manual calibration to provide a more acceptable fit to observed NO₃-N concentrations at the 8 and 75 cm depths.

Water content dependence of nitrification and denitrification were incorporated into the solute models. Water content dependence of reaction rates uses a modified version of the Walker (1974) equation 8: 8 where ω_r is the reaction rate constant (d⁻¹) at the reference water content (θ_ref), ω(θ) is the rate constant (d⁻¹) at the actual water content (θ), and B is a dimensionless solute-dependent parameter. Nitrification and denitrification rates were considered to be optimum at pressure heads of −500 cm and −20 cm, respectively.

Model Initial Conditions and Inverse Data. The initial condition for soil pressure head at the uppermost observation node in each model was set to the corresponding water content that was measured at the beginning of the model period. The rest of the model space was set so that soil water movement was only due to gravitational flow (constant pressure head with depth). Inverse data for water content were input into HYDRUS-1D as the average daily measured volumetric water content of the three replications. Initial NO₃-N concentrations were input into the model based on both soil samples and lysimeter samples taken immediately before corn planting. Average soil water NO₃-N concentrations were determined by taking the average soil NO₃-N level from the three replications and converting it to a concentration using the water content and bulk density of the soil.

Model Performance Criteria. Observed data in the model were average soil water content at the 8 and 23 cm depths and average soil water NO₃-N concentration at the 8 and 75 cm depths. As a plot-scale model, no performance evaluation criteria are currently established in the literature for HYDRUS-1D that evaluate the agreement between predicted and observed data, and therefore determined by the objectives and statistical performance measures selected by the user (Moriasi et al. 2015). Therefore coefficient of determination (R²), PBIAS, and index of agreement (d) were used to assess calibration and validation model performance along with visual inspection of graphical data.

R² describes the proportion of the variance in measured data explained by the model (Moriasi et al. 2007). R² ranges from 0 to 1, and the degree of collinearity between predicted and observed data increases as values approach 1 (Moriasi et al. 2015). PBIAS measures whether the predicted data tend to be smaller or larger than their observed counterparts (Gupta et al. 1999). Positive values indicate an underestimation of water flow or NO₃-N transport while negative values indicate overestimations of the parameters. Index of agreement (d) is the ratio between the mean square error and the potential error with values ranging from 0 to 1. As d values approach 1, the agreement between predicted and observed values increases (Willmott 1984).