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. <EOS>