Over the last two decades, the application of Markov Chain Monte Carlo (MCMC) to Multi-Observable Thermochemical Tomography (MTT) has become increasingly popular for inferring the thermochemical structure of the Earth, thanks to its ability to sample the highly nonlinear posterior probability distribution functions (PDFs) that lack closed-form expressions [1]. The computational overburden of this stochastic approach is exceedingly demanding, as it entails performing many evaluations of complex forward problems, potentially limiting its applicability to real-world scale studies. The challenge is further complicated when attempting to incorporate Stokes-like models (i.e., mantle convection and lithospheric removal) related to the present-day dynamic topography (PDDT), which are fundamental to modeling the inherent nonlinear dynamic effects. For these reasons, standard numerical schemes are often inefficient for solving such problems, and several Reduced-Order Models (ROMs) have been developed to obtain efficient and reliable approximations. Our work aims to address the physical phenomenon in its most general form, accurately capturing all relevant effects and drawing consistent insights at the end of the inversion process. In particular, we built upon the work of [2], introducing higher levels of complexity through relaxing the restrictive assumption of linearity (i.e., assuming non-Newtonian rheology). Subsequently, we developed a hybrid ROM methodology, which leverages a distinctive combination of the latest advancements in scientific machine learning (SciML) and tailors it to develop an efficient ROM capable of reducing the fully nonlinear 3D Stoke’s problem, ensuring that accuracy is not compromised. Furthermore, the new solver can be applied non-intrusively within the inversion process, enabling seamless model updates as the inversion progresses through the parameter space exploration.