Introduction

Rainfed agriculture plays an important role in food production worldwide as it accounts for more than 75% of global cropped area, being responsible for more than 60% of global cereal production (Cassman et al., 2003; Connor & Mínguez, 2012). Sustaining food production by rainfed crops in the years ahead will require productivity gains in resource use (Fischer & Connor, 2018). Specific challenges include estimating the magnitude and thus the value of yield gaps, identifying limiting factors, and implementing profitable and sustainable strategies. Closing yield gaps in rainfed farming, while improving resource use efficiency, is also expected to minimize the expansion of arable land and emissions while reducing other undesirable externalities (Cassman et al., 2003; Snyder et al., 2009). Over the last few decades, new technologies have been evolving within the context of rainfed farming (e.g., precision agriculture), aiming to increase the productivity of resources while offering substantial environmental benefits (Griffin et al., 2018).

Spatial variations in rainfed crop yields are commonly observed worldwide (Florin et al., 2009; Griffin et al., 2020), which reveal, theoretical, opportunities for optimizing resource use (e.g., fertilizer) through precision agriculture (PA). A central concept of PA is the spatially variable rate application (VRA) of fertilizer according to the intra-plot variations of yield (Basso et al., 2013; Bullock & Lowenberg‐DeBoer, 2007; Robertson et al., 2008; Pedersen et al., 2021). In fields of undulating topography, where rainfed crops experience different degrees of water stress due to spatial water variations caused by lateral flows, yield varies spatially within the same field (Halvorson & Doll, 1991; Tenreiro et al., 2022). These variations may imply different nutrient requirements and application rates over fields (Nielsen & Halvorson, 1991; Sadras, 2002). However, although the theoretical reasons for adopting VRA technology in rainfed systems are well accepted, its adoption has not been widespread (Basso et al., 2013; Lowenberg-DeBoer & Erickson, 2019; Robertson et al., 2012).

Considerable attention has been devoted to VRA adoption constraints (Robertson et al., 2008; Schimmelpfennig, 2016; Welsh et al., 2003). Sceptics argue that there is a lack of evidence of the relative advantage of VRA adoption. The principal issues are related to technical obstacles, complexity on the use of equipment and software, problems with data availability and access, operational incompatibility and ambiguity regarding financial benefits (Pathak et al., 2019; Robertson et al., 2007). Nevertheless, current trends in energy prices, transportation costs and impacts on both cereal prices and fertilizers costs (Chowdhury et al., 2021; Deloitte, 2021; EUC, 2021; FAO, 2021; Glauber & Laborde, 2022; Khalfaoui et al., 2021; USDA, 2021) may alter the trade-offs involved in the adoption of VRA in rainfed systems. Generally, there are other advantages attributed to the increase of fertilizer use efficiency that go beyond economic reasons which should not be ignored. Adequate use of fertilizer is not only needed to sustain yields and to increase the efficiency of water and energy use, but it also contributes to minimize the overall emissions throughout the supply chain (Snyder et al., 2009).

Understanding and capitalizing on yield variability is one of the main objectives of PA and the yield gap (YG) is an important concept to be used in this context as it provides a benchmark to explore yield variations (Cassman, 1999). The YG is defined as the difference between the potential (Yp), or water-limited yield (Yw) in the case of rainfed cropping systems, and the actual yield (Ya) achieved by farmers (Fischer, 2015; Loomis & Connor, 1992). However, most studies on YG analysis ignore intra-plot variability (Fischer et al., 2014; Lobell et al., 2009; Schils et al., 2018), which is partially due to data availability constraints (Beza et al., 2017) and to the limitations of crop models to simulate processes such as spatial water distribution (Tenreiro et al., 2020b).

In fields of undulating topography where water flows from higher to lower elevation zones, the assessment of spatial variations of Yw may be combined with yield mapping from combine harvesting to understand how YG varies within a field. Assuming that Ya is achieved under current management practices (i.e.., uniform nutrient application), the magnitude of the YG and its variation over space (i.e., how site-specific Ya relates to Yw and how it varies from upslope to downslope zones) could be taken as an indication of the relative advantage of VRA adoption. The larger the differences in YG among different zones, the greater would be the opportunity for adopting VRA (Robertson et al., 2008).

This study builds upon a methodological framework that assessed both Yw and Ya spatial variations in rainfed wheat fields (Tenreiro et al., 2022), and it explores how the YG varied within fields and from year to year. Intra-plot yield variation was modelled, and management units were determined according to the YG variation. Some economic implications associated with multiple scenarios of VRA adoption were explored through a case study and recommendations were proposed. The present study performs a novel economic analysis on VRA adoption to benefit from the spatial water variations in rainfed wheat systems. The hypothesis that spatial variations in soil water supply may justify a variable fertilization rate over space is tested, considering the following objectives:

  1. 1.

    The identification of economic trade-offs between the scale of cropping and opportunities for VRA adoption;

  2. 2.

    The assessment of return on investments at the farm scale in different scenarios, where the spatial water variations over fields determine yield variability and scope for adopting VRA.

Materials and methods

Experimental conditions and on-farm data collection

The experimental farm is located in Córdoba province, southern Spain (37.8° N, 4.8° W, mean elevation 165 m amsl., mean precipitation 605 mm year−1, range of slope steepness 2–6%), and it has a total arable area of 320 ha (Fig. 1a). At the farm level, crop rotations include autumn and spring sown crops (e.g., winter wheat, rapeseed, sunflower, chickpea), with winter wheat being grown 1–2 times every four seasons. This study is limited to wheat (Triticum durum L.), which is a major crop within the farming system with a relative area share of 0.25–0.3 year−1. The soils are predominantly of clay texture (40–50% clay and 15–22% sand) with high bulk densities (1.65–1.88 g cm3) and moderate depth (1.2–1.6 m). Soil organic matter content is on average 1.05% (± 0.3) and mean pHKCl is 7.4 (± 0.05). More details on soil properties and other geophysical conditions are provided in Tenreiro et al. (2022).

Fig. 1
figure 1

Experimental fields: a Total farm area (320 ha); b Digital Elevation Model (DEM), values expressed in m amsl.; c Soil types according to USDA classification system; d TOPMODEL Topographic Index (TMTI), values are unitless; e Yw zones map (i.e., ‘LIF’ and ‘no-LIF’ zones)

Full size image

Depending on the year conditions and operational aspects at the farm level, sowing dates varied from mid-November to end of December (Table 1). Wheat is sown under direct seeding at an average rate of 200 (± 20) kg ha−1. The current N management system is spatially uniform, following a fertilization plan that consists of two applications per season (one pre-stem elongation and the other near the end of the vegetative stage). Mean N-P-K application rates to wheat are 165-60-0 units per ha year−1. Typical N-fertilizers are DAP, CAN, and Urea, namely diammonium phosphate (18-46-0), calcium ammonium nitrate (24-0-0 + 8% Ca) and urea fertilizer (46% N), respectively. The crop potential yield is used to define the crop N net requirements, which are considered as a proxy for defining the N application rate. Total N rate is estimated according to a simple N-balance, assuming that N inputs are fully recovered in the harvested grain. Total N rate is equally divided among the two applications conducted. More information on crop management is provided in Table 1.

Table 1 Wheat crop management characteristics in different fields on the experimental farm, located in Córdoba, Spain, over six cropping seasons (2015–2021)
Full size table

Lateral inflow (LIF) zone mapping

Spatial water variations caused by lateral flows imply yield variations. Zones with different potential yields, within the same field, were delineated according to lateral inflow (LIF) levels. For zone delineation, the TOPMODEL Topographic Index (TMTI) was used as described in detail by Beven et al. (2021) and expressed as following:

$${\text{TMTI}}_{{({\text{x}},{\text{y}})}} = {\text{ ln}}\left( {\frac{{\text{a}}}{{{\text{tan}}\beta }}} \right)$$
(1)

where a represents the flow accumulation, defined as the upslope contributing area to that point (x, y), and β represents the field slope angle in the same point. While the hydraulic gradient ‘tan β’ is defined with respect to the plan distance, the flow accumulation rate is defined with respect to the plan unit area. a is expressed in m2 and β in degrees.

This investigation scaled up to the farm level the results of Tenreiro et al. (2022), who followed an experimental design that consisted of choosing three different sampling zones per field, where the TMTI values (expressed in m2 per unit contour length) ranged from 0 to 6 with median values varying from 1 to 2 in ‘no-LIF zones’ up to 4–6 in ‘LIF zones’ (Appendix Fig. 6). The Yw simulation results in Tenreiro et al. (2022) are assumed to be dependent on the scale of the experimental fields. To avoid an extrapolation of results to different scales from those in experimental conditions, all points with a TMTI falling outside of the distribution range in experimental conditions (Appendix Fig. 6) were rejected (i.e., points with TMTI > 6 were discarded). For the same reason, only those fields with a predominance of clay texture were selected (Fig. 1). According to Tenreiro et al. (2022), only two zones were significantly different in terms of soil water content variations due to LIF (i.e., zones A and C, Appendix Fig. 6). Therefore, in the present study, only two different zones were considered:

  • ‘LIF zones’ (i.e., downslope zones with significant amount of water supplied through lateral flow taking place from upslope areas of the same field);

  • ‘No-LIF zones’ (i.e., upslope zones where an insignificant amount of water is supplied through lateral flow);

According to the standard deviation of each significantly different zone (Appendix Fig. 6), a TMTI cutting threshold value equal to 5 was considered and each location was classified as following:

  • IF TMTI(x,y) < 5 [point classified as ‘no-LIF zone’];

  • ELSE [point classified as ‘LIF zone’];

Based on the previously described criteria, a subset plot of 92 ha was selected (Fig. 1d, e), from which 76 ha were classified as ‘no-LIF zones’ while the remaining 16 ha were considered to be ‘LIF zones’. Therefore, the selected farm area has an overall share of ‘LIF zones’ equal to 17.4% (Fig. 1e).

According to Tenreiro et al. (2022), ‘LIF zones’ are characterized by significant water supplied through lateral inflow, while ‘no-LIF zones’ are characterized by insignificant lateral inflow. TMTI is described in detail by Beven et al. (2021). TMTI values ranged from zero to five in no-LIF zones and were larger than five in LIF zones (Appendix—Fig. 6).

Intra-plot spatial assessment of yield gaps

The YG was spatially determined at the intra-plot level with a spatial resolution of 10 × 10 m using the following equation:

$${\text{YG}}_{{({\text{x}},{\text{y}})}} = {\text{ Yw}}_{{({\text{x}},{\text{y}})}} {-}{\text{ Ya}}_{{({\text{x}},{\text{y}})}}$$
(2)

where YG, Yw and Ya represent the yield gap, the water-limited yield and the actual yield, respectively. The subscripts (x,y) indicate the corresponding polygon centroid coordinates. Yields are expressed in terms of Mg DM (grain) ha−1.

Two zones were considered according to lateral inflow (LIF) water supply (i.e., LIF and no-LIF zones) and Yw was year-specific, and was assumed to be constant within each zone. Yw values (2016–2021 period) were obtained from Tenreiro et al. (2022) and Ya was determined with the ‘New Holland’ Precision Land Manager (PLM) software, with a spatial resolution of 10 × 10 m, taking as an input the shapefiles generated by the combine harvester monitor (Fendt PLI C 5275).

To obtain LIF zones’ Yw, lateral inflow was simulated as an additional water supply to the effective rainfall using the AquaCrop simulation model (Steduto et al., 2009) as described in Tenreiro et al. (2022). For No-LIF zones’ Yw was simulated, assuming that water inflow takes place only due to vertical infiltration (Tenreiro et al., 2020b). Figure 2 shows seasonal water supply over the six seasons considered. Seasonal precipitation (P) from sowing to harvest varied from 228 to 535 mm and seasonal LIF, which was only considered in ‘LIF zones’, ranged from 30 to 112 mm. Full rainfall information regarding the cropping system analysed in this manuscript is provided in Tenreiro et al. (2022).

Fig. 2
figure 2

Water supply: seasonal precipitation (P) and lateral inflow (LIF) for the experimental dataset. Values are expressed in mm. More information is provided in Tenreiro et al. (2022)

Full size image

Quantifying the relative advantage of VRA adoption

The assessment focused on the relation between the intra- and inter-season economic benefits of VRA adoption, and on the total investment costs in relation to the amortization of equipment. The relative advantage (Robertson et al., 2012) in terms of income risk reduction or environmental benefits, as done by Swinton and Ahmad (1996), was not explored.

Exploring the relative advantage of VRA adoption must dedicate attention to the scope for improvement, not only regarding the target scenario but also the current management system (i.e., the baseline scenario). The baseline scenario is characterized by a uniform nutrient application rate and thus its agronomic performance is captured by the determination of Ya. The target scenario implies VRA adoption, and it is benchmarked by the simulated Yw level. Therefore, the relative advantage is estimated as the difference between the target and the baseline scenario. This was expressed as the differential gross margin (DGM), which is the difference between the gross margin obtained with VRA and without it (Pedersen et al., 2021). Since the gross margin equals revenue minus expenditure, it was computed as a function of grain yield. Grain yield determines both the output revenue and the input costs (associated with fertilizer requirements). The yield difference between the target (Yw) and the baseline scenario (Ya) is defined as the YG.

The relative advantage was measured as the average DGM. The average DGM (€ ha−1) was estimated by averaging the annual ADGM values (€ ha−1 year−1), expressed as:

$$DGM = \mathop \sum \limits_{t = 1}^{T} \frac{{ADGM_{t} }}{T}$$
(3)

where T is the number of seasons (N = 6 years) of the dataset (i.e., 2016–2021), and the subscript t indicates the specific season. ADGM is the difference between the annual differential revenue (ADRt) and the annual differential costs (ADCt), expressed as:

$$ADGM_{t} = ADR_{t} - ADC_{t}$$
(4)

where both ADRt and ADCt are expressed in € ha−1 year−1. Since the DGM represents the differential gross margin, it applies to the differential area, which is defined by the LIF area share. Therefore, both ADRt and ADCt are defined per ha of LIF zones. The ADRt was computed as a function of YG, expressed as follows:

$$ADR_{t} = YGz_{t} \cdot 1000 \cdot \left( {Price_{WHEAT} + LPP} \right) + DPA$$
(5)

where YGzt represents the YG difference among zones (i.e., between LIF and no-LIF zones) for a specific year (t), and the PriceWHEAT corresponds to the grain price (€ kg−1 DM grain). LPP represents a linked production payment in some scenarios (i.e., a subsidy linked to yield and expressed in € kg−1 DM grain), and DPA represents a direct payment on the annual area sown with wheat (€ year−1). DPA does not consider current direct payments, but it explores the introduction of an additional payment on area over current subsidy levels in one scenario. The larger YGzt is, the higher is the expected ADRt. The YGzt was computed as follows:

$$YGz_{t}= \sqrt {\left(\left[\mathop \smallint \limits_{cum.d = 0}^{cum.d = 1} ECDF\left( {YG_{LIF} } \right)\right]-\left[\mathop \smallint \limits_{cum.d = 0}^{cum.d = 1} ECDF\left( {YG_{No - LIF} } \right)\right]\right)^2}$$
(6)

where ECDF(YGLIF) is the ‘empirical cumulative distribution function’ of YG in ‘LIF zones’ and ECDF(YGNo-LIF) is the ‘empirical cumulative distribution function’ of YG in ‘No-LIF zones’. Cum.d delineates the cumulative density, varying from 0 to 1. The 'empirical cumulative distribution function' is the distribution function of a sample measure. In this specific case, this applies to the yield gap distribution over field. It expresses the fraction of yield gap observations which are lower than, or equal to, a specified value. For ECDF computation, the function ‘stat_ecdf’ from the ‘ggplot2’ library in R-studio was used (Wickham, 2007).

The ADCt was estimated by the following equation:

$$ADC_{t} = (YGz_{t} \cdot 1000 \cdot N_{GRAIN} \cdot Price_{N} ) + AC_{VRA}$$
(7)

where NGRAIN is the nitrogen content in grain (expressed in terms of % DM), the PriceN is the price of N (expressed in € per kg N) and the ACVRA is the annual cost (Table 2) associated with VRA use (e.g., technology renting plus external consultancy costs). NGRAIN was set equal to 2.8% of DM (Quemada et al., 2016). PriceN was estimated according to the 2020/2021 fertilizer price index for Diammonium Phosphate (DAP), Calcium Ammonium Nitrate (CAN) and Urea (EUC, 2021), and considering the mean N content of fertilizer which was estimated as the simple average of N content in DAP, Urea and CAN. The baseline PriceN was set equal to 1.10 € kg−1 N.

Table 2 The ten different scenarios considered to analyse the impact of price support policies, extra direct payments, and different market prices on VRA economic relevance
Full size table

Economic modelling and analysis of future scenarios

The capital recovery of VRA adoption was modelled according to 10 future scenarios (Table 2) which considered different prices and payments (i.e., application of extra subsidies). The Net Present Value (NPV) was estimated over a period of 10 years which was assumed to be the total amortization period of VRA equipment (Drabik & Peerlings, 2016; Tozer, 2009). NPV was estimated as follows:

$$NPV = - II_{VAR} + \mathop \sum \limits_{i = 1}^{N} \frac{{DGM_{i} \cdot LIF_{AreaSHARE} }}{{\left( {1 + \Upsilon } \right)^{i} }}$$
(8)

where IIVRA is the initial investment (i.e., the acquisition cost), ϒ is the discount rate, i is the cash-flow year and N represents the total amortization period (10 years). The acquisition costs (Appendix—Table 6) were based on information compiled from different sources (AAEA, 2000; Batte & Ehsani, 2006; Finco et al., 2021; Griffin, 2006; Tozer, 2009), including the acquisition costs of the GPS guiding system, the precision application system RTK, the GPS receiver, the base station, the replicators and the N application controller (Appendix—Table 6). The overall total gain (OTG) was estimated by solving the NPV series over a period of 10 years (N = 10), which was a function of both the discount rate used and the price/payment scenario adopted (Table 2). The Internal Rate of Return (IRR) and the return-on-investment payback time (ROIt) were respectively estimated by solving the following equations:

$$II_{VAR} = \mathop \sum \limits_{i = 1}^{N} \frac{{DGM_{i} \cdot LIF_{AreaSHARE} }}{{\left( {1 + IRR} \right)^{i} }}$$
(9)
$$ROI_{t} = \frac{{II_{VAR} }}{{\left( {DGM \cdot LIF_{AreaSHARE} \cdot Wheat_{area} } \right)}}$$
(10)

where Wheatarea is the annual area sown as wheat (ha). The relation between ROIt and Wheatarea was obtained through regression analysis, for each price/payment scenario.

The impact of price support policies, extra direct payments on crop area and different market prices were analysed according to ten different scenarios (Table 2). The considered ratios for increased price scenarios were defined according to the observed connectedness between fertilizer and grain market prices. By averaging the values reported by Chowdhury et al. (2021) and Khalfaoui et al. (2021), a ratio among product prices of 0.46 was assumed, meaning that a 100% increase in fertilizer prices would increase wheat prices by 46%. The comparison of capital recovery among different scenarios was made by adopting the current situation as the baseline scenario. For each scenario, three different sub-scenarios were considered by setting the discount rate at 2.5%, 5% and 7.5% (Table 2).

The following 10 scenarios were considered:

  1. 1.

    Baseline scenario (S-1) overall ‘LIF-zones’ share within the cropping system equal to 17.4%, according to the farm conditions, current CAP direct payments and product prices.

  2. 2.

    Enhanced LIF area share scenario (S-2) overall ‘LIF-zones’ share within the cropping system increased by 35% in relation to the baseline scenario, current CAP direct payments and product prices.

  3. 3.

    Moderate increased prices scenario (S-3) baseline scenario plus product prices increased by 30% and 66%, respectively for wheat grain and N fertilizer price.

  4. 4.

    Drastically increased prices scenario (S-4) baseline scenario plus product prices increased by 100% and 219%, respectively for wheat grain and N fertilizer price.

  5. 5.

    Price support scenario (S-5) baseline scenario plus the introduction of a direct payment linked to production (DPP) equal to 0.02 € kg−1.

  6. 6.

    Price support plus drastically increased prices scenario (S-6) scenario S-5 plus product prices increased by 100–219% (i.e., wheat grain and N fertilizer price).

  7. 7.

    Additional direct payment on cropped area scenario (S-7) baseline scenario plus the introduction of an additional payment on cropped area equal to 46 € ha−1.

  8. 8.

    Additional direct payment on crop area plus drastically increased prices scenario (S-8) baseline scenario plus the introduction of an additional payment on crop area equal to 46 € ha−1 plus product prices increased by 100–219% (i.e., wheat grain and N fertilizer price).

  9. 9.

    Support VRA investment (S-9) baseline scenario plus the introduction of a support on investment corresponding to 50% of initial investment covered by a subsidy.

  10. 10.

    Support VRA investment plus drastically increased prices scenario (S-10) baseline scenario plus the introduction of a 50% subsidy of initial investment covered by a subsidy, plus product prices increased by 100–219% (wheat grain and N fertilizer prices).

Currently, there is no specific support for cereals in Spain, where the present study was conducted, but the example of what was allocated as crop-specific aid per ha in the cases of oilseeds and legume-crops during the year 2020 was followed (FEGA, 2021). According to existing crop-specific aids, a mean support of 40–55 € ha−1 was assumed. This interval supposes that for the mean yield levels here considered, an additional payment of 0.02 € kg−1 could be applied as LPP, equivalent to a support of 46 € ha−1 as DPA.

The relation between return time on investment (ROIt, expressed in years) and the annual wheat area (expressed in ha) was modelled. Regression analysis was used to estimate model coefficients, several alternative models were tested, including linear, quadratic, logarithmic, exponential and power models. Least Squares Fitting (LSF) and statistical hypothesis testing were respectively used to estimate the regression coefficients and their significance level (stats package in R; Team, 2000). Root mean square error (RMSE) and R2 were used as statistical indicators of performance for the best-model selection. For each scenario, the minimum area for adoption, corresponding to the threshold below which the return on investment takes longer than the amortization of VRA equipment, was also estimated. The null hypothesis was tested with the non-parametric Tukey's range test (HSD-test).

Results

YG analysis

Figure 3 and Table 1 show simulated Yw values (2016–2021), which were obtained from Tenreiro et al. (2022). Measured Ya values ranged from 0.3 to 6.4 Mg DM ha−1 over space and time. The intervals of intra-plot variation for each season are shown in Fig. 3. Ya showed a mean coefficient of variation of 14.7% (± 3.1%) among the six seasons considered.

Fig. 3
figure 3

Yield variations: a Cumulative probability distribution of actual yields (Ya) for each season (2015/2016 to 2020/2021) and within each zone (i.e., solid lines correspond to ‘no-LIF zones’ while dashed lines relate to ‘LIF-zones’); b Cumulative probability distribution of yield gaps (YG) for each season and within each zone. Vertical lines indicate Yw levels obtained from Tenreiro et al. (2022) simulations

Full size image

The probability distribution curves of YG between zones and within each season are also shown in Fig. 3. YG values ranged from 0 to 4.8 Mg DM ha−1. Median YG was 1.7 Mg DM ha−1 with a standard deviation of 0.135 Mg DM ha−1. In relative terms, the average YG was 30.5% of Yw. The YG’s were larger in LIF zones in four out of six seasons (Fig. 3). The yield maps are shown in the Appendix (Figs. 7, 8, 9), and the YG results are synthesized in Table 3.

Fig. 4
figure 4

The negative relation between return time on investment and the annual wheat area. The annual area sown as wheat is subjected to an average LIF area share of 17.4%. Scenarios (a–e) are described in Table 2

Full size image
Table 3 Mean yield values (Ya), yield gaps (YG), and coefficients of variation (CV) in LIF and No-LIF zones
Full size table

The economics of VRA adoption

The capital recovery of VRA adoption was expressed as a function of both the scale of adoption (i.e., annual area sown with wheat) and the financial discount rate. This relation was directly affected by the differential gross margin (DGM), which was computed as the average of Table 4 values.

Table 4 The economics of VRA adoption, under spatial water variations in rainfed wheat systems according to the experimental conditions
Full size table

On average, DGM ranged from 12.1 to 147.5 € ha−1, depending on the scenario considered. Under current conditions (see S-1 and S-2 in Table 4), three out of six seasons (i.e., 2017/2018, 2018/2019, 2019/2020) presented a negative ADGMt, indicating risk of economic inefficiencies associated with VRA adoption for 50% of the years investigated. DGM was negative when ADRt was lower than ADCt. For each year’s specific conditions, there was a great variation of ADGMt, which ranged from -42.3 (baseline scenario) to 405.2 € ha−1 (scenario S-8, i.e., a policy support through an additional direct payment of + 46€ ha−1 on crop area plus high wheat grain and N fertilizer prices; Table 2).

The economics of VRA adoption varied considerably, not only among different years but also for the different scenarios (Table 4). The ADGMt and ADRt showed stronger sensitiveness to the differences between years and economic scenarios than ADCt, which presented a lower range of variation (Table 4). LIF coefficients (i.e., season LIF divided by season precipitation) varied from 11 to 27% according to Fig. 2 values. Those years with higher LIF coefficients (e.g., 2016/2017 and 2020/2021) showed higher ADRt and ADGMt. The opportunity cost showed a trend to decrease with the increase of LIF coefficient, but regressions were not capable of explaining more than 30% of the overall variation (results not shown).

Economic trade-offs for viability of VRA

Model coefficients were obtained through regression analysis (Table 5). From the set of models tested, power models were selected as the best fitting. Power models maximized R2 values and minimized the RMSE. While the R2 values of linear, quadratic, logarithmic and exponential models were respectively < 0.1, 0.1–0.2, 0.2–0.4 and 0.7–0.8, the RMSEs (expressed in years) were larger than 10 for the first three model types and approximated 0.5 for the exponential models. Power models (a*x−1) were the best fitting models with R2 > 0.98 and RMSE < 0.1. Model coefficients varied from scenario to scenario but the statistics of model performance did not change between scenarios. The model coefficient a and the minimum area for adoption of VRA are shown in Table 5 for each of the 10 different scenarios. The minimum area for adoption corresponds to the threshold value below which the return on investment takes longer than the amortization of VRA equipment.

Table 5 Parameters of the best fitting (power) models, the minimum area for adoption of VRA, the mean return on investment time and the mean overall total gain (OTG)
Full size table

There was a negative relation between the return time on investment (ROIt, expressed in years) and the annual area sown as wheat(ha) which was best fitted by a power model. Figure 4 shows the power models that were fitted for each economic scenario. The 10-year overall total gain (OTG) changed significantly from case to case, and according to the different discount rates (Fig. 5). The lower the discount rate, the larger the profitability of VRA investment.

Fig. 5
figure 5

The 10-year overall total gain (OTG) under multiple scenarios (units expressed in thousand €). The OTG was estimated by solving the Net Present Value series over a period of 10 years, which was a function of both the discount rate used and the price/payment scenarios adopted (a–e as described in Table 2)

Full size image

Over a 10-year period and under the baseline scenario, at a discount rate of 2.5%, the OTG is expected to vary from  − 100.9 thousand € (for a total of 50 ha of wheat sown every year) to 310.3 thousand € (in the case of 2000 ha of wheat sown annually). Under the most profitable scenario (S-8), for a median farmer with 100 ha wheat sown per year, the 10-year OTG is forecasted as 125.3, 134.2 and 144.6 thousand €, respectively with a discount rate of 7.5, 5 and 2.5% (Fig. 5 and Supplementary material).

In this case (i.e., annual sown area equal to 92 ha with an overall share of ‘LIF zones’ equal to 17.4%), the internal rate of return (IRR) was negative for almost all cases considered (Table 5 and Fig. 5), indicating a lack of relative advantage, or additional profitability, for all the scenarios and discount rates considered. Over the ten different scenarios explored, and assuming a median case of 92 ha sown, it was observed that the capital return could increase up to 5% in the case of scenario S-8 (Table 5 and Appendix—Table 7), but this would depend on a drastic price increase, or other changes in agricultural policies.

Discussion

Crop yields and yield gaps

Actual wheat yields (Ya) varied within the range of other studies (Padilla et al., 2012; Schils et al., 2018) and exhibited significant differences among zones in two out of six seasons (Table 3). This is associated with yield maps fitting into variograms with irregular spatial structure as both CV and differences among means did not show consistency from year to year. This highlights that temporal instability is an important issue for site-specific management because the agronomic implications of asymmetric spatial variations differ greatly with the crop × year setting, as also discussed by Tenreiro et al. (2020a).

Coefficients of variation were consistently higher in ‘no-LIF zones’, which is in line with the results of Tenreiro et al. (2022), who highlighted that the contribution of LIF to yield spatial variations tends to be stronger in years of relatively low water supply. This indicates that the degree of variation tends to increase with the level of water stress. Observed CVs of Ya were similar to those reported by others (Batchelor et al., 2002; Florin et al., 2009; Whelan & McBratney, 2000).

The CVs of YG were notably larger than CVs of Ya, which is attributed to more sources of variation affecting the process of YG mapping. According to Table 3, YG’s were statistically, and systematically, different between zones. This indicates that the YG appears to be a more precise benchmark (instead of using solely Ya or Yw) for precision management decision-making, because it captures a higher degree of spatial variation and it delivers the magnitude of the expected response under a VRA for each year-specific conditions.

The economics of VRA adoption

The results presented here are conditioned by the specific topographic conditions of the chosen farm. Under different geomorphological conditions (e.g., increased LIF area share), the returns on investment are expected to change. The larger the relative share of LIF zones within the farming system, the lower the minimum area required for VRA adoption (see S-2 in Table 5).

Under current conditions (S1), a relative advantage associated with VRA adoption was computed but only for an annual area sown as wheat larger than 567 ha year−1 (Table 5). This is considerably larger than representative European (arable) farm sizes, which typically range from 4 to 62 ha (Andersen, 2017). For cases with less sown area, the investment costs could still be recovered but it would not be due to the relative additional gain. This means that a farm with a lower annual sown area and currently profitable with uniform N applications, could pay for the VRA investment, but the overall profitability of the farm would decrease if VRA is adopted.

All the alternative scenarios considered accelerated capital recovery (Table 4 and Fig. 4). However, wheat and N prices were the most important factors for VRA viability (Table 5 and Fig. 4). The slope of the 10-year OTG increased notably with both wheat and N prices increasing (Fig. 5). However, this is determined by the N-grain prices’ relation adopted. The most promising scenario, which showed the largest gain on capital recovery (see S-8 in Table 4 and Fig. 5), was the inclusion of an additional direct payment on crop area (+ 46 € ha−1) plus a drastic evolution of both input and output prices.

An additional direct payment on area could also deliver advantages for a wider range of farmers as the minimum annual wheat area for adoption decreases more than 75% under scenarios S-7 and S-8 (Table 5). For the same available budget considered, the introduction of an additional direct payment on area would imply a more robust advantage for VRA adoption than payments linked to production (Tables 4, 5). Payments linked to grain yield would be mostly diluted over the Ya level, which relates to the current management system, and not over the differential YG closing effort that is attributed to VRA adoption. Since current mean YGs are approximately 30% of mean Yws (Table 3), the major part of subsidy support linked to production would not apply to the differential gross revenue, caused by the technological shift, but rather to the actual yield level that is already achieved under the current management (i.e., uniform nutrient application). In this sense, a much wider fraction of financial inflow would be directly attributed to the technological shift in the case of an additional payment on area (Table 4). This indicates that the direct payment of 46 € ha−1 year scenario could be a better option in terms of policy support on investment.

The opportunity of lower initial investment costs associated with lower VRA acquisition costs of equipment was not investigated, but this was partly explored through scenarios S-9 and S-10. These scenarios explored the impacts of a policy support on investment through a support-payment equal to 50% of the initial acquisition costs (Table 2). This strategy impacts the return on investment by decreasing ROIt (Fig. 4e), but achieving significant impacts on OTG requires further price changes (Fig. 5e). In the absence of further price changes, this option does not guarantee incentives for adoption by farmers and this could be a constraint from a strategic point of view.

The results showed that the profitability of VRA investment would respond more to changes in market prices than to policy supports (Tables 4, 5). In the absence of additional policy supports, the minimum area for adoption of VRA decreased substantially under both price evolution scenarios (i.e., respectively 69% and 88% according to Table 5). The evolution of prices that were assumed is likely to turn VRA into a viable technology for a much wider population of farmers, as the minimum area for adoption, according to the considered range of price increases, is expected to decrease from 567 ha to an interval of 67 to 177 ha year−1. This clearly indicates that, given current trends on both wheat and N prices observed in early 2022 (Glauber & Laborde, 2022; Vos et al., 2022), more farmers would be inclined to adopt VRA in local rainfed systems.

The mean OTG values shown in Table 5 were estimated as the simple average of OTG obtained for the following series of sown areas: 50, 100, 250, 500, 1000 and 2000 ha. According to the HSD-test results shown in Table 5, scenarios S-1, S-2, S-5 and S-9 did not significantly differ from each other, which indicates that different market prices are the most important condition to improve VRA economic viability.

The obtained ADGMt fluctuated in a wider range than the values reported by Robertson et al., (2007, 2012). The results ranged from − 42.3 to 405.2 € ha−1 year−1 (Table 4). Nevertheless, the results are expressed in € ha−1 of LIF area. Since the conditions are characterized by a mean LIF area share of 17.4%, the ADGMt values must be extrapolated to the total crop area for a direct comparison with the results of Robertson et al., (2007, 2008). In terms of total crop area, the obtained ADGMt ranged from 7.2 to 70.5 € ha−1 which is more in agreement with the literature. The average DGM was equal to 12.2 € ha−1 year−1 for the baseline-scenario (S-1) and, considering all the scenarios explored, the average DGM was 61.1 € ha−1 year−1, which is very much in line with the range of values reported by Robertson et al., (2007, 2008).

Seasons characterized by a negative ADGMt showed lower annual revenues than costs. Annual costs are often characterized by larger fixed costs than variable costs, because the variable costs are insensitive to farm size and structure (Pedersen et al., 2021). When the annual revenue did not overcome the 90 € ha−1 year−1 threshold (Appendix—Table 6), the net margin was negative. This was observed in three out of six years, indicating a risk of economic losses caused by VRA adoption in half of the years investigated. However, since both Ya and YG were significantly different between zones for those same years, this could still justify the adoption of VRA from both an agronomic and environmental perspective (Mulla & Schepers, 1997; Pathak et al., 2019; Plant, 2001).

From a financial point of view, the viability of VRA is strongly dependent on the annual sown area, which depends on how farmers value the return on capital. Under current price conditions (S-1 and S-2), a convincing IRR was only obtained for annual sown areas above 500 ha year−1 (Appendix—Table 7). The IRR can be profitable for sown areas between 125 and 500 ha year−1, but it would depend on changes to product prices (see S-4, S-7 and S.10 in Appendix—Table 7). The results here presented assume that farmers would purchase both the VRA equipment and the tractor (see Appendix—Table 6). This may not be the case in other regions where renting or subcontracting might be a better option, particularly for small farms. In those cases, different returns on capital could be expected as the initial investment costs would be substituted by additional annual variable costs associated with the renting or outsourcing of the VRA service.

Methodological considerations and practical issues

VRA may increase farmers’ profit by reducing costs or increasing the value of production, because fertilizer rates can be both increased or decreased between differentiated zones. Under rainfed conditions, inter-annual climatic variation leads to considerable asymmetries in crop yield patterns and financial returns on VRA. Some years show advantages for increasing N rates downslope, some benefit from decreasing applications, and others from applying N uniformly. It is assumed that, independent of the year type, the analysis succeeded well in modelling the expected marginal returns on VRA because the YGzt magnitude was computed as a module. This is valid for both years of reduced or increased N rates in LIF zones. In both situations, the differential application rates are captured by the present methodology.

This analysis assumes that N application rates approach the crop net N requirements, considering that most N inputs are recovered in the harvested grain. Considering that mean Ya values range from 3.1 to 4.5 Mg DM ha−1 (Table 3), crop net N-requirements would be 88–126 kg N ha−1. For this N use efficiency was taken into consideration as 25–40% more N is applied on average (Table 2). However, it is considered that this is in line with the ‘characteristic operating space’ for N use efficiency that is found in literature addressing European commercial (cereal) farms (Panel, 2015; Quemada et al., 2020; Silva et al., 2021).

This is an ‘ex-post analysis’ which may be a disadvantage for guiding decisions under new season conditions (Bullock & Lowenberg‐DeBoer., 2007). This study dealt with historical data (i.e., pre-collected information) and obtained from six consecutive seasons which may not be sufficient to capture the entire variability of rainfed systems. In reality, farmers face much greater uncertainty because decisions need to be taken according to specific year conditions which are highly variable. In addition, it is difficult to follow the approach used hereas ‘in-season’ YG assessment because it requires post-harvest yield information.

Therefore, the following practical question arises: how could farmers manage N applications under VRA, when their decisions must be taken without full access to the same kind of information as here presented?

It has been shown by Tenreiro et al. (2022) that the net yield response to LIF varies substantially from year to year. This makes decisions on VRA adoption challenging because farmers face possible asymmetries in crop responses patterns to differential N rates. However, some important guidelines can be proposed from these results. Under the study conditions, the maximum differential YG was 0.83 Mg DM ha−1 (Fig. 3b), corresponding to a maximum differential N rate of ± 23.2 kg N ha−1 (i.e., considering crop N requirements equal to grain uptake and assuming a N concentration equal to 2.8% DM−1 grain, Quemada et al., 2016). This represents a variation of approximately ± 16% N ha−1 between zones. In addition, it was suggested that the cost of opportunity for VRA adoption tends to decrease with the LIF coefficient (i.e., LIF divided by season precipitation). The larger the fraction of LIF contribution to total season water supply, the greater the economic advantages of VRA adoption could be. These results have practical implications for nutrient management in areas of undulating topography. In this sense, the following recommendations are offered:

  1. 1.

    Use VRA for first fertilizer application, applying up to + 8% more N in ‘LIF zones’ in comparison to ‘No-LIF zones’.

  2. 2.

    Adjust that pattern on second top-dressing applications, conducted at dates prior to flowering and according to the following criteria:

    • If the LIF coefficient falls within the top 25% near flowering (i.e., if LIF coefficient is larger than 25%, according to Tenreiro et al. (2022) findings), up to + 8% more N could be applied in ‘LIF zones’ in comparison to ‘No-LIF zones’.

    • If the LIF coefficient falls within the bottom 25% at flowering date (i.e., if LIF coefficient is smaller than 15%, according to Tenreiro et al. (2022) findings), the application pattern could be inversed by lowering N rates to a minimum of -8% less N in ‘LIF zones’, in comparison to ‘No-LIF zones’.

    • For seasons with a LIF coefficient ranging from 15 to 25% at flowering, a plausible recommendation would be to not vary the N rate for top-dressing applications.

The results are largely conditioned by the yield simulations of Tenreiro et al. (2022) and must not be directly extrapolated to other cases without further cautious considerations. For many farms, either small sized or presenting low yield variations among zones, none to minor economic advantages associated with VRA are expected for N management. It is therefore essential to farmers and advisers to take into consideration the scale effects addressed here before promoting a technological shift of this kind.

Conclusions

This study demonstrated how the relative economic advantage of VRA adoption in rainfed wheat systems of undulating topography would change greatly from year to year and from farm to farm. Both farm size (i.e., annual sown area) and topographic structure (influencing the redistribution of water from high to low parts of the fields) impacted the dynamics of investment returns. Considerable effects of scale were observed and the minimum area for adoption varied widely among different economic scenarios. The study suggests that there are economic opportunities for N management through VRA as a strategy for bridging yield gaps at the intra-plot level, which are caused by lateral inflows from high to low elevation parts of a field. In the case considered here, VRA adoption shows, for the current policy-prices scenario, that VRA adoption would have an economic advantage in farms with an annual sown area greater than 567 ha year−1, which is considerably larger than typical European cereal farm sizes. The profitability of adopting VRA is expected to respond largely to future market prices, and, in the absence of additional policy supports, the minimum area for adoption of VRA could decrease to a range of 68–177 ha year−1, depending on the price increase scenario. The effects of policy support on VRA adoption were also investigated with additional payments on crop area being the most promising from both public and private interest perspectives. The combination of further price increases and an additional payment on crop area could lower the adoption threshold down to 46 ha year−1, making VRA technology economically viable for a much wider population of farmers. Over the total amortization period, the (mean) differential gross margin of this case study that is attributed to VRA adoption was 12.2 € ha−1 year−1. Nevertheless, considerable inter-annual variation is expected, and farmers might experience net financial losses in some specific years.