A computer implemented stochastic inversion method for estimating model parameters of an earth model of a subsurface geological volume of interest, the method comprising: a) acquiring at least one geophysical data set that samples a portion of the subsurface geological volume of interest, wherein each geophysical data set defines an acquisition geometry of the subsurface geological volume of interest; b) generating a specified number of boundary-based multi-dimensional models of the subsurface geological volume of interest, said models being defined by model parameters; c) generating forward model responses of the models for each specified acquisition geometry; d) generating a likelihood value of the forward model responses matching the geophysical data set for each specified acquisition geometry; e) saving the model parameters as one element of a Markov Chain for each model; f) testing for convergence of the Markov Chains; g) updating the values of the model parameters for each model and repeating c) to f) in series or in parallel, until convergence is reached; h) deriving probability density functions for each model parameter of the models which form the converged Markov Chains; and i) calculating the variances, means, modes, and medians from the probability density functions of each model parameter for each model to generate estimates of model parameter variances and model parameters for the earth models of the subsurface geological volume of interest which are utilized to determine characteristics of the subsurface geological volume of interest, wherein Markov Chain Monte Carlo sampling is used to generate, through iterations of the updating and repeating which each comprise saving the model parameters, a sequence of saved model parameters for each model that form the Markov Chain concerned. a) acquiring at least one geophysical data set that samples a portion of the subsurface geological volume of interest, wherein each geophysical data set defines an acquisition geometry of the subsurface geological volume of interest; b) generating a specified number of boundary-based multi-dimensional models of the subsurface geological volume of interest, said models being defined by model parameters; c) generating forward model responses of the models for each specified acquisition geometry; d) generating a likelihood value of the forward model responses matching the geophysical data set for each specified acquisition geometry; e) saving the model parameters as one element of a Markov Chain for each model; f) testing for convergence of the Markov Chains; g) updating the values of the model parameters for each model and repeating c) to f) in series or in parallel, until convergence is reached; h) deriving probability density functions for each model parameter of the models which form the converged Markov Chains; and i) calculating the variances, means, modes, and medians from the probability density functions of each model parameter for each model to generate estimates of model parameter variances and model parameters for the earth models of the subsurface geological volume of interest which are utilized to determine characteristics of the subsurface geological volume of interest, wherein Markov Chain Monte Carlo sampling is used to generate, through iterations of the updating and repeating which each comprise saving the model parameters, a sequence of saved model parameters for each model that form the Markov Chain concerned.