Soil physical and hydraulic properties in the Donato stream basin, RS, Brazil. Part 2: Geostatistical simulation

Propriedades físico-hidráulicas do solo na bacia do arroio Donato, RS, Brasil. Parte 2: Simulação geoestatística

ABSTRACT

Geostatistical simulation has been the most promising and used technique for the analysis of uncertainties of soil physical and hydraulic properties, with high spatial heterogeneity. This study carried out a stochastic analysis of saturated hydraulic conductivity (K_sat) and soil water retention curve parameters in the Donato stream basin, located in the municipality of Pejuçara, in the state of Rio Grande do Sul, Brazil, with geographic coordinates between 28º 25’ 34” S and 53º 40’ 30” W, and 28º 24’ 50” S and 53º 41’ 30” W, 590 m of altitude. Soil samples were collected during the period from August to November of 2012. Sequential Gaussian simulation technique was used to generate 100 random fields of each variable. The results showed great uncertainties for K_sat and the parameter α of the soil water retention curve. The uncertainties between the percentiles 5 and 95% for K_sat indicated values from 24 to 44 cm d^-1, and for the parameter α, the uncertainties could be estimated from 0.622 to 1.122 cm^-1.

Key words:
sequential Gaussian simulation; random fields; spatial variability analysis

RESUMO

A simulação geoestatística tem sido a técnica mais promissora e utilizada para a análise de incertezas de propriedades físico-hidráulicas do solo com grande heterogeneidade espacial. Este estudo realizou uma análise estocástica da condutividade hidráulica saturada (K_sat) e dos parâmetros da curva de retenção de água no solo na bacia do arroio Donato, localizada no município de Pejuçara, na região Noroeste do estado do Rio Grande do Sul, com coordenadas geográficas entre 28º 25’ 34” S e 53º 40’ 30” O, e 28º 24’ 50” S e 53º 41’ 30” O, 590 m de altitude. Amostras de solo foram coletadas durante o período de agosto a novembro de 2012. Foi empregada a técnica de simulação sequencial Gaussiana para geração de 100 campos aleatórios de cada variável. Os resultados revelaram consideráveis incertezas para K_sat e o parâmetro α da curva de retenção de água no solo. As incertezas entre os percentis de 5 e 95% para K_sat indicaram valores de 24 a 44 cm d^-1, e para o parâmetro α, as incertezas foram estimadas de 0,622 a 1,122 cm^-1.

Palavras-chave:
simulação sequencial Gaussiana; campos aleatórios; análise de variabilidade espacial

Introduction

Geostatistical techniques have often been used to model uncertainties about the unknown value of a variable. One of the ways to quantify such uncertainties is through geostatistical simulation, which consists of stochastic methods based on obtaining several random fields for spatial distribution of a variable (Goovaerts, 1997).

Geostatistical simulation has been the most promising technique for the analysis of uncertainties. Simpler estimation methods, such as ordinary kriging, are based on formulas with weighted averages, and consequently, result in estimates that have significant levels of smoothing (Yamamoto, 2008). This means that the spatial continuity of the kriged surface is greater than that of the observed data.

In this context, the high spatial variability of soil properties themselves is a strong justification to admit that the results of productivity, water demand for irrigation, evapotranspiration rates, etc. will also vary in space, demonstrating the importance of stochastic models, such as the geostatistical simulation.

Several studies have used geostatistical methods to evaluate the spatial variability of soil physical and hydraulic attributes (Yamamoto, 2008; Santos et al., 2012; Furtunato et al., 2013; Dec & Dörner, 2014; Melo, 2015; Xu et al., 2017). Hu et al. (2007) used the HYDRUS-1D model (Šimůnek et al., 1998) to simulate the movement of water through the soil profile, but applying it to different random fields generated by sequential Gaussian simulation (SGS). Their results justified the use of a stochastic approach that takes into account the soil variability, as well as the associated uncertainties.

Based on the foregoing, this study aimed to perform a stochastic analysis of the saturated hydraulic conductivity (K_sat) and the parameters of the soil water retention curve (RC) in the Donato stream basin, RS, Brazil, presented and discussed in Part 1 (Siqueira et al. 2019).

Material and Methods

The study area comprises the Donato stream basin, located in the municipality of Pejuçara, RS, Brazil, with an area of 1.10 km² and geographic coordinates between 28º 25’ 34” S and 53º 40’ 30” W, and 28º 24’ 50” S and 53º 41’ 30” W, 590 m of altitude.

A regular grid was established covering the whole area of this basin. The soil was sampled at the crossing points of the grid, with regular spacing of 140 m and the most distant samples were collected in a regular spacing of 200 m, making a total of 55 points. These samples were used to determine the soil water retention curve, and its parameters as given by Genuchten (1980), and the saturated hydraulic conductivity. For more information about the study area and soil sampling, see Part 1 (Siqueira et al., 2019) of this research.

The analysis of the spatial variability of the saturated hydraulic conductivity and the soil water retention curve presented in Part 1 (Siqueira et al. 2019) was done by means of geostatistical simulation using the sequential Gaussian simulation (SGS) method.

Geostatistical simulation algorithms assume that a value of the variable will be randomly extracted from the conditional cumulative distribution function (ccdf) at each point in the area, obtained from the kriging mean and variance at that point (Goovaerts, 1997; Deustsch & Journel, 1998):

$$F\left(u;z|\left(m\right)\right)=\text{Prob}\left(Z\left(u\right)\le z|\left(m\right)\right)$$ (1)

where:

u - location of the variable Z(u);

z - cut value of the ccdf; and,

m - information conditioned to the simulated values.

Likewise, different simulated maps (random fields) can be generated together from the random sampling of the ccdfs at all points of the grid:

$$F\left(u_1^{\text{'}},u_2^{\text{'}}\mathrm{,...,}{u}_N^{\text{'}};z_1\mathrm{,...,}{z}_N|\left(m\right)\right)=\text{Prob}\left(Z\left(u_1^{\text{'}}\right)\le z_{\mathrm{1,...,}}Z\left(u_N^{\text{'}}\right)\le Z_N|\left(m\right)\right)$$ (2)

where:

N - number of grid nodes;

u’_i - variable location where there is no observed data;

F - conditional cumulative distribution function of the variable Z at each point u_α, for α = 1, 2, ..., N; and,

Prob - probability of non-exceedance.

The geostatistical softwares SGeMS (Remy et al., 2009) and GSLIB (Deutsch & Journel, 1998) were used, both of public domain.

Initially, all variables were normalized as:

$$Y\left(u\right)=\Phi \left(Z\left(u\right)\right)$$ (3)

where:

Z(u) - variable of interest;

Y(u) - normalized variable; and,

Φ - normalization function.

Then, the directional variograms of the normalized data of each variable were modeled in 8 directions, from 0 to 157.5º, varying by 22.5º. The equation for obtaining the experimental variogram is given by Goovaerts (1997).

$$\gamma \left(H\right)=\frac{1}{2N\left(H\right)}{\displaystyle \sum _{\alpha =1}^{N\left(H\right)}{\left[z\left(u_{\alpha }\right)-z\left(u_{\alpha }+H\right)\right]}^2}$$ (4)

where:

γ(H) - variance function at distance H;

H - distance between pairs of samples;

N(H) - number of pairs of samples separated by the distance H; and,

z(u_α) and z(u_α + H) - samples that compose the pairs separated by H.

Based on the variograms, a block model was generated for the entire area and for each variable through the sequential Gaussian simulation method. The hydraulic conductivity data and the soil water retention curve parameters were initially normalized.

The geostatistical simulation was carried out for a grid of 1500 x 1300 m, considering blocks of 30 x 30 m. A hundred realizations (random fields) of simulated values were generated for each variable in the study area. The simulated values were then back transformed from the Gaussian space into the original sample space.

The validation of the geostatistical simulation performed for each hydraulic variable was conducted respecting the following criteria: a) all generated maps (random fields) generated by the simulation should honor the values of the variable at m sampled locations; b) the histograms of the simulated values should be very similar to the histograms of the original data; and c) the spatial variability model of the variable (variogram) must also be reproduced by the various random fields.

Results and Discussion

Table 1 presents the parameters of the variograms in the directions of greatest and least spatial continuity of K_sat and the parameters of the RC. For all variograms the best fit was obtained by the spherical model.

Table 1
Variogram model parameters of normalized data for all variables

Variable	Continuity	Azimuth (º)	Nugget effect(1)	Sill(1)	Range (m)
ln(Ksat + 1)	Large	67.5	0.050	0.950	462
	Small	157.5			158
n	Large	157.5	0.200	0.800	518
	Small	67.5			196
ln(α)	Large	135.0	0.050	0.950	620
	Small	45.0			158
θsat	Large	112.5	0.125	0.875	273
	Small	22.5			220

• ⁽¹⁾ In this case, the sum of the nugget effect and the sill should be 1 (total variance) due to the data normalization

The variograms of ln (α) and n showed the largest ranges, while the lowest range was obtained for the θ_sat variogram. The number of pairs in the experimental variogram is considered insufficient to adequately represent spatial variation over long distances (longer than 1000 m). Hu et al. (2008) also did not find good spatial dependence for α and K_sat in their study for steep soils, and this may also be a reason for the difficulty of the present study, whose study area is moderately undulating, resulting in a hydrographic basin of only 1 km². In the same way, Paleologos & Sarris (2011) did not obtain good fits of the experimental variograms of K_sat.

In general, the K_sat ranges found in the literature (Hu et al., 2008; Santos et al., 2012) are smaller than those presented here due to the differences in sample spacing and probably due to some scale effect, as discussed by Blöschl & Sivapalan (1995). These authors comment that spatial resolutions are poorer than when working on time scales, and this is one of the biggest obstacles to handling in the field of sampling. However, Deb & Shukla (2012) commented that the hydraulic conductivity presents great spatial variability, but it can result in both short and long spatial ranges.

The parameters of the variograms demonstrate different spatial behaviors of the variables, reflected by different ranges and nugget effects of the fitted models. In addition, the directional variograms of all variables revealed that there is geometric anisotropy of the horizontal variation in the study area.

The different soil covers may also have influenced this process, since they varied during each sampling. One possibility would be to perform samplings throughout the area during the permanence of each soil cover and to compare them individually. However, this was not possible in practice due to the amount of samples collected and the time required to obtain the laboratory soil water retention curves. It is believed that this new configuration would favor the modeling of variograms, since different data populations would be individually considered and more correlated in space.

Sampling errors are also common in this type of variables, since, for example, the samples collected in cylinders can be detached from them and create preferential flow paths, masking the actual value of the variable obtained in laboratory tests, as in the case of determining K_sat with permeameters. This may have been the justification for the high K_sat values obtained.

One map of the 100 random fields of each variable is shown in Figure 1, where it is observed that there is no smoothing effect of kriging, represented by the abrupt color changes from one block to the other.

Figure 1
Maps of some simulations performed by the sequential Gaussian simulation method: random field #84 of ln(K_sat + 1) (A), random field #5 of ln(α) (B), random field #57 of n (C) and random field #29 of θ_sat (D)

As the first validation criterion of the simulations performed, it is assumed that the generated fields (Figure 1) are in agreement with the sampled values and their locations [see Figure 3, Part 1 (Siqueira et al., 2019)]. In this case, it is noted that the high and low values were well represented by the simulation. Figure 2 shows the histograms of the simulated data already transformed into the original sample space of some simulated random fields.

Figure 2
Comparison of observed (left) and simulated (right) data histograms for: ln(K_sat + 1) (A) and random field #84 of ln(K_sat + 1) (B), ln(α) (C) and random field #5 of ln(α) (D), n (E) and random field #57 of n (F), θ_sat (G) and random field #29 of θ_sat (H)

Visually, the shapes of the histograms are very much similar, as well as their statistics. Taking these randomly chosen fields as examples, it is seen that the simulated data statistics are very similar to those of the original data. This is another strong indication that the simulations performed are satisfactory. Considering the statistics of all simulations (all random fields), the results obtained are presented in Table 2.

Table 2
Statistics of all random fields for each back transformed variable

Variable	Mean			Variance
	Minimum	Maximum	All simulations	Minimum	Maximum	All simulations
Ksat (cm d-1)	10.935	33.677	22.931	1705.895	7418.207	3888.108
n (dimensionless)	1.103	1.174	1.1382	0.013	0.023	0.018
α (cm-1)	0.550	1.240	0.8248	1.130	2.530	2.105
θsat (cm3 cm-3)	0.440	0.470	0.4570	0.006	0.009	0.008

These results (illustrated in histograms of Figures 2B, D, F and H), when compared with the statistics presented in Figures 2A, C, E and G, indicate that the statistics of the parameter θ_sat were more accurately reproduced by the 100 simulated random fields, without any smoothing of the mean or change in the data variance. Similar behavior was observed for the random fields statistics of the parameter n. However, for K_sat and the parameter α, significant smoothing of the averages occurred, mainly for the first, which may be due to the existence of large areas in the grid without samples and also to the difficulties in modeling the variograms.

Similar results for K_sat were found by Jang & Liu (2004), whose simulations revealed high spatial variability of K_sat due to the wide variation of the mean and the variance between the simulations obtained. However, geostatistical simulation, on the contrary of kriging, does not seek to minimize the error of estimates (local). Instead, its objective is to reproduce the characteristics of the histograms and the spatial continuity (variogram model), besides honoring the sampled data (Goovaerts, 1999).

The good representativeness of the spatial continuity of the simulations can be checked in Figure 3. Only few simulations of each variable were chosen to facilitate the visualization of the variograms. When compared to the variograms of the original data (black line), it can be seen that the spatial continuity is maintained.

Figure 3
Comparison of original data and some simulated random fields variograms of ln(K_sat + 1) (A), ln(α) (B), n (C) and θ_sat (D)

It is observed that the simulations reproduced well the behavior of the average spatial relationship between pairs of samples, although in the case of K_sat some simulations have underestimated the original data variance, which is in agreement with a smaller variance of all the simulations (Table 2).

As discussed by Yamamoto (2008), the Gaussian simulation guarantees the reproduction of the histograms and semivariograms of the transformed data (normalized or even indicators) when compared with the histograms and semivariograms of the original data. However, the transformation of the simulated data to the original scale will not necessarily be reproduced.

With the use of an exhaustively sampled database, Yamamoto (2008) has shown that, for data with lognormal distribution (as is the case of K_sat), this difference in semivariograms can be even more pronounced. Goovaerts (1997) also comments that the simulation algorithms reproduce the variographic models only on the average of all the realizations (random fields) and that the greatest fluctuations are expected when using the sequential simulation, as used in this study.

As a final parameter to measure the uncertainties of the variables considered here, were chosen the random fields related to the percentiles of 5 and 95% of non-exceedance of the mean values of all random fields. This choice served to represent a more favorable and a more unfavorable scenario of available water in the soil and, consequently, to the development of the plant. Table 3 presents the fields whose averages refer to these percentiles for each variable. The transformations performed to facilitate the simulation of K_sat and α were removed to determine the percentiles, and thus, to obtain the final results of this study.

Table 3
Random fields associated with the percentiles of 5 and 95% of non-exceedance of the means of each variable

Variable	Percentile 5%a	Simulation number	Percentile 95%	Simulation number
Ksat (cm d-1)	24.041	32	44.601	65
n (dimensionless)	1.118	61	1.162	21
α (cm-1)	0.622	45	1.122	28
θsat (cm3 cm-3)	0.446	40	0.470	47

• ^a The means of each simulation (random field) were calculated and the percentiles of all means were obtained

Among the analyzed variables, K_sat and α present the greatest uncertainties, represented by the amplitude between the percentiles. This means that if the entire study area could be sampled with respect to these parameters, it would be possible to find, for example, a situation of K_sat mean values of both 24 and 44 cm d^-1.

These values would represent significant differences of the estimates of other hydraulic or hydrological variables which depend on the parameters analyzed here, such as surface drainage, groundwater drains design for agricultural drainage, availability of water for irrigation, etc. For example, although the effect of K_sat on soil water flow is not direct, this variable is indispensable to estimate the hydraulic conductivity of the unsaturated soil. In addition, K_sat is one of the parameters that are usually calibrated in surface hydrological models, since most of them present great sensitivity to variations in this soil property.

Conclusions

1. The random fields generated by the stochastic method of sequential Gaussian simulation were validated according to the criteria related to the good representativeness of the sample spatial continuity and the reproducibility of their variograms.

2. All simulations performed allowed an evaluation of the global uncertainty of the analyzed variables, which proved to be significant when evaluated by the amplitude between the percentiles, mainly for the saturated hydraulic conductivity and for the parameter α.

3. The difficulties of modeling the variograms may be due to factors such as sample spacing, different soil covers, sampling errors, as well as lower spatial connectivity of some variables.

Literature Cited

• Blöschl, G.; Sivapalan, M. Scale issues in hydrological modeling: A review. Hydrological Processes, v.9, p.251-290, 1995. https://doi.org/10.1002/hyp.3360090305
  » https://doi.org/10.1002/hyp.3360090305
• Deb, S. K.; Shukla, M. K. Variability of hydraulic conductivity due to multiple factors. American Journal of Environmental Science, v.8, p.489-502, 2012. https://doi.org/10.3844/ajessp.2012.489.502
  » https://doi.org/10.3844/ajessp.2012.489.502
• Dec, D.; Dörner, J. Spatial variability of the hydraulic properties of a drip irrigated andisol under blueberries. Journal of Soil Science and Plant Nutrition, v.14, p.589-601, 2014. https://doi.org/10.4067/S0718-95162014005000047
  » https://doi.org/10.4067/S0718-95162014005000047
• Deutsch, C. V.; Journel, A. G. GSLIB: Geostatistical Software Library and User’s Guide. 2.ed. Oxford: Oxford University Press, 1998. 369p.
• Furtunato, O. M.; Montenegro, S. M. G. L.; Antonino, A. C. D.; Oliveira, L. M. M. de; Souza, E. S. de; Moura, A. E. S. S. de. Variabilidade espacial dos atributos físico-hídricos de solos em uma bacia experimental no Estado de Pernambuco. Revista Brasileira de Recursos Hídricos, v.18, p.135-147, 2013. https://doi.org/10.21168/rbrh.v18n2.p135-147
  » https://doi.org/10.21168/rbrh.v18n2.p135-147
• Genuchten, M. T. van. A closed form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal, v.44, p.892-898, 1980. https://doi.org/10.2136/sssaj1980.03615995004400050002x
  » https://doi.org/10.2136/sssaj1980.03615995004400050002x
• Goovaerts, P. Geostatistics for natural resources evaluation. Oxford: Oxford University Press , 1997. 483p.
• Goovaerts, P. Geostatistics in soil science: State-of-the-art and perspectives. Geoderma, v.89, p.1-45, 1999. https://doi.org/10.1016/S0016-7061(98)00078-0
  » https://doi.org/10.1016/S0016-7061(98)00078-0
• Hu, K.; White, R. E.; Chen, D.; Li, B.; Li, W. Stochastic simulation of water drainage at the field scale and its application to irrigation management. Agricultural Water Management, v.89, p.123-130, 2007. https://doi.org/10.1016/j.agwat.2006.12.010
  » https://doi.org/10.1016/j.agwat.2006.12.010
• Hu, W.; Shao, M. A.; Wang, Q. J.; Fan, J.; Reichardt, K. Spatial variability of soil hydraulic properties on a steep slope in the loess plateau of China. Scientia Agricola, v.65, p.268-276, 2008. https://doi.org/10.1590/S0103-90162008000300007
  » https://doi.org/10.1590/S0103-90162008000300007
• Jang, C.-S.; Liu, C.-W. Geostatistical analysis and conditional simulation for estimating the spatial variability of hydraulic conductivity in the Choushui River alluvial fan, Taiwan. Hydrological Processes, v.18, p.1333-1350, 2004. https://doi.org/10.1002/hyp.1397
  » https://doi.org/10.1002/hyp.1397
• Melo, T. M. de. Simulação estocástica dos impactos das mudanças climáticas sobre as demandas de água para irrigação na região Noroeste do Rio Grande do Sul. Porto Alegre: UFRGS, 2015. 133p. Tese Doutorado
• Paleologos, E. K.; Sarris, T. S. Stochastic analysis of flux and head moments in a heterogeneous aquifer system. Stochastic Environmental Research and Risk Assessment, v.25, p.747-759, 2011. https://doi.org/10.1007/s00477-011-0459-7
  » https://doi.org/10.1007/s00477-011-0459-7
• Remy, N.; Boucher, A.; Wu, J. Applied geostatistics with SGeMS: A user's guide. Cambridge: Cambridge University Press, 2009. 284p. https://doi.org/10.1017/CBO9781139150019
  » https://doi.org/10.1017/CBO9781139150019
• Santos, K. S.; Montenegro, A. A. A.; Almeida, B. G. de; Montenegro, S. M. G. L.; Andrade, T. da S.; Fontes Júnior, R. V. de P. Variabilidade espacial de atributos físicos em solos de vale aluvial no semiárido de Pernambuco. Revista Brasileira de Engenharia Agrícola e Ambiental, v.16, p.828-835, 2012. https://doi.org/10.1590/S1415-43662012000800003
  » https://doi.org/10.1590/S1415-43662012000800003
• Šimůnek, J.; Huang, K.; Genuchten, M. T. van. The HYDRUS code for simulating the one-dimensional movement of water, heat and multiple solutes in variably-saturated media. Riverside: U.S. Salinity Laboratory, 1998. 164p.
• Siqueira, T. M.; Louzada, J. A. S.; Pedrollo, O. C.; Castro, N. M. dos R.; Oliveira, M. H. C. de. Soil physical and hydraulic properties in the Donato stream basin, RS, Brazil. Part 1: Spatial variability. Revista Brasileira de Engenharia Agrícola e Ambiental , v.23, p.669-674, 2019. http://dx.doi.org/10.1590/1807-1929/agriambi.v23n9p669-674
  » http://dx.doi.org/10.1590/1807-1929/agriambi.v23n9p669-674
• Xu, Z.; Wang, X; Chai, J.; Qin, Y.; Li, Y. Simulation of the spatial distribution of hydraulic conductivity in porous media through different methods. Mathematical Problems in Engineering, v.2017, p.1-10, 2017. https://doi.org/10.1155/2017/4321918
  » https://doi.org/10.1155/2017/4321918
• Yamamoto, J. K. Estimation or simulation? That is the question. Computational Geosciences, v.12, p.573-591, 2008. https://doi.org/10.1007/s10596-008-9096-8
  » https://doi.org/10.1007/s10596-008-9096-8

Publication Dates

• Publication in this collection
  12 Aug 2019
• Date of issue
  Sept 2019

History

• Received
  19 June 2018
• Accepted
  16 July 2019
• Published
  31 July 2019