The inverse problem in groundwater modeling is often made numerically tractable and computationally practical by estimating only a small fraction of the many unknown system parameters. However, this parsimonious approach restricts the solution of the inverse problem to a pre-determined subspace of the true parameter space. To reflect detailed local variations in hydraulic conductivity or recharge, it may be desirable to estimate a very large number of parameters during calibration, which requires an inversion technique that can accommodate highly parameterized models. The least-squares QR (LSQR) decomposition is an iterative solution method that can solve for many hundreds or thousands of parameters. LSQR has been used successfully in seismic tomography inversion problems. As an iterative method, LSQR can solve sparse and dense inverse problems of the form Ax=b using significantly less computer storage than direct solution methods. We test the applicability of the LSQR method for solving the inverse problem for groundwater flow using a synthetic model and compare results with those obtained using the more commonly employed method, the singular value decomposition (SVD). Parameter sensitivities are calculated using forward differences and the adjoint-state method.