Simulation of variably saturated soil water flow requires the use of pressure head, or soil moisture, or a switching between the two, as the primary variable for solving Richards’ equation. Under unfavorable conditions, such as heterogeneity, rapidly changing atmospheric boundary, or sudden infiltration into dry soils, the traditional non-switching method suffers from numerical difficulties. Solving this problem with a primary variable switching method is less preferred due to the mathematical complexity. While the Picard method is more popular for solving the non-switching models due to its simplicity and stability, two different forms of Richards’ equation are combined into one numerical scheme for switching under specific hydraulic conditions. The method is successfully implemented in a one-dimensional model solved by a Picard iteration scheme. A threshold saturation based on the soil moisture retention relation is used for switching between either form of the Richards’ equation. The method developed here is applicable for simulating variably saturated subsurface flow in heterogeneous soils. Compared with traditional methods, the proposed model conserves mass well and is numerically more stable and efficient.