The spatial variability of soil variables is a critical component of modeling, estimation, prediction and risk assessment in soil science. On one hand, spatial variability must be taken into account for optimal spatial interpolation (e.g., kriging) and risk assessment (e.g., evaluating the probability that the value of a given property is higher than an established threshold); on the other, spatial variability influences the output of physically based models (e.g., rainfall-runoff). As soil variables are usually known at only a small number of experimental locations, their spatial variability must be evaluated using statistical tools such as the spatial covariance (or the semivariogram), which, in turn, are modeled by a few parameters (e.g., nugget variance, sill, range). These covariance parameters, however, can only be inferred with an associated statistical uncertainty. The main objective of this paper is to show that this uncertainty can be assessed by a maximum likelihood approach to inference. Three different methods for obtaining joint confidence intervals of spatial covariance parameters are considered: (i) the Fisher information matrix, (ii) likelihood intervals, and (iii) the likelihood ratio statistic. The authors show that, when expressed as likelihood intervals, the confidence regions provided by the likelihood ratio statistic are especially suitable for applications because of its straightforward calculation. A second objective of the paper is to demonstrate how the uncertainty of the spatial covariance parameters can be included in applications. The interpolation uncertainty is included by Bayesian kriging whereas it is included in risk assessment and physically based models by geostatistical simulation. A case study of the zinc content of soils of the Swiss Jura is used to illustrate the methodology.