Two alternative methods for numerical modeling of variably saturated flow in layered media are developed. The first method is based on a space and time finite difference discretization of the mixed form of Richards' equation. The solution of the resulting equations follows the scheme introduced by Celia et al. The unique features in the proposed solution are treatments of (i) the saturated–unsaturated transition and (ii) the interface between material layers. In the second method, following work by Marinelli and Durnford, Richards' equation is discretized in time. At each time step a system of two first-order equations in space is obtained. These equations are solved with an ordinary differential equation (ODE) shooting method. The ODE solution uses a fixed space grid and first order Euler integration. Both the finite difference and shooting methods are used to model the drainage of an initially saturated one-dimensional three-layer soil column. Grid converged results indicate a close agreement between the finite difference and shooting method predictions.