Three implicit numerical procedures for solving Richards' equation were compared. In some tests non-linear source terms were included. The main focus is on two mass-conserving procedures where the highly non-linear nature of Richards' equation is overcome by substantially different iterative techniques. In the first procedure (P1), which is well described in the literature, the time step is set back when difficulties in obtaining a convergent solution occur and any source terms are included explicitly. The second procedure (P2), which is presented in this paper, uses a fixed time step and convergence is obtained through a gradually increasing underrelaxation technique. Source terms are included implicitly, and in addition, a linear dependency between source terms and the water potential is included directly in the solution. The third procedure (P3) is based on a standard implicit discretization of Richards' equation, which does not facilitate mass conservation. All tests showed that mass conservation can be obtained without the added cost of more complicated programming or a longer calculation time. In challenging tests where fluctuating boundary conditions were specified as known fluxes, the total number of iterations needed in PI to simulate the tests was many orders of magnitude higher than that for P2. In tests including source terms, the more advanced treatment of these in P2 resulted in convergence rates that were significantly higher that for the explicit coupling in P1. The gradually increasing underrelaxation technique and the more advanced treatment of source terms in P2 have proven to be effective approaches, particularly in long-term simulations of water movements in soils where calculation time can be prohibitive. As a secondary benefit, the implementation of P2 is less complicated than P1. (C) 1999 Elsevier Science Ltd. All rights reserved.