Summary Reducing the gap between geophysical inversion and geological interpretation can be achieved by integrating geological modelling into geophysical inversion. For this, we use a generalised, iterative level-set gravity inversion scheme in which geological units are deformed automatically. During the inversion process, a regularisation term is defined using automated geological modelling to account for geological data and principles. This provides model-dependent geological constraints and encourages geological realism throughout inversion. To alleviate the dependence on the starting model and consider the possibility of features unseen by direct observations, an automated geophysical data-driven method is proposed to insert new rock units in the model. Uncertainty quantification is achieved through the null space shuttle algorithm, which is used to generate a series of alternative models that are consistent with geophysical data. This methodology is applied to assess the uncertainties of a pre-existing 3D crustal-scale geological model of the Western Pyrenean orogeny (France, Spain). The area is characterized by a positive gravity anomaly generally attributed to the presence of a shallow mantle body. The impact of variations in shape and density of key crustal and mantle features is investigated. Different scenarios are explored in 3D space to produce a range of viable, relatively simple crustal scale models of the area. This application demonstrates the capability and potential of this approach to evaluate alternative interpretations of geophysical data. The results show the plausibility of scenarios with a shorter subducted Iberian lower crust and a denser Axial Zone than in the pre-existing model.