Statistical Methods for the Analysis of Soil Spatial and Temporal Variability

Ahmed Douaik1, Marc van Meirvenne2 and Tibor Tóth3,4 1Research Unit on Environment and Conservation of Natural Resources, National Institute of Agricultural Research (INRA), Rabat, 2Department of Soil Management, Ghent University, Ghent, 3Department of Soil Science, Research Institute for Soil Science and Agricultural Chemistry, Budapest, 4Rural, Water and Ecosystem Resources Unit, Institute for Environment and Sustainability, European Commission – Joint Research Center, Ispra 1Morocco 2Belgium 3Hungary 4Italy


Introduction
Soils vary in space and time due to the combined effects of physical, chemical, and biological processes that operate at different scales and with different intensities.This variability is caused by soil forming factors, i.e., climate, parent material, time, topography, and vegetation (Jenny, 1980), farming practices (fertilization, irrigation, tillage, etc.), erosion, etc.Each cause may operate independently or in combination with other factors and over a wide range of spatial and/or temporal scales.Spatial scales can manifest at microscale to watershed and more, while temporal scale can span from seconds to decades and longer.These processes and causes create a pattern of nested variability or soil heterogeneity.This scale-dependency means that soil properties may display spatial and/or temporal patterns only over certain distances and not over others.Characterizing spatial variability of soil properties and crop parameters by inventorying them is needed to evaluate the effectiveness of their management.Also, characterizing their temporal variability by monitoring them is important in order to compare different management systems for sustainability and environmental quality.In addition, an understanding of the spatial and/or temporal variability can provide a framework for developing effective sampling schemes for future site management and efficient experimental designs for research approaches.Due to the spatial and/or temporal variability of soil properties, numerous samples need to be taken and the measurements need to be repeated as conditions change or to determine if they are changing.Soil properties are traditionally determined in laboratory from samples taken from field.However, this approach, although more precise, is more time consuming and expensive.Alternatively, soil properties can be easily, quickly, cheaply, and more intensively measured, even indirectly, using terrain attributes derived from a DEM and using emerging soil-landscape mapping technologies like GPS, geophysical techniques (EMI, GPR, TDR, etc.), soil spectroscopy and remote sensors, etc.Even with these auxiliary or proxy data sources, a full inventory of soil resources and properties is non-realistic and impossible.Therefore, soil can be measured at only limited space locations and time instants.Soil properties at non-sampled locations and/or times should be predicted.These predictions and related uncertainties can be obtained by using statistical models for soil spatial and/or temporal variability.In this way, soil variability can be considered as the combination of a systematic (deterministic) and a random (stochastic) components.The former represents the gradual or distinct changes (trends) while the latter involves differences not related to a known cause and those not discerned by the nature or the scale of investigation.The choice of the adequate model depends on the behaviour of the phenomenon of interest (a priori knowledge), the objective of the analysis, the nature of the available data, and, from a practical point of view, the existence of good statistical methods and adequate software packages.Excellent review papers on modelling spatial and temporal variability were published a decade ago (Heuvelink & Webster, 2001;Kyriakidis & Journel, 1999).The latter review focused only on geostatistical methods for modelling particularly the joint space-time variability while the former discussed methods for modelling purely spatial, purely temporal, and joint space-time variability.Our chapter is basically a review of statistical methods that have been used in order to describe and model the spatial and temporal variability of soil resources and properties.Once this variability is characterized and modelled, it can be advantageously used for different practical applications like optimal sampling of soil resources, efficient designing of experiments, interpolation and mapping of soil properties, etc.Therefore, in this chapter, we will present the main statistical methods with a short description of each of them, discuss eventually the possible links between some of them, identify their assumptions, merits and limitations, stress some of their practical applications, and list some relevant bibliographic references.The chapter fills the gap from the Heuvelink & Webster (2001) paper by considering more statistical methods like temporal stability (Douaik et al, 2006;Vachaud et al, 1985) and dynamic spatial variation (Douaik et al, 2007;Lesch et al, 1998), updates the presentation of the methods from the Kyriakidis & Journel (1999) paper by taking into account the research work done, and especially the last developments in pedometrics, during the past decade and adds some examples of practical applications like mapping using the Bayesian maximum entropy (Christakos, 2000;Douaik et al, 2005;2009).After presenting some classical statistical methods (Section 2), we provide an overview of the geostatistical methods including structural analysis (Section 3) and prediction (Section 4).Then, we present the Bayesian maximum entropy, another method of prediction (Section 5).Finally, we discuss the different methods and conclude with some suggestions for choosing among them (Section 6).

Classical statistical methods
Classical statistical methods are basically aspatial and atemporal methods where the spatial and temporal coordinates are ignored.The different approaches are: coefficient of variation (2.1), variance component analysis (2.2), temporal stability (2.3), and dynamic spatial variation (2.4).

Coefficient of variation
It is assumed that the observed value for a soil attribute (z) at any location (s) at a given time instant (t) is: with the population mean and (s,t) random, spatially and temporally uncorrelated error assumed normally distributed with zero mean and variance 2 .The population mean is estimated from a sample: with n: the sample size and the variance, 2 , is estimated by: ( ) Its root square, σ , is the standard deviation.The coefficient of variation is defined as: The coefficient of variation can be calculated for the temporal probability distribution function (PDF) at each space location as well as for the spatial PDF at each time instant.As an example, Van Wesenbeeck & Kachanoski (1988) studied the spatial distribution (from 100 locations) of the surface soil water content throughout the growing season (65 measurement dates).They concluded that there exists a temporal dependence of the spatial variability and a spatial dependence of the temporal variability.They also suggested that random sampling of the locations should be avoided.James et al. (2003) defined the temporal variability of soil moisture content as the coefficient of variation at each site across the whole year (52 time instants) and used an analysis of variance (ANOVA) to test for differences between habitats (three sites) in the temporal variability of soil moisture content.They found that soil moisture varied significantly with time but not between habitats while the interaction between time and habitat was significant which reflects differences between habitats in temporal pattern of soil moisture.The coefficients of variation were significantly different between habitats.Other examples of the use of the coefficient of variation for the assessment of the spatial and temporal variabilities can be found in Ehrenfeld et al. (1997) and Guo et al. (2002).
The standard deviation and, hence, the coefficient of variation does not incorporate spatial and/or temporal information and therefore doesn't provide a representation of the nature of the spatial and/or temporal behaviour of a given soil property.It is considered as an indicator of global variability as opposed to local variability (Carter & Pearen, 1985), the latter is identified using correlogram, covariance function, or variogram which will be presented later in this chapter (sections 3.3 to 3.6).

Variance component analysis
In the variance component analysis approach, the spatial and temporal contributions are considered as factors with random effects in the ANOVA model.In this way the different components of the total variance can be determined.Campbell et al. (1989) adopted the following model, for the log concentrations of a given solute (log c) in subsample l (3) from sampling area k (2) at site j (2) and at time i (6 dates over 20 days): with µ: the grand mean; T: the time effect; S: the site effect; A: the sampling area effect; (A.U): the subsampling within sampling area effect, and : an error term.The temporal and spatial variances were compared to the sampling area variance while the latter was compared to the subsampling within sampling area variance.They found that the temporal variation was far less significant than the variation between the two sites.Also they found that the variation between areas a few meters apart was greater than that between the two sites even though the texture and the management of the sites were quiet different.Other examples of the evaluation of the temporal and spatial variabilities based on the concept of variance component analysis can be found in Van Es (1993) and Van Es et al. (1999).
The variance component analysis is directly related to the coefficient of variation presented in the precedent section.The variance component estimates for each variability source are computed, and then the square root of these estimates is divided by the mean value.The output is exactly the definition of the coefficient of variation.Thus, the variance component analysis is suffering from the same drawbacks as the coefficient of variation approach.

Temporal stability
The concept of temporal stability was first introduced by Vachaud et al. (1985).It is defined as the time invariant association between spatial location and statistical parameters of soil properties.They distinguished mainly two approaches: relative differences and Spearman rank correlation.

Relative differences
Let z ij be the observed value of a soil property at location i (i=1, …, n) and time j (j=1, …, m).The relative differences ij are defined as: with j z : the spatial average for the time instant j.The temporal average of ij is: and its corresponding temporal standard deviation is defined as: A zero value for i δ indicates that the temporal average j z represents the average value over the whole study area at any time.The field average value is overestimated if i δ > 0 while it is underestimated when i δ < 0. A more time stable location will be indicated by a www.intechopen.comsmall value of ( ) i σδ whereas a high value of the latter is an indication of a less time stable location.The i δ values can be plotted graphically against their rank with the corresponding temporal standard deviations.Based on the relative differences, Kachanoski & De Jong (1988) showed that a useful test for temporal stability is the Pearson correlation coefficient between soil property values measured at two consecutive time instants.They also showed that temporal stability exists if the relative differences remain constant between two time instants: 21 ii δδ = .This equation implies that: 1 z and 2 z being spatial averages for the time instants j = 1 and 2, respectively.This equation establishes a linear relationship between the soil property at two different times with an intercept equal to zero and a slope equal to the ratio of the mean value observed at the second time frame to the mean value observed at the first time frame.Consequently, a good test for the temporal stability is the correlation between soil property values measured at consecutive time frames (Kachanoski & De Jong, 1988).Furthermore, the former equation implies that the regression between relative differences from two consecutive time frames should have a zero intercept and a unity slope.This is another way to check the existence of temporal stability of the spatial pattern of a variable.However, the regression between z i2 and z i1 is, in general, of the following form: Using equations ( 9) and ( 10) we expect that, if temporal stability exists, the following conditions are verified: So, another way to check the temporal stability is to test for a zero intercept and a unity slope of the regression between two consecutive time instants.If there is a constant increase or decrease in the soil property, w, in all the locations between two time frames, we expect from Eq. ( 10) that (Kachanoski & De Jong, 1988): Equation ( 12) indicates that if the slope of the regression between values at two consecutive time instants is not significantly different from one, then the regression intercept represents the constant change that has occurred between the two time instants.Based on whether the intercept (b) is significantly different from 0, the slope (a) is significantly different from 1, or both, four possible scenarios can be distinguished (Douaik et al, 2007;Grant et al., 2004): -Scenario 1 (a = 1 and b = 0): there is no net change between the two sampling dates -Scenario 2 (a ≠ 1 and b = 0): the mean soil property level has changed and there is a static (uniform) change in the soil property spatial pattern -Scenario 3 (a = 1 and b ≠ 0): the mean soil property level has not changed and the change in the soil property spatial pattern is dynamic (non-uniform) but, when averaged, this change is not significant -Scenario 4 (a ≠ 1 and b ≠ 0): the mean soil property level has changed and the change in the soil property spatial pattern is dynamic (non-uniform).These refinements of the concept of temporal stability were used by, among others, Da Silva et al. (2001), Douaik (2005), Douaik et al (2006;2007), Ehrenfeldt et al. (1997)), Goovaerts & Chiang (1993), Petrone et al. (2004), Van Wesenbeeck et al. (1988).The concept of relative differences is based on the same statistical parameters (mean and standard deviation) as the two precedent approaches.However, it uses the standardized data (relative differences) for each time instant instead of the original data and the mean and standard deviation are computed for each space location.This concept is based, like the two precedent approaches, on the Gaussian distribution of data and the absence of autocorrelation.Both assumptions must be checked before the use of relative differences in the study of temporal stability.The normality can be checked using the standard statistical tests and the independence can be verified using either a test for autocorrelation like the Moran test or by computing the variogram and determining its range.Normality can be gained via data transformation and independence by the selection of sample locations separated by distances larger than the autocorrelation range.For non normally distributed data, but still non autocorrelated, the Spearman rank correlation can be used.The latter is presented in the next section.

Spearman rank correlation
It refers to the tendency of a soil property, measured at different locations in space, to maintain their relative ranking over time.It is defined as: with R ij and R ik the ranks of z ij observed at location i on time instants j and k, respectively.A value of this coefficient equal to one indicates a perfect temporal stability between time instants j and k, and thus identity of ranks for any location whereas a lack of temporal stability implies that r s = 0.Many researchers used this concept to evaluate the temporal stability of different soil properties, among them, Campbell et al. (2001), Douaik (2005), Farley & Fitter (1999), Guo et al. (2002), Heathman et al. (2009), Hu et al (2011), Reichardt et al. (1993), Si (2002), andWendroth et al. (1999).The Spearman rank correlation coefficient is a non parametric statistics therefore it can be used even for non normally distributed data.However, in the case of autocorrelation, its use is not recommended and it will be cautious to use statistical methods which consider the spatial and temporal nature of the given soil property.In addition, even for independent data, the Spearman rank correlation inform us only about the relative ranking of space locations between two time instants and do not provide a quantification of the temporal change.This quantification can be obtained using either the paired-t test for comparing mean values or the dynamic spatial variation approach.The former is presented in what follows while the latter is discussed in section 2.5.

The
with d the mean of the differences, d i , between the observations for the second and the first time frames: ( ) s d is the standard deviation of the differences: ( ) and n the number of the differences, i.e. the number of locations.
The test statistic t (equation 14) has (n-1) degrees of freedom.The calculated statistic is compared to a tabulated t value with the same degrees of freedom and a probability of error of type I (taken generally equal to 5%).Equivalently, the corresponding probability to the test statistic t is compared to 5% to check if the difference of mean values is significant or not.As an example of application of the paired-t test, Kenny et al. (2002) checked the temporal trend of the mean thickness of the Ap horizon.

Dynamic Spatial variation
The concept of temporal stability handles each soil property separately and, as described before (section 2.3), it can quantify the temporal change even without testing if this change is significant or not.Lesch et al. (1998) discussed another approach that can test for change in the mean value as well as for a dynamic spatial variation between two time instants.Dynamic spatial variation refers to the spatially variable change of a soil property between two time instants, i.e., the change is not the same between two time instants for each location.
To test for the existence of dynamic spatial variation, Lesch et al. (1998) identified different steps.First of all, the pairs of soil data obtained at a first time instant (y and z), from n space locations, are used to establish a calibration equation based on the regression model: The parameters of this regression equation (b 0 and b 1 ) are estimated using the ordinary least squares method.In a second step, the residuals of the model are tested for spatial independence, normality, homoscedasticity and the linearity of the relationship.The three last assumptions can be verified using classical regression tools.The assumption of spatial independence can be checked using the Moran test (Lesch et al., 1995).In the next step, a mixed ANOVA model is set: i=1, 2 are the two time instants and j=1, … , n are the locations.This model is equivalent to: ( ) with 021 dtt =−is the difference between mean values for the second and first time instants and j η independently and identically normally distributed with zero mean and variance equal to ( ) represents the observed values of the first soil property at the first time instant and at locations j (j=1, …, n) and z 2k represents the same property observed at the second time instant at the locations k (k=1, …, m), the difference between the two values is given by Eq. ( 19) with d 0 representing the shift in the average value between the two time instants and representing the dynamic spatial variation.These two quantities can be statistically tested: the first if it is different from zero and the variance of the second if it is higher than zero.Based on whether the first parameter (the mean difference, d 0 ) is significantly different from zero, the second parameter (spatial shift, 2 ) is significantly higher than zero, or both, four possible scenarios can be distinguished (Douaik et al., 2007) : -Scenario 1 (d 0 = 0 and 2 = 0): there is no net change with time in the soil property -Scenario 2 (d 0 ≠ 0 and 2 = 0): the mean soil property level has changed and there is a static (uniform) change in the soil property spatial pattern -Scenario 3 (d 0 = 0 and 2 > 0): the mean soil property level has not changed and the change in the soil property spatial pattern is dynamic (non-uniform); again, when averaged, this change is not significant -Scenario 4 (d 0 ≠ 0 and 2 > 0): the mean soil property level has changed and the change in the spatial pattern of the soil salinity is dynamic (non-uniform).An application of the approach of Lesch et al (1998) can be found in Douaik (2005), Douaik et al (2007) and Tóth et al (2002).The approach of Lesch et al (1998) is based on the determination of a linear regression between two soil properties; hence it assumes that the residuals from the regression model are normally distributed with equality of variances and spatially not autocorrelated.All these assumptions must be checked.If spatial autocorrelation is present in the residuals, the estimates of the regression parameters and their variances are biased and the statistical tests are no longer valid.In this special case, spatial regression and/or geostatistical approaches (sections 3 and 4) are recommended.

Geostatistics: Structural analysis
Structural analysis is concerned with the description and the modelling of the spatial and/or temporal variability using geostatistical tools like the covariance function, the correlogram, or the variogram.The analysis of space-time data requires the definition of a space-time continuum with a coordinate system and a measure of space-time distance and models and techniques which make possible the link between spatio-temporally distributed data (Christakos et al., 2002).The spatial coordinates are defined, in general, in two dimensions s = (s 1 , s 2 ) and the temporal coordinate t along the temporal axis such that the space-time coordinates are p = (s, t).Christakos et al. (2002) defined space-time random field (STRF) model as 'a mathematical construction that rigorously represents the distribution of natural phenomena across space and time.The STRF model provides scientifically meaningful representation, explanation, and prediction of these phenomena in uncertain environments.These uncertainties may be due to measurement errors, heterogeneous databases, erratic fluctuations in the space-time variation of the underlying process, and insufficient knowledge'.An STRF, X(p), is a collection of realizations χ of the space-time distribution of a natural variable.It can be viewed as a collection of correlated random variables x=[x 1 , …, x n ]' at the space-time points p=[p 1 , …, p n ]' and a realization of the STRF at these points is χ=[ χ 1 , …, χ n ]'.It assigns a probability that a realization χ in n dimensions will occur following the multivariate probability distribution function (PDF) of X(p) or, equivalently, its cumulative distribution function (CDF).The PDF and CDF represent a complete characterization of the STRF.A partial, but sufficient, characterization of the STRF is provided by its space-time statistical moments.The first order moment, which is the mean function, expresses trends or systematic structures in space-time.The moment of second order, the covariance function, expresses correlations and dependencies between two different points p and p'.In addition to these two moments, an STRF can be also characterized by its variogram.Christakos and Raghu (1996) reported that some forms of averaging over time or space are done before the analysis of the spatial or temporal random fields.In doing so, the joint space-time variability is not accounted for and there is some loss of information.Initially to consider both space and time variability, some authors resorted to some simplifications.For example, Bilonick (1985) and Egbert & Lettenmaeir (1986) decomposed their STRF in a purely spatial and a purely temporal components whereas Rodriguéz-Iturbe & Eagelson (1987) considered separable models where spatial variations were assumed to be independent of temporal variations.Three main approaches can be distinguished for the analysis of the geostatistical space-time data (Kyriakidis & Journel, 1999;Stein & Sterk, 1999;Stein et al., 1998): -Methods based on vectors of spatial random fields (SRF) for the case of many space locations and few time instants (Bogaert & Christakos, 1997;Papritz & Fluhler, 1994).It does not include the temporal dependence existing between observations and can predict only at the observed time instants.The spatial variability is modelled either through a separable spatial variogram for each time instant, or by a single spatial variogram considering time instants as replicates; -Methods based on vectors of time series (temporal random fields, TRF) for the case of long time series with few space locations (Rouhani & Wackernagel, 1990;Solow & Gorelick, 1986).This approach doesn't take into account the spatial dependence and it predicts only at the observed locations.In a similar way as above, independent TRF or TRF as replicates in space can be considered.-Methods using an STRF, thus considering the joint space-time variability; For the first approach, the focus of the analysis is on smooth interpolated soil attribute maps over specific time instants.The intention is to capture single instantaneous snapshots (a static picture) of the natural process.The objective can be also the comparison of the various www.intechopen.commaps or the detection of the temporal persistence or change in the spatial patterns (Goovaerts & Chiang, 1993;Van Meirvenne et al., 1996).The intention of the second approach is to capture a sequence of successive snapshots at single space locations, which give temporal profiles.Only the third approach includes both the spatial and temporal dependencies so the interpolation is more precise and can be done for unsampled time instants at unsampled space locations.The focus is here on video sequence of successive spatial pictures like a movie.As a fully satisfactory stochastic model should involve explicitly both spatial and temporal aspects, we present in the following sections the STRF concept and, subsequently, some special cases.The spatial and/or temporal variability can be described and modelled using either the covariance function or the variogram; although the latter is far more used in the geostatistical literature than the former and, therefore, only the variogram forms will be presented here.Also, in the case of two different soil properties or for the same soil property but observed at two different spatial locations or time instants, the joint variability can be characterized using the cross covariance function or the cross variogram.We present in section 3.3 the joint STRF while the spatial random field (SRF) is presented in section 3.4 and the temporal random field (TRF) in section 3.5.Finally, the space-time, the spatial, and the temporal cross variograms are discussed in section 3.6.

Joint Space-Time Random Field model
The space time variogram is defined as: It is estimated, from the observed data, by (Stein et al., 1998): x s h t Nh (21) with (,) i j x st and (,) s τ ++ ij xh t are pairs of observations with a spatial distance equal h at the time instant t j , the total number of such pairs is N j (h), and a temporal distance equal , the number of such pairs is N( ).Bilonick (1985) and Stein et al. (1998) are examples of application of the space-time variogram while examples of application of the space-time covariance function include Bennett (1975) for the analysis of population diffusion, Bell (1987) for the study of rainfall, Douaik (2005) and Douaik et al. (2004a) for the analysis of soil salinity.

Separate Spatial Random Field model
To each of the n t time instants corresponds a spatial random field (SRF) X(s, t j ), j=1, … , n t .The spatial variogram, corresponding to each of the above SRF, is defined as: It is estimated by: x st , ( ) ] of data separated by a spatial distance h, at the time instant t j .The spatial variogram (equation 22) is a special case of the space-time variogram (equation 20), i.e. when the temporal lag is zero meaning the same time instant: () The spatial variograms include the spatial dependence but the dependence in time is not accounted for.It is useful in the case of rich data in space and scarce data in time.The spatial variograms can be one, two or three dimensional.

Separate Temporal Random Field model
In a similar way to SRF, for a fixed space location s i , i=1, …, n s , corresponds the independent temporal random fields X(s i , t).
The temporal variogram at the location s i , i=1, …, n s , is: It is estimated as: The temporal dependence is included whereas the spatial dependence is not considered.Thus, it is useful if many data are collected in time and few in space.The temporal variograms are one dimensional.Application of temporal variograms in the soil literature is scarce and Cichota et al (2006) and Franzluebbers et al (2002) are examples of such application.

Space-time, spatial, and temporal cross variograms
The space-time cross variogram is defined, for two soil properties X and Y, as: It is estimated by: In the case of two SRF observed at two time instants t j and t j' , the spatial cross variogram is defined by: It is estimated by: -- τ Nhis defined as (,) τ Nh for a particular fixed temporal lag .By analogy with the spatial cross variogram, the temporal cross variogram for two TRF observed at two spatial locations s i and s i' is defined as: It is estimated by: is defined as (,) τ Nh for a particular fixed spatial lag h.Bishop & Lark (2006), Maas et al. (2010), and Motaghian & Mohammadi (2011) are examples of application of cross-variogram in the soil science.The definition of the spacetime, the spatial, and the temporal cross-variograms requires that data are observed at the same space locations and time instants, at the same space locations for two different time instants or at the same time instants for two different space locations, respectively.If it is not the case, the space-time, the spatial or the temporal pseudo cross-variogram (Papritz et al., 1993) needs to be used.Vanderlinden et al. (2006) and Zhang et al. (1999) are examples of application of pseudo cross-variogram and cokriging in soil science literature.

Geostatistics: Prediction
When a theoretical model is fitted to the experimental (observed) space-time covariance function or variogram, it becomes possible to tackle the problem of predicting attributes at unsampled space locations and/or time instants.Consider the problem of predicting the value of a continuous attribute x at any unsampled space location s 0 and time instant t 0 , x(s 0 , t 0 ), using only the x-data available over the space region S and time domain T, say, n data { } ( , ), 1,..., ; 1,..., Kriging is a family of generalized least squares regression algorithms that allows to predict x(s 0 , t 0 ).The latter is considered, in geostatistics, as a realization of the STRF X(s 0 , t 0 ).The problem of predicting x(s 0 , t 0 ) is defined as (Rouhani & Myers, 1990) with ij λ , the kriging weights for the time instant t j and the space location s ij .This is the general form and there are some simplifications which are related to the kind of STRF considered.We present first the simple cases where we consider only SRF or TRF, and then we present the more general case for the single STRF.

Spatial or temporal kriging
In case of spatial or temporal random field (sections 3.4 and 3.5) we saw that the dependence in data is modelled via a spatial or temporal covariance function or variogram.For example, in the case of an SRF, we are dealing in the space domain and the different kriging algorithms defined for SRF can be applied.Let X(s, t 0 ) denotes the SRF at the time instant t 0 .The predictor X(s 0 , t 0 ), at the space location s 0 , is a linear combination of the n s observations x 1 , …, x ns available at the time instant t 0 : The n s weights are calculated such that X(s 0 , t 0 ) is unbiased and that the variance of the prediction error is minimal.A detailed description of the computation procedure can be found in geostatistical handbooks (Cressie, 1993;Goovaerts, 1997).
Three kriging algorithms can be distinguished (Goovaerts, 1997): -Simple kriging: the mean is known and constant thorough the space region; -Ordinary kriging: the mean is unknown but constant in the sub-domains of the space region; -Kriging with trend model or universal kriging: the mean is unknown, is local and varies smoothly in the sub-domains of the space region.In a similar way, the predictor X(s 0 , t 0 ), at the time instant t 0 , is a linear combination of the n t observations x 1 , …, x nt available at the space location s 0 .The prediction is done for each time instant (space location) separately and independently from the other time instants (space locations).Many applications can be found in the soil science literature: ordinary kriging using variograms converted from a basic pooled spatial variogram on the sampling time instants (Sterk & Stein, 1997;Stein & Sterk, 1999); ordinary kriging and universal kriging (Ettema et al., 1998;Figueira et al., 2001;Liu et al, 2009); and ordinary log-normal kriging (Ettema et al., 2000).

Spatial or temporal cokriging
For the same cases as before (spatial or temporal random fields) more accurate estimates, based on cokriging, can be obtained by using the information available at a precedent time instant as well as the contemporaneous data for spatial prediction (D'Agostino et al., 1998) or information from two different TRFs observed at two different space locations.Cokriging can also be used to estimate the temporal change for a spatially correlated soil property between two time instants using the pseudo cross-variogram (Lark, 2002;Papritz & Fluhler, 1994).Bogaert (1996) showed that the simple cokriging system between n t coregionalized space variables (one for each time instant) is equivalent to the simple space-time kriging.This equivalence was shown for a more general form of the covariance functions (Kyriakidis & Journel, 1999).Bogaert (1996) demonstrated also that space ordinary cokriging system with one unbiasedness condition is equivalent to the space-time ordinary kriging system when it is expressed in terms of covariance functions.Pseudo cross-variograms should replace the cross-variograms if the cokriging system has to be written in terms of variograms.However, the use of the unique unbiasedness constraint is restricted to the second order stationary situations.The author noted also that space-time ordinary kriging is preferable to space ordinary cokriging whenever it is possible because, although the differences in the prediction variances are negligible when data are abundant, it is not the case when some time instants involve limited data.We recall that for SRF, the prediction can be done at unsampled locations only for one or more of the observed time instants.As examples of application, we cite Bogaert &Christakos (1997) andD'Agostino et al. (1998).In the case of TRF, cokriging allows forecasting and hindcasting, the latter is useful in the estimation of missing values.Examples of application can be found in Rouhani et al. (1992a, b).Cokriging was used also to estimate a missing observation at a given time instant by considering the observations at this time instant as the variable of interest and the observations at other time instants as covariables (Stein, 1998).Let a given soil property be sampled at n 1 space locations during the first time instant, 1 () i Xs, i = 1, …, n 1 , and at n 2 (less than n 1 ) space locations later during the second time instant, 2 () j Xs , j = 1, …, n 2 .The kriging estimator of X 2 at an unsampled location s 0 is (Goovaerts, 1997): with 1 i and 2 j are the kriging weights relative to the variables X 1 and X 2 , respectively.The weights must satisfy the following conditions:

Space-time kriging
For the more general single STRF, two kriging algorithms can be distinguished, i.e. two step space-time kriging and anisotropic space-time kriging.They are presented in the next sections followed by other forms of space-time kriging.

Two-step space-time kriging
The prediction of the value of X 0 = X(s 0 , t 0 ) is a linear combination of n t predictors in space.In the first step, at each of the n t time instants t j , j = 1, …, n t , the predictor of the value of X(s 0 , t j ) is determined as (Stein et al., 1994): The kriging weights are determined in the standard manner for kriging in the space domain (Goovaerts, 1997).In the second step X(s 0 , t 0 ) is predicted, using a linear combination of the predictors determined at the first step.This predictor is (Stein et al., 1994): with j are the kriging weights from the second step.The predictor 00 ˆ(,) Xs t explicitly uses spatial variograms that change over time.This algorithm was applied in different situations (Douaik, 2005;Stein et al., 1994Stein et al., , 1998)).

Anisotropic space-time kriging
For this algorithm, the prediction is done in one and unique step capitalizing on a space-time variogram model which incorporates the spatial and temporal dependencies with a space-time anisotropy ratio (Bechini et al., 2000;Rouhani & Myers, 1990;Stein et al., 1994Stein et al., , 1998)).The advantage of the anisotropic over the two step space-time kriging is that the former allows predictions to be made at every point in space and time.Also it avoids the uncertainties related to the predictors determined at the first step for the latter.

Space-time cokriging
De Iaco et al. ( 2004) extended the Linear Model of Coregionalization (LMC) to the spacetime domain using generalized product sum models for the variograms.In addition, they proposed an extension of ordinary cokriging to the space-time domain.Although De Iaco et al. ( 2004) discussed the more general case of several secondary variables; we present here only the case of one variable of interest and one covariable.A linear space-time predictor of X(s, t) at an unsampled point (s 0 , t 0 ) in the space-time domain is: where Λ i (s 0 , t 0 ) are the (2x2) matrices of kriging weights whose elements (,) are the weights assigned to the value of the β th variable, β = 1, 2, at the i th sampled point to predict This corresponds to an extension of ordinary cokriging to the space-time domain.More detail on how to determine the kriging weights (,)  (2004).In their work, they compared space-time ordinary kriging to space-time ordinary cokriging.They found that the correlation coefficients between observed and predicted values were higher when the latter algorithm was used.Also, the cokriging error variances were lower than those from kriging.

Other forms of space-time kriging
Additional forms of space-time kriging exist in the literature.Among these variants, indicator kriging (Bilonick, 1988), kriging with external drift (Snepvangers et al., 2003), and factorial kriging (De Iaco et al., 2003;Rouhani et al., 1992a, b;Van Meirvenne & Goovaerts, 2002) were used.Now after that the classical geostatistical techniques of space-time interpolation have been laid out, we present a versatile and more general method of interpolation, i.e., the Bayesian Maximum Entropy (BME).

Bayesian maximum entropy
For most of the classical geostatistical methods of interpolation, i.e. kriging, the prediction is based solely on the hard data, accurate observed or measured data.However, other data sources can be useful, mostly, for space locations and/or time instants where hard data are missing and they can improve the accuracy of the predictions.Such data can be different types of soft data (interval, probabilistic, etc.), physical laws, moments of higher order, etc. Bayesian Maximum Entropy (BME) (Christakos, 2000;Christakos & Li, 1998) is a recent approach developed for the spatio-temporal mapping of natural processes using uncertain information.It offers the flexibility to incorporate various sources of physical knowledge.This incorporation of physical knowledge bases enables global prediction features and the adoption of probability distributions without a need for any assumptions like for example to be Gaussian.Christakos (2000) defined a knowledge base as 'a collection of knowledge sources relevant to the problem at hand to be involved by a reasoning process aiming at the solution of the problem'.In the BME framework, the total knowledge (K) available regarding a natural process is considered to be formed from two main bases: the general knowledge (G) and the specificatory knowledge (S).The general knowledge represents the knowledge that one has about the distribution of the natural variable to be mapped before any specific data.It encompasses physical laws, statistical moments of any order (including the mean and variogram or covariance function), multipoint statistics, etc.It is said to be general because it does not depend on the specific random field realization at hand.The specificatory knowledge includes the data specific to a given experiment or situation.It refers to a particular occurrence of the natural variable at a particular space location and a particular time instant.It is divided into two main categories depending on the accuracy of the data.The first category is called hard data which are exact measurements of the natural process and considered to be error-free.It encompasses accurate measurements obtained from realtime observation devices, computational algorithms, simulation processes, etc.The second category is called soft data.It involves uncertain observations, empirical charts, assessments by experts, etc.They are incomplete or qualitative data linked to opinions, intuition, etc. Different kinds of soft data are available, often in the form of interval domain data or of a probabilistic nature.To illustrate these notions of hard and soft data, we give here an example.Soil salinity can be assessed by measuring the electrical conductivity of the soil.It is measured either in the laboratory (ECe) or in the field (ECa).The measurement in the laboratory is done with high precision in standard conditions and is directly related to soil salinity.Thus, it is considered as an accurate measurement of the soil salinity and represents the hard data.The apparent electrical conductivity measured in the field (ECa) is only an indirect measurement of soil salinity and still needs to be calibrated with ECe values.Thus, based on a calibration model (frequently a regression model), ECa values are converted into 'predicted' ECe values.As these values are only predictions from a stochastic model, they are entailed with some uncertainty, and are consequently considered as soft data.The predicted ECe values and their corresponding standard deviations allow us to define the two kinds of soft data: the confidence intervals are built which give the interval soft data while the two statistic parameters are used to describe Gaussian PDFs which represent our probabilistic soft data.

The three steps of BME analysis
The BME has a double goal: informativeness (prior information maximization given the general knowledge) and cogency (posterior probability maximization given specificatory knowledge).The BME analysis is done in three main steps of knowledge acquisition, integration and processing: the structural or prior step, the meta-prior step, and the integration or posterior step.The goal of the prior step is the maximization of the information content considering only the general knowledge before any use of the data, which corresponds to the first goal of the BME.This step assumes an inverse relation between information and probability: the more informative an evaluation of a mapping situation is, the less probable it is to occur.This means that if a theory is general and vague, it includes more alternatives so it is more probable however it is less informative.A random field is completely defined by its multivariate probability distribution function (PDF), which forms the prior PDF.The latter should be derived by means of an estimation process that takes into consideration physical constraints under the form of prior information or knowledge.This information is measured, in the context of BME, using the Shannon's entropy function (Shannon, 1948), thus the E in BME.The expected information is given by the entropy function which needs to be maximized subject to the physical constraints provided by prior information or knowledge G, which justifies the M in the BME acronym: where map χ represents the available data and the one to be predicted, Gm a p f () χ is the unknown multivariate PDF associated with the general knowledge G, is a normalization constant, and α µ are Lagrange multipliers.
During the meta-prior step, the specificatory knowledge is collected and organized in appropriate quantitative forms that can be easily incorporated in the BME framework.The available data can be divided into two main types: the hard data and the soft data.The hard data are considered exact and accurate measurements of the natural process while the soft data are indirect and inaccurate measures of the variable of interest.The soft data may be in the form of intervals, probability distribution functions, etc.In the posterior step, the two knowledge bases (G and S) are integrated.The goal is the maximization of the posterior PDF given the total knowledge K, which is the second goal of BME.The prior PDF is updated, by considering the available site-specific knowledge (the data).This updating is performed using a Bayesian conditionalization, thus the B in BME.For interval soft data, the posterior PDF Kk f () χ is: and for probabilistic soft data, the posterior PDF Kk f () χ is: with soft χ the soft data and A is a normalization constant.

Some BME estimates
The posterior PDF (equations 44 and 45), which is not limited to the Gaussian type, describes fully the random field at the estimation point.It provides a complete picture of the mapping situation as well as different estimators and their associated estimation uncertainty.The prediction points lie, in general, on a regular grid and the predictions are used to create space-time maps.The different estimates presented in the following sections are specific to a general knowledge involving the first two statistical moments and a specificatory knowledge encompassing hard data as well as interval and probabilistic soft data.From the posterior PDF, different statistical parameters can be derived like the mode, the conditional mean, and the variance of the prediction error.The conditional mean is: www.intechopen.comwhile the variance of the prediction error is:

Kriging as a special case of BME
When the general knowledge is limited to the mean and covariance/variogram functions and the specificatory knowledge is restricted to the only hard data, the BME posterior PDF is Gaussian (the mean and the mode are equal) and the BME mean estimate is equivalent to the kriging (best minimum mean squared error, MMSE) estimate, which is the conditional mean.When the space-time random field is Gaussian, the kriging estimate becomes linear and optimal among all MMSE estimators.Consequently, BME is a more general interpolation approach and kriging is a special case in limiting situations.

Other space-time interpolation methods
The BME is one of the space-time interpolation methods that allow the incorporation of data with different information support and of different degrees of certainty and the knowledge in the form of physical laws.The data assimilation approach (Swinbank et al., 2003) is another method for such situations.It has two main variants: the Kalman filtering (Kalman, 1960) and the Bayesian data fusion (Fasbender et al., 2009).However, their application in the soil science literature is limited (Bierkens, 2001;Heuvelink et al., 2006), mostly for the latter since it was very recently developed.

Discussion and conclusions
This chapter has the main objective of reviewing the mostly used statistical methods for appraising the spatial and temporal variability of soil properties.There is a large panoply of techniques and, as stated in the introduction, the choice of the adequate model will depend on the behaviour of the phenomenon of interest (a priori knowledge), the objective of the analysis, and the nature of the available data.As a first step, the classical statistical method of coefficient of variation, which ignores the spatial and temporal coordinates, can be used as an exploratory data analysis to get an idea about the variability of data either for a given location, a given time instant or both.The coefficient of variation is considered as an indicator of global variability as opposed to local variability; the latter is identified using, for example, the variogram, a geostatistical tool.Another classical statistical method, the variance component analysis, in which the spatial and temporal contributions are considered as factors with random effects and the joint space-time contribution as interaction in the ANOVA model, can be used to determine their relative contribution to the total variance.Since this approach is directly related to the coefficient of variation, it is suffering from the same drawbacks as the later. www.intechopen.com Still using classical statistical methods, there are opportunities to consider the spatial and temporal dimensions of data using, for example, the concept of temporal stability.It is defined as the time invariant association between spatial location and statistical parameters of soil properties and can be assessed using either the relative differences or the Spearman rank correlation.A good test for the temporal stability is the correlation/regression between soil property values measured at consecutive time frames.A way to check the temporal stability of the spatial pattern is to test for a zero intercept and a unity slope of the regression between two consecutive time instants.If the slope of the regression between values at two consecutive time instants is not significantly different from one, then the regression intercept represents the constant change that has occurred between the two time instants.Based on the significance test of these two regression parameters, four possible scenarios can be distinguished: no net change between the two sampling dates; the mean soil property level has changed and there is a static (uniform) change in the soil property spatial pattern; the mean soil property level has not changed and the change in the soil property spatial pattern is dynamic (non-uniform) but, when averaged, this change is not significant; and the mean soil property level has changed and the change in the soil property spatial pattern is dynamic (non-uniform).The concept of relative differences is based on the same statistical parameters (mean and standard deviation) as the two precedent approaches (coefficient of variation and variance component analysis).This concept is based, like the two precedent approaches, on the assumptions of Gaussian distribution of data and the absence of autocorrelation.Both assumptions must be checked before the use of relative differences in the study of temporal stability.The normality can be checked using the standard statistical tests and the independence can be verified using either a test for autocorrelation like the Moran test or by computing the variogram and determining its range.Normality can be gained via data transformation and independence by the selection of sample locations separated by distances larger than the autocorrelation range.For non normally distributed data, but still non autocorrelated, the Spearman rank correlation can be used.It refers to the tendency of a soil property, measured at different locations in space, to maintain their relative ranking over time.In the existence of autocorrelation, its use is not recommended and it will be cautious to use statistical methods which consider the spatial and temporal nature of the given soil property.In addition, even for independent data, the Spearman rank correlation inform us only about the relative ranking of space locations between two time instants and do not provide a quantification of the temporal change.This quantification can be obtained using the paired-t test for comparing mean values.Thus, the paired-t test, still a classical statistical method, can be used following a temporal stability analysis based on the concept of relative differences.The concept of temporal stability handles each soil property separately and quantifies the temporal change even without testing if this change is significant or not.Dynamic spatial variation approach can test for change in the mean value (mean shift) as well as for uniform/non-uniform spatial variation between two time instants (spatial shift).Dynamic spatial variation refers to the spatially variable change of a soil property between two time instants, i.e., the change is not the same between two time instants for each location.Again, like with the relative differences concept, based on whether the mean difference is significantly different from zero, the spatial shift parameter is significantly higher than zero, or both, four possible scenarios can be distinguished: no net change with time in the soil property; the mean soil property level has changed and there is a static (uniform) change in the soil property spatial pattern; the mean soil property level has not changed and the change in the soil property spatial pattern is dynamic (non-uniform), again, when averaged, this change is not significant; and finally, the mean soil property level has changed and the change in the spatial pattern of the soil property is dynamic (non-uniform).
The dynamic spatial variation approach is based on the determination of linear regression between two soil properties; hence, it assumes that the residuals from the regression model are normally distributed with equality of variances and spatially non autocorrelated.All these assumptions must be checked.If spatial autocorrelation is present in the residuals, the estimates of the regression parameters and their variances are biased and the statistical tests are no longer valid.In this special case, spatial regression and/or geostatistical approaches are recommended.
The concept of temporal stability requires a moderate number of measurements in space and time, but once it is checked, the number of future measurements can be reduced to locations with measurements representing important statistics such as the average value.
Likewise, based on the dynamic spatial variation test, if it was found that there was a static (uniform) spatial variation between two dates, the sampling effort could be limited to the first soil property survey (an easy to measure and cheap variable) for the first date, while the second soil property (a hard and expensive variable) could be done for two dates but at a reduced number of locations.So the joint use of the concept of temporal stability and the temporal mean shift and spatial shift tests could result in a drastically reduced sampling effort.As illustration, Douaik et al. (2006) found that it was sufficient to measure soil salinity at just two locations and nine dates, from initially 20 locations and 19 time instants, to characterize reliably the average field soil salinity of 25 km² study area.The statistical method to use will depend on the data availability and the aim of the study.If data are available only for one of the two soil properties, there is only one possibility for checking the spatial pattern: the concept of temporal stability.For checking the temporal change in the average soil property, there is a choice between the paired-t test and the concept of temporal stability.If, on the other hand, data are available for both soil properties, there are more options: the paired-t test and the concept of temporal stability for checking the temporal change in the average level separately for both soil properties or the temporal mean shift test for checking the temporal change in the average level combining both soil properties.
Regarding the check of the temporal stability of the spatial patterns, the same choices are available except the paired-t test.
For an implicit and formal assessment of the spatial and temporal variability of soil properties, structural analysis is used with the goal of describing and modelling this variability using geostatistical tools like the covariance function, the correlogram, or the variogram.The models of spatial and temporal variability can be used, later, for interpolating between space locations and time instants for mapping soil properties.
Basically two interpolating methods can be used for space-time data: the geostatistical method of kriging and the Bayesian maximum entropy.Three main approaches for the analysis of space-time data can be distinguished: -Methods based on vectors of spatial random fields (SRF) for the case of many space locations and few time instants.It does not include the temporal dependence existing between observations and can predict only at the observed time instants.The spatial variability is modelled either through a separable spatial variogram for each time instant, or by a single spatial variogram considering time instants as replicates; -Methods based on vectors of time series (temporal random fields, TRF) for the case of long time series with few space locations.This approach doesn't take into account the spatial dependence and it predicts only at the observed locations.In a similar way as above, independent TRF or TRF as replicates in space can be considered.-Methods using an STRF, thus considering the joint space-time variability.For the first approach, the focus of the analysis is on smooth interpolated soil attribute maps over specific time instants.The intention is to capture single instantaneous snapshots (a static picture) of the natural process.The objective can be also the comparison of the various maps or the detection of the temporal persistence or change in the spatial patterns.The intention of the second approach is to capture a sequence of successive snapshots at single space locations, which give temporal profiles.Only the third approach includes both the spatial and temporal dependencies so the interpolation is more precise and can be done for unsampled time instants at unsampled space locations.The focus is here on video sequence of successive spatial pictures like a movie.When a theoretical model is fitted to the experimental (observed) space-time variogram, it becomes possible to tackle the problem of predicting attributes at unsampled space locations and time instants.Kriging, a family of generalized least squares regression algorithms, is one of the methods allowing these predictions.In case of spatial (temporal) random field, the dependence in data is modelled via a spatial (temporal) variogram.The prediction is done for each time instant (space location) separately and independently from the other time instants (space locations).More accurate estimates, based on cokriging, can be obtained by using the information available at a precedent time instant as well as the contemporaneous data for spatial prediction or information from two different TRFs observed at two different space locations.Cokriging can also be used to estimate the temporal change for a spatially correlated soil property between two time instants using the pseudo cross variogram.Space-time kriging is preferable to spatial cokriging whenever it is possible because, although the differences in the prediction variances are negligible when data are abundant, it is not the case when some time instants involve limited data.For the more general single STRF, two kriging algorithms can be distinguished: two-step space-time kriging and anisotropic space-time kriging.For the latter algorithm, the prediction is done in one and unique step capitalizing on a space-time variogram model which incorporates the spatial and temporal dependencies with a space-time anisotropy ratio.The advantage of the anisotropic over the two-step space-time kriging is that the former allows predictions to be made at every point in space and time.Also it avoids the uncertainties related to the predictors determined at the first step for the latter.In the case of the availability of spacetime data for more than one soil property, space-time cokriging can be used for improving the interpolation of the soil property of interest by considering its spatial and temporal dependence with one or more auxiliary soil properties.Additional forms of space-time kriging exist in the literature.Among these variants, indicator kriging, kriging with external drift, and factorial kriging were used in soil science.The prediction for most of the geostatistical methods of interpolation, i.e. kriging, is based solely on the accurate observed or measured data called hard data.However, other data sources can be useful, mostly, for space locations and/or time instants where hard data are missing and they can improve the accuracy of the predictions.Such data can be different types of uncertain data (soft data), physical laws, statistical moments of higher order, etc. Bayesian Maximum Entropy (BME) is a recent approach developed for the spatio-temporal mapping of natural processes using uncertain information in addition to accurate measured data.It offers the flexibility to incorporate various sources of physical knowledge enabling global prediction features and the adoption of probability distributions without a need for any assumptions like for example to be Gaussian.The BME has a double goal of informativeness (prior information maximization given the general knowledge) and cogency (posterior probability maximization given specificatory knowledge or the data at hand).The BME probability distribution function describes fully the random field at the prediction point.It provides a complete picture of the mapping situation and different statistical parameters can be derived like the mode, the conditional mean, and the variance of the prediction error.The prediction points lie, in general, on a regular grid and the predictions are used to create space-time maps.In limiting situations (when the general knowledge is limited to the mean and variogram and the specificatory knowledge is restricted to the only hard data), the BME probability distribution function is Gaussian and the BME mean estimate is equivalent to the kriging estimate.Consequently, BME is a versatile and more general method of interpolation and kriging is a special case in limiting situations.Besides BME, there are other space-time interpolation methods allowing the incorporation of data with different information support and of different degrees of certainty and the knowledge in the form of physical laws.The general data assimilation approach is one of these methods, with its two main variants: the Kalman filtering and the Bayesian data fusion.However, their application in the soil science literature is limited, mostly for the latter since it was very recently developed.In conclusion, different statistical methods are presently available for dealing with spatial and temporal variability of soil properties and others may be developed in the future.The behaviour of the phenomenon of interest (a priori knowledge), the objective of the research work, and the nature of the available data should help the practitioner in the choice of the adequate method(s).
data separated by a temporal distance (time period) , at the space location s i .Similarly to the spatial variogram, the temporal variogram (equation 25) is a special case of the space-time variogram (equation 20) when the spatial lag is zero meaning the same spatial location: in De Iaco et al.

Paired-t test Equation
(McClave & Sincich, 2006)e is a constant change or not between two time instants, but we still need to check if the mean values of a soil property between these two times are significantly different or not.To answer this question, the paired-t test(McClave & Sincich, 2006)is used.The test statistic t is computed as: www.intechopen.comthe  th variable,  = 1, 2 , at the unsampled point (s 0 , t 0 ).The predicted value includes two components: