To reduce uncertainties in reconstructed images, geologic information must be introduced in a numerically robust and stable way during the geophysical data inversion procedure. In the context of potential (gravity) data inversion, it is important to bound the physical properties by providing probabilistic information on the number of lithologies and ranges of values of possibly existing related rock properties (densities). For this purpose, we have introduced a generalization of bounding constraints for geophysical inversion based on the alternating direction method of multipliers. The flexibility of the proposed technique enables us to take into account petrophysical information as well as probabilistic geologic modeling, when it is available. The algorithm introduces a priori knowledge in terms of physically acceptable bounds of model parameters based on the nature of the modeled lithofacies in the region under study. Instead of introducing only one interval of geologically acceptable values for each parameter representing a set of rock properties, we define sets of disjoint intervals using the available geologic information. Different sets of intervals are tested, such as quasidiscrete (or narrow) intervals as well as wider intervals provided by geologic information obtained from probabilistic geologic modeling. Narrower intervals can be used as soft constraints encouraging quasidiscrete inversions. The algorithm is first applied to a synthetic 2D case for proof-of-concept validation and then to the 3D inversion of gravity data collected in the Yerrida Basin (Western Australia). Numerical convergence tests show the robustness and stability of the bound constraints that we apply, which is not always trivial for constrained inversions. This technique can be a more reliable uncertainty reduction method as well as an alternative to other petrophysically or geologically constrained inversions based on the more classic “clustering” or Gaussian-mixture approaches.