We consider the least-squares mixed finite element formulation of nonlinear elliptic boundary value problems using a combination of standard conforming finite elements with Raviart-Thomas spaces for the flux. For the solution of the arising nonlinear algebraic least-squares problems suitable variants of the Levenberg-Marquardt method are proposed. In particular, for elliptic problems arising from an implicit time discretization of a variably saturated subsurface flow model, we study an appropriate norm for the definition of the trust-region. We study an inexact version of this approach where the linear least-squares problems in each step are approximately solved by an adaptive multilevel method. For realistic variably saturated subsurface flow problems, the results of computational experiments are presented.