The sensitivity of gravity and magnetic data to deep structures and the broad availability of regional data sets and surveys of high resolution make them suitable for determining detailed three-dimensional (3D) models of the subsurface. However, the sole consideration of gravity and magnetic information cannot properly resolve heterogeneous 3D environments. Advocated to solve this problem, we present an automated refinement technique for three-dimensional multilayer models as conditioned by gravity and magnetic data and by meaningful geometrical and physical constraints. We construct our model by an aggregate of rectangular prisms and aim to estimate their bottom depths, which define the geological layers. We summarize mathematically our concept of refinement in an objective function that includes the misfit to the data, the similitude to an a priori geological-geophysical model, and the smoothness of the relief of the layers. Importantly, our objective function also includes inequality constraints that prevent the superposition of layers and integrate the surface and borehole geology with the multilayer deep model. The objective function is solved using quadratic programming in a stable iterative scheme. The resulting algorithm is tested on synthetic data and applied to crustal and sedimentary basin environments from southern Baja California, Mexico. The assimilation of the geological and geometrical constraints to the inversion process produces models that correlate with the surface geology and reveal the three-dimensional features of the subsurface. (c) 2004 Elsevier B.V. All rights reserved.