A computer-implemented method for ensuring stability of iterative inversion of seismic data to infer a model of at least one physical property of a subsurface region, wherein a model update is computed, using a programmed computer, for a next iteration by optimizing an objective function measuring misfit between the seismic data and model-simulated seismic data, wherein said method comprises: determining when a model update will cause an unstable simulation, and in response to such a determination, using a Projection Onto Convex Sets to find a nearest stable model, wherein said method is characterized in that it further comprises: determining a matrix of elastic constants of the subsurface region, said matrix corresponding to said model of the at least one physical property of the subsurface region; determining whether an updated model will cause an unstable simulation based on whether the matrix of elastic constants corresponding to the updated model is or is not a positive semi-definite matrix; in response to a determination of instability, defining a stability projection operator that will convert the matrix corresponding to the updated model to a positive semi-definite matrix, said stability projection operator being a Projection Onto Convex Sets operator; applying the stability projection operator to the matrix corresponding to the updated model, generating a projected matrix, and adjusting the updated model to correspond to the projected matrix; and using the adjusted model to simulate synthetic seismic data for a next iteration. determining when a model update will cause an unstable simulation, and in response to such a determination, using a Projection Onto Convex Sets to find a nearest stable model, wherein said method is characterized in that it further comprises: determining a matrix of elastic constants of the subsurface region, said matrix corresponding to said model of the at least one physical property of the subsurface region; determining whether an updated model will cause an unstable simulation based on whether the matrix of elastic constants corresponding to the updated model is or is not a positive semi-definite matrix; in response to a determination of instability, defining a stability projection operator that will convert the matrix corresponding to the updated model to a positive semi-definite matrix, said stability projection operator being a Projection Onto Convex Sets operator; applying the stability projection operator to the matrix corresponding to the updated model, generating a projected matrix, and adjusting the updated model to correspond to the projected matrix; and using the adjusted model to simulate synthetic seismic data for a next iteration. determining a matrix of elastic constants of the subsurface region, said matrix corresponding to said model of the at least one physical property of the subsurface region; determining whether an updated model will cause an unstable simulation based on whether the matrix of elastic constants corresponding to the updated model is or is not a positive semi-definite matrix; in response to a determination of instability, defining a stability projection operator that will convert the matrix corresponding to the updated model to a positive semi-definite matrix, said stability projection operator being a Projection Onto Convex Sets operator; applying the stability projection operator to the matrix corresponding to the updated model, generating a projected matrix, and adjusting the updated model to correspond to the projected matrix; and using the adjusted model to simulate synthetic seismic data for a next iteration.