Physical models in the form of partial differential equations represent an important prior for many under-constrained problems.One example of such uses is tumor treatment planning, which heavily depends on an accurate estimate of the spatial distribution of tumor cells in a patient’s anatomy. Medical imaging scans can identify the bulk of the tumor, but they cannot reveal its full spatial distribution. Tumor cells at low concentrations remain undetectable, for example, in the most frequent type of primary brain tumors, glioblastoma. Deep-learning based approaches fail to estimate the complete tumor cell distribution due to a lack of reliable training data. Most existing works therefore rely on physics-based simulations to match observed tumors, providing anatomically and physiologically plausible estimations. However, these approaches struggle with complex and unknown initial conditions and are limited by overly rigid physical models. In this work, we present a novel method that allows balancing data-driven and physics-based cost functions. In particular, we propose a unique discretization scheme that allows quantifying the adherence of our learned spatiotemporal tumor and brain tissue distributions to their corresponding growth and elasticity equations. This quantification, serving as a regularization term rather than a hard constraint, enables greater flexibility and proficiency in assimilating patient data than existing models. We demonstrate improved coverage of tumor recurrence areas compared to existing techniques on real-world data from a cohort of patients. The method brings the potential to enhance clinical adoption of model-driven treatment planning for glioblastoma.