A new framework is developed for efficient implementation of semi-empirical terramechanics models in multibody dynamics environments. In this approach, for every wheel in contact with soft soil, unilateral contact constraints are added for both the normal direction and the tangent plane. The forces associated with the latter, like traction and rolling resistance, are formulated in this approach as set-valued force laws, their properties being determined by deregularization of the terramechanics relations. As shown in the paper, this leads to the dynamics representation in the form of a linear complementarity problem (LCP). With this formulation, stable simulation of rovers is achieved even at relatively large time steps. In addition, a high-resolution height-field (HF) is employed to model terrain-surface deformation and changes in hardening of soil under the wheel. As a result, the multipass effect is captured in the presented approach. In addition, an extensive set of experiments was conducted using a version of the Juno rover (Juno II). The experimental results are analyzed and compared with the model developed in the paper.