Addressing the issue of model error in stochastic geophysical inversion

Solving inverse problems in geophysics and hydrology requires a forward model linking subsurface physical parameters to measured data, which is typically assumed to be perfectly known in the inversion procedure. However, to make the solution of the inverse problem computationally tractable using stochastic methods in order to quantify uncertainties, for example through Markov-chain-Monte-Carlo (MCMC) posterior sampling, fast approximations of the forward model are commonly considered. Use of such “proxy” forward solvers introduces model error into the inverse problem, which can strongly bias posterior statistics, yield highly over-confident estimates of the wrong parameters, and hamper data integration efforts. Here, we present a new methodology for addressing the issue of model error in Bayesian solutions to hydrogeophysical inverse problems that is geared towards the common case where these errors cannot be effectively described (i) globally through some parametric statistical distribution; or (ii) locally based on interpolation between a small number of computed stochastic model-error realizations. Rather than attempting to build a global or local error model, we instead focus on identifying the model-error component of the residual through an orthogonal-projection-based approach. To this end, pairs of approximate and detailed forward solver runs are stored in a dictionary that grows at a specified rate during the stochastic inversion procedure. At each iteration, a local model-error basis is constructed for the current test set of model parameters using the K-nearest neighbouring entries in the dictionary. This basis is then used to identify and separate the model-error component of the residual from the parameter-and measurement-error components before computing the likelihood of the proposed set of model parameters. As an example, we apply our technique to crosshole GPR traveltime tomography for three different subsurface parameterizations of varying complexity: 5 horizontal layers; a 2D medium described by 20 Karhunen-Loève coefficients; and a 20x40 pixel-based parameterization. We generate synthetic data using the eikonal equation (our assumed “full-physics” forward solver), whereas a straight-ray solver (our proxy) is assumed in the inversion procedure. In each case, our developed model-error approach allows us remove posterior bias and obtain a more honest characterization of posterior uncertainty.