Wavelet estimation method, particularly advantageous for full wavefield inversion (“FWI”) of seismic data, that makes use of both the primary and multiple reflections in the data. The inventive method uses an FWI algorithm to generate a subsurface model from primary reflections (101) in a shallow layer before first arrival of multiple reflections (101). The model is then used to simulate multiples (102). The wavelet is subsequently modified (104) such that the simulated multiples closely match the true recorded multiples (103). The simulated multiples may then be subtracted from the measured data (105) thereby creating a deeper top layer of data substantially free of multiples, and the method may then be repeated to extend the subsurface model to a greater depth (106).