A computer-implemented method for avoiding drilling and production hazards in a subsurface reservoir, the method comprising: a. receiving, at one or more processors, 4D seismic data that has differences between the recorded seismic data at time T 1 and the seismic data recorded at time T 2 and at least one reservoir model, the differences caused by differences in subsurface properties which affect P-wave velocity, Vp, and S-wave velocity, Vs, where time T 2 is after time T 1; b. inverting, via the one or more processors, the 4D seismic data using elastic full waveform inversion to get stiffness changes, comprising: iterating, using the elastic full waveform inversion, the P-wave velocities and the S-wave velocities to obtain an inverted Vp perturbation model and an inverted Vs perturbation model; and obtaining the stiffness changes based on the inverted Vp perturbation model and the inverted Vs perturbation model; c. converting, via the one or more processors, the stiffness changes to strain changes using 3 rd -order elasticity relations; d. inverting, via the one or more processors, the strain changes to time-lapse stresses and reservoir pore pressures in a subsurface reservoir using elastostatic Green functions connecting reservoir pore pressure changes to induced strain; e. using the reservoir pore pressure and the time-lapse stress changes to refine the at least one reservoir model to generate at least one refined reservoir model; f. performing, via the one or more processors, seismic imaging using the at least one refined reservoir model to generate a reservoir image; and g. identifying or indicating potential drilling and production hazards based on the reservoir image. a. receiving, at one or more processors, 4D seismic data that has differences between the recorded seismic data at time T 1 and the seismic data recorded at time T 2 and at least one reservoir model, the differences caused by differences in subsurface properties which affect P-wave velocity, Vp, and S-wave velocity, Vs, where time T 2 is after time T 1; b. inverting, via the one or more processors, the 4D seismic data using elastic full waveform inversion to get stiffness changes, comprising: iterating, using the elastic full waveform inversion, the P-wave velocities and the S-wave velocities to obtain an inverted Vp perturbation model and an inverted Vs perturbation model; and obtaining the stiffness changes based on the inverted Vp perturbation model and the inverted Vs perturbation model; iterating, using the elastic full waveform inversion, the P-wave velocities and the S-wave velocities to obtain an inverted Vp perturbation model and an inverted Vs perturbation model; and obtaining the stiffness changes based on the inverted Vp perturbation model and the inverted Vs perturbation model; c. converting, via the one or more processors, the stiffness changes to strain changes using 3 rd -order elasticity relations; d. inverting, via the one or more processors, the strain changes to time-lapse stresses and reservoir pore pressures in a subsurface reservoir using elastostatic Green functions connecting reservoir pore pressure changes to induced strain; e. using the reservoir pore pressure and the time-lapse stress changes to refine the at least one reservoir model to generate at least one refined reservoir model; f. performing, via the one or more processors, seismic imaging using the at least one refined reservoir model to generate a reservoir image; and g. identifying or indicating potential drilling and production hazards based on the reservoir image.