1. A method for determining a discrete physical properties model of a subsurface region, referred to herein as the model or the subsurface model, by iteratively inverting measured geophysical data acquired from the subsurface region, comprising: approximating a Hessian matrix of an objective function, then times a vector, called Hessian times vector, using a computer, with a single forward-wave simulated propagation and a single computation of gradient of the objective function, in a modified subsurface model, thereby requiring only three forward-wave or reverse-wave propagations; wherein said approximation is based on an approximate equation for a Born scattered pressure field, {tilde over (p)} d,F , where the Hessian times vector is approximated by a gradient computation using {tilde over (p)} d,F as an artificial residual; then computing a direction in model parameter space for an update to a current model by multiplying inverse of the Hessian matrix times a gradient of the objective function, wherein the inverse of the Hessian matrix is computed iteratively using a conjugate gradient method in which said approximate Hessian times vector is used to evaluate the Hessian matrix times a model perturbation vector; performing a line search to determine magnitude of the model update using said computed direction; and adding the model update to the current model to form an updated model, and using the updated model for geophysical prospecting.