Decorrelation between s‐wave and bulk sound velocities in the lowermost mantle and explicit density models based on seismic tomography give evidence for non‐thermal lateral density variations in the lowermost mantle. Here we implement such variations in a numerical model of mantle flow driven by density anomalies which are derived from seismic tomography. In the lowermost ≈300 km, we either use both s‐wave and bulk sound velocities to infer both thermal mantle density anomalies and non‐thermal heterogeneities, or we assume that regions with s‐wave speed anomalies below −1% have an additional non‐thermal density anomaly. We predefine the shape of the viscosity profile in the upper mantle, transition zone and lower mantle, but absolute values of viscosity in the lithosphere, upper mantle, transition zone and lower mantle are free parameters. We use geoid, radial heat flux profile, viscosity “Haskell” average, core‐mantle boundary (CMB) excess ellipticity, as well as (optionally) long‐wavelength root mean square (RMS) CMB topography and reliable point estimates of CMB topography as constraints to optimize the model in parameter space in a least squares sense. We are able to obtain a reasonable fit to all data constraints. Computed RMS CMB topography is predominantly long‐wavelength and with 1–1.5 km RMS amplitude somewhat larger than the long‐wavelength component inferred from seismology. Geoid variance reduction is 75 to 83% in our preferred parameter range. Best fit models have a viscosity maximum close to 1023 Pas about 600 km above the CMB, and a viscosity drop near the base of the mantle, corresponding to a thermal boundary layer about 300 km thick with temperature increase from ≈2500 to 3500–4000 K.