Burgess, T.M. & Webster, R. 1980. Optimal interpolation and isarithmic mapping of soil properties. I. The semi‐variogram and punctual kriging. Journal of Soil Science, 31, 315–331

This landmark paper by Burgess & Webster (1980a) signalled a new era in the spatial mapping 2 of the soil. The emergence of pedometrics as a distinct subdiscipline of soil science was a gradual 3 process, and had its roots in earlier studies than this one, but if one publication is to mark the 4 start of pedometrics, then this is it. 5 The publication of this paper, and the successors in the series that it initiated, showed that 6 statistics could do more than just assist the soil surveyor in quantitative prediction from maps 7 based on a legend of discrete units. It showed soil scientists that there was a very different way 8 of mapping the soil, namely through geostatistical interpolation. 9 The impact of the paper on practice has been substantial. Soil scientists have used geosta10 tistical prediction to map the distribution of nutrients or pests at within-field scale for precision 11 agriculture, to delineate contaminated land, to map changes in soil properties over time at the 12 regional scale to map the risk of problems such as salinity and even for forensic inference. 13 The value of geostatistics was soon recognized in soil science. This is, in part, because when 14 one considers the state of the art at the time of its publication it is clear that it fulfilled a need for 15 new methods. At the same time, its impact was not simply due to the method it introduced to 16 soil science, ordinary kriging, but to the interest it stimulated among soil scientists and applied 17 statisticians who went on to develop and improve the methodology in several critical ways. 18

I have told before in a tribute to Danie Krige how geostatistics came to be applied to soil (Webster, 2015), but the hard graft and fortuitous events that led to this first paper on the subject will bear repeating. It began with my appointment as soil chemist by the British Colonial Office to the Government of Northern Rhodesia, now Zambia, in 1957. The government's aim was to evaluate the suitability of land for agricultural development, and it became my job to survey and map soil for that purpose. My training beforehand in Britain was based on what was then conventional technique, namely the identification of distinct classes of soil and their delineation on maps. The notion underlying the procedure, usually implied but less often stated, was that if one could identify from a map the class of soil at a place then one should be able to predict the soil's properties there. Transferring the practice to the soil on the deeply weathered rocks on the Zambian plateaux was problematic, however. The soil varied gradually over the landscape; there were no obvious boundaries between one kind of soil and another, although there were the repetitive sequences related to the topography, for which Milne had coined the term 'catena'. Further, the fairly dense miombo woodland in the north of the country meant that one could rarely see for more than a few tens of metres, and air-photo interpretation had so far been of little help. Survey had to be carried out almost entirely by sampling, and I made my maps from the sparse punctual observations by interpolation. The last was simply by hand and eye, hardly satisfactory in an increasingly quantitative era. I wanted something better.
Enter into Zambia Philip Beckett, another sceptical chemist. He too was interested in prediction of soil conditions between observation sites. During service in the British army he had noted the successes and failures of conventional soil maps in predicting going conditions for military vehicles. Sitting around a paraffin lamp out in the bush in 1960, we returned evening after evening to the question: how can we predict from punctual data with known measures of confidence? We recognized that the problem was essentially statistical and that we needed a statistical solution. The following year I joined Philip at Oxford to pursue the matter in collaboration with the Royal Engineers.
We were not alone in our quest; several other engineers had begun to realize what the problem was and were toying with combinations of classical soil maps and statistical prediction. Morse & Thornburn (1961) sampled at random the classes of agricultural soil maps in Illinois and computed the means and variances of engineering properties within classes so as to predict values at unsampled places. Kantey & Williams (1962) in South Africa took matters a stage further; they made their own engineering soil maps and sampled with the same aims. Beckett and I planned a thorough study along similar lines. We classified and mapped a large part of the Oxfordshire landscape, sampled it to a stratified random design and measured properties of the soil. Then by analysis of variance we assessed our classification for its effectiveness in (i) diminishing the variances within the classes and (ii) predicting values with acceptable precision. We also assessed maps made and sampled by several of our collaborators. We had mixed success. Our map of the Oxford region enabled us to predict the mechanical properties of the soil reasonably well. It predicted relatively poorly the soil's pH and organic matter content, and it was useless for predicting the plant nutrients in the soil. Table 1, from Webster & Beckett (1970), summarizes our results. Nevertheless, even in the most favourable situations there was substantial residual variation for which the classification could not account. It was unlikely to be white noise, 'pathological variation' as Priestley (1981) later called such impossible variation for a continuous variable in continuous space or time. It must be structural. We had not solved the problem arising from gradual change or trend. If we simply drew arbitrary boundaries in those situations then the residuals would be spatially correlated. At about that time petroleum geologists were having some success with trend-surface analysis to map simple smoothly varying geological structures. In our situation, however, it was unsatisfactory because (i) fluctuation in one part of a region affected the fit of the surface elsewhere and (ii) the residuals were correlated so that calculated prediction variances were biased.

Soil properties as stochastic processes
The next significant step came when H.E. Cuanalo arrived in Oxford from Mexico. He pointed out that time-series analysts have similar problems; they treat actuality as realizations of stochastic processes to describe quantitatively fluctuations in time. We should be able to do the same for soil in space. What heresy! It would require a leap of imagination, away from the then current view of soil as the deterministic outcome of certain factors (and predictable if we know what those factors were) to soil as the product of random processes. Looking back, we should not have been surprised. Our knowledge of the soil-forming factors is far from certain, and their interactions are so complex that it is hard to distinguish the outcome from random (Webster, 2000).
Cuanalo set to with vigour to test the feasibility of the approach. He recorded the soil profile at 10-m intervals in 321 pits on a 3.2-km-long transect across the Jurassic scarp-lands of north Oxfordshire. We computed correlograms of the properties from the data and reported them, the first of their kind, in Webster & Cuanalo (1975). All showed strong spatial correlation extending to 200-250 m. This distance corresponded approximately to the average width of the outcrops and to the evident changes in soil. If we filtered out the variation due to the presence of the distinct underlying Jurassic strata we discovered that there was still spatial correlation in the residuals, although with a range of only about 80 m. Figure 1 shows an example, here with variograms rather than correlograms, for the clay content in the subsoil (at ≈ 65 cm deep). The upper sequence of points is of the raw data, and the curve is that of the spherical function fitted to them. The lower sequence is of the residuals after the outcrop means have been filtered out; again, the curve is of the spherical function fitted to them. The spherical function is: where c 0 is the variance at lag zero, the so-called nugget variance, c 1 is the correlated component of variance, and r is the range, the limiting distance of spatial dependence. Table 2 lists the values of the parameters of the curves shown in Figure 1. Seen from today, the situation seems obvious; we had two sources of variation, one from class to class of outcrops, which we might treat as a fixed effect, and the other within classes, which we should treat as random. We needed a mixed model to describe it, as I explained in my tribute to Danie Krige (Webster, 2015) and illustrated with some results.
To some extent north Oxfordshire was a fortunate choice for exploring spatial correlation in soil. We could see qualitatively how the soil varied and chose a sampling interval to accord with this. In the event we took a cautious approach and sampled more closely than necessary to reveal the form and scale of variation. We might have drawn the same conclusion with a coarser sampling interval and less work. But what was one to do where one could not see, as on the Zambian plateaux? How could one discover the spatial scale(s) of variation with modest labour before more intensive sampling to obtain reasonably accurate variograms? I had the opportunity to explore while working with B.E. Butler, the doyen of Australian pedology in the CSIRO. Our first task was to discover the spatial scale(s) on which soil properties varied on the Southern Tablelands of Australia, another region of deep weathering. We sampled to a balanced spatially nested design and estimated the components of variance from a hierarchical analysis of variance, quite unaware that the technique had been proposed 36 years earlier by Youden & Mehlich (1937) and lain almost forgotten in their house journal. We summed the components to form rough variograms and discovered that different soil properties varied on disparate scales in that landscape (Webster & Butler, 1976). Figure 2 shows such variograms of two of the variables. Armed with that information I could compute for the soil's phosphorus content a conventional variogram and map with it, but the grid interval of 180 m was too course for potassium (Webster, 2011). But I anticipate; we did not know how to use the knowledge at the time. The spacings in Figure 2 increase in logarithmic progression. Notice that as the sampling interval shortens (i.e. as the scale becomes increasingly fine, from right to left on the graph) the number of degrees of freedom increases twofold with each step after the first. If one wants more steps on the graph for greater refinement, then doubling the sampling at each stage to maintain balance soon becomes unaffordable. The increased precision at the shorter lag distances is also unnecessary. Margaret Oliver recognized the problem and sacrificed balance and analytical elegance for greater efficiency for studying the soil in the Wyre Forest of England. She and I designed a five-stage hierarchy but without doubling all branches of the hierarchy at the lowest stage, and we programmed Gower's (1962) algorithm to estimate the components of variance (Oliver & Webster, 1987). Shortly afterwards Boag and I devised an extreme form of unbalanced hierarchy with equal degrees of freedom at all stages apart from the first in a study of the distribution of cereal cyst nematodes in soil, and again we estimated the components of variance by Gower's method (Webster & Boag, 1992). We have since replaced Gower's method, which although unbiased is not unique, by the more efficient residual maximum likelihood method reml (Webster et al., 2006). The publications cited above and the more recent one by Lark (2011) should ensure that this efficient, economical way of obtaining a first rough estimate of the variogram in unknown territory is available for any pedometrician's tool kit.
A second topic in my Australian research was to characterize the spatial pattern of gilgai terrain. Gilgais are typically shallow wet depressions a few metres across in otherwise flat plains, and their patterns seem to be regular. I and two technicians sampled the soil to 1 m at regular 4-m intervals on a 1.5-km transect on the Bland Plain of New South Wales (Webster, 1977). The correlograms of several properties appeared wavy, and I transformed them to their corresponding power spectra. Chapter 7 in Webster & Oliver (2007) is an up-to-date account of that research.
But how could we use this intelligence to interpolate? Figure 2 Accumulated variances as proportions of the total variances for soluble phosphorus (P) and soluble potassium (K) estimated by hierarchical analysis of variance from nested sampling at Ginninderra, Australia (data from Webster & Butler, 1976). The degrees of freedom for the distances are printed immediately above the abscissa.

The answer: kriging
Perhaps the most significant event during my sojourn in Australia occurred a week or so before I was due to leave. A complete stranger breezed into my office without a by-your-leave and asked me bluntly: 'What's this kriging?'. I had never heard the term before, and rather than plead complete ignorance I played for time. Who was this intruder? Why did he ask? He was Daniel Sampey, a mining geologist. I let him talk, which he did for about 20 minutes. He told me of a certain Professor Krige and of Georges Matheron, of the theory of regionalized variables and of its application in geostatistics. Then, clearly disappointed that I knew even less than he did, he left as abruptly as he had arrived. His parting shot was that as I was about to return to Britain I should visit Leeds University where mining engineers knew a thing or two. In those 20 minutes I realized that my problem of spatial prediction of soil conditions at unvisited places had been solved, at least in principle, and in general terms I understood how. On my return to Britain I contacted Anthony Royle at Leeds. He amplified what Daniel Sampey had told me, and he generously gave me a copy of his lecture notes on the subject and some references to the literature, including to Matheron's (1965) seminal thesis.
On my return I recruited Trevor Burgess, who was keen to apply his newly won mathematics degree to the task of kriging. Together we turned Matheron's equations into algorithms and algorithms into computer code. From colleagues we obtained data on several soil variables from surveys in mid-Wales and Norfolk and plugged them into our code, first to estimate the variograms and then, having modelled the variograms, to interpolate by ordinary kriging. Finally, we contoured the kriged predictions to map the variables. We submitted our papers (Burgess & Webster, 1980a,b) describing our achievement to the Journal of Soil Science, and they appeared shortly afterwards. They were the first to describe for soil scientists the variogram as we know it today and the first to display maps of soil properties made by kriging.
On re-reading those papers for this issue of the journal I was surprised to discover how deep was our understanding at that time. There in black and white are accounts of the underlying theory and mathematics. They have been the basis of my teaching of geostatistics ever since. More worryingly, I have had to explain them time and time again in critical reviews of papers sent to me by the editors of learned journals. I wish that authors would read those papers before simply pressing buttons on computer packages and copying uncritically the results into their scripts.
I return to my reflections, and to more good fortune for us as soil scientists. Soil is almost the ideal medium for practical geostatistics. It forms a continuous mantle over large parts of the earth's land surface. Access is easy over much of that, so that sampling at the working scale of the individual field, farm or estate can be cheap. We are not confined to tunnels underground, as gold-miners are; our targets do not move, as fish and the weather do. Some of the soil's most important properties, such as its pH, concentrations of the major plant nutrients and trace elements, salinity and pollutant heavy metals, are also cheap to measure nowadays: we need not be short of data for these variables. With so many data we can estimate spatial covariance functions and variograms accurately. The statistical distributions of these properties are in most instances 'well-behaved' in that they are either close to normal (Gaussian) or to lognormal. In the latter situations simple transformation to logarithms makes analysis straightforward and efficient. Further, although the laws of physics must be obeyed as soil is formed, the numerous processes that operate and have operated in combination over many millennia to form the present-day soil have produced a complexity that is indistinguishable from random (Webster, 2000). So, we can often treat it as the outcome of random processes without harming our professional reputations.
It is a small step from there to assume that the soil variables, transformed if necessary, are intrinsically stationary and that the variogram is a sufficient summary of the spatial structure. Ordinary kriging follows, and from the 1980s onwards it has become the workhorse of geostatistics in surveys not only of soil itself but also in the related fields of agronomy, pest infestation and pollution; there are hundreds of examples of its application described in the literature. It has proved to be valuable in modern precision agriculture in particular (see Oliver (2010) and the topics therein).
There have been advances and refinements, of course: science and technology do not stand still. One in particular seems to me especially important; it is the treatment of geographic trend, or 'drift' as it is usually known among geostatisticians, and the mathematically related 'external drifts' (i.e. related variables). If there is drift then the assumption of intrinsic stationarity is violated. The kriging itself requires no more than an augmentation of the ordinary kriging system, which Matheron (1969) called 'universal kriging'. The difficulty, and a serious impediment for a couple of decades, was obtaining a valid estimate of the variogram of the random component of the variation (i.e. of the residuals from the drift). The way to overcome this difficulty has been to use likelihood methods, as pointed out by Stein (1999), and residual maximum likelihood (reml) in particular. Lark and colleagues have been at the forefront of this development to obtain what they call the empirical best linear unbiased predictor, or 'e-blup' (Lark et al., 2006). I leave Lark et al. (2019) to expand on this in their commentary.

Addenda
While reflecting I mention one other way in which I was fortunate. All of the popular valid variogram functions, apart from the unbounded linear model, contain non-linear distance parameters, and these cannot be estimated by ordinary least-squares regression. Some, such as the exponential and power functions, can be re-parameterized so that they are linear. Others such as the spherical and related functions cannot; they must be estimated by numerical approximation, and doing that requires expertise in numerical analysis. Rothamsted had that expertise; Ross (1987) had written his program, MLP, for non-linear estimation, and we soil scientists used it to advantage for fitting models to experimental variograms (McBratney & Webster, 1986). The algorithms are incorporated in GenStat, Rothamsted's general statistical program, now in its 19th release (Payne, 2018), to provide the facilities for estimating and modelling spatial covariances completely under the control of the practitioner and with transparent monitoring of the processes. These facilities include the choice of steps, bins and maximum lags, the robust estimators and variable weighting according to the expected values. They include the linear model of coregionalization for two or more random variables. GenStat also contains the algorithms to analyse spatial data by reml and so cope with trend and external drift. It has been a boon to those of us at Rothamsted and the scientists elsewhere who have sought our help and collaboration.
And what about the future? I shall not be drawn; it is in the hands of my successors, and Lark et al. (2019) in their commentary have dealt admirably with my sentiments. Theirs is the future.