A robust, efficient numerical solution of variably saturated flows in three dimensions will be presented. The solution algorithm was developed by applying the mixed Lagrangian-Eulerian (LE) and Finite Element (FE) approach to Richard's equation. The Lagrangian-Eulerian approach with accurate and robust particle tracking algorithms has been considered a better resolution for sharp front problems. However, the implementation of the algorithm to the advective form of the Richards equation has encountered difficulties resulting from the flux boundary conditions. Therefore, it was decided to develop a mixed LE and FE approach to circumvent the problem. In this approach, LE method was applied to interior nodes and the FE method to incoming-flux-boundary nodes. The goal was to develop numerical tools that can efficiently analyze the flow behaviors in saturated/unsaturated porous media. Application of the hybrid approach to a test problem shows that it is more robust and efficient over the conventional finite element methods.