Physics-informed neural networks (PINN) have recently become attractive to solve partial differential equations (PDEs) that describe physics laws. By including PDE-based loss functions, physics laws such as mass balance are enforced softly. In this paper, we examine how well mass balance constraints are satisfied when PINNs are used to solve the resulting PDEs. We investigate PINN's ability to solve the 1D saturated groundwater flow equations for homogeneous and heterogeneous media and evaluate the local and global mass balance errors. We compare the obtained PINN's numerical solution and associated mass balance errors against a two-point finite volume numerical method and the corresponding analytical solution. We also evaluate the accuracy of PINN to solve the 1D saturated groundwater flow equation with and without incorporating hydraulic head as training data. Our results showed that one needs to provide hydraulic head data to obtain an accurate solution for the head for the heterogeneous case. That is, only after adding data at some collocation points the PINN methodology was able to give a solution close to the finite volume and analytical solutions. In addition, we demonstrate that even after adding hydraulic head data at some collocation points, the resulting local mass balance errors are still significant at the other collocation points compared to the finite volume approach. Tuning the PINN's hyperparameters such as the number of collocation points, epochs, and learning rate did not improve the solution accuracy or the mass balance errors compared to the finite volume solution. Mass balance errors may pose a significant challenge to PINN's utility for applications where satisfying physical and mathematical properties is crucial.