Many natural phenomena and processes remain poorly understood and analytically intractable. One such example is the contact dynamics between non-Newtonian materials, such as soft soils and wheels. This challenge becomes particularly critical in scenarios where the use of large wheels, tracks, or tires is constrained, as is often the case in planetary rover missions. To better anticipate a rover's behavior for different terrain conditions it may encounter, researchers have developed a wide range of terramechanical models. These models differ significantly in fidelity, computational cost, and domain of applicability. Given the lack of explicit analytical formulations in the form of differential equations, machine learning–based surrogate modeling is well-suited for approximating the behavior of complex dynamical systems without requiring. Existing studies on wheel–soil contact dynamics encompass both real-world experiments and a broad spectrum of simulation approaches, ranging from fast empirical models to computationally intensive discrete element methods. By combining information from heterogeneous data sources, from low-fidelity, inexpensive simulations to high-fidelity but costly and rare measurements or simulations, multi-fidelity machine learning models can predict both accurately and computationally efficiently. This work presents a multi-fidelity Gaussian Process surrogate model for predicting wheel–soft soil interactions, specifically the forces and torques acting on a moving wheel. In addition, the study investigates strategies for enforcing physical constraints within the machine learning model. This is achieved by embedding physical knowledge through dataset construction and experimental design, and by adopting a hybrid modeling approach in which the machine learning component is restricted to learning discrepancies between different fidelity levels. The proposed approach yields surrogate models that exhibit physically consistent behavior in free-running dynamical simulations and outperform traditional terramechanical models in terms of the accuracy–computational efficiency trade-off. Furthermore, the inclusion of uncertainty quantification enables the assessment of model robustness and enhances the interpretability of the predictions.