Subsurface multiphysics simulations are challenging because they involve coupled flow, transport, mechanics, heat transfer, heterogeneous permeability (uncertain parameters), fracture evolution, and long-time dynamics. Reliable computation for these systems requires more than accuracy alone: numerical methods must preserve local conservation, respect physical bounds, handle uncertainty, and remain efficient for nonlinear high-dimensional models. This talk presents a unified research program on structure-preserving and data-driven methods for subsurface applications. I will discuss enriched Galerkin finite element methods for thermo-poro-elasticity systems with local mass conservation; flux-corrected transport methods for maximum-principle preservation; data assimilation methods for Biot and two-phase Darcy flow; and scientific machine learning approaches for inverse permeability reconstruction. I will also describe recent work on heterogeneous permeability modeling, phase-field fracture, and neural-network extrapolation of long-time dynamics. The central theme is that robust subsurface computation benefits from combining partial differential equations structure, numerical analysis, and data-driven modeling. Thus,  the unifying message is not ‘classical numerics versus SciML’; it is how to combine them so that each compensates for the other’s weaknesses.