{'fraction_non_alphanumeric': 0.04090739497052459, 'fraction_numerical': 0.014133273636295799, 'mean_word_length': 4.785356587808669, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Lithological Tomography with the Correlated Pseudo-Marginal Method

We consider lithological tomography in which the posterior distribution of (hydro)geological parameters of interest is inferred from geophysical data by treating the intermediate geophysical properties as latent variables. In such a latent variable model, one needs to estimate the intractable likelihood of the (hydro)geological parameters given the geophysical data. The pseudo-marginal method is an adaptation of the Metropolis-Hastings algorithm in which an unbiased approximation of this likelihood is obtained by Monte Carlo averaging over samples from, in this setting, the noisy petrophysical relationship linking (hydro)geological and geophysical properties. To make the method practical in data-rich geophysical settings with low noise levels, we demonstrate that the Monte Carlo sampling must rely on importance sampling distributions that well approximate the posterior distribution of petrophysical scatter around the sampled (hydro)geological parameter field. To achieve a suitable acceptance rate, we rely both on (1) the correlated pseudo-marginal method, which correlates the samples used in the proposed and current states of the Markov chain, and (2) a model proposal scheme that preserves the prior distribution. As a synthetic test example, we infer porosity fields using crosshole ground-penetrating radar (GPR) first-arrival travel times. We use a (50x50)-dimensional pixel-based parameterization of the multi-Gaussian porosity field with known statistical parameters, resulting in a parameter space of high dimension. We demonstrate that the correlated pseudo-marginal method with our proposed importance sampling and prior-preserving proposal scheme outperforms current state-of-the-art methods in both linear and non-linear settings by greatly enhancing the posterior exploration.


Introduction
Geophysical investigations are rarely performed with the sole aim of inferring distributed subsurface models of geophysical properties. Rather, the underlying motivation is often to gain knowledge and constraints on other properties (e.g., permeability, clay fraction or mineral composition) and state variables (e.g., water saturation, salinity, temperature) of interest. Geophysical inverse theory has traditionally focused on assessing the resolution and uncertainty of inferred geophysical properties (e.g., Parker 1994;Menke 2018;Tarantola 2005;Aster et al. 2018), while interpretation procedures in terms of properties or state variables of interest have received less attention. This is changing in hydrogeophysics (Binley et al. 2015), for instance, where it is now well-established that dedicated inversion approaches are needed when using geophysical data to gain knowledge about hydrogeological properties and state variables (e.g., Kowalsky et al. 2005). For example, when inferring hydraulic conductivity by observing geophysical observables sensitive to water content or salinity during a tracer test experiment (Linde & Doetsch 2016). However, these considerations have general validity and relevance for exploration and more fundamental geophysical studies. In a mantle context, for instance, one example concerns the inference of thermo-chemical constraints from seismological observations as reviewed by Zunino et al. (2016).
Multiple inversion frameworks have been proposed that combine hydrogeological and geophysical data in order to build predictive hydrogeological models (e.g., Ferré et al. 2009;Linde & Doetsch 2016). A critical aspect of such frameworks relates to how geophysical properties (sensed by geophysical data) are linked to hydrogeological target properties and variables of interest through petrophysical (rock physics) relationships. Brunetti & Linde (2017) distinguish between three sources of uncertainty related to petrophysical relationships: model uncertainty, parameter uncertainty and prediction uncertainty. While the first two refer to uncertainty in the choice of the appropriate petrophysical model and its parameter values, the latter is related to scatter and bias around the calibrated petrophysical model. In hydrogeophysical inversion studies targeting hydrogeological properties or state variables of interest, we note that the petrophysical relationship is often assumed to be perfect (deterministic) with known or unknown parameter values (e.g., Lochbühler et al. 2014;Kowalsky et al. 2005). However, ignoring petrophysical prediction uncertainty and its spatial correlation patterns results in bias, too narrow uncertainty bounds and overly variable hydrogeological parameter estimates (Brunetti & Linde 2017). Unfortunately, analytical solutions to such inverse problems are available only when considering linear forward models and petrophysical relationships under the assumption of Gaussian distributions (Tarantola 2005;Bosch 2004). Geophysical applications, however, often involve non-linear physics and non-linear petrophysical relationships (e.g., Mavko et al. 2009).
Inversion approaches that account for petrophysical prediction uncertainty are often based on a twostep procedure: geophysical properties are first estimated using deterministic gradient-based inversions and then converted into parameters of interest using uncertain petrophysical relationships (e.g., Chen et al. 2001;Mukerji et al. 2001;Gonzalez et al. 2008;Grana & Della Rossa 2010;Shahraeeni & Curtis 2011). The results of such a two-step approach can be misleading if neglecting the spatially-varying and typically much lower resolution of smoothness-constrained geophysical inversion models compared with the scale at which petrophysical relationships are developed (core or borehole logging scale) (Day-Lewis et al. 2005). Furthermore, with such an approach it is next to impossible to ensure that the geophysical inversion accounts for the prior constraints on the (hydro)geological target variable (Ferré et al. 2009) and physical constraints such as conservation of mass, continuity and momentum. Moreover, for a deterministic inversion setting, Bosch (2004) showed that with a non-linear petrophysical relation, the two-step approach is an inherent approximation (Bosch 2004).
As an alternative to the two-step approach, coupled inversions directly target hydrogeological properties by inversion of geophysical data (e.g., Hinnell et al. 2010;Kowalsky et al. 2005). They are often formulated within a Bayesian framework whereby one seeks to characterize the posterior probability density function (PDF) of hydrogeological parameters θ given geophysical data y. Since it is often impossible to sample directly from the posterior PDF p(θ|y) of interest, Markov chain Monte Carlo (MCMC) methods, such as the Metropolis-Hastings method (MH;Hastings 1970;Metropolis et al. 1953), are used. Since the intermediate variable, the geophysical property X, connecting observations and target variables is unobservable (latent), one speaks of a latent variable model. In this study, we consider a setup where the latent geophysical property is given by X = F(θ) + ε P , with θ → F(θ) representing the deterministic component of a petrophysical relationship and ε P the petrophysical prediction error. Assuming an integrable and centered petrophysical prediction error ε P , F(θ) stands for the expected value of the latent variable X. The geophysical data is given by Y = G(X) + ε O with x → G(x) denoting the geophysical forward solver and ε O describing the observational noise.
For a latent variable model as the one described above, the likelihood of observing the geophysical data given the proposed hydrogeological parameters, p(y|θ) = p(y, x|θ)dx, is often intractable. In the present context, this implies that the integral has an unknown or non-existing analytical form, which makes the direct implementation of the MH and related algorithms impossible. One way to circumvent this difficulty is to instead infer the joint posterior PDF (θ, x) → p(θ, x|y) of the hydrogeological and geophysical parameters from which p(θ|y) is readily obtained by marginalization. Lithological tomography as introduced by Bosch (1999) pioneered such an approach to estimate the joint posterior by combining geophysical data, geological prior knowledge and uncertain petrophysical relationships. Within lithological tomography, pairs of the target and latent variables are proposed using marginal sampling of θ and conditional sampling of X. Then, these pairs are accepted or rejected with p(y|θ, x), used in the acceptance ratio of the MH algorithm (where p(y|θ, x) = p(y|x) is valid for our latent variable model). In Bosch (1999), the conditional PDF p(x|θ) to sample X is given by a multivariate Gaussian distribution based on a suitable petrophysical relationship. In practice, this is achieved by adding brute force Monte Carlo realizations of the petrophysical prediction error ε P to the output of F(θ) at each iteration of the MCMC chain (i.e., Bosch et al. 2007). Linde et al. (2017) suggest that such an implementation is inefficient when considering large geophysical datasets with high signal-to-noise ratios and significant petrophysical uncertainty. The reason is that brute force Monte Carlo sampling of the petrophysical prediction error using p(x|θ) induces high variability in the values taken by the likelihood function p(y|θ, x), even for the same θ, which could lead to prohibitively low acceptance rates even in the limiting case when the MCMC model proposal scale for θ goes to zero. Brunetti & Linde (2017) proposed an alternative approach to sample from the joint posterior PDF p(θ, x|y). In their method referred to herein as full inversion, the petrophysical prediction error ε P is parameterized and treated as the other unknowns within the MH algorithm. That is, the MH proposal mechanism draws new realizations of both the target variable θ and the petrophysical prediction error ε P , which combined also lead to a realization of the latent variable X used to calculate the likelihood function p(y|θ, x). Brunetti & Linde (2017) presented a convincing performance of the full inversion approach with clear improvements in efficiency compared with the original formulation of lithological tomography by Bosch (1999). Nonetheless, the full inversion method suffers from high dimensionality, and the strong (posterior) correlation between ε P and θ makes standard MCMC inversions inefficient (e.g., Deligiannidis et al. 2018).
In this study, we evaluate an inversion method targeting directly the marginal posterior p(θ|y) by approximating the intractable likelihood p(y|θ) = p(y|θ, x)p(x|θ)dx. In the pseudo-marginal (PM) method introduced by Beaumont (2003) and studied by Andrieu & Roberts (2009), the true likelihood is replaced with a non-negative unbiased estimator resulting in a MH algorithm sampling the same target distribution as when using the true likelihood. In their work, Beaumont (2003) and Andrieu & Roberts (2009) use an unbiased likelihood estimator based on Monte Carlo averaging over samples of the latent variable. In our setting with the latent variable X = F(θ) + ε P , we note that the original lithological tomography approach of Bosch (1999) is closely related to the pseudo-marginal method. In the original lithological tomography method targeting the joint posterior PDF p(θ, x|y), the MCMC chains store the conditional draws of the latent variables together with the target variables, and the target posterior PDF p(θ|y) is obtained by marginalization. The PM method applied with one draw of the latent variable leads to equivalent results in terms of the marginal posterior PDF. In the PM method, the draws of the latent variable are not stored but only used to estimate the likelihood p(y|θ). Using only one sample of the latent variable in the PM method typically leads to impractically-low acceptance rates due to the high variability of the ratio of log-likelihood estimators. To achieve an efficient algorithm, the standard deviation of the log-likelihood estimator needs to be around 1.2-1.5 (Doucet et al. 2015), which is ensured by increasing the number of samples and applying importance sampling. schemes. In the context of state-space models, the number of Monte Carlo samples used in the likelihood estimator needs to increase linearly with the number of observations, which becomes impractical in data-rich applications (Deligiannidis et al. 2018). To obtain low-variance log-likelihood ratio approximations with a smaller number of Monte Carlo samples, Deligiannidis et al. (2018) introduced the correlated pseudo-marginal (CPM) method by which the draws of latent variables used in the denominator and numerator in the likelihood ratio are correlated. Both the PM and CPM methods are general in that they allow for non-linear and non-Gaussian assumptions, but their implementation and applicability in data-rich high-dimensional geophysical settings remain untested.
Inferring hundreds or thousands of parameters with a MH algorithm is challenging as the number of iterations needed for convergence grows with the number of target parameters (e.g., Robert et al. 2018). To ensure adequate performance in such settings, it is crucial to equip the algorithm with a well-working proposal scheme. In the context of Gaussian random fields with high dimension, Cotter et al. (2013) demonstrated that standard random walk MCMC algorithms leads to strong dependence on the discretization of the target field and highly inefficient algorithms. Their proposed solution lies in preserving the prior PDF within the proposal scheme such that the acceptance probability of model proposals only depends on the likelihood ratio. This type of proposal schemes was explored in geophysics by Mosegaard & Tarantola (1995), in what is often referred to as the extended Metropolis algorithm. In a high-dimensional target space, the extended Metropolis approach still needs an efficient model proposal scheme (Ruggeri et al. 2015). Following Brunetti & Linde (2017), we use the adaptive multi-chain algorithm DREAM (ZS) (DiffeRential Evolution Adaptive Metropolis using an archive of past states) by Laloy & Vrugt (2012), which is widely used in various geophysical inversion studies (e.g., Bikowski et al. 2012;Rosas-Carbajal et al. 2014;Hunziker et al. 2017). We adapt herein the DREAM (ZS) 's formulation in order to accommodate prior-preserving model proposals.
As an exemplary problem, we consider inference of high-dimensional multi-Gaussian porosity fields using crosshole ground-penetrating radar (GPR) first-arrival travel times. We consider both a linear straight-ray solver, to enable comparisons with analytical solutions, and a more physically-based non-linear eikonal solver. We compare the results obtained by our prior-sampling-based proposal and importance-sampling-based implementation of the (correlated) pseudo-marginal method with standard model proposals and without importance sampling. Furthermore, we compare against the original lithological tomography formulation, full inversion and MCMC inversions that simply ignore the presence of petrophysical prediction uncertainty. With these examples, we will demonstrate that our implementation of the CPM method is outperforming the other inversion methods by greatly enhancing the posterior exploration. This paper is structured as follows. Section 2 introduces the methodology by discussing Bayesian inference in the context of high-dimensional settings, presenting the inversion approaches considered and the tools employed for performance assessment. Section 3 presents the two test examples with linear and non-linear physics. The results and wider implications are discussed in Section 4, followed by conclusions in Section 5.

Methodology
The methodology section starts by introducing the considered latent variable model (Section 2.1), followed by general considerations concerning Bayesian inference and MCMC in high-dimensional settings (Section 2.2). The correlated pseudo-marginal method and our IS procedure are introduced in Section 2.3 and baseline methods used for comparative purposes are presented in Section 2.4. Finally, Section 2.5 presents the performance assessment metrics used to evaluate the results.

Latent variable model
We consider a latent variable model where the unobservable variable X = (X 1 , X 2 , ..., X L ) is related to the d target parameters θ = (θ 1 , θ 2 , ...θ d ) and the T measurements y = (y 1 , y 2 , ..., y T ). We write for G : R L → R T and F : R d → R L with errors ε O and ε P . In our setting, x → G(x) describes the physical forward solver with ε O denoting the observational noise and θ → F(θ) represents the petrophysical relationship with ε P denoting the petrophysical prediction error (PPE). We assume both errors to be Gaussian such that the distribution of X|θ can be represented with the PDF p(x|θ) = ϕ L (x; F(θ), Σ P ) and the one of Y |θ, X with the PDF p(y|θ, x) = ϕ T (y; G(x), Σ Y ), with the notation ϕ M (·; µ, Σ) denoting the PDF of a M -variate Normal distribution with mean µ and covariance matrix Σ.

Bayesian Inference with Markov Chain Monte Carlo
In Bayes' theorem, the posterior probability density function (PDF) p(θ|y) of the model parameters θ given the measurements y is specified by with the prior PDF p(θ) of the model parameters, the likelihood function p(y|θ) and the evidence p(y).
Generally, there is no analytical form of the posterior PDF. If the posterior PDF can be evaluated pointwise up to a normalizing constant, MCMC methods can be used to generate posterior samples. The basic idea of MCMC algorithms is to construct a Markov chain with the posterior PDF of interest as its stationary distribution (see e.g., Robert & Casella 2013). MCMC algorithms iteratively propose new values for the states of the Markov chain that are accepted or rejected with a prescribed probability. One foundational MCMC algorithm is Metropolis-Hastings (MH;Metropolis et al. 1953;Hastings 1970). It proceeds as follows at iteration j: First, using the model proposal density q(·|θ (j−1) ), a new set of states θ (j) is proposed. Then, the acceptance probability, is calculated and the proposed If the proposed θ (j) is rejected, the old state of the chain is kept and θ (j) = θ (j−1) .
Within the MH algorithm, we need to evaluate the likelihood function θ → p(y|θ) in order to compute the acceptance probability. In our latent variable model (see Section 2.1), the likelihood is given by, and the integral has generally no analytical form. In Sections 2.3, 2.4.2 and 2.4.3, we present three methods to circumvent the difficulties of an intractable likelihood function.
We use a high-dimensional pixel-based parameterization of the target field, θ = µ θ + Σ θ 1/2 Z, where Z is a D 2 -dimensional random vector consisting of i.i.d. standard-normal distributed variables. To infer the target field, we need to estimate the Z-variables. Similar to Ruggeri et al. (2015), we do not apply any further dimensionality reduction of the parameter space beyond the discretization (in contrast with, for instance, Brunetti & Linde (2017) who used the dimensionality reduction approach of Laloy et al. (2015)). This is done to avoid distorted posterior PDF estimates that may arise in response to a reduction of the parameter space. Furthermore, we seek to evaluate performance in a challenging high-dimensional setting with thousands of unknowns.
When inferring model parameters with the MH algorithm, it is crucial to choose the model proposal scale well. If the model proposal steps are too large, the acceptance rate is low and the Markov chain needs many iterations until convergence. If the step-width is too small, the exploration of the parameter space is very slow and the Markov chain will similarly need many iterations until convergence (see Section 2.5 for the assessment of convergence). To deal with this challenge of tuning the proposal scale of each model parameter, we use the adaptive multi-chain algorithm DREAM (ZS) (DiffeRential Evolution Adaptive Metropolis using an archive of past states) by Laloy & Vrugt (2012) for which details can be found in Appendix A.
MCMC algorithms generally suffer from the curse of dimensionality as the number of iterations needed for convergence increases with the number of target parameters (e.g., Robert et al. 2018). In the context of Gaussian random fields, Cotter et al. (2013) show that MCMC methods based on standard random walk proposals lead to strong dependencies on the discretization of the target field and to inefficient algorithms when employed in high dimensions. For a given proposal scale, refining the grid representing the random field leads to a decreasing acceptance rate with zero as the limiting value for an infinite number of unknowns. To make MCMC algorithms robust to discretization and maintain a reasonable stepsize when inferring thousands of unknowns, they propose model proposal schemes such as the pCN (preconditioned Crank-Nicholson) that preserve the prior PDF. For a target variable Z with a Standard-Normal prior, the proposal of a standard random walk method is given by Z (j) = Z (j−1) +γζ, with γ being the step size and ζ ∼ N (0, 1), respectively. Instead, the pCN proposal scheme uses Z (j) = 1 − γ 2 Z (j−1) + γζ, ensuring that Z (j) remains standard-normally distributed. Cotter et al. (2013) show that proposal schemes preserving the prior PDF lead to (1) algorithms that mix more rapidly and (2) the convergence being insensitive to the discretization of the target field. We note that the idea of defining a model proposal scheme preserving the prior distribution was proposed more than 25 years ago in geophysics by Mosegaard & Tarantola (1995). This approach is often referred to as the extended Metropolis algorithm and has mainly been explored in the context of inversion with complex geostatistical prior models (a detailed description of the method can be found in Hansen et al. (2012)). Defining a proposal density q(·|θ (j−1) ) such that the MCMC algorithm samples the prior PDF in the absence of data implies that q(θ (j−1) |θ (j) ) holds true, with the implication that the MH acceptance-ratio of Equation (3) is reduced to the likelihood ratio, The extended Metropolis approach still needs an efficient model proposal scheme (Ruggeri et al. 2015), which is why we use DREAM (ZS) in this work. In the case of a Gaussian-distributed prior, the standard DREAM (ZS) proposal scheme does not generate samples that preserve the prior distribution. In order to adapt extended Metropolis to DREAM (ZS) , we rely on a transformation of the variables to the Uniform space (details in Appendix A). This transformation makes it possible to create a proposal mechanism which unites (1) the efficiency of the DREAM (ZS) proposals with (2) the robustness of the prior-preserving proposals. In what follows, our proposal scheme using the uniform transform will be referred to as prior-sampling DREAM (ZS) proposals, while the the standard proposal scheme of DREAM (ZS) will be referred to as standard DREAM (ZS) proposals. We stress that both priorsampling DREAM (ZS) and standard DREAM (ZS) target the same posterior PDF, but the former is expected to be more efficient.

Pseudo-marginal method
Beaumont (2003) shows that a MH algorithm using a non-negative unbiased estimator of the likelihood samples the same target distribution as when using the true likelihood. He exploits this property by estimating the likelihood in Equation (4) on the basis of Monte Carlo averaging over samples of the latent variable X. Andrieu & Roberts (2009) adopt this approach in their pseudo-marginal (PM) method and provide a theoretical analysis of the scheme. When one brute force Monte Carlo sample of the latent variable is drawn in each MCMC iteration without importance sampling (c.f., the original lithological tomography by Bosch (1999); see Section 2.4.2), the algorithm is likely to suffer from a low acceptance rate due to the high variability of the log-likelihood estimator. This is due to the fact that a likelihood estimator given by p(y|θ, X) takes very different values depending on the draw of the latent variable X, even for the same θ. This occurs as the scatter (ε P ) has a strong effect on the data response, and hence, the likelihood. To improve the efficiency, Beaumont (2003) and Andrieu & Roberts (2009) use many samples drawn by importance sampling (IS; e.g. Owen & Zhou 2000). Consequently, they propose the following unbiased estimator of the likelihood p(y|θ), where X n i.i.d ∼ m(·|θ) for n = 1, 2, ..., N with m(·|θ) being the importance density function. More details about the importance sampling procedure will follow in Section 2.3.3.

Correlated pseudo-marginal method
For the PM method to be efficient, the number of samples N used in the likelihood estimator (Eq. (7)) should be selected such that the variance of the log-likelihood ratio estimator is low enough (Doucet et al. 2015). If it is too high, the algorithm will suffer from an impractically low acceptance rate. In the state-space model context, this implies that N needs to scale linearly with T leading to a computational cost of order T 2 at every MCMC iteration, which can be prohibitively expensive for large T (Deligiannidis et al. 2018). To reduce the computational cost, Deligiannidis et al. (2018) introduced the correlated pseudo-marginal (CPM) method by which the draws of latent variables used in the denominator and numerator of the likelihood ratio estimators are correlated. The underlying idea is that the variance of a ratio of estimators is lower if they are positively correlated (Koop 1972). Assuming that the latent variable X is standard-normal distributed, the CPM method proposes (in iteration j) a realization of the n-th latent variable draw by means of pre-conditioned Crank-Nicholson proposals, The assumption that the latent variable has a standard-normal distribution hardly limits the general applicability of the CPM method, since there exist transformations from numerous distributions that will allow proposals to act on Gaussian distributions (e.g. Chen et al. 2018; Section 2.3.3). We stress that if the proposed θ (j) with X (j) n is rejected by the CPM algorithm, we keep X Compared to standard MCMC algorithms, the CPM method requires two additional parameters: the latent variable sample size N and the correlation parameter ρ. To achieve optimal performance, the parameters should be chosen such that the variance of the log-likelihood ratio estimator for a fixed target variable θ, takes values between 1.0 and 2.0 in regions with high probability mass (Deligiannidis et al. 2018).
(y|θ) refer to the likelihood estimators (Eq. (7)) obtained with the accepted latent variable of iteration j − 1 and the proposed (and not necessarily accepted) latent variable of iteration j, that is, the likelihood estimators used in the acceptance ratio of the MH algorithm. In order to choose the parameter values, we first fix the number of samples N at a value that is smaller than the number of available parallel processors. Then, we evaluate different ρ and estimate corresponding values of V ar(R) for a fixed θ in a region with high posterior probability mass (e.g., chosen based on initial MCMC runs).

Importance sampling procedure
For high-dimensional problems with large data sets exhibiting high signal-to-noise ratios, it is necessary to use importance sampling when drawing samples of latent variables to be used within the likelihoodestimator (Eq. (7)). This is a consequence of the integrand p(y|θ, x) in Equation (4) having a peak in a region of X having small probability under p(x|θ). Importance sampling proceeds by sampling from a so-called importance distribution given by the PDF x → m(x|θ) that preferentially generates samples with high p(y|θ, x)p(x|θ). Furthermore, the support of the importance distribution must include all values x, for which p(y|θ, x)p(x|θ) > 0 (Owen & Zhou 2000). It holds, leading to the unbiased importance sampling estimate of the likelihood given in Equation (7). To ensure minimal variance of the estimator, we seek x → m(x|θ) to be nearly proportional to x → p(y|θ, x)p(x|θ) as recalled in Owen & Zhou (2000) referring to the results of Kahn et al. (1953). Since p(x|θ, y) ∝ p(y|θ, x)p(x|θ), it is sensible to base the importance density on x → p(x|θ, y).
Within a latent variable model with a non-linear physical forward solver (Section 2.1), we can not derive the exact expression for p(x|θ, y). Here, we derive local approximations of this posterior by relying on linearization. To do so, we use a linearization of the map x → G(x) around x lin = F(θ lin ) + ε P lin based on a first-order expansion, with J x lin being the Jacobian matrix of the forward solver corresponding to x lin . Ideally, x lin should be given by a realization of the latent variable similar to the one the algorithm is currently exploring. By approximating p(y|θ, and, applying p(x|θ) = ϕ L (x; F(θ), Σ P ) and the relationships between marginal and conditional Gaussians out of Bishop (2006) given in Appendix B, we get, for an approximation of p(x|θ, y). To incorporate importance sampling within the CPM method, we need to correlate the draws of latent variables. To achieve this, we rely on the fact that a realization of the latent variable X can be generated with µ IS + Σ IS 1/2 Z P , where Z P is standard Gaussian distributed in R L . Using this representation, we can correlate the (standard-normal distributed) Z Pvariables using Equation (8).

Baseline inversion methods
We present now the inversion approaches used for comparison with the CPM method. These include a method ignoring the petrophysical prediction errors and two approaches (original formulation of lithological tomography without importance sampling and full inversion) accounting for the PPEs by inferring the joint posterior PDF p(θ, x|y) of the target and latent variables. An overview of all inversion methods (including CPM) is given in Table 1.

Ignore petrophysical prediction errors
This inversion method (no PPE) ignores the presence of petrophysical prediction errors in the MH algorithm. For the latent variable model introduced in Section 2.1, this results in an approximation of the likelihood function with the Gaussian PDFp(y|θ) = ϕ T (y; G(F(θ)), Σ Y ), where the forward response G(F(θ)) is simulated without accounting for PPEs. The method is included in the comparison as it is commonly used in practice as discussed by Brunetti & Linde (2017).

Lithological Tomography
One way to consider PPEs while circumventing the difficulty of an intractable likelihood function is to infer the joint posterior PDF (θ, x) → p(θ, x|y) of the hydrogeological and geophysical parameters. Lithological tomography (Bosch 1999) pursues this strategy and uses a factorization of the joint posterior PDF as p(θ, x|y) ∝ p(θ)p(x|θ)p(y|θ, x), where p(y|θ, x) = p(y|x) is valid for our setting. To sample from this posterior PDF, Bosch (1999) proceeds as follows: First, realizations from the joint prior of θ and X are created by marginal sampling of θ and conditional sampling of X. Then, the pairs of model proposals are accepted or rejected with p(y|x), used in the acceptance ratio of the MH algorithm. In practice, this means that brute force Monte Carlo realizations (no importance sampling) of the petrophysical prediction error ε P are added to the output of the petrophysical relationship F(θ). For our latent variable model, this results in an approximation of the likelihood function witĥ p(y|θ) = ϕ T (y; G(x), Σ Y ), where the latent variable X = F(θ) + ε P is obtained with a draw of ε P from the multivariate Gaussian with PDF ϕ L (·; 0, Σ P ).

Full Inversion
The full inversion approach infers the joint posterior PDF by treating the latent variables analogously to the other unknowns. In the context of our latent variable model (Section 2.1), this means that in iteration j of the MH, not only a new θ (j) but also a new ε P (j) is proposed by the algorithm's proposal scheme. Then the likelihood function p(y|θ, (2017) applied full inversion to infer porosity fields by inversion of crosshole GPR first-arrival travel times, that is, to a setting similar to ours. For the parametrization of the porosity field of interest, they used a spectral representation combined with the dimensionality reduction approach of Laloy et al. (2015). Brunetti & Linde (2017) achieved convincing results and improvements compared to standard lithological tomography without importance sampling (Section 2.4.2). Nevertheless, full inversion is expected to suffer from high dimensionality and strong correlation among the latent and target variables as the two sets of variables are treated as being independent within the proposal scheme (e.g., Deligiannidis et al. 2018).

Performance assessment
To assess the performance of the different inversion approaches, we primarily focus on the exploration of the posterior PDF. The reason for this will become clear in the results section (Section 3).
To declare convergence, we use theR-statistic of Gelman & Rubin (1992) that compares the withinchain variance with the between-chain variance for the second half of the MCMC chains. The general convention is that convergence is declared once this statistic is smaller or equal to 1.2 for all model parameters. Since we deal with a high-dimensional parameter space with thousands of unknowns, we relax this condition slightly and declare convergence if 99 % of the parameters satisfy this criterion. When an algorithm is considered convergent, we compare the resulting posterior samples with those of the other approaches.
For the test case with linear physics in Section 3.2, we compare the results with the analytical solution of the posterior PDF p(θ|y). For these comparisons, we use histograms and the Kullback-Leibler divergence (KL -divergence; Kullback & Leibler 1951). The KL -divergence between two PDFs p 1 (·) and p 2 (·) is defined as, To obtain the PDF of the estimated posterior, we can use the MCMC samples to either (1) make a kernel density estimate or to (2) estimate the mean and variance for a Gaussian approximation (Krueger et al. 2016). Here we use the second option since the posterior is Gaussian. If the PDFs p 1 (·) and p 2 (·) are Gaussians with p 1 = N (µ 1 , σ 2 1 ) and p 2 = N (µ 2 , σ 2 2 ), the expression of the KL-divergence reduces to, A KL-divergence of zero indicates that the two PDFs are equal and it increases as the distributions diverge from each other.
For the test example with non-linear physics in Section 3.3, there is no analytical solution to compare with. Hence, we compare the estimated posterior distribution with a single value (the known true porosity at each pixel). We achieve this by applying so-called scoring rules (Gneiting & Raftery 2007) assessing the accuracy of a predictive PDF θ →p(θ) with respect to a true value θ. Scoring rules are functions that assign a numerical score for each prediction-observation pair (p, θ), with a smaller score indicating a better prediction. They assess both the statistical consistency between predictions and observations (calibration) and the sharpness of the prediction. We use the logarithmic score (logS; Good 1952) defined by logS(p, θ) = − logp(θ) that is related to the Kullback-Leibler divergence (Gneiting & Raftery 2007). As for the linear case, we use the MCMC samples to obtain a Gaussian approximation of the estimated posterior PDF. The logarithmic score favours predictive PDFs under which the true value has high probability. We supplement this measure with two simpler ones: the number of pixels in which the true porosity value was in the range of the posterior samples and the standard deviation of the estimated posterior PDF.
We also consider the acceptance rates (AR) and the integrated autocorrelation time (IACT). We aim for an acceptance rate of 15% -30% as proposed by Vrugt (2016). The IACT of the chain {θ (j) ; j = 1, 2, ...} Corr(θ (1) , θ (1+l) ). In practice, the estimated autocorrelation for large values of l is noisy such that we need to truncate the sum. Following Gelman et al. (2004), we truncate the sum when two successive autocorrelation estimates are negative. We renounce from discussing the CPU time as it depends strongly on the chosen forward model and discretization as well as on other parameters pertaining to the computing equipment.

Results
We consider the problem of inferring the porosity distribution using crosshole GPR first-arrival travel times. We first address a test case with linear physics (straight-rays) to allow for comparison with analytical solutions and then one with non-linear physics (eikonal solver) to address a more challenging and physically-based setup. Our examples are synthetic and the water-saturated porosity field is described by a multi-Gaussian random field.

Synthetic data generation
Our considered subsurface domain is 7.2 m × 7.2 m and we use 25 equidistant GPR transmitters located on the left side and 25 receivers on the right side of the model domain, resulting in 625 first-arrival travel times. The transmitter-receiver layout is depicted in Figure 1c. As introduced in Section 2.2.1, we assume the porosity field to be a Gaussian random field GRF (µ θ (·), C θ (·, ·)). We use µ θ (·) = 0.39 and an exponential covariance function C θ (·, ·). For the latter, we use a sill of 2e −4 and geometric anisotropy where the main, horizontal direction has an integral scale of 4.5 m and the integral scale ratio between the horizontal and vertical direction is 0.13. We use a (50 × 50)-dimensional pixel-based parameterization of the porosity field; the true synthetically generated field is shown in Figure 1a. Note that porosity is a positive quantity bounded between zero and one while a Gaussian prior distribution has a full support. The Gaussian prior is used here to ensure an analytical solution in the linear physics case. Given the presented mean and the sill, it is extremely unlikely that a porosity value outside the physical boundaries is generated. In other settings, one could use a transform of the porosity (e.g., as in Bosch 2004) or choose a bounded distribution.
To predict the dielectric constant κ, we use the complex refractive index model (CRIM; Roth et al. 1990), Figure 1: (a) Porosity field θ, (b) PPE field ε P , (c) slowness field x with transmitter-receiver layout, (d) dependency of slowness on porosity obtained without (line) and with (scatter) PPE and (e) noisecontaminated first-arrival travel times y for the linear and the non-linear forward solver corresponding to the true synthetic model.
where κ w and κ s are the dielectric constants of water [81] and mineral grains [5], respectively. The resulting slowness field (which in our case is the latent variable X) depicted in Figure 1c is given by, where c is the speed of light in vacuum [0.3 m/ns]. This specifies the petrophysical relationship to be linear with θ → F(θ) = 1 We add a petrophysical prediction error (PPE) ε P that is a realization of a centred GRF over a regular 2D grid of size 50 × 50. We are assuming that the PPE field (depicted in Figure 1b) has an exponential covariance function C P (·, ·) with a sill of 2.1e −2 and the same correlation structure as the porosity field. The dependency of the slowness on the value of the porosity and the PPE is indicated in Figure 1d. Finally, the resulting 625 GPR first -arrival travel times are calculated with (i) a linear (straight-ray) forward solver referred to as G s and (ii) a non-linear (eikonal) forward solver referred to as G e (the time2D solver of Podvin & Lecomte (1991)), such that, with i.i.d. centered normal observational noise ε O with standard deviation of 1 ns. The two sets of traveltimes are depicted in Figure 1e.

Inversion settings and prior assumptions
All considered inversion methods (Sections 2.3 and 2.4) are implemented with prior-sampling and standard DREAM (ZS) proposals using the same parameter settings of the DREAM (ZS) algorithm with four MCMC chains running in parallel. For the prior on porosity, we use the Gaussian PDF p(θ) = ϕ 2500 (θ; µ θ , Σ θ ) assuming the mean µ θ and covariance structure Σ θ to be known (the same values as for the data generation). Using a pixel-based parameterization of the field, we infer the 2500-dimensional vector Z defining the porosity by θ = µ θ + Σ θ 1/2 Z, with Z having a multivariate standard-normal prior PDF. The full inversion has to estimate another 2500 Z P -variables for the PPE field leading to a total of 5000 inferred parameters. For the PPE ε P we also use a Gaussian prior PDF p(ε P ) = ϕ 2500 (θ; 0, Σ P ) with known covariance structure Σ P , leading to a Gaussian prior PDF for the slowness field (for fixed porosity) given by p(x|θ) = ϕ 2500 (x; F(θ), Σ P ). For the likelihood function, we assume that the 625-dimensional vector describing the observational noise ε O has a Gaussian distribution with zero mean and diagonal covariance matrix Σ Y ; the standard deviation is assumed to be 1 ns as in the data generation process.

Linear physics
To enable comparisons of the inferred posterior PDFs with the analytical solution for p(θ|y), we first consider the case of linear physics. Then, with J s being the Jacobian (i.e., forward operator) of the linear forward solver. The analytical posterior PDF can be derived as detailed in Appendix B. Figure 2a shows the posterior mean and Figures 2b -2d depict three draws from the posterior distribution.
When employing the PM and CPM method in this setting of large datasets with low noise, it is crucial to use a well-chosen importance sampling for the latent variable. As introduced in Section 2.3.3, it is sensible to use x → p(x|θ, y) as a basis for the importance density. As long as we are in the linear Gaussian case, we can derive the analytical expression for this posterior (Appendix B), resulting in a zero-variance importance sampling density (Owen & Zhou 2000). Since it then does not make sense to use multiple importance density samples (the importance weights are constant), we combine in this linear case importance sampling with PM using N =1 (original lithological tomography algorithm enhanced with importance sampling that we will hereafter refer to as LithTom IS). We note that using the exact formula for the importance sampling corresponds to having access to the exact likelihood p(θ|y). The use of larger N is considered in Section 3.3 for the case of non-linear physics. This linear setting for which analytical solutions are available serves mainly (1) to demonstrate the necessity of a well-working importance sampling distribution, (2) to investigate the exploration capabilities of MCMC-based inversion approaches that estimate the intractable likelihood using Monte Carlo samples (lithological tomography, PM and CPM methods) and (3) to compare the performances of the prior-sampling and standard DREAM (ZS) proposal mechanisms. Figure 3 presents the estimated posterior means of the porosity field obtained when applying the no PPE (Fig. 3a), the full inversion (Fig. 3b) and the LithTom IS ( Fig. 3c) with standard DREAM (ZS) proposals, as well as for LithTom IS with prior-sampling DREAM (ZS) proposals (Fig. 3d). These are the cases for which we reached convergence of the chains. The porosity field obtained with the inversion ignoring PPEs has, as expected (Brunetti & Linde 2017), a higher variance. Visually, all other estimates are very similar in terms of structure and magnitude with respect to the analytical posterior mean in Figure 2a. The estimated posterior mean of LithTom IS with the prior-sampling DREAM (ZS) proposals has a slightly lower variance than for standard DREAM (ZS) proposals. The ARs (Table 2) for standard DREAM (ZS) proposals are the highest for LithTom IS, while the method ignoring PPEs and full inversion have lower ARs. Classical lithological tomography without importance sampling leads to an AR of less than 0.1 % such that, in practice, it unfeasible to reach convergence. Applying the CPM method without IS for N =50 and ρ=0.95 also results in an only slightly larger AR (roughly 0.2 %), thereby, highlighting the need for importance sampling for the considered problem. Since less than 5 % of the parameters converged after 200'000 iterations, we renounce from showing further results for the CPM and PM method without IS. The method ignoring PPEs and the full inversion using prior-sampling DREAM (ZS) proposals suffer from very low ARs and did not reach convergence after 200'000 iterations. Table 2 shows the number of iterations needed for the 99 th percentile of the parameters'R-statics to be below 1.2. It also shows the IACTs of the cell in the very middle of the porosity field for all inversion approaches reaching convergence within 200'000 MCMC iterations. We observe that the iterations needed for convergence and the IACT of the LithTom IS method with prior-sampling DREAM (ZS) proposals are the lowest. Figure 3e shows the evolving log-likelihood values. When ignoring PPEs or performing the full inversion, the chains converge to much higher log-likelihoods than for the LithTom IS method. This is expected as they rely on the likelihood p(y|θ, x) (where X = F(θ) + ε P , with ε P = 0 for the algorithm ignoring PPEs), while LithTom IS estimates p(y|θ) = p(y|θ, x)p(x|θ)dx. This example highlights that LithTom IS broadens the likelihood function. Figures 3f and 3g show the prior probabilities (logarithmically transformed) for the posterior samples obtained with the three different inversion approaches using the two alternative proposal schemes. We observe that the LithTom IS method using prior-sampling DREAM (ZS) proposals (Fig. 3f) is the only approach for which the prior probability of the true porosity field is sampled. All other methods and proposal scheme combinations sample porosity fields with higher prior probabilities than the true field (black solid line). Practically speaking, this implies for these cases that none of the posterior samples are close to the true model. full inversion confirm the significant improvements of the exploration capabilities of this approach. An even better representation of the analytical posterior was obtained with the LithTom IS approach when using prior-sampling DREAM (ZS) proposals. This is indicated by the histogram in Figure 4d and by a further two-fold decrease of the KL-divergence in Figure 5d. An overview of the mean KL-divergences is given in Table 2.
This linear example has been used to show that importance sampling and prior-preserving proposal schemes are essential to obtain meaningful results in our considered high-dimensional setting. For this example, one can get accurate results using LithTom IS alone. The next section dealing with the non-linear case will serve to demonstrate the benefits of the CPM method in non-linear settings.

Non-linear physics
We now consider a non-linear test case in which the 625 arrival times are generated with the eikonal 2D traveltime solver time2D of Podvin & Lecomte (1991) such that, Given the non-linear physics, the likelihood function p(y|θ) is intractable and there is no analytical expression for the posterior PDF p(θ|y) to compare with. The same applies for the PDF p(x|y, θ)  that we previously used for the importance sampling of the latent variable X. Hence, as importance sampling distribution we rely on the approximation of the PDF p(x|y, θ) introduced in Section 2.3.3. For x lin = F(θ lin ) + ε P lin = 1 c √ κ s + ( √ κ w − √ κ s )θ lin + ε P lin , we use the last state of the porosity field for θ lin and the previous importance sampling mean µ IS for ε P lin . To decrease computational resources, we only update the linearization every 100 MCMC iterations. Since the expression is approximate, we further inflate the importance sampling covariance matrix Σ IS by multiplying Σ Y with a factor. After initial testing, we found that 1.2 yielded the best performance. Figure 6 depicts the dependence of the variance of the log-likelihood ratio estimator R (Eq. (9)) on the correlation parameter ρ for N = 1, N = 10 and N = 50 samples of the latent variable X (with θ being fixed at a region with high posterior probability mass). Figure 6a depicts estimates when drawing the realizations of the latent variable proportionally to its prior distribution p(x|θ) and Figure 6b for the case where the latent variable is sampled with importance sampling. The two plots highlight three fundamental aspects of the CPM method in our geophysical setting. First, it is crucial to use a well chosen importance sampling for the latent variable draws, since for a correlation of, say, ρ = 0, the variance of the log likelihood ratio estimator can be reduced from values between 10 000 and 1 000 000 (using sampling from prior) to values between 3 and 31 (using importance sampling). Second, increasing the number of draws of latent variables (N ) decreases the variance of the log-likelihood ratio estimator further and, third, this is also achieved by increasing the amount of correlation (ρ) used for two subsequent draws of latent variables. The variance for ρ = 1 is equal to zero for all parameter settings (as we use the same values for X (j−1) and X (j) ). Without importance sampling, we could still obtain a variance of the log-likelihood ratio estimator between 1 and 2 as recommended by Deligiannidis et al. (2018), but with the need of a very high N or a ρ very close to 1. In practice, this would either result in excessively high computational costs or slow mixing in the draws of the latent variables.
Due to the high variances displayed in Figure 6a and since the pseudo-marginal approaches without importance sampling have already proven to be highly inefficient in the linear case (Table 2), we now restrict ourselves only to CPM implementations involving IS. In stark contrast to the linear case, the LithTom IS approach (N = 1, ρ = 0) leads to a highly inefficient algorithm, as the variance of R around 30 is much higher than the upper recommended threshold of 2.0. For the CPM method, we set the number of samples to 10 and the correlation to ρ = 0.95 as this values leads to a variance of the log likelihood ratio estimator in-between 1.0 and 2.0. The autocorrelation of one cell of the latent variable field is given by Corr(X 1 , X 1+l ) = ρ l for lag l with the correlation mechanism of Equation (8), such that for ρ = 0.95 roughly 100 (accepted) iterations are needed to draw an independent realization of the latent variable. In practice, the decorrelation will be slower as we only move on with accepted proposals (Section 2.3.2) .
The results for both DREAM (ZS) proposal schemes are shown in Figure 7 and Table 3. For the estimates of the posterior mean of the porosity field ( Fig. 7a-7d), we observe similar results as in the linear case: Using prior-sampling DREAM (ZS) proposals results in a porosity field estimate with lower variance and using the method ignoring PPEs (Fig. 7a for standard DREAM (ZS) proposals) leads to higher variance. The highest acceptance rate is obtained with applying the CPM IS method using standard DREAM (ZS) proposals (Table 3) and the acceptance rates for prior-sampling DREAM (ZS) proposals are lower. The LithTom approach with IS has an AR of less than 1 % and would, therefore, require far more than 200'000 iterations to converge. Trace plots of the evolving log-likelihood values are shown in Figure 7e. As expected and in agreement with the linear test case (Fig. 3e), the (a) (b) Figure 6: Variance of the log likelihood ratio estimator R = log p (j) for the non-linear test example and θ fixed at a region with high posterior probability mass as a function of ρ (used to correlate the latent variables X (j) and X (j−1) as in Equation (8)) for N = 1, N = 10 and N = 50 samples of the latent variable X; the realizations of the latent variable are drawn (a) from the prior p(x|θ) and (b) with importance sampling. The black lines delimit the range between 1.0 and 2.0 recommended by Deligiannidis et al. (2018). methods converge to different values. As in the linear case, we find that CPM IS with prior-sampling DREAM (ZS) proposals is the only case providing posterior samples that match the prior probability of the true porosity field (Fig. 7f and 7g).  Table 3. We find that the method that ignores PPEs fails to sample a range of values including the true porosity value in more than 10% of the pixels and has a mean estimated posterior standard deviation that is up to 50 % smaller than the other methods. The CPM IS method generates posterior samples with ranges that include, in more than 99 % of the pixels, the true porosity value with the percentages obtained using prior-sampling DREAM (ZS) proposals being even higher. Finally, the full inversion does not sample the true porosity value in almost 3% of the pixels and has a reduced mean estimated posterior standard deviation by up to 40 % compared to the CPM IS method. We also note that the IACT of the CPM methods are the lowest (Table 3).

Discussion
This study showed clearly that the correlated pseudo-marginal (CPM) method, which accounts for petrophysical prediction uncertainty within the estimate of the likelihood function p(y|θ), combined with importance sampling (IS) and prior-sampling MCMC proposals leads to a broader exploration of the target posterior p(θ|y) than the other presented combinations of inversion methods and proposal schemes. The CPM method is an exact and general method, but it needs in the considered highdimensional setting an efficient importance sampling and prior-sampling proposals to work well even for the case of linear physics. In the linear setting (with available analytical solutions for the PDFs), the CPM method using importance sampling performs well using only one uncorrelated sample of the PPE (LithTom IS). In absence of importance sampling, even a high number of samples N and correlation ρ could not prevent the algorithm from being highly inefficient (Table 2). We find that the exploration of the posterior PDF is much improved when using the LithTom IS approach compared with full inversion (Fig. 4 and Fig.  5). Although theR-statistic of Gelman & Rubin (1992) suggests that the full inversion algorithm (using standard DREAM (ZS) proposals) has converged, we demonstrate a significant underestimation of the posterior standard deviation and posterior samples with far too high prior probabilities compared with the true model ( Fig. 3f and 3g). Indeed, the full inversion's high acceptance rate (for standard DREAM (ZS) proposals) may be mainly a consequence of local exploration combined with an adaptive MCMC expanding its archive. This (1) points out that Gelman-Rubin'sR-statistics and the acceptance rate are insufficient metrics to assess the performance of an adaptive MCMC algorithm such as DREAM (ZS) and (2) highlights issues with over-fitting when using adaptive MCMC. Indeed, Robert et al. (2018) warn against using adaptive MCMC methods without due caution as adaptations to the proposal scheme can lead to algorithms relying too much on previous iterations, thereby, excluding parts of the parameter space that have not yet been explored.
The need for a well-chosen importance sampling distribution is also demonstrated for the non-linear setting by analysing the variances of the log-likelihood ratio estimator (Fig. 6). This analysis also confirmed the strong influence of N and ρ. Since the importance sampling distribution is no longer showing the number of iterations needed for the 99 th percentile of the parameter'sR-statics to be below 1.2 (or the percentage of parameters with aR-statistics below 1.2 if the the algorithm did not converge), the percentage of pixels in which the true porosity value lies within the range of posterior samples (θ true ), the mean logarithmic score (logS), posterior standard deviation (Post SD) and the integrated autocorrelation time (IACT) for the cell in the very middle of the porosity field θ. The CPM IS method was evaluated with the parameter choice of N = 10 and ρ = 0.95. proposals) is better than the full inversion by observing that (1) the range of the posterior samples includes more often the true porosity value while (2) the logarithmic score is lower and (3) the mean estimated posterior standard deviation is higher (Table 3).

Method
We recommend to work in the full parameter space whenever possible such that any distortions in the posterior estimations due to model reductions can be avoided. The presented adaptive priorpreserving proposal scheme (prior-sampling DREAM (ZS) proposal) is developed in the spirit of the extended Metropolis algorithm of Mosegaard & Tarantola (1995) and the pCN proposal of Cotter et al. (2013). It is a simple correction of the standard DREAM (ZS) proposal that (1) makes the algorithm robust to the choice of the discretization of the target field and (2) maintains its capabilities to sample efficiently in complex high-dimensional parameter spaces. We find that the prior-sampling DREAM (ZS) proposals lead to an enhanced exploration of the posterior PDF and a stable AR (Tables  2 and 3). Indeed, the CPM IS approach with prior-sampling proposals is the only one generating samples with a prior probability comparable to the one of the true porosity field (Figs. 3 and 7). Due to dependencies between latent and target variables, the full inversion with prior-sampling DREAM (ZS) proposals suffers from a very low acceptance rate as the method does not allow for large proposal steps. This dependency is bypassed by the CPM IS, allowing larger steps for a given AR. In general, combinations of adaptive Metropolis and pCN-proposals are referred to as DIAM (dimension independent adaptive Metropolis) proposals and were introduced by Chen et al. (2016). Another way to increase the efficiency of the pCN proposal was proposed by Rudolf & Sprungk (2018) with the so-called generalized pCN-proposal (gpCN), in which the proposal scheme is tuned to have the same covariance as the target posterior distribution.
We emphasize that this study only considers synthetic data. We demonstrate that all but our method of choice (CPM IS with prior-sampling DREAM (ZS) proposals) have severe problems in exploring the full posterior distribution even in this well-specified setting. A field demonstration of CPM IS with prior-sampling DREAM (ZS) proposals is a natural next step. Furthermore, our entire study remains within Gaussian assumptions for the target field, petrophysical prediction uncertainty and observational noise. In the presented results, we deal only with weak non-linearity in our forward operator and assume the petrophysical relationship to be linear. In the future, it would be useful to consider test cases involving stronger non-linearity, be it through a higher variance of the slowness field or a non-linear petrophysical relationship. Stronger non-linearity would affect the accuracy of the first-order expansion used to derive the importance sampling distribution for the CPM method, implying that the approximations would become less accurate. This could lead to a decrease of efficiency that could be counter-acted by using larger N or ρ. An important topic for future research would be to develop and assess importance sampling schemes that do not rely on Gaussian assumptions. Potential starting points could be efficient importance sampling by Richard & Zhang (2007) or multiple importance sampling introduced by Veach & Guibas (1995) and popularised by Owen & Zhou (2000).
In agreement with Brunetti & Linde (2017), we find that ignoring petrophysical prediction uncertainty leads to biased estimates and too tight uncertainty bounds. While the need for a method accounting for PPEs grows with increasing integral scale of the target field (Brunetti & Linde 2017), the ratio of the variances of the PPE, the target variable and the observational noise also influences the results. The need for a well-working importance sampling for CPM grows with increasing petrophysical prediction uncertainty and decreasing observational noise. At the same time, large petrophysical prediction uncertainty leads to a flattened likelihood function p(y|θ), thereby, decreasing the variance of the likelihood estimators (assuming a well-working importance sampling) and, therefore, enhancing the efficiency of the algorithm. Our present work focuses on petrophysical prediction uncertainty for a known covariance model, but it would be possible to expand this to an unknown covariance model, an uncertain petrophysical model or uncertain model parameters.

Conclusions
We consider lithological tomography in which geophysical data are used to infer the posterior PDF of target (hydro)geological parameters. In such a latent variable model, the geophysical properties play the role of latent variables that are linked to the properties of interest through petrophysical relationships exhibiting significant scatter. Compared with the original formulation of lithological tomography that does not consider importance sampling, we make the approach more applicable to high dimensions (thousands of unknowns) and large data sets with high signal-to-noise ratios. To account for the intractable likelihood appearing in the Metropolis-Hastings algorithm in this setting, we explore the correlated pseudo-marginal (CPM) method using an importance sampling distribution and priorsampling proposals. For the latter, we adapt the standard (adaptive) proposal scheme of DREAM (ZS) with a prior-sampling approach, leading to a further improvement in exploration compared with standard model proposals when dealing with high-dimensional problems. We find that our implementation of the CPM method outperforms standard lithological tomography and the full inversion approach, which parameterizes and infers the posterior petrophysical prediction uncertainty. For a linear test example, the mean KL-divergence with respect to the analytical posterior can be reduced by 99 % by our implementation of the CPM method (even without using correlations) compared with full inversion. In the case of non-linear physics, we reduce the mean logarithmic score with respect to the true porosity field by up to 33 % compared with the full inversion method. The CPM method is generally applicable and accurate, but it requires a well-working importance sampling distribution (presently based on Gaussian random field theory) to be efficient. Future work with the CPM method could consider field data applications, more non-linear physics and non-linear petrophysical relationships as well as relaxing the assumptions of Gaussian random fields. Furthermore, the method's use in coupled hydrogeophysical inversions involving hydrogeological flow and transport models would be of interest.
A DREAM algorithms and prior-sampling proposals To perform a high-dimensional inversion with the MH algorithm, one needs a well-working proposal scheme. To deal with this challenge, Ter Braak (2006) introduced an adaptive random walk MH algorithm named Differential Evolution Markov chain (DE-MC). This method runs C Markov chains in parallel, where at each iteration j, the C different realizations of the model parameters define a population {Z (j) c ; c = 1, 2, ..., C}, which is used to guide new model proposals. For chain c, two chains (denoted as a and b) are drawn without replacement from the remaining set of chains. Then, the algorithm proposes a new state for the c-th chain with, where γ denotes the jumping rate and ζ is a draw from N (0, s 2 ) with a small standard deviation s used to ensure that the resulting Markov chain is irreducible. By accepting or rejecting the resulting proposals with the MH-ratio of Equation (3), a Markov chain with the posterior PDF as its stationary distribution is obtained (Proof in Vrugt et al. 2008a). This leads to an algorithm which is automatically adapting the scale and the orientation of the proposal density along the way to the stationary distribution, allowing it to provide efficient sampling on complex, high-dimensional, and multi-modal target distributions. Based on the DE-MC, Vrugt et al. (2008b) introduced the adaptive multi-chain MCMC algorithm called DREAM (DiffeRential Evolution Adaptive Metropolis). It enhances the efficiency of DE-MC by applying subspace sampling (only randomly selected dimensions of the model parameter are updated) and outlier chain correction. An excellent overview of the theory and application of the DREAM algorithm is given by Vrugt (2016). For our case study, we use the extended version DREAM (ZS) introduced by Laloy & Vrugt (2012), as its proposal scheme using an archive of past states leads to further improved convergence and posterior exploration.
To adapt extended Metropolis to DREAM (ZS) , we rely on a transformation of the variables to the Uniform space. In our case study with Gaussian target variable Z ), we hence force the algorithm to use a proposal scheme that samples from the prior PDF.

B Analytical posterior PDF and importance density for linear physics
Assuming linear physics and petrophysics, it is possible to derive an analytical expression for the posterior PDF p(θ|y) of the porosity (or other variable of interest). We consider here both relationships being linear without intercept (G(X) = J s X and F(θ) = J p θ), however, an intercept (as the one used for F(θ) in our test case; Section 3.1.1) is easily included. For the 2D grid of the porosity θ and the latent variable X, we use the following prior PDFs: p(θ) = ϕ D 2 (θ; µ θ , Σ θ ), p(x|θ) = ϕ L (x; J p θ, Σ P ).
To derive the (in this case) tractable likelihood p(y|θ), we use a standard result about marginal and conditional Gaussians (Bishop 2006):

Lemma 1 Marginal and Conditional Gaussians
Assume a marginal Gaussian distribution for X ∈ R L and a conditional Gaussian distribution for Y ∈ R T given X in the form p(x) = ϕ T (x; µ, Λ −1 ), p(y|x) = ϕ T (y; Ax + b, L −1 ), with ϕ T (·; µ, K) denoting the PDF of the T -variate Normal distribution with mean µ and covariance matrix K. Then, the marginal distribution of Y and the conditional distribution of X given Y are given by where Σ = (Λ + A T LA) −1 .
Using the prior on the latent variable X and the Gaussian likelihood p(y|x, θ) = ϕ 625 (y; J s x, Σ Y ), we get with Equation (22), p(y|θ) = ϕ T (y; J s J p θ, Σ Y + J s Σ P J T s ).
Subsequently, the analytical form of the posterior p(θ|y) is derived with Equation (23), the prior on porosity and the expression of the likelihood p(y|θ) from the last equation: p(θ|y) = ϕ D 2 θ; µ θ|Y , Σ θ|Y , For the case with linear physics, the importance density p(x|θ, y) = ϕ L (x; µ IS , Σ IS ) introduced in Section 2.3.3 is an exact expression for p(X|θ, y) and the IS mean and covariance matrix reduce to: