# Downscaling FWI images for reducing structural uncertainty

Giusi Ruggiero, Paul Cupillard, Guillaume Caumon

► To cite this version:

Giusi Ruggiero, Paul Cupillard, Guillaume Caumon. Downscaling FWI images for reducing structural uncertainty. Geophysical Journal International, 2026, 244 (3), pp.1-15. 10.1093/gji/ggaf528. hal-05484981

HAL Id: hal-05484981

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

Submitted on 20 Feb 2026

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 CC BY 4.0 - Attribution - International License
# Downscaling FWI images for reducing structural uncertainty

Giusi Ruggiero,<fnref n="1" /> Paul Cupillard and Guillaume Caumon <fnref n="1" />,<fnref n="2" />

Accepted 2025 December 8. Received 2025 October 25; in original form 2025 June 20

## SUMMARY

Despite ongoing improvements in the accuracy and efficiency of seismic interpretation, resolving small-scale geological structures remains challenging due to inherent limitations. Indeed, seismic imaging techniques are based on frequency band-limited data, and therefore they provide smooth representations of the real complex geology. This means that, within the resolution limit of seismic data, there is room for alternative valid interpretations at the small-scale. We suggest the use of a downscaling (or inverse homogenization) approach to make quantitative inferences and reduce uncertainty on small-scale geological structures. In particular, from a smooth velocity model obtained through the full waveform inversion (FWI) technique, the downscaling inversion aims to recover all the finer scale models compatible with it. In this context, the non-periodic homogenization serves as forward operator to construct smooth long-wavelength equivalent models. This approach can effectively be applied to specific areas of interest where seismic imaging is deemed reliable but the resolution is insufficient to support decision making. In this study we present an application in which downscaling is used to properly detect a fault and estimate small-scale P-wave velocity variations. In particular, the approach uses geostatistical simulation to generate velocity model realizations in the pre-faulting state and a kinematic modelling approach to generate fault displacement. The inverse problem is solved with a Bayesian stochastic framework in which fault-related and stratigraphic uncertainties are jointly incorporated into the inversion process. The methodology, validated on a synthetic example, is shown to be a reliable approach to reduce uncertainty on both fine-layering velocity and fault structure. Despite the additional challenges to be considered for real-case studies, we believe this target-oriented approach represents an efficient and cost-effective strategy to locally reduce FWI interpretation uncertainty.

**Key words:** Bayesian inference; Inverse theory; Numerical modelling; Crustal imaging; Seismic anisotropy; Fractures, faults, and high strain deformation zones.

## 1 INTRODUCTION

Mapping fault systems in the earth's subsurface and estimating fault throw is vital for numerous scientific and practical applications. Indeed, it is well known that fault geometry, distribution and petrophysical properties primarily control crustal fluids migration and resistance to rupture. This has major implications in subsurface earth resources applications, such as hydrocarbon exploration and production (e.g. T. Manzocchi et al. 2010), radioactive waste and CO₂ storage (e.g. B. Dockrill & Z.K. Shipton 2010) and hydrogeological and geothermal systems (e.g. J. Fairley 2009). Hence, understanding fault structures is essential to ensure the integrity and performance of projects involving structurally complex reservoirs (e.g. F.L. Richards et al. 2015; G. Yielding 2015), as well as to mitigate the environmental impact and risks associated with them, such as gas leakage or induced seismicity (e.g. K.A. Gilmore

et al. 2022). Beyond resource exploration, high-resolution studies of faults are essential for understanding faults growth, evolution and interactions in earthquake rupture modelling, seismic hazard assessment studies and geodynamics (e.g. J.J. Holmes et al. 2021; Z. Zhang et al. 2022).

In this context, seismic reflection data have been widely used to investigate fault zone structures. In particular, faults are mainly detected through the mapping of reflection discontinuities in seismic migrated sections (e.g. B. Wernicke 1995). Nowadays, these interpretations are also derived from velocity models obtained by full waveform inversion (FWI), in addition to migrated images (A. Alghuraybi et al. 2024; D. Kumar & R. Ali 2024; H. Dinh et al. 2024). FWI is an advanced method of seismic imaging that aims to recover high-fidelity physical property models of the subsurface using the full recorded wavefield, in contrast to traditional velocity modelling workflows based only on first-arrivals [see J. Virieux &

<footnote n="1">¹ Université de Lorraine, CNRS, GeoRessources, F-54000 Nancy, France. E-mail: giusi.ruggiero@univ-lorraine.fr, paul.cupillard@univ-lorraine.fr, guillaume.caumon@univ-lorraine.fr</footnote>

<footnote n="2">²Institut Universitaire de France, Paris, France</footnote>
S. Operto (2009) for an overview]. Conventionally, it is solved as an optimization problem where the solution is a best-fitting model that minimizes an objective function based on data derivatives and a regularization term (e.g. R.-E. Plessix 2006).

Despite being a robust and powerful inversion approach, the ill-
posed nature of the inverse problem, which stems from insufficient
data coverage, limited data bandwidth and measurement uncertain-
ties, allows for the presence of distinct plausible models, leading to
geological interpretations uncertainty.

Several techniques have been developed for reducing the inverse
model's interpretation ambiguity through the introduction of con-
straints during the inversion. For instance, some studies propose to
parametrize the FWI problem using a level-set method to represent
the boundary of a geologic body (e.g. A. Kadu et al. 2017; J.B.
Muir & J. B. Tsai 2020; T. Dahlke et al. 2020; J.B. Muir et
2022). More recently, the use of deep learning to incorporate prior
information on geological structures or to directly learn the map be-
tween data and model domains has been investigated (e.g. Z. Zhang
et al. 2018; Y. Li et al. 2021; Z. Zhang & T. Alkhalifah 2022; F.
Wang & T. Alkhalifah 2024). However, these methods often rely on
geological structures that are too simple, or they lack uncertainty
estimations, making it difficult to assess the solution confidence.
Alternatively, various approaches of Bayesian probabilistic FWI
have been introduced to incorporate geologic prior knowledge and
to assess uncertainty on the solution (e.g. S. Singh et al. 2018; G.
Visser et al. 2019; L. Mosser et al. 2020; O. Kolbjørnsen et
2020; X. Zhao & A. Curtis 2024). Fully nonlinear Bayesian FWI,
which avoids Gaussian assumptions on the prior and posterior dis-
tributions, can be solved using sampling techniques such as Markov
Chain Monte Carlo (MCMC) (R.M. Neal 1993; W.R. Gilks et
1995; S. Brooks et al. 2011). For these techniques, the main limita-
tion is the increased computational demand associated with running
a large amount of forward models to evaluate the numerous samples
needed to describe the solution space. Reducing computational cost
often involves reducing the number of model parameters, leading
to a loss of resolution. For this reason, alternative methods need to
be considered to ensure the feasibility of uncertainty assessments in
practical applications, where reliable information at the small scale
is required.

Among all the limitations that make it hard to obtain the desired
high-resolution model with FWI, we focus our attention on the lim-
ited frequency-band of seismic data. Due to physical and computa-
tional limitations, only a limited seismic bandwidth can be consid-
ered in practice. Given a maximum frequency present in the seismic
data, the wavefield exhibits a corresponding minimum wavelength.
Seismic resolution is limited to one quarter of this wavelength ver-
tically and even less laterally (R. Sheriff 1985; A.R. Brown 2011;
D.A. Herron 2011). This defines an approximate lower practical
constraint on seismically resolvable fault offsets (T.S. Faleide et al.
2021; V. Dimmen et al. 2023), meaning that within the limits of
resolution of seismic data, there is room for alternative valid in-
terpretations. Hence, if we assume that the global minimum of the
objective function is achieved, deterministic FWI provides only a
smooth representation of the real complex structures, which is not
suited for the detection of discontinuities in the observed properties
at the small scale (T. Bodin et al. 2015; Y. Capdeville & L. Métivier
2018), such as those induced by faults with throws below vertical
seismic resolution.

The objective of this paper is to propose an approach to reduce
uncertainty on small-scale geological structures from FWI mod-
els by combining elastic and structural parameters. The approach
relies on the fact that a long-wavelength equivalent model of the

earth obtained via deterministic FWI is similar to a smooth model
obtained from the application of the homogenization operator to
a smaller scale heterogeneous model (e.g. G.E. Backus 1962; Y.
Capdeville & J.-J. Marigo 2007; T.H. Jordan 2015; P. Cupillard &
Y. Capdeville 2018). Indeed, in the context of wave propagation, the
homogenization operator smooths out the heterogeneities of scale
smaller than a given minimum wavelength. The result is, therefore,
a smooth equivalent medium (the medium ‘seen’ by the wavefield),
such that the limited frequency waveforms computed in the true
and homogenized medium are the same up to a controlled error (Y.
Capdeville & J.-J. Marigo 2007).

This study uses an inverse homogenization approach to locally
downscale a smooth model obtained through FWI (Section 2). Follow-
ing the ideas of T. Bodin et al. (2015), N. Hedjazian et al. (2021)
and T. Santos et al. (2024), the inverse homogenization (or down-
scaling inverse problem) uses the elastic homogenization as forward
operator (Section 3.1). In this frame, the smooth inverted model is
interpreted as a probability distribution over possible small-scale
models in order to reduce uncertainty on small-scale stratigraphic
and fault structures. For this purpose, the methodology integrates
a geostatistical technique for estimating small-scale velocity varia-
tions and a stochastic fault modelling technique for reducing struc-
tural uncertainty (Section 3.2). This allows to explicitly encode
structural and stratigraphic information into the downscaling in-
verse problem (Section 3.3). A schematic and visual representation
of the proposed methodology is shown in Fig. 1. The benefits of the
proposed method are illustrated in Section 4 on a synthetic problem
and results are discussed in Section 5 before concluding our work.

## 2 DOWNSCALING AS AN INVERSE HOMOGENIZATION PROBLEM

Downscaling generally refers to the process of inferring small-scale
information from large-scale variables. Within the scope of our
work, the downscaling serves to derive models of the earth elastic
properties at small scale, given an estimation of its smooth effective
elastic properties. In the context of seismic wave propagation, there
are several different approaches used to compute effective proper-
ties (e.g. V. Grechka 2003; B. B. Engquist et al. 2009; Y. Capdeville
et al. 2010b; K. Gao et al. 2015; T.H. Jordan 2015; A. Fichtner &
S. M. Hanasoge 2017). Here, we apply the non-periodic homoge-
nization developed by Y. Capdeville et al. (2010a, b) and L. Guillot
et al. (2010). This method is based upon the two-scale homogeniza-
tion theory for periodic media (E. Sánchez-Palencia 1980), which
involves an explicit scale separation between large and small scales
and the introduction of two space variables, one for the large scales
and one for the small scales. The extension to non-periodic media
was proposed by Y. Capdeville & J.-J. Marigo (2007), Y. Capdeville
et al. (2010b), L. Guillot et al. (2010) and successfully applied in
seismology to upscale general elastic media in 1-D (Y. Capdeville
et al. 2010a), then in 2-D (Y. Capdeville et al. 2010b) and finally
in 3-D (P. Cupillard & Y. Capdeville 2018). The method is based
on the minimum wavelength of the wavefield to separate the scales
and makes it possible to compute an effective medium, valid up to a
given maximum source frequency, from a given fine scale descrip-
tion of a medium (Y. Capdeville et al. 2020).

In the FWI context, the homogenization has already been used by
several authors (A. Fichtner et al. 2013; Y. Capdeville et al. 2013;
M. Afanasiev et al. 2016; Y. Capdeville & L. Métivier 2018). In
A. Fichtner et al. (2013), the non-periodic homogenization is used
to construct smooth long-wavelength equivalent models to apply a
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Overview of the proposed methodology. Downscaling is cast into a Bayesian framework where information coming from the observed data, a smooth effective medium derived from HFWI, is combined with prior knowledge on the small-scale elastic and structural properties to retrieve the posterior distribution of all possible small-scale models that fit the data set.</figcaption></figure>

full waveform inversion method that integrates data across a wide range of spatial scales and frequencies. In the work of M. Afanasiev et al. (2016), the homogenization is applied to constrain FWI gradient updates in the layered medium case, ensuring that fine-scale model adjustments preserve the long-wavelength characteristics of the seismic model. Y. Capdeville & L. Métivier (2018) extend results from the stratified media case of Y. Capdeville et al. (2013), and numerically show that FWI based on limited frequency band data can retrieve, at best, the smooth homogenized version of the reference model. Hence, they suggest to invert for a smooth effective medium through the homogenized FWI (HFWI), where the FWI solution is homogenized after each iteration or only after the last iteration. The authors also show that incorporating homogenization constraints can improve the robustness and convergence of FWI by effectively guiding the solution space.

However, the smooth effective models produced by HFWI leave
possible ambiguities for geological interpretations. For this rea-
son, N. Hedjazian et al. (2021) propose to solve the seismic imag-
ing problem as a two-stage inversion: first perform a deterministic
waveform inversion where the aim is to recover a macroscale ef-
fective elastic medium, based on HFWI; then carry out a second
inverse problem, the downscaling or inverse homogenization, to
properly interpret this macroscale model and resolve small-scale
cavities. Along the same lines, T. Santos et al. (2024) downscale
tomographic images of marble cake mantle models by using a gen-
erative adversarial network (GAN)-based sampler to generate prior
models.

This study contributes to the ongoing research on the downscal-
ing as an inverse homogenization problem by testing the method-
ology in the 2-D case of a fault characterization problem with
a fine-scale stratified medium. The main contribution is the for-
mulation of a downscaling Bayesian stochastic framework that
jointly incorporates stratigraphic and fault-related uncertainty into
the inversion process. Prior knowledge on the fine-scale elas-
tic and fault-related properties are provided by, respectively, a
Gaussian random field and a parametric model of the fault ge-
ometry and of the associated near-field displacement. The pos-
terior distribution of all possible fine-scale models that fit the
data set (i.e. a macroscale effective medium that can be obtained

through the HFWI method) provides a quantification of uncertainty
(Fig. 1).

A significant feature of the methodology is that it can be applied
to a selected subset of a large data set (red frame in Fig. 1), as-
sumed to have a relevant impact on decision making, thus reducing
the computational demand that arises in large domains. This can
be attributed to the fact that the data is an homogenized/effective
medium, and thus, the effective properties at one location depend
only on the fine-scale properties in the neighbourhood (N. Hedjazian
et al. 2021). The inclusion of a priori information about small-scale
geological structures, the low dimensionality of the model space and
a computationally affordable forward modelling which allows the
use of global search approaches, makes the proposed downscaling
strategy an efficient methodology to be used, post-FWI, for reducing
structural interpretation uncertainty.

3 METHODOLOGY

3.1 Forward modelling: elastic homogenization

This section provides an introduction and summary of the main prin-
ciples behind the homogenization of complex deterministic media,
with no intrinsic scale separation, for wave propagation problems
[see Y. Capdeville et al. (2020) for a detailed overview]. In this
context, the main purpose of the homogenization is to determine
the relevant fine-scale properties to be smoothed, so that the propa-
gating wavefield in the homogenized medium is the same as in the
original, fine-scale medium up to a controlled accuracy. In contrast
to the classical periodic two-scale homogenization which is limited
to periodic media (E. Sánchez-Palencia 1980), the non-periodic
homogenization is applied to media that do not present any natu-
ral scale separation, introducing an artificial scale according to the
minimum wavelength as an input parameter.

Assuming that in the far-field, the wavefield has a minimum
wavelength $\lambda_{\min} = V_{\min}/f_{\max}$, the scale of details that can be found
in the homogenized medium with respect to $\lambda_{\min}$ is controlled by
the adimensional parameter $\varepsilon_0 = \lambda_0/\lambda_{\min}$. The input constant $\lambda_0$
sets the separation between macroscopic (x) and microscopic (x/$\varepsilon_0$)
scales. Therefore, choosing $\varepsilon_0 < 1$ yields an effective medium that
<figure><img src="image_3.png" /><figcaption>**Figure 2.** (a) Original isotropic fine-scale layered medium shown in terms of the P-wave velocity $v_p$. (b) On the left, the effective medium obtained by homogenizing it at $\lambda_0 = 15$ m and represented with the first coefficient $C_{11}^*$. The medium does not contain fine scales smaller than $\lambda_0$. On the right, the total amount of anisotropy of the effective medium calculated with eq. (5). Low anisotropy values are related to weak velocity contrasts.</figcaption></figure>

includes some heterogeneity scales smaller than $\lambda_{min}$, and $\varepsilon_0 > 1$ an effective medium where some scales larger than $\lambda_{min}$ are removed.

For any elastic property of the medium or physical quantity which depends on both scales, such as the displacement $\mathbf{u}(\mathbf{x}, \mathbf{x}/\varepsilon_0, t)$, the homogenization aims at finding its effective version, $\mathbf{u}^*(\mathbf{x}, t)$, that only depends on the macroscopic scale $\mathbf{x}$. The effective displacement $\mathbf{u}^*$, solution of the effective wave equation with a scale separation parameter $\varepsilon_0$, can be found as the leading order term of the asymptotic expansion in $\varepsilon_0$ of the true displacement $\mathbf{u}$:

$$ \mathbf{u}\left(\mathbf{x}, \frac{\mathbf{x}}{\varepsilon_0}, t\right) = \mathbf{u}^*(\mathbf{x}, t) + \varepsilon_0 \chi\left(\mathbf{x}, \frac{\mathbf{x}}{\varepsilon_0}\right) : \mathbf{\epsilon}(\mathbf{u}^*)(\mathbf{x}, t) + \mathbf{O}(\varepsilon_0^2). \quad (1) $$

To the leading order in $\varepsilon_0$, the displacement does not depend on the microscopic scale $\mathbf{x}/\varepsilon_0$. To the first order, the displacement depends on the microscopic scale, but only through a 'site effect', that is the effect on the recorded displacement of structures local to the receiver. This effect is expressed by the tensor product of the strain $\mathbf{\epsilon}(\mathbf{u}^*)$ and the first-order corrector $\chi$ in eq. (1).

The effective wave equation is obtained from the classical wave equation by replacing the density $\rho$, elasticity tensor $\mathbf{c}$, stress tensor $\sigma$ and source term $\mathbf{f}$ by their effective quantities:

$$ \rho^* \ddot{\mathbf{u}}^* - \nabla \sigma^* = \mathbf{f}^*, \quad (2) $$

$$ \sigma^* = \mathbf{c}^* : \mathbf{\epsilon}(\mathbf{u}^*). \quad (3) $$

The non-periodic elastic homogenization $\mathcal{H}$ projects any fine-scale elastic medium $\mathbf{m} \in \mathcal{M}$ to its homogenized version $\mathbf{m}^*$ in the effective, finite dimensional model space $\mathcal{M}^*$:

$$ \mathbf{m}^* = (\rho^*, \mathbf{c}^*) = \mathcal{H}(\rho, \mathbf{c}). \quad (4) $$

Building the upscaling operator $\mathcal{H}$ requires to solve an elastostatic problem with periodic boundary conditions referred to as the cell problem. While this problem can be solved analytically for layered media, a numerical method needs to be applied for general media, based on finite elements (Y. Capdeville et al. 2010b; P. Cupillard & Y. Capdeville 2018) or on an iterative FFT scheme (Y. Capdeville et al. 2015). In this work, we choose the second approach, the Fast Fourier Homogenization, to represent the forward problem.

The resulting effective elastic medium is smooth (i.e. it does not contain scales smaller than $\lambda_0$) and fully anisotropic (6 independent elastic parameters in 2-D, and 21 parameters in 3-D), as fine-scale isotropic heterogeneities lead to anisotropic effective heterogeneities (G.E. Backus 1962). Fig. 2 shows an example of effective medium and corresponding anisotropy for a fine-scale layered

medium. The amount of anisotropy is defined as

$$ \frac{\| \mathbf{c}^{\text{iso}} - \mathbf{c}^* \|_2}{\| \mathbf{c}^{\text{iso}} \|_2}, \quad (5) $$

where $\mathbf{c}^{\text{iso}}$ is the orthogonal projection of $\mathbf{c}^*$ on the isotropic tensor subspace (J.T. Browaeys & S. Chevrot 2004; F.I. Fedorov 2013). In the following, we omit $\rho^*$ and we refer to the effective medium as the effective anisotropic elasticity tensor $\mathbf{c}^*$.

## 3.2 Parametrization: structural-velocity models

Because this downscaling study is concerned with understanding fault-related uncertainty, a central aspect of our method is the definition of geologically meaningful model parameters that can be inferred in a probabilistic framework. This section describes the methodology employed to stochastically simulate models of faulted stratigraphy in the 2-D case.

Inspired by the work of G. Godefroy et al. (2018), who extend the work of F. Georgsen et al. (2012) and G. Laurent et al. (2013), fault displacement is produced by using a direct kinematic modelling approach. As compared to physics-based geodynamic modelling, kinematic modelling is purely geometric and does not account for stress transfer, making it more computationally efficient. This represents a significant advantage for the integration into a probabilistic setting. In this approach, the discontinuous displacement of geologic horizons in the region surrounding the fault is expressed by a numerical fault operator. Before applying the fault operator, we simulate a continuous horizon geometry and we define at each point in the computational domain the elastic properties ($v_p, v_s$ and $\rho$) in the pre-faulting space using piecewise constant values in each layer and Gaussian random fields (GRF). Realizations of GRF are produced by geostatistical simulation (P. Goovaerts 1997; C. Lantuéjoul 2001) and can be used to reflect property variations at different scales, from local (small-scale heterogeneities) to regional (large-scale trends). Moreover, geostatistical simulation allows for integrating different data sources, such as well logs, seismic tomography and rock physics constraints, ensuring that the generated spatial distribution of the modelled properties honors the data (e.g. T.M. Hansen et al. 2006; D. Grana et al. 2023).

### 3.2.1 Velocity modelling

The first step is to parametrize the model with a regular Cartesian grid $\mathbf{X} = \{\mathbf{x}_i \in \mathbb{R}^2 | i = 1, \dots, N\}$, and subdivide the domain into $J$ initially flat rock units separated by $J - 1$ internal horizons. The
<figure><img src="image_4.png" /><figcaption>**Figure 3.** (a) 1-D macro-layered velocity field, obtained by defining a constant velocity in each macro layer. (b) Fine-layered velocity field, generated as described in eq. (6). (c) Corresponding vertical velocity profiles.</figcaption></figure>

P-wave velocity within the domain $v_p$ is decomposed into a macro-layered velocity field $w$, defined by a homogeneous velocity $w_j$ in each rock unit $j$ and a residual field simulated by generating samples from an anisotropic Gaussian random field $Z$ of zero mean and unit variance (Fig. 3). The macro-layered velocity model is traditionally built by defining discontinuities from structural interpretations and assigning velocity to layers based on tomography or average interval velocities from vertical seismic profiles data. The residual field $Z$ is added to the macro-layered velocity to generate thin-layering structures inside each major layer, which may also be constrained from well log data. For this study, $Z$ is obtained by turning bands simulation (TBS) (G. Matheron 1973; J.R. Carr 2002), performed using the gstlearn Python package from MINES Paris (gstlearn 2023). The total fine-layered velocity at $x$ can be thus expressed as

$$v_p(\mathbf{x}) = \sum_{j=1}^{J} w_j \mathbf{1}_j(\mathbf{x}) + \alpha Z(\mathbf{x}), \qquad (6)$$

where $\mathbf{1}_j(\mathbf{x})$ is the indicator function equal to 1 if $x$ belongs to unit $j$ and 0 elsewhere, and $\alpha$ is a scalar factor tuning the intensity of the random field $Z(\mathbf{x})$.

The characterization of the residual field $Z$ is governed by the parameters defining the variogram—such as correlation range and anisotropy—and describing its spatial variability. The specific parametrization of the simulation depends on the geological settings and on data available. For this work, we assume a flat stratigraphy in the geological model, so we consider only vertical changes in the velocity, which are controlled by the vertical correlation length $l_y$ of $Z$, and we neglect lateral velocity variations within each stratigraphic layer (Fig. 3).

Downscaling aims to reconstruct the 1-D fine-scale velocity features of $v_p$. In order to reduce the model space dimensions and to facilitate an efficient sampling, we consider a subset of points $s = \{s_k | k = 1, \dots, K\}$ where $v_p$ is inferred, as done in the pilot point method (G. de Marsily et al. 1984). The velocity field in the whole domain is then obtained with the geostatistical simulation

conditioned on the values $v_p(s)$. The choice of the $s$ grid locations is dependent on the vertical spatial variability of the velocity field defined by $l_y$. For the simplified test of this study, we use a set of deterministic values for the variogram parameters, including $l_y$. Hence, we assume that this information is known from other data sources. However, in a realistic scenario, these variables are subject to uncertainty, so they could be considered as part of the unknowns (e.g. P. Wang et al. 2022). The geostatistical model parametrization used in this study will be detailed in Section 4.

### 3.2.2 Fault displacement

Once the fine-scale velocity field $v_p$ has been defined within the medium, faulting is simulated by applying the kinematic fault operator to the set of points $X$. Applying the fault operator results in the mapping of each point’s location pre- and post-faulting and it is used in this work to generate models with continuous elastic properties that are offset by faults. Each fault is defined by a location, spatial extent of the fault-induced deformation, dip, slip and slip type (i.e. inverse or reverse fault). Following G. Godefroy et al. (2018), the kinematic model is characterized by a discontinuous displacement field defined within an elliptic shape centred on the fault (Fig. 4a). The spatial extent of the displacement is given by the ellipse’s axes semilength, controlling the distance between the centre and the tip of the fault along the fault line, $a$, and in the orthogonal direction, $b$ (Fig. 4a). The location of the fault in the domain and the fault dip are given by the coordinates of the centre $(x_c, y_c)$ and the rotation angle $\phi$ of the ellipse (Fig. 4a). The displacement magnitude is controlled by a maximum amount of slip $D_{max}$ assumed to be at the fault centre and two parametric profiles defining a displacement parallel to the fault line and discontinuous across it, decreasing smoothly from the centre to the boundary of the ellipse (Fig. 4b). These fault-related parameters $\theta = (x_c, y_c, a, b, \phi, D_{max})$ are treated as unknowns in the downscaling inverse problem and will be sampled to retrieve an ensemble of possible fault models, consistent with the FWI image.
<figure><img src="image_5.png" /><figcaption>Figure 4. (a) Elliptical fault space. The major-axis of the ellipse defines the dip direction, while the minor-axis defines the direction orthogonal to the fault. (b) Parametric profiles defining the displacement inside the elliptical fault space along the two directions (modified from G. Godefroy et al. 2018).</figcaption></figure>

As a final step, the post-faulting elastic values are interpolated bilinearly onto the regular grid X before forward modelling with the non-periodic elastic homogenization.

## 3.3 Bayesian inference with MCMC

In this study, the objective of the downscaling is to quantify uncertainty on both P-wave velocity $v_p$ and fault structure $\mathcal{V}$ of a fine-scale elastic medium $\mathbf{m}$, given a smooth effective medium $c_{\text{obs}}^*$. Hence, the model is defined as a structural-velocity model consisting of two components $\theta = \{\mathcal{V}, v_p\}$, where $\mathcal{V}$ are variables describing the structural model (e.g. the location, dip and slip of the fault), and $v_p$ describes the vertical fine-layered velocity profile.

For a given fine-scale elastic medium $\mathbf{m}(\theta)$, the corresponding effective medium $c^*$ is computed through the forward operator $\mathcal{H}$ as $c^* = \mathcal{H}(\mathbf{m}(\theta))$, where $\mathcal{H}$ is the non-periodic elastic homogenization operator discussed above. The observed effective medium $c_{\text{obs}}^*$ refers to a smooth and homogenized FWI model, achieved via HFWI, that gathers the information present in the seismic data and does not contain scales smaller than the homogenization wavelength $\lambda_0$ (Y. Capdeville & L. Métivier 2018; N. Hedjazian et al. 2021).

For incorporating prior knowledge and quantifying uncertainty on $\theta$, a Bayesian formalism is adopted to solve the downscaling inverse problem. The objective in Bayesian inference is to determine the posterior distribution $p(\theta | c_{\text{obs}}^*)$ of the parameter set $\theta$ according to the well-known Bayes equation:

$$p(\theta | c_{\text{obs}}^*) = \frac{p(c_{\text{obs}}^* | \theta)p(\theta)}{p(c_{\text{obs}}^*)}, \quad (7)$$

where $p(c_{\text{obs}}^* | \theta)$ is the likelihood function which quantifies how the model $\mathcal{H}(\mathbf{m}(\theta))$ fits the observations $c_{\text{obs}}^*$, $p(\theta)$ denotes the prior distribution of model parameters $\theta$, and $p(c_{\text{obs}}^*)$ denotes the prior distribution of measured data that is usually regarded as a normalizing constant. The posterior distribution in eq. (8) expresses the probability that $\theta$ provides a fine scale elastic medium $\mathbf{m}(\theta)$

whose effective medium is $c_{\text{obs}}^*$. Assuming the distribution of errors in the observed data to be multivariate-Gaussian distributed with zero mean and covariance matrix $\Sigma$, the likelihood is expressed as:

$$p(c_{\text{obs}}^* | \theta) \propto \exp \left( -\frac{1}{2} (c_{\text{obs}}^* - \mathcal{H}(\mathbf{m}(\theta)))^T \Sigma^{-1} (c_{\text{obs}}^* - \mathcal{H}(\mathbf{m}(\theta))) \right). \quad (8)$$

Prior information is generally inferred from other independent assessments, such as well log samples and rock-physics measurements. For the structural model, we choose uninformative uniform prior distributions for each fault-related parameter in $\mathcal{V}$, reflecting the assumption of little prior knowledge on the fault structure. The initial $\mathcal{V}$ will be randomly chosen from the prior distribution. For the velocity model, we assume a prior knowledge coming from indirect data on the macro-layered velocity field $\mathbf{w}$. The initial values of $v_p$ at the pilot point locations $\mathbf{s}$ are set equal to $\mathbf{w}$. Combining the prior and the likelihood we obtain the posterior probability distribution in eq. (7), which is then used to quantify the uncertainty of the parameters of interest $\theta = \{\mathcal{V}, v_p\}$ (A. Tarantola 2005).

### 3.3.1 Sampling algorithm

To circumvent the computationally challenging estimation of $p(c_{\text{obs}}^*)$ in eq. (7), a popular method for characterizing $p(\theta | c_{\text{obs}}^*)$ is a sampling algorithm known as Markov Chain Monte Carlo (MCMC). This stochastic process starts from a given initial guess of the model parameters, it samples repeatedly in the model space according to a proposal distribution, and eventually reaches the stationary regime sampling the posterior distribution.

For this study, we make use of an adaptive Metropolis sampling method (H. Haario et al. 1999), which combines the traditional Metropolis-Hastings (N. Metropolis et al. 1953; W.K. Hastings 1970), with the capacity to tune step sizes continuously while approaching the target distribution. At each iteration $i$, the MCMC algorithm generates trial moves from the current state $\theta^{i-1}$ to the
<figure><img src="image_6.png" /><figcaption>**Figure 5.** (a) Reference model parameters $\theta = \{\vartheta, v_p\}$, consisting of the structural model $\vartheta$ and the 1-D velocity model $v_p$. (b) Fine-scale elastic medium $\mathbf{m}(\theta)$ used as reference medium for the synthetic application. The elastic medium is illustrated here in terms of P-wave velocity $v_p$. The elastic properties $v_s$ and $\rho$, required for the elastic homogenization, are computed from $v_p$ as explained in Section 4.1. (c) First component of the effective elasticity tensor $\mathbf{c}^* = \mathcal{H}(\mathbf{m}(\theta))$ obtained by homogenizing the reference fine-scale medium with $\lambda_0 = 20$ m.</figcaption></figure>

proposal state $\theta^i$, with the probability of the proposal being accepted depending on the current sampled model only. Each sampled vector $\theta^i$ corresponds to a fine-scale medium realization $\mathbf{m}(\theta^i)$ (Figs 5a and b). The forward simulation (elastic homogenization) gives the effective elastic properties of these realizations $\mathcal{H}(\mathbf{m}(\theta^i))$ (Fig. 5c), allowing for the computation of the likelihood function (eq. 8). Specifically for the implementation in this study, the Markov chain consists of two possible updates:

(i) The velocity model: one point location $s_k \in \mathbf{s}$ is randomly selected and the velocity $v_p(s_k)$ is updated according to a Gaussian proposal distribution $q(v_p(s_k))$. To simplify the model and reduce computational complexity, the velocity at the selected pilot point is perturbed independently, without regard for the velocity values of other pilot points (this aspect is discussed in Section 5.1). The velocity in the whole medium $\mathbf{m}(\theta^i)$ is simulated conditioned to the whole set of pilot points $v_p(s)^i$ values. In this way, the Markov chain explores the parameter space by adjusting the values at the points $\mathbf{s}$, which in turn influences the simulated values at other locations.

(ii) The structural model: all the fault-related parameters $\vartheta = (x_c, y_c, a, b, \phi, D_{\max})$ are perturbed according to a Gaussian proposal distribution $q(\vartheta)$. By doing so, the fault structure is updated in the fine-scale elastic medium $\mathbf{m}(\theta)^i$.

The new candidate sample $\theta^i$ is accepted following the Metropolis-Hastings rule:

$$ \alpha(\theta^i | \theta^{i-1}) = \min \left[ 1, \frac{p(\mathbf{c}_{\text{obs}}^* | \theta^i) p(\theta^i) q(\theta^{i-1} | \theta^i)}{p(\mathbf{c}_{\text{obs}}^* | \theta^{i-1}) p(\theta^{i-1}) q(\theta^i | \theta^{i-1})} \right]. \quad (9) $$

When the proposal distribution is symmetric, that is $q(\theta^{i-1} | \theta^i) = q(\theta^i | \theta^{i-1})$, as in the case of a Gaussian distribution, and the prior is non-informative, the MH acceptance probability simplifies to the ratio of likelihoods of the proposed state $p(\mathbf{c}_{\text{obs}}^* | \theta^i)$ and current state $p(\mathbf{c}_{\text{obs}}^* | \theta^{i-1})$ of the chain. In practice, if a random number drawn uniformly between 0 and 1 does not exceed $\alpha(\theta^i | \theta^{i-1})$, then the proposal state $\theta^i$ is accepted as the new state of the Markov chain; otherwise, it is rejected and the new state is $\theta^{i-1}$.

# 4 SYNTHETIC APPLICATION

## 4.1 Set-up

The performance of the presented method is assessed on a synthetic data application. A model involving a single fault that offsets a fine-layered stratigraphy is used to represent the reference earth model (Fig. 5b). This model is parametrized with a regular grid where each cell has a 1 m × 1 m size and a total domain of 180 m × 100 m. To build the reference fine-scale elastic medium $\mathbf{m}(\theta_{\text{ref}})$, we first define a constant P-wave velocity in each macro-layer $w_j$, and then the covariance or variogram model describing the spatial structure of the Gaussian random field Z employed to simulate the fine-layering velocity model $v_p$ within the domain. Because we aim at reproducing a horizontal stratigraphic orientation, we select the variogram model parameters to account for anisotropy: the horizontal direction is chosen as the primary orientation where spatial correlation is stronger. Thus, we choose a Gaussian covariance function with vertical correlation length $l_y = 5$ m and a horizontal to vertical anisotropy ratio of 500 : 1, typical of deepwater or deltaic depositional environments (M.J. Pyrcz & C.V. Deutsch 2014). Following this, fault displacement is applied as described in Section 3.2.2. The S-wave velocity $v_s$ and density $\rho$ properties, also needed for the forward step, are calculated from $v_p$ using the fixed ratio $v_p/v_s = 1.8$ (this hypothesis is discussed in Section 5.1) and Gardner's relation $\rho = 0.31 \mathbf{v}_p^{0.25}$ (G. Gardner et al. 1974).

In this synthetic application, we approximate the observation (the smooth HFWI model) with the homogenized reference medium $\mathcal{H}(\mathbf{m}(\theta_{\text{ref}}))$, to which we add spatially correlated Gaussian noise $\eta$ to increase the complexity of the inverse problem. $\mathcal{H}(\mathbf{m}(\theta_{\text{ref}}))$ is computed from $v_p$, $v_s$ and $\rho$ by applying the homogenization operator with a scale separation parameter $\lambda_0 = 20$ m. As so, no spatial variations of wavelength smaller than $\lambda_0$ are kept in the smooth effective properties. For the correlated noise $\eta \sim N(0, \Sigma)$, we choose a Gaussian covariance function with characteristic length scales $l_x = l_y = 20$ m equal to $\lambda_0$, in order to reproduce possibly structured correlations that may be present in an elastic model obtained from HFWI, and noise variance equal to 4 per cent of the
<figure><img src="image_7.png" /><figcaption>**Figure 6.** The six components of the effective elasticity tensor $\boldsymbol{c}_{\text{obs}}^*$ representing the observation in the synthetic application. These are obtained by homogenizing the reference model at $\lambda_0 = 20$ m and by adding spatially correlated Gaussian noise $\eta$.</figcaption></figure>

mean value of the effective elastic properties. The resulting tensor $\boldsymbol{c}_{\text{obs}}^* = \mathcal{H}(\mathbf{m}(\boldsymbol{\theta}_{\text{ref}})) + \eta$ is shown in Fig. 6. As can be observed, the off-diagonal terms $\boldsymbol{c}_{13}^*$ and $\boldsymbol{c}_{23}^*$ of the effective properties, which give information about effective medium anisotropy, are non-zero only around the fault and the high-contrast velocity layer in the central part of the model. The small-scale layering is not producing effective anisotropy since we choose a cutoff wavelength $\lambda_0$ higher than the vertical spatial variation of the fine-scale velocity $l_y$.

Because the true noise structure is generally unknown in practice, to test the inversion under more realistic and imperfect conditions, we compute an approximate data covariance matrix $\boldsymbol{\Sigma}$ instead of using the exact covariance of the generated noise $\eta$ (an idealized scenario). In particular, $\boldsymbol{\Sigma}$ considers errors in the six elastic coefficients of the data vector $\boldsymbol{c}_{\text{obs}}^* = (\boldsymbol{c}_{11}^*, \boldsymbol{c}_{12}^*, \boldsymbol{c}_{13}^*, \boldsymbol{c}_{22}^*, \boldsymbol{c}_{23}^*, \boldsymbol{c}_{33}^*)$ are correlated, but does not account for spatial correlation between data points, meaning that errors are considered independent for each location.

The choice of the pilot points $\mathbf{s}$ where $\mathbf{v}_p$ is sampled significantly affects the feasibility and efficiency of the method. For this application, we choose to use points regularly spaced along the vertical direction with a spacing between two successive points $\Delta s = 2$ m, approximately half the vertical correlation length $l_y$, to control the recovery of the reference fine-scale velocity heterogeneities through the inversion process. The choice of this parametrization will be further discussed in Section 5.1.

## 4.2 Results

We apply the MCMC approach described in Section 3.3.1 to sample posterior models. To verify the convergence, we run four MCMC parallel chains starting from different initial models, shown in Fig. 9(a). The convergence plots of the four runs are shown in Fig. 7, which displays the negative log-likelihood, equivalent of the misfit function in Bayesian inference, and the acceptance rate of all the realizations. We discard the first 50 000 realizations, considered part of the burn-in period, and compute the posterior probability based on the remaining 180 000 realizations in the stationary regime, for which the acceptance rate is around 10 per cent. MCMC convergence is typically assessed by examining posterior distribution statistics across iterations (S.P. Brooks & A. Gelman 1998). A widely used diagnostic is the Gelman–Rubin statistic $\hat{R}$ (A. Gelman & D.B. Rubin 1992), which provides an estimate of how well the Markov chains have mixed and whether they have converged to the target distribution (typically when $\hat{R} < 1.1$). For all the fault-related parameters $\theta$, the estimated $\hat{R}$ values were close to 1, indicating convergence. For the velocities defined at the pilot points $\mathbf{s}$, convergence is inherently slower because only one pilot point is perturbed at each iteration. This behaviour is reflected by the small decreasing trend of the log-likelihoods after the burn-in (Fig. 7b). Consequently, although not all pilot point velocities achieved a Gelman–Rubin statistic below 1.1, we chose to terminate the sampling as a computational compromise.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** (a) Top graph: Negative log-likelihood functions of four parallel MCMC chains for number of iterations. The burn-in phase is considered to finish after 50 000 iterations when the functions reach a plateau and start sampling the posterior. Bottom graph: Acceptance rate of four parallel MCMC chains for number of iterations. As the number of iterations increases, the acceptance rates stabilize around 0.1. (b) Zoomed-in view of the two graphs post burn-in.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Inversion results of four MCMC chains. (a) Velocity profiles $v_p$ extracted at $x = 50$ m. The blue line is the reference fine-scale velocity profile. The black and red lines show the mean of the prior and posterior distributions. Red-shaded region marks the 90 per cent confidence interval (5th and 95th percentile range), capturing the most probable range of velocity values. (b) Marginal posterior densities of the fault-related parameters $\vartheta$ from four MCMC chains (shown in different shades of blue) after burn-in. Black dotted line indicates the true values from the reference model in Fig. 5(b). The minimum and maximum values of the uniform priors for each fault-related parameter are also indicated.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Example of starting and posterior models for different MCMC chains. (a) The initial fault parameters $\vartheta$ are randomly sampled from the uniform prior distributions, whereas the initial velocity model $v_p$ corresponds to the macroscale velocity field $w$. (b) Posterior models show a fault that matches the target structure (Fig. 5b), and only differs on some fine-scale velocity layers.</figcaption></figure>

The posterior distributions of the fault parameters $\vartheta$ are shown in Fig. 8(b). Despite using large non-informative priors, also given in Fig. 8(b), the reference value for all the fault-related parameters $\vartheta_{\text{ref}}$ (marked by the dashed line on the plots) lies within the posterior distribution. This indicates that data, as reflected through the likelihood function, have sufficiently informed the posterior, even in presence of noise $\eta$. In particular, the posterior densities suggest how the likelihood function tightly constrains the posterior for some fault parameters, leading to sharp and concentrated distributions and low standard deviation in the posterior samples. This is a sign that the model is confident in these estimates. In contrast, data are less informative for other parameters, such as for the semi-major *a* and semiminor *b* axis of the ellipse, defining the spatial extent of the fault-induced deformation. This is reasonable because the fault displacement is assumed to be maximum at the centre of the fault and decreases towards the tips of the elliptic fault zone. Hence, a small variation of the ellipse semi-axis lengths does not produce a significant change in the data. In addition to a reduction of uncertainty in the structural model, there is also a reduction of uncertainty in the velocity model. The inverted velocity is expressed as the posterior mean computed with posterior samples of the four chains. Fig. 8(a) shows velocity with depth at $x = 50$ m. The prior mean velocity reflects the macro-layered velocity (Fig. 3a), while the mean of the posterior velocity profiles follows the shorter wavelength trend of the reference model's velocity profile. However, the posterior mean velocity does not align with the reference profile at all depths. Some regions—more affected by data

noise—exhibit higher uncertainty, as indicated in Fig. 8(a) by a larger confidence interval, whereas others are better constrained. An animated figure provided as supplementary material illustrates how the simulated velocity profile evolves during MCMC sampling, shown at intervals of 1000 iterations. Examples of posterior fault-velocity models for different MCMC chains are shown in Fig. 9(b).

The mean, standard deviation and relative error maps are computed and displayed in Fig. 10. The relative error is calculated as the difference between the reference model and the mean model, divided by the standard deviation. These posterior statistics are confirming how the fault model is successfully recovered with low uncertainty, whereas a relatively large uncertainty is associated to the fine-layered velocity, especially in the upper and lower areas of the domain. Overall, these results demonstrate the ability of the methodology to resolve uncertainty on the unknown small-scale features in the synthetic case, in spite of the spatially correlated noise added to the reference homogenized model. The accuracy of the results stems from the fact that this correlated noise has a smaller horizontal variogram range than the prior horizontal range of the residual velocity field. Consequently, the proposed model samples manage to filter out part of the noise at the expense of a lower data fit. A complementary explanation for the obtained accuracy comes from the relatively strong control of the elastic parameters provided by the deterministic macro-layered velocity. Indeed, the continuous velocity residual field tends to preserve the strong anisotropy at the macro-layer boundaries.
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Posterior statistics post burn-in: mean, standard deviation and relative error of the small-scale medium $m(\theta)$ obtained from posterior samples of 4 MCMC chains. The fault structure has been successfully recovered and well-constrained, whereas higher uncertainty is observed for some fine-scale velocity layers.</figcaption></figure>

# 5 DISCUSSIONS

The proposed downscaling approach aims at capturing small-scale geological structures from a smooth effective model that gathers the information present in the seismic data. For the specific implementation adopted here, downscaling integrates geostatistical and kinematic fault modelling techniques to jointly account for stratigraphic and structural (fault-related) uncertainty.

This method is related to a series of previous works that attempt to bridge the gap between FWI models and geological interpretation. Indeed, it is well known that FWI focuses purely on matching wave physics, not necessarily producing geologically realistic results (M. Warner et al. 2013). A reason is that FWI outputs are constrained by data bandwidth and partial illumination, providing at best smooth representations of the real earth small-scale structures. Hence, these models are prone to interpretation ambiguities when it comes to interpreting geological structures associated to sharp discontinuities. Bridging this gap through joint inversion, regularization with geological priors and integrated modelling workflows remains an active area of research (e.g. V.C. Tsai et al. 2023; J. Giraud et al. 2024).

In the presented workflow, an homogenized FWI solution, assumed to provide the effective properties of the real earth model, is used to inform small-scale geological models with a probabilistic downscaling inversion, accounting for uncertainty. The presented downscaling approach can be used for all applications involving small-scale anomalies or heterogeneities causing discontinuities in elastic properties. Small-scale variations are mapped into effective anisotropy in the homogenized medium. This information is then exploited in the downscaling process to find its finer scale characteristics. Geological applications include, for example, imaging the subsurface to detect faults, fluids, karsts, tunnels or cavities.

This paper shows an application to a fine-layering medium characterized by fault displacement. Although the method has shown good performance in the synthetic test, capturing both fine-layering stratigraphy and fault-related variables, we are aware of some assumptions and simplifications made in the application setup as well

as of the additional challenges to address in the case of real data application. These aspects are further discussed below.

## 5.1 Parametrization

To parametrize the fault model, we adopt a kinematic operator where fault deformation is modelled in a symmetric ellipse centred on the fault line (as in F. Georgsen et al. 2012; G. Godefroy et al. 2018).

For simplicity, we test the approach with an isolated normal fault inducing a displacement symmetric and parallel to the fault line. However, this fault model can account for asymmetric displacement and more complex fault geometries, such as listric faults, by using asymmetric displacement profiles (Fig. 4b) (G. Godefroy et al. 2018). The downscaling inversion can also be tested in the case of several faults where fault displacements are applied sequentially and the relationship between faults is controlled by the sequence order (N. Cherpeau et al. 2010). In such case, the number of faults in the model and the order in which they truncate each other, for instance, can be considered uncertain and more sophisticated MCMC algorithm, such as rj-MCMC (e.g. M. Sambridge et al. 2006), or more informative priors could be adopted to address these problems (e.g. L. Holden et al. 2003; N. Cherpeau & G. Caumon 2015; O. Aydin & J.K. Caers 2017; F. Taty Moukati et al. 2025).

Regarding the velocity model, for the sake of simplicity, we decide to use a reference model in which the velocity varies only along the vertical direction to reproduce a horizontal stratigraphy and we derive the S-wave velocity from P-wave velocity values assuming a fixed ratio of $v_p/v_s = 1.8$. When only $v_p$ is known, a $v_p/v_s$ ratio of roughly 1.7 – 1.8 is commonly assumed in numerical simulations since many crustal rocks exhibit a Poisson’s ratio near $v \approx 0.25$, which corresponds to $v_p \approx 1.73 v_s$ (R.E. Sheriff & L.P. Geldart 1995; D. Komatitsch & J. Tromp 1999; H. Gercek 2007). Alternatively, empirical $v_p - v_s$ relationships derived from lab data, well logs or field measurements are adopted, such as the one proposed by J.P. Castagna et al. (1985) for clastic sediments or by T.M. Brocher (2005) for crustal rocks. Making this ratio spatially variable using an additional variable at pilot points could be envisioned for more
realism, but this option was not chosen here to harness the number of
model parameters and reduce the computational cost. The fine-scale
P-wave velocity is simulated with a geostatistical approach. Hence,
its spatial variability is defined by the vertical correlation length of
the variogram model. Having an accurate estimate of the geosta-
tistical variables describing the variogram model is determinant in
recovering the fine-scale velocity through downscaling inversion. In
this synthetic test, we assume that the variogram model is known,
which is a significant hypothesis. In a real case study, a common
practice is to fit a chosen model to an empirical variogram computed
using the available borehole data. In this case, the approach can be
easily extended to model lateral heterogeneities provided there is
sufficient data to constrain the horizontal correlation length. If in-
stead data are not available or insufficient to reliably constrain the
variogram model, variables defining the variogram (such as range,
sill and anisotropy) could also be treated as unknowns and estimated
by the inversion.

To infer the velocity profile with MCMC, we use a pilot point parametrization strategy, wherein the inversion focuses on a reduced set of velocities defined at selected point locations. The velocity field in the whole domain is then obtained with the geostatistical simulation conditioned on the velocity values sampled at the pilot point locations. This approach reduces the effective parameter space explored during MCMC sampling, while approximately honoring the spatial correlation structure. In this scheme, a critical choice is the number and spacing of pilot points, which affects the trade-off between model accuracy and computational feasibility. Indeed, including as many pilot points as possible for calibration leads to a good fit of observed data but increases the computational burden. In a stationary domain, the choice of the pilot points spacing should be guided by the spatial correlation structure of the velocity field (i.e. the variogram range) (e.g. J.J. Gómez-Hernández et al. 1997; A. Alcolea et al. 2006; R. Aissat et al. 2023). For the synthetic application of this work, we define regularly spaced pilot points with a spacing of approximately half the variogram range to adequately capture the spatial heterogeneity of the velocity field. For simplicity, in the MCMC algorithm we adopt a proposal scheme where the velocity at pilot points is perturbed independently (see Section 3.3.1). Because these points are spatially correlated (Δs < l<sub>y</sub>), independent perturbations could violate the assumed spatial structure of the variogram model if the perturbations are too large or uncorrelated with nearby points. We have assessed the limited impact of this simplification in our experiments by comparing the empirical variogram of conditional simulations from the posterior samples with the theoretical (assumed) model (see the online supplementary material for additional information). As an alternative, pilot points values could be sampled from the conditional probability distribution defined by kriging and the estimation variance computed from the neighbouring pilot points (X. Freulon & C. de Fouquet 1993).

## 5.2 Homogenization and geological heterogeneities scales

A further significant aspect to take into consideration is that the
described methodology is dependent on the cutoff wavelength λ₀
used to homogenize the elastic properties. In the synthetic test we
approximate the observation (homogenized FWI model) with the
homogenized reference model using λ₀ = 20 m, to which we add
Gaussian noise. Choosing λ₀ as four times bigger than the vertical
spatial correlation range of the velocity field, and considering a spa-
tially correlated Gaussian noise are meant to increase the complexity
of the inverse problem for the synthetic case. Indeed, for shorter λ₀,

geological heterogeneities at smaller scales would be preserved, and
the solution would be better constrained by the observations. Never-
theless, the fine-scale structures of the target medium are recovered.
This can be attributed to the significant effective anisotropy ob-
served at the fault position (Fig. 6), providing valuable information
to the downscaling procedure when constraining this structure. The
effective anisotropy is enhanced by the presence of a high-contrast
velocity layer in the reference model which is offset by the fault.
In addition, the use of the same geostatistical parametrization to
simulate fine-layering velocity in both predictions and observations
facilitates the accurate recovery of the velocity field.

In the case of real data application, the choice of λ₀ for homogenizing the FWI solution and the proposed fine-scale models is crucial. This wavelength has to be carefully chosen according to the minimum resolvable wavelength λmin in FWI. Theoretically, λmin is part of the unknowns but in practice, it can be computed approximately with the maximum frequency in the waveform data and the presumed minimum seismic velocity in the medium. Assuming that FWI can reach λmin/2 to λmin/4 resolution, all the heterogeneities at smaller scales would be poorly resolved by the seismic data. Hence, they could be smoothed out by choosing λ0 = λmin/2, that is choosing ε0 = 0.5 (see Section 3.1). This smoothing would ensure that the original and homogenized models produce similar wavefields.

An adaptive approach with progressive refinement of the ho-
mogenization wavelength λ₀ and pilot points grid could also be
adopted in a complex real case-study. In such case, the MCMC in-
version would start with a large λ₀ (smoother effective properties)
and a coarse pilot points grid (fewer points where velocity is in-
ferred) and, as the algorithm converges, gradually decrease λ₀ and
refine the pilot points grid. This multiscale strategy would allow
to progressively increase the resolution and the level of small-scale
details as the posterior distribution becomes better constrained, thus
improving convergence efficiency.

5.3 Data uncertainty

An additional difficulty to be considered for real data application is
the definition of data uncertainty. In the synthetic application pre-
sented above, we compute an approximate data covariance matrix
considering that errors in the six independent elastic coefficients
obtained by homogenization are correlated with each other. How-
ever, we do not account for spatial correlation between data points.
This approximation has been done for simplicity although spatial
correlation induced by the homogenization (variables are spatially
correlated on a length scale λ₀) should be accounted for by using a
full covariance matrix. For real data application, the estimation of
uncertainty in the observation, the HFWI solution model, is not triv-
ial. In particular, errors in HFWI can derive from theoretical errors
on the homogenization operator *H* and from errors in the FWI solu-
tion, which may come from noise in the data, lack of data coverage,
convergence towards a local minimum in the waveform inversion,
etc. How to compute uncertainty associated to the FWI solution
obtained through a local optimization approach is a widespread
subject of research in seismic tomography. In several works, this is
obtained by computing local uncertainties around the global mini-
mum associated to this solution through the inverse Hessian matrix,
which approximates the posterior covariance matrix (A. Tarantola
2005). Because an exact evaluation of the posterior covariance re-
quires the storage and inversion of Hessian-sized matrices, some
form of rank reduction are often used to generate tractable posterior
covariance estimates (e.g. T. Bui-Thanh et al. 2013; H. Zhu et al. 2016; Q. Liu & D. Peter 2019). However, in HFWI, the updated model is projected into the space of effective models $M^*$ using the homogenization operator after each iteration (or only at the final one). This means that we should account for uncertainties in $M^*$, which is challenging. A possible strategy to approximate the data covariance matrix would be, for instance, that of N. Hedjazian et al. (2021), who compute uncertainties performing multiple HFWI inversions for multiple noise realizations on the synthetic data. As such, this approximation only considers errors associated to the noise in waveforms and mapped to the HFWI solution through the inversion. A more accurate estimate of data uncertainties should be thoughtfully addressed in future work.

## 5.4 Computational aspects

The presented approach, in which forward seismic modelling is replaced by the homogenization technique, represents an efficient alternative to costly stochastic FWI. For large-scale FWI problems, the computational cost of solving the forward problem remains prohibitively expensive. In this work, elastic homogenization is performed with an FFT-based solver on a regular grid, which has the great advantage of being computationally affordable. The cost of one forward modelling for the 2-D model used in this study is less than one second on a laptop equipped with 20 cores. This makes the implementation of a MCMC approach feasible to sample the ensemble of possible solutions.

In general, MCMC sampling is known to become computationally intensive for high-dimensional parameter spaces due to the curse of dimensionality (A. Curtis & A. Lomax 2001). Nevertheless, our downscaling approach can be used as a target-oriented uncertainty characterization framework by applying it on a subset of a larger scale inverted model. This allows to further spare computational time by reducing the dimension of the model space.

## 6 CONCLUSION

We presented a downscaling approach, also referred to as inverse homogenization, to quantify stratigraphic and structural uncertainty in the context of FWI. The objective of the downscaling is to derive an ensemble of small-scale models of the earth elastic properties, knowing an estimation of its macroscale effective properties, which we assume can be derived by homogenizing a FWI result. With a synthetic application on a faulted fine-layered velocity model, in which uncertainty about fault model properties and seismic P-wave velocity is addressed with Bayesian inference, we assessed the feasibility and the robustness of the methodology and demonstrated its benefits for practical, real-world applications. In particular, the proposed uncertainty quantification framework allows for 'localized' inference: it can be locally applied to a large-scale FWI data set to downscale a relevant geologic target assumed to have an impact on decision making. Furthermore, it allows the integration of prior geological knowledge and makes use of simplified forward physics (the elastic homogenization), which accelerates the forward predictions and makes the methodology computationally tractable, even with uninformative priors. With ongoing advancements in this research, the methodology will be suited to all applications requiring a detailed characterization of the stratigraphy and fault structure at the small-scale, such as well placement optimization, sedimentary basin analysis, reservoir monitoring or seismic hazard assessment.

## ACKNOWLEDGMENTS

The authors would like to thank Y. Capdeville for providing the Fast Fourier Homogenization code employed for forward modelling. We would also like to thank Y. Capdeville, T. Bodin, L. Métivier and P. Thore for constructive discussions about the subject. This work was performed in the frame of the RING project (http://ring.geo ressources.univ-lorraine.fr/) at Université de Lorraine. We would like to thank for their support the industrial and academic sponsors of the RING-GOCAD Consortium managed by ASGA. The High Performance Computing resources were partially provided by the EXPLOR centre hosted by the University de Lorraine.

## SUPPORTING INFORMATION

Supplementary data are available at *GJIRAS* online.
*suppl_data*

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

## DATA AVAILABILITY

This study uses synthetic data only. All the codes used will be shared on reasonable request to the corresponding author.

## REFERENCES

Afanasiev, M., Boehm, C., May, D. & Fichtner, A., 2016. Using effective medium theory to better constrain full waveform inversion, in *78th EAGE Conference and Exhibition 2016*, Vol. **2016**, pp. 1–3, European Association of Geoscientists & Engineers(online).

Aissat, R., Pryet, A., Saltel, M. & Dupuy, A., 2023. Comparison of different pilot point parameterization strategies when measurements are unevenly distributed in space, *Hydrogeol. J.*, **31**(8), 2381–2400.

Alcolea, A., Carrera, J. & Medina, A., 2006. Pilot points method incorporating prior information for solving the groundwater flow inverse problem, *Adv. Water Resour.*, **29**(11), 1678–1689.

Alghuraybi, A., Bell, R.E., Jackson, C.A., Sim, M. & Jin, S., 2024. Full-waveform inversion as a tool to predict fault zone acoustic properties, *Geophys. Prospect.*, **72**, 3168–3183. doi:

Aydin, O. & Caers, J.K., 2017. Quantifying structural uncertainty on fault networks using a marked point process within a Bayesian framework, *Tectonophysics*, **712**, 101–124.

Backus, G.E., 1962. Long-wave elastic anisotropy produced by horizontal layering, *J. geophys. Res.*, **67**(11), 4427–4440.

Bodin, T., Capdeville, Y., Romanowicz, B. & Montagner, J.-P., 2015. Interpreting radial anisotropy in global and regional tomographic models, in *The Earth's Heterogeneous Mantle: A Geophysical, Geodynamical, and Geochemical Perspective*, pp. 105–144, eds Khan, A. & Deschamps, F., Springer.

Brocher, T.M., 2005. Empirical relations between elastic wavespeeds and density in the earth's crust, *Bull. seism. Soc. Am.*, **95**(6), 2081–2092.

Brooks, S., Gelman, A., Jones, G. & Meng, X.-L., 2011. *Handbook of Markov chain Monte Carlo*, CRC Press.

Brooks, S.P. & Gelman, A., 1998. General methods for monitoring convergence of iterative simulations, *J. Comput. Graph. Stat.*, **7**(4), 434–455.

Browaeys, J.T. & Chevrot, S., 2004. Decomposition of the elastic tensor and geophysical applications, *Geophys. J. Int.*, **159**(2), 667–678.

Brown, A.R., 2011. *Interpretation of Three-Dimensional Seismic Data*, American Association of Petroleum Geologists.

Bui-Thanh, T., Ghattas, O., Martin, J. & Stadler, G., 2013. A computational framework for infinite-dimensional Bayesian inverse problems Part I: the
linearized case, with application to global seismic inversion, *SIAM J. Sci. Comput.*, **35**(6), A2494–A2523.

Capdeville, Y. & Marigo, J.-J., 2007. Second order homogenization of the elastic wave equation for non-periodic layered media, *Geophys. J. Int.*, **170**(2), 823–838.

Capdeville, Y. & Métivier, L., 2018. Elastic full waveform inversion based on the homogenization method: theoretical framework and 2-D numerical illustrations, *Geophys. J. Int.*, **213**(2), 1093–1112.

Capdeville, Y., Guillot, L. & Marigo, J.-J., 2010a. 1-D non-periodic homogenization for the seismic wave equation, *Geophys. J. Int.*, **181**(2), 897–910.

Capdeville, Y., Guillot, L. & Marigo, J.-J., 2010b. 2-D non-periodic homogenization to upscale elastic media for P–SV waves, *Geophys. J. Int.*, **182**(2), 903–922.

Capdeville, Y., Stutzmann, E., Wang, N. & Montagner, J.-P., 2013. Residual homogenization for seismic forward and inverse problems in layered media, *Geophys. J. Int.*, **194**(1), 470–487.

Capdeville, Y., Zhao, M. & Cupillard, P., 2015. Fast Fourier homogenization for elastic wave propagation in complex media, *Wave Motion*, **54**, 170–186.

Capdeville, Y., Cupillard, P. & Singh, S., 2020. An introduction to the two-scale homogenization method for seismology, in *Advances in Geophysics*, Vol. **61**, pp. 217–306, eds Moseley, B. & Krischer, L., Elsevier.

Carr, J.R., 2002. *Data Visualization in the Geosciences*, Prentice-Hall Upper Saddle River, NJ.

Castagna, J.P., Batzle, M.L. & Eastwood, R.L., 1985. Relationships between compressional-wave and shear-wave velocities in clastic silicate rocks, *Geophysics*, **50**(4), 571–581.

Cherpeau, N. & Caumon, G., 2015. Stochastic structural modelling in sparse data situations, *Petrol. Geosci.*, **21**, 233–247.

Cherpeau, N., Caumon, G. & Lévy, B., 2010. Stochastic simulations of fault networks in 3D structural modeling, *C. R. Géosci.*, **342**(9), 687–694.

Cupillard, P. & Capdeville, Y., 2018. Non-periodic homogenization of 3-D elastic media for the seismic wave equation, *Geophys. J. Int.*, **213**(2), 983–1001.

Curtis, A. & Lomax, A., 2001. Prior information, sampling distributions, and the curse of dimensionality, *Geophysics*, **66**(2), 372–378.

Dahlke, T., Biondi, B. & Clapp, R., 2020. Applied 3D salt body reconstruction using shape optimization with level sets, *Geophysics*, **85**(5), R437–R446.

de Marsily, G., Lavedan, G., Boucher, M. & Fasamino, G., 1984. Interpretation of interference tests in a well field using geostatistical techniques to fit the permeability distribution in a reservoir model, in *Geostatistics for Natural Resources Characterization (NATO Advanced Study Institute)*, pp. 831–849, eds Verly, G., David, M., G. Journel, A. & Marechal, A., Springer.

Dimmen, V., Rotevatn, A. & Lecomte, I., 2023. Imaging of small-scale faults in seismic reflection data: insights from seismic modelling of faults in outcrop, *Mar. Petrol. Geol.*, **147**, 105980.

Dinh, H., Latter, T., Townsend, M., Grinde, N., Hogden, S., Robb, N., Ak- sland, M. & Bertrand, A., 2024. High-frequency FWI imaging: repurpos- ing seismic data for imaging shallow hazards, *First Break*, **42**(8), 55–63.

Dockrill, B. & Shipton, Z.K., 2010. Structural controls on leakage from a natural CO₂ geologic storage site: Central Utah, USA, *J. Struct. Geol.*, **32**(11), 1768–1782.

Engquist, B., Holst, H. & Runborg, O., 2009. Multi-scale methods for wave propagation in heterogeneous media, *Commun. Math. Sci.*, **9**, 33–56.

Fairley, J., 2009. Modeling fluid flow in a heterogeneous, fault-controlled hydrothermal system, *Geofluids*, **9**(2), 153–166.

Faleide, T.S., Braathen, A., Lecomte, I., Mulrooney, M.J., Midtkandal, I., Bugge, A.J. & Planke, S., 2021. Impacts of seismic resolution on fault interpretation: Insights from seismic modelling, *Tectonophysics*, **816**, 229008.

Fedorov, F.I., 2013. *Theory of Elastic Waves in Crystals*, Springer Science & Business Media.

Fichtner, A. & Hanasoge, S.M., 2017. Discrete wave equation upscaling, *Geophys. J. Int.*, **209**(1), 353–357.

Fichtner, A., Trampert, J., Cupillard, P., Saygin, E., Taymaz, T., Capdeville, Y. & Villasenor, A., 2013. Multiscale full waveform inversion, *Geophys. J. Int.*, **194**(1), 534–556.

Freulon, X. & de Fouquet, C., 1993. Conditioning a Gaussian model with inequalities, in *Geostatistics Tróia'92, Vol. 1*, pp. 201–212, ed. Soares, A., Springer.

Gao, K., Chung, E.T., Gibson, R.L., Jr, Fu, S. & Efendiev, Y., 2015. A numerical homogenization method for heterogeneous, anisotropic elastic media based on multiscale theory, *Geophysics*, **80**(4), D385–D401.

Gardner, G., Gardner, L. & Gregory, A., 1974. Formation velocity and density—the diagnostic basics for stratigraphic traps, *Geophysics*, **39**(6), 770–780.

Gelman, A. & Rubin, D.B., 1992. Inference from iterative simulation using multiple sequences, *Stat. Sci.*, **7**(4), 457–472.

Georgsen, F., Rœ, P., Syversveen, A.R. & Lia, O., 2012. Fault displacement modelling using 3D vector fields, *Comput. Geosci.*, **16**, 247–259.

Gercek, H., 2007. Poisson's ratio values for rocks, *Int. J. Rock Mech. Mining Sci.*, **44**(1), 1–13.

Gilks, W.R., Richardson, S. & Spiegelhalter, D., 1995. *Markov Chain Monte Carlo in Practice*, CRC Press.

Gilmore, K.A., Sahu, C.K., Benham, G.P., Neufeld, J.A. & Bickle, M.J., 2022. Leakage dynamics of fault zones: experimental and analytical study with application to CO₂ storage, *J. Fluid Mech.*, **931**, A31.

Giraud, J., Ford, M., Caumon, G., Ogarko, V., Grose, L., Martin, R. & Cupillard, P., 2024. Geologically constrained geometry inversion and null-space navigation to explore alternative geological scenarios: a case study in the western Pyrenees, *Geophys. J. Int.*, **239**(3), 1359–1379.

Godefroy, G., Caumon, G., Ford, M., Laurent, G. & Jackson, C. A.-L., 2018. A parametric fault displacement model to introduce kinematic control into modeling faults from sparse data, *Interpretation*, **6**(2), B1–B13.

Gómez-Hernánez, J.J., Sahuquillo, A. & Capilla, J., 1997. Stochastic simulation of transmissivity fields conditional to both transmissivity and piezometric data—I. theory, *J. Hydrol.*, **203**(1–4), 162–174.

Goovaerts, P., 1997. *Geostatistics for Natural Resources Evaluation*, Oxford Univ. Press.

Grana, D., de Figueiredo, L. & Mosegaard, K., 2023. Markov chain Monte Carlo for seismic facies classification, *Geophysics*, **88**(3), M131–M143.

Grechka, V., 2003. Effective media: a forward modeling view, *Geophysics*, **68**(6), 2055–2062.

gstlearn, 2023. *Geostatistics and Machine Learning Toolbox*. Available at: https://gstlearn.org>.

Guillot, L., Capdeville, Y. & Marigo, J.-J., 2010. 2-D non-periodic homogenization of the elastic wave equation: SH case, *Geophys. J. Int.*, **182**(3), 1438–1454.

Haario, H., Saksman, E. & Tamminen, J., 1999. Adaptive proposal distribution for random walk metropolis algorithm, *Comput. Stat.*, **14**, 375–395.

Hansen, T.M., Journel, A.G., Tarantola, A. & Mosegaard, K., 2006. Linear inverse Gaussian theory and geostatistics, *Geophysics*, **71**(6), R101–R111.

Hastings, W.K., 1970. Monte Carlo sampling methods using Markov chains and their applications, *Biometrika*, **57**, 97–109. doi:

Hedjazian, N., Capdeville, Y. & Bodin, T., 2021. Multiscale seismic imaging with inverse homogenization, *Geophys. J. Int.*, **226**(1), 676–691.

Herron, D.A., 2011. *First Steps in Seismic Interpretation*, Society of Exploration Geophysicists.

Holden, L., Mostad, P., Nielsen, B.F., Gjerde, J., Townsend, C. & Ottesen, S., 2003. Stochastic structural modeling, *Math. Geol.*, **35**, 899–914.

Holmes, J.J., Driscoll, N.W. & Kent, G.M., 2021. High-resolution 3D seismic imaging of fault interaction and deformation offshore San Onofre, California, *Front. Earth Sci.*, **9**, 653672.

Jordan, T.H., 2015. An effective medium theory for three-dimensional elastic heterogeneities, *Geophys. J. Int.*, **203**(2), 1343–1354.

Kadu, A., Van Leeuwen, T. & Mulder, W., 2017. Parametric level-set full-waveform inversion in the presence of salt bodies, in *SEG Technical
*Program Expanded Abstracts* 2017, pp. 1518–1522, Society of Exploration Geophysicists. Available at: https:// homepage.tudelft.nl/ a1k53/Papers/Conferences/SEG2017_Kadu.pdf

Kolbjørnsen, O., Buland, A., Hauge, R., Røe, P., Ndingwan, A.O. & Aker, E., 2020. Bayesian seismic inversion for stratigraphic horizon, lithology, and fluid prediction, *Geophysics*, **85**(3), R207–R221.

Komatitsch, D. & Tromp, J., 1999. Introduction to the spectral element method for three-dimensional seismic wave propagation, *Geophys. J. Int.*, **139**(3), 806–822.

Kumar, D. & Ali, R., 2024. Integration of FWI-derived reflectivity in seismic interpretation, *Leading Edge*, **43**(12), 836–842.

Lantuëjoul, C., 2001. *Geostatistical Simulation: Models and Algorithms, no. 1139*, Springer Science & Business Media.

Laurent, G., Caumon, G., Bouziat, A. & Jessell, M., 2013. A parametric method to model 3D displacements around faults with volumetric vector fields, *Tectonophysics*, **590**, 83–93.

Li, Y., Alkhalifah, T. & Guo, Q., 2021. Target-oriented time-lapse waveform inversion using deep learning-assisted regularization, *Geophysics*, **86**(4), R485–R495.

Liu, Q. & Peter, D., 2019. Square-root variable metric based elastic full-waveform inversion—Part 2: Uncertainty estimation, *Geophys. J. Int.*, **218**(2), 1100–1120.

Manzocchi, T., Childs, C. & Walsh, J., 2010. Faults and fault properties in hydrocarbon flow models, *Geofluids*, **10**(1–2), 94–113.

Matheron, G., 1973. The intrinsic random functions and their applications, in *Advances in Applied Probability (Vol. 5, Issue 3)*, pp. 439–468, ed. Taylor, P., Cambridge Univ. Press

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**(6), 1087–1092.

Mosser, L., Dubrule, O. & Blunt, M.J., 2020. Stochastic seismic waveform inversion using generative adversarial networks as a geological prior, *Math. Geosci.*, **52**(1), 53–79.

Muir, J.B. & Tsai, V.C., 2020. Geometric and level set tomography using ensemble Kalman inversion, *Geophys. J. Int.*, **220**(2), 967–980.

Muir, J.B., Clayton, R.W., Tsai, V.C. & Brissaud, Q., 2022. Parsimonious velocity inversion applied to the Los Angeles Basin, CA, *J. geophys. Res.: Solid Earth*, **127**(2), e2021JB023103.

Neal, R.M., 1993. *Probabilistic inference using Markov chain Monte Carlo methods*. Available at: https://www.cs.columbia.edu/~blei/fogm/2025F/r eadings/Neal1993.pdf

Plessix, R.-E., 2006. A review of the adjoint-state method for computing the gradient of a functional with geophysical applications, *Geophys. J. Int.*, **167**(2), 495–503.

Pyrcz, M.J. & Deutsch, C.V., 2014. *Geostatistical Reservoir Modeling*, Oxford Univ. Press.

Richards, F.L., Richardson, N.J., Bond, C.E. & Cowgill, M., 2015. Interpretational variability of structural traps: implications for exploration risk and volume uncertainty, *Geol. Soc. (Special Publications)*, **421**(1), 7–27.

Sambridge, M., Gallagher, K., Jackson, A. & Rickwood, P., 2006. Trans-dimensional inverse problems, model comparison and the evidence, *Geophys. J. Int.*, **167**(2), 528–542.

Sánchez-Palencia, E., 1980. *Non-Homogeneous Media and Vibration Theory* (Lecture Note in Physics, Springer-Verlag), Vol. **320**, pp. 57–65, Springer.

Santos, T., Bodin, T., Soulez, F., Ricard, Y. & Capdeville, Y., 2024. Refining tomography with generative neural networks trained from geodynamics, *Geophys. J. Int.*, **238**(3), 1676–1695.

Sheriff, R., 1985. *Aspects of Seismic Resolution: Chapter 1*. *AAPG Memoir*, Vol. **39**, American Association of Petroleum Geologists.

Sheriff, R.E. & Geldart, L.P., 1995. *Exploration Seismology*, Cambridge Univ. Press.

Singh, S., Tsvankin, I. & Naeini, E.Z., 2018. Bayesian framework for elastic full-waveform inversion with facies information, *Leading Edge*, **37**(12), 924–931.

Tarantola, A., 2005. *Inverse Problem Theory and Methods for Model Parameter Estimation*, SIAM.

Taty Moukati, F., Stoica, R.S., Bonneau, F., Wu, X. & Caumon, G., 2025. From fault likelihood to fault networks: stochastic seismic interpretation through a marked point process with interactions, *Math. Geosci.*, **57**(1), 115–151.

Tsai, V.C., Huber, C. & Dalton, C.A., 2023. Towards the geological parametrization of seismic tomography, *Geophys. J. Int.*, **234**(2), 1447–1462.

Virieux, J. & Operto, S., 2009. An overview of full-waveform inversion in exploration geophysics, *Geophysics*, **74**(6), WCC1–WCC26.

Visser, G., Guo, P. & Saygin, E., 2019. Bayesian transdimensional seismic full-waveform inversion with a dipping layer parameterization, *Geophysics*, **84**(6), R845–R858.

, Kitanidis, P & Caers, J., 2022. Hierarchical Bayesian inversion of global variables and large-scale spatial fields, *Water Resources Res.*, **58**(5), e2021WR031610.

Wang, F. & Alkhalifah, T., 2024. Geological and well prior assisted full waveform inversion using conditional diffusion models, in *85th EAGE Annual Conference and Exhibition (including the Workshop Programme)*, Oslo, Norway, Vol. **2024**, pp. 1–5, European Association of Geoscientists & Engineers.

Warner, M. et al., 2013. Anisotropic 3D full-waveform inversion, *Geophysics*, **78**(2), R59–R80.

Wernicke, B., 1995. Low-angle normal faults and seismicity: A review, *J. geophys. Res.: Solid Earth*, **100**(B10), 20 159–20 174.

Yielding, G., 2015. Trapping of buoyant fluids in fault-bound structures, *Geol. Soc. London, Special Publ.*, **421**(1), 29–39.

Zhang, Z. & Alkhalifah, T., 2022. Regularized elastic full-waveform inversion using deep learning, in *Advances in Subsurface Data Analytics*, pp. 219–250, eds Bhattacharya, S. & Di, H., Elsevier.

Zhang, Z., Alkhalifah, T., Naeini, E.Z. & Sun, B., 2018. Multiparameter elastic full waveform inversion with facies-based constraints, *Geophys. J. Int.*, **213**(3), 2112–2127.

Zhang, Z., Deng, Y., Qiu, H., Peng, Z. & Liu-Zeng, J., 2022. High-resolution imaging of fault zone structure along the creeping section of the Haiyuan fault, NE Tibet, from data recorded by dense seismic arrays, *J. geophys. Res.: Solid Earth*, **127**(9), e2022JB024468.

Zhao, X. & Curtis, A., 2024. Physically structured variational inference for Bayesian full waveform inversion, *J. geophys. Res.: Solid Earth*, **129**(11), e2024JB029557.

Zhu, H., Li, S., Fomel, S., Stadler, G. & Ghattas, O., 2016. A Bayesian approach to estimate uncertainty for full-waveform inversion using a *priori* information from depth migration, *Geophysics*, **81**(5), R307–R323.