Changes in soil hydraulic properties with depth are common in field soils and need to be described in numerical simulations. Water content is undesirable to calculate Darcian flux since it is not continuous across boundaries between different soil layers. A general form of the one-dimensional water content–based Richards equation, which was first derived by Hills et al., is adopted here to simulate one-dimensional unsaturated flow in gradational soils by adding a correction term. For cell-centered grid, several algorithms are presented to find the correction term for the two nodes in the vicinity of the layer interface. For vertex-centered grid, the dyadic values of water content at the interface are represented by one sole water content value through the introduction of composite soil water curve. The proposed algorithms are implemented in an iterative model, and then the model is tested with several synthetic cases in gradational soil and layered soil. The water content–based Richards equation leads to superior numerical performance in terms of mass conservation, accuracy, and efficiency over the mixed form Richards equation, especially for the flow in relatively dry soil and with coarse grid. In layered soil, it is found that the algorithm based on vertex-centered grid is the most robust among all the methods discussed here.