Nitrogen Fertilizer Management Based on Site-Specific Maximized Profit and Minimized Pollution

Study Site

Seven fields in Montana representing different climatic and soil regions were selected from a larger population of fields for which data is being collected for the OFPE project led by Dr. Bruce Maxwell (Fig. 1). The benefit of using producers who have prior experience and participation in site-specific management with OFPE is their experience in research management on their fields as well as our access to their historical data. Fields were selected based on the crop grown in 2020 and 2021 (winter wheat) and on the quality and promptness of data collection by the producers. Additionally, producers covered a three different cooperator farms, generating a dataset with measurements in different geographic locations.

The response variables of interest were crop productivity (yield in bushels per acre) and quality (grain protein percent), both of which were gathered from monitors mounted in the cab of the farmer’s combine harvesters. Yield was collected about every three seconds while protein was measured about every ten seconds. Yield is related to price received for wheat based on volume (bu), however protein in winter wheat is judged based on a nonlinear scale. On typical years, protein content above 11.5% is rewarded with a 2 cent per half point protein increase on the base price received, however if protein is below 11.5%, the price received is decreased by 8 cents per half point protein down to 9.5%. Thus, the combination of quantity and quality of the winter wheat crop influences farmer net returns. Farm equipment used in this study varied between farmers, with one using John Deere and the other using Case. All yield monitors are calibrated every spring by the farmers, according to their respective manufacturing instructions. Grain protein content was measured with Next Instrument’s CropScan 3300H near infrared monitor (product citation) which was installed, maintained, and calibrated by Triangle Ag Crop Consultants (website citation). Both yield and protein data are subjected to proprietary cleaning practices according to the software of the machine used to make the measurement. In addition to the yield, protein, and as-applied N data that were collected from the machines on the field, remotely sensed covariate data from open sources was gathered (Table 2). These data were gathered or derived from Google Earth Engine (Gorelick et al. 2017) and aggregated together at the locations of the observations in the yield and protein datasets.

Table 2. Table of covariate data types gathered from Google Earth Engine to enrich the crop yield and protein datasets gathered from on-farms.

Crop Response Modeling

The first objective was achieved by characterizing the response and variation in yield and protein data to N fertilizer application using multiple statistical methods such as regression, Bayesian, and machine learning approaches. These models were fit using the large amounts of data collected from yield harvesters and protein analyzers mounted on the combine, and from aggregated 10m scale covariate data aggregated from the earth observation satellites.

1) Modified Non-Linear Sigmoid Model (Beta Function)

The first model type was parametric and modified from a non-linear sigmoidal model called the beta function (Yin, Goudriaan, Lantinga, Vos, & Spiertz, 2003). Following the literature, a model that assumes an asymptotic nature of crop response curves to additional N fertilizer was tested here (Anselin, Bongiovanni, & Lowenberg-DeBoer, 2004; Reynolds, Drury, Phillips, Yang, & Agomoh, 2021; Yin et al., 2003). Initial tests compared the beta function to multiple forms of non-linear models including: hyperbolic, Richards, Gompertz, Weibull, and  logistic forms. Selection of the beta function (Yin et al., 2003) over these other non-linear models was based on three factors; 1) the beta function yielded the lowest mean RMSE across fields for both yield and protein (Table 4), 2) the beta function required the least locking of parameters to facilitate convergence, and 3) the beta function captured the asymptotic response and a downturn of responses at high experimental N fertilizer rates.

Table 4. Results of the mean RMSE across the seven fields of yield and protein from each of the non-linear models tested with a randomly split training dataset with 70% of observations across three years for model fitting and a hold-out dataset with the other 30% of observations for predicting responses with new data and calculating RMSE.

The beta model was modified by incorporating an intercept term to represent the response of the crop at zero N-rate. Additionally, a Gaussian spatial correlation structure on the UTM coordinates of observations to account for the spatial dependence of responses was included (De Bastiani, Mariz de Aquino Cysneiros, Uribe-Opazo, & Galea, 2015; Diggle & Tawn, 1998; Mardia & Marshall, 1984; Xia, Ding, & Wang, 2008). The beta function was fit using generalized non-linear least squares regression through the nlme package in R (Pinheiro, Bates, DebRoy, Sarkar, & R Core Team, 2021). Bottom-up feature selection was performed where predictors that did not contribute to a decrease in AIC by two units were omitted from the model. As crop responses differed across fields, not all beta functions fit for yield or protein for a given field and set of training years have the same predictors. All final models for each field were selected from the following beta function used for yield (kg ha^-1) or grain protein concentration (%) forecasting, with a Gaussian process spatial correlation structure, and the predictors that influence the intercept term (α);

(1)

Where R is yield (kg ha^-1) or grain protein concentration (%), α is the fit minimum yield at a zero N-rate, β is the fit maximum yield (asymptote) at high N-rates, N is the as-applied N-rate in kg ha^-1, δ₂ is the fit N-rate where β is approached, δ₁ is the fit N-rate at the inflection point (center of the upturn of the curve), Gaussian(x, y) is the Gaussian process spatial correlation structure, and ε is random error assumed to be normal and independent N(0, σ²). The intercept term, α, was assumed to be influenced by all topographic variables and information from prior years;

(2)

Where α₀ - α₁₅ were coefficients for the intercept and covariate data listed in Table 3.

2) Generalized Additive Model (GAM)

Staying in the non-linear realm but relaxing assumptions on the shape of the model for crop responses, a generalized additive models (GAMs) was tested. Generalized additive models are a blend of semi-parametric generalized linear models with additive smoothing terms for modeling non-linear relationships between covariates and the response (Hastie & Tibshirani, 1986; Wood, 2017; Pinilla & Negrín, 2021). Despite their flexibility, the use of GAMs for modeling crop responses has been limited (Chen, O’Leary, & Evans, 2019; Estes et al., 2013; Joshi, Kazula, Coulter, Naeve, & Garcia y Garcia, 2021; Liew & Forkman, 2015). The GAMs for grain yield and grain protein concentration used gamma family distributions and thin plate splines with shrinkage for the basis functions on most of the covariates in Table 3. A gamma family distribution was selected to ensure that estimates and predictions were above zero in accordance with the reality of grain yield and grain protein concentration measurements. Shrinkage allowed estimated degrees of freedom to fall to zero, combining the model fitting and feature selection process. The only covariates that did not use a thin plate basis function were the spatial UTM coordinates, which were included as smoothing terms with a Gaussian process basis function to account for spatial autocorrelation (Gotway & Stroup, 1997; Guisan, Edwards, & Hastie, 2002; Holland, De Oliveira, Cox, & Smith, 2000; Zuur & Camphuysen, Kees, 2012). To account for non-constant variance found in initial model fits, a log link function was used. Models were fit using the mgcv package in R (Wood, 2003; Wood, Pya, & Säfken, 2016). As crop responses differed across fields, not all GAMs for yield or protein for a given field and set of training years have the same predictors but selected from the following initial GAM with a gamma family distribution and log link function;

(3)

Where R is the expected yield (kg ha^-1) or grain protein concentration (%) with a gamma log link function, I is the intercept,  f₁ is the tensor product of the x and y coordinate (x, y) with a Gaussian process basis function,  f₂ is the tensor product of centered aspect in the E/W direction (aspect_cos) and aspect in the N/S direction (aspect_sin) with a thin plate shrinkage spline basis function,  f₃ – f₁₇ are smoothing functions fit with maximum likelihood and thin plate shrinkage splines as the basis function for covariates listed in Table 3, and ε is the random error, assumed to be normal and independent N(0, σ²).

3) Non-Linear Bayesian Regression Model (Bayes NLR)

Two Bayesian models were also tested. Bayesian approaches are commonly applied for modeling crop responses (Besag & Higdon, 1999; Huang, Huang, Yu, Ni, & Yu, 2017; Jiang, He, Kitchen, & Sudduth, 2009; Lawrence et al., 2015; Nandram, Berg, & Barboza, 2014; Ramsey, 2020; J. C. Wang et al., 2012) but are also used for models that influence policy such as yield gap and potential (Prost, Makowski, & Jeuffroy, 2008) and commodity prices (Drachal, 2019). The first Bayesian model used was adapted from Lawrence, Rew, & Maxwell (2015). Two modifications were made to fit this dataset and scenario, including replacing electrical conductivity data with clay content, based on the correlation between the two properties (McBratney, Minasny, & Whelan, 2005). Second, the precipitation term used as a predictor for crop responses was constrained in this model to precipitation summed across the previous year, as only this data would be available when forecasting is required to make management decisions in an upcoming year. Finally, because this data represented scattered points across a field rather than a lattice grid structure in Lawrence et al., simultaneous spatial auto regression (SAR) was used rather than conditional auto regression (Ver Hoef, Hanks, & Hooten, 2017; Wall, 2004). The base of this model was identical to Lawrence et al. (2015) where grain yield and protein were assumed to follow a normal distribution;

(4)

Where R was the grain yield (kg ha^-1) or grain protein concentration (%) that followed a normal distribution with a variance, σ_ε², that followed an inverse gamma distribution;

(5)

 And the mean crop response, μ, was defined by a non-linear crop model;

(6)

Where N was the as-applied N-rate (kg ha^-1), prec_py and claycontent are defined in Table 3, and z was the SAR spatial random effect;

(7)

Where η had a mean E(η) ≡ 0 and Cov(η) ≡ σ_z²I. I was the identity matrix, and B was the n x n spatial dependence matrix between n_i and n_j  where i ≠ j and i, j = 1, ..., n  (Hooten, Ver Hoef, & Hanks, 2019).

All priors used were the same as in Lawrence et al., where the prior for σ_ε² followed a gamma distribution with 0.01 used for both shape parameters in (5). For β terms in (6), truncated normal distributions, TN(0, 1000), were used for all terms except where β_max = N(0, 1000). All Bayesian methods were executed in R using the Stan based brms package (Bürkner, 2021; 2018; 2017).

4) Bayesian Multiple Linear Regression Model (Bayes MLR)

The second type of Bayesian approach tested was a multiple linear regression (MLR) model. Based on initial exploration of the data (Hegedus & Maxwell, 2022), a normal distribution as the base of a SAR MLR model was used (Anselin et al., 2004; Hooten et al., 2019; Kissling & Carl, 2008; Ver Hoef, Petersen, Hooten, Hanks, & Fortin, 2018);

(8)

Where R was the grain yield (kg ha^-1) or grain protein concentration (%), σ_ε² was identical in formation and priors as (5), and;

(9)

where β was the coefficient vector (β₀, ..., β_p-1) with p the number of predictors from Table 3, X is the n x p covariate matrix of covariate vectors, and z was the SAR random effect. To avoid overfitting, top-down feature selection was applied and omitted covariates where dropping them from the model reduced the AIC by more than two units. Thus, Bayesian MLR models vary in their covariates between fields and splits but are all fit from the same initial model;

(10)

Where predictor variables are covariates from Table 3 with uniform priors on coefficients β₀ - β₁₆ (Gelman, 2006), and the SAR spatial random error, z, defined in (7). As with the non-linear modeling, the Stan based brms package in R was used for Bayesian analyses (Bürkner, 2021; 2018; 2017).

5) Random Forest Regression (RF)

Random forest (RF) regression is an increasingly popular non-parametric machine learning method that has been used to model spatial data (Georganos et al., 2021; Jing, Yang, Yue, & Zhao, 2016; Li, Heap, Potter, & Daniell, 2011; Shi et al., 2015), specifically crop response data (Everingham, Sexton, Skocaj, & Inman-Bamber, 2016; Jeong et al., 2016; Lawes, Oliver, & Huth, 2019; Mariano & Mónica, 2021; Marques Ramos et al., 2020; Paccioretti et al., 2021; L. Wang, Zhou, Zhu, Dong, & Guo, 2016), where crop responses are estimated by the ensemble of regression trees using binary splits of observations based on covariate data. The ranger package in R was used for fitting and generating predictions (Wright & Ziegler, 2017). The predictors included in the RF models for yield and protein are in Table 3. To avoid over-fitting, top-down feature selection was performed where covariates were omitted if dropping them from the model resulted in an AIC decrease of over two units. The number of trees (ntrees) and the number of covariates sampled at each node (mtry) were optimized during the fitting process. To account for spatial effects, UTM coordinates were included as a variable definition following spatial random forest regression practices in the literature (Janatian et al., 2017; Langella, Basile, Bonfante, & Terribile, 2010; Walsh, Kreakie, Cantwell, & Nacci, 2017; Y. Wang et al., 2017).

NUE Sampling and Calculation

Table 5. Collaborator farm ID (Fig. 1), field names, field size (hectares), crop harvest history, and the year sampled. The letter identifier for the Farm corresponds to Figure 1.

The time and money required for detailed soil sampling are limitations on the ability of most farmers to assess N transport and soil characteristics within their fields. Based on the work of Sigler and John et al. in Montana small grain dryland agroecosystems, the thickness of fine textured soils, or simplified to soil water storage, is assumed to be a soil characteristic that influences variation of NUE within fields. If precipitation inputs exceed crop water uptake and soil water storage, nitrate leaching is hypothesized to be the dominant pathway of N loss because nitrate that has not been taken up by the crop or retained in the soil will be transported below the crop rooting depth in well-drained soil. If N limitation of dryland small grain agroecosystems is a function of leaching loss with excess N application, then soil water storage (WHC) will dictate yield capacity when N fertilizer application exceeds crop requirements. Thus, soils with high WHC will contribute to higher yields and less nitrate leaching of excess fertilizer because the water is retained in the root zone.

To account for the variation in NUE explained by soil water storage, a high spatial density estimate of WHC was created. If soil water storage is a function of the depth of fine textured soils, then consistency in productivity over time at a specific point in a field could reveal insight into WHC at that point. While no studies have found have linked historical NDVI imagery to NUE or WHC, Araya et al. used NDVI to estimate available water content, indicating the promise NDVI has for estimating soil properties (Araya, Lyle, Lewis, & Ostendorf, 2016). The implications for identifying areas of a field at risk of N loss from freely available remotely sensed data benefits farmers by providing information on where NUE is potentially low, which can inform site-specific N fertilizer management decisions.

The fields used in this study were classified into WHC zones that were used to stratify soil and plant biomass sampling. Normalized Difference Vegetation Index (NDVI) was derived from satellite imagery from 2000 – 2021. Landsat imagery was collected on 30m x 30m raster grids, a coarser resolution than satellites such as Sentinel-2 however, was chosen for the length of imagery history available (Table 6). Landsat 7 data was only used for one year due to a failure in the Scan Line Corrector causing a gap-filled zig-zag pattern during data collection.

Imagery was downloaded from Google Earth Engine (GEE), where image processing was performed to correct images for clouds and other atmospheric anomalies (Gorelick et al., 2017). The cloud computing platform of GEE was taken advantage of to calculate NDVI as;

(11)

where ‘NIR’ is the near-infrared wavelength (~865 nm) and ‘Red’ correspond to red light wavelength (~665 nm). Landsat data was downloaded for each of the boundaries of the individual farms, from which images were clipped to the boundary of each field. Although data was collected from 2000 – 2021, years in which average NDVI across each field fell below 0.3 were omitted because these represent fallow years in which no crops were grown, and any productivity was due to weeds (Fig. 3).

Fig. 3. Example of the data reduction process in Field 2 where fallow years were omitted based on mean NDVI across the field for a given year. Fallow years (mean NDVI < 0.3) are shown in the figure on the left, and mean NDVI without fallow years is shown on the right. While these figures show an example for one field, this process was repeated for all fields sampled.

For each field, the rasterized NDVI for each year was stacked, and WHC was classified based on each pixel’s mean and coefficient of variation (CV) to produce a single map with classifications of high, medium, and low WHC. High, medium, and low categories for mean and CV for each pixel were delineated into the three groups via the Jenks optimization method for dividing the data (Jenks, 1967). A matrix the mean and CV groups for NDVI were then derived to estimate hypothesized WHC zones (Table 7). If an area of a field is consistently productive over time, the crop is acquiring adequate water and N to maximize growth across varying annual precipitation. If the water is held in the soil, then so too is nitrate, so these areas are hypothesized to have a high WHC and a low risk of N loss. If the inverse is the case and a crop is consistently unproductive over time, the crop does not acquire enough water and N to maximize productivity under varying annual precipitation. The WHC at these points is expected to be low due to a shallow depth of fine textured soils, resulting in a high risk of N loss and low NUE.

Table 7. Classification scheme of WHC zones based on mean NDVI and the CV. Each pixel is classified based on the mean and CV across the 20 years of data collected.

Soil samples were taken at 75 locations in each field twice per year, in March and August, randomly stratified on the classification for WHC and the experimentally varied N rates (Fig. 4). Soil samples were taken at depths of 0-15, 15-30, 30-60, 60-90, and from 90cm to as deep as possible. The depth of the deepest sample corresponded to either gravel content, rock layer or the limit of the probe (110 cm). The maximum depth at each location was recorded during sampling. Soil samples were obtained by a truck mounted hydraulic probe in March (pre-fertilization) and August (post-harvest). Soil samples were taken within 7 days prior to fertilization and within 7 days post-harvest. Winter wheat biomass samples were hand harvested in June at plant maturity at each of the sampling locations. Samples were transported in coolers on ice from the fields back to MSU where soil and plant samples were weighed wet, oven dried at 50 degree C for 7-14 days, and then weighed dry to determine water content. Plant samples were ground using an UDY Cyclone Sample Mill (UDY Corporation, Fort Collins, CO, USA). Soil samples were milled to break up peds and aggregates and sieved to remove rocks and roots to homogenize samples with a Dynacrush DC-5 soil crusher (Custom Laboratory, Hodlen, MO, USA). Aliquots were taken from each sample (plant and soil) and analyzed for total N with a combustion analyzer (Costech analytical technologies Inc., Valencia, CA, USA) in the Environmental Analytical Laboratory at MSU. Costech measurements for soil samples from March and August, Costech measurements for plant biomass samples, and as-applied N fertilizer rates were co-located by sample point, along with the derived WHC zone.

Fig. 4. Water holding capacity classifications and sampling locations for each field.

A partial mass balance approach was used to determine an estimate for efficiency for each sample location. Measured soil percent TN and bulk density were used to calculate soil TN, which was converted to kg ha^-1. Crop N uptake was derived from the percent TN in the plant aboveground biomass and the closest crop yield measurement, typically a maximum distance away of 5m. Nitrogen taken up by the plant or left in the root zone after harvest were regarded as economically beneficial (efficient) uses of N because N in the crop contributes to net-return in the season the crop was grown, while N left in the root zone has the potential to be taken up by future crops and contribute to net-returns for farmers in upcoming years. However, knowledge of NUE requires knowledge about more aspects of the N cycle. Nitrogen use efficiency was estimated using soil samples prior to fertilization and post-harvest, plant biomass samples at plant maturity, and estimates of other inputs to derive the NUE. Nitrogen use efficiency was thus defined as outputs deemed efficient over the total amount of N inputs;

(12)

where N_crop is N taken up by the crop measured from plant biomass samples, N_t is the amount of N in the soil after harvest measured by soil sampling, N_t-1 is the amount of N in the soil prior to N fertilizer application in the spring measured by soil sampling, and N_inputs are;

(13)

where N_dep is N deposited from the atmosphere, N_min is mineralized N, N_fert is the inputs of N fertilizer gathered from the farmer’s as-applied N fertilizer maps, and N_fix is N fixed from the atmosphere. Due to the nature of the crop-fallow rotation there was no N fixation, so N_fix was assumed to be zero. Total N deposition was estimated to be 2.92 kg N ha^-1 year^-1 based on the mean total N deposition from 2000 – 2018 recorded at the closest CASTNET station at Glacier National Park (U.S. EPA, 2018). Pathways recognized in which N can be made unavailable to plants are volatilization, leaching, immobilization, runoff, and subsurface lateral flow, however the focus of this work was on quantifying how much N was partitioned to productive outputs rather than quantifying the amount of N lost via each output pathway.

However, the equation of NUE contains one unknown, the N production via mineralization within fields (N_min). Due to the limitations of the mass balance approach measuring total N, mineralized N was modeled using a regression model derived from a review of multiple mineralization studies on the relationship between mineralization rate and percent total nitrogen (Vigil, Eghball, Cabrera, Jakubowski, & Davis, 2002);

(14)

Where N_min is mineralized N in kg ha^-1, and TN is the percent total nitrogen in the soil. Mineralization of soil organic N (SON) to plant available forms is a function of organic matter, moisture, and temperature, however equation 14 uses only percent soil total N to estimate mineralized N, introducing uncertainty into mineralized N predictions (Vigil et al., 2002). Soil organic N is typically most prevalent in the upper surface of soils, thus mineralization occurs predominately in the soil surface (Cassman & Munns, 1980). Crop N uptake is also greatest in the surface of soils due to mineralization of N and the surficial application of fertilizer N and deposition of atmospheric N. Due to most processes that contribute to NUE occurring predominately in the surface soil, only soil samples from the top 30cm were used to calculate and model NUE. A final equation for calculating NUE of the top 30cm of soil was;

(15)

where all terms are defined above. This equation was used to calculate NUE at all 75 points within each field for both years. Sample points that did not have corresponding measurements at each depth in both March and August were omitted, yielding observations at 67, 69, 60, and 75 out of 75 sample locations in BS 2020, BSE 2021, WH 2020, and MS 2021, respectively.

NUE Modeling

To extrapolate estimates of NUE to other portions of a field and other dryland winter wheat fields, a generalizable model for estimating NUE at a subfield scale was developed with low-cost covariate data. The calculated NUE from equation 15 was the response variable used for developing a generalized model of NUE for dryland, small grain, conventionally managed Montana farms. The model intensions were to estimate subfield NUE for a given year without requiring producers to take detailed soil samples or provide parameters for complex N transport models.

Development of a general model for predicting subfield NUE relied on open source remotely sensed covariates gathered from Google Earth Engine (GEE), which were assumed to influence NUE (Table 2). All remotely sensed data were gathered from open-source GEE data repositories (Gorelick et al., 2017) and aggregated with soil sample data. Information from each covariate was extracted to each soil sample location to form the analysis dataset (Hegedus & Maxwell, 2022b; 2022c). Temporal data were gathered and delineated by year and constrained in range to the point in time that a farmer needs to make decisions on N fertilizer management of dryland winter-wheat in Montana (Hegedus & Maxwell, 2022c). Although harvest data were collected for each field, these data would not be available to make predictions at the point of a decision on N fertilizer management and were not included as covariate information for modeling NUE. The covariates selected for modeling were based on assumptions about environmental variables that influence the components of NUE and observed relationships from evaluation of observed NUE (Fig. 5). The green circles are measurements directly used to calculate NUE, besides atmospheric deposition of N. The lines leading from the green circle to NUE (blue) are thickest to illustrate the strength of the relationship between these variables and the actual NUE indicating that NUE predictions would be most accurate if these variables were measured directly. However, soil sampling is time and money intensive, so modeling NUE using data available to farmers increases a farmer’s ability to use data about their fields to make efficient management recommendations for their fields. Orange circles represent types of data relevant to each component of the calculation for NUE (green). Note that soil characteristics, topographical features, and weather all influence both mineralization of nitrogen and crop uptake at any location within a field. Grain yield and grain protein concentration are gathered from farms, but as farmers needs to make recommendations on nitrogen fertilizer in the spring prior to harvest, these data would not be available at a farmer’s decision point and not used in the model. The yield and protein circles are colored black because this data is not available, however yellow circles indicate remotely sensed covariate data that may benefit predictions of NUE (Fig. 5).

Figure 5. Conceptual model of the development of a general NUE model for Montana dryland agroecosystems. Green circles are variables that are used to directly calculate NUE, thus measuring these variables would provide the best characterization of true NUE. Orange circles are general data types that influence the components used to calculate NUE (green). Yellow circles are data types organized by general data type (orange), that are measured from remote sensing sources. The thickness of connecting arrows represents the strength of the relationship between the variable and NUE, where a variable is more directly related to NUE with thick arrows and more indirectly related with a thin arrow. The strength of a relationship between data and NUE decreases as distance between data and the blue NUE circle increases. Yield and protein would be beneficial data types for predicting crop N uptake and subsequent NUE but are not available data when a farmer models NUE to derive optimal N fertilizer rates.

A set of five candidate models were tested on their ability to predict NUE. Candidate models were chosen based on their potential for predictive and not inferential modeling. Model candidates were selected to represent simple to complex approaches with varying degrees of flexibility and assumptions. The models tested were 1) a baseline simple linear regression (SLR) model that only included as-applied nitrogen as a covariate, 2) a multiple linear regression model (MLR) using available covariates, 3) a non-linear model using an exponential decay function with as-applied N (NLR), 4) a generalized additive model (GAM), and two machine learning approaches, 5) a random forest regression model (RF) and 6) a support vector regression model (SVR). The exponential decay function was selected based on exploratory analysis of the relationship between NUE and N sources. Random forest regression is an ensemble tree approach where covariates are randomly sampled in each tree and used to classify observations (Breiman, 2001). Support vector regression maps observations to the feature space generated by the covariate data (Lau & Wu, 2008; Smola & Schölkopf, 2004).

Calculated NUE exhibits spatial dependence between neighbors, as NUE values are more similar for observations near each other than for observations with a greater distance between them. Spatial covariance structures were incorporated into the models in varying ways, from a Gaussian process spatial correlation structure used in the SLR, MLR, and NLR, (De Bastiani, Mariz de Aquino Cysneiros, Uribe-Opazo, & Galea, 2015; Diggle & Tawn, 1998; Mardia & Marshall, 1984; Xia, Ding, & Wang, 2008), a Gaussian basis function on spatial coordinates in the GAM (Gotway & Stroup, 1997; Guisan, Edwards, & Hastie, 2002; Holland, De Oliveira, Cox, & Smith, 2000; Zuur & Camphuysen, Kees, 2012), to common machine learning approaches for spatial data with the RF and SVR (Janatian et al., 2017; Langella, Basile, Bonfante, & Terribile, 2010; Walsh, Kreakie, Cantwell, & Nacci, 2017; Wang et al., 2017).

All models except the SLR and NLR were subjected to feature selection during the fitting process to optimize model performance. Top-down AIC based feature selection was performed for the MLR and GAM models, where a reduction of two AIC units with the removal of a covariate justified its omission. As AIC is not appropriate for machine learning methods, top-down RMSE based feature selection was performed for the RF and SVR models. Feature selection was based on any reduction in RMSE to determine whether withholding a covariate improved model performance. All models, except the SLR and NLR that only use as-applied N fertilizer, began feature selection with a full model that included the WHC classification, and the covariates illustrated in the yellow circles of Fig 5.

The initial candidate models were subjected to a hold-one field out cross validation approach to initially assess the ability of the model to predict NUE. Each model was independently trained on data from three out of four observed fields, and performance was tested by calculating the root-mean square error (RMSE) between the predicted and observed NUE in the withheld field (Fig. 6). The RMSE across the four test datasets were averaged for each model and compared to identify the two candidate models that generated the lowest mean RMSE. The two models selected from the hold-one field out cross validation approach were then empirically assessed using 5x2 cross validation (Dietterich, 1998). All the data from the fields were pooled and randomly split 50/50 into a training and test dataset. Both models were trained on the training dataset and RMSE was calculated based off predicted and observed NUE in the test dataset. Then the training and testing datasets were swapped (2 folds), and the process was repeated to generate four error metrics and calculate the variance of the difference in the predictions. Beginning with splitting the data 50/50, the process is repeated five times (5 folds), and the error estimates and variance from the five repetitions were used to calculate the 5x2 t-statistic (Dietterich, 1998). Assuming a t-distribution with five degrees of freedom, evidence against the null hypothesis that there was no difference in RMSE between observations and predictions from the two models was evaluated with a two-tailed significance test (α = 0.05).

Figure 6. Hold-one field out cross validation scheme where one field is held out as a test dataset while three fields are used to train models in each of the four cases.

Optimized Site-Specific N Fertilizer Rates 

To derive optimized site-specific N fertilizer recommendations, field-specific models for predicting yield and protein were used in simulations to derive optimum rates under economic uncertainty. Simulation was used to develop site-specific N management that accounts for the tradeoff in net-return and NUE (SSOPT). A 10m grid was generated for each field, to which remotely sensed data were aggregated to the grid centroids to create the simulation datasets. The field-specific crop responses models trained on observed data were then used to predict grain yield and grain protein at every point in each simulated dataset for every N rate from 0 to 168 kg N ha^-1 .

At every N rate for every point, the yield and protein predictions were used to calculate net-return;

(16)

where NR is the net-return ($ ha^-1) received and a function of the product of the yield (kg ha^-1) and the price received (P), minus the cost of the applied input (CA) multiplied by the as-applied rate of the input (AA), the fixed costs (FC) associated with production ($ ha^-1) that do not include the input, and the cost per hectare of the site-specific application (ssAC). As Montana farmers receive a premium/dockage to the base price received for wheat based on grain protein concentration, the price received for calculating net-return was (Hegedus & Maxwell, 2022b; 2022c);

(17)

where P is the final price received ($ kg^-1), Bp is the base price received ($ kg^-1), B0pd is the intercept of the protein premium/dockage function set by the grain elevator, B1pd is the coefficient on the grain protein concentration (protein, %) and B2pd is the coefficient on the squared protein term.

In addition to simulating the net-return, NUE was predicted at every point for every N rate to account for the tradeoff in net-return and NUE. The optimum nitrogen fertilizer rate for any given point in the field was then identified as the rate that maximized profit and efficiency. Under the assumption that nitrogen use efficiency cannot exceed 100%, 1 – NUE represents the inefficient proportion of nitrogen available for crops. Desiring maximization of net-return and minimization of nitrogen use inefficiency, the optimum N fertilizer rate at each point was the rate at which the distance between net-return and nitrogen inefficiency was greatest (Fig. 7). The FFOPT net-return was based on the N rate that maximized net-return and minimized inefficiency when it was uniformly applied across the field.

Fig. 7. Conceptual diagram of the optimization process for finding the nitrogen fertilizer rate at each point that maximizes net-return and minimizes nitrogen use inefficiency, calculated as 1 – nitrogen use efficiency.

Uncertainty in economic conditions for an upcoming year was addressed through a Monte Carlo simulation where economic conditions from previous years were randomly selected at each iteration and used to calculate net-return (Hegedus et al., 2022). Management recommendations for the SSOPT method are produced for farmers by taking the mean optimized N rates from every point in the field across the simulation. As farmers are also uncertain in the weather conditions of upcoming years, optimized N fertilizer rates were developed under four different weather conditions, where the year from 2000 – 2021 representing the most recent year OFPE was applied (Recent), most average precipitation year (Average), driest (Dry), and wet (Wet) were simulated for each field.