The modelling of unsaturated-saturated flow in soil leads to a first-order system of time-dependent nonlinear partial differential equations for the flux and the hydraulic potential. The least-squares finite element approach is applied to the system resulting from fully implicit time discretization using lowest-order Raviart-Thomas elements for the flux and standard conforming $P_1$ elements for the potential. We solve the resulting nonlinear algebraic least-squares problems by a Gauss-Newton method involving the solution of a sequence of linear least-squares problems at each time step. The properties of the least-squares functional as a local error estimator and the effect of the resulting locally refined triangulations on the accuracy and on the performance of the Gauss-Newton method are studied. Finally, we show numerical results obtained for some transient two-dimensional unsaturated-saturated subsurface flow problems from hydrology applications.