Variably saturated subsurface flow is modeled with a scalar, highly nonlinear equation. After spatial and temporal discretization, the nonlinearities in this equation are often solved by application of Newton's method within each time step. For large-scale problems, the Newton method is combined with a preconditioned Krylov method for solution of the linear systems within each nonlinear iteration. In this talk, we discuss both software and mathematical issues related to transitioning the variably saturated flow module within the PNNL subsurface simulation code, STOMP, from use of a block ILU preconditioner to use of a multigrid preconditioner. Results related to tuning the multigrid preconditioner will be shown as well as weak scaling studies.