Uncertainty quantification of geological model parameters in 3D gravity inversion by Hessian informed Markov chain Monte Carlo

Geological modeling has been widely adopted to investigate underground geometries. How- 12 ever, modeling processes inevitably have uncertainties due to scarcity of data, measurement 13 errors and simplification of modeling method. Recent developments in geomodeling methods 14 have introduced a Bayesian framework to constrain the model uncertainties by consider- 15 ing additional geophysical data into the modeling procedure. Markov chain Monte Carlo 16 (MCMC) methods are normally used as tools to solve the Bayesian inference problem. To 17 achieve a more efficient posterior exploration, advances in MCMC methods utilize deriva- 18 tive information. Hence, we introduce an approach to efficiently evaluate second-order 19 derivatives in geological modeling and introduce a Hessian-informed MCMC method, the 20 generalized preconditioned Crank-Nicolson (gpCN), as a tool to solve the 3D model-based 21 gravity Bayesian inversion problem. The result is compared with two other widely applied 22 MCMC methods, random walk Metropolis-Hasting and Hamiltonian Monte Carlo, on a 23 synthetic three-layer geological model. Our experiment demonstrates that superior perfor- 24 mance is achieved by the gpCN, which has the potential to be generalized to more complex 25 models. 26


INTRODUCTION
In many geoscience applications, inversion methods are used to estimate subsurface prop-27 erties (e.g., structure, density and porosity) from observed geophysical data. Conventional 28 geophysical inversion aims to find the best-fit parameter sets that minimize the error be-29 tween observed geophysical data and simulation results. However, in practical cases, ob-

PRELIMINARIES
In this section, we briefly introduce the forward model, including the construction of the 131 structural geological model and the gravity forward simulations, and formulate the Bayesian 132 inverse problem. We review two commonly used MCMC algorithms namely the RMH 133 algorithm and HMC algorithm, which were implemented in our study as a comparison.
134 Implicit modeling 135 To generate the 3D geological models, we applied a implicit surface representation method, Solving the co-kriging system based on the interface variance and covariance information 156 and the above constraints provides a continuous scalar field. A 2D scalar field example Gravity method, to date still one of the most widely used geophysical techniques in ex-162 ploration (Nabighian et al., 2005) can be used to integrate geophysical information into 163 geomdeling procedures. One of the most popular methods to include geophysical data is 164 through inversion (Tarantola, 2005 where x,y, and z are the Cartesian components of the targeting prism, ρ is density and 172 r stands for the Euclidian distance from the center of the cell to the receiver points. In this 173 work we use virtual receivers on the ground surface, which is usually the top of the model.

174
By summing up the attractions of all the cells, the gravity at that receiver point can be 175 simulated.

176
The resolution of the model, which reflects the size of the cells, has an impact on the simulation. Low-resolution cells introduce more uncertainties into the model. Additionally, receivers close to the boundary of the model will have boundary effects. Therefore, it is necessary to further extend the model to minimize the boundary effects. A model with a small cell size as well as a large model extent would certainly be ideal; however, our computational resources are limited. The attraction force decays exponentially with respect to distance r from the gravitation equation: Hence, instead of using a regular grid (which means all cells have the same size), we apply 177 a irregular grid (termed as centered grid in this paper) (see also (Güdük et al., 2021,p.5)).

178
To simplify the calculation, we make this grid isotropic around a center position where the 179 gravity response is calculated (illustrated in Figure 2). Both the height and widths of the  Bayesian inference starts from Bayes' theorem, where p (d obs | m) is known as the likelihood of the data d obs given model parameters and C is the covariance matrix of the proposal distribution.

212
The random-walk method has limitations in cases of high dimensionality and high correlations. Gradient information is therefore employed to assist in the posterior exploration.
The most popular gradient-based method, the HMC algorithm, is designed to draw independent samples, and therefore efficiently explore the state space (Duane et al., 1987; Neal, 1993; Chen et al., 2014; Betancourt, 2017). HMC introduces auxiliary momentum variables r, and therefore we can write the Hamiltonian of a particle as: where K and V are kinetic and potential energies respectively and are defined as: Here, the mass matrix M is often set as the identity matrix (Chen et al., 2014)

213
Hamiltonian dynamics can be simulated by following Hamilton's equations with an 214 artificially introduced time variable t: The HMC algotithm requires that Hamilton's equations are integrated using a symplec-216 tic integrator. In practice, the "leapfrog" integrator is often used. The algorithm can be Simulate discretization of Hamiltonian dynamics in Eq. 1: Propose (m cand , r cand ) = (m n , r n ) RMH is simple to implement and gradient-free, while HMC uses gradient information 221 to obtain low autocorrelated samples; however, both methods are popular and widely ap-222 plied. In this study, we will compare the efficiency of the generalized preconditioned Crank-

METHODS
In this work, we conducted an end-to-end procedure from the model construction to the 225 geophysical simulation and finally to Bayesian inference. We first generated a geological  of Nadam is given as follows: where α is the step size.b t andv t are the bias-corrected first and second momentum  Using methods such as finite difference (FD) is not only computationally costly but also 273 can suffer from numerical inaccuracy. Therefore, we adopted the automatic differentiation 274 (AD) technique, which is widely applied in the field of artificial intelligence and is critical 275 to the success of training neural networks. Here we briefly introduce how AD works and how gradient and higher-order derivatives can be evaluated efficiently using AD.

277
Consider an arbitrary function y = F (x), where F is the mapping function maps the to the cost function y = y 1 , y 2 , ..., y j , (y ∈ R D 1 ).

279
The first-order partial derivative, also known as the Jacobian matrix, is given as follows: Jacobian matrix can be represented by iteratively applying the chain rule: approach is summarized in the following.

289
To evaluate the gradient in the forward mode AD, a tangent vector v is selected at the 290 evaluation point x, For example, ∂(F ) i ∂(x 1 ) can be evaluated by a tangent vector v = (1, 0, . . . , 0), v ∈ R i . Thus, 292 the forward mode AD provides the directional derivative.

293
In contrast, backward mode AD evaluation is based on a fixed dependent variable, and the quantity of interest is the adjoint

RESULTS
When the available data are limited, uncertainties are inevitable in 3D geological models. 313 We can constrain the uncertainties by additional geophysical data. To configure such a 314 problem in a Bayesian inference framework as described above, both the prior information 315 and the likelihood must be expressed in terms of a probability distribution. As an example 316 we constructed a 3D geological model using the implicit method described above. The   (Table 1), while preserving the computa-381 tional efficiency, even though we only used three leapfrog steps in the experiments for HMC.

382
Although we demonstrate only a simple example, we can expect a better performance of 383 the gpCN in a higher dimension and more correlated case.

384
The main reason for the superior performance of gpCN is the highly correlated posterior 385 space. If we use the sample of HMC as a benchmark and plot the adopted proposal distri-386 bution used in the gpCN (Figure 9), we can see that the proposal distribution captured the    Finally, we represent the uncertainties following the information entropy method intro-389 duced by Wellmann and Regenauer-Lieb (2012). In Figure 10, we compared the uncer-   can benefit greatly from Hessian information, even for high-dimensional parameters.

CONCLUSION
In summary, this study extended the previous development of stochastic geological modeling 459 methods . We used the automatic differential technique implemented in TensorFlow to 460 allow second-order derivative information to be efficiently evaluated in geological models. 461 We applied the recently developed Hessian-informed MCMC, the generalized preconditioned