Lithospheric density structure can be constructed from seismic tomography, gravity modelling, or using both data sets. The different approaches have their own uncertainties and limitations. This study aims to characterize and quantify some of the uncertainties in gravity modelling of lithosphere densities. To evaluate the gravity modelling we compare gravity-based and seismic velocity-based approaches to estimating lithosphere densities. In this study, we use a crustal model together with lithospheric isostasy and gravity field observations to estimate lithosphere densities. To quantify the effect of uncertainty in the crustal model, three models are implemented in this study: CRUST1.0, EuCrust-07 and a high-resolution P-wave velocity model of the British Isles and surrounding areas. Different P-wave velocity-to-density conversions are used to study the uncertainty in these conversion methods. The crustal density models are forward modelled into gravity field quantities using a method that is able to produce spherical harmonic coefficients. Deep mantle signal is assumed to be removed by removing spherical harmonic coefficients of degree 0–10 in the observed gravity field. The uncertainty in the resulting lithosphere densities due to the different crustal models is ±110 kg m−3, which is the largest uncertainty in gravity modelling. Other sources of uncertainty, such as the VP to density conversion (±10 kg m−3), long-wavelength truncation (±5 kg m−3), choice of reference model (<±20 kg m−3) and Lithosphere Asthenosphere Boundary uncertainty (±30 kg m−3), proved to be of lesser importance. The resulting lithosphere density solutions are compared to density models based on a shear wave velocity model. The comparison shows that the gravity-based models have an increased lateral resolution compared to the tomographic solutions. However, the density anomalies of the gravity-based models are three times higher. This is mainly due to the high resolution in the gravity field. To account for this, the gravity-based density models are filtered with a spatial Gaussian filter with 200 km half-width, which results in similar density estimates (±35 kg m−3) with the tomographic approach. Lastly, the gravity-based density is used to estimate laterally varying conversion factors, which correlate with major tectonic regions. The independent gravity-based solutions could help in identifying different compositional domains in the lithosphere, when compared to the tomographic solutions.@en