In recent years, the application of machine learning methods in the derivation of physical governing equations has gained significant attention. This has become increasingly relevant due to the growing complexity of problems that are difficult to fully comprehend. Instead of driving solely by data, this study incorporated the conservation of mass into its framework. To ensure the physical rationality of the derived equations, dimensional analysis was incorporated into the algorithm. This facilitated establishing connections between physical parameters and each term of the target equation, ensuring the validity and reliability of the resulting equation. To enhance the interpretability of the resulting partial differential equations (PDEs), we analyzed and compared the results obtained from sparse regression, multi-objective optimization, and then proposed a sequential identification method, namely PHY-PDE. To validate this approach, the identified PDEs were rigorously tested against groundwater-related equations, specifically the Darcy’s equation and the advection–diffusion equation. Additionally, various scenarios involving parametric models, unknown or missing information, and different levels of noisy data were considered. The complexity of the resulting PDE was found to be directly proportional to the inputted information. Furthermore, a polynomial regression method was employed to address the noisy interruption, yielding satisfactory results for noise levels of up to approximately 45%. This innovative approach significantly contributes to PDEs identification under varying conditions, ensuring a more physically grounded outcome.