This paper discusses a hybrid optimization method for FWI, whereby the l-BFGS and truncated Newton methods are used for FWI iterations. The latter method utilizes a matrix-free, finite difference-based formulation to compute Hessian-vector products, yielding model updates that are preconditioned by the inverse Hessian. The hybrid implementation achieves model resolution comparable to using a truncated Newton method for all iterations, but at significantly lower cost. Synthetic results demonstrate that the hybrid method delivers similar results as those expected with the truncated Newton method, and improved model reconstruction over gradient-based methods. Application to 3D North Sea data demonstrates uplift in imaging of near-surface velocity channels.