Analysis
of Proximal and Remotely Sensed Spectral Data for the Prediction of Soil
Properties
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
Biography
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)
INSEAD
Boulevard de
Constance
Fontainebleau
77300
France
Tel: 33 (0) 1 60 72 44 56
Fax: 33 (0) 1 60 74 55
64
Email: Benjamin.Warr@insead.edu
Website: http://benjamin.warr.free.fr/
Acknowledgements
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.
Preface
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 everlarger spatial scales, and improved spatiotemporal 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’ handeye 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 nondestructive 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 handeye 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
Characterising soil variation.
1
The study area and survey data. 7
Ordinary and indicator kriging. 13
MultiGaussian Random Functions. 17
Conditioning multiGaussian Random Functions. 18
Quantitative classification. 25
Interpolation and simulation. 27
Structural analysis and interpolation. 27
Predicting and simulating the geometry of soil
horizons. 38
Kriging continuous attributes with positive skew.. 40
Kriging categorical attributes to delineate zones of presence. 41
Structural analyses of horizon attributes. 42
Comparing prediction methods. 45
Simulating the insitu variables to characterize spatial uncertainty. 45
Comparing prediction bias using a simulated dataset 46
Insitu prediction using the ‘optimal’ method. 48
Kriging and simulation measures of uncertainty. 50
Soil
colour and spectral reflectance.
55
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
Psychophysical colour spaces: Hue, saturation, lightness. 70
Spectral reflectance measurements. 71
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
Multivariate spatial predictions of soil
colour. 88
Issues when kriging multivariate colour representations. 89
Nonuniformity of scale and the property of additivity. 90
Nonuniformity of scale and nonlinear multivariate correlation. 90
Compositional data: Constant sum constraints and linear dependence. 91
Compositional kriging for constrained predictions. 91
Remotely sensed observations of the soil
106
Sensor systems and platforms. 109
Multispectral optical sensors. 109
Synthetic Aperture Radar (SAR)  Imaging and Applications. 114
Sensor noise and image statistics. 117
Fourier analysis using the 2D Fast Fourier Transform (FFT). 118
Kriging the intrinsic component 121
Kriging a secondorder stationary component 122
Study 1: Comparing Fourier analysis and Factorial Kriging analyses. 122
Applying Factorial Kriging. 122
Structural Analysis of Multitemporal SAR Imagery. 122
Structural Analysis of a SPOTXS Image. 123
Factorial Kriging of Multitemporal SAR Imagery. 124
Factorial Kriging of SPOTXS 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
Structural analysis and variogram fitting. 133
Cleaning Thresholded Ratio Images Results. 138
Study 3: Applying FK to improve soil property – image derived variable correlations. 142
Structural Analysis of Landsat TM Imagery. 143
Factorial Kriging of Landsat TM Imagery. 144
Calibrated soil reflectance models for the
prediction of soil properties. 144
Study 1: Partial least squares regression. 144
Introduction: Problems with ordinary least squares regression. 144
Collinearity and Overfitting. 144
The Partial Least Squares (PLS1) algorithm.. 144
Methods: Calibrating the PLS model 144
Study 2: Multivariate spatial prediction methods. 144
Introduction: maximising the use of available data. 144
Regression and residual kriging. 144
Fuzzy Classification Mapping. 144
Kriging with an external drift 144
Methods: scenarios of data availability. 144
Table of Figures
Figure 1. 2. Bourke Irrigation Area. Origin (bottom left): 37500 E, 6655000 N. 9
Figure 1. 4. Natural eucalyptus woodland vegetation typical cover for the Red soil. 20
Figure 1. 16. Ordinary kriged (OK) map for (a) sand (b) silt and (c) clay content. 34
Chapter
2
Chapter
3
Figure 3. 1. Diagrammatic representation of an electric and magnetic field. 59
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. 8. (a) Mapping of XYZ tristimulus values to (b) an x y Chromaticity Diagram. 68
Figure 3. 12. Schematic diagram of a double beam spectrometer. 71
Chapter
4
Figure 4. 6. CIE L^{*} a b estimates vs. Munsell estimates (a) Hue, (b) Value and (c) Chroma 100
Chapter
5
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. 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. 19. July Longrange Components  20, 50, 100m Grid. 131
Figure 5. 20. Ratio Image Variogram Models and Filtered Image Models. 134
Figure 5. 21. Filtered Ratio Image qq Plots (black), 1:1 bisector (grey ). 135
Figure 5. 22. JulyAugust Ratio Images. 136
Figure 5. 23. Effects of Connected Components Removal 139
Figure 5. 24. Landsat TM anisotropic variogram models for (a) band 1, (b) band 2. 143
Chapter
6
Figure 6. 3. PLS1 weighting vectors 1 to 4. 144
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 decisionmaking 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 costeffective 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 sitespecific 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 costeffective 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 timeconsuming 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 costtime 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 ‘nondestructive’ measurement technique. Proximal sensing (PS) and remote sensing (RS) technologies also provide nondestructive 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 threeway tradeoff 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 semiarid 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 spatialspectralcoverage tradeoffs, 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 highresolution 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 (PalaciosOrueta 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 unmixing 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 evergreater quantities of insitu 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 highdimensionality 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 highdimension spectrometer data is replaced by a lowdimension 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.
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 HapexSahel 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 northwest 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 900km^{2} 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 (010 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 ‘micropippette 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 (~23cm). Horizon thickness data is a censored variable. From the same cores soil colour measurements were made of horizons are at four depths (020cm, 3040cm, 7080cm and 90100cm) 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. 
MethodsThe rest of this introductory chapter describes the theoretical basis and application of ‘standard’ nonspatial 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 ClassificationTraditionally 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 soillandscape 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 decisionmaking 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, m_{ij} is the membership of an individual i in class j and d^{2} 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 cutoff 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 soillandscape models are developed to provide predictive information about the spatial distribution of soil classes D_{k}, over the study domain D. Laboratory analyses provide statistics that describe the average attribute values z_{k}(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 R_{E}V’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).
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 zdata into normal score yvalues,
Equation 1.
4
_{}
where f is a normal score transform that relates the pquantiles z_{p} and y_{p}_{ }of the two distributions (Rivoirard, 1994). The variogram of the y (x_{a}) can then be recalculated using Equation 1. 3. Typically the y(x_{a}) variogram model g_{Y}(h), has a smaller relative nugget effect and greater spatial continuity over short lag distances than its g_{Z}(h) equivalent, because of the reduction of influence of extreme values.
Indicator variables,
I(x_{a} ;z_{k}), are generated by applying an
indicator transform to a continuous variable,
z(x_{a}), at selected thresholds, z_{k} which are
often chosen to be equal to the decile values of the cumulative
distribution,
Equation 1. 5
_{}
The indicatortransformed 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 z_{k} = 0. The
result is a categorical indicator variable equal to one if the horizon
s_{k}, is present and zero if
not,
Equation 1. 6 _{}_{}
The experimental variogram of an (categorical) indicator
variable I(x ; z_{k}) 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 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 x_{a}, using the surrounding n data {z(x_{a}), a = 1,…n}. Kriging is optimal for the estimation of quasistationary 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 l_{a}(x) is the weight assigned to the data z(x_{a}) and m(x_{a}) is the mean (Goovaerts, 1997). Simple kriging (SK) requires that the mean m (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 l_{a} 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 nonbias condition a zero sum expression including a Lagrangian L(u) parameter 2m_{OK} 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
x_{0}. 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 (1980ac) 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) Nonbiased, 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
equiprobable 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 multiGaussian 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 singlepoint ccdfs _{}. Each estimate when considered individually is best (in the leastsquares 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 equiprobable realisations. In contrast to kriging, simulation randomly selects a value z(x¢_{j}), j = 1,…, N, at location x, from the Npoint 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 nonsequential methods. A brief description of the concepts underlying the most commonly used sequential methods follows.
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
X_{n}, which have Gaussian statistical distributions
~N (0,1) have favourable properties of stability through addition
((X_{1} +
X_{2 }= X_{3} is N (0,1)) and
through linear combination (lX_{1} + lX_{2 }= X_{3 } is N (0,1)), and
this can be generalised,
Equation 1. 15
_{}.
The central limit theorem infers that as _{}, any combination of independent random variables
X_{1}, X_{2},…, X_{n}, 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 ~N (0,Var(X)). The multiGaussian RF also
has favourable spatial properties. The twopoint 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.
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 nonstationary, because SK is an exact interpolator. To simulate independent stationary residuals the following algorithm is used.
1)
Nonconditional simulation of RF,
Y(x) giving Y_{S}(x) at all
locations.
2)
Krige
using the Y_{S}(x) that coincide with data
locations as new samples to estimate a new variable
Y_{S}(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 Y_{SC}(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 Y_{i}(x) within a neighborhood surrounding a location, to give Y(x)^{K} at all locations.
3) Simulate a variable ~N (0,1) and multiply this by s_{K}, and add to either the conditioning data value or Y(x)^{K}.
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 pseudoprobabilities 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 
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 010cm 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.
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 lowlying 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 coarsesandy 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 lightcoloured (5 YR/V/C), finetextured ‘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 lightcoloured 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 semiarid 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 buildup 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 MurrayDarling 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 lowlying 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 010cm
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 largegrained 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 onfarm flood and drain irrigation or when the MurrayDarling 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 oneway 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.
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. 
Table
1. 1.
Summary of oneway 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, d^{2}(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. Insitu soil type allocation based on expert knowledge is a useful and cheap practice. The eyebrain 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.
Histograms for each of the insitu 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. 12hj) 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).
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 nondirectional 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 shortrange 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 shortrange 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 insitu 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 postprocessed 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. 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 shortrange 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 (ultisolsalfisols) corresponding to the southwestnortheast 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)
Fe^{* }_{OK} 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 shortrange 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 northeast 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 MurrayDarling 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.
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 (CO_{2}). 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 ZiTong, 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
multiGaussian 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 preexisting 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 insitu data and (ii) then to investigate the uncertainty associated with the prediction of the soil horizons lateral and vertical extent by geostatistical simulation.
The
ordinary kriging (OK) system is derived under the assumption that the mean
m(x_{a}) is constant but unknown within a local neighbourhood
W(x). Therefore the local mean must be estimated for each
z(x_{a}). 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 multiGaussian kriging (MG) and indicator
kriging (IK). Both require a nonlinear (back) transform of the data prior
to and after kriging. The first step of multiGaussian kriging is the
normalisation of the z(x_{a}) to generate a variable with a Gaussian distribution. This
facilitates inference of the mean when kriging. The kriging estimates
Y^{*}_{SK} (x_{a}) are then backtransformed 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 z_{k} 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
highfrequency shortrange 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 nonzero zones based on
the predicted values. Indeed if any cutoff 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.
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 s_{k} 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 viceversa. 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.
A
sensible approach to the estimation of local horizon thickness should involve
twosteps; 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 insitu
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 crossvalidation, 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 nonparametric sis algorithms are discussed.
The insitu variables used were presence (or absence) of: 1) diffuse CO_{3}^{2+} (D), 2) nodular CO_{3}^{2+} (N) and 3) gypsum (G). We refer to the categorical indicator variables as I(x ; D) for the diffuse CO_{3}^{2+}, I(x ; N) for the nodular CO_{3}^{2+} and I(x ; G) for gypsum. The horizon thickness values are D_{T}, or N_{T,} or G_{T}, 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 
Table 2. 1. Variogram parameters of insitu horizon indicators (D, N, G), the normal score transforms of insitu horizon thickness (D_{Tg}, N_{Tg}, G_{Tg}) and the simulated transect thickness T_{SIM} and lateral extent P_{SIM} (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 colocated with diffuse carbonates (r = 0.45). Their spatial codependence 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 multiGaussian 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 cutoff 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 normalscore transformed insitu
thickness variable, was then used to simulate a variable
T_{SIM g} to 500
points along the transect using the 30 samples to constrain the sequential
Gaussian algorithm. The Gaussian simulated values were backtransformed to
the raw space of the insitu diffuse carbonate horizon thickness. A
single realisation T_{SIM}, 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 multiGaussian kriging. Thickness
T_{SIM}, was also indicatortransformed 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 cutoffs applied to diffuse horizon
thickness, z_{k} (a) z_{1} = 0.06 (b)
z_{2} = 0.15 (c) z_{3} = 0.32 (d)
z_{4} = 0.51 (e) z_{5} = 0.67 (f)
z_{6} = 0.83 (g) z_{7} = 1.13 (h)
z_{8} = 1.27 (i)
z_{9} = 1.36.
(REMOVED) Finally a categorical indicator was defined to create
the simulated horizon lateral extent P_{SIM}, 
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 setaside for validation purposes.
In a
first step, OK, MG and IK were used to provide estimates of simulated horizon
thickness, T_{SIM}. 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 nonzeroes.In a second step, the simulated horizon lateral extent
P_{SIM}, 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(x ; P_{SIM})]^{*},
Equation 2. 5
_{}
The
improvement in the quality the predicted thickness T_{SIM}
obtained by using the prior prediction of the probability of horizon presence
P_{SIM} to delineate the zones where the thickness value is
valid, was compared by the predictions against the true the exhaustive simulated
dataset.
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 (Etype 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 postprocessed to find the
Etype uncertainty measures.
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 (be) are all smoother than reality (Figure 2. 4a). Both the prediction transects and error histograms
(Figure
2. 4ad) 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. 
 
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 nonzero 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 P_{SIM} and local thickness T_{SIM} were severely overestimated 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)
Z_{SIM} – Z^{*}_{SIM OK}, (b)
Z_{SIM} – Z^{*}_{SIM IK}, (c)
Z_{SIM} – Z^{*}_{SIM MG}, (d)
Z_{SIM} – 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
P_{SIM}_{.
}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 nonzeroes. · 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
z_{k} = 0.5. Figure
2. 4.
Simulated transect, (a) true simulated thickness, T_{SIM}, (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. 

The optimal twostep technique was applied to predict the lateral extent and local horizon thickness for the insitu data. Each of the alternative methods described earlier were applied to the insitu 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. 5ac) 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 sidebyside thick horizon estimates.
Table
2. 3.
Insitu
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 Etype mean for each horizon
type.
Figure 2. 5di (above) and
Figure 2. 6af (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. (ac) sGs indicator simulation standard
deviation, (df) 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 postprocessing 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 Etype 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 25 % 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 nonlinear 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
postprocessing.  
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 knockon 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 twostep 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 insitu 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