Spatial modelling of soil hydraulic properties integrating different supports

Journal of Hydrology 511 (2014) 1–9

Contents lists available at ScienceDirect

Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol

Spatial modelling of soil hydraulic properties integrating different supports Ana Horta a,⇑, Maria João Pereira b, Maria Gonçalves c, Tiago Ramos d, Amílcar Soares b a

Faculty of Agriculture and Environment, The University of Sydney, Sydney, NSW 2006, Australia Centre for Natural Resources and Environment, Instituto Superior Técnico, Technical University of Lisbon, Lisbon, 1049-001 Lisbon, Portugal c Estação Agronómica Nacional (EAN), INIA, Instituto Nacional de Recursos Biológicos, Quinta do Marquês, Av. República, 2784-505 Oeiras, Portugal d CEER – Biosystems Engineering, Institute of Agronomy, Technical University of Lisbon, Tapada da Ajuda, 1349-017 Lisbon, Portugal b

a r t i c l e

i n f o

Article history: Received 13 July 2013 Received in revised form 28 December 2013 Accepted 12 January 2014 Available online 21 January 2014 This manuscript was handled by Andras Bardossy, Editor-in-Chief, with the assistance of Uwe Haberlandt, Associate Editor Keywords: Geostatistics Soil hydraulic properties Data integration Block simulation Uncertainty

s u m m a r y Modelling soil–water interactions provides important outputs for agriculture management and environmental monitoring. Most of the existing models rely on soil hydraulic properties (SHP) as input data. Geostatistical approaches based on stochastic simulations provide the spatial distribution of SHP and the uncertainty attached to their estimates. Frequently, SHP are measured in different data supports thus one needs to guarantee the integration of these different supports to provide a coherent and reliable model. One possible solution is to use block sequential simulation (Liu and Journel, 2009). Our work presents an application of this algorithm to map and quantify the spatial uncertainty of total porosity. Based on the simulation outputs and on the comparison with direct sequential simulation that does not account for the multiple supports of the data, we concluded that block sequential simulation should be used to produce reliable input data to dynamic soil–water models. Ó 2014 Elsevier B.V. All rights reserved.

1. Introduction Intensive agricultural activities in many regions of Europe are affecting the quality of main environmental resources such as soil and groundwater. It is of crucial importance to identify these environmental problems and to propose solutions to minimise the potential risks for public health and the ecosystems. With this purpose, models are developed to design irrigation schemes and to describe soil–water regimes and nutrient dynamics useful for environmental monitoring (Verhagen, 1997; Alphen et al., 2001). Research on this subject has produced diverse models towards quantifying and integrating the most important physical, chemical and biological processes active in the unsaturated zone of soils. The application of these models is fully dependent on the availability and quality of input data (namely, soil hydraulic properties (SHP)) that will determine the accuracy of modelling results (Wösten et al., 2001). However, lack of easily accessible, representative and accurate SHP is recognised as a limitation to soil–water models reliability (Wösten et al., 1999). The reasons behind the lack of quality of soil hydraulic data are mostly related with the

⇑ Corresponding author. Tel.: +61 0286271050. E-mail address: [email protected] (A. Horta). 0022-1694/$ - see front matter Ó 2014 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jhydrol.2014.01.027

techniques used to obtain direct measurements of these properties. The majority of these techniques remain relatively time consuming and therefore costly (Schaap et al., 1998; Minasny and McBratney, 2002). The cost-effectiveness of obtaining SHP can be improved by using indirect methods to predict hydraulic properties using more easily, widely available, routinely, or cheaply measured properties (Minasny and McBratney, 2002). Methods developed for this purpose use the so-called pedotransfer functions (PTFs) (Bouma, 1989). By using regression or data-mining models, PTFs relate hydraulic properties to easier to measure soil attributes such as soil texture, organic matter content and/or other data routinely measured in soil surveys. However, reliability of pedotransfer estimates should be carefully analysed since PTFs are based on general data sets and validation with ‘true’ field data is often lacking. Also, PTFs are specific for certain edafo-climatic conditions and soil family types, hence extrapolation to other areas beyond their geographical dataset should be regarded carefully (Minasny et al., 1999). Concerning the accuracy of SHP predictions, generally measured in terms of the uncertainty attached to the result, many PTFs have been developed without including any uncertainty calculation (McBratney et al., 2011). One possible way to overcome PTFs limitations and to provide an accurate estimate for SHP is to use a spatial inference model

2

A. Horta et al. / Journal of Hydrology 511 (2014) 1–9

able to integrate the spatial variability of SHP (Wösten et al., 2001). This can be implemented using geostatistical models. Geostatistical estimation of SHP, or any other soil property, is about predicting the attribute at unsampled locations. Few studies have previously adopted this approach for SHP mapping. The work by Sinowski et al. (1997) shows an example of how geostatistical estimation (kriging) can be used to generate maps of water retention curves (integrating PTF information). The results obtained are a good example of the advantages of using geostatistics for the quantification of SHP. Also, the benefits of using geostatistical models to predict any other soil property have been recognised by several authors (for example, Goovaerts (1999) and Webster (2000)). In this paper we present a geostatistical approach that goes beyond previous studies in that it provides an estimate of the spatial distribution of SHP together with an evaluation of the uncertainty associated with the result, while accounting for the different data support. When modelling uncertainty we face two main decisions: whether to model uncertainty locally or spatially, and how to account for the different supports of the data and of the model discretization. Goovaerts (2001) recommends choosing the type of uncertainty modelling guided by practical criteria. Local uncertainty modelling is the most traditional approach and it amounts to assigning a confidence interval to the estimate, or, at most, to build a local uncertainty distribution from which to derive an estimate and its associated uncertainty. However, when the aim of the study is to use the estimates as input to another model for further prediction (such as would be the case here in which we propose to use the estimates as input to a subsurface water flow model) the combination of a map with local estimates and a map of local variances is of no use (Gómez-Hernández and Wen, 1998). For this reason, we opt for a stochastic approach to spatial uncertainty modelling by generating multiple, equally likely realizations of the SHP from their joint distribution. Later, the postprocessing of these realizations through the flow model will provide a model of uncertainty about the target result. The SHP considered in this study is total porosity (Ptotal). The data available for our study consists in measurements of Ptotal made for samples collected in the same soil layer but not at the same depth hence the measured values are not reported for the same support. This usually happens when soil surveys are performed in different conditions such as different surveys or sampling protocols. Commonly, the different data supports are ignored and all data are treated as point data. Alternatively, a change of support can be applied to obtain an interpolated value for a larger support, a method largely applied in geostatistics and known as point-to-area interpolation or block kriging (Isaaks and Srivastava, 1989); or measurement values can be considered as average SHP values that should be downscaled to a smaller support. Several authors (Gotway and Young, 2002, 2005; Kyriakidis, 2004; Goovaerts, 2008) have proposed to use kriging to predict point values from block data, an approach referred to as ‘area-to-point’ kriging. We wish to account for the data at their sampling support without the need for upscaling or downscaling. This can be accomplished using block sequential simulation (BSSIM) as proposed by Liu and Journel (2009). In sequential simulation, each realization is generated by sequentially visiting each node and generating a value at each location. At any given location, a local probability distribution function is built after solving a kriging system. This probability distribution must account for the conditioning data and for previously simulated values. BSSIM is able to incorporate block data, as well as point data, as conditioning data, by using pointto-point, area-to-point, point-to-area and area-to-area covariances. In addition, BSSIM uses the DSS algorithm (Journel, 1994; Soares, 2001) to build the local probability distribution in original space

of the variable without the need of any transformation, and it is free of any distribution assumption on the attribute. In this paper we demonstrate how to spatially characterize the uncertainty of SHP accounting for different supports using BSSIM. In order to show the importance of properly accounting for the different data supports, the BSSIM results are compared with those resulting from the application of direct sequential simulation considering all data as point values. 2. Materials and methods 2.1. Soil database and study area The soil data used in this work was extracted from a Portuguese database (PROPSOLO, Ramos et al., 2007) comprising 347 georeferenced (Lisboa Hayford Gauss IGeoE projection with datum Lisboa Hayford) soil profiles, collected in several locations across Portugal, from 1977 to 2011. During this 30-year period, profiles were collected for different public and private projects conducted by EAN (Estação Agronómica Nacional), a Portuguese research centre for agriculture and soil science development. Some of the data included in PROPSOLO were provided to HYPRES (Hydraulic Properties of European Soils; Wösten et al., 1999), a European database for soil hydraulic data. Each sampling campaign complied with Portuguese and international procedures (Cardoso and Fernandes, 1972; FAO, 2006). For each soil profile, a qualitative description was made, each horizon was sampled, and soil properties were determined in the laboratory. All soil hydraulic properties were measured on undisturbed soil cores collected in the soil horizons/layers of the different soil profiles included in the database (Ramos et al., 2006). For this paper, 46 soil profiles were chosen in an area located in the South of Portugal, with circa 3800 km2 (Fig. 1). The major soil groups in this area are predominantly Luvisols, Cambisols, and Vertisols. The criteria for selecting these soil profiles included the number and the proximity of sampling points as well as the homogeneity of soil type and land use. The selected area was converted from rainfed agriculture to irrigation over the last two decades. To sustain the irrigation projects a greater number of soil profiles were collected than in other Portuguese regions. For our study, we selected only the topsoil layer for each location. Topsoil depths varied between 10 cm and 48 cm. Total porosity (Ptotal) was obtained from the maximum holding capacity of 100 cm3 undisturbed soil cores in volumetric basis. Several 100 cm3 samples were taken from each core, and their average value was provided as Ptotal for the given core. For the 46 topsoil samples collected in the study area, the measured Ptotal varies between 0.33 cm3/cm3 and 0.64 cm3/cm3, its histogram is shown in Fig. 2. The spatial distribution of Ptotal is displayed in Fig. 1. The variability observed for Ptotal can be explained by the effects of local soil texture, soil structure, organic matter content, particle dispersion, soil crusting, changes in the concentration and ionic composition of the soil solution, microbiological activity, and the loading (stress) history on the soil profile (van Genuchten and Šimu˚nek, 1996; Gupta et al., 2006; Strudley et al., 2008). In the next sub-section, we describe the geostatistical approach proposed not only to estimate Ptotal but to address two issues: how to deal with different sampling supports and how to provide an uncertainty measure regarding the variability of Ptotal in the study area and not only at specific locations. 2.2. Block sequential simulation The geostatistical model proposed in our work uses the block sequential simulation (BSSIM), an algorithm included in the BGeost

A. Horta et al. / Journal of Hydrology 511 (2014) 1–9

3

Fig. 1. Study area and measured values for total porosity.

mean=0.427; variance=0.005

8 6 0

2

4

Frequency

10

12

Ptotal

0.30

0.35

0.40

0.45

0.50

0.55

0.60

0.65

Ptotal

Fig. 2. Histogram and main statistics of total porosity data measured for the study area.

package developed by Liu and Journel (2009) for the integration of data with different support volumes. By support we refer to the measure (length, area or volume) usually associated with the physical representation of the data values. Most variables in geostatistics are usually associated with ‘‘points’’ that represent a location in space. However, sometimes we need to differentiate the support in which the variable is observed. That is often due to the scale at which the variable is measured. For example, precipitation can be estimated at a coarse spatial resolution by global climate models or can be measured at a fine scale resolution in hydrometric stations. Some variables are only available at the coarse scale (classically termed ‘‘block

support’’ in geostatistics), for example, health data or socioeconomic data provided at the county level. In the BGeost package both point and block support data are considered simultaneously. This is different from the common change of support approach. In a typical change of support application, one aims to upscale or downscale the data (Li et al., 2011a,b), depending if the initial datum is provided in point or block support. The kriging algorithm can be used to perform this upscaling or downscaling since it was developed to accommodate different spatial supports for both the data and the predicted unit (Journel and Huijbregts, 1978; Goovaerts, 2008). Point-to-area kriging (Isaaks and Srivastava, 1989) is used to provide a value representative of a block when only point data are available. On the other hand, area-to-point kriging (Gotway and Young, 2002, 2005; Kyriakidis, 2004; Goovaerts, 2008) is used to disaggregate block into point data. Recently, Goovaerts (2011) proposed a geostatistical approach to incorporate point and areal data for spatial interpolation using two methods. One relies on area-to-point kriging whereas the other combines point and areal data in the same kriging system using a similar approach as proposed by Liu and Journel (2009). By using the BGeost approach we can account for the data in their original sampling support without the need for upscaling or downscaling. In our case study, we have a soil hydraulic property (total porosity) measured in samples that were collected in the same soil layer but with different depths depending on the sampling location. Instead of assuming that all data have been taken over the same constant depth for the soil layer, we chose to work with the true sampling depths and analyse the importance of accounting for the different measurements supports to map the spatial variability and uncertainty. BGeost comprises three algorithms: BKRIG (block kriging estimation), BSSIM (block sequential simulation) and BESIM (block error simulation). For this work, we chose BSSIM to obtain an

4

A. Horta et al. / Journal of Hydrology 511 (2014) 1–9

estimate of the spatial distribution of total porosity (Ptotal) and the uncertainty associated. We opted for a stochastic simulation approach to obtain multiple, equally likely realizations of Ptotal. Each Ptotal realization can be processed afterwards using a transfer function (for example, a flow simulator) where each realization will provide a unique value for each model response. The distribution of the response values corresponding to the input realizations provides a model of uncertainty about the target result. Also, the set of alternative realizations provides a visual and quantitative measure of spatial uncertainty about the spatial distribution of Ptotal (Goovaerts, 1997). BSSIM is a simulation algorithm based on the sequential simulation principle (Gómez-Hernández and Cassiraga, 1994; Goovaerts, 1997), i.e., it visits sequentially all nodes of the simulation grid and assigns a value to each node by drawing from a probability distribution function. This probability distribution function is built from the solution of a kriging system in which both the conditioning data, as well as the previously simulated node values in a neighbourhood are used as data to compute both the kriging estimate and the kriging variance. Two key points of BSSIM as compared with other sequential simulation algorithms are: (i) in order to use the block data as conditioning data, the covariances between block supports and between block and point supports must be accounted for, and (ii) by using the fundamental principle of direct sequential simulation, the local conditional distribution is built from the prior experimental distribution: there is no gaussianity implied, and for this reason data can be treated in their original space without any need to transform them. Block data are input to the algorithm by specifying their value and the number of discretizing points that they cover. In our case, since data are taken up to 48 cm deep, we decided to consider that the point support is of 15 cm; block data will cover one, two or three points depending if their support is up to 15, 30 or 45 cm deep, respectively. This information will be internally used by the algorithm to compute the block-to-block and the block-topoint covariances together with the formal relationship between point and block support data, which, in the BGeost algorithms is a linear average relationship:

BðV a Þ ¼

1 jV a j

Z

0

La ðZðu0 ÞÞdu ; 8a

ð1Þ

kðuÞ ¼ K 1 k

ð2Þ

1

where K is the data covariance matrix and k is the vector of datato-unknown covariances.K1and k are built as (Liu and Journel, 2009):

K¼

– For each realization: For each location u along a random path visiting all nodes to be simulated, i. search for the conditioning data consisting of point and block data, and previously simulated values; ii. compute local block-to-block, block-to-point, point-to-block and point-to-point covariances; iii. build and solve the mixed-scale kriging system, which provides the local kriging mean and variance; iv. define the local conditional cumulative distribution function (ccdf) with its mean and variance given by the kriging estimate and variance as described in the DSSIM implementation; v. draw a value from the local ccdf defined in (iv) and add the simulated value to the data set; – repeat the previous steps to generate another realization. BSSIM honours the experimental data, both point and block data, and reproduces the point histogram and main spatial continuity patterns as modelled by the experimental variogram. In the end, we obtain point support realizations for the distribution of Ptotal in the topsoil layer. In the next section, we demonstrate an application of BSSIM for the spatial characterization of uncertainty of Ptotal accounting for different supports. The version of BSSIM used is the one implemented in the BGeost package built into the Stanford Geostatistical Modelling Software, version 2.1 (Remy et al., 2008).

3. Results and discussion

Va

where Z(u) is the point support attribute, Va is the volume support, B(Va) is the block value and La is a linear function. At each node u in the realization, the kriging system that is built results in the following kriging weights k(u)

"

is important since the linear relation between point and block data occurs in original space, and will not be honoured if a non-linear transformation of the data is performed (for instance trying to make the transformed data follow a Gaussian distribution). For more details about the direct sequential simulation see, Soares (2001) and Horta and Soares (2010). The BSSIM algorithm can be summarized as follows:

C PP0

C PB

C tPB

C BB0

#

" k¼

C PP0 C BP0

# ð3Þ

where C denotes the point (P) covariance submatrix, C denotes a covariance submatrix involving block (B) data and P0 is the point to be estimated. The construction of matrices K and k as functions of the pointto-point covariance, the block-to-point covariance and the blockto-block covariance is necessary when block and point data are mixed for the simulation of a given support. In addition, using the DSS paradigm (the simulated value is drawn from a local distribution built from an interval of the global distribution defined by the simple kriging mean and variance computed at each node), the algorithm can work in the original space without the need of any (Gaussian) transformation. This last point

3.1. Spatial distribution of porosity To apply the proposed methodology the first step is to define the support of each datum according to the sampling depth of the topsoil layer. The smallest support, i.e., the point support in the BSSIM algorithm, is 15 cm; all other data will have a larger support which will be discretized in elements of 15 cm for the purpose of computing the area-to-point covariances and the area-to-area covariances. The linear operator used to relate point and block data is the arithmetic average. For the 46 locations (see Fig. 3), only 4 sampling points were classified as point data. The remaining 42 points were divided in 24 blocks straddling two points, and 18 blocks straddling three points. The high Ptotal values (above the mean) are mostly associated with block data (only one location defined as point presents high Ptotal). Next, we need the covariance function of the point data, which can be derived from the variogram function. The model for the variogram is based on the experimental variogram computed with the 46 data values. For the purpose of computing the variogram, all data were considered as point support data. The model finally used is an exponential isotropic variogram with a range of 20 km and a sill of 0.047 after manually fitting the experimental variogram (Fig. 4). These data and parameters were used to perform BSSIM over a grid of 235  133  3 elements, each element of size

A. Horta et al. / Journal of Hydrology 511 (2014) 1–9

5

Fig. 3. Block and point data defined for the study area.

300 m  300 m  0.15 m. Only the topmost layer was simulated. A set of 30 realizations was obtained representing the spatial distribution of Ptotal on the top layer. Fig. 5a shows an example of two equiprobable realizations whereas Fig. 5b presents the average of the 30 realizations obtained using BSSIM. This average image suggests that high values of porosity (exceeding the mean and the upper quantile of the experimental histogram) are likely to be found in the East/Southeast quadrant of the study area. Low values (below the mean) are preferentially located in two poles in the Northwest/Southwest of the study area. As in every simulation output, each BSSIM realization should reproduce the experimental variogram (Fig. 4) and histogram (Fig. 2).

Ptotal histograms of simulated values presented some differences with respect to the experimental histogram that are worth exploring to better understand BSSIM performance. As an example, the histograms of two equiprobable BSSIM realizations are shown in Fig. 6a and b. For the Ptotal simulations, the histograms tend to a normal distribution, i.e., the most populated class contains values closer to the mean than what was observed in the experimental histogram. This can be explained because the classes of low values (below the mean) in the experimental histogram represent mostly clusters of low values that bias the experimental count for those classes, giving them a larger frequency than what they should have. During the simulation procedure, these biased classes of the experimental

Fig. 4. Experimental variogram for total porosity (fitted with an exponential model with a 20 km range).

6

A. Horta et al. / Journal of Hydrology 511 (2014) 1–9

Fig. 5a. Spatial distribution of total porosity (2 realizations obtained using BSSIM).

– the porosity value that is exceeded by 90% of the simulated values (0.1 quantile) (Fig. 7c) and the porosity value that is exceeded by 10% of the simulated values (0.9 quantile) (Fig. 7d).

Fig. 5b. Spatial distribution of total porosity (average of 30 realizations obtained using BSSIM).

histogram tend to be spatially declustered, explaining the shift observed in the histograms of simulated values. 3.2. Spatial uncertainty characterization The average map presented before (Fig. 5b) provides the main patterns of the spatial heterogeneity of Ptotal honouring the spatial variability model given by the experimental histogram and variogram. Spatial uncertainty was quantified using the following statistics computed from the distribution of values obtained through the ensemble of realizations at each grid node: – the variance of the simulated values (Fig. 7a); – the probability of exceeding the experimental Ptotal mean (0.42) (Fig. 7b),

The variance map (Fig. 7a) shows low uncertainty for high and low Ptotal estimates displayed in our average map. An exception is made for the South/Southeast part of the study area where available data is sparse resulting in higher uncertainty. High Ptotal values are expected only near the sampling locations with high values, as indicated by the probability map in Fig. 7b. This assessment is confirmed by the 0.1 quantile map (Fig. 7c) that identifies high Ptotal values with a 90% chance of being exceeded in the same zones. This information combined with the previous variance map allows for a better definition of high Ptotal areas, especially in what concerns the South/Southeast area. The lowest values presented in Fig. 7d (0.9 quantile map) indicate areas where there is a 90% chance that the true porosity will be below it, that is, areas where Ptotal values must be certainly small. These maps are just some of the statistics that can be derived from the ensemble of realizations to analyse the uncertainty about the attribute. In addition, these maps can serve as input to a model, so that the analysis of the ensemble of responses obtained can serve to measure the uncertainty in model responses.

3.3. Importance of integrating different supports To show the importance of properly accounting for the different data supports, we performed a comparison between the results obtained by BSSIM and the ones obtained using another sequential simulation algorithm assuming that all experimental values have the same (point) support. For this purpose the direct sequential simulation (DSS) algorithm was used. DSS used the Ptotal

A. Horta et al. / Journal of Hydrology 511 (2014) 1–9

30000

BSSIM 1

mean=0.439; variance=0.003

15000 0

5000

10000

Frequency

20000

25000

(a)

7

0.30

0.35

0.40

0.45

0.50

0.55

0.60

0.65

0.60

0.65

Ptotal BSSIM 2 mean=0.437; variance=0.005

15000 10000 0

5000

Frequency

20000

25000

(b)

0.30

0.35

0.40

0.45

0.50

0.55

Total_porosity_2

Fig. 6. (a and b) Histogram of the simulated values obtained for one BSSIM realization.

measured in the topsoil layer disregarding its support. The variogram parameters used were the same as the ones applied in BSSIM. This comparison was based in the average maps shown in Fig. 8(a and b) (ensemble averages of the 30 realizations generated by BSSIM and DSS) and in the probability maps of Fig. 8(c and d). Fig. 8(c and d) represents, at each location, the probability that porosity is larger than the upper quartile of the experimental histogram (0.48) as determined by the count of realizations which, at that location, exceed that threshold, the value 1 represents probability values higher than 0.5. The comparison between the average maps indicates that areas with low Ptotal are larger in the DSS map, namely in the Northwest/Southwest area. On the other hand, areas with higher Ptotal are larger in the BSSIM map as shown in the maps in Fig. 8(a and b). There is also a substantial difference on the estimates surrounding the point A on the Southwest area. This distribution of Ptotal estimates is also confirmed by both histograms of the simulated values presented in Fig. 9(a and b). The analysis of the histogram shows how both average images differ in terms of the distribution of Ptotal: DSS concentrates most of Ptotal estimates in the lower classes of the histogram whereas BSSIM accumulates more values in histogram classes that exceed the upper quantile value.

Fig. 7. Uncertainty evaluation using four different statistics: variance (a), probability above the mean (b), 0.1 quantile (c), 09. quantile (d).

8

A. Horta et al. / Journal of Hydrology 511 (2014) 1–9 BSSIM mean=0.441; variance=0.001

10000 0

5000

Frequency

15000

(a)

0.30

0.35

0.40

0.45

0.50

0.55

0.60

0.65

0.60

0.65

Ptotal

10000 0

5000

Frequency

(b)

15000

DSS mean=0.427; variance=0.001

0.30

0.35

0.40

0.45

0.50

0.55

Ptotal

Fig. 9. Comparison between the histogram of simulated values obtained by BSSIM (a) and DSS (b).

These results raise the question of whether accounting for different data supports overestimates the porosity content or reproduces better the variability given by the data. The answer, in short, is that considering block support data as point data will always underestimate the variance of the final point realizations. This underestimation of the variance, for a skewed distribution such as that of porosity, will result in an underestimation of the high values. The reason for this underestimation of the variance is that block values (being averages of point values) belong to a population with a smaller variance than the point values, for this matter, there is a need to ‘‘inflate’’ the variance of the conditioning data if they are treated as point data. By properly accounting for the block support, through the point-to-block covariance, the block kriging equations will correct the variance when estimating point values from block data. In conclusion, we consider that integrating different supports provides a more realistic reproduction of the variability given by experimental data. Therefore, BSSIM outputs are more reliable when combined with the uncertainty assessment presented before. 4. Conclusions Fig. 8. Comparison between average maps generated by BSSIM (a) and DSS (b) outputs and probability maps for total porosity estimates higher than the upper quantile for BSSIM (c) and DSS (d) (1 represents a probability higher than 0.5).

In this paper we propose a geostatistical approach using stochastic simulation to spatially characterize the uncertainty of soil hydraulic properties. The use of a stochastic simulation algorithm

A. Horta et al. / Journal of Hydrology 511 (2014) 1–9

is particularly important when the estimated values are to be used by other prediction models. In the case of soil hydraulic properties, the estimates are to be used as input to hydrological models. Stochastic simulation will not only provide realistic realizations of the spatial distribution of the hydraulic properties of interest, but also a global measure of uncertainty, instead of local uncertainty estimates (as given by kriging-type algorithms). The soil hydraulic property used in our case study was total porosity (Ptotal). The data available comprised measurements of Ptotal made for samples collected in the same soil layer but at different depths and different supports. The different sampling supports are accounted for by using the block sequential simulation algorithm (BSSIM) developed by Liu and Journel (2009). This algorithm treats simultaneously point and block data using a framework that uses kriging to combine point and block covariances and stochastic simulation to compute an estimate at unsampled locations. The output of BSSIM consists in point support realizations of Ptotal in the topsoil layer, from which estimates of the main patterns of variability of total porosity can be derived, together with different measures of uncertainty around these estimates. A comparison with another stochastic simulation algorithm that does not take into account the data support was performed to show the importance of accounting for the support of the data. Based on the results obtained for this case study, we conclude that it is important to account for different supports to characterize the spatial variability of porosity since the output reproduces the true variability in the data. This is achieved by assuming that block data are linear averages of point data; that block data can be represented by a set of aligned points according to specific intervals; and that each estimate is conditioned by both block and point using mixed covariance values. An algorithm as BSSIM implements all these premises plus uses a direct simulation approach that allows using the data in their original space. One should keep in mind that this approach cannot be applied on attributes that do not average linearly, such as permeability; however, it could be applied to a point transformation, as is often the case of permeability, for which log-permeability can be considered that averages linearly. For our case study, the proper integration of the values measured at different depths in a coherent and unique model allowed for a better definition of areas with high porosity that otherwise could have been underestimated. Finally, the spatial modelling of soil properties using a geostatistical approach as the one proposed here is an alternative to the use of pedotransfer functions. Acknowledgements First author acknowledges Professor Jaime Goméz-Hernández for his outstanding support in reviewing this paper. This work was possible due to the financial support of the Portuguese Foundation for Science and Technology (FCT) through the postdoctoral scholarship SFRH/BPD/70315/2010. References Alphen, B., Booltink, H., Bouma, J., 2001. Combining pedotransfer functions with physical measurements to improve the estimation of soil hydraulic properties. Geoderma 103 (1–2), 133–147. Bouma, J., 1989. Using soil survey data for quantitative land evaluation. Adv. Soil Sci. 9, 177–213. Cardoso, J.C., Fernandes, J.F., 1972. Normas para a observação e descrição de perfis e para a colheita de amostras, second ed. Serviço de Reconhecimento e de Ordenamento Agrário, Secretaria de Estado da Agricultura, Ministério da Economia, Lisboa.

9

FAO, 2006. Guidelines for soil description. In: World Soil Resources Reports, fourth ed. Food and Agriculture Organization of the United Nations, Rome. Gómez-Hernández, J., Cassiraga, E., 1994. Theory and Practice of Sequential Simulation. In: Armstrong, M., Dowd, P.A. (Eds.), Geostatistical Simulations. Kluwer Academic Publishers, Dordrecht, pp. 111–124. Gómez-Hernández, J., Wen, X.-H., 1998. To be or not to be multi-Gaussian? A reflection on stochastic hydrology. Adv. Water Resour. 21 (1), 47–61. Goovaerts, P., 1997. Geostatistics for Natural Resources Evaluation. Oxford University Press, New York. Goovaerts, P., 1999. Geostatistics in soil science: state-of-the-art and perspectives. Geoderma 89, 1–45. Goovaerts, P., 2001. Geostatistical modelling of uncertainty in soil science. Geoderma 103, 3–26. Goovaerts, P., 2008. Kriging and semivariogram deconvolution in presence of irregular geographical units. Math. Geosci. 40 (1), 101–128. Goovaerts, P., 2011. A coherent geostatistical approach for combining choropleth map and field data in the spatial interpolation of soil properties. Eur. J. Soil Sci. 62 (3), 371–380. Gotway, C.A., Young, L.J., 2002. Combining incompatible spatial data. J. Am. Stat. Assoc. 97 (459), 632–648. Gotway, C.A., Young, L.J., 2005. Change of support: an inter-disciplinary challenge. In: Renard, P., Demougeot-Renard, H., Froidevaux, R. (Eds.), GeoENV V— Geostatistics for Environmental Applications. Springer, Berlin, pp. 1–13. Gupta, S.D., Mohanty, B.P., Köhne, J., 2006. Soil hydraulic conductivities and their spatial and temporal variations in a vertisol. Soil Sci. Soc. Am. J. 70, 1872–1881. Horta, A., Soares, A., 2010. Data integration model to assess soil organic carbon availability. Geoderma 160 (2), 225–235. Isaaks, E.H., Srivastava, R.M., 1989. Applied Geostatistics. Oxford University Press, New York. Journel, A.G., 1994. Modeling Uncertainty: Some Conceptual Thoughts. In: Dimitrakopoulos, R. (Ed.), Geostatistics for the Next Century. Kluwer, Dordrecht, pp. 30–43. Journel, A.G., Huijbregts, C.J., 1978. Mining Geostatistics. Academic Press, New York. Kyriakidis, P.C., 2004. A geostatistical framework for area-to-point spatial interpolation. Geogr. Anal. 36 (3), 259–289. Li, L., Zhou, H., Gómez-Hernández, J.J., 2011a. A comparative study of threedimensional hydraulic conductivity upscaling at the macrodispersion experiment (MADE) site, Columbus air force base, Mississippi (USA). J. Hydrol. 404 (3–4), 278–293. Li, L., Zhou, H., Gómez-Hernández, J.J., 2011b. Transport upscaling using multi-rate mass transfer in three-dimensional highly heterogeneous porous media. Adv. Water Resour. 34 (4), 478–489. http://dx.doi.org/10.1016/ j.advwatres.2011.01.00. Liu, Y., Journel, A., 2009. A package for geostatistical integration of coarse and fine scale data. Comput. Geosci. 35 (3), 527–547. McBratney, A., Minasny, B., Tranter, G., 2011. Necessary meta-data for pedotransfer functions. Geoderma 160 (3–4), 627–629. Minasny, B., McBratney, A., 2002. The efficiency of various approaches to obtaining estimates of soil hydraulic properties. Geoderma 107 (1–2), 55–70. Minasny, B., McBratney, A., Bristow, K., 1999. Comparison of different approaches to the development of pedotransfer functions for water–retention curves. Geoderma 93 (3–4), 225–253. Ramos, T.B., Goncalves, M.C., Martins, J.C., van Genuchten, M.T., Pires, F.P., 2006. Estimation of soil hydraulic properties from numerical inversion of tension disk infiltrometer data. Vadose Zone J. 5 (2), 684–696. Ramos, T.B., Gonçalves, M.C., Martins, J.C., Pires, F.P., 2007. PROPSOLO – Base de dados georreferenciada de propriedades do solo. In: Proceedings ‘‘II Congresso Nacional de Rega e Drenagem’’. Remy, N., Boucher, A., Wu, J., 2008. Applied Geostatistics with SGeMS: A Users Guide. Cambridge University Press, New York. Schaap, M.G., Leij, F.L., van Genuchten, M.Th., 1998. Neural network analysis for hierarchical prediction of soil hydraulic properties. Soil Sci. Soc. Am. J. 62, 847– 855. Sinowski, W., Scheinost, A., Auerswald, K., 1997. Regionalization of soil water retention curves in a highly variable soilscape, II. Comparison of regionalization procedures using a pedotransfer function. Geoderma 78 (3–4), 145–159. Soares, A., 2001. Direct sequential simulation and cosimulation. Math. Geol. 33 (8), 911–926. Strudley, M.W., Green, T.R., Ascough II, J.C., 2008. Tillage effects on soil hydraulic properties in space and time: state of the science. Soil Till. Res. 99, 4–48. van Genuchten, M.Th., Šimu˚nek, J., 1996. Evaluation of pollutant transport in the unsaturated zone. In: Rijtema, P.E., Eliáš, V. (Eds.), Regional Approaches to Water Pollution in the Environment. Kluwer Acad. Publ, Dordrecht, pp. 139– 172. Verhagen, J., 1997. Site specific fertiliser application for potato production and effects on N-leaching using dynamic simulation modeling. Agric. Ecosyst. Environ. 66 (2), 165–175. Webster, R., 2000. Is soil variation random? Geoderma 97, 149–163. Wösten, J., Lilly, A., Nemes, A., Le Bas, C., 1999. Development and use of a database of hydraulic properties of European soils. Geoderma 90 (3–4), 169–185. Wösten, J., Pachepsky, Y., Rawls, W., 2001. Pedotransfer functions: bridging the gap between available basic soil data and missing soil hydraulic characteristics. J. Hydrol. 251 (3–4), 123–150.