SUMMARY Using forward mantle convection models starting at 140 Ma, and assimilating plate reconstructions as surface velocity boundary condition, we predict present-day mantle structure and compare them with tomography models, using geoid as an additional constraint. We explore a wide model parameter space, such as different values of Clapeyron slope and density change across 660 km, density and viscosity of the thermochemical piles at the core–mantle boundary (CMB), internal heat generation rate, and model initiation age. We also investigate the effects of different strengths of a weak layer below 660 km and weaker asthenosphere and slabs. Our results suggest that slab structures at different subduction zones are sensitive to the viscosity of the asthenosphere, strength of slabs, values of Clapeyron slope and the density and viscosity of the thermochemical piles, while different internal heat generation rates do not affect the slab structures. We find that with a moderately weak asthenosphere ($10^{20}$ Pa·s) and strong slabs, the predicted slab structures are consistent with the tomography models, and the observed geoid is also matched well. Moreover, our models successfully reproduce the degree-2 structure of the lower mantle beneath Africa and the Pacific, also known as Large Low Shear Velocity provinces (LLSVPs). A moderate Clapeyron slope of −2.5 MPa K−1 at 660 km aids in slab stagnation while higher values result in massive slab accumulation at that depth, ultimately leading to slab avalanches. We also find that the convective patterns in the thermal and thermochemical cases with slightly denser LLSVPs are similar, although the geoid amplitudes are lower for the latter. However, with more dense LLSVPs, the slabs cannot perturb them and no plumes are generated. Plumes arise as thermal instabilities from the edges of the LLSVPs, when cold and viscous slabs perturb them. While our predicted plume locations are consistent with the observed hotspot locations, matching the plume structures in tomography models is difficult. These plumes are essential in fitting the finer features of the observed geoid. In longer-duration models, more voluminous subducted material reaches the CMB, which tends to erode the LLSVPs significantly, and yields a poor fit to the observed geoid. Our results suggest that with the presence of a thin, moderately weak layer below 660 km, a slightly dense LLSVP, and Clapeyron slope of −2.5 MPa K−1, the velocity anomalies in seismic tomography and the long-wavelength geoid can be matched well. One of the limitations of our models is that the assimilated plate motion history may be too short to overcome arbitrary initial conditions effects. Also, assimilated true plate velocities in our models may not represent the true convective vigour of the Earth.