A two-dimensional Galerkin finite element model for water flow in variably saturated soil is presented. A fourth-order Runge-Kutta time integration method is employed which allows use of time steps at least 2 times greater than for a traditional finite difference approximation of time derivatives. For short total simulation times computer execution costs for the Runge-Kutta method are greater than for the finite difference approximation due to the start up cost of the Runge-Kutta method, but for longer simulation times the Runge-Kutta method requires considerably less computational effort even when automatic time-step adjustment is used with the finite difference procedure. A comparison of the method of influence coefficients and 2 × 2 Gaussian integration to compute element matrices indicates that the influence coefficient method reduces total execution time to 60% of that required for numerical quadrature. Computed pressure heads using the influence coefficient method and numerical integration are found to be in close agreement with each other even under conditions of highly non-linear soil properties in a heterogeneous domain. Fluxes computed by the two methods are also generally in close agreement except under extremely non-linear conditions when some deviations were observed at short simulation times.