Space-time mapping of soil salinity using probabilistic bayesian maximum entropy

The mapping of saline soils is the first task before any reclamation effort. Reclamation is based on the knowledge of soil salinity in space and how it evolves with time. Soil salinity is traditionally determined by soil sampling and laboratory analysis. Recently, it became possible to complement these hard data with soft secondary data made available using field sensors like electrode probes. In this study, we had two data sets. The first includes measurements of field salinity (EC_a) at 413 locations and 19 time instants. The second, which is a subset of the first (13 to 20 locations), contains, in addition to EC_a, salinity determined in the laboratory (EC_2.5). Based on a procedure of cross-validation, we compared the prediction performance in the space-time domain of 3 methods: kriging using either only hard data (HK) or hard and mid interval soft data (HMIK), and Bayesian maximum entropy (BME) using probabilistic soft data. We found that BME was less biased, more accurate and giving estimates, which were better correlated with the observed values than the two kriging techniques. In addition, BME allowed one to delineate with better detail saline from non-saline areas.     

BME prediction of continuous geographical properties using auxiliary variables

Using auxiliary information to improve the prediction accuracy of soil properties in a physically meaningful and technically efficient manner has been widely recognized in pedometrics. In this paper, we explored a novel technique to effectively integrate sampling data and auxiliary environmental information, including continuous and categorical variables, within the framework of the Bayesian maximum entropy (BME) theory. Soil samples and observed auxiliary variables were combined to generate probability distributions of the predicted soil variable at unsampled points. These probability distributions served as soft data of the BME theory at the unsampled locations, and, together with the hard data (sample points) were used in spatial BME prediction. To gain practical insight, the proposed approach was implemented in a real-world case study involving a dataset of soil total nitrogen (TN) contents in the Shayang County of the Hubei Province (China). Five terrain indices, soil types, and soil texture were used as auxiliary variables to generate soft data. Spatial distribution of soil total nitrogen was predicted by BME, regression kriging (RK) with auxiliary variables, and ordinary kriging (OK). The results of the prediction techniques were compared in terms of the Pearson correlation coefficient (r), mean error (ME), and root mean squared error (RMSE). These results showed that the BME predictions were less biased and more accurate than those of the kriging techniques. In sum, the present work extended the BME approach to implement certain kinds of auxiliary information in a rigorous and efficient manner. Our findings showed that the BME prediction technique involving the transformation of variables into soft data can improve prediction accuracy considerably, compared to other techniques currently in use, like RK and OK.     

Assimilation of fuzzy data by the BME method

Modern spatiotemporal geostatistics provides a powerful framework for generation of predictive maps over a spatiotemporal domain by accounting for general knowledge to define a space of plausible events and then restricting this space of plausible events to be consistent with available site-specific knowledge. The Bayesian maximum entropy (BME) method is one of the most widely used modern geostatistics methods. BME results from assigning probabilities of plausible events based on general knowledge through information maximization and then applying operational Bayesian conditionalization that can explicitly assimilate stochastic representations of various uncertain (soft) data bases. The paper demonstrates that fuzzy data sets can be indirectly assimilated by BME through a two-step process: (a) reinterpretation of the fuzzy data as probabilistic through a generalized defuzzification procedure, and (b) efficient assimilation of the probabilistic results of generalized defuzzification by the BME method. A numerical demonstration involves site-specific probabilistic results obtained from the generalized defuzzification of a simulated fuzzy data set and general knowledge that includes the spatial mean trend and correlation structure models. The parameters of these models can be inferred from the hard data equivalent values of the probabilistic results. Accordingly, details of inference based on probabilistic soft data are also considered.     

A BME solution of the inverse problem for saturated groundwater flow

In most real-world hydrogeologic situations, natural heterogeneity and measurement errors introduce major sources of uncertainty in the solution of the inverse problem. The Bayesian Maximum Entropy (BME) method of modern geostatistics offers an efficient solution to the inverse problem by first assimilating various physical knowledge bases (hydrologic laws, water table elevation data, uncertain hydraulic resistivity measurements, etc.) and then producing robust estimates of the subsurface variables across space. We present specific methods for implementing the BME conceptual framework to solve an inverse problem involving Darcys law for subsurface flow. We illustrate one of these methods in the case of a synthetic one-dimensional case study concerned with the estimation of hydraulic resistivity conditioned on soft data and hydraulic head measurements. The BME framework processes the physical knowledge contained in Darcys law and generates accurate estimates of hydraulic resistivity across space. The optimal distribution of hard and soft data needed to minimize the associated estimation error at a specified sampling cost is determined. This work was supported by grants from the National Institute of Environmental Health Sciences (Grant no. 5 P42 ES05948 and P30ES10126), the National Aeronautics and Space Administration (Grant no. 60-00RFQ041), the Army Research Office (Grant no. DAAG55-98-1-0289), and the National Science Foundation under Agreement No. DMS-0112069.     

Space–time forecasting of groundwater level using a hybrid soft computing model

Forecasting of space–time groundwater level is important for sparsely monitored regions. Time series analysis using soft computing tools is powerful in temporal data analysis. Classical geostatistical methods provide the best estimates of spatial data. In the present work a hybrid framework for space–time groundwater level forecasting is proposed by combining a soft computing tool and a geostatistical model. Three time series forecasting models: artificial neural network, least square support vector machine and genetic programming (GP), are individually combined with the geostatistical ordinary kriging model. The experimental variogram thus obtained fits a linear combination of a nugget effect model and a power model. The efficacy of the space–time models was decided on both visual interpretation (spatial maps) and calculated error statistics. It was found that the GP–kriging space–time model gave the most satisfactory results in terms of average absolute relative error, root mean square error, normalized mean bias error and normalized root mean square error.     

Radial basis function neural network for hydrologic inversion: an appraisal with classical and spatio-temporal geostatistical techniques in the context of site characterization

This paper investigates three techniques for spatial mapping and the consequential hydrologic inversion, using hydraulic conductivity (or transmissivity) and hydraulic head as the geophysical parameters of concern. The data for the study were obtained from the Waste Isolation and Pilot Plant (WIPP) site and surrounding area in the remote Chihuahuan Desert of southeastern New Mexico. The central technique was the Radial Basis Function algorithm for an Artificial Neural Network (RBF-ANN). An appraisal of its performance in light of classical and temporal geostatistical techniques is presented. Our classical geostatistical technique of concern was Ordinary Kriging (OK), while the method of Bayesian Maximum Entropy (BME) constituted an advanced, spatio-temporal mapping technique. A fusion technique for soft or inter-dependent data was developed in this study for use with the neural network. It was observed that the RBF-ANN is capable of hydrologic inversion for transmissivity estimation with features remaining essentially similar to that obtained from kriging. The BME technique, on the other hand, was found to reveal an ability to map localized lows and highs that were otherwise not as apparent in OK or RBF-ANN techniques.     

Space–time forecasting using soft geostatistics: a case study in forecasting municipal water demand for Phoenix, Arizona

Managing environmental and social systems in the face of uncertainty requires the best possible forecasts of future conditions. We use space–time variability in historical data and projections of future population density to improve forecasting of residential water demand in the City of Phoenix, Arizona. Our future water estimates are derived using the first and second order statistical moments between a dependent variable, water use, and an independent variable, population density. The independent variable is projected at future points, and remains uncertain. We use adjusted statistical moments that cover projection errors in the independent variable, and propose a methodology to generate information-rich future estimates. These updated estimates are processed in Bayesian Maximum Entropy (BME), which produces maps of estimated water use to the year 2030. Integrating the uncertain estimates into the space–time forecasting process improves forecasting accuracy up to 43.9% over other space–time mapping methods that do not assimilate the uncertain estimates. Further validation studies reveal that BME is more accurate than co-kriging that integrates the error-free independent variable, but shows similar accuracy to kriging with measurement error that processes the uncertain estimates. Our proposed forecasting method benefits from the uncertain estimates of the future, provides up-to-date forecasts of water use, and can be adapted to other socio-economic and environmental applications.     

Accounting for the uncertainty in the local mean in spatial prediction by Bayesian Maximum Entropy

Bayesian Maximum Entropy (BME) has been successfully used in geostatistics to calculate predictions of spatial variables given some general knowledge base and sets of hard (precise) and soft (imprecise) data. This general knowledge base commonly consists of the means at each of the locations considered in the analysis, and the covariances between these locations. When the means are not known, the standard practice is to estimate them from the data; this is done by either generalized least squares or maximum likelihood. The BME prediction then treats these estimates as the general knowledge means, and ignores their uncertainty. In this paper we develop a prediction that is based on the BME method that can be used when the general knowledge consists of the covariance model only. This prediction incorporates the uncertainty in the estimated local mean. We show that in some special cases our prediction is equal to results from classical geostatistics. We investigate the differences between our approach and the standard approach for predicting in this common practical situation.     

Spatial prediction of categorical variables: the Bayesian maximum entropy approach

 Being a non-linear method based on a rigorous formalism and an efficient processing of various information sources, the Bayesian maximum entropy (BME) approach has proven to be a very powerful method in the context of continuous spatial random fields, providing much more satisfactory estimates than those obtained from traditional linear geostatistics (i.e., the various kriging techniques). This paper aims at presenting an extension of the BME formalism in the context of categorical spatial random fields. In the first part of the paper, the indicator kriging and cokriging methods are briefly presented and discussed. A special emphasis is put on their inherent limitations, both from the theoretical and practical point of view. The second part aims at presenting the theoretical developments of the BME approach for the case of categorical variables. The three-stage procedure is explained and the formulations for obtaining prior joint distributions and computing posterior conditional distributions are given for various typical cases. The last part of the paper consists in a simulation study for assessing the performance of BME over the traditional indicator (co)kriging techniques. The results of these simulations highlight the theoretical limitations of the indicator approach (negative probability estimates, probability distributions that do not sum up to one, etc.) as well as the much better performance of the BME approach. Estimates are very close to the theoretical conditional probabilities, that can be computed according to the stated simulation hypotheses.     

Defining probability of saturation with indicator kriging on hard and soft data

In humid, well-vegetated areas, such as in the northeastern US, runoff is most commonly generated from relatively small portions of the landscape becoming completely saturated, however, little is known about the spatial and temporal behavior of these saturated regions. Indicator kriging provides a way to use traditional water table data to quantify probability of saturation to evaluate predicted spatial distributions of runoff generation risk, especially for the new generation of water quality models incorporating saturation excess runoff theory. When spatial measurements of a variable are transformed to binary indicators (i.e., 1 if above a given threshold value and 0 if below) and the resulting indicator semivariogram is modeled, indicator kriging produces the probability of the measured variable to exceed the threshold value. Indicator kriging gives quantified probability of saturation or, consistent with saturation excess runoff theory, runoff generation risk with depth to water table as the variable and the threshold set near the soil surface. The probability of saturation for a 120 m × 180 m hillslope based upon 43 measurements of depth to water table is investigated with indicator semivariograms for six storm events. The indicator semivariograms show high spatial structure in saturated regions with large antecedent rainfall conditions. The temporal structure of the data is used to generate interpolated (soft) data to supplement measured (hard) data. This improved the spatial structure of the indicator semivariograms for lower antecedent rainfall conditions. Probability of saturation was evaluated through indicator kriging incorporating soft data showing, based on this preliminary study, highly connected regions of saturation as expected for the wet season (April through May) in the Catskill Mountain region of New York State. Supplementation of hard data with soft data incorporates physical hydrology of the hillslope to capture significant patterns not available when using hard data alone for indicator kriging. With the need for water quality models incorporating appropriate runoff generation risk estimates on the rise, this manner of data will lay the groundwork for future model evaluation and development.     

Compositional Bayesian indicator estimation

Indicator kriging is widely used for mapping spatial binary variables and for estimating the global and local spatial distributions of variables in geosciences. For continuous random variables, indicator kriging gives an estimate of the cumulative distribution function, for a given threshold, which is then the estimate of a probability. Like any other kriging procedure, indicator kriging provides an estimation variance that, although not often used in applications, should be taken into account as it assesses the uncertainty of the estimate. An alternative approach to indicator estimation is proposed in this paper. In this alternative approach the complete probability density function of the indicator estimate is evaluated. The procedure is described in a Bayesian framework, using a multivariate Gaussian likelihood and an a priori distribution which are both combined according to Bayes theorem in order to obtain a posterior distribution for the indicator estimate. From this posterior distribution, point estimates, interval estimates and uncertainty measures can be obtained. Among the point estimates, the median of the posterior distribution is the maximum entropy estimate because there is a fifty-fifty chance of the unknown value of the estimate being larger or smaller than the median; that is, there is maximum uncertainty in the choice between two alternatives. Thus in some sense, the latter is an indicator estimator, alternative to the kriging estimator, that includes its own uncertainty. On the other hand, the mode of the posterior distribution estimator, assuming a uniform prior, is coincidental with the simple kriging estimator. Additionally, because the indicator estimate can be considered as a two-part composition which domain of definition is the simplex, the method is extended to compositional Bayesian indicator estimation. Bayesian indicator estimation and compositional Bayesian indicator estimation are illustrated with an environmental case study in which the probability of the content of a geochemical element in soil being over a particular threshold is of interest. The computer codes and its user guides are public domain and freely available.     

Selecting best mapping strategies for storm runoff modeling in a mountainous semi‐arid area

Accurate runoff and soil erosion modeling is constrained by data availability, particularly for physically based models such as OpenLISEM that are data demanding, as the processes are calculated on a cell‐by‐cell basis. The first decision when using such models is to select mapping units that best reflect the spatial variability of the soil and hydraulic properties in the catchment. In environments with limited data, available maps are usually generic, with large units that may lump together the values of the soil properties, affecting the spatial patterns of the predictions and output values in the outlet. Conversely, the output results may be equally acceptable, following the principle of equifinality. To studyhow the mapping method selected affects the model outputs, four types of input maps with different degrees of complexity were created: average values allocated to general soil map units (ASG1), average values allocated to detailed map units (ASG2), values interpolated by ordinary kriging (OK) and interpolated by kriging with external drift (KED). The study area was Ribeira Seca, a 90 km² catchment located in Santiago Island, Cape Verde (West Africa), a semi‐arid country subject to scarce but extreme rainfall during the short tropical summer monsoon. To evaluate the influence of rainfall on runoff and erosion, two storm events with different intensity and duration were considered. OK and KED inputs produced similar results, with the latter being closer to the observed hydrographs. The highest soil losses were obtained with KED (43 ton ha^{? 1} for the strongest event). To improve the results of soil loss predictions, higher accurate spatial information on the processes is needed; however, spatial information of input soil properties alone is not enough in complex landscapes. The results demonstrate the importance of selecting the appropriate mapping strategy to obtain reliable runoff and erosion estimates. Copyright © 2013 John Wiley & Sons, Ltd.     

Estimating threshold-exceeding probability maps of environmental variables with Markov chain random fields

Estimating and mapping spatial uncertainty of environmental variables is crucial for environmental evaluation and decision making. For a continuous spatial variable, estimation of spatial uncertainty may be conducted in the form of estimating the probability of (not) exceeding a threshold value. In this paper, we introduced a Markov chain geostatistical approach for estimating threshold-exceeding probabilities. The differences of this approach compared to the conventional indicator approach lie with its nonlinear estimators—Markov chain random field models and its incorporation of interclass dependencies through transiograms. We estimated threshold-exceeding probability maps of clay layer thickness through simulation (i.e., using a number of realizations simulated by Markov chain sequential simulation) and interpolation (i.e., direct conditional probability estimation using only the indicator values of sample data), respectively. To evaluate the approach, we also estimated those probability maps using sequential indicator simulation and indicator kriging interpolation. Our results show that (i) the Markov chain approach provides an effective alternative for spatial uncertainty assessment of environmental spatial variables and the probability maps from this approach are more reasonable than those from conventional indicator geostatistics, and (ii) the probability maps estimated through sequential simulation are more realistic than those through interpolation because the latter display some uneven transitions caused by spatial structures of the sample data.     

Stability analysis in nonstationary spatial covariance estimation

Space deformation modelling and estimation techniques based on Multidimensional Scaling (MDS) methods play an important role in nonparametric approaches to the covariance structure analysis of the spatiotemporal processes underlying environmental studies. Since any related procedure depends on the planar MDS representation, the stability of the estimated dispersion, together with the determination of the most influential stations in the estimation of the dispersion space, are important issues that must be analysed before performing the final mapping. In this paper, stability analysis, both in terms of the MDS model and of the variogram function, as well as concerning the derivation of kriging interpolation estimates, is addressed using a special analytical jackknife procedure. Furthermore, the influence of each station in the solution given is assessed, thus providing relevant information regarding not only the MDS procedure but also the interpolation process and the variogram estimation of the spatial dispersion.     

The Bayesian maximum entropy method for lognormal variables

The Bayesian maximum entropy (BME) method can be used to predict the value of a spatial random field at an unsampled location given precise (hard) and imprecise (soft) data. It has mainly been used when the data are non-skewed. When the data are skewed, the method has been used by transforming the data (usually through the logarithmic transform) in order to remove the skew. The BME method is applied for the transformed variable, and the resulting posterior distribution transformed back to give a prediction of the primary variable. In this paper, we show how the implementation of the BME method that avoids the use of a transform, by including the logarithmic statistical moments in the general knowledge base, gives more appropriate results, as expected from the maximum entropy principle. We use a simple illustration to show this approach giving more intuitive results, and use simulations to compare the approaches in terms of the prediction errors. The simulations show that the BME method with the logarithmic moments in the general knowledge base reduces the errors, and we conclude that this approach is more suitable to incorporate soft data in a spatial analysis for lognormal data.     

Sequential kriging and cokriging: Two powerful geostatistical approaches

A sequential linear estimator is developed in this study to progressively incorporate new or different spatial data sets into the estimation. It begins with a classical linear estimator (i.e., kriging or cokriging) to estimate means conditioned to a given observed data set. When an additional data set becomes available, the sequential estimator improves the previous estimate by using linearly weighted sums of differences between the new data set and previous estimates at sample locations. Like the classical linear estimator, the weights used in the sequential linear estimator are derived from a system of equations that contains covariances and cross-covariances between sample locations and the location where the estimate is to be made. However, the covariances and cross-covariances are conditioned upon the previous data sets. The sequential estimator is shown to produce the best, unbiased linear estimate, and to provide the same estimates and variances as classic simple kriging or cokriging with the simultaneous use of the entire data set. However, by using data sets sequentially, this new algorithm alleviates numerical difficulties associated with the classical kriging or cokriging techniques when a large amount of data are used. It also provides a new way to incorporate additional information into a previous estimation.     

ESTIMATING INFILTRATION PARAMETERS FROM BASIC SOIL PROPERTIES

Infiltration data were collected on two rectangular grids with 25 sampling points each. Both experimental grids were located in tropical rain forest (Guyana), the first in an Arenosol area and the second in a Ferralsol field. Four different infiltration models were evaluated based on their performance in describing the infiltration data. The model parameters were estimated using non-linear optimization techniques. The infiltration behaviour in the Ferralsol was equally well described by the equations of Philip, Green–Ampt, Kostiakov and Horton. For the Arenosol, the equations of Philip, Green–Ampt and Horton were significantly better than the Kostiakov model. Basic soil properties such as textural composition (percentage sand, silt and clay), organic carbon content, dry bulk density, porosity, initial soil water content and root content were also determined for each sampling point of the two grids. The fitted infiltration parameters were then estimated based on other soil properties using multiple regression. Prior to the regression analysis, all predictor variables were transformed to normality. The regression analysis was performed using two information levels. The first information level contained only three texture fractions for the Ferralsol (sand, silt and clay) and four fractions for the Arenosol (coarse, medium and fine sand, and silt and clay). At the first information level the regression models explained up to 60% of the variability of some of the infiltration parameters for the Ferralsol field plot. At the second information level the complete textural analysis was used (nine fractions for the Ferralsol and six for the Arenosol). At the second information level a principal components analysis (PCA) was performed prior to the regression analysis to overcome the problem of multicollinearity among the predictor variables. Regression analysis was then carried out using the orthogonally transformed soil properties as the independent variables. Results for the Ferralsol data show that the parameters of the Green–Ampt and Kostiakov model were estimated relatively accurately (maximum R² = 0.76). For the Arenosol, use of the second information level together with PCA produced regression models with an R² value ranging from 0.38 to 0.68. For the Ferralsol, most of the variance was explained by the root content and organic matter content. In the Arenosol plot, the fractions medium and fine sand explained most of the observed variance.     

Classifying wave forecasts with model-based Geostatistics and the Aitchison distribution

This paper proposes a non-parametric method of classification of maps (i.e., variable fields such as wave energy maps for the Western Mediterranean Sea) into a set of D typical regimes (calm, E-, SW- or N/NW-wind dominated storms, the 4 synoptic situations more often occurring in this region). Each map in the training set is described by its values at P measurement points and one of these regime classes. A map is thus identified as a labelled point in a P-dimensional feature space, and the problem is to find a discrimination rule that may be used for attaching a classification probability to future unlabelled maps. The discriminant model proposed assumes that some log-contrasts of these classification probabilities form a Gaussian random field on the feature space. Then, available data (labelled maps of the training set) are linked to these latent probabilities through a multinomial model. This model is quite common in model-based Geostatistics and the Gaussian process classification literature. Inference is here approximated numerically using likelihood based techniques. The multinomial likelihood of labelled features is combined in a Bayesian updating with the Gaussian random field, playing the role of prior distribution. The posterior corresponds to an Aitchison distribution. Its maximum posterior estimates are obtained in two steps, exploiting several properties of this family. The first step is to obtain the mode of this distribution for labelled features, by solving a mildly non-linear system of equations. The second step is to propagate these estimates to unlabelled features, with simple kriging of log-contrasts. These inference steps can be extended via Markov-chain Monte Carlo (MCMC) sampling to a hierarchical Bayesian problem. This MCMC sampling can be improved by further exploiting the Aitchison distribution properties, though this is only outlined here. Results for the application case study suggest that E- and N/NW-dominated storms can be successfully discriminated from calm situations, but not so easily distinguished from each other.     

Estimating flow resistance in steep slope rills

Recent research recognized that the slope of 18% can be used to distinguish between the ‘gentle slope’ case and that of ‘steep slope’ for the detected differences in hydraulic variables (flow depth, velocity, Reynolds number, Froude number) and those representatives of sediment transport (flow transport capacity, actual sediment load). In this paper, using previous measurements carried out in mobile bed rills and flume experiments characterized by steep slopes (i.e., slope greater than or equal to 18%), a theoretical rill flow resistance equation to estimate the Darcy-Weisbach friction factor is tested. The main aim is to deduce a relationship between the velocity profile parameter Γ, the channel slope, the Reynolds number, the Froude number and the textural classes using a data base characterized by a wide range of hydraulic conditions, plot or flume slope (18%–84%) and textural classes (clay ranging from 3% to 71%). The obtained relationship is also tested using 47 experimental runs carried out in the present investigation with mobile bed rills incised in a 18%—sloping plot with a clay loam soil and literature data. The analysis demonstrated that: (1) the soil texture affects the estimate of the Γ parameter and the theoretical flow resistance law (Equation 25), (2) the proposed Equation (25) fits well the independent measurements of the testing data base, (3) the estimate of the Darcy-Weisbach friction factor is affected by the soil particle detachability and transportability and (4) the Darcy-Weisbach friction factor is linearly related to the rill slope.