Current generation LSMs have several shortcomings, including the lack of robust numerical methods to solve discretized equations, monolithic software design that hinders testing of a process representation in isolation from other components, and absence of a flexible coupling framework to solve tightly couple multi-physics problems. While the LSM community vigorously engages in accessing the accuracy of the conceptual model to represent the real system (validation) through numerous Model Intercomparison Projects, the fidelity of the numerical implementation of the conceptual model (verification) is seldom evaluated. In this work, we develop a numerically robust standalone Multi-Physics Problem (MPP) library for solving terrestrial processes at global scales with support for flexible coupling strategies. The MPP library uses the Method of Manufactured Solutions (MMS) to verify the numerical implementation when the exact solution is unknown. The underlying numerical engine of the MPP library is PETSc, whose DMComposite subclass provides an interface for assembling individual physics components as the part of a global matrix for solving tightly coupled multi-physics problems. Physics formulations supported by the MPP library include the variably saturated Richards equation and a thermal advection-diffusion model. The MMS technique was applied to verify the MPP library for a range of single and multi-physics problems that were solved in one or three dimensions. The numerical library solves variably saturated flow that includes coupled soil and plant (root, xylem) hydraulics and a thermal advection-diffusion model. For all problems, the error between the simulated and manufactured solution reduced as mesh size decreased and expected order of accuracy metric was obtained. The MPP library has been coupled with the E3SM Land Model.