Variational data assimilation is a technique for combining measured data with dynamical models. It is a key component of Earth system state estimation and is commonly used in weather and ocean forecasting. The approach involves a large-scale generalized nonlinear least-squares problem. Solving the resulting sequence of sparse linear subproblems requires the use of sophisticated numerical linear algebra methods. In practical applications, the computational demands severely limit the number of iterations of a Krylov subspace solver that can be performed and so high-quality preconditioners are vital. In this paper, we introduce variational data assimilation from a numerical linear algebra perspective and review current solution techniques, with a focus on the challenges that arise in large-scale geophysical systems.