Physics-informed neural networks (PINNs) integrate physical principles into machine learning, finding wide applications in various scientific and engineering fields. However, solving nonlinear hyperbolic partial differential equations (PDEs) with PINNs presents challenges due to inherent discontinuities in the solutions. This is particularly true for the Buckley-Leverett (B-L) equation, a key model for multiphase fluid flow in porous media. In this paper, we demonstrate that PINNs, in conjunction with Welge's construction, can achieve superior precision in handling the B-L equations in different scenarios including one shock and one rarefaction wave, two shocks connected by a rarefaction wave traveling in the same direction, and two shocks connected by a rarefaction wave traveling in opposite directions. Our approach accounts for variations in fluid mobility, fluid solubility, and gravity effects, with applications in modeling 1D water flooding, polymer flooding, gravitational flow, and CO2 injection into saline aquifers. Additionally, we applied PINNs to inverse problems to estimate multiple PDE parameters from observed data, demonstrating robustness under conditions of slight scarcity and up to 5% impurity of labeled data as well as shortages in collocation data.