Newton-Krylov-Schwarz methods solve nonlinear equations by using Newton's method with a Schwarz domain decomposition preconditioned Krylov method to approximate the Newton step. In this talk we will discuss the design and implementation of Newton-Krylov-Schwarz solvers in the context of the implicit temporal integration on an unstructured three-dimensional spatial mesh of Richards' equation for flow in the unsaturated zone. The issues include nonsmooth nonlinearities, construction and efficient implementation of the coarse-mesh problem, and temporal integration. The second part of the talk, given by E. W. Jenkins, will discuss the convergence theory for the method.