Introduction

Soil spatial variability has been studied for several decades (Burrough, 1993; Heuvelink & Webster, 2001), and it has been established that within-field variability is a primary source of variable crop yield response across the field (Earl et al., 2003). Field-scale variability in soil properties and processes greatly affect soil moisture and nutrient availability, which influence crop growth and development, and ultimately impacts crop yield (Warrick & Gardner, 1983). Soil spatial variability influences soil microbial diversity (Juma, 1993; Naveed et al., 2016) and soil health status (Adhikari et al., 2021; Zebarth et al., 2019) leading difference in yield across a given field provided uniform inputs (e.g., fertilizer, seeding rate). Moreover, interactions among yield limiting factors such as soil, topography, and weather and management are also responsible for non-uniform yield across a field (Kravchenko et al., 2005). Precision agriculture (PA) is a tool that may be used to manage soil and crop spatial variability by applying inputs in a manner that is site-specific to optimize yield, environmental benefits, and sustainability (Adhikari et al., 2009; Bongiovanni & Lowenberg-Deboer, 2004; Oliver et al., 2013). Application of PA principles to soil and water conservation at farm or landscape scale is even more pronounced when optimization in conservation practices and farm profitability are desired, and this can be achieved through precision conservation application (Capmourteres et al., 2018; McConnell & Burger, 2011). Positive effects of field-scale application of conservation practices has long been reported (Her et al., 2016).

Variability in crop yield occurs spatially and temporally, but only with several years’ worth a geospatially explicit yield data can these be fully understood. However, temporal yield patterns are not always spatially consistent over years; part of a field could produce more in one year but not the same amount in a different year (Cox & Gerard, 2007). Porter et al. (1998) reported greater interannual variability than spatial variability based on 10-year analysis of plot-scale corn yield. Interannual yield variability is influenced by several factors and therefore, a single-year yield information is insufficient for PA decisions (Eghball & Varvel, 1997). For example, potential management zones derived from a single year yield data are not reliable, since the data likely do not reflect underlying soil-yield relationship properly, as this relationship is highly dependent on weather and management (Cox et al., 2006).

Spatio-temporal yield variability is visualized with yield maps derived from point yield measurements, and yield mapping is one of the primary goals of PA applications. Point yield measurements can be upscaled to a continuous yield map using different interpolation techniques such as inverse distance weighting, splines, natural neighbor, and geostatistical techniques that are most common in PA (Oliver, 2010). Ordinary kriging (OK), one of the widely used geostatistical techniques in PA (Frogbrook et al., 2002; Oliver, 2010), relies on a precise characterization of autocorrelation among observations using a variogram. A variogram measures the mean degree of dissimilarity (semivariance) between samples over locations, and thus can describe autocorrelation among observations at specified distances (Deutsch & Journel, 1999). In general, the semivariance increases with increasing distance between observations until an upper bound is achieved. The value of semivariance at which the variogram plateaus is called the sill, the distance at which the sill is reached is called range, and the semivariance at zero distance is called nugget. The nugget is an estimate of the residual error or spatially uncorrelated error. The observations located within the range are spatially dependent or autocorrelated, while observations away from this distance are not.

Repeated measures of yield data help producers and crop advisors identify production trends and accumulate information about temporal yield stability. Cropping systems are stable when the variance component of yield is small. Producers prefer unstable (large variance) high mean yield compared to stable (small variance) low yield (Piepho, 1998). Thus, yield stability could be identified using mean yield and coefficient of variation (CV) (Francis & Kannenberg, 1978; Piepho, 1998). Maps delineating zones within a field that produce either low or high yield (stable zone) over the production period separated from unstable zones is one option to evaluate precision management. Long and Ketterings (2016) divided the field into yield stability zones or quadrants (Q1, Q2, Q3, and Q4) by using mean yield and CV of multi-year yield data. Based on this zonation, Q2 and Q3 would be unstable zones due to higher CV, whereas Q1 and Q4 would be stable zones with lower CV. This approach has been successfully tested in other studies as well, for example, Kharel et al. (2019) used mean and standard deviation of temporal yield data to delineate yield stability zones in New York dairy farms. Kharel et al. (2019) also evaluated whether such division is reasonable using percent yield variability explained by each zone and reported that the developed zones alone could explain 35–56% yield variability.

In PA, management decisions are implemented and evaluated based on gross margin from the fields, and developing a proper enterprise budget and record-keeping are important tools for efficient farm management. Gross margin can be calculated as the difference between variable and fixed cost (cost related to seeds and fertilizer, chemicals, equipment, etc.,) and income (revenue collected from selling the product) (Nix, 1995). Developing a gross margin map of the field is more desirable in PA or precision conservation applications as the map gives monetary value from each small unit within the field, directly facilitating easy comparisons between fields and yield stability zones. Such maps can be used to identify potential areas to set-aside within each field using economic justification (Blackmore, 2000), thus providing potential precision conservation applications. The primary goal of this study is to investigate within-field yield stability zones and gross margins within those zones as a tool to implement precision conservation decisions. Specific objectives were to (i) map within-field corn yield variability using geostatistics, (ii) delineate yield stability zones using spatial–temporal yield data, and (iii) map and evaluate gross margin across fields and zones for precision conservation decisions.

Material and methods

Study site

The study was conducted in nine corn fields located near Riesel, Texas (31° 28′ 16.2″ N, 96° 53′ 06.5″ W) (Fig. 1, left), which are managed by the USDA-Agricultural Research Service, Grassland, Soil and Water Research Laboratory. The fields are representative of fertile Blackland Prairies in north central Texas, and cover a total area of 56.3 ha, where the largest field is of 9.55 ha (Field 16-A), and the smallest field is 3.86 ha (Field SW-16). Soils in these fields are predominantly grouped into Houston Black clay Soil Series and classified as Vertisols (fine, smectitic, thermic, udic Haplustert) according to Soil Taxonomy (Soil Survey Staff, 2010). These smectitic shrink-swell soils have a typical particle size distribution of 17% sand, 28% silt, and 55% clay with low permeability when wet, but show high infiltration rates when dry due to preferential flow associated with soil cracks (Arnold et al., 2005). Prior to 2017, these fields were managed differently with corn, wheat, and oats as main crops in rotation, and leguminous cover crop as green manure or for harvest as hay. After 2018, all fields were planted with corn as a main crop with cover crop, and reduced tillage in selected fields. Reduced tillage, and contour terracing are the major conservation practices applied in these fields. The site is characterized by a humid subtropical climate with an average annual temperature of 19.6 °C, and average annual precipitation of 940 mm, respectively.

Fig. 1
figure 1

Location of investigated fields in Riesel, Texas (Left); corn yield measurements (A–C) and predicted maps (D–F) from the three selected fields Y-6, Y-8, and W-13 for year 2018, 2019, and 2020, respectively

Full size image

Corn yield measurements

Corn yield was recorded for three consecutive years (2018, 2019, and 2020) using a combine harvester equipped with an Ag Leader yield monitor system. Harvest equipment has a swath width of about 6 m that harvests 8 corn rows and the average speed during harvesting was about 6.5 km h−1. At each measurement location, geographic coordinates were recorded using Differential Global Positioning System. Raw data from the harvester were cleaned using Yield Editor software (Sudduth & Drummond, 2007) to correct, for example, for flow delays or slow combine velocity at the beginning and end field passes. All statistical and geospatial data analyses were performed on cleaned data.

Field economic data

Field economic data were collected for each field during 2018, 2019, and 2020 crop period to maintain farm budget for cost–benefit analysis. Detailed farm management data for corn including farm operations such as tillage, planting, harvesting, pest and weed control, crop yield and selling price, fertilizer, seed, and chemical purchase records, and fuel, machinery and labor were properly maintained by the farm. As there was no variable rate technology applied, all these records were reported as uniform application across individual fields, and were tabulated as per field or per hectare basis.

Statistical analysis

Statistical analysis included data cleaning and performing descriptive statistical analysis of yield data measurements from 2018, 2019, and 2020 for each field together with zonal statistics, and data visualization for graphics and reporting. These analyses were performed in JMP (SAS Institute Inc., 2009), and in R software program (R Development Core Team, 2008).

Geospatial modeling

Yield mapping

Georeferenced point yield measurements from 2018, 2019, and 2020 were used to generate continuous corn yield maps of each field using OK. An omni-directional experimental variogram of yield measurement with 20 lags with varying lag distance depending on the field was computed to which three theoretical variogram models (Spherical, Exponential and Gaussian), were fitted in VESPER program (Minasny et al., 2005). Among the three models tested, the best fitted model was selected based on Akaike Information Criteria with the lowest value for the best model. The variogram parameters of the best fitted model were than used to map yield across the field at a grid resolution of 5 m × 5 m in ArcGIS (ESRI, 2012). Unlike in the Spherical model where the variogram reaches the sill, the range of Exponential and Gaussian model was approximated at which the variogram value was 95% of the sill because the latter two models reach its sill asymptotically (Isaaks & Srivastava, 1989).

Delineating yield stability zones

The kriged maps from the study period for each field were used to generate yield stability zones. The procedure included the following steps:

  1. (1)

    Derive the map of mean, and CV from 2018, 2019, and 2020 yield maps,

  2. (2)

    Export the mean, and CV values for each pixel, and plot them against each other in a scatter plot,

  3. (3)

    Using the mean of yield, and CV as reference values, divide the pixels into quadrants or zones,

  4. (4)

    Map the zone as representative of each quadrant as yield stability zone.

Figure 2 shows an example on how the stability zones were derived for field Y-8, and changes in yield and CV across zones. Mean yield, and CV in this field during the study period were 5.8 Mg/ha and 40.1%, respectively, and these values were used to divide the field into zones. Zone A corresponds to those pixels having yield higher than the mean but CV lower than its mean; Zone B is where the pixels have yield and CV values higher than their corresponding means; Zone C has CV higher than mean but yield lower than mean, and Zone D has both values lower than corresponding mean values.

Fig. 2
figure 2

Dividing the fields into yield stability zones (Zone A, Zone B, Zone C, and Zone D) using mean yield and coefficient of variations (CV) for 2018, 2019, and 2020 (left) for a selected field Y-8, and corresponding yield and CV values in four zones (right). The vertical and horizontal dotted line in the left figure represent the corresponding mean value, and the dot in the box plot in the right figure is the mean

Full size image

Assessment and mapping of gross margin

Gross margin from each field was calculated by deducting fixed and variable cost from annual revenue, and the map of gross margin of a field was compiled by calculating it for each grid within the field. It was assumed that the income from a field varies across the field due to changes in yield but the fixed and variable costs remain the same since inputs were applied uniformly across each field. The gross margin for each grid was calculated as Eq. (1) and was mapped for the entire field in ArcGIS (ESRI, 2012):

$${GM}_{i}=\left({Y}_{i}\times CSP\right)-(VC+FC)$$
(1)

where GMi ($/ha) is gross margin at grid i, Yi (Mg/ha) is the yield from that grid, CSP ($/Mg) is a crop selling price, VC ($/ha) and FC ($/ha) are the variable cost and fixed cost for a field. The Yi × CSP is the income from selling the grain, VC covering the cost related to seed, fertilizer, chemicals, fuel, machinery labor, and repair and maintenance, FC is for equipment investment and machinery depreciation per field. If any treatments such as seed, fertilizer, or herbicides in the field varied spatially, VC and FC should be calculated for each grid. Similarly, if the CSP is not uniform, for example, if it changes with grain quality, it should also be calculated for each grid.

The gross margin map for each field was then intersected with corresponding stability zone map to calculate and compare gross margins across stability zones, and assessment results were then used to develop a precision conservation plan. Among different conservation plans that could be applied, we plan to take these unstable and less productive areas out of crop production and possibly transfer them into some type of permanent cover.

Results and discussion

Yield measurement

Descriptive statistics for corn yield from nine fields during the study period are listed in Table 1. In general, mean yield in 2018 was lower than 2019 or 2020 in all fields. Field Y-8 gave the highest yield of 2.8 Mg/ha (± 0.6 SD) for 2018 followed by field 6–12 (2.4 Mg/ha ± 0.5 SD), and Y-10 had the lowest yield of all (1.8 Mg/ha ± 0.5 SD). During the study period field Y-6 had the greatest variability (CV = 44.7% in 2019), and field W-13 had the least variability (CV = 16.1% in 2020). Comparing yield between 2019 and 2020, fields 6–12, W-12, Y-13, and Y-8 had greater yield in 2020, but fields 16-A, Y-6, Y-10, and W-13 had greater yield in 2019. However, yield from SW-16 was similar for 2019 and 2020, but the CV was higher (41%) for 2019 than 2020 (27.8%). In general, the northern part of the field Y-6 and eastern part of Y-8 had lower yield compared to the rest of the field. For the field W-13, the southern part had lower yield than the rest of the field, specially in 2019.

Table 1 Descriptive statistics of corn yield measurements and predicted maps for 2018, 2019, and 2020 [SD standard deviation, CV coefficient of variation]
Full size table

Yield maps

Table 2 summarizes variogram parameters used in OK mapping. Among the three variogram models tested, Spherical and Exponential models were the best models to characterize yield measurements autocorrelation, which is common as also reported in Frogbrook et al. (2002). The nugget variance which represents the unstructured variance that was not modeled by the variogram ranged between 0.04 and 2.40, respectively, from field Y-8 and 6–12 in 2019. A maximum range (418 m) of the variogram was reported from Y-8 in 2020, and minimum (20 m) from SW-16 and W-12 in 2019, and 2020, and these range values determine the extent of patches created due to autocorrelation in the predicted maps (Frogbrook et al., 2002). Similarly, a maximum sill of the variogram was reported for Y-6 in 2019, followed by 2020 for the same field, while the lowest sill of all was observed for W-13 in 2018.

Table 2 Variogram parameters of corn yield mapping using ordinary kriging
Full size table

Overall, all fields had lower yield in 2018 compared to the yield from 2019 and 2020, the highest yield (11.7 Mg/ha ± 1.4 SD) in all years was mapped from W-12 in 2020, and the lowest yield (1.7 Mg/ha ± 0.3 SD) was from Y-10 in 2018. Year 2018 was the driest year of all directly impacting the yield (Adhikari et al., 2021). Figure 1, D to F display the yield maps of three selected fields Y-6, Y-8, and W-13 as a reference. Mean yield from Y-6 and W-13 was greatest in 2019, whereas Y-8 had the greatest yield in 2020 compared to the yield from 2018 or 2019. While looking at the yield distribution during the three years in each field, a comparable yield distribution pattern was observed in almost all fields. For example, the northern part of Y-6 and eastern part of Y-8 had lower yield compared to the rest of the fields, with the pattern being more evident for 2019 and 2020. For W-13, the southern part of the field, in general, had lower yield compared to the rest of the field for all years.

From the yearly yield map of 2018, 2019, and 2020 in each field, mean yield and associated CV were calculated, and the resulting maps are displayed in Figs. 3 for the selected fields Y-6, Y-8, and W-13. Overall mean yield from field 16-A was 5.11 Mg/ha (± 0.9 SD), where the mapped values ranged from 1.20 to 7.82 Mg/ha. Field W-12 had the greatest overall mean yield (8.17 Mg/ha ± 0.9 SD) among all fields, and field Y-10 had the lowest overall mean yield (3.64 Mg/ha ± 0.4 SD). Similarly, overall CV for the entire study period ranged from 39.0 (Y-10) to 53.9% (W-12). The maps also revealed some distinct yield patterns during the study. For example, the drainage pattern from a draw in the northwest portion of field 16-A presented with low yield and low CV. The northern part of Y-10 has terraces running east–west with alternate high and low yield values but with higher CV. Similarly, the eastern part of Y-8, and northern part of Y-6 are known to have seeps, which produced low yield accompanied with low CV from those areas. This showed that some parts of the field were stable in terms of its production-either high or low yield-while the other parts were very unstable during the study period.

Fig. 3
figure 3

Average corn yield from the three selected fields for 2018, 2019, and 2020 (A), and associated coefficient of variation (B)

Full size image

Yield stability zone

Yield stability zones were identified using mean and CV of yield for the duration of the study from each field (Fig. 4). Few other studies such as Blackmore et al. (2003) and Kharel et al. (2019) have used SD instead of CV for identifying yield stability zones. However, Francis and Kannenberg (1978) and Piepho (1998) preferred CV to represent yield variability as applied in this study. We identified and delineated four distinct zones in each field with reference to how stable the yield was during those 3 years. Zone A and Zone B generally produced a greater yield than the mean in which Zone B was unstable due to its higher CV values relative to Zone A. Both Zone C, and Zone D produced less than the field mean yield, but Zone D is more stable due to lower CV. So, Zone A and Zone D are classified as stable zones, whereas Zone B and Zone C are unstable zones (Francis & Kannenberg, 1978; Kharel et al., 2019). While looking at individual field maps, for example, most of W-12 had greater than mean yield with high CV indicating that the portion of the field produced higher yield but were unstable (Zone B). The middle portion of field 16-A was stable with greater than mean yield, and was thus classified as Zone A. Similarly, a major portion of Field Y-8 in the east, Y-6 in the north, and W-13 in the south were classified as Zone D, that was stable in producing low yields during the study period.

Fig. 4
figure 4

Yield stability zones derived from corn yield data for 2018, 2019, and 2020 [Zone A: High yield and stable; Zone B: High yield and unstable; Zone C: Low yield and unstable; and Zone D: Low yield and stable]

Full size image

Table 3 shows the extent of each stability zone for the study fields. Over 42% of total study area was covered by Zone B, 29% by Zone D, while Zone A and Zone C represented about 14% each. The areal coverage and extent of these unstable zones can be strategic to apply precision conservations decision (Berry et al., 2005; Gelder et al., 2008). While looking at individual fields, most fields had greater extent of Zone B compared to the other zones, except for Y-10, and W-13 which had a maximum area, 29.9%, and 42.9%, respectively, under Zone D.

Table 3 Area covered in hectare (ha) by yield stability zones in different fields
Full size table

Gross margin map

Figure 5 shows gross margin maps for all fields based on yield and farm economics data from 2018, 2019, and 2020. The gross margin values ranged between − $693 and $775/ha across the study fields. Four fields (Y-6, W-13, W-12, and 6–12) resulted in positive gross margins, whereas the remaining fields (16-A, SW-16, Y-10, Y-13, and Y-8) returned negative gross margins. Field W-12 had the highest gross margin average ($541/ha) followed by 6–12 with a margin of $161/ha, and Y-10 had the lowest margin (− $257/ha) of all. Overall, the distribution of gross margin across fields followed the spatial pattern of predicted yield as the map was heavily dependent on yield (Blackmore, 2000). The northeastern part of Y-8, southern part of W-13, and northwest part of Y-6 had low or negative gross margin relative to the rest of the respective field. Similarly, most of W-12, 6–12 and Y-13, and eastern part of SW-16 had a relatively greater gross margins, whereas gross margin from Y-10 was a mixture of high and low values in regular patterns based on terrace placement. Specific to field 16-A which had a gross margin of − 166/ha, a 35–40-m wide channel running north–south towards the northwest corner of the field had a large negative margin compared to the rest of the field.

Fig. 5
figure 5

Average gross margin across fields based on gross margin maps of 2018, 2019, and 2020

Full size image

Evaluating gross margin by yield stability zones, we observed that all zones from fields W-12, W-13, and 6–12 had a positive margin, whereas all zones in 16-A, and Y-10 had negative margin (Fig. 6). A maximum positive margin among all fields was reported for Zone B from W-12, and the minimum (most negative margin) for Zone D of Y-8. Zones A and B had positive gross margins for all fields except for 16-A, and Y-10 where gross margins were negative, whereas Zones C and D had negative gross margins in all fields except for W-12, 6–12, and W-13. It has been reported that the gross margin maps can be an important tool in PA, as they assess the spatial variability of economic return and assist farmers to manage their farms more efficiently (Bazzi et al., 2015). Such maps can be used as a basis for precision conservation decisions with economic justification (Blackmore, 2000), as on-farm adoption of conservation practices depends on yield and economic impact, among other factors (Perry-Hill & Prokopy, 2014).

Fig. 6
figure 6

Average gross margin by yield stability zones. Error bars represent standard deviation

Full size image

Precision conservation plan and cost saving

A precision conservation plan option for the study area was to leave Zone D from selected fields out of production leaving it be taken over by the native plant species. This would have multiple benefits such as promoting native plant species and pollinators, could act as a natural grass water ways trapping soil and nutrient loads from upstream, and build-up soil organic matter naturally, among others. Implementation of precision conservation plans using a management unit approach has also been suggested by Gelder et al. (2008); however, the approach was slightly different than the one we adopted in this study. Table 4 lists average input cost by stability zones during 2018, 2019, and 2020. Average input cost in Zone D from different fields ranged between $929 and $2,201, and by adopting the precision conservation plan in each field, a total of $13,877 could have been saved in the past 3 years.

Table 4 Average input cost in US dollar ($) from yield stability zones in 2018, 2019, and 2020
Full size table

As a part of long-term research, information from this study was used to adopt precision conservation in some of these fields. Fields SW-16 and W-13 were selected to reduce inputs (seed and fertilizer) for the contiguous areas of Zone D to 60% of the full rates used for corn crops during 2018–2020. This would reduce costs by roughly $450 and $570 for fields SW-16 and W-13, respectively. Further, Y-8 and Y-13 were selected to remove Zone D from row crop production to be replaced with permanent cover, reducing annual input costs in these two fields by more than $3,200 combined. While field Y-6 does have a large contiguous area that would be a prime target for such precision conservation, we chose to not alter management of this field to serve as a control to the other fields.

Conclusions

Three years (2018–2020) of yield data from nine corn fields in Texas Blackland Prairie soils were used to map within-field yield variability, as well as gross margin variations using geostatistical techniques. Each field was divided into yield stability zones based on 3-year average yield and associated CV. Farm economic data from each field was used to quantify gross margin, and the values were mapped across the fields for spatial assessments. Yield stability zones were then used as baseline management units for which gross margin values were calculated per zone basis, and compared within, and across fields for potential precision conservation decisions.

There was large variation in corn yield measurements during 2018, 2019, and 2020 in the study area. Further, yields in 2018 were poor due to a very dry growing season. Mean yield during the study period was greatest for field W-12 (8.17 Mg/ha) and lowest for field Y-10.

This study demonstrated that yield stability zones could be efficiently derived from the mean yield and CV using several years of yield data. A major portion (57%) of the study area was classified as an unstable zone with high CV and either high or low yield. Zone D (low yield with low CV) represented nearly 29% of the area. Gross margins ranged between − $693 and $775/ha. Overall, almost all fields under zones A and B had positive gross margins, and zones C and D had negative gross margins.

For the next step of this research, inputs will be reduced in Zone D in fields SW-16 and W-13 and Zone D will be removed from production for fields Y-8 and Y-13. Based on previous years this will result in reduced input costs for these fields of $450 to $2 200. We will continue to monitor yields and gross margins from these fields and will also monitor environmental impacts of these precision conservation management options.

Overall, this study highlighted a methodology to identify yield stability zones and to map gross margins across corn fields in Texas Blackland soils and proposed a suitable precision conservation plan with economic justification. This approach could be an effective conservation tool helping farmers save money and contribute to environmental benefits. We believe that similar methodology could be applied in other regions as well for yield stability assessment as conservation benefits including economic return.