Obtaining reliable estimates of the subsurface state distributions is essential to determine the location of, e.g., potential nuclear waste disposal sites. However, providing these is challenging since it requires solving the problem numerous times, yielding high computational cost. To overcome this, we use a physics-based machine learning method to construct surrogate models. We demonstrate how it produces physics-preserving predictions, which differentiates it from purely data-driven approaches.