1. A method comprising: defining, with a computer, a virtual subsurface box having sources and receivers located on its edges or faces, wherein the virtual subsurface box represents a selected part of Earth's subsurface; estimating, with a computer, Green's functions by estimating direct arrival times from locations of the sources to locations of the receivers; updating, with a computer, the Green's functions by one of iterative optimization that minimizes error in source field reconstruction, a series expansion, or iterative wavefield extrapolation; and running a simulation, with a computer, wherein the simulation simulates seismic waves generated by the sources being reflected off of at least two subsurface reflectors, which causes simulated internal multiple reflections within the box from the optimized Green's functions, and subtracting the simulated internal multiples from seismic data collected by actual receivers, wherein the seismic data includes internal multiple reflections, and the subtracting reduces the internal multiple reflections in the seismic data; generating, with a computer, a subsurface image from the seismic data with the simulated internal multiples subtracted therefrom.