Inverse problems involve estimating parameters from physical measurements; for example, hydraulic tomography is a method to estimate hydraulic parameters related to the subsurface from pressure measurements obtained by a series of pumping tests. For many inverse problems, the solution of a shifted system of linear equations is the major computational bottleneck. By exploiting the shift-invariance property of Krylov subspaces, only a single Krylov basis is computed and the solution for multiple shifts can be obtained at a cost that is similar to the cost of solving a single system. We have developed flexible Krylov solvers for shifted systems and proposed fast methods for solving large-scale inverse problems based on these Krylov solvers. We demonstrate the performance gains (up to 20x speedup) on data sets taken from hydraulic tomography.