The purpose of this work is to test different computational algorithms for unsaturated flow for accuracy and robustness by comparing computed results in a finite element program with analytical solutions. Because real-world problems are complex, testing codes for accuracy is often difficult. This is particularly true for flow in the vadose zone where Richards' equation is highly nonlinear. Recently, however, Tracy (Tracy WRRJ 2006) [1] (Tracy JHYD 2007) [2] has derived analytical solutions for a box-shaped flow region that is initially dry until water is applied to the top of the region. Two-dimensional and three-dimensional versions of these solutions for both steady-state and transient flow are available to be used in the testing process. Numerical precision and nonlinear solver robustness were investigated for varying degrees of nonlinearity by varying the Gardner α parameter. As α was increased, three ways of modeling relative hydraulic conductivity inside individual finite elements and two versions of the nonlinear solver were tested using three different ways to measure the error. The results of these tests are given in this paper.