We compared three methods for discretizing the storage term of the Richards equation: the traditional chain rule expansion approach (Method A), a mass conservative approach based on a mixed form of the Richards equation (Method B), and another mass conservative approach using a chord-slope approximation for the specific moisture capacity (Method C). The results of three test problems indicated that Method A could not achieve a perfect global mass balance even if the iteration number in the Picard iteration scheme was large enough to bring a perfect solution convergence. Both Methods B and C successfully produced perfect mass balances. Method C produced step-by-step decreases in the global mass balance error as the Picard iteration level increased, which corresponded well with step-by-step decreases in the solution convergence error. Method B produced more accurate mass balances than Methods A and C for every Picard iteration level; when this method was used, the global mass balance error became negligible before the solution converged. Analytical evaluation of Method B revealed the mechanisms for removing the mass balance error. As the difference between the matric pressure head, ψ, lessens between the previous and current Picard iteration levels, the water retention curve becomes more linear in the region bounded by the two ψ values. As a result, the difference in the water content between two consecutive Picard iteration levels is accurately approximated by using the difference in the two ψ values, which results in a remarkable reduction in the mass balance error, allowing Method B to produce better results than Method C.