Unsaturated flow problems in porous media often described by Richards’ equation are of great importance in many engineering applications. In this contribution, we propose a new numerical flow approach based on isogeometric analysis (IGA) for modeling the unsaturated flow problems. The non-uniform rational B-spline (NURBS) basis is utilized for spatial discretization whereas the stable implicit backward Euler method for time discretization. The nonlinear Richards’ equation is iteratively solved with the aid of the Newton–Raphson scheme. Owing to some desirable features of an efficient numerical flow approach, major advantages of the present formulation involve: (a) numerical oscillation at the wetting front can be avoided or facilitated, simply by using either an h-refinement or a lumped mass matrix technique; (b) higher-order exactness can be obtained due to the nature of the IGA features; (c) the approach is straightforward to implement and it does not need any transformation, e.g., Kirchhoff transformation or filter algorithm; and (d) in contrast to the Picard iteration scheme, which forms linear convergences, the proposed approach can however yield quadratic convergences by using the Newton–Raphson method for solving resultant nonlinear equations. Numerical model validation is analyzed by solving a three-dimensional unsaturated flow problem in soil, and its derived results are verified against analytical solutions. Numerical applications are then studied by considering three extensive examples with simple and complex configurations to further show the accuracy and applicability of the present IGA.