Assessing sampling designs for determining fertilizer practice from yield data Computers and Electronics in Agriculture

Many farmers sample their soil to measure the concentrations of plant nutrients, so as to decide how much fertilizer to apply. Now that fertilizer can be applied at variable rates farmers want to know whether maps of nutrient concentration made from grid samples or of ﬁeld subdivisions (zones within their ﬁelds) are merited: do such maps lead to greater proﬁt than would a single measurement on a bulked sample for each ﬁeld when all costs are taken into account? We have examined the merits of grid-based and zone-based sampling strategies over single ﬁeld-based averages using continuous spatial data on wheat yields at harvest in six ﬁelds in southern England and simulated concentrations of phos- phorus (P) in the soil. We have taken into account current prices of wheat, P fertilizer and sampling and laboratory analysis. Variograms of yield provide guides for sampling. We show that where variograms have large variances and long effective ranges grid-sampling and mapping are feasible and have large probabilities of being cost-effective. Where effective ranges are short, sampling must be dense to reveal the spatial variation and be expensive, and variable-rate application of fertilizer is likely to be impracti-cable and almost certainly not cost-effective. We found zone-based sampling was less likely to be cost effective in a similar situation when the management zones were poorly correlated to P concentrations. Copyright (http://creativecommons.org/licenses/by/4.0/).


Introduction
We have known for more than 150 years that shortages of phosphorus and potassium in the soil limit crop growth. Thousands of experiments have been done to estimate the responses of crops to additions of these nutrients and to calculate the needs for fertilizers. Farmers now want to use these results to vary their applications within fields. Perhaps surprisingly, there are few reports linking variation in the concentrations of these elements in the soil to yields within individual fields on commercial farms. There are examples, however, where positive correlations were found for cereals (Frogbrook et al., 2002) and pastures (McCormick et al., 2009;Serrano et al., 2011). Many farmers in the United Kingdom sample their soil every four years to measure the nutrients, in particular phosphorus (P) and potassium (K), in the soil so as to decide how much fertilizer to apply to their crops. Sampling is often done at points on a 'W' shape across each field, and then individual samples are bulked before analysis in the laboratory (PDA, 2011). Even though this sampling configuration does not follow the principles of design-based statistics, it has been widely adopted by farmers as the results are generally no less accurate than those obtained from stratified random sampling (Marchant et al., 2012). By bulking the sample, however, all information on the variation of the nutrients across the field is lost, and so any local deficiency or excess is obscured. If a farmer wants to map the variation in nutrients across a field, so that fertilizer could be adjusted spatially, for example, then he or she should ideally sample the soil on a grid (with perhaps some additional points at closer spacings) and measure the nutrient content in each sample of soil separately (Mallarino and Wittry, 2004;Sawchik and Mallarino, 2007;Fu et al., 2013). Kriging, which makes best use of such data, could then be used to map the variation in nutrients (Kravchenko, 2003;Webster and Oliver, 2007). To krige, however, one needs an accurate estimate of the variogram or covariance function for the variable of interest, and for that at least 100 measurements are needed (Oliver and Webster, 2014). This creates a problem because measuring the concentration of P or K of each sample in the laboratory is costly.
Often the reason for the variation in yield across a field will be obvious to the farmer. For example, the farmer might know that a particular part of the field is prone to drought and that this is a major cause of the variation in yield. If the farmer suspects a local nutrient deficiency then it would make sense to divide the land into management zones and estimate the nutrient status for each zone separately.
So, which of these three sampling approaches, commonly used by farmers and advisors, should be adopted in any particular situation to apply P and K fertilizer spatially? It is a question taxing agricultural advisers who want to advise farmers on best practice for within-field sampling for plant nutrients-see, for example, Oliver and Kerry (2013), Mylavarapu and Wonsok (2014) and Hawkins et al. (2016). The grid-based approach should result in the most accurate prediction of fertilizer requirement (Mallarino and Wittry, 2004;Sawchik and Mallarino, 2007). But the money saved by varying the application of the fertilizer locally within fields to match the requirement of the crop might be less than the cost of sampling and measurement (Fleming et al., 2000;Mallarino and Wittry, 2004). The balance of the two, and therefore the merit of the approach, depend on the magnitude and variation of the nutrient content of the soil (Sawchik and Mallarino, 2007). These two variables can be determined accurately only from measurements made on samples. The variation can, however, be assessed indirectly from crop yields. Many farmers are already monitoring yields as they harvest their crops, and variation in the data they record is in many instances a reflection of the variation in the availability of the nutrients in the soil (Stafford et al., 1999;Diker et al., 2004;Flowers et al., 2005). If the variation in yield is small then it is unlikely that the nutrients vary substantially. We know that factors other than nutrient supply can cause large variations in yield. Nevertheless, nutrient supply does dominate yield variation in many cases, and for present purpose we proceed on that assumption. In such circumstances Lark et al. (2003) proposed metrics based on the variogram of yield data (which captures the magnitude of variation of the yield) to assess the scope for variable rate management. Their approach was to use the metrics as factors in a decision tree designed to determine the potential for variable rate management.
In this context we aimed (i) to compare the merits of measuring plant nutrients by three sampling schemes (whole fields, zones within fields and grids) and (ii) to assess the extent to which yield maps might be used to determine the most cost-effective sampling. Which of the above sampling approaches is suitable for a given situation depends on a farmer's profit margin over the cost of fertilizer and soil sampling. It is not possible to do such a comparison in the field because the test requires perfect knowledge of how the nutrients vary across the field, we therefore resorted to simulation. In the approach presented here we simulated the variation in nutrients across fields from geostatistical models of the nutrients and used these to test the cost effectiveness of each sampling scheme. We modelled the associated yield variation for each realization and explored the use of the metrics of the yield variogram to decide which sampling strategy was likely to be most cost-effective.

Data
We collated yield data, denoted y, from monitors on board combine harvesters and measurements of extractable P, denoted z, in the soil for six fields on a farm near Newbury, England (Table 1). Soils were medium to heavy textured with slight to moderate stoniness. The yields were of winter wheat from the seasons 2001 to 2011 recorded at approximately 20-m intervals by the monitor. The measurements of Olsen extractable P (sodium bicarbonate extract at pH 8.2) were made on a 100-m grid at 24-36 locations across each field. Fig. 1 displays yield maps for the six fields for a single year.

Determining management zones
Several methods have been devised for creating management zones (Oliver and Webster, 1989;Fleming et al., 2000;Diker et al., 2004;Flowers et al., 2005;Zhang et al., 2009). The multivariate technique of Dray et al. (2006) has most recently been applied successfully by Peralta et al. (2015) for wheat farming and by Córdoba et al. (2016) for grain cropping more generally.
We created management zones from the yield data using a spatially smoothed version of a fuzzy k-means classification devised by Lark (1998) (see also Milne et al., 2012). The data for the classification consisted of yields of wheat for p years at the n nodes of a square grid at intervals of 10 m. We denote the grid coordinates as x x 1 ; x 2 f gand the yields in the p years as y 1 ðxÞ; y 2 ðxÞ; . . . ; y p ðxÞ. From these data we created a classification by a 'hard' k-means algorithm. We standardized each of the y j ; j ¼ 1; 2; . . . ; y p to zero mean and variance of 1, denotedỹ j . For this method we choose k, the number of classes. We divided the whole set of the standardized data into that number of classes in such a way as to minimize the trace or determinant of the within-classes variance-covariance matrix. Each grid node then belonged to one and only one of the k classes, and in general it resembled other members of its class more than the members of the other classes.
For fuzzy k-means classification we first define a measure of dissimilarity, d, between an individual node i and a class q. A convenient measure is the Euclidean distance in the vector space: whereỹ ij is the standardized yield at node i in the jth year, andỹ jq is the mean ofỹ in class q in that year.
Nomenclature x x 1 ; x 2 f g spatial coordinates in two dimensions y yield of crop y standardized yield y r realized yield y 0 target yield z quantity of phosphorus, P z Ã realization of z z fert quantity of fertilizer P z soil initial quantity of P in the soil z total z fert þ z soil k the number of classes in the k-means classification variogram parameters c 0 nugget variance c 1 variance of spatially correlated structure a distance parameter Costs G wheat price of grain, assumed to be £150 t À1 G fert price of P fertilizer, assumed to be £0.31 kg À1 G sample cost of soil analysis, assumed to be £5 sample À1 For this method we assume that each node belongs to some degree to every class, and we create a classification by minimizing a pooled 'belongingness', say b: in which u iq is the degree of membership of node i to class q, and x is the fuzzyness parameter. The membership across all classes must sum to 1: The parameter x must lie between 1 (in which case we obtain a hard classification) and 2. We set x ¼ 1:25 to create our classifications.
As above, we needed to choose k. We did this by experimenting with several values between 2 and 5 (the most that a farmer is likely to distinguish). For each class we then computed the normalized classification entropy nðkÞ, proposed by Dunn (1977): We then ploted nðkÞ against k and identified the value of k at which nðkÞ falls below the overall trend. This was the value we chose. The  procedure is analogous to that in hard k-means classification in which the trace of the within-classes variance-covariance matrix or its determinant is plotted against k. The next step in the zonation is to smooth the classes. It turned out that the distributions of the memberships of the nodes were strongly bimodal, and so, following Lark (1998), we transformed them with the symmetric log-ratio to unimodal distributions. We smoothed the transformed membership (denotedũ iq ) using a weighted average of the transformed memberships in circular neighbourhoods, R, of radius r: Like the original memberships, the transformed memberships must lie in the range 0-1 and must sum to 1. This means that the weights in R must sum to 1. The weights were derived from the variogram, Eq. (6). We can write a simple bounded model in general as In this equation c 0 is a spatially uncorrelated variance, the 'nugget variance', corresponding to white noise and c 1 is the spatially correlated component of variance; their sum, c 0 þ c 1 , constitute the 'sill'. The function f ðhÞ defines the form of the variogram and contains a distance parameter. The weights are obtained as where h ij is the separation in distance and direction, the lag, between nodes i and j. Note that only the functional form of the variogram and its distance parameter affect the weights; neither c 0 nor c 1 do so. The neighbourhood R defines the region over which the membership values are smoothed, unless the variogram reaches its sill within it. In the latter case the effective range of the variogram defines the smoothing region. The farmer, of course, must have a hard classification; for practical management he must have each position in the field belonging to one class and one class only. So the final stage in the zonation is therefore to assign each node to the class for which its smoothed membership is greatest.
Note that the size of R affects the results. The larger it is the greater is the smoothing. If R is small then the classification is likely to be too fragmented; if it is too large then the memberships will be smoothed too much and the final classes not sufficiently homogeneous. Lark (1998) proposed a coherency index to identify an appropriate radius for R defined as where g a is the proportion of pairs of nodes within a distance a that belong to the same class, and w q is the proportion of nodes that belong to class q. The larger is the value of H, the more spatially coherent are the classes. We chose a ¼ 10 ffiffiffi 2 p m (the distance between two points on the diagonal of our 10-m grid) so that we were effectively comparing each node with its neighbours along the rows, columns and diagonals on the grid.

Modelling phosphorus in soil
For each of the six fields, we used the measurements, z, of P to create realizations of P concentrations with plausible means and spatial variations that could differ from one management zone to another. Our data came from well managed fields, but we wanted our simulated values of P to limit the yield. Therefore we scaled the measured data by 0.5 before fitting the models of spatial variation for all fields except Mantclose. We used a similar approach to that described by Marchant et al. (2012). First we standardized the measurements, z i ; i ¼ 1; 2; . . . ; n, to have a variance = 1 by dividing by the standard deviation, s, of the data to give values z i ; i ¼ 1; 2; . . . ; n. Then we characterized the mean and spatial variation within each of the k zones of the field by fitting a linear mixed model to the transforms,z: where M is an n Â k fixed effects design matrix which permits the mean concentrations to differ between zones. If the ith value ofz is in zone j then the element in ith row of the jth column of M is 1, otherwise it is 0. The vector b is of length k and contains the coefficients of the fixed effects (i.e. the mean concentration within each zone). The component g $ UNð0; VÞU where Nð0; VÞ is a vector of spatially correlated random residuals with a normal distribution with zero mean and covariance matrix V, and U is a diagonal matrix this is a zone-dependent scaling factor of the variance). If we assume second-order stationarity, so that the covariance function exists and can be obtained from the variogram parameters, then these and the fixed effects coefficients b can be estimated simultaneously by residual maximum likelihood, REML, (Patterson and Thompson, 1971). We assumed that the spatial variation is represented by an isotropic exponential variogram model: in which c 0 and c 1 are the nugget and spatially correlated components of the variance, as mentioned above, and a is a distance parameter. The function approaches its maximum, c 0 þ c 1 , asymptotically, and the distance 3a is often taken to be the effective range of the spatial correlation (see Webster and Oliver, 2007). In the discussion below we shall often refer to the 'effective range' with this meaning.
We simulated values for P on a 10 m Â 10 m grid across each field using the Cholesky decomposition technique, also known as lower-upper or LU technique (Webster and Oliver, 2007). We used the variogram model that we fitted to the transformed data to create an t Â t covariance matrix C and scaled this for each zone inde- where t is the number of simulated points on the 10 m Â 10 m grid and b U is a diagonal matrix where element b Uði; iÞ ¼ r j when the ith simulated value is in zone j. This was then decomposed into its lower and upper triangular form where The simulated values, z Ã , are then given by where g is a t Â 1 vector of random numbers drawn from a standard normal distribution, and M sim is a t Â k design matrix. If the ith datum is in zone j then the element in the ith row of the jth column of M sim is 1, otherwise it is 0. Fig. 2 shows examples of the simulated values of P. For ease of calculating the yield response (see Section 2.5) we converted our simulated P values to mg kg À1 by assuming the soil has a bulk density of 1.1 g cm À3 .
To simulate field measurements with larger variances than in the observed data we also scaled the values of r j and simulated values of P. In each additional zone, concentrations of P were modified at four scales (0.5, 1, 2, 3) resulting in four sets of runs for fields with two zones and 16 sets of runs for fields with three zones.

Modelling yield and quantifying its spatial variation
For each realization of simulated phosphorus, z Ã , we simulated the associated yields (y) using yield response models for P (see Section 2.5) and computed experimental variograms from simulated yield values by the method of moments: where yðx j Þ and yðx j þ hÞ are the simulated values at positions x j and x j þ h separated by the lag h, which in the isoptropic case is the scalar distance h, and mðhÞ is the number of comparisons at that lag. By changing h we obtained the ordered set to which we fitted exponential models, Eq. (10), by weighted least-squares approximation.

Yield response model
The yield response model for P was derived by Marchant et al. (2012) from published data (Syers et al., 2008;Johnston, 2005;Milford and Johnston, 2007;Johnston and Goulding, 1988). For every 1 kg of P added in fertilizer, we assumed that 0.18 kg is available to the crop. We also assumed that this addition is contained in the top 30 cm of soil and the soil has a bulk density of 1.1 g cm 3 .
This means that an addition of 1 kg P ha À1 leads to an increase in the concentration of this layer of j ¼ 0:054 mg kg À1 . Thus the total nutrient available after addition of a quantity of fertilizer z fert was where z soil is the initial concentration of the nutrient in the soil.
The yield response to added nutrients was modelled by where y r is the realized yield, y 0 is target yields and A and B are parameters. We set our target yield at 8.8 t ha À1 . The model parameters were A ¼ 1:33 and B ¼ 0:68 (see Fig. 3).

Sampling strategies
In each of the simulated fields, we compared three sampling schemes to see which would result in a treatment map that gave the greatest profit. The sampling schemes we considered were (i) a W-shaped design across the whole field, (ii) a zone-based scheme with W-shaped designs within each zone and (iii) a grid-based scheme with samples taken on a 100 m Â 100 m grid. Each W-shaped design comprised 10 sampling points. In practice soil samples taken according to schemes (i) and (ii) are bulked before analysis resulting in either a single value for each field or each zone within a field. To simulate these we calculated the average of the nutrient values from the sample points on each Wshaped design. In grid-based designs the aim is to map the variation in a nutrient so that fertilizer rates can be adjusted accordingly. Therefore these samples are not bulked before laboratory analysis; rather a measurement is made on each sample. The number of sampling points depends on the size of the field, but typically there are too few for prediction by kriging (the fields we considered here vary from 20 to 30 ha in size). The reason is that a minimum of about 100 points is needed for an accurate estimate of the variogram. Therefore inverse distance weighting is often used for the predictions instead of kriging, but these predictions are likely to be sub-optimal because practitioners do not know the spatial scale(s) of the variation.

Evaluation of sample designs
The profit from applying any of the strategies we have set out can be assessed simply as the difference between the cost of the measurements plus that of the fertilizer and price of the crop that is produced, in this case wheat. This is the net gain, which we denote D for unit area, thus: where G wheat is the price of the grain, which we assumed to be £150 t À1 ; G fert is the price of fertilizer P, assumed to be £0.31 kg À1 ; n is the number of individual soil samples analyzed in the laboratory, each costing G sample ¼ £5, and y and z fert are the yield and quantity of fertilizer added, respectively.
We estimated the nutrient concentration in the soil at each location on the 10 m Â 10 m grid using each of the three sampling schemes. This resulted in a single estimate for the whole-field sampling scheme, an estimate for each zone for the zone-based scheme and a spatially varying estimate for the grid-based scheme. We denote these estimates of z soil by b z field , b z zone and b z grid . Using these estimates with Eqs. (14)-(16) we calculated the amount of fertilizer that should be added (z fert ), noting that this value is not the true optimum as it is based on the estimated nutrient supply not the true nutrient supply, z soil .
For each realization, i, of the fields we calculated the profit margin under each sampling scheme. We computed the difference in profit margin given by the zone-based, D zone ðiÞ, and grid based, D grid ðiÞ, schemes compared with that of the field-based D field ðiÞ scheme. We also computed the excess fertilizer applied under the zone-and grid-based schemes and compared this with the field-based scheme.

Using metrics of variation to guide sampling strategies
For each set of simulations, we used multiple linear regression to see how much of the variation in D grid À D field and D zone À D field could be explained by the distance parameter of the yield variogram a and the variance parameter c 1 . More importantly we wanted to compute the probability that the grid-or zone-based sampling strategies were more profitable than the field-based strategy for given parameters of a and c 1 . For each field we fitted the model to the data, where 'scheme' is 'zone' or 'grid'. The assumption underlying the model is that the residuals are normally distributed about the mean prediction and that the standard error, s obs for predicting a single observation is given by where VðbÞ is the covariance function for the parameter estimates b fb 0 ; b 1 ; b 2 ; b 3 g and s 2 mse is the mean square error. From this we calculated the probability that D j À D field > 0 for any given combination of a and c 1 . We did a similar analysis for excess fertilizer.

Results
Based on the normalized classification entropy, two or three management zones were identified for each field (Table 2). Table 2 lists the parameters of the models fitted to the nutrient measurements, Eq. (9), and the fitted variograms are shown in Fig. 4. The sill variances, c 0 þ c 1 , vary from 0.09 to 2.23, and the effective range from 63.0 m to 522.6 m showing substantial differences in the spatial structure of P between the fields. Yield / t ha -1 4 5 6 7 8 9 Fig. 3. Response of yield to added P to a soil with an Olsen extractable P of 2.5 mg kg À1 . Table 2 The number of management zones derived for each field and the parameters of the linear mixed model (Eq. (9)) fitted to the measured values of P. In each case r1 ¼ 1.

Comparison of P estimated by the sampling schemes
The range of the simulated field averages for P (Table 3) were consistent with those reported by analytical laboratories for arable fields in England and Wales (NRM Laboratories, 2017) and within the measured range (Table 1). For each realization we computed the difference in profit margin given by the zone-based (D zone ) and grid based (D grid ) schemes compared with that of the  field-based (D field ) scheme. Thus, the differences are (D grid À D field ) and (D zone À D field ) for the grid-based and zone-based values, respectively. Table 3 reports the means and standard errors for these differences, which are also shown in Fig. 5. In all cases the grid-based estimates gave larger profits than did the zone-based estimates. This was largely because the cluster classes were not significant factors in explaining the variation in the nutrients (see supplementary information). The smallest differences between the zone-based estimate and the grid-based one were for the fields with the largest differences in mean concentrations of P between zones (Gravelpit and Roodhill). Bassdown had the smallest mean of D zone À D field , and this is likely to result from both the small difference in mean values of P between the zones, sbð1Þ; sbð2Þ and sbð3Þ, and the short-range variation ( Table 4 The percentage variance accounted for and model parameters for the multiple linear regression models (Eq. (17)) fitted to the Dzone À D field and D grid À D field data.
Field name Model parameters Percentage variance accounted for Zone sampling range 63.0 m) which is smaller than the distance across feasible management zones. The largest mean profits were for Limes and Mantclose which have a long spatial structure (Fig. 4) with values of P in a treatable range (Table 3 and Fig. 3). The distributions of D zone À D field are more symmetrically distributed than for D grid À D field for all fields except for Roodhill, (as illustrated in Fig. 5). The distributions of D grid À D field are positively skewed. For Roodhill and Gravelpit a large number of realizations had values of b z field and b z zone that were not limiting (99% and 66% respectively). This resulted in recommendations of no fertilizer application, and so D zone À D field is simply the difference in sampling costs. These correspond to the large peaks in the distributions of D grid À D field .
We also computed the difference in excess fertilizer applied when estimates were based on the zone-based (b z zone ) and grid based (b z grid ) sampling schemes compared with the field-based (b z field ) scheme. We define excess fertilizer as the amount applied over and above that which would have been applied if we had perfect knowledge of the true variation in P across the field. Table 3 reports the means and standard error for these differences (b z zone À b z field and b z grid À b z field ).
There was no consistent pattern to the mean responses. In some fields (Easton Right, Gravelpit and Roodhill) the field-based sampling resulted in less of an excess than the zone-based sampling with positive differences and in some fields (Bassdown, Gravelpit and Roodhill) the field-based sampling resulted in less of an excess than the grid-based sampling (Table 3). Similarly there was no consistent pattern between grid-and zone-based sampling.
3.2. Assessing the extent to which yield maps can be used to predict the most appropriate sampling scheme The parameters for the models fitted by multiple linear regression are listed in Table 4  for. For Gravelpit, Mantclose and Roodhill the regression model was fitted to the subset of realizations where b z field or b z zone were limiting. The models fitted to D zone À D field explained very little of the variation in the data. The variation in D grid À D field was better explained, although the value for Bassdown were still small. The realizations for this site were generated from the model with a small effective range and the nugget to sill ratios of the realizations were in general larger than other sites (at least 37% larger, data not shown) indicating a relatively large component of unstructured variance.
Figs. 6a and 6b shows the probability that D zone and D grid are larger than D field . In all cases the probability increases with both c 1 and the effective range (3a). This is also true of D zone À D field for Gravelpit and Roodhill, but these sites are dominated by simulations where b z field and b z zone where not limiting, and so the effects are negligible.
The multiple linear regressions showed that the effective range and c 1 explained very little of the variation in the excess fertilizer data, b z zone À b z field and b z grid À b z field (Table 5).

Discussion and conclusion
We have investigated, through simulation, the costeffectiveness of three sampling strategies commonly used to guide fertilizer recommendations. We aimed to see if the variation captured in yield monitor data could be used to determine which sampling strategy would be best in any given situation. The first option was the conventional single bulked sample for a whole field. The second was to delineate spatially coherent management zones from the yield data, sample the soil of each zone independently (zone-based sampling) and apply fertilizer (P in our case) in accord with the P content of the soil measured. The third option was to estimate the magnitude of the variation of P across the field, and judge from this how fertilizer might be applied at a variable rate. This required grid-based sampling.
Despite the expected increases in yield in response to additions of fertilizer, we know that there are situations in which there are inverse relations between nutrient concentration and yields (see, for example, Lake et al., 1997, for wheat andCox et al., 2003, for soya beans). It seems that in those situations the larger crops have extracted more of the nutrients and depleted the stock in the soil.
We have shown that the advantages of using grid-and zonebased sampling strategies over field-based ones vary from field to field. In our simulations, on average, the grid-based sampling performed better than the zone-based sampling. This was largely because of the variation in the concentration of P, and the mean concentrations, did not differ sufficiently from one zone to another (even after scaling down the P measures used to model the variation in soil P by 0.5). This is likely to be so for many arable fields in the UK because so much P has been added to them in the last 150 years. Other investigators (Mallarino and Wittry, 2004;Flowers et al., 2005;Sawchik and Mallarino, 2007) similarly found that grid-based sampling was better than the zone-based approach for estimating the likely responses of field crops to added fertilizer. Factors such as changes in water availability (due to variation in drainage or soil texture) or local infestations of weeds or disease will often be the main cause of large variations. In practice, many farmers would be able to explain the observed differences between proposed management zones and so be able to predict whether the zones were determined by differences in nutrient availability. This would provide valuable information on whether zone-based sampling was sensible.
The maps in Fig. 6a and 6b show that the probability of gridbased sampling's being more profitable increases with both increases in effective range and in c 1 . Larger values of c 1 imply large differences in nutrients, and a farmer might wish to apply fertilizer differentially in accord. This is feasible in practice only if the effective range is also large, and for two reasons. One is the difficulty of varying the application at a fine scale; the other is the cost of sampling and soil analysis on grids fine enough to map the concentration of the nutrient in the soil.
Notice that the probability maps in Fig. 6a and 6b change from field to field. Note also that they are based on current prices of wheat and fertilizer and sampling costs, which are subject to variation, and that our models might have introduced bias in profitability. Indeed, farmers will typically measure several soil variables at once (e.g. P, K, Mg and pH). Soil sampling may reveal anyone of these variables (or a combination of) to be the dominant cause of yield variation. Accurate quantification of the benefits of one sampling approach over another are therefore difficult, and so the absolute values shown in Figs. 6a and 6b should not be applied in other contexts, although, the principles hold. As variation in yield at scales appropriate for management increases, so does the likelihood that grid-based sampling will be more profitable than a single field estimate.
Farmers and their advisors are building increasingly large data sets of crop yields, soil properties and previous management decisions. They should find these data valuable for managing their fields more effectively both to improve productivity and to reduce unnecessary inputs. The principles described here could be integrated into a decision-support system to help farmers decide when Table 5 The percentage variance accounted for and model parameters for the multiple linear regression models (Eq. (17)) fitted to theẑzone Àẑ field andẑ grid Àẑ field .

Field name
Model parameters Percentage variance accounted for variable-rate application of fertilizer might be cost effective and if so what sampling strategy to apply. This is the first step, then the farmer must decide how to vary the fertilizer spatially. The reliable prediction of how much fertilizer to apply depends not only on the quantification of the soil nutrient supply but also the yield potential of the crop and the efficiency of nutrient uptake by the crop (Kindred et al., 2015). Predicting yield potential and the efficiency of nutrient uptake are substantial topics of research that can involve sensor technology, empirical analysis and simulation modelling (Dhital and Raun, 2016;van Wart et al., 2013;Kindred et al., 2015). In our study we followed the practices of agricultural advisors in the UK and considered zone-based and grid-based sampling schemes. The grid-based sampling scheme was on a 100-m grid, but in practice the variogram could be further used to optimize the grid sampling according to the cost of samples and the apparent dominant scale of spatial variation in the yield data.