# Transdimensional joint inversion of flow and well log data using a cascaded Metropolis sampler on a layer-cake model

Julien Herrero, Guillaume Caumon, Thomas Bodin, Jeremie Giraud

► To cite this version:

Julien Herrero, Guillaume Caumon, Thomas Bodin, Jeremie Giraud. Transdimensional joint inversion of flow and well log data using a cascaded Metropolis sampler on a layer-cake model. Geoenergy Science and Engineering, 2025, 246, pp.213605. 10.1016/j.geoen.2024.213605. hal-04873511

HAL Id: hal-04873511

https://hal.univ-lorraine.fr/hal-04873511v1

Submitted on 8 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Transdimensional joint inversion of flow and well log data using a cascaded Metropolis sampler on a layer-cake model

Julien Herrero<sup>a,*</sup>, Guillaume Caumon<sup>a,b</sup>, Thomas Bodin<sup>c</sup>, Jeremie Giraud<sup>a,d</sup>

## ARTICLE INFO

**Keywords:**

History matching

Inverse problem

Markov chain Monte Carlo

Uncertainty quantification

Reservoir modeling

## ABSTRACT

History matching is a type of inverse problem for subsurface geore source evaluation that adjusts the geological and petrophysical properties of a reservoir from dynamic data. Here, we aim to infer the geological structures and the associated permeability values of a model in a joint inversion using observed production data (pressure and fluid saturation) and permeability log data. To determine the appropriate level of complexity in the unknown model, we introduce a transdimensional Monte Carlo method applied to a one-dimensional layered reservoir model. In contrast to classical Bayesian modeling in which permeability values are predicted in a fixed reservoir geometry, the number of layers is treated as an unknown variable in the inversion. Layers remain horizontal for simplicity and are characterized only by their number, thickness, interface depth, and constant isotropic permeability values. Our method employs a cascaded Metropolis sampler based on a two-step sequential acceptance criterion in the Markov chain to reduce computational costs. The likelihood probability of well log data is first quickly computed as in a standard regression problem, and flow simulations are performed only on candidate models accepted by this first regression step to further reduce the uncertainty. Numerical flow simulations are solved on a mesh conformal to horizons, which is locally updated at each iteration. Synthetic tests on a simple reference model highlight the capability of the transdimensional cascaded Metropolis algorithm, which jointly uses flow and well log data, to recover a parsimonious representation of the subsurface. Results show that the proposed approach is able to successfully locate major geological discontinuities, quantify uncertainty, and spatially capture the level of information brought by different datasets in the joint inversion process. This suggests potential applicability to real data, more complex reservoir geometries, and other physics.

## Introduction

Subsurface exploration for estimation or storage of georesources involves significant scientific, financial, environmental and societal decisions that require studies based on numerical simulations to explain and predict available observations. In the context of all reservoir issues for the safe exploitation of hydrocarbons, CO₂ sequestration, recovery of geothermal resources, or gas storage, numerical simulations are essential for making predictions and establishing scenarios to inject or produce these different reservoir fluids. This leads to a number of inverse problems where observations are used to determine the values of the variables which describe the system.

In the case of a reservoir model or an aquifer defined by a set of layers, it is common to update the knowledge about rock property

fields (e.g., rock permeability) by performing history matching of production data (Oliver et al., 2008). History matching aims to estimate reservoir model parameters and their uncertainty from production data (e.g., flow rate, pressure, fluid compositions measured at the wells) given the forward relation between the parameters and the data. In a history matching problem, the number of model parameters to be estimated is generally larger than the number of observed data to be matched, and the relationship between the production data and the model variables is complex and non-linear. The history matching problem is therefore ill-posed, and numerous strategies have been proposed in the last decades to make the problem tractable (e.g., Christie et al., 2006; Oliver and Chen, 2011; Rwechungura et al., 2011; Liu and Durlofsky, 2020). Among these, stochastic Bayesian approaches are well suited to address the non-linearity and ill-posedness of the

<footnote><sup>a</sup> Université de Lorraine, CNRS, GeoRessources, F-54000 Nancy, France</footnote>

<footnote><sup>b</sup> Institut Universitaire de France, France</footnote>

<footnote><sup>c</sup> Laboratoire de Geologie de Lyon: Terre, Planetes et Environnement, CNRS, ENSL, Univ Lyon 1, 69100, Villeurbanne, France</footnote>

<footnote><sup>d</sup> Centre for Exploration Targeting (School of Earth Sciences), The University of Western Australia, 35 Stirling Highway, Perth, Australia</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: julien.herrero@univ-lorraine.fr (J. Herrero).</footnote>
problem (e.g., Gómez-Hernández et al., 1997; Barker et al., 2001). Bayesian formulation also allows for integration of diverse data types, each carrying its own uncertainties and sensitivities (Bodin et al., 2012b).

This study aims to provide probabilistic estimates of the geometry and permeability values of a 1D layered and non-periodic medium using a Markov chain Monte Carlo (MCMC) algorithm and two different sets of reservoir measurements in a joint inverse problem: Static (well logs of permeability) data and dynamic (history matching) data (Section 1). The common reservoir modeling practice of using a prescribed number of layers of fixed geometry neglects geological interpretation uncertainty (Lallier et al., 2012; Wellmann and Caumon, 2018) and can lead to forecasting biases (Carrera et al., 2005). Therefore, we allow the number and the geometry of layers of the reservoir model to self-adapt over MCMC iterations. This approach increases the degrees of freedom in the problem, which allows for the joint characterization of both the geometric and petrophysical uncertainties in the subsurface inverse problem. Under this framework, the complexity of the model is constrained only by reservoir data, eliminating the need for restrictive geological assumptions in the prior probability. To sample the unknown numbers of geological layers, we employ the reversible jump Markov chain Monte Carlo (RJMCMC) algorithm, which is a global search method based on a transdimensional Markov chain (Section 2). This algorithm was originally presented by Geyer and Møller (1994) and Green (1995) as an extension of the standard Metropolis-Hastings MCMC sampler (Metropolis et al., 1953; Hastings, 1970) to variable dimension problems.

To date, only a few studies have considered transdimensional settings for porous medium problems. Mondal et al. (2010) pioneered the use of RJMCMC for uncertainty quantification of flows in heterogeneous porous media. They reconstruct the shape of channel boundaries combining a level-set representation of these boundaries and a stochastic two-point permeability model constrained by well log values of permeability. Channel boundaries are parameterized with piecewise linear functions of variable dimension involving point addition, deletion, or vertical move. Jiménez et al. (2016) introduce a “smart” pilot point method, where not only the values but also the number and optimal locations of these points are determined through a transdimensional Bayesian inversion, demonstrating its potential for high-resolution hydrogeological model calibration. Somogyvári et al. (2017) use a transdimensional inversion to reconstruct 2D fracture network geometries by using tracer tomography synthetic data directly injected from wells in the fractures. The model parameters are the fracture locations and lengths, and fractures are added or deleted during the reversible jump iterations. This method was then extended to 3D by Ringel et al. (2021). In geophysics, since the seminal work of Malinverno and Leaney (2000) for the inversion of zero-offset vertical seismic data, transdimensional methods have been applied several times to infer the number of layers in the case of 1D (horizontal) multilayer models (e.g., Malinverno, 2002; Dettmer et al., 2010; Bodin et al., 2012b; Seillé and Visser, 2020) or 2D models adding a variable parameter of interface slope (Visser et al., 2019).

Our work builds on these approaches: Similar to the geophysical applications mentioned above, we implement a variable number of geological layers with a piecewise constant parameterization on layer interfaces to build layer-cake models. Similar to Mondal et al. (2010), we solve an inverse problem where both the model geometry and permeability property are updated using flow data and well log permeability information. The joint inversion is done with a sequential RJMCMC sampler based on the cascaded Metropolis algorithm described by Tarantola (2005) (page 182, chapter 6). Although cascaded Metropolis samplers have been previously used in geosciences (e.g., by Manassero et al. (2021) for magnetotelluric data inversion), their application in a transdimensional context is novel. It is especially convenient in our context as the forward problem associated with well logs is computationally cheap as compared to the flow simulation.

Hence, this problem is evaluated first while the more time-consuming fluid flow simulation is only performed if the first acceptance criterion is met. We also use an unstructured mesh model representation to provide the needed degree of geometric accuracy during flow simulation. Because the mesh discretization must stay conformal to layer boundaries to ensure the stability of the forward problem throughout transdimensional iterations, we apply the recent local mesh updating method of Legentil et al. (2022).

To illustrate the benefits of the methodology, we propose a synthetic case study on a layer-cake model (Section 3). Additionally, we present different ways to exploit the wealth of information provided by the transdimensional posterior distribution. Uncertainties are estimated not only for permeabilities at each depth but also for layer locations. Furthermore, we propose to quantify the level of information contributed by each dataset along depth (i.e., information gain) using a Kullback-Leibler divergence between the prior and posterior distributions.

This application allows us to evaluate the practical convergence of the method and to check the parsimonious nature of the RJMCMC sampler (Section 4), which converges towards the minimum number of parameters that explain the data, accounting for uncertainties (MacKay, 2003; Gallagher et al., 2009). Section 4 discusses the implications of the results for reservoir modeling, and their extension to more realistic 2D and 3D cases.

# 1. Transdimensional inversion for a generic layered model

## 1.1. Bayesian inference

The aim of Bayesian inference is to quantify the a posteriori probability distribution, or posterior distribution, which is the probability density function (PDF) of the model parameters given the observed data (Smith, 1991). This formalism relies on Bayes' theorem (Bayes, 1763), that poses the data integration problem and formalizes the use of prior knowledge to give the posterior distribution $P(m|d_{obs})$ as

$$P(m|d_{\text{obs}}) = \frac{P(m) \cdot P(d_{\text{obs}}|m)}{P(d_{\text{obs}})}, \qquad (1)$$

where $a|b$ means a given (or conditional on) b and is the probability of event a when the event b is known, and:

* $P(m)$ is the a priori, or prior probability density of the model m, that reflects what we know about the parameters of interest of the model m before measuring the data $d_{obs}$. A Bayesian formulation allows for the integration of any prior knowledge of the model, provided that this information can be translated into a probability distribution (Scales and Snieder, 1997). Thus, all inferences from the data are made relative to this prior distribution.

* $P(d_{obs}|m)$ is the likelihood function, a conditional PDF which is the probability of observing the measured data given a particular model. This function quantifies the agreement between simulated and observed data. Simulated data d are synthetic measurements predicted from the model parameters as $d = g(m)$. Observed data $d_{obs}$ are contaminated with errors $\epsilon$, related to model predictions as $d_{obs} = g(m) + \epsilon$, where $\epsilon$ includes both measurement and model errors (Tarantola, 2005). The predictions are typically compared to the observed data through a misfit function $\phi$ measuring the square of the Mahalanobis distance (Mahalanobis, 2018) to take into account correlations between observed and predicted data as

$$\phi(m) = (g(m) - d_{\text{obs}})^T C_d^{-1} (g(m) - d_{\text{obs}}), \qquad (2)$$

where $C_d$ is the covariance matrix of modeling and measurement errors.
<figure><img src="image_2.png" /><figcaption>Fig. 1. Layered model parameterization used in this paper. The geometry is defined in 1D by the combination of piecewise constant basis functions with a set of interface depth values delimiting geological layers of variable thickness. Each layer owns a constant value of permeability.</figcaption></figure>

A common choice is to define the likelihood distribution as a Gaussian function:

$$P(\mathbf{d}_{\text{obs}}|\mathbf{m}) = \frac{1}{\sqrt{|C_d|(2\pi)^n}} \exp\left\{-\frac{1}{2}\phi(\mathbf{m})\right\}, \quad (3)$$

where n is the number of data points.

This Gaussian assumption is justified by the composition of the noise from diverse error sources, leading to an asymptotically Gaussian distribution, a consequence of the central limit theorem, assuming finite variance errors.

* $P(\mathbf{d}_{\text{obs}})$ is a normalization factor integrated over the model space. This term is called the total evidence and is very expensive to compute. MCMC methods circumvent this computation, so that Eq. (1) is often written without evidence as a proportionality relationship:

$$P(\mathbf{m}|\mathbf{d}_{\text{obs}}) \propto P(\mathbf{m}) \cdot P(\mathbf{d}_{\text{obs}}|\mathbf{m}). \quad (4)$$

Typically, this posterior probability cannot be calculated analytically, leading to approximation via sampled models. Sampling schemes aim to generate an ensemble of models, whose density in parameter space approximates the posterior probability distribution (Mosegaard and Sambridge, 2002). This posterior distribution represents the system after prior knowledge has been updated by data. It is common to derive the optimal or average model or to establish marginal probability distributions for model parameters from the available model set. Spatially, model parameter uncertainties can also be visualized by an error map computed from the samples.

Bayesian inference is extensible to transdimensional settings in which the number of parameters is poorly known. In the case of a layered model that incorporates an unknown number of geological layers $n_L$, Bayes' theorem is rewritten as (Sambridge et al., 2006)

$$P(\mathbf{m}, n_L | \mathbf{d}_{\text{obs}}) \propto P(\mathbf{m}, n_L) \cdot P(\mathbf{d}_{\text{obs}} | \mathbf{m}, n_L) \\ \propto P(\mathbf{m} | n_L) \cdot P(n_L) \cdot P(\mathbf{d}_{\text{obs}} | \mathbf{m}, n_L), \quad (5)$$

where the number of layers $n_L$ is a variable parameter describing model dimension and following its own prior discrete PDF $P(n_L)$.

To estimate the likelihood and prior probabilities, and get realizations of the posterior probability density related to the Bayesian formalism expressed above, the definition of three components is imperative: The model vector $\mathbf{m}$, the observed field data $\mathbf{d}_{\text{obs}}$, and the model predictions $\mathbf{g}(\mathbf{m})$.

## 1.2. Model parameterization

We use the same formalism as Sambridge et al. (2013) to introduce the transdimensional approach with the proposed parameterization $\mathbf{m}$ of the problem. We consider a one-dimensional conformal model composed of $n_L$ geological horizontal layers delimited by a set of $n_L - 1$ interface locations. The depths of these interfaces $z_i$ are defined between $z_{min}$ and $z_{max}$, the fixed depth of the upper and lower horizons, respectively. This model defines the permeability field $k$ by

$$k(z) = \sum_{i=1}^{n_L} k_i \psi_i(z), \quad (6)$$

corresponding to the physical quantity of interest to be constrained by the data at all possible depth locations $z$ and for which the dimension (i.e., the total number of parameters) is of unknown size. The model is expanded in terms of a variable number of piecewise constant basis functions $\psi_i(z)$ representing the geometry of the model: The $i$th basis function $\psi_i$ is unity in the $i$th layer and zero elsewhere, and is defined with an indicator function $1$ as

$$\psi_i(z) = \begin{cases} 1 & \text{if } z \in [z_i, z_{i+1}) \\ 0 & \text{otherwise.} \end{cases} \quad (7)$$

In Eq. (6), $k_i$ is the coefficient of the $i$th basis function of the targeted petrophysical property, i.e., the effective permeability of the layer $i$. For simplicity, we only consider isotropic permeability values.

Fig. 1 displays the key elements defining the model geometry. Each geological layer is characterized by a constant permeability value $k_i$, a thickness value $h_i$ and bounded by interfaces of depth values $z_i$ and $z_{i+1}$ for a 1D depth coordinate system. The model $\mathbf{m}$ is, therefore, defined by a parameter vector of size $2n_L - 1$ that is $\mathbf{m} = \{\mathbf{z}, \mathbf{k}\}$ with $\mathbf{z} = [z_2, z_3, \dots, z_{n_L}]^T$, $\mathbf{k} = [k_1, k_2, \dots, k_{n_L}]^T$. Note that the values in $\mathbf{z}$ are sorted by increasing depth to ensure the consistency of indices in Eqs. (6) and (7).

## 1.3. Prior probability distribution

Considering the interface random vector $\mathbf{z}$ as independent of the permeability random vector $\mathbf{k}$ conditionally to the number of layers $n_L$, $P(\mathbf{m}|n_L)$ equates to $P(\mathbf{z}|n_L) \cdot P(\mathbf{k}|n_L)$. Hence, the transdimensional prior probability writes

$$\begin{aligned} P(\mathbf{m}, n_L) &= P(\mathbf{m}|n_L) \cdot P(n_L) \\ &= P(\mathbf{z}|n_L) \cdot P(\mathbf{k}|n_L) \cdot P(n_L). \end{aligned} \quad (8)$$

To describe all possible positions of interfaces in the system, we adopt the “grid trick” approach from Bodin and Sambridge (2009) by discretizing the field with an underlying finite 1D grid of $n_z$ cells for mathematical convenience, so that no more than one horizon lies in a cell, hence the number of horizons is finite. This corresponds to a minimal layer thickness value $h_{min} = (z_{max} - z_{min})/n_z$. Then, the uniform probability of horizon depths for $n_L$ layers (or $n_L - 1$ interfaces) is

$$P(\mathbf{z}|n_L) = \left[ \frac{n_z!}{n_L!(n_z - n_L + 1)!} \right]^{-1}. \quad (9)$$

Regarding the prior distribution for permeabilities, it comprises a set of independent constant values in each layer defined by an uninformative uniform distribution as

$$P(\mathbf{k}|n_L) = \begin{cases} \prod_{i=1}^{n_L} \frac{1}{\Delta_k} & \text{if } k_i \in K \\ 0 & \text{otherwise,} \end{cases} \quad (10)$$

where $\Delta_k = (k_{max} - k_{min})$ and $K = \{k_i \in \mathbb{R} | k_{min} \le k_i \le k_{max}\}$.

The prior discrete probability of the number of model layers $P(n_L)$ is distributed uniformly as

$$P(n_L) = \begin{cases} \frac{1}{n_L} & \text{if } n_L \in I \\ 0 & \text{otherwise,} \end{cases} \quad (11)$$
where $\Delta_{n_L} = (n_{Lmax} - n_{Lmin})$ and $I = \{n_L \in \mathbb{N} | n_{Lmin} \le n_L \le n_{Lmax}\}$.

Combining Eqs. (8), (9), (10), and (11), the prior probability becomes

$$P(\mathbf{m}, n_L) = \begin{cases} \frac{n_L!(n_z-n_L+1)!}{n_z!(\Delta_k)^{n_L}(\Delta_{n_L})} & \text{if } n_L \in I, k_i \in K \\ 0 & \text{otherwise.} \end{cases} \quad (12)$$

As shown in Eq. (4), combining $P(\mathbf{m})$ with the information provided by the data in the likelihood $P(\mathbf{d}_{\text{obs}}|\mathbf{m})$ allows the evaluation of the posterior distribution for any given model $\mathbf{m}$.

## 1.4. Data

Our second focus lies with the observed data $\mathbf{d}_{\text{obs}}$. We use two distinct datasets in a complementary way during inversion: Well log data $\mathbf{k}_{\text{obs}}$, and flow data $\mathbf{f}_{\text{obs}}$ obtained during reservoir exploitation as $\mathbf{d}_{\text{obs}} = \{\mathbf{k}_{\text{obs}}, \mathbf{f}_{\text{obs}}\}$. In this paper, we use synthetic data to explore the properties and check the results of the proposed methodology.

### 1.4.1. Well logs of permeability

The first dataset is composed of a collection of points making up one or more vertical well logs. These static reservoir data provide $n_k$ well points of localized and regularly sampled rock permeability values $\mathbf{k}_{\text{obs}} = [k_{\text{obs}1}, k_{\text{obs}2}, \dots, k_{\text{obs}n_k}]^T$ along the vertical z-axis.

The influence of well log data on parameter inference can be controlled by the noise model. Following Bodin (2010) and geostatistical literature (e.g., Chiles and Delfiner, 2012), we propose a correlated noise model, as correlated data errors provide less information than independent errors. Indeed, independent data provide unique and separate pieces of information about the system, while correlated data are partly redundant, thus do not add as much new knowledge (Tarantola, 2005). This configuration of a noise describing both independent measurement errors and intrinsic spatial correlation of the rocks is suitable since spatial correlations exist in reality (Chiles and Delfiner, 2012). However, as layer interfaces correspond to permeability jumps in our model, we consider spatial correlation only within the same layer.

To introduce realistic uncertainty into permeability estimates, we apply a random noise following a log-normal distribution to the permeability data $\mathbf{k}_{\text{obs}}$. This approach involves initially converting the permeability to a logarithmic scale. A Gaussian noise is then added to data according to the distribution $N(\log(k_i), \sigma_{\log(k_{\text{obs}})})$ and a correlation factor $r_{\mathbf{k}_{\text{obs}}}$, where $k_i$ is the permeability within the current layer $i$ and $\sigma_{\log(k_{\text{obs}})}$ is the (assumed constant) standard deviation of the Gaussian noise of $\log(k_{\text{obs}})$.

The synthetic well log data are sampled from a reference layer-wise constant permeability model by randomly generating a vector in each layer. For example, in a five-layer reference model, five independent datasets containing spatially correlated well log samples are produced per well. For a more complete explanation of the construction of synthetic well data, refer to Appendix A.

### 1.4.2. Flow data

In the context of history matching for reservoir uncertainty quantification, flow data, also known as production data, encompass time-varying responses that are non-linearly related to static reservoir attributes such as facies, porosity, and permeability. Notable among these dynamic data are water-cut $wc$, a ratio of water produced compared to the volume of total liquid produced, and well bottomhole pressure $P$, which reflects the reservoir ability to push fluids to the surface (Dake, 2001).

Two observed flow data vectors are considered in this work: $\mathbf{f}_{\text{obs}} = \{wc_{\text{obs}}, P_{\text{obs}}\}$ with $wc_{\text{obs}} = [wc_{\text{obs}1}, wc_{\text{obs}2}, \dots, wc_{\text{obs}n_f}]^T$ and $P_{\text{obs}} = [P_{\text{obs}1}, P_{\text{obs}2}, \dots, P_{\text{obs}n_f}]^T$ for $n_f$ well measurements for each property. We assume an independent and identically distributed noise for both properties. To effectively manage the constraints of the water-cut property, which is a ratio defined between 0 and 1 with values frequently

at zero, we employ a logarithmic transformation defined as $\log(wc/(1-wc))$. The water-cut data are thus represented in the transformed scale to add a truncated Gaussian noise following $N(\log(wc_i/(1-wc_i))), \sigma_{\log(wc_{\text{obs}})}$ if $wc_i > 0$ and set to 0 if $wc_i = 0$, where $wc_i$ is the water-cut measurement at time step $i$ and $\sigma_{\log(wc_{\text{obs}})}$ is the standard deviation of the independent Gaussian noise applied to the log-transformed data. After applying the noise, the water-cut data are transformed back to their original physical space for visualization and interpretation. The pressure data noise is introduced from $N(P_i, \sigma_{P_{\text{obs}}})$, where $P_i$ is the well pressure measurement at time step $i$ and $\sigma_{P_{\text{obs}}}$ is the independent Gaussian noise level.

The synthetic flow data are sampled using the reference model mentioned in the previous part, with the associated forward problem related to fluid transfers (Section 1.5.2).

### 1.5. Forward problem

A total number of $n = n_k + 2n_f$ data points are available since the history matching properties are characterized by two vectors of size $n_f$. We relate the model vector $\mathbf{m} = \{\mathbf{z}, \mathbf{k}\}$ with $n$ observable data $\mathbf{d} = \{\mathbf{k}, \mathbf{wc}, \mathbf{P}\}$ by $\mathbf{d} = \mathbf{g}(\mathbf{m}) = [\mathbf{g}_k(\mathbf{m}), \mathbf{g}_{wc}(\mathbf{m}), \mathbf{g}_P(\mathbf{m})]^T$ denoting predicted data for all data types.

#### 1.5.1. Set of layer permeability values

The measured permeability well logs $\mathbf{k}_{\text{obs}}$ can be used to constrain the model permeability property in the medium. There are numerous geostatistical estimation methods or simulation methods to extrapolate well log data into a discretized model (Pyrcz and Deutsch, 2014), but the transfer of well data at model resolution is a source of approximations. Therefore, this paper considers $\mathbf{k}_{\text{obs}}$ as data to invert in a Bayesian sense, jointly with flow data. Predictions correspond to the vector of model layer permeabilities $\mathbf{k}_{\text{pred}} = \mathbf{g}_k(\mathbf{m}) = \Psi \cdot \mathbf{k}$, where $\Psi$ is the $n_L \cdot n_k$ matrix of coefficients $\psi_{ij} = \psi_i(z_j)$ (Eq. (7)) and $z_j$ is the depth of the sample $k_j$. Hence, the forward problem simply consists of collecting the permeability values of the layers from the geological model. Fig. 2 is an example of a reservoir model with a set of constant permeability values in the layers. In this representation, the data associated with the model come from a water injector well and an oil producer well, which are used for reservoir exploitation and history matching, and a permeability log associated with the injector well.

#### 1.5.2. Flow simulation and remeshing

The relationship between flow data and reservoir permeability is a function of many factors, including fluid properties, well control (rate, pressure, perforations), and reservoir geometry and properties (porosity, facies) (Oliver and Chen, 2011). The reservoir model requires the inclusion of state variables, such as water saturation to reflect the fraction of pore space occupied by water. The temporal evolution of these variables follows Darcy's law and mass conservation principles, and is inherently dependent on permeability (Oliver et al., 2008).

The forward step solving porous flow problem equations to compute history matching property predictions $\{wc, P\} = [\mathbf{g}_{wc}(\mathbf{m}), \mathbf{g}_P(\mathbf{m})]^T = ([wc_1, wc_2, \dots, wc_{n_f}]^T, [P_1, P_2, \dots, P_{n_f}]^T)$ is commonly called flow simulation. The evolution through time of reservoir variables can be numerically computed, in accordance with boundary and well conditions. Although the model parameterization is 1D (Section 1.2), the flow simulation is performed within a 2D (or 3D) field to capture the spatial heterogeneities and complexities of a reservoir.

In this paper, we run a two-dimensional incompressible two-phase flow simulation, based on conservation law and Darcy's law, discretized with CVFEM, using the RINGFlow code (Anquez et al., 2020), further detailed in Appendix B. The water injector and oil producer are considered as punctual elements, each located at a single point within the reservoir. The water injector is positioned at the lower-left corner, and the oil producer is at the upper-right corner. Both are set as boundary conditions in the simulator.
<figure><img src="image_3.png" /><figcaption>Fig. 2. The model **m** is constrained by both static data k<sub>obs</sub> and dynamic data f<sub>obs</sub> from field wells. The distance between well log measurements and model permeability values is measured for each layer. The layer geometry is discretized by a conformal triangular mesh for 2D flow simulation providing production data.</figcaption></figure>

The spatial discretization of the reservoir models is done with a
triangular mesh conformable to the horizons (Fig. 2). To efficiently
maintain the mesh conformity throughout geometric model perturba-
tions performed during inversion, we use the local updating strategy
for triangular meshes of Legentil et al. (2022).

1.6. Combining datasets using a joint likelihood probability

Datasets of different nature can be integrated in a joint inversion.
Capitalizing on the inherent property of Bayesian methods, the joint
inversion is particularly suitable for a robust and holistic interpretation
of complex systems (Vozoff and Jupp, 1975).

Assuming all data types are independent, the likelihood of the
observed field data **d**<sub>obs</sub> = {**k**<sub>obs</sub>, **w**<sub>obs</sub>, **P**<sub>obs</sub>} writes

$$
P(\mathbf{d}_{\text{obs}} | \mathbf{m}) = P(\mathbf{k}_{\text{obs}} | \mathbf{m}) \cdot P(\mathbf{w}_{\text{obs}}, \mathbf{P}_{\text{obs}} | \mathbf{m}). \quad (13)
$$

In this Bayesian framework, the noise levels and data covariance naturally act as weighting factors in the joint likelihood: A dataset with a lower noise level is given more weight and vice versa. These elements dictate the confidence assigned to each dataset, balancing their integration based on reliability and internal consistency (Bodin et al., 2012b). To evaluate the joint probability in Eq. (13), we must specify likelihood functions for each respective dataset.

1.6.1. Well log data likelihood function

Because the well log data noise is correlated, the data fit $\phi_k$ is defined as the Mahalanobis distance (Eq. (2)) between layer permeability values $g_k(\mathbf{m})$ and vectors of well log data $\mathbf{k}_{\text{obs}}$. Thus, for a given model layer, we sum the distance between each well log point that belongs to this layer and the current layer permeability. The process is repeated over all model layers, and this piecewise constant regression provides a single fit-to-data value for the whole model:

$$
\phi_k(\mathbf{m}) = \phi_k(\mathbf{z}, \mathbf{k}) = ((\mathbf{\Psi} \cdot \mathbf{k}) - \mathbf{k}_{\text{obs}})^T \mathbf{C}_{\mathbf{k}}^{-1} ((\mathbf{\Psi} \cdot \mathbf{k}) - \mathbf{k}_{\text{obs}}). \quad (14)
$$

The complete description of the piecewise correlated noise model,
leading to the calculation of the inverse covariance matrix, of the misfit
function, and of the likelihood probability is presented in Appendix C.
This parameterization requires two factors: The estimated noise stan-
dard deviation σ<sub>k</sub> and a correlation factor r<sub>k</sub>. For simplicity, we assign
a constant value to these factors, identical to the true field parameters
σ<sub>log(k<sub>obs</sub>)</sub> and r<sub>k<sub>obs</sub></sub> associated with observed well log data. From Eq. (3),
the resulting expression for the likelihood PDF is

$$
P(\mathbf{k}_{\text{obs}} | \mathbf{m}) = \frac{1}{\sqrt{|C_{\mathbf{k}}| (2\pi)^{\eta_{\mathbf{k}}}}} \exp \left\{ -\frac{1}{2} \phi_{\mathbf{k}}(\mathbf{m}) \right\}. \quad (15)
$$

1.6.2. Flow data likelihood function

Contrary to well log data, flow data errors are assumed independent,
with constant noise standard deviation σ<sub>wc</sub> and σ<sub>p</sub>. Therefore, the
covariance matrix is diagonal and the misfit functions read

$$
\phi_{wc}(\mathbf{m}) = \left\| \frac{g_{wc}(\mathbf{m}) - w\mathbf{c}_{obs}}{\sigma_{wc}} \right\|^2, \quad \phi_{p}(\mathbf{m}) = \left\| \frac{g_{p}(\mathbf{m}) - \mathbf{P}_{obs}}{\sigma_{p}} \right\|^2. \qquad (16)
$$

The final evaluation of likelihood probability for reservoir produc-
tion curves is simply a multiplication of both likelihood expressions
from Eqs. (3) and (16):

$$
P(\mathbf{w}_{\text{obs}}, \mathbf{P}_{\text{obs}} | \mathbf{m}) = \frac{1}{(2\pi)^{\eta_{\text{f}}}} (\sigma_{\text{wc}} \sigma_{\text{P}})^{\eta_{\text{f}}} \exp \left\{ -\frac{1}{2} [\phi_{\text{wc}}(\mathbf{m}) + \phi_{\text{P}}(\mathbf{m})] \right\}. \quad (17)
$$

2. Implemented reversible jump algorithm

2.1. Theoretical framework

Two significant challenges arise when characterizing the full pos-
terior density P(m|d<sub>obs</sub>). First, the posterior density function cannot
be expressed in a simple analytical form, only allowing its estimation
through evaluations at varied positions in the model space. Second, as
the dimension of the model space increases, the number of models to
test becomes huge, making the complete characterization of the poste-
rior or uniform random sampling impractical: This is commonly called
the “curse of dimensionality” (see Fernández-Martínez and Fernández-
Muñiz, 2020, and references therein). For history matching problems,
these difficulties are paramount as the model space to explore is very
large, due to the large number of unknown parameters related to the
small amount of data.

Markov chain Monte Carlo sampling, a flexible and popular tech-
nique, helps (partly) overcome these problems. This iterative stochas-
tic method generates samples from the posterior distribution and is
applicable even when the distribution is non-Gaussian (e.g., Doucet
and Wang, 2005). MCMC methods are frequently applied to fixed-
dimensional geophysical inverse problems (Sambridge and Mosegaard,
2002).

For inference on both model parameters and dimensionality (i.e.,
when the number of parameters is variable), the reversible jump MCMC
sampler pioneered by Geyer and Møller (1994) and Green (1995)
presents a generalization of the MCMC technique. This methodol-
ogy enables sampling or “jumping” across spaces of varying dimen-
sions, extending the settings of the well-known Metropolis-Hastings
algorithm (Metropolis et al., 1953; Hastings, 1970).

This method builds a sequence of sampled models in a Markov chain
of order 1, which means that each generated model is a perturbation of
the previous sample. The chain starts at a random point, drawn from
the prior distribution, and a proposal probability distribution q(m'|m)
is used to generate a new random model m' from the current model
m. The prime used in mathematical expressions denotes the model
state after being randomly perturbed. The vector lengths of the current
model m and the proposed model m' may not necessarily be identical in
a transdimensional problem, as the perturbed model could incorporate
a different number of parameters.

In all MCMC methods, models are accepted or rejected based on
the acceptance probability α(m'|m). This probability governs the model
transition, with the chain either moving to the candidate $\mathbf{m}'$ if $\alpha(\mathbf{m}'|\mathbf{m})$ is higher than a randomly drawn number $u \sim U(0, 1)$ or staying with the current model $\mathbf{m}$ if not. The acceptance probability is given by

$$ 
\begin{aligned}
\alpha(\mathbf{m}'|\mathbf{m}) &= \min[1, \text{prior ratio} * \text{likelihood ratio} * \text{proposal ratio} * |J|] \\
&= \min \left[ 1, \frac{P(\mathbf{m}')}{P(\mathbf{m})} \cdot \frac{P(\mathbf{d}_{\text{obs}}|\mathbf{m}')}{P(\mathbf{d}_{\text{obs}}|\mathbf{m})} \cdot \frac{q(\mathbf{m}|\mathbf{m}')}{q(\mathbf{m}'|\mathbf{m})} \cdot |J| \right] \\
&= \min \left[ 1, \frac{P(\mathbf{m}'|\mathbf{d}_{\text{obs}})}{P(\mathbf{m}|\mathbf{d}_{\text{obs}})} \cdot \frac{q(\mathbf{m}|\mathbf{m}')}{q(\mathbf{m}'|\mathbf{m})} \cdot |J| \right],
\end{aligned} 
\quad (18) $$

where $J$ is the Jacobian matrix accounting for model dimension changes. In our parameterization, which adheres to fully discretized parameter spaces with non-continuous probability distributions, $|J| = 1$, indicating a volume-preserving transformation (Denison et al., 2002). After an initial burn-in period, corresponding to the number of iterations before which the random walk is assumed to have reached stationarity, the chain asymptotically converges to the posterior distribution $P(\mathbf{m}|\mathbf{d}_{\text{obs}})$ (Tierney, 1994; Green, 2003).

When the dimensionality of $\mathbf{m}$ and $\mathbf{m}'$ remains constant, the proposal distribution $q(\mathbf{m}'|\mathbf{m})$ can be achieved by modifying the $i$th component of $\mathbf{m}$ with a random variable $\varepsilon$ sampled from a normal distribution $N(0, \sigma_{\varepsilon})$ as

$$ \mathbf{m}' = \mathbf{m} + \varepsilon. \quad (19) $$

Transition probabilities for a jump between dimensions, expressed by the perturbation $q(\mathbf{m}'|\mathbf{m})$, are a more delicate aspect of transdimensional methods. The design of these transitions significantly influences the convergence rate of the algorithm but not the inversion outcomes. Hence, while proposals can be flexibly defined, a suitable mathematical definition is necessary for efficient algorithm dynamics.

## 2.2. Proposal distribution for the layered reservoir case

We consider four possible types of perturbations of the current model $\mathbf{m}$ at each step of the Markov chain. These include a Gaussian perturbation of the permeability parameter, and geometric modifications comprising moves, births, and deaths of layer interfaces (Fig. 3(a)). Combining all these perturbations, it is possible to generate a broad range of configurations and permeability models with variable degrees of freedom to explore the model space. The selection and application of the perturbation type from $\mathbf{m}$ to $\mathbf{m}'$ follow this scheme:

1. At each even chain step, a layer permeability value is randomly selected with probability $1/n_L$ and changed (Fig. 3(b)) according to a Gaussian distribution with a standard deviation $\sigma_{\varepsilon_k}$ as $k'_i = k_i + \varepsilon_k$, $\varepsilon_k \sim N(0, \sigma_{\varepsilon_k})$. If the new permeability value falls outside the predefined prior range, the proposed model is rejected and the current model is counted twice in the ensemble solution.

2. At each odd chain step, the perturbation concerns the model geometry, with each type of perturbation being randomly selected and having an equal probability of 1/3:

*   **Move:** Randomly changes the positions of layer interfaces (Fig. 3(c)). One interface $i$ is selected among the existing set with probability $1/(n_L - 1)$ and its position is moved with a standard deviation $\sigma_{\varepsilon_z}$ by $z'_i = z_i + \varepsilon_z$, $\varepsilon_z \sim N(0, \sigma_{\varepsilon_z})$. This perturbation allows for potential layer swaps as interfaces are sorted by increasing depth after the perturbation. If the new position is outside boundaries or if some layer thickness $h_i = z_{i+1} - z_i$ is less than the minimal allowed layer thickness $h_{\min}$, the proposed model is rejected, while the current model is retained again.

*   **Birth:** Creates a new layer in the model (Fig. 3(d)). A location along the model $z$-axis that is not already occupied is selected from the uniform prior distribution $U(z_{\min}, z_{\max})$ to create an interface. This results in the formation of a new layer ($n_L \to n_L + 1$) with a thickness value that

<figure><img src="image_4.png" /><figcaption>Fig. 3. (a) Current model at the $i$th iteration of RJMCMC. Four different model perturbations are allowed to generate a candidate model. (b) Change of permeability in a randomly selected layer. (c) Move of a randomly selected layer interface. (d) Birth of a new interface and sampling of the permeability value in the new layer. (e) Death of a randomly selected interface and merge of permeability between neighboring layers.</figcaption></figure>

must be greater than $h_{\min}$ to prevent preliminary rejection. A random choice determines if the new layer is located above or below the inserted interface. The new permeability value $k'_{n_L+1}$ inside the layer affected by a birth is computed from a Gaussian proposal distribution based on the permeability value of the previous layer at birth location $k_i$ and a standard deviation $\sigma_{\varepsilon_b}$, that is $k'_{n_L+1} = k_i + \varepsilon_b$, $\varepsilon_b \sim N(0, \sigma_{\varepsilon_b})$.

*   **Death:** Deletes an existing layer (Fig. 3(e)). The death step randomly selects an already existing interface from the current set with probability $1/(n_L - 1)$ and deletes it to cause the removal of one model layer ($n_L \to n_L - 1$). This perturbation acts as a reversal of the birth process, merging the deleted layer with the permeability value of the neighboring layer selected randomly from above or below the removed interface, to set the post-death permeability value $k'_j$.

All the proposal standard deviation terms defined above characterize the perturbation intensity that must be small enough to accept a sufficient number of proposed models but high enough to ensure effective algorithm convergence (Brooks et al., 2003). In this implementation, the proposal standard deviation values are automatically recalculated with a correction factor every 100 iterations if the acceptance rate for a perturbation type is too high (> 30%) or too low (< 10%) to have an adaptive proposal that prevents the algorithm from becoming computationally prohibitive. The selection process for a specific perturbation could also be optimized to enhance convergence. For instance, permeability perturbations do not require geometric changes, thus eliminating the need for local model remeshing when evaluating the flow data likelihood. Consequently, their evaluation is less expensive compared to other perturbations and could be performed more frequently.
<figure><img src="image_5.png" /><figcaption>Fig. 4. Flowchart of the transdimensional cascaded Metropolis algorithm.</figcaption></figure>

To evaluate the acceptance probability (Eq. (18)), we still need to compute the probability ratios (proposal ratio, prior ratio and like- lihood ratio) between the proposed model and the current model. The acceptance criteria for the different perturbations are the same as described by Bodin and Sambridge (2009), see Appendix D.

2.3. Cascaded Metropolis

The Bayesian framework allows for the incorporation of various
data types into the model likelihood probability. Although the flow
equations are connected to the permeability field through a non-linear
relationship and require a high computational cost for the numerical
simulation to be evaluated, the problem related to the well log data
is a computationally inexpensive piecewise regression. It is therefore
convenient to first evaluate the likelihood probability of well data given
the considered model P(kobs|m) and calculate the porous flow forward
problem to assess P(wcobs, Pobs|m) only when necessary.

This can be achieved by considering the available data sequentially
with the cascaded Metropolis, as described by Tarantola (2005). This
algorithm is an extension of the classic Metropolis sampler that subdi-
vides the acceptance–rejection criterion for the candidate model into as
many criteria as needed. Here, we use an acceptance criterion for each
dataset: α<sub>well_log</sub>(m'|m) for the well logs of permeability and α<sub>flow</sub>(m'|m)
for the history matching data. The complete process for the resulting
cascaded Metropolis sampler is illustrated in Fig. 4 and can be described
as follows:

0. At the starting point of the chain, generate the initial model **m**,
evaluate its likelihood, and set it as the current model.

1. At iteration *i*, perturb the current model **m** to propose a candidate **m'** (Section 2.2).

2. Solve the forward problem related to well log data (Section 1.5.1).

3. Determine $P(k_{\text{obs}}|\mathbf{m}')$ (Section 1.6.1).

4. Draw a number $u \sim U(0, 1)$ and evaluate $\alpha_{\text{well\_log}}(\mathbf{m}'|\mathbf{m})$. If $\mathbf{m}'$ is accepted ($u \le \alpha$), proceed to step 5; otherwise, return to step 1.

5. Solve the forward problem related to flow data: Remesh the model if the perturbation is geometric, and perform the flow simulation (Section 1.5.2).

6. Determine $P(wc_{\text{obs}}, \mathbf{P}_{\text{obs}}|\mathbf{m}')$ (Section 1.6.2).

7. Draw a number $u \sim U(0, 1)$ and evaluate $\alpha_{\text{flow}}(\mathbf{m}'|\mathbf{m})$. If $\mathbf{m}'$ is accepted ($u \le \alpha$), replace $\mathbf{m}$ with $\mathbf{m}'$ and, in all cases, return to step 1.

This cascaded algorithm allows for significant computational time
savings without affecting the convergence of the chain towards the
posterior distribution.

2.4. Parallel computing and temperature

When introducing flow data into the problem, the non-linearity
and non-uniqueness make convergence more difficult to reach due
to the existence of multiple local minima. To overcome convergence
issues during the joint inversion, parallel Markov chains can be used
to increase the number of realizations (Rosenthal, 2000). Starting from
different initial states, these chains sample the model space indepen-
dently, and post burn-in results are merged after inversion in only one
collection of model space samples.
<figure><img src="image_6.png" /><figcaption>Fig. 5. (a) Simulated reference model (black) that aims to be recovered by the transdimensional method, and noisy well log data k<sub>obs</sub> (gray). (b) Water saturation map obtained from the flow simulation on the reference model at half the injected pore volume. (c) Water-cut data wc<sub>obs</sub> generated with independent (truncated at 0) Gaussian noise applied on reference model values. (d) Well pressure data P<sub>obs</sub> generated with independent Gaussian noise applied on reference model values.</figcaption></figure>

To further accelerate convergence, a temperature *T* is also introduced in the likelihood distribution of (normally distributed and independent) flow data (Eq. (16)). In this technique inspired by simulated annealing (Kirkpatrick et al., 1983), the system starts with a high temperature. This initial state relaxes the constraints imposed by the data term and encourages exploration by artificially increasing the data variance, yielding more variability in the models that are accepted. *T* gradually cools down to converge towards the desired solution, acting as a gradually decreasing coefficient during the burn-in phase:

$$ \phi_{wc}(\mathbf{m}) = \left\| \frac{g_{wc}(\mathbf{m}) - \mathbf{w}c_{obs}}{T \cdot \sigma_{wc}} \right\|^2, \quad \phi_P(\mathbf{m}) = \left\| \frac{g_P(\mathbf{m}) - \mathbf{P}_{obs}}{T \cdot \sigma_P} \right\|^2, \quad (20) $$

where $T_i = f T_{i-1}$ is the temperature for the current iteration *i* and *f* is the magnitude of the decrease while $T > 1$.

## 3. Examples

### 3.1. Constructing a synthetic dataset

We evaluate the properties of the joint inversion method using synthetic datasets of well logs and production curves for history matching. To create these data, a reference model geometry is drawn from a uniform distribution for *n*<sub>L</sub> and *z* values.

For our first tests of RJMCMC convergence, we build a simple 2D geological model with unitless spatial boundaries for *x* and *z* between 0 and 1, and a minimum thickness value $h_{min} = 0.03$ to avoid micro-layers

(which would raise remeshing challenges and may be unrealistic). Furthermore, the model adopts a layer-cake structure, which simplifies the sequence of interfaces into straight lines expressed by a single *z* coordinate. This makes the problem one-dimensional, and the *x* coordinate only plays a role in the forward flow problem.

The reference model of the reservoir field is represented in Fig. 5(a) (black) with a one-dimensional depth section. It contains seven non-uniform layers with a thick layer at the bottom, cut by a very thin layer, and overlaid by thinner layers at the top. The associated isotropic permeability field of this reference model is randomly computed from a uniform distribution and represented on a logarithmic scale log(*k*) to compare different orders of magnitude within ranges that yield realistic permeability contrasts, from 0.1 mD (10<sup>-16</sup> m<sup>2</sup>) for the thinnest layer (e.g., shale) to 1000 mD (10<sup>-12</sup> m<sup>2</sup>) for the reservoir top (e.g., clean, well-sorted sandstone).

For simplicity, a single well is simulated from the reference permeability field with $n_k = 50$ regularly spaced sample points k<sub>obs</sub>, each having a noise standard deviation $\sigma_{k_{obs}} = 3$ mD which becomes $\sigma_{log(k_{obs})} = \log(3) = 1.0986$ in log(*k*). The constant correlation coefficient between two samples is $r_{k_{obs}} = 0.85$, corresponding to an exponential variogram range of 6.2 (see Appendix C). The geological model can be seen as an effective representation of the permeability of the medium while the well log measures the spatial correlations of this permeability property inside layers. For instance, the thick layer located between 0.5 and 0.9 could represent an alternation of sandstone and sand-clay strata
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** (a) Results of the RJMCMC inversion of well log data, on 30,000 post burn-in iterations through a depth section. The initial state of the transdimensional Markov chain (blue), the reference model (black), the noisy well log data (gray), and the average model from all accepted realizations of the model space (red) ± one standard deviation (dashed green) are plotted for the log-permeability property. (b) Posterior probability density to have an interface at each depth represented in 1D. (c) Evolution of convergence properties (normalized well log misfit function, MAE of log-permeability between the current model and the reference model, and global acceptance rate) during post burn-in inversion. (d) Evolution of the number of layers and histogram of the posterior ensemble $P(n_L|d_{obs})$ during post burn-in inversion.</figcaption></figure>

as shown in the well log representation of Fig. 5(a) (gray), depicted by a single constant value of log-permeability in the reference model. Conversely, the very thin layer located towards the reservoir base contains only one point, and should therefore not be easily detectable from well log data only.

A flow simulation is performed to generate synthetic observed production data. The domain is initially filled with oil ($S_o = 1$) and discretized in 2D by an initial coarse synthetic mesh of 184 triangles (this number may vary during transdimensional iterations based on remeshing outcomes) to optimize computation time, as the discrete system in flow simulation must be solved at each triangle node. We inject water at the bottom left corner and produce oil and water at the opposite corner. The porosity is initialized as a constant across the field and is set to 0.2. The water is injected with a rate of 0.1 pore volume (PV = porosity * total volume of the rock) per day up to 1 PV. Fig. 5(b) shows the water saturation property at 0.5 PV of water injection. The simulation comprises $n_f = 80$ time steps for each property $w_c$ and $P_{obs}$, represented respectively in Fig. 5(c) and Fig. 5(d). Independent Gaussian noise $\sigma_{\log(w_c)} = 0.5$ is added to the log-transform of $w_c$ values (Section 1.4.2) and $\sigma_{P_{obs}} = 25$ (dimensionless for simplicity) is added to $P_{obs}$ values. These noise magnitudes are chosen based on the range of output reference values from our simulator for these properties.

## 3.2. RJMCMC results for a variable number of layers

This section presents the results of the transdimensional inversion using the RJMCMC approach, first considering well log data only. Note that we have also conducted comparisons with a fixed-dimensional MCMC using the same dataset, highlighting that choosing a fixed number of geological layers, as often done in reservoir grids, significantly influences the inversion results (see Appendix E). Subsequently, we contrast these results with a joint inversion that incorporates both

flow and well log datasets combining the RJMCMC algorithm and the cascaded Metropolis sampler.

For both log-only and joint inversions, the width of the uniform prior distribution is $n_L \sim U_d(1, 15)$ for layers and $k \sim U(0.08, 1200)$ mD for layer permeability values based on the permeability ranges of the reference model (slightly extrapolated to account for uncertainty). The starting point of the Markov chain is a model randomly drawn from these prior distributions. The initial Gaussian perturbations for the different proposals have standard deviations $\sigma_{\epsilon_k} = 100$ mD, $\sigma_{\epsilon_z} = 0.05$, and $\sigma_{\epsilon_b} = 1000$ mD.

### 3.2.1. Inversion of well log data

The results of the well log regression problem are presented in Fig. 6(a) after 40,000 iterations including a burn-in period of 10,000 iterations (25% of the total run length). Starting from a random initial model (blue) and noisy well log data (gray), we estimate the log-permeability contrasts of the reference model (black) by inversion. The resulting average model (red) is obtained by an average of log-permeability values of all accepted post burn-in realizations. This averaged solution appears as a smoothed curve, characterized by significant variability not only in log-permeability values but also in layering. The geometry of the reference model is globally captured, since at least five distinct layers are clearly visible. The range defined by one standard deviation of the posterior distribution on the log-permeability $\sigma_{\log(k)}$ (dashed green) corroborates this observation, as it includes a significant portion of the reference model. However, at some depths, the posterior mean remains significantly distant from the reference. Specifically, the thin layer located at $z = 0.9$ is not adequately identified, even though a small permeability decrease is observed. The mean absolute error (MAE) is calculated on the log-permeability to quantify the distance between the average model and the reference model (0.91), and the same distance relative to the standard deviation of the realizations ($0.85 \sigma_{\log(k)}$). Fig. 6(b) provides a 1D estimation of the posterior probability distribution for the presence of layer interfaces
<figure><img src="image_8.png" /><figcaption>Fig. 7. (a) Results of the cascaded RJMCMC joint inversion of well log and flow data, on 150,000 post burn-in iterations of 30 independent Markov chains, merged into one result, through a depth section. (b) Posterior probability density to have an interface at each depth represented in 1D, merged from the 30 Markov chains. (c) Histogram of the posterior ensemble $P(n_L|d_{obs})$ during post burn-in inversion. (d) Main settings and outputs of the joint inversion.</figcaption></figure>

at each depth. The positions of the interfaces in the reference model are detected, except for the thin layer of low permeability.

Fig. 6(c) shows the convergence of the Markov chain after burn-in during the inversion of well log data. The sampling behavior aligns with the expected characteristics of a Markov chain that has reached convergence towards the posterior solution: The data misfit, as well as the MAE of log-permeability between the current model and the reference, which stabilized during the burn-in period, remain steady after the burn-in phase. The same observation is made for the global acceptance rate, which stays around 15%. Upon inspecting the acceptance rates for each specific perturbation, we note that the changes of layer permeability are accepted approximately 13% of the time, the move acceptance rate remains around 20%, while the birth and death acceptance rates are both approximately 15%. These percentages are consistent with effective RJMCMC sampling (Brooks et al., 2003). Fig. 6(d) presents the evolution of the number of layers $n_L$, and the histogram of the sampled posterior discrete distribution $P(n_L|d_{obs})$, after the burn-in phase. This number regularly oscillates around the reference value ($n_L = 7$), underscoring the natural parsimony property of the RJMCMC algorithm.

### 3.2.2. Joint inversion of well log and flow data

The joint inversion is computed with the cascaded sampler (Section 2.3), and using parallel Markov chains and a temperature factor $T$ to help overcome convergence issues (Section 2.4). A total of 30 parallel chains are run over 200,000 iterations each. We start with $T = 120$ with a decrease factor $f = 0.9999$ at each iteration, meaning that 47,873 iterations are needed for the temperature to decrease to $T = 1$, which we round to 50,000 iterations for the burn-in period.

The merged results in Fig. 7(a) from parallel Markov chains show several effects: The log-permeability values of the reference model layers are better recovered by the posterior average (red) than in Fig. 6(a). In particular, the very thin layer of low permeability at $z = 0.9$ is partly recovered thanks to the flow data. However, this layer is poorly resolved as the average is thicker than in the reference. The interval of one standard deviation (dashed green) is sometimes reduced, meaning

lower local variability in sampling (for instance at $z = 0.4$), but is quite similar overall to the one obtained by considering only well log data. The MAE of log-permeability between the reference model and the average of the ensemble (0.86) is slightly lower than that observed after well log inversion (0.91), but is higher relative to the standard deviation of the realizations (1.06 $\sigma_{log(k)}$) compared to the error of well log inversion (0.85 $\sigma_{log(k)}$). This suggests that, on average across the entire domain, the samples are closer to the reference model, but they are relatively more dispersed compared to their own variability than when inverting for well log only. The posterior distribution for the location of interfaces according to depth (Fig. 7(b)) provides a high probability value for all actual layer contacts (even if the approximations for the thin layer mentioned above are observed again), indicating the value of information added by production data to the problem.

Due to the non-linearity of the inverse flow problem, the consideration of flow information makes the chain convergence more difficult. Hence, the burn-in period is strongly variable from one chain to the next, and the arbitrarily chosen value of 50,000 iterations is sometimes insufficient to cover the number of iterations required for convergence. Although the global acceptance rate remains roughly the same (14%), the move rate decreases to 10% and birth and death rates decrease drastically and drop to 1%. On the contrary, the acceptance rate for permeability changes increases to approximately 24%. Another notable point is the introduction of a higher number of layers than in the previous inversion. Exceeding the reference number of 7 layers, the sampling is less parsimonious (Fig. 7(c)). A potential explanation for the difficulties of chain convergence, and the excessive generation of layers may lie in the coarse mesh used for flow simulations, thus producing a theoretical error that is not fully accounted for in the likelihood function (Section 4 includes more discussion on this aspect).

To quantify the computational resources required, the results of the joint inversion were obtained through parallel processing on a standard workstation powered by an Intel Xeon E3-1240 v6 CPU running at 3.70 GHz (4 cores, 8 threads), resulting in an average computational time of approximately 2.5 s per sample for a given chain. It should be noted that the computation of a single sample for the well log
<figure><img src="image_9.png" /><figcaption>Fig. 8. (a) (b) (c) Separate and joint inversion results through a depth section. (d) Discrete Kullback-Leibler divergence between prior and posterior permeability distributions for flow-only, log-only, and joint transdimensional inversions.</figcaption></figure>

inversion takes about 1 ms, so the average time for a joint inversion sample is substantially reduced by the cascaded approach when fast early rejections of models occur at the first step. The table in Fig. 7(d) summarizes the main information about the joint inversion.

### 3.2.3. Two complementary data types

In addition to the previous results, we run an inversion based only on flow data, computed under the same conditions as the joint inversion (Fig. 8(a)). As expected, flow data appear to only provide information about large scales with significant uncertainty, and become more precise when there are substantial fluid flow rate variations. In contrast, well logs (Fig. 8(b)) impose strong local constraints on the model but may miss some thin features. Utilizing both data types jointly (Fig. 8(c)) enhances the reliability of the results.

To quantitatively assess the contribution of each dataset on the outcomes, we compute a discrete Kullback-Leibler divergence (KLD) between prior and posterior log-permeability distributions along depth for the log-based, flow-based, and joint inversions (Fig. 8(d)). The KLD is a measure of how one probability distribution $P$ diverges from a second expected probability distribution $Q$. In other words, it quantifies the “distance” between two probability distributions in terms of value of information. For discrete distributions, the KLD is defined as

$$D_{KL}(P \parallel Q) = \sum_{i} P(i) \log \left( \frac{P(i)}{Q(i)} \right), \qquad (21)$$

where in our case $P$ is the posterior distribution and $Q$ the prior distribution.

Here, a larger value in KLD means a larger information gain relative to the prior. While seeming moderately informative when used alone, flow data still provide additional information relative to the prior distribution, supplementing well log data at some specific depths, especially within the two layers with very low permeability. This is because the simulator is particularly sensitive to high permeability contrasts. However, the exclusive use of well log data brings slightly more information within the thick central layer. This observation can be attributed to the pronounced uncertainty of flows in this layer, which slightly reduces the divergence from the uniform prior distribution, as well as to a likely trade-off between flow and well log information. Nevertheless, the joint use of data produces information that is better spatially constrained, combining the fine detection of layer permeability changes provided by well log data, and the highlighting of strong flow rate variations from flow data.

## 4. Discussion

### 4.1. Numerical diffusion in the flow simulation

The results of the joint inversion have shown a drastic decrease in the acceptance rates of interface births and deaths, and an average number of layers per model that is higher than the expected value. Finally, some strata, such as the thin layer of low permeability (at $z = 0.9$), are identified but poorly resolved.

These effects may be related to the coarse discretization of the domain. Whereas the flow data noise level in the likelihood function was chosen consistently with the noise added in the reference data, it may underestimate the numerical errors due to limited mesh resolution. A minor perturbation in the model can result in a distinct remeshing outcome, leading to a significantly different solution and a substantial modeling error. As this theoretical modeling error is currently excluded from the likelihood function, the noise magnitude is underestimated, resulting in overfitting and an inflated number of layers compared to the actual scenario (Sambridge et al., 2013).

Trials with simulations on finer meshes have highlighted notable differences in results, suggesting a numerical diffusion effect. A change for a finer mesh could potentially correct or mitigate these errors, but at the cost of increased computational time. As an alternative, the modeling error could be estimated by conducting simulations with various mesh resolutions, or using a Hierarchical Bayes framework (Section 4.3.3), and then incorporated into the likelihood function. This assumption also paves the way for upscaling strategies. For instance, an initial simulation could be computed on a coarse mesh, followed by a more refined mesh-based simulation only if the acceptance criterion is met in the preliminary one (Mondal et al., 2010).

### 4.2. Convergence

After the complete run of the 30 Markov chains in the joint inversion, some differences are still present between chains, especially in the lower part of the model as observed in Fig. 9(a). Fig. 9(b) presents the evolution of all the misfit functions of flow data after the burn-in period. It shows that the convergence is not complete for all chains as some oscillate at relatively large misfit levels. The observed difficulties in achieving convergence can be attributed to the underestimation of modeling errors (Section 4.1), but also to the disparate natures and sensitivities of the two datasets. Static data
<figure><img src="image_10.png" /><figcaption>Fig. 9. (a) Results of the cascaded RJMCMC joint inversion (average model of permeability ± one standard deviation) for the 30 independent Markov chains. (b) Flow misfit functions for the 30 Markov chains. (c) Gelman–Rubin convergence diagnostic with post burn-in PSRF values of well log, flow, and joint data misfit functions.</figcaption></figure>

relate to model parameters through a piecewise constant regression, whereas the relationship between reservoir permeability and dynamic data is highly non-linear, likely leading to multiple local misfit minima. Consequently, convergence is substantially determined by noise parameters and by the sensitivities of each data type. Whereas convergence is mathematically proven, the Markov chain could in practice spend significant time sampling around a local maximum of the posterior PDF (i.e., minimized misfit value), failing to adequately explore the global model space.

The Gelman–Rubin convergence diagnostic (Gelman and Rubin, 1992) on the misfit functions of well log data, flow data, and joint data, is very useful for evaluating convergence at different stages of the sampling process. The Gelman–Rubin diagnostic assesses multiple Markov chain convergences by comparing within-chain to between-chain variances to calculate the Potential Scale Reduction Factor (PSRF) $\hat{R}$:

$$ \hat{R} = \sqrt{\frac{\frac{T-1}{T}W + \frac{1}{T}B}{W}} \qquad (22) $$

where T is the chain length, B is the between-chain variance (the variance of the chain means), and W is the within-chain variance (the average of the variances of each chain). The numerator represents an estimate of the posterior variance for the considered parameter, based on the samples. The PSRF indicates a high degree of convergence to the target distribution when it decreases close to 1. Conventionally, a cut-off value of 1.2 indicates good convergence (Brooks and Gelman, 1998).

Fig. 9(c) displays the outcome of this statistic for the misfit functions at each post burn-in iteration. The result appears coherent because the PSRF values decrease over time, indicating reducing variances among the chain misfit functions. The PSRF for the well log misfit stabilizes faster than for the flow and joint misfits. This is actually explained by the use of the cascaded Metropolis and the temperature in the likelihood evaluation of flow data: During first steps, all chains rapidly converge towards the same result because the flow data have almost no impact. However, as the temperature gradually decreases, the flow data bring new information and lead to divergences in the speed of convergence. This effect is also reflected in the rise of the

PSRF for the well log misfit in the later post burn-in iterations (> 80,000 iterations), because the system slightly degrades the fit to the well log data to improve the fit to the flow data. Although the curve shapes are similar, their values do not reach the cut-off of 1.2, but rather remain above 1.5. This suggests the presence of local minima or artifacts within the models, potentially causing permanent differences between chains (visible in Fig. 9(a)). These results are common in non-linear problems where convergence is challenging. Generally, it is better to interpret the trends of the curves rather than using absolute cut-off values (Somogyvári and Reich, 2020), thus the chains can be considered as having converged.

To improve the low acceptance rate of birth and death proposals and enhance convergence, a conceivable solution could be the use of parallel tempering (e.g., Scalzo et al., 2019), adaptive MCMC (e.g., Lykkegaard, 2022), or gradient MCMC methods (e.g., Zhao and Sen, 2021). The gradient in flow and transport problems is not directly available but could be computed from adjoint state (Ackerer et al., 2014; Delay et al., 2019).

### 4.3. The “inverse crime” and extension to real-world scenarios

We acknowledge the “inverse crime,” wherein the same forward and noise models are employed both for generating synthetic data and for the subsequent inversion. Additionally, the inversion is applied to a one-dimensional layer-cake model with constant properties within each layer. Given that subsurface reservoirs and aquifers are typically complex, heterogeneous, and non-Gaussian (Zhou et al., 2014), these simplified layer-cake models fail to capture all the geological and physical complexities observed in reality, where more sophisticated parameterizations are required. Moreover, significant errors in measurement noise and model representation can accumulate with real-world scenarios. This section discusses the limitations of this approach and explores ways to account for more realistic use cases.

#### 4.3.1. Geological complexity

In the modeling of geological concepts with reservoir and hydrogeological data, uncertainty is omnipresent. The model structure is not known accurately, data contain errors, and modeling errors commonly
occur at the prediction stage. Realizing a realistic incorporation of geological information is challenging, as boundaries of model formations are uncertain, and main structural and stratigraphic features such as faults and fractures, paleochannels, or depositional patterns, may remain unknown (Carrera et al., 2005). Geostatistical approaches and suitable parameterizations have been developed to represent the small-scale heterogeneities of the field (Linde et al., 2015). While optimization methods commonly used in industry are efficient (e.g., Doherty and Hunt, 2010), they do not produce sufficient quantification of uncertainties, rely on regularization schemes to stabilize the solution, and typically provide only one best model (Moore and Doherty, 2006).

Stochastic simulations, such as sequential data assimilation (e.g., Aanonsen et al., 2009) and Monte Carlo simulations, offer a better way to account for parameter uncertainty (Zhou et al., 2014). However, stochastic simulations are generally not handled optimally by reservoir engineers because the history matching problem is ill-posed, the number of equivalent solutions is large, and the geometric model used is subject to geological interpretation. Furthermore, in real-world applications, it is common for the model to be mis-specified, meaning that the true solution may not lie within the parameter space defined by the chosen model. This occurs when the parameterization of the model does not fully capture the complexities or variability of the real system, often leading to biased predictions and an underestimation of uncertainties (Doherty and Welter, 2010; Finsterle and Zhang, 2011).

Recent developments are moving towards the use of high-dimensional parameterizations and complex prior information to enhance the geological realism of posterior model realizations (Hansen et al., 2012; Heße et al., 2019). However, the complexity of underground flow transport forward problems calls for convenient parameterizations with a reduced number of parameters (Linde et al., 2017). Instead, our approach considers geological uncertainties through a transdimensional inverse problem where the number of geological objects is variable, with a piecewise constant parameterization that involves a few model parameters. It allows for the inference of both physical and geometrical parameters with sufficient degrees of freedom, where no informative prior knowledge is necessary, reducing the risk of overfitting to a mis-specified model. By allowing the model structure to vary, the transdimensional method enables a more complete exploration of the parameter space and quantifies uncertainty in both the parameter values and the model structure, thus offering a more comprehensive view of possible system configurations (Sambridge et al., 2006). This approach marks a departure from the classical sequential reservoir modeling approach, where layers are first defined by experts and then left untouched during history matching. In contrast, the visualization of the posterior horizon depth distribution (Fig. 7(b)) and Kullback-Leibler divergence (Fig. 8(d)) appears as a rigorous way to evaluate the contribution of well log and flow data to geological interpretation. Consequently, the proposed methodology shows significant potential for geologically informed numerical well test interpretation in non-periodic and non-stationary media.

Nonetheless, the underlying one-dimensional geological model, consisting of horizontal layers, does not accurately capture the complexities typically found in real-world scenarios. This choice of representation is questionable at the scale of subsurface reservoirs and even in the well neighborhood, which may include many angular unconformities due to erosion and non-deposition. For more realism, Visser et al. (2019) and Herrero et al. (2023) propose to extend the parameterization of the layered model by adding an angle parameter for layer interface slope. The dipping angle enables the handling of geometric uncertainties associated with complex geological structures in a two-dimensional problem. This approach, grounded in seismic interpretation and well logging data, could help to address uncertainties in well correlation where the shape of sedimentary wedges is unknown (Baville et al., 2022). Another promising avenue is the use of level-set methods to implicitly define layer shapes, offering greater flexibility in modeling complex subsurface reservoir geometries (Caumon et al., 2012; Wellmann and Caumon, 2018), capitalizing on Legentil et al. (2022) for local mesh updating.

### 4.3.2. Physical complexity

Regarding the physics employed, it is dictated by simplifying assumptions: Fluid displacement depends on a single parameter in the inverse problem, rock permeability, assumed to have a constant value within layers. However, natural porous formations are inherently heterogeneous, and spatial variability influences flow and transport processes (Dagan, 2012). This assumption could be relaxed by considering vertical and lateral trends within each layer, permeability anisotropy, and spatial heterogeneity to produce more realistic petrophysical properties.

Effective medium models, as discussed in the context of anisotropic permeability in laminated systems (Zhan et al., 2022), offer a way to account for heterogeneities. This is particularly relevant when interpreting well log data, where the vertical resolution is often insufficient to capture fine laminations, potentially leading to errors if isotropy is assumed. These effective medium models could be integrated into the transdimensional inversion process to account for key anisotropic effects. Effective medium approaches simplify the problem by averaging properties over certain scales but may not capture localized variations that can be critical for multiphase flow and transport processes. Therefore, a fundamental challenge remains in appropriately choosing between explicit and effective (equivalent) representations of subsurface heterogeneities. Finding the right balance between model complexity and computational feasibility is also crucial for successful inversion of real-world data.

In complement to effective anisotropic permeability, explicit random fields could be used in each layer to generate realistic heterogeneities. However, integrating such random fields in the transdimensional inverse problem appears as computationally challenging. A more feasible approach could be to represent the permeability within each layer using spatially varying functions, such as linear or higher-order gradients. For instance, introducing a lateral linear gradient in permeability within each layer would involve adding one additional parameter per layer to capture the lateral variation of permeability. Alternatively, the permeability within layers could be modeled using polynomial functions, where the degree of the polynomial controls the complexity of the spatial variation. This strategy can provide flexibility in capturing more complex permeability variations while moderately increasing the number of parameters. Other petrophysical properties, such as rock porosity, could also be included in the inverse problem to enhance inversion results and model realism.

### 4.3.3. Noise estimation

In the synthetic case study presented in Section 3, we used the reference noise parameters in the likelihood function (Eq. (3)). However, in real-world scenarios, the true noise parameters are unknown, and their estimation is necessary. If the noise covariance does not accurately reflect the real noise in the system, the model may adjust its parameters to compensate for this misrepresentation, leading to biased or incorrect parameter estimates and an underestimation of uncertainty (Oliver and Alfonzo, 2018). In transdimensional inversion, this effect can be amplified since there is a trade-off between model complexity (i.e., the inferred number of layers) and the noise model. Underestimating the noise leads to the generation of too many layers, overfitting the data, while overestimating the noise results in overly smooth models with too few parameters (Bodin et al., 2012b; Sambridge et al., 2013).

Several approaches can be employed to address the problem of noise variance estimation. In many cases, noise parameters can be estimated using Maximum Likelihood Estimation, where the noise is treated as a parameter to be optimized along with the model parameters (Carrera and Neuman, 1986). Another method is cross-validation, in which the dataset is divided into training and validation subsets and the noise level is adjusted based on how well the model fits unseen data. Empirical Bayesian methods offer an alternative in which the variance of the residuals between observed and predicted data is used to infer the level of noise without needing to specify it a priori (Huang et al., 2006).
Regarding the covariance of errors, one common approach is to use a Gaussian Process or variogram inference. This technique is widely used in geostatistics, where noise is assumed to have a structured spatial correlation, and the model fit accounts for both measurement noise and intrinsic spatial variability (Chiles and Delfiner, 2012).

A common issue with these methods is that they tend to focus primarily on measurement error, without clearly differentiating it from modeling error. However, a parameterization that fails to accurately represent the geological reality of the heterogeneous medium can lead to significant model errors, as discussed in Section 4.3.1. A simple and effective solution to address this limitation is Hierarchical Bayes. This approach enables the simultaneous inference of model parameters and noise characteristics directly from the data by treating noise amplitude, spatial correlation, or both, as unknown hyper-parameters (Bodin et al., 2012a). By allowing the data to directly inform the noise model, Hierarchical Bayes can capture complex noise structures, such as spatial correlations, and account for uncertainties related to measurement errors and model approximations. However, this approach increases the number of parameters to infer, expanding the model space to sample, and introduces new parameter perturbations in the RJMCMC proposal function, which can make convergence more challenging in complex problems.

#### 4.3.4. How to invert for full physics-based simulation models?

In summary, transdimensional inversion of flow data yields promising results, but for simplicity the study is confined to a 1D parameterization. By combining the previous points (Sections 4.3.1, 4.3.2, and 4.3.3), we consider that the methodology developed in this work can be extended to 2D or 3D parameterizations.

However, making this workflow applicable to such real-world cases remains challenging. Full-physics simulations in 3D subsurface models typically require several hours to run, rendering the execution of tens of thousands of samples impractical without optimization. Furthermore, more complex parameterizations will increase the number of geometric and petrophysical parameters in the inverse problem, which is likely to make convergence more complex due to the curse of dimensionality and the non-linearities of the considered physics.

Two ways can be considered to reduce computational costs. The first is to optimize the sampling process of the RJMCMC algorithm, reducing the number of iterations needed for convergence, as discussed in Section 4.2. The second strategy, which complements the first, is to reduce the computational time of the forward problem. In this context, proxy models and surrogate techniques are particularly relevant solutions. Proxy and surrogate models serve as simplified representations of complex simulations, designed to approximate the output of full-scale simulations with significantly reduced computational costs. These models are useful in iterative processes, such as history matching, where repeated evaluations of flow simulations would be prohibitively expensive. For instance, Asher et al. (2015) provide a review of surrogate models in groundwater modeling, focusing on reducing computational costs while maintaining adequate accuracy. Dachanuwattana et al. (2019) include a proxy-based history matching workflow for unconventional oil reservoirs within an MCMC framework, while Ma et al. (2022) utilize a convolutional network model for surrogate-based history matching to reduce the cost of production forecasting. Similarly, Tang et al. (2022) introduce a deep-learning-based surrogate model that combines flow and geomechanical processes for efficient CO₂ sequestration simulation.

With these solutions, future studies could extend the transdimensional approach to diverse application fields such as CO₂ storage (Cho and Jun, 2021; Legentil et al., 2023) or geothermal forecasting (Dashti et al., 2023), both in 2D and 3D contexts.

### 4.4. Comparison with existing history matching techniques

In the context of reservoir characterization and history matching, ensemble-based data assimilation techniques such as the Ensemble Kalman Filter (EnKF) (Evensen, 2009) and Ensemble Smoother with Multiple Data Assimilation (ES-MDA) (Emerick and Reynolds, 2013) have shown their ability to handle high-dimensional problems. These methods update an ensemble of random field realizations to match observed data, providing a practical approach for uncertainty quantification in large-scale reservoir simulations.

In comparison, the primary drawback of the transdimensional inversion method proposed in this work is its computational intensity. In MCMC algorithms, the need to run a large number of forward simulations to adequately sample the model space can be a significant limitation when dealing with complex, high-dimensional models. Moreover, the inclusion of birth and death steps in the RJMCMC sampler increases the number of iterations required for the chain to converge to the equilibrium distribution.

However, ensemble-based methods often rely on assumptions that may limit their applicability. A key limitation is the assumption of Gaussianity in the posterior distribution of parameters and the linearity of the model responses, which may not hold in complex, non-linear systems. When the Gaussian error assumption is not satisfied, the performance of ensemble methods can decline, sometimes resulting in ensemble collapse, even with larger ensemble sizes (DeChant and Moradkhani, 2012). Additionally, these methods typically operate within a fixed model structure and do not account for uncertainties in model dimension or topology. While some advancements have attempted to address structural uncertainties by dynamically updating reservoir geometry associated with horizons (Seiler et al., 2010b) or faults (Seiler et al., 2010a.), they still operate within a fixed model dimension and topology and do not allow for changes in the number of layers or their configuration.

Although transdimensional inversion has not yet been applied in conjunction with ensemble methods, it complements these fixed-dimensional approaches by using adaptive parameterization and variable model dimensions. This enables rigorous quantification of structural uncertainties. Furthermore, the transdimensional method does not rely on multi-Gaussian assumptions for model parameters and states, which provides more flexibility in representing complex, non-Gaussian variables. This is particularly relevant in our context, where permeability is often governed by discrete (hence non-Gaussian) geological features such as layers, fractures and other geological bodies.

Thus, the methodology presented in this work pursues different objectives but could serve as a complementary approach to ensemble-based methods. Combining both approaches could be a promising direction by integrating the adaptive parameterization of transdimensional methods with the computational efficiency of EnKF or ES-MDA within each layer.

## 5. Conclusion

The complexity of forward problems in subsurface porous media calls for suitable parameterizations with a reduced number of parameters to solve the inverse problem. We choose in this work to address the question of stratigraphic uncertainties by acting directly on the number of geological layers. While usually a fixed model geometry is considered to perform Bayesian inversion of reservoir data, we propose a transdimensional joint inversion. This formulation, based on RJMCMC and local remeshing, allows the geometry to vary by dimension jumps in order to discover the number of layers needed to explain available reservoir observations. Although the computational cost is higher compared to other history matching methods, such as ensemble-based approaches, this method offers the advantage of adaptive model parameterization and addresses the uncertainty about the level of complexity needed to fit the data. The model parameters are
constrained both by well logs (with a simple regression problem) and by flow data (through a history matching problem), which are sequentially integrated into a joint inversion using a cascaded Metropolis sampler. The transdimensional cascaded sampler can be extended to integrate other types of data.

The first application of this method to a one-dimensional layered porous model with non-homogeneous layer thicknesses demonstrates its ability to recover the main geological discontinuities from per- meability well logs when the appropriate level of noise is applied. Furthermore, incorporating flow data into the problem improves hori- zon detection, reduces variability, and allows for the spatial assessment of the value of various types of information. This work paves the way for reservoir and aquifer modeling with realistic parameteriza- tions that enable the exploration of model parameter uncertainties and dimensionality.

### CRediT authorship contribution statement

**Julien Herrero:** Writing – original draft, Visualization, Validation, Software, Resources, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **Guillaume Caumon:** Writing – review & editing, Validation, Supervision, Methodology, Funding acquisition, Conceptualization. **Thomas Bodin:** Writing – review & editing, Vali- dation, Supervision, Methodology, Conceptualization. **Jeremie Giraud:** Writing – review & editing, Supervision, Conceptualization.

### Funding sources

This work was performed in the frame of the RING project (http://ring.georessources.univ-lorraine.fr/) supported by the RING Consor- tium managed by ASGA. No sponsor was involved in the study design, analysis and interpretation of the data. Jeremie Giraud was funded by the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 101032994.

### Declaration of competing interest

The authors declare that they have no known competing finan- cial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

We would like to thank for their support the industrial and academic sponsors of the RING-GOCAD Consortium managed by ASGA. The authors are also thankful to Mustapha Zakari for his help with the flow simulator RINGFlow, and to Capucine Legentil who created the LUMOS software for local mesh updating.

## Appendix A. Generating well log noise vector

The synthetic noise vector of permeability is generated separately in each layer by sampling from a normal distribution defined with a zero mean and covariance matrix of errors $C_{kobs}$. See Appendix C for the calculation of $C_{kobs}$. To sample from this distribution, we diagonalize the covariance matrix $C_{kobs}$ using Singular Value Decomposition (SVD):

$$C_{kobs} = USV^T. \qquad (A.1)$$

Given $C_{kobs}$ symmetry and positive definiteness, $U$ and $V^T(=U^{-1})$ perform the bijective transformation and its inverse for a coordinate system change, where $S$ represents a diagonal matrix in the decompo- sition. In the transformed coordinate system, the normal distribution becomes uncorrelated, simplifying the generation of a random vector $S\epsilon$, where $\epsilon = [\epsilon_1, \dots, \epsilon_n]^T$ is a vector consisting of independent standard normal random variables (or Gaussian white noise). The noise vector is then reverted to the physical coordinate system with

$$\Sigma = US\epsilon. \qquad (A.2)$$

## Appendix B. Flow simulator

In this paper, we use a Control Volume Finite Element Method (CVFEM) scheme implemented in the RINGFlow code (Anquez et al., 2020).

We consider a two-phase, isothermal, incompressible, Newtonian and immiscible system (Aarnes et al., 2007) consisting of a pure water phase ($\omega$) and a hydrocarbon phase ($\alpha$). Rock is assumed to be incomp- ressible. The fraction volume phase occupied in a pore volume is its saturation $S_\alpha$ ($\alpha = \omega, \alpha$). The pressure drop between the two phases $P_{cow}(P_{cow} = P_\omega - P_\omega)$ is the capillary pressure. For simplicity, we assume filled pores ($S_\omega + S_\alpha = 1$), no capillary pressure ($P_\omega = P_\alpha = P$) and we neglect gravity effects.

Within each phase, continuity equations give

$$\frac{\partial(\phi \rho_a S_\alpha)}{\partial t} + \nabla \cdot (\rho_a u_a) = M_\alpha \quad \alpha = \omega, \alpha, \qquad (B.1)$$

where $S_\alpha$ is the phase saturation, $\phi$ is the rock porosity, $\rho_a$ is the phase density (kg.m$^{-3}$), $u_a$ is the Darcy velocity (m.s$^{-1}$), $M_\alpha$ is the sink or source term (kg.m$^{-3}$.s$^{-1}$).

Assuming the rock and the two fluid phases are incompressible, Eq. (B.1) reduces to

$$\phi \frac{\partial(S_\alpha)}{\partial t} + \nabla \cdot (u_\alpha) = \tilde{M}_\alpha \quad \alpha = \omega, \alpha, \qquad (B.2)$$

where $\tilde{M}_\alpha = M_\alpha/\rho_\alpha$ is the volumetric source term (kg.s$^{-1}$).

We pose $\lambda_\alpha$ the phase mobility as

$$\lambda_\alpha = - \frac{k k_{ra}}{\mu_\alpha}, \qquad (B.3)$$

where $k$ corresponds to the rock permeability (m$^2$), $k_{ra}$ is the fluid relative permeability, and $\mu_\alpha$ is the phase viscosity (Pa.s$^{-1}$).

Summing up the two phases using Eq. (B.2) and Darcy's law, we obtain the pressure equation:

$$\nabla \cdot (-\lambda_T \nabla P) = M, \qquad (B.4)$$

where $M = \tilde{M}_\omega + \tilde{M}_\alpha$ and $\lambda_T = \lambda_\omega + \lambda_\alpha$ is the total mobility defining the total velocity:

$$u_T = -\lambda_T \nabla P. \qquad (B.5)$$

Water continuity (Eq. (B.2)) combined with the pressure (Eq. (B.4)) leads to

$$\phi \frac{\partial(S_\omega)}{\partial t} + \nabla \cdot (f_\omega u_T) = \tilde{M}_\omega, \qquad (B.6)$$

where $f_\omega = \lambda_\omega/\lambda_T$ is the water fractional mobility. We observe that, in this formulation, the water velocity is a fraction of the total velocity.

In what follows, we use a simple analytical model for saturation- dependent relative permeabilities:

$$\begin{cases} k_{rw} = S_\omega^2 \\ k_{ro} = (1 - S_\omega)^2. \end{cases} \qquad (B.7)$$

Finally, we set no flow boundary condition which means that the boundaries of the spatial domain are impervious.

In hydrocarbon recovery, this equation is generally solved numeri- cally using the finite-volume method (Aziz and Settari, 1979; Schäfer, 2006), while the finite-element method is common in hydrology.

Our numerical model has been developed for 2D triangular meshes. Following the CVFEM-DFM approach of Monteagudo and Firoozabadi (2004) and the discretization scheme of Karimi-Fard and Firoozabadi (2001), RINGFlow integrates equations of pressure (Eq. (B.4)) and saturation (Eq. (B.6)) on node centered median dual control volumes.

The key elements of the numerical formulation are (Anquez et al., 2020):

* The input permeability field function $k(x)$, in Darcies, is the 2D anisotropic tensor $[k_{xx} \quad k_{xz} \quad k_{zz}]$ computed for each model triangle.
* The pressure is interpolated in each triangle using a Lagrange P1 finite element interpolation.

* The saturation arguments of the fractional function $f_w$ are discretized spatially using an upwind scheme.

* Only water is injected and both oil and water are produced (the water production source term is given by $f_w$ times the pressure equation production term).

The discrete saturation and pressure equations are integrated over time $t$ with a time step $\Delta_t$. Pressure and saturation are coupled with the Implicit Pressure, Explicit Saturation (IMPES) method.

## Appendix C. Correlated noise model and likelihood evaluation for well log data

We employ a covariance function $C_k$ of permeability data errors expressed in terms of a noise correlation (symmetric, diagonal constant) matrix and a noise magnitude (Agostinetti and Malinverno, 2010) to generate the correlated noise model. Assuming the noise is stationary, meaning constant magnitude and correlation over space, $C_k$ can be presented within each layer as

$$C_k = \sigma_k^2 \begin{bmatrix} 1 & c(h_1) & c(h_2) & \dots & c(h_{n_{kL}}) \\ c(h_1) & 1 & c(h_1) & \dots & c(h_{n_{kL}-1}) \\ c(h_2) & c(h_1) & 1 & \dots & c(h_{n_{kL}-2}) \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ c(h_{n_{kL}}) & c(h_{n_{kL}-1}) & c(h_{n_{kL}-2}) & \dots & 1 \end{bmatrix}, \quad (C.1)$$

where $\sigma_k^2$ denotes the constant noise variance (the magnitude of data noise) for well log data, and $c(h_i)$ is the noise correlation between two points separated by a distance $h_i$ in the current layer crossed by $n_{kL}$ points. We parameterize the covariance matrix using an exponential variogram model. The covariance function is given by

$$C_k(h_i) = \sigma_k^2 \exp \left( -\frac{h_i}{a} \right), \quad (C.2)$$

where $a = -1/\log(r_k)$ represents the range of the variogram for a constant correlation factor $r_k$ that describes the correlation between two samples separated by one unit distance, and the noise variance $\sigma_k^2$ is the sill of the variogram. This formulation characterizes the Gaussian noise with an exponential correlation structure. The correlation factor $r_k$, defined between 0 and 1, implies independent data noise as it approaches 0 (but cannot be exactly 0), and perfect correlation as it approaches 1.

To evaluate the misfit function (Eq. (2)), one must find the inverse of the covariance matrix $C_k^{-1}$. A significant benefit of the parameterization based on the exponential correlation function is that $C_k^{-1}$ has convenient analytical forms expressed by $\sigma_k^2$ and $r_k$. It can be demonstrated that the inverse of $C_k$ is a symmetric tridiagonal matrix (Malinverno and Briggs, 2004; Bodin, 2010):

$$C_k^{-1} = \frac{1}{\sigma_k^2(1 - r_k^2)} \begin{bmatrix} 1 & -r_k & 0 & \dots & 0 & 0 \\ -r_k & 1 + r_k^2 & -r_k & \dots & 0 & 0 \\ 0 & -r_k & 1 + r_k^2 & \dots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \dots & 1 + r_k^2 & -r_k \\ 0 & 0 & 0 & \dots & -r_k & 1 \end{bmatrix}. \quad (C.3)$$

The memory requirement for the inverse covariance matrix advantageously scales with the total number of points in a layer $n_{kL}$, and computation of the misfit function calls for order $n_{kL}$ operations.

To evaluate the likelihood probability defined in Eq. (3), the determinant of the covariance matrix $|C_k|$ is also required additionally to the misfit function. Again, Malinverno and Briggs (2004) provide

an analytical formula for the determinant, leveraging the tridiagonal property of the inverse covariance matrix:

$$|C_k|_L = \sigma_k^{2n_{kL}} \cdot (1 - r_k^2)^{n_{kL}-1}. \quad (C.4)$$

This covariance function yields stable analytical solutions for $C_k^{-1}$ and $|C_k|$, making it well suited for the implementation of a Hierarchical Bayes algorithm, which allows $\sigma_k$ and $r_k$ to be expressed as variables to be determined during the inversion. This can be useful when using real data as these values are typically unknown during inversion. For simplicity in this paper, the expression of the estimated $\sigma_k$ and $r_k$ remains constant, and the chosen values correspond to the parameters of synthetic well log data $\sigma_{kobs}$ and $r_{kobs}$. Accounting for the piecewise discretization of the layered model, the covariance determinant for the whole model is

$$\begin{aligned} |C_k| &= \sum_{i=1}^{n_L} |C_k|_i \\ &= \sigma_k^{2n_k} \cdot (1 - r_k^2)^{n_k-1}. \end{aligned} \quad (C.5)$$

This implies that two distinct models share the same covariance matrix determinant, provided $\sigma_k$ and $r_k$ are constant, and the final expression for the likelihood in the context of exponential correlated noise parameterization is given by

$$P(\mathbf{k}_{\text{obs}}|\mathbf{m}) = \sum_{i=1}^{n_L} \frac{1}{\sqrt{|C_k|_i (2\pi)^{n_k L}}} \exp\left\{-\frac{1}{2}\phi_{k_i}(\mathbf{m})\right\}, \quad (C.6)$$

that simplifies using Eqs. (14) and (C.5) as

$$P(\mathbf{k}_{\text{obs}}|\mathbf{m}) = \frac{1}{\sqrt{|C_k|(2\pi)^{n_k}}} \exp\left\{-\frac{1}{2}\phi_{\mathbf{k}}(\mathbf{m})\right\}. \quad (C.7)$$

## Appendix D. Evaluation of complete acceptance probability

Determining the acceptance probability $\alpha(\mathbf{m}'|\mathbf{m})$ (Eq. (18)) involves computing three ratios: The prior (Eq. (12)) ratio $P(\mathbf{m}')/P(\mathbf{m})$, the likelihood (Eq. (3)) ratio $P(\mathbf{d}_{\text{obs}}|\mathbf{m}')/P(\mathbf{d}_{\text{obs}}|\mathbf{m})$, and the proposal ratio $q(\mathbf{m}|\mathbf{m}')/q(\mathbf{m}'|\mathbf{m})$. Since the data noise and correlation factor remain constant for each generated model, the likelihood ratio for $n$ data can be simplified as follows:

$$\frac{P(\mathbf{d}_{\text{obs}}|\mathbf{m}')}{P(\mathbf{d}_{\text{obs}}|\mathbf{m})} = \frac{\exp\left\{-\frac{1}{2}\phi(\mathbf{m}')\right\}}{\exp\left\{-\frac{1}{2}\phi(\mathbf{m})\right\}}, \quad (D.1)$$

where $\phi$ is defined by Eqs. (14) and (16).

The proposal ratio is perturbation type-dependent. Neither the move nor the modification of layer permeability values causes a dimension change. Hence, the proposal ratio is one, as the forward and reverse proposal probabilities are symmetrical, i.e., $q(\mathbf{m}'|\mathbf{m}) = q(\mathbf{m}|\mathbf{m}')$. This reverts to the classical Metropolis-Hastings algorithm, where a simple perturbation from the proposal distribution is often employed, as in Eq. (19). Similarly, the prior ratio is unity. Therefore, the acceptance or rejection criterion simplifies when the model dimension remains fixed, with only the likelihood ratio remaining:

$$\alpha_k(\mathbf{m}'|\mathbf{m}) = \alpha_m(\mathbf{m}'|\mathbf{m}) = \begin{cases} \min \left[ 1, \frac{\exp\left\{-\frac{1}{2}\phi(\mathbf{m}')\right\}}{\exp\left\{-\frac{1}{2}\phi(\mathbf{m})\right\}} \right] & \text{if } k_i \in K \\ 0 & \text{otherwise,} \end{cases} \quad (D.2)$$

where $K = \{k_i \in \mathbb{R} | k_{\min} \le k_i \le k_{\max}\}$.

Calculating the acceptance probability becomes more complex when considering birth and death Gaussian proposal ratios. These computations require determining the reverse step probabilities, essentially assessing the probability of deleting an interface at a specific location when a birth occurs and vice versa for the death process. Once determined, the proposal ratio is combined with prior and likelihood ratios.
<figure><img src="image_11.png" /><figcaption>Fig. E.10. Results of the MCMC inversion of well log data with a fixed number of layers, on 30,000 post burn-in iterations through a depth section. Three average models of log-permeability ± one standard deviation are represented from accepted realizations with respectively 5, 10, and 30 layers.</figcaption></figure>

The detailed computation of this combination of probability ratios has been expounded by Bodin and Sambridge (2009), hence, we only present the final formulas for birth acceptance and death acceptance.

After simplifications, the expression for a birth reduces to

$$ \alpha_b(\mathbf{m}'|\mathbf{m}) = \begin{cases} \min \left[ 1, \frac{\sigma_{\epsilon_b} \sqrt{2\pi}}{\Delta_k} \cdot \exp \left\{ \frac{(k'_{n_L+1}-k_i)^2}{2\sigma_{\epsilon_b}^2} - \frac{\phi(\mathbf{m}')-\phi(\mathbf{m})}{2} \right\} \right] & \text{if } (n_L+1) \in I, \\ 0 & \text{otherwise,} \end{cases} \quad (D.3) $$

where $\Delta_k = (k_{max} - k_{min})$, $I = \{n_L \in \mathbb{N} | n_{Lmin} \le n_L \le n_{Lmax}\}$, $i$ denotes the layer in the current model where the interface birth takes place, $k_i$ is the previous permeability value at birth location, $k'_{n_L+1}$ is the new permeability value, and $\sigma_{\epsilon_b}$ is the proposal standard deviation for birth perturbation.

For the death step, we invert the prior ratio, yielding

$$ \alpha_d(\mathbf{m}'|\mathbf{m}) = \begin{cases} \min \left[ 1, \frac{\Delta_k}{\sigma_{\epsilon_b} \sqrt{2\pi}} \cdot \exp \left\{ -\frac{(k'_{n_L}-k_i)^2}{2\sigma_{\epsilon_b}^2} - \frac{\phi(\mathbf{m}')-\phi(\mathbf{m})}{2} \right\} \right] & \text{if } (n_L-1) \in I, \\ 0 & \text{otherwise,} \end{cases} \quad (D.4) $$

where $i$ denotes the removed layer from the current model, $j$ corresponds to the merged layer “containing” the deleted interface location in the proposed model, $k_i$ is the initial permeability value at death location, and $k'_{n_L}$ is the permeability value in the reshaped layer after the death.

## Appendix E. MCMC results for a fixed number of layers

In georesource exploration, a model geometry is generally constructed by discretizing a grid uniformly to capture contrasts of the reservoir parameters under investigation. This grid is based on a fixed number of petrophysical parameters (e.g., cell permeability) and geometrical features (e.g., geological layers, faults). One interesting aspect that can be tested is the ability to approximate the reference model (Fig. 5(a)) in the case of a traditional MCMC sampler with a fixed number of layers. To recover the reference model, we invert for permeability using only well log data. The permeability values of the layers are initially drawn from a prior distribution $k \sim U(0.08, 1200)$ mD, and then updated during the random walk by the proposal distribution.

We compare three inversions including different numbers of layers (5, 10, and 30 layers) and permitting variations in layer thickness by enabling interface “move” operations as described in Section 2.2. For each inversion, 40,000 iterations incorporating a burn-in period of 10,000 are performed. In the 1D depth section presented in Fig. E.10, we observe that, despite the inclusion of the appropriate noise parameters in the likelihood function, the average model of log-permeability built with 5 layers (green) is under-constrained, lacking sufficient parameters to capture critical information. A configuration of 10 layers (blue) yields more satisfactory results, which makes sense given that this number is higher than but close to the reference number of layers (7 layers), thereby allowing for a relatively good approximation of the reference model. Finally, the 30-layer setup (orange) tends to overfit the noise within the data, rendering it unrealistic. These findings underscore the importance of selecting a suitable number of layers to balance model complexity against the risk of overfitting or oversimplification. In comparison, the transdimensional inversion (Fig. 6(a)) correctly recovers most of the structures of the reference model, while avoiding making an arbitrary choice on the number of layers, instead allowing the algorithm to decide on an appropriate value for this number.

In conclusion, the comparative analysis between fixed-dimensional and transdimensional MCMC inversions illustrates the advantages of the RJMCMC algorithm in the context of our layered model. By determining an appropriate (parsimonious) number of parameters, these methods not only enhance convergence but also give results that are more geologically plausible, while simultaneously providing a measure of uncertainty regarding the number, position, and thickness of the layers.

## Data availability

Data and codes associated to this paper are accessible to the sponsors of the RJMCMC Consortium and upon request to the authors.
## References

Aanonsen, S.I., Nævdal, G., Oliver, D.S., Reynolds, A.C., Vallès, B., 2009. The Ensemble Kalman Filter in reservoir engineering-a review. SPE J. 14, 393-412. http://dx.doi.org/10.2118/1172742-PA.

Aarnes, J.E., Gimse, T., Lie, K.A., 2007. An introduction to the numerics of flow in porous media using MATLAB. In: Geometric Modelling, Numerical Simulation, and Optimization. Springer, pp. 265-306. http://dx.doi.org/10.1007/978-3-540-68783-2_9.

Ackerer, P., Trottier, N., Delay, F., 2014. Flow in double-porosity aquifers: Parameter estimation using an adaptive multiscale method. Adv. Water Resour. 73, 108-122. http://dx.doi.org/10.1016/j.adwwatres.2014.07.001.

Agostinetti, N.P., Malinverno, A., 2010. Receiver function inversion by trans-dimensional Monte Carlo sampling. Geophys. J. Int. 181, 858-872. http://dx.doi.org/10.1111/j.1365-246X.2010.04530.x.

Anquez, P., Zakari, M., Caumon, G., 2020. Comparing three DFN Simplification Strategies for Two-Phase Flow applications. In: ECMOR XVII 2020. pp. 1-21. http://dx.doi.org/10.3997/2214-4609.202035112.

Asher, M.J., Croke, B.F., Jakeman, A.J., Peeters, L.J., 2015. A review of surrogate models and their application to groundwater modeling. Water Resour. Res. 51, 5957-5973. http://dx.doi.org/10.1002/2015WR016967.

Aziz, K., Settari, A., 1979. Petroleum Reservoir Simulation, vol. 476, Applied Science Publishers.

Barker, J.W., Cuypers, M., Holden, L., 2001. Quantifying uncertainty in production forecasts: Another look at the PUNQ-S3 problem. SPE J. 6, 433-441. http://dx.doi.org/10.2118/74707-PA.

Baville, P., Apel, M., Hoth, S., Knaust, D., Antoine, C., Carpentier, C., Caumon, G., 2022. Computer-assisted stochastic multi-well correlation: Sedimentary facies versus well distality. Mar. Pet. Geol. 135, 105371. http://dx.doi.org/10.1016/j.marpetgeo.2021.105371.

Bayes, T., 1763. LII. An essay towards solving a problem in the doctrine of chances. By the late Rev. Mr. Bayes, F. R. S. communicated by Mr. Price, in a letter to John Canton, A. M. F. R. S. Philos. Trans. R. Soc. Lond. 53, 370-418. http://dx.doi.org/10.1098/rstl.1763.0053.

Bodin, T., 2010. Transdimensional Approaches To Geophysical Inverse Problems. (Ph.D.) thesis. Australian National University, http://dx.doi.org/10.25911/5d5feebb08e4.

Bodin, T., Sambridge, M., 2009. Seismic tomography with the reversible jump algorithm. Geophys. J. Int. 178, 1411-1436. http://dx.doi.org/10.1111/j.1365-246X.2009.04226.x.

Bodin, T., Sambridge, M., Rawlinson, N., Arroucau, P., 2012a. Transdimensional tomography with unknown data noise. Geophys. J. Int. 189, 1536-1556. http://dx.doi.org/10.1111/j.1365-246X.2012.05414.x.

Bodin, T., Sambridge, M., Tkalcíč, H., Arroucau, P., Gallagher, K., Rawlinson, N., 2012b. Transdimensional inversion of receiver functions and surface wave dispersion. J. Geophys. Res.: Solid Earth 117, http://dx.doi.org/10.1029/2011JB008560.

Brooks, S.P., Gelman, A., 1998. General methods for monitoring convergence of iterative simulations. J. Comput. Graph. Stat. 7, 434-455. http://dx.doi.org/10.1080/10618600.1998.10474787.

Brooks, S.P., Giudici, P., Roberts, G.O., 2003. Efficient construction of reversible jump Markov chain Monte Carlo proposal distributions. J. R. Stat. Soc. Ser. B: Stat. Methodol. 65, 3-39. http://dx.doi.org/10.1111/1467-9868.03711.

Carrera, J., Alcolea, A., Medina, A., Hidalgo, J., Slooten, L.J., 2005. Inverse problem in hydrogeology. Hydrogeol. J. 13, 206-222. http://dx.doi.org/10.1007/s10040-004-0404-7.

Carrera, J., Neuman, S.P., 1986. Estimation of aquifer parameters under transient and steady state conditions: 1. Maximum likelihood method incorporating prior information. Water Resour. Res. 22, 199-210. http://dx.doi.org/10.1029/WR022i002p00199.

Caumon, G., Gray, G., Antoine, C., Titeux, M.O., 2012. Three-dimensional implicit stratigraphic model building from remote sensing data on tetrahedral meshes: theory and application to a regional model of La Popa Basin, NE Mexico. IEEE Trans. Geosci. Remote Sens. 51, 1613-1621. http://dx.doi.org/10.1109/TGRS.2012.2207727.

Chiles, J.P., Delfiner, P., 2012. Geostatistics: Modeling Spatial Uncertainty. In: Wiley Series in Probability and Statistics, vol. 713, John Wiley & Sons.

Cho, Y., Jun, H., 2021. Estimation and uncertainty analysis of the CO₂ storage volume in the Sleipner field via 4D reversible-jump Markov-chain Monte Carlo. J. Pet. Sci. Eng. 200, 108333. http://dx.doi.org/10.1016/j.petrol.2020.108333.

Christie, M., Demyanov, V., Erbas, D., 2006. Uncertainty quantification for porous media flows. J. Comput. Phys. 217, 143-158. http://dx.doi.org/10.1016/j.jcp.2006.01.026.

Dachanuwattana, S., Yu, W., Sepehrnoori, K., 2019. An efficient MCMC history matching workflow using fit-for-purpose proxies applied in unconventional oil reservoirs. J. Pet. Sci. Eng. 176, 381-395. http://dx.doi.org/10.1016/j.petrol.2019.01.070.

Dagan, G., 2012. Flow and Transport in Porous Formations, second ed. Springer, Berlin, Heidelberg, http://dx.doi.org/10.1007/978-3-642-75015-1.

Dake, L.P., 2001. The Practice of Reservoir Engineering (Revised Edition). Elsevier.

Dashti, A., Korzani, M.G., Geuzaine, C., Egert, R., Kohl, T., 2023. Impact of structural uncertainty on tracer test design in faulted geothermal reservoirs. Geothermics 107, 102607. http://dx.doi.org/10.1016/j.geothermics.2022.102607.

DeChant, C.M., Moradkhani, H., 2012. Examining the effectiveness and robustness of sequential data assimilation methods for quantification of uncertainty in hydrologic forecasting. Water Resour. Res. 48, http://dx.doi.org/10.1029/2011WR011011.

Delay, F., Badri, H., Ackerer, P., 2019. Heterogeneous hydraulic conductivity and porosity fields reconstruction through steady-state flow and transient solute transport data using the continuous adjoint state. Adv. Water Resour. 127, 148-166. http://dx.doi.org/10.1016/j.adwwatres.2019.03.014.

Denison, D.G., Holmes, C.C., Mallick, B.K., Smith, A.F., 2002. Bayesian methods for nonlinear classification and regression. In: Wiley Series in Probability and Statistics, vol. 386, John Wiley & Sons.

Dettmer, J., Dosso, S.E., Holland, C.W., 2010. Trans-dimensional geoacoustic inversion. J. Acoust. Soc. Am. 128, 3393-3405. http://dx.doi.org/10.1121/1.3500674.

Doherty, J.E., Hunt, R.J., 2010. Approaches to highly parameterized inversion-A guide to using PEST for groundwater-model calibration. 2010-5169, US geological survey. http://dx.doi.org/10.3133/sir20105169.

Doherty, J., Welter, D., 2010. A short exploration of structural noise. Water Resour. Res. 46, http://dx.doi.org/10.1029/2009WR008377.

Doucet, A., Wang, X., 2005. Monte Carlo methods for signal processing: a review in the statistical signal processing context. IEEE Sig. Process. Mag. 22, 152-170. http://dx.doi.org/10.1109/MSP.2005.1550195.

Emerick, A.A., Reynolds, A.C., 2013. Ensemble smoother with multiple data assimilation. Comput. Geosci. 55, 3-15. http://dx.doi.org/10.1016/j.cageo.2012.03.011.

Evensen, G., 2009. The ensemble Kalman filter for combined state and parameter estimation. IEEE Control Syst. Mag. 29, 83-104. http://dx.doi.org/10.1109/MCS.2009.932223.

Fernández-Martínez, J.L., Fernández-Muñiz, Z., 2020. The curse of dimensionality in inverse problems. J. Comput. Appl. Math. 369, 112571. http://dx.doi.org/10.1016/j.cam.2019.112571.

Finsterle, S., Zhang, Y., 2011. Error handling strategies in multiphase inverse modeling. Comput. Geosci. 37, 724-730. http://dx.doi.org/10.1016/j.cageo.2010.11.009.

Gallagher, K., Charvin, K., Nielsen, S., Sambridge, M., Stephenson, J., 2009. Markov chain Monte Carlo (MCMC) sampling methods to determine optimal models, model resolution and model choice for Earth Science problems. Mar. Pet. Geol. 26, 525-535. http://dx.doi.org/10.1016/j.marpetgeo.2009.01.003.

Gelman, A., Rubin, D.B., 1992. Inference from iterative simulation using multiple sequences. Statist. Sci. 7, 457-472. http://dx.doi.org/10.1214/ss/1177011136.

Geyer, C.J., Møller, J., 1994. Simulation procedures and likelihood inference for spatial point processes. Scand. J. Stat. 21, 359-373.

Gómez-Hernández, J.J., Sahuquillo, A., Capilla, J.E., 1997. Stochastic simulation of transmissivity fields conditional to both transmissivity and piezometric data— I. Theory. J. Hydrol. 203, 162-174. http://dx.doi.org/10.1016/S0022-1694(97)00098-X.

Green, P.J., 1995. Reversible jump Markov chain Monte Carlo computation and Bayesian model determination. Biometrika 82, 711-732. http://dx.doi.org/10.1093/biomet/82.4.711.

Green, P.J., 2003. Trans-dimensional Markov chain Monte Carlo. In: Highly Structured Stochastic Systems, vol. 27, Oxford University Press, pp. 179-198.

Hansen, T.M., Cordua, K.S., Mogesaard, K., 2012. Inverse problems with non-trivial priors: efficient solution through sequential Gibbs sampling. Comput. Geosci. 16, 593-611. http://dx.doi.org/10.1007/s10596-011-9271-1.

Hastings, W.K., 1970. Monte Carlo sampling methods using Markov chains and their applications. Biometrika 57, 97-109. http://dx.doi.org/10.1093/biomet/57.1.97.

Herrero, J., Caumon, G., Bodin, T., Lacheux, M., 2023. Transdimensional sampling of two-dimensional layered geological models with variable slope: a proof of concept. In: 5th EAGE Conf. Pet. Geostat. EAGE, pp. 1-5. http://dx.doi.org/10.3997/2214-4609.202335016.

Heße, F., Comunian, A., Attinger, S., 2019. What we talk about when we talk about uncertainty. Toward a unified, data-driven framework for uncertainty characterization in hydrogeology. Front. Earth Sci. 7, 118. http://dx.doi.org/10.3389/feart.2019.00118.

Huang, C.F., Gerstoft, P., Hodgkiss, W.S., 2006. Uncertainty analysis in matched-field geoacoustic inversions. J. Acoust. Soc. Am. 119, 197-207. http://dx.doi.org/10.1121/1.2139075.

Jiménez, S., Mariethoz, G., Brauchler, R., Bayer, P., 2016. Smart pilot points using reversible-jump Markov-chain Monte Carlo. Water Resour. Res. 52, 3966-3983. http://dx.doi.org/10.1002/2015WR017922.

Karimi-Fard, M., Firoozabadi, A., 2001. Numerical simulation of water injection in 2D fractured media using discrete-fracture model. In: SPE Annu. Tech. Conf. Exhib. SPE. pp. SPE-71615. http://dx.doi.org/10.2118/71615-MS.

Kirkpatrick, S., Gelatt, Jr., C.D., Vecchi, M.P., 1983. Optimization by simulated annealing. Science 220, 671-680. http://dx.doi.org/10.1126/science.220.4598.671.

Lallier, F., Caumon, G., Borgomano, J., Viseur, S., Fournier, F., Antoine, C., Gen-tilhomme, T., 2012. Relevance of the stochastic stratigraphic well correlation approach for the study of complex carbonate settings: Application to the Malam-paya buildup (Offshore Palawan, Philippines). Geol. Soc. Lond. Spec. Publ. 370, 265-275. http://dx.doi.org/10.1144/SP370.12.

Legentil, C., Pellerin, J., Cupillard, P., Froehly, A., Caumon, G., 2022. Testing scenarios on geological models: Local interface insertion in a 2D mesh and its impact on seismic wave simulation. Comput. Geosci. 159, 105013. http://dx.doi.org/10.1016/j.cageo.2021.105013.
Legentil, C., Pellerin, J., Raguenel, M., Caumon, G., 2023. Towards a workflow to evaluate geological layering uncertainty on CO₂ injection simulation. Appl. Comput. Geosci 18, 100118. http://dx.doi.org/10.1016/j.acags.2023.100118.

Linde, N., Ginsbourger, D., Irving, J., Nobile, F., Doucet, A., 2017. On uncertainty quantification in hydrogeology and hydrogeophysics. Adv. Water Resour. 110, 166-181. http://dx.doi.org/10.1016/j.advwatres.2017.10.014.

Linde, N., Renard, P., Mukerji, T., Caers, J., 2015. Geological realism in hydrogeological and geophysical inverse modeling: A review. Adv. Water Resour. 86, 86-101. http://dx.doi.org/10.1016/j.advwatres.2015.09.019.

Liu, Y., Durlofsky, L.J., 2020. Multilevel strategies and geological parameterizations for history matching complex reservoir models. SPE J. 25, 081-104. http://dx.doi.org/10.2118/193895-PA.

Lykkegaard, M.B., 2022. Multilevel Delayed Acceptance MCMC with Applications To Hydrogeological Inverse Problems (Ph.D. thesis). University of Exeter, United Kingdom.

Ma, X., Zhang, K., Zhang, J., Wang, Y., Zhang, L., Liu, P., Yang, Y., Wang, J., 2022. A novel hybrid recurrent convolutional network for surrogate modeling of history matching and uncertainty quantification. J. Pet. Sci. Eng. 210, 110109. http://dx.doi.org/10.1016/j.petrol.2022.110109.

MacKay, D.J., 2003. Information Theory, Inference and Learning Algorithms. Cambridge university press.

Mahalanobis, P.C., 2018. On the generalized distance in statistics. Sankhyä: Indian J. Stat. Ser. A (2008-) 80, S1-S7.

Malinverno, A., 2002. Parsimonious Bayesian Markov chain Monte Carlo inversion in a nonlinear geophysical problem. Geophys. J. Int. 151, 675-688. http://dx.doi.org/10.1046/j.1365-246X.2002.01847.x.

Malinverno, A., Briggs, V.A., 2004. Expanded uncertainty quantification in inverse problems: Hierarchical Bayes and empirical Bayes. Geophysics 69, 1005-1016. http://dx.doi.org/10.1190/1.1778243.

Malinverno, A., Leaney, S., 2000. A Monte Carlo method to quantify uncertainty in the inversion of zero-offset VSP data. In: SEG Tech. Prog. Exp. Abstr. 2000. SEG, pp. 2393-2396. http://dx.doi.org/10.1190/1.1815943.

Manassero, M.C., Afonso, J.C., Zyserman, F.I., Jones, A., Zlotnik, S., Fomin, I., 2021. A reduced order approach for probabilistic inversions of 3D magnetotelluric data II: joint inversion of MT and surface-wave data. J. Geophys. Res.: Solid Earth 126, e2021JB021962. http://dx.doi.org/10.1029/2021JB021962.

Metropolis, N., Rosenbluth, A.W., Rosenbluth, M.N., Teller, A.H., Teller, E., 1953. Equation of state calculations by fast computing machines. J. Chem. Phys. 21, 1087-1092. http://dx.doi.org/10.1063/1.1699114.

Mondal, A., Efendiev, Y., Mallick, B., Datta-Gupta, A., 2010. Bayesian uncertainty quantification for flows in heterogeneous porous media using reversible jump Markov chain Monte Carlo methods. Adv. Water Resour. 33, 241-256. http://dx.doi.org/10.1016/j.advwatres.2009.10.010.

Monteagudo, J., Firoozabadi, A., 2004. Control-volume method for numerical simulation of two-phase immiscible flow in two-and three-dimensional discrete-fractured media. Water Resour. Res. 40, http://dx.doi.org/10.1029/2003WR002996.

Moore, C., Doherty, J., 2006. The cost of uniqueness in groundwater model calibration. Adv. Water Resour. 29, 605-623. http://dx.doi.org/10.1016/j.advwatres.2005.07.003.

Mosegaard, K., Sambridge, M., 2002. Monte Carlo analysis of inverse problems. Inverse Problems 18, R29-R54. http://dx.doi.org/10.1088/0266-5611/18/3/201.

Oliver, D.S., Alfonzo, M., 2018. Calibration of imperfect models to biased observations. Comput. Geosci. 22, 145-161. http://dx.doi.org/10.1007/s10596-017-9678-4.

Oliver, D.S., Chen, Y., 2011. Recent progress on reservoir history matching: a review. Comput. Geosci. 15, 185-221. http://dx.doi.org/10.1007/s10596-010-9194-2.

Oliver, D.S., Reynolds, A.C., Liu, N., 2008. Inverse Theory for Petroleum Reservoir Characterization and History Matching. Cambridge University Press, Cambridge, http://dx.doi.org/10.1017/CBO9780511535642.

Pyrcz, M., Deutsch, C., 2014. Geostatistical Reservoir Modeling. Oxford University Press.

Ringel, L.M., Jalali, M., Bayer, P., 2021. Stochastic inversion of three-dimensional discrete fracture network structure with hydraulic tomography. Water Resour. Res. 57, e2021WR030401. http://dx.doi.org/10.1029/2021WR030401.

Rosenthal, J.S., 2000. Parallel computing and Monte Carlo algorithms. Far East J. Theor. Stat. 4, 207-236.

Rwechungura, R., Dadashpour, M., Kleppe, J., 2011. Advanced history matching techniques reviewed. In: M.E. Oil Gas Show Conf. SPE, http://dx.doi.org/10.2118/142497, SPE-142497-MS.

Sambridge, M., Bodin, T., Gallagher, K., Tkalčić, H., 2013. Transdimensional inference in the geosciences. Philos. Trans. R. Soc. A: Math. Phys. Eng. Sci. 371, 20110547. http://dx.doi.org/10.1098/rsta.2011.0547.

Sambridge, M., Gallagher, K., Jackson, A., Rickwood, P., 2006. Trans-dimensional inverse problems, model comparison and the evidence. Geophys. J. Int. 167, 528-542. http://dx.doi.org/10.1111/j.1365-246X.2006.03155.x.

Sambridge, M., Mosegaard, K., 2002. Monte Carlo methods in geophysical inverse problems. Rev. Geophys. 40, http://dx.doi.org/10.1029/2000RG000089, 3-1-3-29.

Scales, J.A., Snieder, R., 1997. To Bayes or not to Bayes? Geophysics 62, 1045-1046. http://dx.doi.org/10.1190/1.6241045.1.

Scalzo, R., Kohn, D., Olierook, H., Houseman, G., Chandra, R., Girolami, M., Cripps, S., 2019. Efficiency and robustness in Monte Carlo sampling for 3-D geophysical inversions with Obsidian v0.1.2: Setting up for success. Geosci. Model Dev. 12, 2941-2960. http://dx.doi.org/10.5194/gmd-12-2941-2019.

Schäfer, M., 2006. Computational Engineering: Introduction To Numerical Methods, vol. 321, Springer, http://dx.doi.org/10.1007/978-3-030-76027-4.

Seiler, A., Aanonsen, S., Evensen, G., Lia, O., 2010a. An elastic grid approach for fault uncertainty modelling and updating using the Ensemble Kalman filter. In: SPE Europec EAGE Conf. Exhib. SPE, http://dx.doi.org/10.2118/130422-MS, SPE-130422-MS.

Seiler, A., Aanonsen, S.I., Evensen, G., Rivenaes, J.C., 2010b. Structural surface uncertainty modeling and updating using the ensemble Kalman filter. SPE J. 15, 1062-1076. http://dx.doi.org/10.2118/125352-PA.

Seillé, H., Visser, G., 2020. Bayesian inversion of magnetotelluric data considering dimensionality discrepancies. Geophys. J. Int. 223, 1565-1583. http://dx.doi.org/10.1093/gji/gga391.

Smith, A.F.M., 1991. Bayesian computational methods. Philos. Trans. R. Soc. Lond. Ser. A: Phys. Eng. Sci. 337, 369-386. http://dx.doi.org/10.1098/rsta.1991.0130.

Somogyvári, M., Jalali, M., Jimenez Parras, S., Bayer, P., 2017. Synthetic fracture network characterization with transdimensional inversion. Water Resour. Res. 53, 5104-5123. http://dx.doi.org/10.1002/2016WR020293.

Somogyvári, M., Reich, S., 2020. Convergence tests for transdimensional Markov chains in geoscience imaging. Math. Geosci. 52, 651-668. http://dx.doi.org/10.1007/s11004-019-09811-x.

Tang, M., Ju, X., Durlofsky, L.J., 2022. Deep-learning-based coupled flow-geomechanics surrogate model for CO₂ sequestration. Int. J. Greenh. Gas Control 118, 103692. http://dx.doi.org/10.1016/j.ijggc.2022.103692.

Tarantola, A., 2005. Inverse problem theory and methods for model parameter estimation. Soc. Ind. Appl. Math. http://dx.doi.org/10.1137/1.9780898717921.

Tierney, L., 1994. Markov chains for exploring posterior distributions. Ann. Stat. 22, 1701-1728. http://dx.doi.org/10.1214/aos/1176325750.

Visser, G., Guo, P., Saygin, E., 2019. Bayesian transdimensional seismic full-waveform inversion with a dipping layer parameterization. Geophysics 84, R845-R858. http://dx.doi.org/10.1190/geo2018-0785.1.

Vozoff, K., Jupp, D.L.B., 1975. Joint inversion of geophysical data. Geophys. J. Int. 42, 977-991. http://dx.doi.org/10.1111/j.1365-246X.1975.tb06462.x.

Wellmann, F., Caumon, G., 2018. 3-D Structural Geological Models: Concepts, Methods, and Uncertainties. In: Advances in Geophysics, vol. 59, Elsevier, pp. 1-121. http://dx.doi.org/10.1016/bs.agph.2018.09.001.

Zhan, Q., Shetty, S., Boyd, A., Liang, L., 2022. A general anisotropic effective medium model for laminated sequence. SPE J. 27, 2318-2333. http://dx.doi.org/10.2118/209600-PA.

Zhao, Z., Sen, M.K., 2021. A gradient-based Markov chain Monte Carlo method for full-waveform inversion and uncertainty analysis. Geophysics 86, R15-R30. http://dx.doi.org/10.1190/geo2019-0585.1.

Zhou, H., Gómez-Hernández, J.J., Li, L., 2014. Inverse methods in hydrogeology: Evolution and recent trends. Adv. Water Resour. 63, 22-37. http://dx.doi.org/10.1016/j.advwatres.2013.10.014.