Precision landings on the lunar surface are essential for accessing scientifically valuable sites and enabling future lunar base operations. Achieving such accuracy demands high-fidelity representations of the spacecraft dynamics during early GNC development phases. This work addresses two major challenges in modeling lunar lander dynamics: 1) the accurate capturing of variable-mass dynamics due to rapid propellant consumption, and 2) the characterization of sloshing-induced perturbations, which are strongly dependent on the fill level and the lander's acceleration. Expanding upon a rigid body dynamical model, in this paper we derive a set of nonlinear equations of motion that incorporates both variable-mass effects and sloshing dynamics using a pendulum-based approximation parameterized for spherical tanks. The formulation is modular and scalable, allowing a straightforward integration of multiple tanks and propulsion systems. The model’s validity is demonstrated through a simplified lunar lander scenario showing very close alignment with published data on sloshing mode characteristics.