Analysis of Proximal and Remotely Sensed Spectral Data for the Prediction of Soil




Doctoral Thesis

Soil Science Department and Geography Department of the University of Reading





Benjamin Warr

Supervisors: Dr. M.A. Oliver (Soil Science) and Dr. K.H. White (Geography)



Date of Submission 30/03/2002


Benjamin Warr was born in Zambia in 1972, and has lived in the UK since the age of five. He completed his first degree in Environmental Soil Science in 1997, which he passed with Honours. Studies for his PhD started immediately on winning a three year University of Reading PhD Scholarship Award and subsequently a 2 year Leverhulme Study Abroad Scholarship. Since 1997 Benjamin has been living in France and studying at the Centre de Géostatistique, at the Ecole des Mines, Fontainebleau, the birthplace of geostatistics, where he has completed the Cycle de Formation Specialisée en Géostatistique (C.F.S.G.), a Masters course in the theory and application of geostatistics. Whilst completing his thesis Benjamin has learnt to speak French, worked for Earth Resource Management Services (ERMS) at TOTAL FINA ELF as a geostatistical engineer and for Elsevier Scientific Publications. For the past two years he has been employed at INSEAD in Fontainebleau, where he continues to research subjects of environmental concern.


Contact details:

Benjamin Warr

Centre for the Management of Environmental Resources (CMER)


Boulevard de Constance




Tel: 33 (0) 1 60 72 44 56

Fax: 33 (0) 1 60 74 55 64




Five years have passed since the inception of this thesis and many realities of both my personal and professional life have changed or developed. Over nine years have passed since I went ‘back to school’, to take control of an overactive imagination and desire to understand. This journey did not and could not have started alone and I would like to thank old friends for providing me with the courage to make the leap. Others who have provided continual support include my family, in particular my mother. Also, Dr. Margaret Oliver for all her patience and thoughtful advice. It was Margarets’ idea to send me to France to study with the best and for this I am most grateful. As a result I have built a life in France and a profession that should see me warm and fed. I would like to thank Dr Margaret Armstrong who invited me to visit the center and whose example encouraged me to stay in France. I would also like to thank Dr. Kevin White, Dr. Balwant Singh, Professor Alex McBratney and Dr. Inakwu Odeh, who welcomed me in Australia to do my fieldwork. In particular, thanks go to Odeh for accompanying me while surveying and passing on his experience of the soils of the region. On a personal note, I want to thank the many friends and colleagues who have supported me and this work, because without them this thesis would never be finished and life wouldn’t be the same.


Pedometrics research has been eloquently defined as the ‘development of systems and methods which can provide us with increasingly accurate measurements or estimates of the properties of the soil with improved spatial and temporal coverage’ Wielemaker et al. (2001). The data and methods of analyses currently being developed by the pedometrics community are enabling environmental scientists form many disciplines to develop and test increasingly accurate crop and environmental models at ever-larger spatial scales, and improved spatio-temporal resolution (Addiscott 1993, Curran 1993). The copious information provided by such models is vital to the quest for sustainable development, as is the measured data, in it’s own right, for monitoring purposes. Yet each advance brings with it questions and problems. The three major tasks facing the spatial analyst today are related to measurement, data preparation and model formulation.

Measurements made in the laboratory have become increasingly important for classifying the soil. They represent a distinct improvement upon the ‘ancient’ hand-eye methods because they are quantitative and objective. It is therefore possible to ascertain a measure of the uncertainty inherent to the measurement process and for different people to replicate and validate the results ceteris paribus. Recent advances in measurement technology have provided two contrasting yet similar techniques that provide data about the soil. These include proximal non-destructive and remote observational measurement. The data provided by these sensors is of varied spatial and spectral resolution, accuracy and format. It often requires much manipulation, to first identify those variables with explanatory power and then to prepare these data for inclusion into spatial models. These data are very different from but are related to one another and in turn to the more traditional hand-eye techniques developed by our predecessors and still widely in use and of considerable importance today.

Yet despite these advances, collecting soil samples and measuring soil properties using ‘standard’ practices remains a time consuming and expensive process. The number of samples collected is therefore limited often by financial constraints to the detriment of the possible resolution and accuracy / precision of the results. Therefore quantitative techniques to spatially manipulate sparse sample data have developed concurrently. The quantitative and statistical methods that are currently available offer many possibilities for the inclusion of diverse data types into the spatial model. The questions posed at the inception of this thesis included,


1.           how can the data be prepared for inclusion in spatial models?

2.           How can modern measurements be related to existing or trandtional measures?

3.           Which methods are best suited to incorporate the data into the spatial model to improve the predictions of soil properties?


These three themes echo the major tasks of measurement, data preparation and model formulation, the central themes of this thesis and much pedometrics research.

Table of Contents


Chapter 1. 1

Characterising soil variation. 1

Introduction. 2

The study area and survey data. 7

Methods. 10

Classification. 10

Interpolation. 11

The variogram.. 12

Indicator coding. 13

Ordinary and indicator kriging. 13

Simulation. 16

MultiGaussian Random Functions. 17

Conditioning multiGaussian Random Functions. 18

Simulating indicators. 18

Analyses and results. 21

The conceptual model 21

Quantitative classification. 25

Interpolation and simulation. 27

Structural analysis and interpolation. 27

Conclusions. 35


Chapter 2. 38

Predicting and simulating the geometry of soil horizons. 38

Introduction. 39

Kriging continuous attributes with positive skew.. 40

Kriging categorical attributes to delineate zones of presence. 41

Methods and Analyses. 41

Structural analyses of horizon attributes. 42

Comparing prediction methods. 45

Simulating the in-situ variables to characterize spatial uncertainty. 45

Results. 46

Comparing prediction bias using a simulated dataset 46

In-situ prediction using the ‘optimal’ method. 48

Kriging and simulation measures of uncertainty. 50

Conclusions. 53


Chapter 3. 55

Soil colour and spectral reflectance. 55

Introduction. 56

Theory. 59

Electromagnetic energy interactions with the soil 59

Colour, the perception of EM energy. 63

An appearance basis: the Munsell colour notations. 64

A psychophysical basis of additive primaries: the CIE and RGB colour systems. 65

CIE tristimulus coordinates. 67

CIE trichromacity coordinates. 67

CIE L* a b system.. 68

The RGB system.. 69

Psychophysical colour spaces: Hue, saturation, lightness. 70

Methods and analysis. 71

Spectral reflectance measurements. 71

Results. 74

Comparing colour space representations and the effects of measurement quality. 76

Observed colour: coordinate representations. 77

Observed colour: changes with depth. 78

Spectrally derived colour: Transformations. 81

Spectrally derived colour: Chromatic representations. 81

Correlations between colour space variables and soil physical and chemical properties  82

Improving the fuzzy classification using quantitative colour measurements. 83

Conclusions. 85


Chapter 4. 88

Multivariate spatial predictions of soil colour. 88

Introduction. 89

Theory. 89

Issues when kriging multivariate colour representations. 89

Non-uniformity of scale and the property of additivity. 90

Non-uniformity of scale and non-linear multivariate correlation. 90

Compositional data: Constant sum constraints and linear dependence. 91

Compositional kriging for constrained predictions. 91

Methods. 94

Analyses and results. 95

Structural analysis. 95

Kriging quality. 98

Representation quality. 101

Conclusions. 104


Chapter 5. 106

Remotely sensed observations of the soil 106

Introduction. 107

Sensor systems and platforms. 109

Multispectral optical sensors. 109

Synthetic Aperture Radar (SAR) - Imaging and Applications. 114

Interpreting SAR Imagery. 115

Sensor noise and image statistics. 117

Methods and analyses. 118

Image Model 118

Fourier analysis using the 2-D Fast Fourier Transform (FFT). 118

Factorial Kriging. 120

Kriging the intrinsic component 121

Kriging a second-order stationary component 122

Data Reduction. 122

Study 1: Comparing Fourier analysis and Factorial Kriging analyses. 122

Applying Factorial Kriging. 122

Structural Analysis of Multi-temporal SAR Imagery. 122

Structural Analysis of a SPOT-XS Image. 123

Factorial Kriging of Multi-temporal SAR Imagery. 124

Factorial Kriging of SPOT-XS Imagery. 125

Applying Fourier Analysis. 125

Qualitative comparison of the methods and results. 127

Data reduction and upscaling using Factorial Kriging. 130

Study 2: Speckle filtering SAR image change detection images. 131

SAR Change Detection Images. 131

Analysis and results. 133

Structural analysis and variogram fitting. 133

Ratio Image Filtering. 134

Cleaning Thresholded Ratio Images Results. 138

Study 3: Applying FK to improve soil property – image derived variable correlations. 142

Methods and analysis. 143

Structural Analysis of Landsat TM Imagery. 143

Factorial Kriging of Landsat TM Imagery. 144

Conclusions. 144


Chapter 6. 144

Calibrated soil reflectance models for the prediction of soil properties. 144

Introduction. 144

Study 1: Partial least squares regression. 144

Introduction: Problems with ordinary least squares regression. 144

Collinearity and Over-fitting. 144

The Partial Least Squares (PLS-1) algorithm.. 144

Methods: Calibrating the PLS model 144

Results. 144

Study 2: Multivariate spatial prediction methods. 144

Introduction: maximising the use of available data. 144

Co-kriging. 144

Regression and residual kriging. 144

Fuzzy Classification Mapping. 144

Kriging with an external drift 144

Methods: scenarios of data availability. 144

Analyses and results. 144

Conclusions. 144


Table of Figures

Chapter 1

Figure 1. 1. Schematic representation of the spatial-spectral-coverage trade-offs, when measuring electromagnetic radiation at the Earths’ surface. 4

Figure 1. 2. Bourke Irrigation Area. Origin (bottom left): 37500 E, 6655000 N. 9

Figure 1. 3. Cotton fields after harvest (June), showing the permanent ridge and furrow planting scheme. 20

Figure 1. 4. Natural eucalyptus woodland vegetation typical cover for the Red soil. 20

Figure 1. 5. Typical Red soil surface, showing a small termite hole. Note the sandy texture and intense red colour. 21

Figure 1. 6. Cotton soil surface, showing cracking and accumulation of salts at the surface, evidenced by the dominance of salt resistant scrub plants. 22

Figure 1. 7. A Transition soil surface, cracking like the Cotton soil, but with iron in its oxidized form like the Red soil. 22

Figure 1. 8. Red soil profile (showing depth of approx. 1.50 m), with duplex clay horizon staring at 110 cm depth. 23

Figure 1. 9. Well drained transition soil profile (showing depth of approx. 1.50 m), showing cracking in the upper 1 m. 23

Figure 1. 10. Box and whisker plots showing the distribution of values for the 0-10cm sample for each soil type. (A = Levee, C = Cotton, R = Red, T = Transition). The width of the boxes is proportional to the number of samples allocated to each soil type, the notch indicates the mean, horizontal bars +/- 2 standard deviations, and the dots possible outliers. 24

Figure 1. 11. Fuzzy classification of sample data (a) partition plot and (b) silhouette plot. The symbols in (a) indicate the closest crisp allocation; r Red, ¡ Cotton. 26

Figure 1. 12. Histograms of soil chemical properties: (a) pH, (b) ECe (Sm-1), (c) organic matter (%), (d) Fe (mg / kg) (e) sand, (f) silt and (g) clay content, (h) diffuse, (i) nodular and (j) gypsum horizon thickness. 28

Figure 1. 13. Variograms of soil chemical properties: (a) pH, (b) ECe (Sm 1), (c + g) organic matter (%), (d) Fe (mg / kg) (e) sand, (f) silt and clay (g) content (units = %), (h) fuzzy class members 1 and 2 (which have identical variograms). 30

Figure 1. 14. Ordinary kriged (OK) map (left column) and sequential gaussian simulation set (n = 500) standard deviation (right column), for pH (a + b), ECe (c + d), organic matter (e + f). 32

Figure 1. 15. Ordinary kriged (OK) map (left side) and sequential Gaussian simulation set (n = 500) standard deviation (right side) for iron content (a + b) respectively for Red soil member 1 (c + d) and Cotton soil member 2 (e + f). 33

Figure 1. 16. Ordinary kriged (OK) map for (a) sand (b) silt and (c) clay content. 34

Figure 1. 17. Maps of the OK estimation standard error for (a) FeOK and (b) pH*OK, illustrating the inutility of expressing variable specific spatial uncertainty using the estimation variance, when the sample locations are sparse and irregularly spaced. 35


Chapter 2

Figure 2. 1. Experimental (dots) and modelled (thick line) variograms for (a) diffuse carbonate indicator, (b) diffuse carbonate thickness, (c) diffuse carbonate thickness normal scores, (d) nodular carbonate indicator, (e) nodular carbonate thickness, (f) nodular carbonate thickness normal scores, (g) gypsum indicator, (h) gypsum thickness, (i) gypsum thickness normal scores. Simulating a dataset with realistic characteristics. 43

Figure 2. 2. Indicator variograms for selected cut-offs applied to diffuse horizon thickness, zk (a) z1 = 0.06 (b) z2 = 0.15 (c) z3 = 0.32 (d) z4 = 0.51 (e) z5 = 0.67 (f) z6 = 0.83 (g) z7 = 1.13 (h) z8 = 1.27 (i) z9 = 1.36. 44

Figure 2. 3. Simulated transect thickness estimate error histograms, (a) ZSIM – Z*SIM OK, (b) ZSIM – Z*SIM IK, (c) ZSIM – Z*SIM MG, (d) ZSIM – Z*SIM OPT. 47

Figure 2. 4. Simulated transect, (a) true simulated thickness, TSIM, (b) T*SIM OK, (c) T*SIM MG, (d) T*SIM IK, (e) T*SIM OPT. The points represent the position of the data used to krige. 48

Figure 2. 5. (a ‑ c) Horizon thickness estimate (units: meters) (d ‑ f) MG kriging estimation variance and (g ‑ i) sGs thickness E-type mean for each horizon type. 50

Figure 2. 6. (a-c) sGs indicator simulation standard deviation, (d-f) sis simulation standard deviation for each horizon type. 51

Figure 2. 7. Scatterplots of simulation post-processing results for the diffuse carbonate horizon – (a) sis mean thickness vs. sis standard deviation, (b) sGs mean thickness vs. sGs standard deviation, (c) sis mean vs. sis standard deviation, (d) sGs indicator coded mean vs. sGs indicator coded standard deviation. 52


Chapter 3

Figure 3. 1. Diagrammatic representation of an electric and magnetic field. 59

Figure 3. 2. A division of the electromagnetic spectrum into regions, showing their frequency in Hertz and potential for proximal (and remote) sensing of various soil properties, including potassium (K), soil colour (SC), soil organic matter, clay content including iron (CC), soil nitrogen (SN), cation exchange capacity (CEC), soil water (SW), and soil structure – soil strength measurements (SS). (Adapted from Viscarro-Rossel & McBratney 1998.) 60

Figure 3. 3. Possible interactions of electromagnetic energy with a material. 61

Figure 3. 4. A possible reflectance spectra for the visible part of the electromagnetic spectrum. 63

Figure 3. 5. (a) Munsell Hue (b) Munsell Value,(c) Munsell Chroma (colour wheel) and (d) the complete Munsell colour system, showing the non-uniformity of the Chroma axis. 64

Figure 3. 6. Additive colour combinations within the RGB colour space. Illustrates how different real RGB primaries generate different colours when combined in the same proportions. 65

Figure 3. 7. (a) CIE colour matching functions and (b) spectral power distribution for the Standard Illuminant C. Required for transformation of reflectance curves into colour coordinates. 66

Figure 3. 8. (a) Mapping of XYZ tristimulus values to (b) an x y Chromaticity Diagram. 68

Figure 3. 9. The CIE xy Chromaticity diagram (a) and the CIELUV Chromaticity diagram (b). Each line represents a colour difference of equal proportion. 68

Figure 3. 10. (a) A simple 3-D representation of the CIELAB system and (b) a 2-D representation of a CIELAB colour plane of equal Luminance. 69

Figure 3. 11. (a) the RGB gamut in 3-D and (b) a 2 dimensional RGB gamut showing that linear combinations of the primaries produces the range of possible colours. 70

Figure 3. 12. Schematic diagram of a double beam spectrometer. 71

Figure 3. 13. Examples of spectra with median loading on factors 1, typical of a Red soil (a) and factor 2 , typical of a Cotton soil (b) [solid black line], with interquartile range [dashed line]. 74

Figure 3. 14. Screeplots (a) number of wavelengths = 1065, (b) number of wavelengths = 252 equivalent to AVIRIS sampling interval and FWHM. 75

Figure 3. 15. Samples plotted in the space of component scores (a) components 1 and 2, (b) components 2 and 3, both show a partition of according to soil type (as indicated by letters: C = cotton soil, R = red, T = transition and A = levee), and (c) components 3 and 4, that does not. 75

Figure 3. 16. Eigenvectors of the first four principal components of the reflectance curves (a-d) number of wavelengths = 1065, (e-f) number of wavelengths = 252 equivalent to AVIRIS sampling interval and FWHM. 76

Figure 3. 17. Scatterplots of Munsell Hue, Value and Chroma and CIE L* a b Hue, Luminance and Chroma analogs, showing the 1:1 relation between the colour space representations. 77

Figure 3. 18. Box and whisker plots showing the changes in observed Munsell Hue, Value and Chroma with depth for 4 soil layers (0-20 cm,.30-40 cm, 70-80 cm and 90-100 cm). The width of the boxes is proportional to the number of observations in each soil type (A = Levée, C = Cotton, R = Red and T = Transition soil types) 78

Figure 3. 19. Box and whisker plots showing the changes in observed red, green and blue coordinates with depth for 4 soil layers (0-20 cm,.30-40 cm, 70-80 cm and 90-100 cm). The width of the boxes is proportional to the number of observations in each soil type (A = Levée, C = Cotton, R = Red and T = Transition soil types). 79

Figure 3. 20. A comparison of colour coordinates derived from observations made in the field and those calculated from spectral reflectance measured under laboratory conditions. 80

Figure 3. 21. MLR predicted values and simultaneous confidence intervals for (a) sand content (%), (b) total iron content (mg / kg), and (c) pH. Soil types are indicated by the shape and shade of the symbols (square=cotton, triangle=transition, round=red). 83

Figure 3. 22. Correlations between fuzzy class memberships from spectral colour and laboratory data (y-axis) and only laboratory data (x-axis), with soil types allocations designated by shade and shape (square=cotton, triangle=transition, round=red). 83

Figure 3. 23. Clusterplots and silhouette plots for two fuzzy classifications performed using (a + b) selected chemical and physical analyses and observed soil colour, and (c + d) the same data with spectrally derived soil colour. 84


Chapter 4

Figure 4. 1. Plot of the spectrally derived xy colour coordinates, showing the soil type allocation. Note that the origin of the chromaticity diagram (a) corresponds to the chromatic point of the standard Illuminant C. The arrows indicate the dominant wavelength of selected samples, and the horizontal and vertical lines possible cut-offs for a rapid colour based classification of the samples into the two modal soil types. 81

Figure 4. 2. Variograms of (a-c) Munsell hue, value and chroma, (d-f) CIE L* a b hue, luminance and chroma, (g-i) red, green blue. Variables calculated from in-situ observation and Munsell coding. 96

Figure 4. 3. Variograms of (a-c) Munsell hue, value and chroma residual from linear regression estimate, (d-f) CIE L* a b hue, luminance and chroma, (g-i) red, green blue. Variables calculated from reflectance spectra. 97

Figure 4. 4. Predicted values vs. sample values of the observed colour values for (a-c) Munsell Hue, Value and Chroma, (d-f) CIE L* a b Hue, Luminance and Chroma both co-kriged and (g-h) rgb trichromaticity coordinates by compositional kriging. 98

Figure 4. 5. Predicted values vs. sample values of the spectrally derived colour values for (a-c) Munsell Hue, Value and Chroma, (d-f) CIE L* A b Hue, Luminance and Chroma both co-kriged and (g-h) rgb trichromaticity coordinates by compositional kriging. 99

Figure 4. 6. CIE L* a b estimates vs. Munsell estimates (a) Hue, (b) Value and (c) Chroma  100

Figure 4. 7. (a-c) Cmp K estimates of xyz and (d-f) the backtransformed XYZ values plotted vs. sample values. 100

Figure 4. 8. Kriged maps of observed colour: CK estimates of Munsell Hue, Value and Chroma (a-c), and CIE L* a b Hue, luminance and Chroma (d-f). 102

Figure 4. 9. Kriged maps of spectrally derived colour: CK Munsell Hue, Value and Chroma (a-c), and CIE L* a b Hue, luminance and Chroma (d-f). 102

Figure 4. 10. False colour composites of soil colour created using the observation derived rgb estimates by (a) cokriging and (b) compositional kriging spectrally derived rgb estimates by (c) cokriging and (d) compositional kriging. 104


Chapter 5

Figure 5. 1. Diagram showing relative atmospheric radiation transmission of different wavelengths. 109

Figure 5. 2. Landsat TM image of the Bourke Study area. The region covers 900km2 (30km by 30km). Reference system details: UTM-55S, Units: meters, Top-Left Corner: 375050, 66859500, Top-Right Corner: 405050, 6685950, Bottom-Right Corner: 405050, 6655950, Bottom-Left Corner: 375050, 6655950. 111

Figure 5. 3. Hyperspectral imaging capabilities of AVIRIS. 113

Figure 5. 4. Transmission and receipt of scattered signal. 114

Figure 5. 5. ERS SAR image of the Bourke study area. 115

Figure 5. 6. Landscape effects on mean backscatter and texture. 116

Figure 5. 7. Chi-Square probability density function with 2n degrees of freedom, (exponential distribution), typical of high and low SAR microwave backscatter regions . High backscatter region (mean = 200), low (mean= 50). 117

Figure 5. 8. Space and frequency domain images (256 by 256 pixels): (a) original image, and the corresponding frequency domain image (b), edited frequency domain image (u, v coordinates) (c) and the original image after FFT filtering (d). 119

Figure 5. 9. (a) SAR nested variogram model, (b) component models. 120

Figure 5. 10. SPOT image (a) nested variogram model, (b) component structures. 124

Figure 5. 11. A comparison of target and retrieved components of the SAR image obtained using a 5 by 5 and a 7 by 7 pixel diameter moving window. 124

Figure 5. 12. A comparison of target and retrieved (a) short-range and (b) long-range components of the Landsat-TM image band. 125

Figure 5. 13. Fourier raw and filtered magnitude images, (a) original SAR FFT magnitude image, and (b) filtered SAR magnitude image (low pass filter), all 256 by 256 in extent. 126

Figure 5. 14. (a) original SPOT FFT magnitude image, and (b) high pass filtered SPOT magnitude image (c) low pass filtered SPOT image, all 512 by 512 pixels in extent. 126

Figure 5. 15. Annotated landcover map of the region, 400 km2, derived from d’Herbes & Valentin, 1997. Colour version with legend available from 127

Figure 5. 16. ERS-1 SAR 25 km2 August 1992 subscene (a) Raw image, (b) FK long-range component image (components Y3 + Y4 estimated). (c) FFT long-range component image (high pass filter applied) and (d – f) SPOT XS 2 June 1992 400 km2 subscene (d) Raw image, (e) FK long-range component image, (f) FFT long-range component image. 128

Figure 5. 17. SPOT XS 2 June 1992 400 km2 subscene (a) FK short-range component image (b) FFT short-range component image. 129

Figure 5. 18. FK filtered and upscaled long-range component SPOT XS-2 image with a pixel size of 100 m by 100 m   129

Figure 5. 19. July Long-range Components - 20, 50, 100m Grid. 131

Figure 5. 20. Ratio Image Variogram Models and Filtered Image Models. 134

Figure 5. 21. Filtered Ratio Image q-q Plots (black), 1:1 bisector (grey ). 135

Figure 5. 22. July-August Ratio Images. 136

Figure 5. 23. Effects of Connected Components Removal 139

Figure 5. 24. June to August: Thresholded at +/- 3.5dB: Long-range Component Images Filtered Using Different Morphological Filters. 140

Figure 5. 25. June to August: Thresholded at +/- 6.5 dB Long-range Component Images Filtered Using Different Morphological Filters. 141

Figure 5. 23. LANDSAT TM images of the Bourke study area (a) NDVI (masked NDVI<0), (b) band 1 raw image (masked NDVI<0). 142

Figure 5. 24. Landsat TM anisotropic variogram models for (a) band 1, (b) band 2. 143

Figure 5. 25. LANDSAT TM images of the Bourke study area (a) FK denoised band 1 (b) FK LRC band 1 (c) band 2 FK denoised, (d) band 2 FK long-range components (LRC), (e) PC 3 of FK denoised bands 1-7, (f) PC 3 of FK LRC bands 1-7. 144


Chapter 6

Figure 6. 1. (a) Sum of square of the Residual matrix Fk indicating the presence of an outlier and (b) the row sum of squares of the Residual Ek  showing the increasing amount of information of explanatory importance extracted from the spectral data with the inclusion of additional PLS factors. Each horizontal trace represents the row sum on inclusion of another factor from top to bottom (example variable – pH). 144

Figure 6. 2. PLS-1 regression coefficients for optimal model (red dashed line) for (a) iron content, (b) sand), (c) ECe (d) organic carbon (e) pH and (f) clay content. Also plotted are the average regression coefficients (solid black line) and + / - 2 standard deviations, calculated from 100 from models developed by jack-knifing using randomly selected training sets (n = 88) and a test sets (n = 30). 144

Figure 6. 3. PLS-1 weighting vectors 1 to 4. 144

Figure 6. 4. PLS-1 component scores, T one and two, of the optimal PLS-1 model for each variable (a) iron content, (b) sand), (c) ECe (d) organic carbon (e) pH and (f) clay content. The shade and shape of the points correspond to the soil types (square: cotton, triangle (up): red, triangle (down): transition). 144

Figure 6. 5. Correlation between true and predicted values for (a) iron, (b) sand, (c) ECe, (d) organic carbon (OC %) , (e) pH and (f) clay. Predictions correspond to the PLS-1 model developed using training set n = 88 and a test set (n = 30). 144

Figure 6. 6. Variograms of PLS-1 estimate error residuals for (a) iron content (mg/kg), (b) sand (%), (c) ECe (Sm-1), (d) organic carbon content (%), (e) pH and (f) clay content. 144

Figure 6. 7. Maps of predicted soil properties. OK (a) iron (mg/kg) and (b) sand content (%), PLS-B2 estimates for (c) iron (mg/kg) and (d) sand content (%). 144


Chapter 1

Characterising soil variation


To achieve environmental sustainability today’s land management interventions must remain effective far into the future. The past fifty years of environmental research has provided modelling tools that can help evaluate the potential behaviour of the soil under different management scenarios. Models that couple the behaviour of the soil, vegetation and atmosphere are no longer novelties and are capable of predicting soil processes accurately for many management purposes when sufficient data is available. Simulation models have been incorporated into systems to support the decision-making process and show considerable potential to help define management strategies that are robust to climatic and economic uncertainty and soil variability (Addiscott, 1993). Yet, the widespread uptake of agricultural models as management tools has been poor (Keating and McCown, 2001). The potential of modelling experiments to expedite landscape management has not been fully developed because of a lack of spatial data to parametize and validate the models. A lack of accurate data to describe the spatial and temporal variations of soil properties is perhaps the primary factor limiting the intensity and quality of human interventions in land management systems. Land managers fortunate to be the focus of academic research have benefited from the learning process that includes stages of data collection, model development, simulation and importantly validation. Those not involved are likely to remain unconvinced of their utility. Cooperation during model implementation and validation are the most effective means to convince a land or production manager of the potential benefits they may provide. The successful diffusion of such practices requires soil survey and sample measurement methods capable of providing data accurate enough to structure the soil model with a precision and a spatial resolution fit for the purposes envisaged. The need for good quality soil data cannot be avoided. Misperceptions of the structure and dynamics of a system are the cause of inaccurate simulation results and consequently the reason for poor management decisions (Moxnes, 2000). The lack of a cost-effective means of providing necessary soils data can also mean that models are poorly parametized. The simulation results may suffer. In contrast to climatic, crop or management variables, which vary continuously during simulation, the future behaviour of the soil and consequently all related components of the system are conditioned by the initial data used to parametize the model. An accurate representation of the soil is therefore an essential component of a successful modelling exercise. If the results are poor land managers will remain unconvinced and the diffusion of such beneficial practices will be slowed and/or reduced. The impacts of today’s management practices on future productivity will not be evaluated and almost inevitably the quality of land management practices will suffer. The expense of collecting sufficient site-specific or spatial data lies at the heart of the problem of delivering the services that are needed to increase land productivity sustainably.

Efforts to provide cost-effective solutions to map the properties of the soil are the focus of this thesis. In the following chapters methods of measuring and predicting soil properties will be described and tested. Soil data has traditionally included only qualitative visual observations and expensive time-consuming chemical and physical analyses performed on samples transported to the laboratory. The advent of modern measurement technologies and analytical methods has meant that new data are available at little extra cost. Soil survey data can now commonly include traditional and modern data types for example,

The cost of collection, the information content and format of the data provided by each measurement method differs significantly. Measurements made by observation are invaluable to provide a first conceptual model of soil spatial variability during the initial stages of soil survey. Physical and chemical analyses are then necessary to validate the conceptual model and provide sample averages for these properties for the region and by soil type. This is a ‘traditional’ approach to soil survey and the spatial prediction of soil properties that suffers from several drawbacks:

·                       visual observations are subjective measurements and are poorly reproducible,

·                       soil survey and traditional laboratory analyses methods are time consuming and expensive,

·                       traditional chemical and physical analyses are destructive

·                       the sample density and coverage are limited by cost-time restraints,

Modern measurement technologies and statistical analyses offer a means of overcoming some of these problems. The most widely used of these technologies exploit the properties of electromagnetic (EM) radiation and its characteristic interactions with those of soil, permitting measurements of the quantities of reflected energy to be used as diagnostic of soil physical and chemical properties. The same characteristic interactions allow an observer to distinguish between soils. When colour matching the human eye is measuring the intensity of light energy arriving at the eye from distant objects, and then processing these stimuli with the brain. Observation is a ‘remote’ and ‘non-destructive’ measurement technique. Proximal sensing (PS) and remote sensing (RS) technologies also provide non-destructive measurements of distant objects by measuring the EM radiation coming from the object. Sensors are capable of measuring EM radiation with precision, over a much larger area and range of wavelengths than an observer can. But at present there is a three-way trade-off between these advantages (Figure 1. 1). The spectral, spatial and temporal resolution of different sensor systems are fixed system parameters that have to date limited the utility of ‘operationally’ available remotely sensed data. Exhaustive and large area spatial coverage landsurface measurements are available from satellite remote sensing devices such as the Landsat Thematic Mapper. But the poor spatial and spectral resolution of much of the commercially available optical remotely sensed data to date, are drawbacks that have implications for nearly all subsequent analyses.

The processes of isolating soil related information from such imagery are fraught with problems. Regardless of the precision of the sensor, vegetation often hides the soil completely in humid environments. While it is possible to interpret the presence or absence of a certain vegetation community with a given soil type it is not desirable to do so. The observations are not direct and there are no reasons why a given vegetation community should exist exclusively on a precise soil type. For this reason many of the applications of RS data to soil are restricted to arid and semi-arid environments (Escafadel et al. 1993), or to certain periods of the year when the soil is bare. Where the soil is entirely covered throughout the year remote sensing technologies provide only partial or surrogate information. Other concerns raised by the use of remotely sensed imagery include the possible error introduced by georectification and positional inaccuracies. Atmospheric absorption may hide important soil information carried by specific bands within these absorption regions. Poor temporal coverage may also render the data useless for the proposed purposes.

Figure 1. 1. Schematic representation of the spatial-spectral-coverage trade-offs, when measuring electromagnetic radiation at the Earths’ surface.

Airborne hyperspectral data appears to offer a potential solution to many of these problems. Capable of being flown at any altitude they can provide high-resolution data that is only weakly affected by atmospheric absorption. Measurements can also be made with a temporal frequency dictated by the objectives of the study, or to collect data at specific times when the soil is bare. Still, the measurement provided by these sensors may not be ‘optimal’, because the sensors were simply not designed specifically for soil studies. Georectification problems can be exacerbated by the movements of the aircraft and the proximity of the scene.

The spatial complexity of the images and the simple fact that there are so many individual measurements across many wavelengths may exacerbate the process of isolating soil related information from features of little interest, for example manmade features or simply shadow.

There are very few studies that have successfully used remotely sensed data of any type to directly predict the value of continuous soil properties (Palacios-Orueta et al 1999).

The use of RS data as soft ancillary information is more common (Odeh & McBratney, 2000). The majority of studies using such imagery use either classification or spectral un-mixing methods to provide abstract measures that relate to the landsurface as a whole, such as the percentage soil and vegetation cover or a degree of class similarity. The complexity of the landsurface and contamination of the soil reflectance signal by more dominant reflectors are such that remotely sensed imagery has perhaps been less useful than originally thought initially possible. When modelling at a coarse scale, soil parameter values identified from surrogate and often abstract information may suffice. However more detailed studies will require spatially distributed continuous parameters to drive soil processes. For this robust calibration and regression methods are required. All of these issues, which are discussed in greater detail throughout this thesis, suggest that the collection of ever-greater quantities of in-situ data cannot be avoided and is indeed necessary. A cheaper alternative to standard chemical and physical analyses is required. Spectral reflectance measurements made directly of the soil / landsurface or on samples in the laboratory offer this potential (Escafadel et al, 1993).

A compromise, which combines the cheap and rapid measurement processes offered by modern laboratory and remote sensor technologies with more standard chemical and physical analyses might represent the cost effective solution. Interpolation methods exist to spatialize the useful information derived from the spectral reflectance data measured at visited locations, using soft exhaustive data to improve the spatial estimates. Remotely sensed data is likely to have a greater correlation with data of a similar nature (electromagnetic) and this virtue should provide a basis for the interpretation of the relationships between the sensed signal and the soil in physical terms.

The progress in measurement technologies permits a flexibility that can provide solutions to certain of the issues raised above, but significant possibilities and questions remain to be explored and answered (see Box 1). The rest of this thesis is devoted to identifying solutions for the inclusion of measurements of the reflectance properties of the soil, into models for the prediction of soil properties.

Spectral reflectance measurements made on soil samples will effectively remove any confusing vegetation effects, but how should this multivariate data be used for soil property prediction; what are the alternatives to classification? What methods exist to reduce the dimensionality of this data while keeping the important soil information?

Data reduction is required for various reasons discussed later, not least, because the sensor systems where not specifically designed for soil studies and many of the measurements made may be redundant. Can the process of data reduction lead to the identification of simpler, cheaper measurements or representations of the soil that are better suited to soil studies or with which soil scientists are more familiar? How can the predictive potential of remotely sensed data be realised within this framework? What improvements in prediction quality can measurement technologies provide? Progressively four types of data including, qualitative observations of horizonation and colour, spectral reflectance measurements collected in the laboratory, colour measurements derived from these and remotely sensed data from operational satellites were explored for their information content and ability to predict soil properties. In the first two chapters standard methods were applied using traditional data types (A) to provide (non)-spatial predictions of soil properties. In the optimal situation of data availability explored in the final chapter laboratory spectral reflectance measurements and remotely sensed data augmented traditional chemical and physical analyses and observations made on samples collected in the field. The main body of this thesis involved the search for methods to overcome calibration problems posed by the high-dimensionality of the diffuse spectral reflectance measurements to enable their inclusion in (spatial) regression models. As well as testing multivariate calibration methods for hyperspectral data, methods were sought to reduce the dimensionality of reflectance data while retaining useful information content. Soil colour has always been an important diagnostic attribute and one of the most important characteristics used to differentiate soil in most systems of classification and soil survey (Stoner & Baumgardner 1980, Kingham 1998). The historical use of soil colour, as a diagnostic and predictive tool and its future role as a means of communicating reflectance information, are the theme of an entire chapter. Colour measurements were also tested as an alternative to hyperspectral data, reflecting a situation where the high-dimension spectrometer data is replaced by a low-dimension colour representation. Colour and remotely sensed data were also tested for their utility as a means to augment the amount of spatial information available and increase the precision of the final estimates.

In the rest of this first chapter the study area and traditional soil measurements, including qualitative soil observations and ‘standard’ destructive chemical and physical analyses are presented and the results of what could be considered ‘standard’ soil survey methods evaluated. While the end objectives of modern soil survey should inevitably lead to quantitative predictions of continuous soil properties the benefits to be gained from the process of observation and conceptualisation should not be disregarded as they indicate the correct framework and methodologies with which to take the analyses further.

The study area and survey data

Aridity was perhaps the most important criteria when identifying likely places to carry out this research. Soils under a low rainfall regime are less densely covered by vegetation, and therefore the remotely sensed signal is less contaminated by vegetation reflectance characteristics. A study site in southwest Niger, the focus of the Hapex-Sahel experiments (Prince et al. 1995) was the first location chosen for this reason and because of the abundance of freely available remotely sensed imagery. However it was not possible to complete the necessary fieldwork there.

The fieldwork for the main body of this study was completed in southwest Australia. The Bourke region is approximately 700 km north-west of Sydney, New South Wales Australia. Bourke is located in the Eromanga Basin, a Jurassic to Cretaceous sedimentary basin that is part of the Great Artesian Basin of New South Wales. The Eromanga Basin is characterised by large expanses of low relief and extensive alluvial systems. The varied depositional environments of the Cainozoic era provided the sediments creating dunes and levées that have since undergone varied pedogenesis (Senior et al. 1978, Kingham 1998). The natural vegetation cover includes eucalyptus forest with a dense shrub understory or tall grass, which in many areas is extensively grazed by sheep. Cleared unimproved grassland extends onto the plains, which in recent years have seen the development of isolated agricultural enterprises. Cotton is the principle crop grown on the plains. In some areas citrus fruit and vineyards have been recently established.

The soil survey that was completed during the Australian winter of 1999 covered an area of 900km2 and lasted one week. The number of samples collected and the density of the sampling scheme were designed to reflect a standard reconnaissance level soil survey. One hundred and fifty sample locations were selected using a stratified random sampling scheme (Figure 1. 2). Of the 150 locations chosen, 110 were visited and at each location at least one soil sample was collected.

At each location visited a series of qualitative observations were made and recorded. These included vegetation structure (forest, scrub, grassland, crop) and moisture and erosion status and soil colour. At 20 locations a second sample was collected approximately 500m away from the first in a direction defined at random. Each sample was bulked from four smaller cores of the uppermost soil layer (0-10 cm) and collected over an area 20m by 20m. This area is approximately equivalent to the landsurface represented by a single Landsat image pixel. If the measurements have a similar support the correlations between them can be considered valid and is likely to be more significant.

A total of 130 surface soil samples were collected of which 110 underwent a series of standard chemical and physical analyses (see Box 2, overpage). Iron content was measured by extraction in 1M hydrochoric acid extraction and organic matter / carbon by the Walkley Black method at a commercial laboratory. The acidity/alkalinity (pH) and electrical conductivity (ECe) were measured by electrode techniques in a 1:5 soil: deionised water solution.

Soil texture was measured using the ‘micro-pippette method’, using only very small quantities of soil, (see Appendix for details). At 100 of the locations it was possible to drill soil cores to a depth of 1 m, from which the presence and thickness of diffuse and nodular carbonate horizons were identified using 0.1 M  HCL. The solution was poured the length of the core and any effervescence noted. Gypsum accumulations, recognised as occupying greater than 10 % of the soil matrix, were identified using the method outlined by Verheye & Boyadgiev (1997).

Reliant on visual identification alone it is likely that in some cores accumulations were not observed where they actually existed. These techniques are not capable of identifying with great detail very thin horizons. There is a detection limit below which the horizon thickness has to be considered zero (~2-3cm). Horizon thickness data is a censored variable. From the same cores soil colour measurements were made of horizons are at four depths (0-20cm, 30-40cm, 70-80cm and 90-100cm) by comparing moistened soil with Munsell colour chips (colour matching). Details of the proximal and remotely sensed measurements and data are provided in subsequent chapters. The specific combinations of the various types of data used for each analysis as well as the validation methodology will be outlined in the methods description specific to each analysis.

Figure 1. 2. Bourke Irrigation Area. Origin (bottom left): 37500 E, 6655000 N.


The rest of this introductory chapter describes the theoretical basis and application of ‘standard’ non-spatial and spatial regression methods for the characterisation and prediction of soil variability. The specific objectives were to

1.                  Describe and test the validity of a conceptual soil classification

2.                  Compare the basis of this classification with its quantitative counterpart

3.                  Provide estimates and maps of important soil properties

4.                  Investigate the spatial variability of the soil properties


Traditionally soil survey and the resulting stratification of the soil into more or less homogeneous units have been based solely on expert knowledge that informed an inductive logic. The framework most often adopted during the development of conceptual soil-landscape models, was developed by Jenny (1941). It states that the soil is a function of climate, parent material, topography, biota. Each of these conditioning factors varies more or less continuously over the surface of the Earth and so the properties of the soil vary in response.

A conceptual model of soil spatial variability developed in the field is based on subjective interpretations of observations, and represents an informal rule based decision-making process. These models represent the first step in formalising the available knowledge about the soils of a region. Often the concepts underlying the model and the methodology required to allocate a soil using these concepts are hard to communicate to others. The model is personal to the soil surveyor. Diagnostic classification systems, for example Soil Taxonomy (Soil Survey Staff, 1975) formalise the rules of soil classification and standardise the description of soil variability, on the basis of the values of diagnostic soil properties. Formal systems are often criticised for being unwieldy, and providing an unrealistic hierarchical representation of soil variability (Webster, 1968). Considerable amounts of soil data, collected at all depths in the profile can be required to correctly allocate a soil to a class. Often this data is unavailable, and under such circumstances informal conceptual models or clustering methods must be adopted.

Quantitative methods of soil classification provide a robust alternative to conceptual models, a means of consistently and quantitatively describing the similarities and dissimilarities between samples using any number of variables that describe as many properties that can be measured (Oliver & Webster, 1990). Individuals can be represented in ‘n’-dimensional character space using various measures to describe the distance between individuals and groups (Gordon, 1987). Many iterative clustering algorithms exist to allocate samples to groups by seeking to minimise the within group dissimilarity of these measures.

Yet, by allocating samples to crisp groups numerical cluster analysis and formal diagnostic classification systems fail to represent the continuous nature of soil spatial variability. Conventional methods assign each location to a single group, regardless of the degree of similarity it shares with locations in other groups. They offer no means of representing the degree of (dis)similarity between two soils. The mathematics of fuzzy sets has provided some useful mathematical tools that correspond closely to our ‘human’ perception of difference (similarity). Zadeh (1965) is responsible for developing the concepts and mathematics behind fuzzy set theory that permits a continuous representation of the ‘degree of similarity’ between samples. Fuzzy classification (FC) methods are valued for their ability to express

a)      generality – when a single concept applies to a variety of situations

b)      ambiguity – a single concept embraces several subconcepts and,

c)      vagueness – there are no precise boundaries

The FC algorithm minimises the objective function J(M, C),

Equation 1. 1         

where c is the number of classes, n is the number of samples, mij is the membership of an individual i in class j and d2 is a measure of dissimilarity (de Gruijter & McBratney, 1988). The main difference between Equation 1. 1 and a standard clustering objective function is the fuzziness exponent j (1 < j < µ). When crisp classifying the value of j is fixed at 1. For FC the value is usually fixed within the range 1.5 – 2.

Burrough et al. (1989, 1992) and McBratney et al. (1992) originally applied fuzzy classification to soil data. It has since proved popular because of its ability to represent soils within continuous character space that reflects the spatial characteristics and our conceptual understanding of how soils vary spatially. Burrough (1989) wrote that “too much flexibility leads to anarchy and too much rigidity causes conflict.” Fuzzy set theory has provided quantitative methods that enable the expression of both flexibility and rigidity and the human concept of vagueness. Where conventional classifications either exclude or include a soil from each class by associating class boundaries with a membership cut-off value, a fuzzy classification describes the relationship between a sample and all classes on a scale of zero to one (zero = no similarity or shared characteristics). The FC algorithm has since been modified (de Gruijter & Mc Bratney, 1988) to distinguish modal soils, from soils that share properties with many groups (intergrades) and those soils that bear no similarities to any of the class centroids (extragrades).


Classification methods are useful when there are measurements of many soil attributes, and this information must be synthesised into a format, which facilitates its interpretation and communication. Increasingly, methods are required to predict the value of a single soil property at unvisited locations, instead of providing a more general description of what the soil is like everywhere.

The objective of soil mapping is to characterise this variability and to provide spatial predictions at unvisited locations. To do this a model of spatial variability is required. A general model of spatial variability of a soil attribute Z assumes some form of structural component m requiring a functional representation, and an error component e,

Equation 1. 2          .

This model is very flexible. Traditionally the functional form of Z(x) was purely conceptual, difficult to communicate and irreproducible. In recent years mathematical and statistical methods have been applied to the problem and various quantitative models have been developed. They formalise a conceptual understanding of how soil varies spatially and provide a toolbox to characterise soil variability and provide predictions at unvisited locations.

During soil survey, empirical and conceptual soil-landscape models are developed to provide predictive information about the spatial distribution of soil classes Dk, over the study domain D. Laboratory analyses provide statistics that describe the average attribute values zk(x) for the class, k = 1 ,…, K present at each location x. The mean m(x), variance s(x) and prediction error e(x) are assumed to be constant within each class. This model is called the discrete (or categorical) model of spatial variation (Heuvelink, 1993) and it is often used to represent the results of a crisp classification. It is the format most often used with Geographical Information Systems, GIS (Burrough, 1989a) because it is simple to associate class statistics with maps in chloropleth or vectorised form. Yet, the assumptions that are made about the structure of the soil spatial variability are unrealistic; soil spatial variation is continuous (Webster, 1968). The categorical model of spatial variation is appropriate for the representation of discontinuous spatial variation, but is often adopted due to a lack of soil observations and sample data (Heuvelink, 1993). The way in which we represent the soil is a compromise between reality and the availability of data that provides useful and reliable information.

The traditional methods of extrapolation or ‘free survey’ that used observations of topography, vegetation and geology to delineate regions of homogeneous soils have largely been superseded by quantitative interpolation methods, for the same reason that quantitative methods now dominate the realm of classification. Tessellation, inverse distance measures, nearest neighbour interpolation, spline functions and trend surfaces have all been applied to the problem (Burrough, 1989). Each method has its drawbacks, for example, local interpolators use data inefficiently and global trend surface techniques fail to predict well at small scales. None of these methods uses the values of the data themselves to determine the spatial weighting and none can provide any measure of the estimation uncertainty (Webster & Oliver, 1990).

Geostatistics, the practical application of Regionalised Variable Theory (Matheron, 1965) developed as a response to these failings and has since evolved to provide a flexible toolbox of statistical methods to explore and predict the spatial, temporal and multivariate structure of data. The Theory of REV’s provides a continuous model of spatial representation. The data z at locations x made over a domain and describing a soil attribute are regionalised random variables z(x) which, within a probabilistic framework, are considered to be just one realisation of a random function Z(x). A random function is an infinite family of random variables at each location over the region of interest. Using the data values measured on sparse samples z(x), geostatistical modelling seeks to characterise the important parameters of the random function Z(x), the structural component of the model of soil spatial variability and e(x) the model error (Equation 1. 2).

The variogram

Geostatistics uses the variogram as a tool to describe the spatial pattern of an attribute and subsequently for prediction and simulation. The experimental variogram of a continuous variable, g*(h), is calculated as half the average squared difference between data pairs,

Equation 1. 3         

where N(h) is the number of data pairs within a given class of distance and direction. The experimental variogram is then modelled using an authorised variogram function (Journel & Huijbregts, 1978). The accuracy with which the variogram parameters are estimated, and hence the variable predicted, may be improved by the reduction of the influence of a small number of large values (Goovaerts, 1997). The spread of the sample histogram can be reduced by transforming the original z-data into normal score y-values,

Equation 1. 4         

where f  is a normal score transform that relates the p-quantiles zp and yp of the two distributions (Rivoirard, 1994). The variogram of the (xa) can then be recalculated using Equation 1. 3. Typically the y(xa) variogram model gY(h), has a smaller relative nugget effect and greater spatial continuity over short lag distances than its gZ(h) equivalent, because of the reduction of influence of extreme values.

Indicator coding

Indicator variables, I(xa ;zk), are generated by applying an indicator transform to a continuous variable, z(xa), at selected thresholds, zk which are often chosen to be equal to the decile values of the cumulative distribution,

Equation 1. 5         

The indicator-transformed values describe the spatial distribution of sets with random location and shape. Therefore the lateral extent of a soil horizon can be described by a categorical indicator variable defined at the threshold zk = 0. The result is a categorical indicator variable equal to one if the horizon sk, is present and zero if not,

Equation 1. 6         

The experimental variogram of an (categorical) indicator variable I(; zk) is computed in exactly the same way as that for continuous variables as

Equation 1. 7         

It describes the frequency with which two locations a vector, h, apart belong to different categories ¾ k = 1,…,K. This is also the probability of passing from one category to another for a given separation vector.

Ordinary and indicator kriging

Ordinary kriging (OK) and variants of this method are today the most commonly used interpolation methods to estimate the value of an attribute z at unsampled locations xa, using the surrounding n data {z(xa), a = 1,…n}. Kriging is optimal for the estimation of quasi-stationary continuous linear and additive variables from sparse sample data (Burgess & Webster, 1980a, b). Kriging systems are based on the same basic estimator Z*(x),

Equation 1. 8         

where la(x) is the weight assigned to the data z(xa) and m(xa) is the mean (Goovaerts, 1997). Simple kriging (SK) requires that the mean (x) is known and constant across the entire domain to be estimated. For the majority of studies, this assumption is unrealistic. Ordinary kriging (OK) permits the relaxation of this constraining assumption of strict stationarity; the mean is considered constant only within a given local neighborhood W(x) of variable size surrounding the location to be estimated. For the estimate to be unbiased, from Equation 1. 8, it is evident that the spatial weights la should sum to one,

Equation 1. 9         

because this ensures that the mean error is equal to zero

Equation 1. 10       

The individual weights are then chosen so as to minimize the estimation variance , using the variogram to determine the spatial weighting under constraint that the estimator is unbiased Equation 1. 9. The estimation variance is

Equation 1. 11       

To minimise the estimation variance with the addition of the non-bias condition a zero sum expression including a Lagrangian L(u) parameter 2mOK is added to the equation for the estimation variance (Equation 1. 11). The ordinary kriging system (OK) is obtained by setting to zero the partial first derivatives with respect to each l and the Lagrange parameter

Equation 1. 12       

The estimation variance equals,

Equation 1. 13       

The estimation variance is a measure of the local uncertainty and is independent of the data values, that is, it describes the estimation error at a single location irrespective of the surrounding data and the value at the location x0. It can be used to define the optimal grid spacing (Webster & Oliver 1990) and a map of the estimation variance is usually presented with each prediction. Thorough details and properties of the ordinary kriging method can be found in textbooks by Goovaerts (1997) and Wackernagel (1998). The series of papers by Burgess & Webster (1980a-c) are a particularly good introduction to the techniques applied to soil data. Bierkens and Burrough (1993) provide a description of the indicator kriging system, which is discussed in greater detail in the next chapter.

Geostatistical methods have flourished because of the flexible way in which models can be developed that are specific to the data that are being analysed. The constraint on the sum of the OK kriging weights is a good example of this. We cannot reasonably assume that the mean value is uniform across the domain neither do we have complete knowledge of the locally varying mean. The weight constraint incorporates this uncertainty into the model. So typically, for a given estimate the ordinary kriging estimation variance is larger than the simple kriging estimation variance. The choice of which system to use should be defined by expert knowledge of the data and the domain under investigation. The advantageous properties of kriging over other methods of interpolation include:

a)      Efficient use of the data to define the spatial weighting and to predict.

b)      Non-biased, exact estimates with minimum estimation variance.

c)      Estimates can be made for isolated locations or areas of any size using a change of support model (Burgess & Webster 1980b).

d)      Takes into account the sampling design and redundancy of information.

e)      Local uncertainty is expressed through the kriging variance.

Kriging is used to provide optimal unbiased estimates and a measure of the estimation error, not specifically to investigate the uncertainty caused by the combination of spatial variability and incomplete sampling. The predictions of soil properties across large regions from small sparsely sampled datasets can be subject to considerable uncertainty, and a study of this uncertainty may reveal as important as that of prediction itself. Geostatistical simulation methods have been developed to generate many equi-probable realisations (maps) of the RF, which have the same spatial structure as the observed variable, regardless of the density or distribution of the sample data. The spatial uncertainty information that they can provide should represent a useful addition to the prediction and estimation variance map.


When sample data are sparsely distributed, the case for most reconnaissance soil surveys at a regional scale, considerable uncertainty is likely to be associated with each estimate. The estimation variance will be larger for samples that are more isolated, and for variables that have a large relative nugget variance. The rate with which the estimation variance increases with increasing distance from a sample location is determined by the structure of the variogram model (with the same range parameter). Clearly, the rate of increase will be slower when kriging using a Gaussian variogram that is continuous towards the origin, than when kriging using an exponential variogram.

The kriging estimation variance and estimate of a multi-Gaussian variable can be combined to provide conditional probability distributions  that describe the local uncertainty of the estimate at location x, using the surrounding sample data, n. It is then possible to simulate local uncertainty by drawing random numbers from this conditional distribution. But no measure of the spatial uncertainty for a series of M locations can be derived from the single-point ccdfs . Each estimate when considered individually is best (in the least-squares sense), but the map of the local predictions “may not be best in a wider sense” (Goovaerts, 1997). Most importantly, the estimation variance does not describe the spatial uncertainty associated with possible patterns of values.

Simulation methods seek to reproduce global statistics such as the histogram or variogram in contrast to kriging, which seeks to optimise the local error variance (Goovaerts 2000). Simulation methods can be used to investigate spatial uncertainty with fine resolution from sparse samples. When simulating uncertainty information is provided by repeated computation of many equi-probable realisations. In contrast to kriging, simulation randomly selects a value z(x¢j), j = 1,…, N, at location x, from the N-point conditional cumulative distribution function (ccdf) that models the joint (hence spatial) uncertainty at the N locations surrounding the simulated point,

Eq. 1. 14                 

This process is repeated for each location x¢a, under constraints that include exactitude, global reproduction of the variogram and histogram of the random function. This does not mean that each individual realisation has the same histogram and variogram as the constraining data. Individually they can deviate considerably. However on average, if the realisations are considered together they should be reproduced. The individual realisations can also be analysed. Patterns that exist in many are considered likely to occur and those that appear irregularly rejected as improbable. The standard deviation can be calculated across all the realisations for each node to provide a measure of the spatial uncertainty prevalent at each node. A map of this variable could provide useful uncertainty information, for example, to identify regions of transition, where soil properties change rapidly and inversely those regions of homogeneous soil.

If the values at individual locations are the only focus of interest, simple kriging in a multiGaussian framework can provide similar uncertainty information. For local uncertainty investigation this is the perhaps the most efficient method, however with the speed today’s computers attain and the large storage devices that are available there is no reason why these techniques should not be superseded by simulation methods which can provide additional spatial uncertainty information that can be used in a number of ways. In contrast to kriged maps of continuous variables, thresholds can be applied to the realisations of a continuous variable because truncation of a single realisation of a random function generates a random set described by an indicator. The realisation sets can be used as inputs for further modelling activities to elucidate the effects of spatial uncertainty on the model outputs (Addiscott, 1993).

Many geostatistical methods exist to simulate both continuous and categorical (indicator) variables (Journel 1996, Goovaerts 1997). Turning Bands (TB), LU decomposition and sequential Gaussian simulation (sGs) methods can reconstruct the variability of realisations of continuous random functions. Sequential indicator simulation (sis) and truncated Gaussian (TG) methods are capable of generating realisations of random sets (Zhu & Journel, 1993). Both types of simulation method can be made conditional to sample data (Lantuéjoul, 1995). The methods can be further divided into two groups, sequential and non-sequential methods. A brief description of the concepts underlying the most commonly used sequential methods follows.

MultiGaussian Random Functions

The objective of simulation is to derive conditional cumulative distribution functions (ccdfs) for each location. This is simplified if all the ccdfs have the same analytical expression and are described by few parameters, for example the mean and variance. To do so a model that describes the complete multivariate distribution of the random function (RF) must be assumed. The assumption of this model gives rise to the name parametric simulation (Goovaerts, 1997). Combinations of variables Xn, which have Gaussian statistical distributions ~(0,1) have favourable properties of stability through addition ((X1 + XX3 is (0,1)) and through linear combination (lX1 + lXX3  is (0,1)), and this can be generalised,

Equation 1. 15        .

The central limit theorem infers that as , any combination of independent random variables X1, X2,…, Xn, with zero mean and unit variance has zero mean and unit variance,

Equation 1. 16       

Equation 1. 17        .

The variance of a combination of Gaussian random variables remains constant if divided by the square root of n,

Equation 1. 18       

Therefore it is possible to simulate and combine many independent random variables each having the same mean and variance as the original simulated Gaussian random variable ~(0,Var(X)). The multiGaussian RF also has favourable spatial properties. The two-point distribution of any pair of RVs is normal and determined by the covariance function (Goovaerts, 1997). When simulating to many locations the multiGaussian random function will have the same variogram as the individual realisations if the variograms g(h)i used to constrain the selection of random values were the same.

Conditioning multiGaussian Random Functions

Simple kriging lies at the heart of the majority of geostatistical simulation algorithms, because for multiGaussian random variables, the simple kriging error is uncorrelated with either the value of the estimated point or other estimated points, because the simulated values are linear combinations of the original independent random variables (Equation 1. 15). However the residual error is non-stationary, because SK is an exact interpolator. To simulate independent stationary residuals the following algorithm is used.

1)      Non-conditional simulation of RF, Y(x) giving YS(x) at all locations.

2)      Krige using the YS(x) that coincide with data locations as new samples to estimate a new variable YS(x)K, which satisfies,

Equation 1. 19       

3)      Combine the simulated residual and an initial kriging estimate of Y(x), Y(x)K,

Equation 1. 20     ,

This process generates the normal simulated variable YSC(x) that is conditioned to the sample locations and variogram model, and honours the normal score transformed sample values at the sample locations.

The sGs method uses a moving window method and simulates points one after another, using already simulated values to update the neighbourhood and thereby contribute when simulating subsequent nodes. A simple overview of the process:

1)      A random path across the domain is selected.

2)      Krige the Yi(x) within a neighborhood surrounding a location, to give Y(x)K at all locations.

3)      Simulate a variable ~(0,1) and multiply this by sK, and add to either the conditioning data value or Y(x)K.

Simulating indicators

The basic idea behind the sequential indicator simulation algorithm is the use of indicator kriging to provide surrogate probabilities. Indicators for locations can be generated by drawing uniform values on a scale, U[0,1] and assigning either 1 or 0 depending upon the probability of occurrence p,

Equation 1. 21       

The probability of occurrence depends on the values at surrounding locations. The conditional probability of occurrence is equal to the kriged estimate of an indicator. Indicator based simulation algorithms use the indicator kriging estimate as a conditional probability (but one which can take values >1 and <0). If the indicator function is considered as a random set A,

Equation 1. 22        .

The value of the indicator of this set known at the data points. The conditional probability of an indicator given the surrounding data , then the simulated indicator is simply . The most commonly used indicator simulation algorithm is sequential indicator simulation (sis) Goovaerts (1997), that uses a random path and a moving window method similar to the sGs method. Simulated values are used to update the neighborhood information.

The use of indicators methods to provide posterior conditional probabilities can be criticized because the kriged indicators are pseudo-probabilities that can exceed the [0,1] interval. Order relation violations may also occur, when the cumulative probability at larger quantiles is smaller than at lower quantiles. In practice the order relations are usually solved by simply reordering and inadmissible values are typically set to the maximum permissible values.

Figure 1. 3. Cotton fields after harvest (June), showing the permanent ridge and furrow planting scheme.

Figure 1. 4. Natural eucalyptus woodland vegetation typical cover for the Red soil

Analyses and results

Conceptual classification, quantitative fuzzy classification and geostatistical interpolation/simulation, were applied using the ‘traditional data types’ (see box 1), which included visual observations of soil colour and chemical and physical analyses completed on 110 bulked samples of the upper 0-10cm soil layer collected in the field and transported to the laboratory. A total of 80 samples were used to calibrate the models while 30 were used for validation purposes.

Conceptual and quantitative fuzzy classification results were compared qualitatively to determine whether they represented a similar division of the samples. A quantitative comparison was not feasible because information about the profile properties of the soil unavoidably entered into the conceptual model. Only surface samples were included in the quantitative classification to permit comparison of prediction results with those from other methods presented later.

Kriging and simulation were used to provide estimate and spatial uncertainty maps of the continuous soil properties and the fuzzy classification memberships. These results were used to provide a preliminary indication of the precision that could be expected by applying current ‘operational’ techniques for soil mapping. Later they shall be compared with those provided by more elaborate methods, and when additional data is available at validation locations.

The conceptual model

Observations made in the field, interpreted in conjunction with the limited existing information for the region outlined earlier were the basis for a conceptual classification. The photographs on the subsequent pages summarise some of the observational information that provided the basis for this classification. Alluvial sediments in the low-lying plains have been transformed into grey (10 YR 5/2) clay soils, Vertosols (Isbell, 1996), used for irrigated cotton production and extensive grazing (Figure 1. 3, Figure 1. 6). Raised areas of higher relief (ancient dunes or levees) (~30 m) are covered by a mixture of fine- and coarse-sandy red or reddish (2.5 YR/V/C) soils, Rudosols (Isbell, 1996) which are thought to be of aeolian and fluvial origin (Northcote 1979). These soils are uncultivated, and support a natural forest vegetation (Figure 1. 4, Figure 1. 5). In the USDA (1975) classification these two soils would be allocated to the Vertisol and Ultisol classes respectively. These names are used interchangeably throughout the thesis.

Figure 1. 5. Typical Red soil surface, showing a small termite hole. Note the sandy texture and intense red colour.

Figure 1. 6. Cotton soil surface, showing cracking and accumulation of salts at the surface, evidenced by the dominance of salt resistant scrub plants.

A light-coloured (5 YR/V/C), fine-textured ‘intergrade’ soil was observed on the slope breaks or the transition zones between the raised levees and low lying alluvial formation (Figure 1. 7). These transitional or intergrade soil for the region are described as Sodosols, Chromosols or Calcarosols in Isbell’s (1996) Australia Soil Classification, Alfisols in the USDA classification and Luvisols within the FAO system. Close to the present day river channel is a light-coloured fine silty coloured soil that has formed on the river levées. Gypsic horizons are found at depth in both Vertosols and Rudosols.

The A horizon of the Vertosols commonly have diffuse carbonates, while nodular carbonates are more common in the Rudosols. The sand, gypsum and the calcium carbonate are probably of similar aeolian origin. The region offers particular advantages for soil survey methods reliant on remote observations because of its semi-arid sparse natural vegetation cover. Also the land management practices associated with cotton farming leave the soil surface visible from space for long periods throughout the year.

Each sample location was allocated while surveying to one of four soil classes; Cotton, Red, Levée or Transition classes, which correspond to the four soil types described above. Soil colour and texture were the key diagnostic attributes. Figure 1. 5 shows the distinct colour and sandy texture of the surface of a typical Red soil.

Figure 1. 7. A Transition soil surface, cracking like the Cotton soil, but with iron in its oxidized form like the Red soil.

In contrast Figure 1. 6 shows the surface of a typical uncultivated Cotton soil. There is evidence of cracking caused by the shrink swell processes operating in this vertisol. The bright light colour is evidence of the build-up of salts at the surface. These two soil types can be considered as being modal for the region. The Cotton soil exists in the alluvial plains, the Red soil on the sandy dunes of an ancient landscape. The Cotton soil is homogeneous having uniform parent material, age and having undergone similar pedogenesis. The Red soil is more variable. The origins of the sandy dune materials is likely to have been more divers (aeolian and fluvial) and deposition is likely to have occurred during different periods over time.

Figure 1. 8. Red soil profile (showing depth of approx. 1.50 m), with duplex clay horizon staring at 110 cm depth.

Figure 1. 9. Well drained transition soil profile (showing depth of approx. 1.50 m), showing cracking in the upper 1 m.

The Transition soil class is true to its name. Figure 1. 7 shows an example of the Transition soil surface. Its colour is ‘intermediate’ between the Red and Cotton soils and often samples allocated to this class were located along the boundary (or transition zone) between the two modal soil types. Relatively few samples of the Levee soil were collected. Samples allocated to this class were situated close to the present day river channel. In the conceptual model, these soils are young soils, with variable texture profiles that result from the heterogeneous depositional environment of the Murray-Darling River. Statistics for the laboratory analyses of chemical and physical analyses are presented in the appendix. The range of values observed for clay content are clearly in error. The cotton soils on the low-lying alluvial plains are vertisols. Hand texture analyses completed in the field indicated a clay content greater than 30% for these samples. The error was caused by the laboratory measurement method.

Figure 1. 10. Box and whisker plots showing the distribution of values for the 0-10cm sample for each soil type. (A = Levee, C = Cotton, R = Red, T = Transition). The width of the boxes is proportional to the number of samples allocated to each soil type, the notch indicates the mean, horizontal bars +/- 2 standard deviations, and the dots possible outliers.

The processes of physical and chemical disaggregation did not successfully disperse and separate the clay particles from within larger aggregates. A more vigorous procedure should have been used to breakdown the strong aggregates caused by salinity and the presence of carbonates and probably gypsum. Nevertheless the remaining reliable chemical and physical laboratory analyses provided a means of validating the conceptual model and testing its predictive ability.

Figure 1. 10 shows the statistical distribution of each measured variable for each soil type allocation. These plots show that the conceptual classification groups soils with similar properties. On average Cotton soils have larger silt, clay, organic matter, iron electrical conductivity and pH. The Red soil is sandier and has characteristically strong hue (intensity of colour) that distinguishes this soil from the dull grey Vertisol. Despite the strong colour of the Red soil, that would suggest intuitively that on average this soil has the larger iron content, it has less iron on average than the Cotton soil. The Red soil is sandier and the iron that is present coats the predominant large-grained sand particles and is therefore far more apparent. In contrast the Cotton soil is a vertisol with high clay content. It is regularly flooded, either as part of the on-farm flood and drain irrigation or when the Murray-Darling River breaks its banks. Iron oxides in the Cotton soil are reduced and intimately mixed with organic matter from the sediments contained in the river water, giving the soil a blue – grey colour.

The transition soil forms an intergrade, with ‘intermediate’ sand and iron contents. On average cotton soils have higher organic matter content, pH and ECe than the red soil. These variables all show a wide dispersion of values about the soil type means. The results of a one-way ANOVA (Table 1. 1), using soil type to estimate the mean attribute value for each soil type confirmed that the conceptual model was capable of estimating the mean values of the pH, iron content, silt and sand content with a high degree of significance.Clay is poorly predicted because of measurement error. The distinction on the basis of soil type was less clear for electrical conductivity and organic matter content. They do not appear to affect the colour of the soil or reflect those properties that do in such a way that a human observer notices. On the other hand the sand and iron content contribute significantly to determining soil colour (Stoner et al. 1980), which is easily incorporated into a mental model of soil variability and was therefore diagnostic when allocating soils in the field.

Quantitative classification

During fuzzy classification the primary choice lies in determining the optimum number of classes (if the fuzzy membership coefficient is fixed). The conceptual model developed during sampling had four classes, two of which dominated the landscape, the Cotton and Red classes.

The Transition soil was identified as an intergrade and the Levée soil (poorly represented by the sample dataset) as resembling more the Cotton and Transition soil than the Red soil. Tests to define an appropriate fuzzy classification of the sample dataset indicated that only two classes could be adequately characterised. Figure 1. 11a plots the samples in the space of the first two transformed components of variation and from which two clear groupings are evident. Figure 1. 11b is a silhouette plot.



Explained variance


(Pr > F)

Silt content (%)



Sand content (%)



Iron content (mg/kg)





6.43x 1010***

Electrical conductivity (Sm-1)



Organic carbon (%)



Clay (%)



Table 1. 1. Summary of one-way ANOVA using soil type (L, C, R, T) to predict soil properties. Level of significance: *0.1, **0.05, ***0.001.

Each horizontal bar corresponds to one sample, the width of the bar is calculated by

Equation 1. 23        ,

Where a(i) is the average dissimilarity between an individual and all other points of the cluster to which i belongs. Then for all clusters C, d2(i,C) is equal to the average dissimilarity of i to all observations of C. The smallest of these is denoted as b(i), and is interpreted as the dissimilarity between i and the cluster to which it has the second largest resemblance. The overall average silhouette width is the average of s(i) over all observations i and can be used to indicate the most suitable number of clusters. The larger the value the more reasonable is the choice of k. Figure 1. 11 shows two groups: the group to the right correspond to the Cotton soil, with 85 % of samples correctly assigned when in the field, the left hand group correspond to the Red soil (75 % same class). Several samples that bridge the gap formed by the two overlapping circles defining the fuzzy class boundaries (membership = 0). Another group of samples can be identified from the silhouette plot (Figure 1. 11b) having negative values of s(i). These samples bear little resemblance to either of the two classes.

Figure 1. 11. Fuzzy classification of sample data (a) partition plot and (b) silhouette plot. The symbols in (a) indicate the closest crisp allocation; r Red, ¡ Cotton.

The former correspond to samples that were allocated to the Transition soil class, which can indeed be considered as an intergrade soil. The majority of the samples with negative s(i) values were allocated to the Levée soil group when in the field. Given that the fuzzy classification is generated using the sample data, the sample data should correlate with one or more of the fuzzy membership values. Unfortunately, the fuzzy memberships from the classification presented above were only very weakly correlated with the individual soil properties and regression was not feasible.

The total variance explained by the fuzzy classification is small (14.18%, see Figure 1. 11) because the division failed to account for the variance caused by the most statistically dispersed variables electrical conductivity and organic carbon content. Nevertheless, these results are testament to the validity of the initial subjective allocation to these classes and quantify the physical basis for the allocations. In-situ soil type allocation based on expert knowledge is a useful and cheap practice. The eye-brain combination is capable of adding considerable value to each sample observation and measurement and also isolating sources of noise that contribute little to the explanatory power of the model. The fuzzy classification mimicked this process quantitatively.

Interpolation and simulation

Structural analysis and interpolation

Histograms for each of the in-situ soil variables are presented in Figure 1. 12. They divide the data into two distinct groups. None of the variables has a distribution that could be described as normal, but the variables describing horizon thicknesses (Figure 1. 12h-j) are very positively skewed. The ordinary kriging estimator makes no assumptions about the shape of a variables statistical distribution, but OK estimates may suffer from exaggerated systematic bias caused by skew (Saito & Goovaerts, 2002). The positive skew of the horizon thickness variables is caused by the discontinuous spatial cover of the horizon attribute. A detailed analysis of these ‘discontinuous’ horizon attributes was completed separately so that various kriging methods could be tested and an appropriate method selected. The results of this analysis are presented in chapter 2. For the other variables, variograms were calculated in four directions with a lag spacing of, 2500m for 8 lag distance using the entire dataset (Figure 1. 13).



Model Parameters









a1 (m)

Iron content (mg/kg)






Organic carbon (%)






ECe (Scm-1)












Sand and silt content (%)






Clay content (%)






Fuzzy class memberships 1 + 2






Table 1. 2. Variogram models of soil properties of the Bourke region. For an example of the notation, Stable (2) refers to a stable model with shape parameter with value 2

When kriging the dataset was divided into training and validation sets, however the small size of the dataset meant that no samples could be set aside for variogram analysis, without introducing concerns about variogram stability. Predictions are always subject to uncertainty, because of approximations in the model, and error in the measured values the model parameters and variables. Variogram estimation requires at least 100 data points (Webster & Oliver, 1992). Estimating the variogram using just the validation points would have introduced an extra source of error that risked hiding the effects of the estimator on the quality of the predictions.


The experimental variograms for each variable showed no evidence of geometric or zonal anisotropy, but all were erratic so a single non-directional variogram was modelled for each variable (Figure 1. 13). The fitted models fell into one of two groups. The first included the attributes pH, ECe and clay content, which have one short-range structure that does not exceed 5000 m and a large relative nugget effect (in particular clay content because of measurement error Figure 1. 13g). The range of the structural component of the second group, which included iron and organic carbon content, and (fuzzy) class membership was of the range 8000 – 10000 m.

These parametizations suggest that short-range processes operating across the landscape are responsible for the chemical properties of the soil, and that these are in turn related to the clay content of the soil. This makes physical sense, because clay is the chemically active component of the solid soil fraction which controls to a large extent the pH buffering and salinity status of the soil.

However this interpretation must be considered carefully given the measurement error of the clay content data. Sand and iron content appear to vary across larger spatial scales that reflect the geomorphological characteristics of the region, the alluvial plains and raised dunes. They share the same spatial characteristics as the regional soil classification.

Each of the raw variables was transformed into a zero mean Gaussian variable with unit variance for the purposes of simulation. Authorized models were fitted to each of these variables normal score experimental variograms. The models closely resembled those of the corresponding raw variables but as expected having lower nugget variance.

Figure 1. 12. Histograms of soil chemical properties: (a) pH, (b) ECe (Sm-1), (c) organic matter (%), (d) Fe (mg / kg) (e) sand, (f) silt and (g) clay content, (h) diffuse, (i) nodular and (j) gypsum horizon thickness.


Figure 1. 13. Variograms of soil chemical properties: (a) pH, (b) ECe (Sm 1), (c + g) organic matter (%), (d) Fe (mg / kg) (e) sand, (f) silt and clay (g) content (units = %), (h) fuzzy class members 1 and 2 (which have identical variograms).

The variogram models of the raw variables were used to predict each in-situ variable to the set of validation points. The same models were used to krige the prediction and estimation variance maps, using the entire dataset.

The validation dataset was used to calculate the mean square error and the correlation between the true and predicted values, both measures of the quality of the fit.

Kriging was performed in a moving neighbourhood that was large enough to ensure that all points in the domain could be estimated. The variograms of the normal score transforms of the soil variables were then used to simulate using sequential Gaussian algorithm (sGs). Tests were completed to identify the minimum number of simulations acceptable (50 to 5,000).

The test results changed little above 100 realisations. One hundred simulations of each variable were generated and the results post-processed to calculate the simulation standard deviation. Various parametizations of the moving window were tested. Finally only a maximum of five simulated data values were included in each window. The quality of the OK estimates were evaluated using the validation statistics (Table 1. 3) calculated at the 30 validation sample locations. The root mean square error is expressed both in the original units, but also as a percentage of the average value of the sample data. The correlation value measures the degree of systematic prediction bias, caused by the overestimation of low values and underestimation of high values.

The sparse nature of the dataset inevitably meant that the prediction quality would leave much to be desired. In the absence of many constraining data the predicted values tended towards the mean value. The most accurate and least biased predictions were for sand, silt and iron content and pH, the worst for organic carbon, ECe, and clay content. Again the quality of the latter is the poorest, further proof of measurement error, if needed.


Iron content




Organic carbon (%)


Sand/Silt content (%)

Clay content (%)

Root mean square error














% error







Table 1. 3. Validation dataset prediction mean square error and correlation statistic between true and predicted values. % error = rmse / variance


Prediction maps for each variable are presented in Figure 1. 14 through figure 1.17. As the variogram models indicated the prediction maps were to fall into two groups, a function of the range of the fitted variogram model. The poorly predicted variables ECe and organic carbon (OC) had short-range variograms (~4000 m). Their prediction maps failed to reveal a strong regionalisation. In contrast the prediction maps of iron content, pH, texture and fuzzy class memberships, modeled using longer range variogram structures (8000 – 10000 m), show a clear regional scale spatial structure.

The density of the sampling is sufficient to characterize the regional scale variability of these attributes. The fuzzy class membership predictions show three regions with high Red soil memberships (ultisols-alfisols) corresponding to the southwest-northeast sandy dunes and ancient fluvial levees. (member 1, figure 1.15c d). A sinuous swath of vertisols soils (member 2) traverses the study region.

The prediction maps for sand and silt content (figure 1.17a and b) resemble the fuzzy class membership maps. The clay content maps must be evaluated with caution because of measurement error. The sand and silt maps provide a clear indication of a swath of soils with high silt content (and consequently high clay content), corresponding to the high class 2 membership predictions.

                                    (a)                                                      (b)

                                 (c)                                                         (d)

                                  (e)                                                          (f)

Figure 1. 14. Ordinary kriged (OK) map (left column) and sequential gaussian simulation set (n = 500) standard deviation (right column), for pH (a + b), ECe (c + d), organic matter (e + f).


(a)                                                                  (b)

                               (c)                                                                      (d)

                              (e)                                                                       (f)

Figure 1. 15. Ordinary kriged (OK) map (left side) and sequential Gaussian simulation set (n = 500) standard deviation (right side) for iron content (a + b) respectively for Red soil member 1 (c + d) and Cotton soil member 2 (e + f).

This sinuous feature describes the location of depressions where fines have accumulated and where vertisols have formed.

Estimation variance maps for the best and worst predicted variables are illustrated in Figure 1. 16. The estimation variance map provided no more information other than that which is evident from a straightforward diagram showing the sample locations combined with a knowledge of the variogram model. Absolutely no difference can be observed between the two maps.

They are an indication of sample density and as previously discussed the choice of variogram, its parametization and the sample design, but provide no variable specific uncertainty information, only model specific uncertainty information.In contrast the simulation maps for each variable differ considerably and can be interpreted very intuitively. The maps indicate the regions where uncertainty is prevalent, caused by the heterogeneity of sample values.

The uncertainty values are independent of the sample location distribution or density. The patterns of values at the constraining sample locations determine the structure of the variable specific spatial uncertainty.

                                   (a)                                                               (b)

Figure 1. 16. Maps of the OK estimation standard error for (a) FeOK and (b) pH*OK, illustrating the inutility of expressing variable specific spatial uncertainty using the estimation variance, when the sample locations are sparse and irregularly spaced.

The simulation standard deviation maps of the fuzzy class memberships show how region is dominated by vertisols (member 2). The high degree of certainty of the presence of this soil is illustrated by the uniformly low simulation standard deviation values (figure1.15f), with spatially uniform the cotton soils.

This finding, suggests that cotton production could be extended to similar regions with a strong likelihood of finding similar soils.

Figure 1. 17. Ordinary kriged (OK) map for (a) sand (b) silt and (c) clay content.

In contrast the spatially heterogeneous and localized nature of the Red soils was apparent from figure 1.15d. This map shows only two isolated regions to the northeast of the region where member 1 (red soil) uncertainty is low. Figure 1.15e and f show low values of uncertainty to the center of the study area corresponding to the center of the region dominated by the Cotton soil.

In Figure 1. 14a c d the simulation standard deviation maps reveal a heterogeneous short-range structure. This indicates that the sampling strategy is insufficiently dense to capture much of the spatial variability of these properties. The reasons for the variability are not regional, but field scale and related to land management practices and landcover type. Indeed, these soil properties are likely to be significantly modified by very local activities and soil processes, in particular the waterlogging, flooding and draining on the cotton fields or local salinisation close to raised water storage tanks.


The purpose of this first chapter was to provide an introduction to the study area and to investigate the data. A reminder of the first two objectives outlined earlier :-

1.                  Describe and test the validity of a conceptual soil classification

2.                  Compare the basis of this classification with its quantitative counterpart

A considerable amount of information and data was gathered during just a single week of soil survey. The process of visiting all of the sample locations, augering and recording observations of the soil proved very informative. It became evident that three soil types dominate the landsurface; large homogeneous vertisol plains, raised dunes with less unfertile ultisols with large tracts of more fertile alfisols forming the transition between them. Ranges of soil physical and chemical properties strongly reflected this division, when not contaminated by measurement error. Unfortunately the clay content values could not be trusted. The conceptual classification of Red, Cotton and Transition soils is by no means a formal one, but it serves its purpose. It is a shorthand means of communicating that which has been learnt during the survey. As more data is collected and information synthesised about the soils of the study area, these expressions can be replaced by more detailed and precise allocations within formal systems. Importantly though, these results confirm that experience and observation can combine to permit the development of empirical and conceptual soil classification models, which have a predictive capability for those soil properties that clearly distinguish the soils of a region.

Lacking sufficient data to allocate each sample accurately to a formal classification system, it was very important to validate the underlying concepts of this rough and ready model using a quantitative technique. Fuzzy classification methods permitted the quantification of a conceptual model. The results of this classification strongly reflected the original division of soils made in the field. The same variables were diagnostic for both the conceptual and quantitative models. Both models distinguished between soil types primarily on the basis of their colour, their iron and sand content. Colour was the main diagnostic soil property, because the relative proportions of sand and iron content play a large role in determining it.

The next objectives were to:

1.                  provide estimates and maps of important soil properties and

2.                  investigate the spatial variability of the soil properties.

The results of any soil classification must be transferred to a spatial representation to fulfil their predictive role. Ordinary kriged estimates represent the ‘control’ in this study against which the relative merits of other methods will be evaluated. Crisp soil classifications stratify the landscape into uniform homogeneous regions. As discussed, this representation contrasts with our understanding of the soil as a continuous phenomenon, but could represent an acceptable level of simplification, depending upon the scale of the soil survey. That which may be perceived as continuous change at one scale of observation, could accurately be represented at a lower resolution as a discontinuity without much loss of generality. However, experience has shown that when a crisp classification is transposed into the landscape the spatial representation is fragmented and bears little resemblance to that which we intuitively know to be correct. By reducing the fragmentation in character space caused by inappropriate exclusion from a given class, spatial contiguity in the spatialisation of a classification can be generated (McBratney et al. 1992).

The fuzzy classification membership maps clearly show the expanses of Red soils, formed on dunes to the north-east and west, also the isolated occurrence of a not dissimilar soil to the center west of the study area. The large homogeneous expanse of Cotton soil, almost exclusively exploited for cotton production, is identified as running in a swath from the northwest to the southeast. This zone lies perpendicular to the present course of the Murray-Darling River, whose presence is indicated by a thin swath of low membership 2 values passing from the northeast to the southwest.

Mapping the continuous fuzzy class memberships has provided a quantitative method to mimic the minds ability to express the vague human concepts that underlie conceptual models of soil spatial variability. The procedure of combining the many univariate measurements into a multivariate representation has isolated the regional components of spatial variation. The process has filtered out the numerical fluctuations associated with the variables that describe the individual soil attributes, and provided a continuous numerical representation of their communality. As a result, the maps of the kriged the fuzzy class memberships were less affected by individual sample variation that are the exception, but by the features common to all samples, global features of the dataset.

Finally, simulation provided a useful alternative to the estimation variance map; the simulation standard deviation maps provided spatial information that was a) specific to each variable, b) a function of the chosen variogram model and, c) the relative spatial distribution of the values. These maps permitted the identification of regions where the value of the attribute could be expected to be relatively constant, and the contrasting transition zones of rapid change and spatial heterogeneity. Where spatial uncertainty was greatest more samples are required to provide a uniform prediction quality. The value of the simulation standard deviation could be used to choose the location of these samples.

If anything this chapter has shown how pedological observations can provide very useful spatial information at little cost. It remains to explore the information content of the horizon observations that were not included in the fuzzy classification of surface samples. Preliminary investigations of this data suggested that alternative methods would be required to provide accurate spatial predictions. The next chapter addresses these issues in further detail.

Chapter 2

Predicting and simulating the
geometry of soil horizons


Pedological observations have traditionally formed the bedrock of soil survey and shall continue to do so as long as resources are scarce or until sufficiently precise and versatile, cost effective technological alternatives are developed. During the reconnaissance stage of soil survey the focus of attention has to be on providing an overview with which to create the framework to describe and further investigate the soil spatial variability. Pedological observations are essential to provide important ancillary (spatial) information for inclusion within the conceptual model. When cost constraints limit the use of detailed analyses to the upper layers of the soil observational measurements may be the only way to provide a first overview of the profile properties of the soil. They may also provide information invaluable in its own right. Soil variation has long been described with reference to diagnostic horizons, the presence or absence of which determines the allocation to a given soil type in formal classification systems. Within the horizon the soil properties (concentrations, density etc) are assumed to be uniform. This represents a reasonable approximation for studies at coarser scales, typical of GIS , regional or global studies if care is taken to accurately represent the aerial and vertical extent of the horizon correctly.

Concerns about global climate change require estimates of various components of atmospheric and terrestrial carbon. This is particularly important with respect to terrestrial inorganic carbon, which is a major repository of global carbon dioxide (CO2). Global and regional estimates of inorganic carbon, particularly soil carbonates are particularly unreliable and in many cases do not take into account the spatial variability of the occurrence or the amount of carbonates in the soil horizons or examine the spatial variability of the predictions. Gypsum on the other hand is important for soil management. Worldwide estimates of the extent of soils containing gypsum range from 90 million (FAO, 1993) to 207 million ha (Eswaran and Zi-Tong, 1991). These estimates are as unreliable as those of terrestrial carbonates. There is, therefore, the need to design appropriate approaches to provide an inventory of soil carbonate and gypsum, particularly at the regional level, which could be aggregated to give estimates at the national, continental and global scale. Of specific concern in the Bourke region are the possible effects that carbonates and gypsum may have on soil stability with the introduction of irrigation.

It is increasingly important to provide uncertainty information concerning these predictions that can be propagated through soil, crop and climate models (Addiscott, 1993). Uncertainty arises in any situation where there is a lack of complete knowledge. Spatial prediction uncertainty is caused both by a lack of precision when collecting data and incomplete knowledge, an unavoidable characteristic of data provided by any measurement process and particularly prevalent when the data is collected across a large area from few sparsely distributed point samples.

Accurate spatial prediction of soil horizon thickness implicitly requires accurate prediction of the locations where a horizon is present. Soil horizon thickness can be expected to lie within a limited range of values, whereas the true horizon lateral extent can differ considerably from that indicated by the samples collected only at a few points across a region. Therefore the uncertainty associated with the prediction of the horizons lateral extent is likely to be far greater than that of predicted horizon thickness. Geostatistical interpolation and simulation techniques can be used to provide spatially distributed predictions and uncertainty measures, but the choice of which technique to use should be defined by our understanding of the pedology of the region and the data themselves. For example, ordinary kriging alone, may not be the most suitable technique to estimate soil properties that are not spatially continuous across the region, at the scale of observation, and hence require a prior delineation of the zones within which the property can be reasonably considered to be continuous (Voltz & Webster, 1990).

Soil horizon thickness is a continuous variable often having a positively skewed distribution. This is caused by many zero values, where no soil horizon is present, and few extreme values. Ordinary kriging (OK) predictions of positively skewed data are often biased, caused by the underestimation of large values and overestimation of small values. Several geostatistical techniques have been developed to overcome the prediction bias caused by positive skew, these include multi-Gaussian kriging (MG) and indicator kriging (IK). A third method consists of using the kriged prediction of an indicator of horizon presence to delineate regions where the soil horizon is most likely to exist, and which can be used to mask the continuous estimate. Zero values of the continuous variable delineate the region where no horizon is present.

At a regional scale of mapping and typically within GIS the boundary between regions of horizon presence and absence are as discontinuities. The prediction of local soil horizon thickness could consist of an initial prediction of the regions where the soil horizon is present, and within which it can reasonably assumed to be continuously spatially distributed, followed by an optimal estimation of the continuous variable within these zones. In some respects this is a similar approach to that suggested by Voltz & Webster (1990), the main difference being that the data themselves are used to delineate the regions of horizon presence and absence instead of a pre-existing soil map. The benefits in terms of prior data requirements are evident, in particular for regions where no prior observations or map delineations exist.

Any statistical prediction is subject to uncertainty (Journel, 1989), but this is exaggerated when the number of samples is small and the distance between samples large; a feature typical of regional scale reconnaissance soil surveys. By recreating sample variability simulation techniques can provide information describing the spatial uncertainty of the attributes across the study area additional to the local measures provided by kriging, such as the estimation variance. Again, there are many geostatistical algorithms available to simulate both continuous and categorical variables and the results that they provide can differ (Journel & Isaaks 1984, Goovaerts, 2001).

In this chapter, methods to estimate soil horizon attributes thickness and lateral extent using a simulated dataset that shared the same statistical and spatial properties exhibited by a real dataset were compared. The objectives were: (i) to apply the method that performed best to the real in-situ data and (ii) then to investigate the uncertainty associated with the prediction of the soil horizons lateral and vertical extent by geostatistical simulation.

Kriging continuous attributes with positive skew

The ordinary kriging (OK) system is derived under the assumption that the mean m(xa) is constant but unknown within a local neighbourhood W(x). Therefore the local mean must be estimated for each z(xa). This estimate is not robust to the presence of extreme values. Systematic prediction bias can be severe for positively skewed data. A single extreme value can impact the estimate of the local mean and cause under and overestimation of the tails of the z(x) histogram. This is known as the ‘smoothing effect’. Smoothing increases with increasing skew, relative nugget effect and with decreasing sample density (Goovaerts, 1997).

Techniques have been developed to reduce systematic prediction bias caused by skew, namely multi-Gaussian kriging (MG) and indicator kriging (IK). Both require a non-linear (back) transform of the data prior to and after kriging. The first step of multi-Gaussian kriging is the normalisation of the z(xa) to generate a variable with a Gaussian distribution. This facilitates inference of the mean when kriging. The kriging estimates Y*SK (xa) are then back-transformed to provide the estimate Z*MG(x).

Indicator kriging involves variogram modelling and kriging for the suite of indicators defined using equation 1. The prediction Z*IK(x) is obtained by multiplying the estimate at each threshold zk by the threshold value and summing the results for all K. Nevertheless, because kriging is a form of weighted spatial averaging OK, MG and IK can all be expected to produce a smoothing effect of the true Z(x) histogram. When kriging the high-frequency short-range variability represented by the nugget effect is filtered. Extreme and erratic values are not generated other than at data points. As the sample density decreases and the predictions tend toward the mean value the smoothing effect can be expected to increase. Therefore none of these interpolation methods can be used to classify zero and non-zero zones based on the predicted values. Indeed if any cut-off is applied to the kriged estimate to provide an estimate of the areal cover of the horizon it too will be biased. Instead an estimate of the probability of the threshold value occurring is required.

Kriging categorical attributes to delineate zones of presence

Consider a categorical indicator defined on the basis of presence or absence. The sample mean of the indicator is a prior estimate of the probability of the category’s presence. The kriged estimate of the categorical indicator represents its conditional expectation, which in turn equals the conditional probability that the category sk is present,

Equation 2. 1         

where |n stands represents the conditioning of the prediction by the surrounding n data surrounding a location x. Because the estimate is can be interpreted as a probability it is possible to threshold and thereby to indicator code the kriged estimate,

Equation 2. 2         

The resulting indicator variable delineates the most probable zones of horizon presence and absence. In the absence of prediction bias a reasonable classification criterion would ensure that the estimated lateral extent should equal the global sample proportion. For an unbiased estimate, values smaller than 0.5 indicate a greater likelihood of absence than presence and vice-versa. If available other sources of information could be used to inform this choice, or thresholds can be applied at various values to provide an indication of the uncertainty of the predicted aerial extent.

Methods and Analyses

A sensible approach to the estimation of local horizon thickness should involve two-steps; an initial delineation of the zones where a horizon is most likely to be present and an optimal kriging of horizon thickness within these zones. This hypothesis was tested using a simulated dataset because the in-situ dataset consisted of relatively few sparsely located samples (100 auger observations). Under such circumstances the comparison of prediction methods on synthetic data with similar statistical properties as the real data represented a useful alternative to cross-validation, for hypothesis testing and investigating uncertainty caused by the choice of technique. The construction of a simulated dataset necessarily started with a detailed analysis of the true sample data to ensure that the simulated data had similar statistical and spatial properties.

The second part of the analysis involved a comparison of the tools that were available to provide information describing spatial uncertainty. In the first instance the measures provided by kriging and simulation methods were evaluated and compared. Finally the difference between the measures of uncertainty of horizon lateral extent provided by parametric sGs and non-parametric sis algorithms are discussed.

Structural analyses of horizon attributes

The in-situ variables used were presence (or absence) of: 1) diffuse CO32+ (D), 2) nodular CO32+ (N) and 3) gypsum (G). We refer to the categorical indicator variables as I(D) for the diffuse CO32+, I(N) for the nodular CO32+ and  I(G) for gypsum. The horizon thickness values are DT, or NT, or GT, respectively, with subscript T indicating the thickness. The sample data indicate that, the diffuse carbonate horizon is more prevalent than horizons with carbonate nodules.

Gypsum horizons are relatively uncommon, being present at only 20 % of the sites. However, the lack of an explicitly defined test method and the reliance solely on visual recognition of its presence as crystals coatings found within the profile probably led to an underestimation of gypsum.

The mean depths of the upper depth limits followed the expected pattern with diffuse carbonates typically occurring closer to the surface than nodular carbonates, and the latter forming above gypsum concretions



Variogram Model


Structure 1

Structure 2






Practical range

Diffuse (D)






Nodular (N)






Gypsum (G)




































Table 2. 1. Variogram parameters of in-situ horizon indicators (D, N, G), the normal score transforms of in-situ horizon thickness (DTg, NTg, GTg) and the simulated transect thickness TSIM and lateral extent PSIM (discussed later).

The diffuse carbonate horizon often extends into the nodular carbonate horizon as expected. The correlation between indicators suggests that nodular calcium is often co-located with diffuse carbonates (r = 0.45). Their spatial co-dependence is weak and intrinsic. Gypsum shows a very weak negative intrinsic correlation with diffuse and nodular carbonates (r = - 0.09). Intrinsic spatial correlation occurs when the properties have the same basic variogram structure and proportional sill variance.

Variogram models of the thickness variables for each horizon are shown in Table 2. 1 and Figure 2. 1. Each was modelled using a combination of a nugget effect and an exponential model. The nugget effect accounted for over half of the modelled semivariance for each horizon type, being largest for the nodular horizon and smallest for the diffuse carbonate horizon. Gypsum and diffuse carbonate horizon exponential models had a similar long range (15300 m and 24280 m respectively), in contrast to that of nodular carbonates, which was approximately half this distance.

For the purposes of multi-Gaussian kriging and simulation, each thickness variable was normal score transformed and the variogram fitting procedure repeated. The resulting models all had smaller relative nugget effect and greater spatial continuity as evidenced by the longer range of the fitted exponential odels.

Figure 2. 1. Experimental (dots) and modelled (thick line) variograms for (a) diffuse carbonate indicator, (b) diffuse carbonate thickness, (c) diffuse carbonate thickness normal scores, (d) nodular carbonate indicator, (e) nodular carbonate thickness, (f) nodular carbonate thickness normal scores, (g) gypsum indicator, (h) gypsum thickness, (i) gypsum thickness normal scores. Simulating a dataset with realistic characteristics.

Variograms were also fitted to the (categorical) indicator variables of each horizon. Again each was modelled using a combination of a nugget effect and exponential model. Interestingly the exponential model ranges of both diffuse and gypsum indicators were approximately half the length of the corresponding thickness variables, while that of the nodular carbonate indicator was twice the length of its thickness counterpart. This suggests that nodular carbonates are widespread but extend to very different depths over short distances.

In contrast diffuse carbonates and gypsum are localised and dispersed. For each variable, the variograms of indicators defined at increasingly large decile values showed a typical destructuration effect (Figure 2. 2). For a cut-off greater than 1 m the models were almost pure nugget effect.

Thirty values were drawn at random from a Gaussian distribution N (0, 1) at equidistant locations along a hypothetical 25 km long transect. The variogram model , of the normal-score transformed in-situ thickness variable, was then used to simulate a variable TSIM g to 500 points along the transect using the 30 samples to constrain the sequential Gaussian algorithm. The Gaussian simulated values were back-transformed to the raw space of the in-situ diffuse carbonate horizon thickness. A single realisation TSIM, was then selected having similar sample variograms and histogram as the diffuse carbonate horizon (a).For estimation purposes a variogram , was adjusted to the chosen realisation to be used for ordinary kriging. A second variogram , was fitted to the Gaussian transform of the same realisation, to be used for multi-Gaussian kriging. Thickness TSIM, was also indicator-transformed at values corresponding to the deciles of the cumulative distribution. Variograms were fitted to each indicator, necessary for indicator kriging.

Figure 2. 2. Indicator variograms for selected cut-offs applied to diffuse horizon thickness, zk (a) z1 = 0.06 (b) z2 = 0.15 (c) z3 = 0.32 (d) z4 = 0.51 (e) z5 = 0.67 (f) z6 = 0.83 (g) z7 = 1.13 (h) z8 = 1.27 (i) z9 = 1.36. (REMOVED)

Finally a categorical indicator was defined to create the simulated horizon lateral extent PSIM,

Equation 2. 3         

and a variogram model, , was adjusted to this variable to be used for kriging the indicator of presence / absence (categorical indicator kriging). From the 500 simulated values thirty were selected to be used for kriging. The remaining 470 were set-aside for validation purposes.

Comparing prediction methods

In a first step, OK, MG and IK were used to provide estimates of simulated horizon thickness, TSIM. The predictions T*SIM OK, MG, IK , were compared against the true transect to determine which showed the least systematic bias. The results were also used to highlight the error that accrues if any of these methods is used to classify zeroes and non-zeroes.In a second step, the simulated horizon lateral extent PSIM, was estimated by kriging the indicator of horizon presence . The conditional probability of horizon presence provided by the kriged categorical indicator , was thresholded to provide an estimate of the simulated horizon lateral extent,

Equation 2. 4         

This estimate was used to mask the kriged estimate of simulated horizon thickness by calculating the product of the thickness estimate T*SIM(x) and the indicator coded estimate of horizon presence [I(; PSIM)]*,

Equation 2. 5         

The improvement in the quality the predicted thickness TSIM obtained by using the prior prediction of the probability of horizon presence PSIM to delineate the zones where the thickness value is valid, was compared by the predictions against the true the exhaustive simulated dataset.

Simulating the in-situ variables to characterize spatial uncertainty

Two runs each consisting of 500 simulation realisations generated: (a) constrained sGs of soil horizon thickness following a normal score transform, and (b) constrained sis of horizon presence. The local neighbourhood was conditioned to use no greater than 25 % simulated values when drawing newly simulated values. The results were then post processed. The mean and standard deviation were calculated for each node across the realisations (E-type estimate). In addition, each realisation of horizon thickness, simulated using sGs, was indicator coded according to Equation 2. 3. This provided a second set of realisations of the indicator of horizon presence, which was subsequently post-processed to find the E-type uncertainty measures.


Comparing prediction bias using a simulated dataset

Figure 2. 3 and Figure 2. 4 show the true and prediction error histograms and transects, Table 2. 2 the error statistics. They provided adequate information with which to identify the ‘most appropriate’ kriging method to use with the real data. The kriged transects (b-e) are all smoother than reality (Figure 2. 4a).

Both the prediction transects and error histograms (Figure 2. 4a-d) show that this is caused primarily by the overestimation of zero values, where no horizon is present and to a lesser extent by the underestimation of large values. Of the OK, IK and MG estimates both OK and IK predicted the presence of a horizon along the entire transect.

Table 2. 2. Statistics of simulated, estimated and error transects. Lateral extent error statistics represent the % of the transect misclassified.









Extent (%)

Simulated statistics
























Estimate statistics








































Error statistics








































The MG estimate predicted a horizon presence along 90 % of the transect (Figure 2. 4b). In reality the horizon is present along only 59 %. As expected the OK estimate, with no prior transformation of the data, was the most biased. Surprisingly the IK estimate was equally biased. We suggest that a lack of non-zero indicator values for certain thresholds was responsible (Saito & Goovaerts, 2000).

Clearly, where no mask was imposed on the thickness estimates, the lateral extent of the horizon PSIM and local thickness TSIM were severely over-estimated Figure 2. 4 e shows the optimal estimate obtained by masking the least biased of the horizon thickness estimates, provided by MG kriging.

Figure 2. 3. Simulated transect thickness estimate error histograms, (a) ZSIM – Z*SIM OK, (b) ZSIM – Z*SIM IK, (c) ZSIM – Z*SIM MG, (d) ZSIM – Z*SIM OPT.

A comparison of Figure 2. 4 c and d shows that imposing the mask has reduced the slight overestimation of true zero values. As a result the estimate of the simulated horizon lateral extent , matches exactly the true lateral extent PSIM.  Moreover we can see from Figure 2. 4e that the location of the zones of horizon absence approximate well the reality (Figure 2. 4a). The validation results were conclusive:

·                           MG kriging was the least systematically biased of the three continuous estimators.

·                           OK, MG and IK should NOT be used to classify zeroes and non-zeroes.

·                           Kriging the conditional probability of horizon presence, thresholding the result and using this to mask the kriged thickness reduced the estimation bias and delineated zones of horizon absence.

·                           In the absence of other information the classification criterion should be reproduction of the sample proportions. In the absence of prediction bias, this can be achieved by applying a threshold to the categorical indicator estimates at a value zk = 0.5.

Figure 2. 4. Simulated transect, (a) true simulated thickness, TSIM, (b) T*SIM OK, (c) T*SIM MG, (d) T*SIM IK, (e) T*SIM OPT. The points represent the position of the data used to krige.

In-situ prediction using the ‘optimal’ method

The optimal two-step technique was applied to predict the lateral extent and local horizon thickness for the in-situ data. Each of the alternative methods described earlier were applied to the in-situ data, to validate the conclusions drawn from the tests using simulated data. These tests completely confirmed the previous conclusions. MG kriging was used to estimate the horizon thickness and categorical IK (cIK) to predict the conditional probability of horizon presence. Thresholding the cIK estimates at a value of 0.5 did not reproduce the sample lateral extent. Instead the categorical indicators of horizon presence were thresholded at quantile values equal to 1 – sample proportion, before multiplying the two predictions. Maps of the predictions (Figure 2. 5a-c) show in black the regions where no horizon is expected to exist. Predicted horizon thickness is typically greatest towards the centre of zones of presence. However, because the mask is estimated independently from the thickness estimate, boundaries of the zones of presence exist side-by-side thick horizon estimates.

variable (units:metres)




Std. Dev.

Extent (%)

Diffuse Carbonates






DT (sample)












Nodular Carbonates






NT (sample)


















GT (sample)












Table 2. 3. In-situ sample and optimal estimate statistics.

The small estimate standard deviations (Table 2. 3) show that kriging has significantly smoothed the thickness estimate histogram. Also, the small values of the mean thickness estimates are evidence of the underestimation of large values and the overestimation of zero values. This effect was best illustrated by scatter plots of validation estimates for each sample location and by the nodular carbonate thickness estimate. Nodular carbonate thickness had the largest relative nugget effect, testament to the number of isolated large thickness values surrounded by zero values, which were subsequently severely underestimated. The smoothing effect was also severe for gypsum thickness because the sample was dominated by zero values (sample proportion 20 %).

The prediction maps confirmed our initial hypotheses of the spatial distribution of carbonate and gypsum horizons. The diffuse and nodular carbonate horizons are more common than that of gypsum. The thickest nodular and diffuse carbonate horizons are collocated to the northeast and west of the study area in Red soils lying on raised levees. The diffuse and nodular carbonate horizon extend into the region dominated by the Cotton soil, to the south of the study area. The horizons in these areas is thin and considerable uncertainty is associated with the predicted thickness values there.

Figure 2. 5. (a ‑ c) Horizon thickness estimate (units: meters) (d ‑ f) MG kriging estimation variance and (g ‑ i) sGs thickness E-type mean for each horizon type.

Kriging and simulation measures of uncertainty

Figure 2. 5d-i (above) and Figure 2. 6a-f (below) show the MG estimation standard deviation and simulation standard deviation for each horizon type and each simulation method. All three have been presented to show how they relate to one another, their shared information content.

Figure 2. 6. (a-c) sGs indicator simulation standard deviation, (d-f) sis simulation standard deviation for each horizon type.

The standard deviation maps derived from the simulated realisations exhibit patterns that are specific to each variable because they reflect the likelihood of finding particular spatial patterns of values. Hence, in zones where the sample values were homogeneously small (zero) or large the spatial uncertainty is low. For zones where the sample values are both small and large, the predicted uncertainties are larger. Exactly the same information could have been provided by multiplying the estimation variance map by the predicted horizon thickness value at each location, because the single point cumulative conditional distribution function (ccdf) is fully specified by the mean and variance for a multiGaussian random variable. Using these maps we can delineate regions based on two criteria, (a) whether we expect there to be a horizon present, using the kriging / simulation conditional probability and (b) on the basis of the degree of certainty in the estimate, via the simulation standard deviation. The individual realisations could also be used as inputs to spatial models.

Figure 2. 7. Scatterplots of simulation post-processing results for the diffuse carbonate horizon – (a) sis mean thickness vs. sis standard deviation, (b) sGs mean thickness vs. sGs standard deviation, (c) sis mean vs. sis standard deviation, (d) sGs indicator coded mean vs. sGs indicator coded standard deviation.

Despite the similarities between the results provided by each simulation method subtle differences exist. For example, the E-type simulated mean lateral extents of the sis realisations approximated well the sample estimates. However the means of the indicator coded sGs realisations, for each horizon type, were consistently 2-5 % larger than the corresponding sample estimates (Table 2. 4). Also the 5 % and 95 % quantile values were systematically larger than those provided by sis. However the simulated range of uncertainty was systematically smaller when using sGs to simulate lateral extent (-4.73 to 7.05 % of the total area).

These differences are caused not just by differences in the variogram models used but importantly by the influence of the sample values when sequentially and randomly drawing values using the simulation model of spatial uncertainty (equation 9). With reference to diffuse carbonates, but representative of all three horizon types Figure 2. 7 a and b show the simulation standard deviation vs. the simulated mean thickness.

There is a non-linear correlation between the sGs thickness simulation mean and standard deviation. In contrast, the sis indicator conditional probability shows no correlation with the estimated value of thickness. The value of thickness has no role to play when simulating an indicator (and when estimating the indicator variogram) and therefore has absolutely no effect on the conditional probability and the standard deviation generated by post-processing.


Unts : km2

Diffuse Carbonates


Nodular Carbonates




Sample estimate







sGs-sis (%)



sGs-sis (%)



sGs-sis (%)

5 %










50 %










95 %




















Table 2. 4. Comparison of the 5, 50 and 95 % quantiles of simulated lateral extent for sis and sGs methods. The difference value SGS ‑ SIS is expressed as a percentage of the total area of the region.

This is not the case when indicators are provided by simulating thickness. Value plays a role when simulating the local thickness and thus the lateral extent and therefore has the knock-on effect of modifying the estimated conditional probabilities and standard deviations of the indicator variable defined thereon.

A comparison of Figure 2. 7 c and d shows that the maximum conditional probabilities are greater for the sGs method. Points with the largest conditional probabilities, with values close to one, have the lowest standard deviations almost equal to zero.

The points with the largest conditional probabilities are in regions of homogeneous thick horizons. This reflects the degree of certainty of there being a horizon present. Yet this effect is not equally expressed for zones of horizon absence.

The conditional probabilities are not equivalent and the simulation local standard deviations not as low (conditional probability @ 0.2, standard deviation @ 0.4) as might be expected for regions that are just as homogeneous as their counterparts. Lower uncertainty in regions of homogeneous thick horizon presence seems logical, however this must be equalled by a similar reduction in uncertainty for regions of homogeneous horizon absence. The sis minimum conditional probability, in regions of homogeneous horizon absence, is associated with a lower standard deviation than that provided by the sGs result (conditional probability @ 0.1, standard deviation @ 0.3). This is a value equivalent to that of equally homogeneous regions of horizon presence (conditional probability @ 0.9, standard deviation @ 0.3).


Simulation provides a useful tool with which to validate the choice of prediction method, in particular when the sample data are sparse. From statistics calculated on the sample dataset, an exhaustive reality can be generated that shares the same spatial properties as the sample dataset. Various prediction methods and hypotheses can be tested against the simulated reality. The tests validated our choice of a two-step kriging method for the estimation of the lateral and vertical extent of soil horizons. If the horizon attribute is not present everywhere, it is not continuous, and the prediction quality benefits from a prior delineation of the regions where the attribute is most likely to be observed. Nevertheless, the in-situ predicted thickness values were systematically biased and the estimates of both thickness and lateral extent, subject to much uncertainty given the scarcity of data.

By simulating the lateral and vertical extent we were able to generate simulation standard deviation maps that provided important measures describing the spatial uncertainty of the aerial extent of the horizons. The optimal kriging prediction, using a multiGaussian thickness estimate and estimation variance could provide similar local uncertainty information, but simulation realisations can be used to investigate spatial uncertainty and therefore represent a richer source of possible information for future use. Objections to their use based on computer memory and processing time can be rejected.

By modelling the spatial uncertainty both