The flow of water through variably saturated subsurface media is commonly modeled by Richards' equation, a nonlinear and possibly degenerate partial differential equation. Due to the nonlinearities, this equation is difficult to solve analytically and the literature reveals dozens of papers devoted to finding numerical solutions. However, the literature also reveals a lack of two important research topics. First, no a priori error analysis exists for one of the discretization schemes most often used in discretizing Richards' equation, cell-centered finite differences. The expanded mixed finite element method reduces to cell-centered finite differences for the case of the lowest-order discrete space and certain quadrature rules. Expanded mixed methods are useful because this simplification occurs even for the case of a full coefficient tensor. There has been no analysis of expanded mixed methods applied to Richards' equation. Second, no results from parallel computer codes have been published. With parallel computer technology, larger and more computationally intensive problems can be solved. However, in order to get good performance from these machines, programs must be designed specifically to take advantage of the parallelism. We present an analysis of the mixed finite element applied to Richards' equation accounting for the two types of degeneracies that can arise. We also consider and analyze a two-level method for handling some of the nonlinearities in the equation. Lastly, we present results from a parallel Richards' equation solve code that uses the expanded mixed method for discretization.