A method and system for generating images of a subsurface region of interest. In general, one embodiment of the present invention includes establishing boundary conditions utilizing seismic data and initial conditions which include excitation from source locations in an earth model. Source wavefields are then propagated forward through the earth model to a maximum time, and saved at a plurality of checkpoints sparsely in time and also corresponding boundary values of the source wavefields at each time step are saved. Source wavefields are also propagated backward through the earth model from the maximum time utilizing the plurality of checkpoints when available and the saved boundary values at each time step. Receiver wavefields are propagated backward concurrently through the earth model from the maximum time.; Imaging conditions are applied at selected time steps to both the backward propagated source wavefields and receiver wavefields and those wavefields are utilized to generate images related to the subsurface region.