Interpolation and inverse modeling problems are ubiquitous in environmental sciences. In many applications, the parameters being estimated or mapped have physical constraints, such as nonnegativity (e.g. concentration, hydraulic conductivity), solubility limits, censored data (e.g. due to dry wells or detection limits), and other physical boundaries or missing data. Geostatistical interpolation and inverse modeling techniques have often been applied for estimating such parameters, but these methods typically cannot enforce physical constraints. This paper describes a statistically rigorous and computationally efficient Gibbs sampler, a Markov chain Monte Carlo technique, based on an a priori truncated Gaussian distribution model, which allows for multiple and variable physical constraints to be enforced within a geostatistical framework. Sample interpolation and inverse modeling applications confirm that estimates, uncertainty bounds and conditional simulations reflect the specified constraints, leading to conclusions that are more consistent with the underlying conceptual model, and provide a more accurate measure of the posterior uncertainty of the parameters being estimated. In addition, especially in inverse modeling applications, a posteriori confidence bounds are narrower even in areas where constraints are not imposed. The method is applicable in multiple dimensions, for data with or without measurement error, and with any variogram model.