A method, comprising: performing, with a computer, a first full wavefield inversion process on input seismic data that includes free surface multiples, wherein the first full wavefield inversion process is performed with a free-surface boundary condition imposed on a top surface of an initial subsurface physical property model, and the first full wavefield inversion process generates a final subsurface physical property model; predicting, with the computer, subsurface multiples with the final subsurface physical property model; wherein the method further includes (a) removing, with the computer, the predicted subsurface multiples from the input seismic data and preparing multiple-free seismic data, and performing, after the removing, a second full wavefield inversion process on the input seismic data with the predicted subsurface multiples removed therefrom, wherein the second full wavefield inversion process is performed with an absorbing boundary condition imposed on a top surface of an initial subsurface physical property model, and the second full wavefield inversion process generates a multiple-free final subsurface physical property model, or (b) performing, with the computer, a second full wavefield inversion process on the input seismic data, wherein the second wavefield inversion process uses an objective function that only simulates primary reflections, the objective function being based on the predicted subsurface multiples, and the second full wavefield inversion process generates a multiple-free final subsurface physical property model; and using the multiple-free final subsurface physical property model as an input to an imaging or velocity model building algorithm, or in interpreting a subsurface region for hydrocarbon exploration or production. performing, with a computer, a first full wavefield inversion process on input seismic data that includes free surface multiples, wherein the first full wavefield inversion process is performed with a free-surface boundary condition imposed on a top surface of an initial subsurface physical property model, and the first full wavefield inversion process generates a final subsurface physical property model; predicting, with the computer, subsurface multiples with the final subsurface physical property model; wherein the method further includes (a) removing, with the computer, the predicted subsurface multiples from the input seismic data and preparing multiple-free seismic data, and performing, after the removing, a second full wavefield inversion process on the input seismic data with the predicted subsurface multiples removed therefrom, wherein the second full wavefield inversion process is performed with an absorbing boundary condition imposed on a top surface of an initial subsurface physical property model, and the second full wavefield inversion process generates a multiple-free final subsurface physical property model, or (b) performing, with the computer, a second full wavefield inversion process on the input seismic data, wherein the second wavefield inversion process uses an objective function that only simulates primary reflections, the objective function being based on the predicted subsurface multiples, and the second full wavefield inversion process generates a multiple-free final subsurface physical property model; and (a) removing, with the computer, the predicted subsurface multiples from the input seismic data and preparing multiple-free seismic data, and performing, after the removing, a second full wavefield inversion process on the input seismic data with the predicted subsurface multiples removed therefrom, wherein the second full wavefield inversion process is performed with an absorbing boundary condition imposed on a top surface of an initial subsurface physical property model, and the second full wavefield inversion process generates a multiple-free final subsurface physical property model, or (b) performing, with the computer, a second full wavefield inversion process on the input seismic data, wherein the second wavefield inversion process uses an objective function that only simulates primary reflections, the objective function being based on the predicted subsurface multiples, and the second full wavefield inversion process generates a multiple-free final subsurface physical property model; and using the multiple-free final subsurface physical property model as an input to an imaging or velocity model building algorithm, or in interpreting a subsurface region for hydrocarbon exploration or production.