11institutetext: Centre for Astrochemical Studies, Max Planck Institute for
Extraterrestrial Physics, Giessenbachstraße 1, 85748 Garching, Germany
11email<EMAIL_ADDRESS>22institutetext: LERMA, CY Cergy Paris
Université, Observatoire de Paris, PSL University, Sorbonne Université, CNRS,
LERMA, F-95000 Cergy, France

# Role of NH3 binding energy in the early evolution of protostellar cores

S. Kakkenpara Suresh 1122 O. Sipilä 22 P. Caselli 22 F. Dulieu 11

###### Abstract

Context. NH3(ammonia) plays a critical role in the chemistry of star and
planet formation, yet uncertainties in its binding energy (BE) values
complicate accurate estimates of its abundances. Recent research suggests a
multi-binding energy approach, challenging the previous single-value notion.

Aims. In this work, we use different values of NH3 binding energy to examine
its effects on the NH3 abundances and, consequently, in the early evolution of
protostellar cores.

Methods. Using a gas-grain chemical network, we systematically vary the values
of NH3 binding energies in a model Class 0 protostellar core and study the
effects of these binding energies on the NH3 abundances.

Results. Our simulations indicate that abundance profiles of NH3 are highly
sensitive to the binding energy used, particularly in the warmer inner regions
of the core. Higher binding energies lead to lower gas-phase NH3 abundances,
while lower values of binding energy have the opposite effect. Furthermore,
this BE-dependent abundance variation of NH3 significantly affects the
formation pathways and abundances of key species such as HNC, HCN, and CN. Our
tests also reveal that the size variation of the emitting region due to
binding energy becomes discernible only with beam sizes of 10 arcsec or less.

Conclusions. These findings underscore the importance of considering a range
of binding energies in astrochemical models and highlight the need for higher
resolution observations to better understand the subtleties of molecular cloud
chemistry and star formation processes.

###### Key Words.:

Astrochemistry – ISM: molecules – Radiative transfer – ISM: abundances –
Stars: protostars

##  1 Introduction

Since its detection (Cheung et al. 1968), NH3 has been ubiquitously observed
in a variety of environments like molecular clouds (Irvine et al. 1987),
prestellar cores (Crapsi et al. 2007), the galactic centre (Winnewisser et al.
1979), diffuse clouds (Liszt et al. 2006), galaxies (Sandqvist et al. 2017;
Gorski et al. 2018), star-forming regions (Fehér et al. 2022), comets (Poch et
al. 2020), and planet-forming disks (Salinas et al. 2016). It is one of the
six major molecules found in interstellar ices in the solid form (Boogert et
al. 2015a). It serves as an important tracer in the interiors of dense,
starless cores where common tracers like CO and CS are depleted from the gas
phase onto dust grains (Caselli et al. 1999; Tafalla et al. 2002) due to
temperatures as low as $\sim$ 6K (Crapsi et al. 2007; Pagani et al. 2007) and
number densities between 104 and 106 cm-3 (Keto & Caselli 2010). Although NH3
is affected by freeze-out (Caselli et al. 2022), NH3 is selectively tracing
dense material, as it is not abundant in the molecular cloud surrounding dense
cores (unlike CO and CS).

An important parameter that influences the grain surface abundance of NH3 (and
other species) is the binding energy. It serves as a measure of the strength
of the interaction between the species and an adsorbate or grain surface.
Additionally, once the species is adsorbed onto dust grains, the binding
energy plays a crucial role in governing both its diffusion across and
residence time on the grain surface. These factors, in turn, significantly
influence the physical and chemical conditions governing the evolution of
star-forming regions. For example, in the vicinity of a low-mass protostar,
the binding energy determines the distance from the central protostar at which
sublimation fronts emerge due to the desorption of volatile species from ice
grain surfaces. Furthermore, the location of these sublimation fronts, also
known as ”snowlines” in protoplanetary disks, directly shapes the composition
of the planets forming within protoplanetary disks. Hence, precise knowledge
of the binding energies is critical for an accurate interpretation of
observational data and reliable prediction of abundances through simulations.

Earlier studies, such as Hama & Watanabe (2013), Penteado et al. (2017) and
Wakelam et al. (2017), provided a singular binding energy value for NH3 on
different types of water ice. However, recent research has challenged this
notion, suggesting that NH3 exhibits a distribution of binding energies.
Experiments conducted by He et al. (2016) revealed a surface coverage-
dependent distribution of binding energies for NH3 at sub-monolayer coverages.
Additionally, Kakkenpara Suresh et al. (2024) observed a distribution of
binding energies on crystalline ice (3780K – 4080K) and compact amorphous
solid water ice (3630K – 5280K) at monolayer coverages. Theoretical
calculations by Ferrero et al. (2020), Tinacci et al. (2022), and Germain et
al. (2022) have also supported this conclusion. The notion of a binding energy
distribution is logical given that the structure of ice on dust grains is
amorphous, resulting in the formation of diverse and unique adsorption sites.
This idea is further substantiated by Bovolenta et al. (2020), who obtained a
Gaussian-like distribution of binding energies for hydrogen fluoride on
amorphous solid water (ASW).

In their numerical study, Grassi et al. (2020) have illustrated that employing
a binding energy distribution allows molecules to occupy higher energy binding
sites, thereby increasing their residence time on the grain surface. This
prolonged residence time enhances their availability to react with other
molecules, even at dust temperatures that conventionally exhibit limited or no
reactivity. Building upon their findings and our previous work in Kakkenpara
Suresh et al. (2024), the present study evaluates the influence of
incorporating a range of binding energy values for NH3 in astrochemical
models, in particular the impact of binding energy on the abundance
distribution of NH3. We also investigate the possible effect that NH3
abundance variations may have on the distributions of other molecules
chemically linked to NH3. We do not expect the measured values of binding
energy to have an impact on the chemistry in prestellar cores due to their low
temperatures. Hence, in the present work, we focus our attention on
protostellar cores. We focus on the very early stages of the evolution of a
protostellar core, and use a model of the well-studied class 0 protostellar
core, IRAS 16293-2422, as a template for our simulations.

This work is organised in the following way. In Section 2, we outline the
physical and chemical model of the protostellar core used. Section 3 presents
our results, and Section 4 delves into the implications of these findings
along with potential avenues for future research. A concise summary of our
work is provided in Section 5.

##  2 Model

###  2.1 Physical Model

The core is characterized using the model by Crimier et al. (2010) (Fig. 1)
where the radial density profile follows as n(H2) $\propto$ r-1.8 within 6900
AU from the centre, where n(H2) is the H2 number density. The source exhibits
a pronounced temperature and density gradient. In this model, the temperature
of the gas rises significantly towards the interior due to both gas
compression during collapse and radiation emitted by the protostar. We assume
that the dust temperature and the gas kinetic temperature are equivalent. The
analytical approach involves dividing the core into concentric shells, and the
final results are combined radially. A two-phase (gas + ice) chemical model
with a dust grain radius = 0.1 $\mu$m is employed, where the entire ice layer
covering the grain is available for desorption.

The abundance profiles are obtained in two steps. First, the initial
conditions corresponding to the parent cloud are obtained by running a single-
point simulation with Tdust = Tgas = 10K, n(H2) = 104 cm-3, grain radius =
0.1$\mu$m, cosmic-ray ionization rate, $\zeta$ = 1.3$\times$e-17 s-1, and
visual extinction AV = 10 mag. Adsorption, desorption, and photodesorption are
also considered, and the simulation is allowed to proceed until a time, t =
106 years, consistent with the previous work of Brünken et al. (2014) and
Harju et al. (2017) on the source. At this stage, the abundances of all the
species are extracted. These are then used as the initial abundances to trace
the evolution of the protostellar core. Finally, the abundances of this core
are extracted after an evolutionary time of t = 104 years. The choice of time
is arbitrary; sub-structure (e.g., a protostellar disk) is expected to form
within the inner regions in a timescale of 104 \- 105 years which is not taken
into account by the present static physical model. Hence, we chose an early
time step (104 years) to obtain the abundances, giving an estimate of the
initial chemical conditions of the forming small-scale structures.

Figure 1: (Top) The H2 number density and (bottom) the temperature
distribution assumed to model IRAS 16293-2422.

###  2.2 Chemical Model

The chemical evolution of the core is monitored by the gas-grain chemical
code, pyRate, discussed in Sipilä (2012); Sipilä et al. (2010). The model is
pseudo-time dependent, i.e., we track the chemical evolution assuming a static
core. pyRate employs the rate equation method for computing the molecular
abundances. The gas-phase network is constructed based on the
kida.uva.2014111https://kida.astrochem-tools.org/networks.html network
(Wakelam et al. 2015) with modifications to incorporate deuterated species and
spin state chemistry (Sipilä et al. 2015a, b). Here we employ a large chemical
network that contains a combined total of over 75,000 reactions in the gas
phase and on grain surfaces, with which we simulate the chemistry of NH3 and
its deuterated forms. In addition, the KIDA network is also used to study the
effect of proton transfer reactions on the abundances of NH3 (Section 4.2).

Initially, the species are assumed to be atomic, except for H2 and HD, and the
initial ortho/para H2 ratio is set to 10-3 consistent with a spin temperature,
$T_{spin}$ $\sim$ 20K (Brünken et al. 2014; Crabtree et al. 2011). Choosing
this specific value for the ratio of ortho/para H2 corresponds to the
assumption that the spin-state ratio has had the time to undergo
thermalization prior to the formation of the core. The initial abundances of
the species are provided in Table 1. The binding energies for NH3 on water ice
surfaces are sourced from Kakkenpara Suresh et al. (2024). Binding energies of
various other species on water ice are obtained from Garrod & Herbst (2006)
and Sipilä (2012). The simulated abundances of NH3 represent the combined
abundances of the ortho and para forms of the molecule.

Species | Abundance
---|---
H2 | 5.00 $\times$ 10-1
He | 9.00 $\times$ 10-2
HD | 1.60 $\times$ 10-5
C+ | 1.20 $\times$ 10-4
N | 7.60 $\times$ 10-5
O | 2.56 $\times$ 10-4
S+ | 8.00 $\times$ 10-8
Si+ | 8.00 $\times$ 10-9
Na+ | 2.00 $\times$ 10-9
Mg+ | 7.00 $\times$ 10-9
Fe+ | 3.00 $\times$ 10-9
P+ | 2.00 $\times$ 10-10
Cl+ | 1.00 $\times$ 10-9
Table 1: Initial chemical abundances with respect to total H nuclei, nH

##  3 Results

###  3.1 Abundance Profiles

Motivated by the experimental results of Kakkenpara Suresh et al. (2024), we
investigated how variations in NH3 binding energy impact the chemical
composition of a protostellar core. For this purpose, we conducted a series of
simulations where we systematically varied only the binding energy of NH3 in
each model and traced the resulting chemical evolution with time within the
core. We extracted the abundances of NH3 at t = 104 years representing the
protostar at its early stages of formation.

Figure 2 presents the abundances of NH3 at this time in the gas-phase (solid
lines) and on grain surfaces (dashed lines) with respect to the radial
distance from the centre to the outer edges of the core. The gas-phase
abundance of NH3 decreases moving radially inwards toward the centre of the
core before increasing rapidly by several orders of magnitude. The abundance
profile shows clear differentiation based on the value of binding energy
employed. The region of high NH3 gas-phase abundance appears closer to the
centre with the radius of the desorption zone varying between 150 and 300 AU
with increasing binding energy value used. Investigation of the reaction rates
at this time step revealed that the primary contributor influencing NH3 gas-
phase abundances is thermal desorption from dust grains. As the binding energy
increases, NH3 remains on the grains until a higher temperature is reached at
which point it acquires sufficient thermal energy for desorption.

Figure 2: Radial abundances of NH3 in gas phase (solid lines) and on grain
surfaces (dashed lines) at 104 years. The binding energies of NH3 used in each
model (colours) are displayed in the lower right corner.

(a)

(b)

Figure 3: Variation in the (a) gas-phase and (b) grain abundances of NH3, HNC,
HCN and CN with binding energy. The colour scheme for the lines follows the
same as in Fig. 2

Notable related effects are only seen in species like HNC, CN, HCN etc (Fig.
3). We find that their abundance profiles vary within the same spatial zone
where NH3 abundances vary. The effects are directly tied via chemical
reactions to the variations in the NH3 abundances. This effect is further
discussed in section 4.1. Furthermore, we analyzed the relative abundances of
key volatile compounds found in ices, namely NH3, CO, CO2, CH4, and CH3OH with
respect to H2O. We compared our predictions with the reported values in
Boogert et al. (2015b), but direct comparisons are challenging due to
differences in the sources and uncertainties related to the parameters
employed for calculating ice abundances in the cited study. Further details of
these comparisons are available in appendix A.

###  3.2 Column Density maps of p-NH3

To assess the potential observational impact of the NH3 abundance variations,
we have simulated NH3 column density maps individually for each BE value. Here
we consider p-NH3 only so that we can compare the column density maps against
the simulated emission of the (1,1) transition of this molecule (see Sect.
3.3). To generate the column density maps, we employ the radial abundances of
p-NH3 obtained in the previous section as input, which are then interpolated
to create a two-dimensional map of the cloud core. The abundances are
convolved to a beam size of 2” for the object situated at 120 pc to simulate
interferometric observations (the distance to IRAS 16293-2422, although this
can be applied to other similar sources by simply changing the distance
value). The radial distribution of column density, denoted as N, is calculated
using the formula:

${\rm N=\sum_{i}p\times X_{i}\times n(H_{2})}$

(a)

(b)

(c)

Figure 4: Column density map of p-NH3 without the envelope (see text) for
binding energy (a) 3870K (b) 4080K and (c) 5280K. The angular size simulated
by each pixel for each map is given in the top left corner of Fig. 4(c).

(a)

(b)

(c)

Figure 5: Integrated intensity map of p-NH3 (1,1) in the core for binding
energy (a) 3870K (b) 4080K and (c) 5280K with envelope

Here, $\Sigma$ represents the summation across i elements of the core model, p
is the path length through an element i in cm, Xi is the abundance of the
species (here, p-NH3) in i with respect to H2, and n(H2) is the volume density
of i in mol cm-3. Figure 4 shows the 2D maps of column densities of p-NH3
within a 300 AU radius where the binding energy-dependent variations are
apparent. The radius of the desorption zone decreases as the binding energy
increases. The column density as displayed in Fig. 4 traces the part of the
distribution that could only be observed with specific NH3 lines, for example,
the (1,1) inversion line. The column density maps also indicate the extent of
desorbed NH3 that will be available in the gas phase.

###  3.3 Radiative transfer studies of NH3 (1,1) transition

Subsequently, one-dimensional radiative transfer modelling of the 23 GHz (1,1)
rotational-inversion transition of p-NH3 is carried out to verify the
observability of the binding energy-dependent variation in the radius of the
desorption zone using the non-LTE radiative transfer code LOC (Juvela 2020).
The first step in this analysis is to simulate the cloud that contains the
protostellar core and its envelope. Given that our source is a Class 0 object,
the core represents a deeply embedded source surrounded by an envelope. The
methodology employed to determine NH3 abundances within the envelope closely
follows the two-step approach outlined for the core in Section 2.1. First, we
obtain the abundances corresponding to the parent cloud by running a single-
point simulation for t = 106 years under the same conditions as the parent
cloud, similar to that described in Section 2.1. After extracting the NH3
abundances at this stage, we run a second single-point simulation, this time
using the specific physical conditions for the envelope. These conditions
include an envelope thickness set to 0.1 pc (5 $\times$ 1017 cm) assuming a
visual extinction, AV, = 5 mag, n(H2) = 104 m-3 and T = 10 K. The modelling is
done taking into consideration the 18 hyperfine components of p-NH3. The
velocity resolution is set to $\sim$ 0.835 km/s such that it is high enough to
distinguish the hyperfine components. The spectra are obtained for 1000 lines
of sight across the radius of the whole object. Each spectrum is then
convolved to a synthesized beam of 2”. The intensity of the transition from
these spectra is integrated and interpolated to create a two-dimensional
intensity map (Fig. 5) assuming spherical symmetry.

The collisional and radiative rate coefficients for the p-NH3 line simulations
were taken from the LAMDA database (Schöier et al. 2005). The provided
coefficients (collisional data from Danby et al. (1988)) do not resolve the
hyperfine structure, and hence we assumed that the coefficients for the
individual hyperfine components were distributed according to LTE. To
investigate the potential effect of hyperfine-resolved collisional rate
coefficients on our results, we ran another series of line simulations
adopting instead the set of collisional rate coefficients presented recently
by Loreau et al. (2023). The hyperfine-resolved radiative transition
frequencies and Einstein A coefficients were derived from data in the CDMS
(Endres et al. 2016). These two approaches led only to small differences in
simulated lines, and in what follows we present the results of simulations
carried out using the data originating in LAMDA.

In each model, we observe that the peak intensity is centrally concentrated
and declines going outwards. The location and size of this zone are in good
agreement with the simulations of the column densities (Fig. 4), despite
optical thickness effects that were missed by calculating the column density
directly from the simulated abundances. Predominantly, the intensity
originates from the inner zone due to elevated NH3 gas-phase abundance, as
indicated by the abundance profiles (Fig. 2). Moreover, the size of the
emitting region diminishes with higher binding energy, attributed to a
reduction in gas-phase NH3 concentrations associated with increasing binding
energy. Similar results are obtained for radiative transfer models using the
(2,2) and (3,3) lines.

Convolutions of the spectra for the (1,1) transitions at larger beam sizes are
also conducted to evaluate the observability of this variation in the size of
the emitting region. However, these effects become discernible only with beam
sizes of 10 arcsec or smaller. This highlights the need for higher-resolution
observations to detect such subtleties in the star-formation process. When
convolved with a 6 arcsec beam, similar to in the observations by Mundy et al.
(1990), we obtain brightness temperatures lower than but within a factor of
two compared to their reported value of 16.5 K for the NH3 (1,1) emission
towards the brightest regions of their source - IRAS 16293-2422. They report
NH3 emissions arising from a ring-like region of 3000-4000 AU with a
resolution of 6 arcsec (960 AU). In contrast, the BE-related effects in the
model used in the present work originate from a region smaller than 300 AU,
indicating that higher resolution observations will be necessary to draw
important conclusions regarding effects due to NH3 binding energy towards this
source. Additionally, we ran radiative transfer models assuming higher
velocity resolution (0.15 km/s) for the lines based on observations by Crapsi
et al. (2007) using the Very Large Array (VLA) and found no difference in our
results. Using the higher resolution requires more computational time but does
not impact the results, so we decided to retain the lower velocity resolution
for our models.

##  4 Discussion

###  4.1 Influence of NH3 binding energy on abundance profiles of chemically
related species

Interestingly, notable effects are observed on other species and their
abundance profiles within the same zone where NH3 abundances vary. A few
examples are demonstrated in Fig. 3. The abundance profiles of these species
vary depending on the chosen binding energy of NH3. To understand this, let us
consider the example of HNC (Fig. 6). In the inner regions of the model with
binding energy = 3780K, HNC is efficiently produced through the reaction
between HCNH+ and NH3, yielding HCN and HNC. HCNH+ originates from HCN through
the interaction with H${}_{3}^{+}$ and can revert to HCN via reactions between
HCNH+ and H2CO (and HCNH+ \+ NH3). Hence, the connection between HCN and HNC
follows the sequence HCN $\rightarrow$ HCNH+ $\rightarrow$ HNC, where the
first step involves H${}_{3}^{+}$, and the second step involves NH3.

Figure 6: Main formation and destruction pathways of HNC as predicted by the
chemical model for the two extreme binding energy (BE) values. In the lower
figure (BE = 5280 K), the dashed arrow represents an alternative pathway for
the formation of HNC through H-abstraction from H2CN. This pathway becomes
prominent when the gas-phase abundance of NH3 is low, thereby inhibiting the
formation via the reaction between NH3 and HCNH+.

In contrast, in the model with binding energy = 5280K, the gas phase abundance
of NH3 is low, which effectively closes off the HCNH+ \+ NH3 route. HNC is
being created instead via H-abstraction from H2CN (dashed arrow in Fig. 6),
but at approximately $\sim$10$\%$ of the rate of the NH3 route of the previous
model; HNC formation is severely inhibited by the lack of available NH3. This
can be seen in Fig. 3 where the shape of the HNC abundance profile follows
that of NH3. Simultaneously, the rate of the reaction HCNH+ \+ H2CO
$\rightarrow$ HCN + H2COH+ increases (as the HCNH+ destruction channel with
NH3 is absent). These findings support the importance of further
investigations into the role of NH3 binding energy in the abundances of
various species in the protostellar core.

###  4.2 Effect of Proton Transfer reactions on NH3 abundances

Taquet et al. (2016) demonstrated that methanol abundances are significantly
enhanced through proton transfer reactions involving NH3, as compared to
simulations lacking NH3. Subsequently, we introduced the following reaction

CH_3OH_2^+ + NH_3 ⟶CH_3OH_gas

derived from the network outlined in their research, which contributes to the
synthesis of dimethyl ether or methyl formate. We then ran a new chemical
simulation with this reaction included, wherein only NH3 binding energies were
varied to investigate potential effects on methanol abundances. For this test,
we used the KIDA network (Wakelam et al. 2015) in place of our full deuterium
and spin-state containing networks, for two main reasons: 1) to check if our
overall conclusions on the binding energy-dependent NH3 desorption region
remain unaffected regardless of the chemical network used (i.e., that the
presence of deuterium and spin states does not affect the conclusions to a
significant degree); and 2) to search for potential effects on some complex
organic molecules (COMs) that are not included in our fiducial chemical
network (for example, CH3OCH3, CH3CH2OH, etc). The grain-surface network used
in this test is the same as our fiducial one, but with deuterium and spin
states removed.

In the investigation of Taquet et al. (2016), protonation by NH3 extended the
survival time of methanol from 104 years to 105 years. In our study, methanol
remained for up to 106 years even without their reaction integrated into the
KIDA network. Upon inclusion of their reaction, our simulations validated
their results, showing an increase in methanol abundances of up to two orders
of magnitude for a specific binding energy value. However, this enhancement
became noticeable only beyond 106 years. Analysis of abundances for different
binding energy values at a specific time step revealed marginal variations in
methanol abundances, even at 106 years, where the most substantial difference
emerged between simulations with and without the proton-transfer reaction. The
discrepancy in the evolutionary time in our model versus that of Taquet et al.
(2016) regarding when the effect on CH3OH becomes apparent is very likely due
to different simulation parameters. We have not attempted to duplicate their
results. Nevertheless, the comparison indicates that proton-transfer reactions
can be very important and should be explored in more detail in later
simulations.

##  5 Conclusions

NH3 is an important molecule widely observed in various astronomical sources
and on dust grain surfaces. However, uncertainties pertaining to its binding
energy values prevent an accurate determination of its abundance and role in
the chemistry of star and planet formation. The binding energy is a key
parameter determining the abundance and chemistry of a species in the ISM.
Recent studies have challenged previous notions of a unique binding energy
value for NH3 proposing instead a multi-binding energy approach.

Drawing on previously established experimental work by Kakkenpara Suresh et
al. (2024), we incorporated multiple NH3 binding energy values derived from
these studies into gas-grain chemical networks to examine their impact on NH3
abundance within a protostellar core. We conducted multiple simulations,
systematically varying only the NH3 binding energy, to discern its influence
on the abundance profiles of NH3 and other species. Our findings reveal a
distinct dependence of NH3 abundance profiles on the binding energy employed,
particularly in the inner warm regions of the model. This variability extends
its influence on other key species, including HCN, HNC, CN, wherein the NH3
abundance dictates the preferred pathway for their formation. On the contrary,
in proton transfer reactions involving NH3, expected to enhance methanol
formation, the abundance variation of NH3 due to binding energy does not
appear to be a significant contributing factor.

Simulation of column density maps for p-NH3 reveals that the size of the
desorption region diminishes as the binding energy increases. These findings
align with radiative transfer studies on the (1,1) inversion line of NH3,
where we added an envelope to our physical model to examine absorption and
emission effects. In these studies, we observed that the peak intensity is
centrally concentrated and decreases outward. Additionally, the intensity
diminishes with higher binding energy due to a reduction in gas-phase NH3. Our
results highlight the importance of considering diverse binding energies in
astrochemical models, providing a refined understanding of molecular cloud
chemistry and star formation processes.

For future work, it is crucial to refine our understanding by further
exploring the impact of varied binding energies on complex organic molecules
(COMs) within astrochemical models. Additionally, investigations into the
spatial distribution and temporal evolution of species influenced by multi-
binding energy approaches would contribute valuable insights. The present
results apply at very early times in the core evolution when the substructure
is still absent and, hence, probe the very early stages of forming
protostellar systems. Consequently, there is a critical need to enhance the
physical model, incorporating more precise representations of observed
substructures within protostellar cores. This improvement is essential for
providing a more accurate and detailed description of the underlying phenomena
in such sources.

###### Acknowledgements.

The authors acknowledge the financial support of the Max Planck Society, CY
Initiative of Excellence (grant ”Investissements d’Avenir” ANR-16-IDEX-0008),
the Programme National ”Physique et Chimie du Milieu Interstellaire” (PCMI) of
CNRS/INSU with INC/INP co-funded by CEA and CNES, a funding programme of the
Region Ile de France.

## References

  * Boogert et al. (2015a) Boogert, A. A., Gerakines, P. A., & Whittet, D. C. 2015a, Annual Review of Astronomy and Astrophysics, 53, 541
  * Boogert et al. (2015b) Boogert, A. A., Gerakines, P. A., & Whittet, D. C. 2015b, Annual Review of Astronomy and Astrophysics, 53, 541
  * Bovolenta et al. (2020) Bovolenta, G., Bovino, S., Vöhringer-Martinez, E., et al. 2020, Molecular Astrophysics, 21, 100095
  * Brünken et al. (2014) Brünken, S., Sipilä, O., Chambers, E. T., et al. 2014, Nature, 516, 219
  * Caselli et al. (2022) Caselli, P., Pineda, J. E., Sipilä, O., et al. 2022, The Astrophysical Journal, 929, 13
  * Caselli et al. (1999) Caselli, P., Walmsley, C., Tafalla, M., Dore, L., & Myers, P. 1999, The Astrophysical Journal, 523, L165
  * Cheung et al. (1968) Cheung, A. C., Rank, D. M., Townes, C. H., Thornton, D. D., & Welch, W. J. 1968, Phys. Rev. Lett., 21, 1701
  * Crabtree et al. (2011) Crabtree, K. N., Indriolo, N., Kreckel, H., Tom, B. A., & McCall, B. J. 2011, The Astrophysical Journal, 729, 15
  * Crapsi et al. (2007) Crapsi, A., Caselli, P., Walmsley, M. C., & Tafalla, M. 2007, Astronomy & Astrophysics, 470, 221
  * Crimier et al. (2010) Crimier, N., Ceccarelli, C., Maret, S., et al. 2010, Astronomy & Astrophysics, 519, A65
  * Danby et al. (1988) Danby, G., Flower, D. R., Valiron, P., Schilke, P., & Walmsley, C. M. 1988, MNRAS, 235, 229
  * Endres et al. (2016) Endres, C. P., Schlemmer, S., Schilke, P., Stutzki, J., & Müller, H. S. P. 2016, Journal of Molecular Spectroscopy, 327, 95
  * Fehér et al. (2022) Fehér, O., Tóth, L. V., Kraus, A., et al. 2022, The Astrophysical Journal Supplement Series, 258, 17
  * Ferrero et al. (2020) Ferrero, S., Zamirri, L., Ceccarelli, C., et al. 2020, The Astrophysical Journal, 904, 11
  * Garrod & Herbst (2006) Garrod, R. & Herbst, E. 2006, Astronomy & Astrophysics, 457, 927
  * Germain et al. (2022) Germain, A., Tinacci, L., Pantaleone, S., Ceccarelli, C., & Ugliengo, P. 2022, ACS Earth and Space Chemistry
  * Gorski et al. (2018) Gorski, M., Ott, J., Rand, R., et al. 2018, The Astrophysical Journal, 856, 134
  * Grassi et al. (2020) Grassi, T., Bovino, S., Caselli, P., et al. 2020, Astronomy & Astrophysics, 643, A155
  * Hama & Watanabe (2013) Hama, T. & Watanabe, N. 2013, Chemical reviews, 113, 8783
  * Harju et al. (2017) Harju, J., Sipilä, O., Brünken, S., et al. 2017, The Astrophysical Journal, 840, 63
  * He et al. (2016) He, J., Acharyya, K., & Vidali, G. 2016, The Astrophysical Journal, 825, 89
  * Irvine et al. (1987) Irvine, W., Goldsmith, P., Hjalmarson, Å., Hollenbach, D., & Thronson, H. 1987, Reidel Dordrecht, 561
  * Juvela (2020) Juvela, M. 2020, Astronomy & Astrophysics, 644, A151
  * Kakkenpara Suresh et al. (2024) Kakkenpara Suresh, S., Dulieu, F., Vitorino, J., & Caselli, P. 2024, Astronomy & Astrophysics, 682, A163
  * Keto & Caselli (2010) Keto, E. & Caselli, P. 2010, Monthly Notices of the Royal Astronomical Society, 402, 1625
  * Liszt et al. (2006) Liszt, H., Lucas, R., & Pety, J. 2006, Astronomy & Astrophysics, 448, 253
  * Loreau et al. (2023) Loreau, J., Faure, A., Lique, F., Demes, S., & Dagdigian, P. J. 2023, MNRAS, 526, 3213
  * Mundy et al. (1990) Mundy, L. G., Wootten, H., & Wilking, B. A. 1990, Astrophysical Journal, Part 1 (ISSN 0004-637X), vol. 352, March 20, 1990, p. 159-166., 352, 159
  * Pagani et al. (2007) Pagani, L., Bacmann, A., Cabrit, S., & Vastel, C. 2007, A&A, 467, 179
  * Penteado et al. (2017) Penteado, E., Walsh, C., & Cuppen, H. 2017, The Astrophysical Journal, 844, 71
  * Poch et al. (2020) Poch, O., Istiqomah, I., Quirico, E., et al. 2020, Science, 367, eaaw7462
  * Salinas et al. (2016) Salinas, V. N., Hogerheijde, M. R., Bergin, E. A., et al. 2016, Astronomy & Astrophysics, 591, A122
  * Sandqvist et al. (2017) Sandqvist, A., Hjalmarson, Å., Frisk, U., et al. 2017, Astronomy & Astrophysics, 599, A135
  * Schöier et al. (2005) Schöier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., & Black, J. H. 2005, A&A, 432, 369
  * Sipilä (2012) Sipilä, O. 2012, Astronomy & Astrophysics, 543, A38
  * Sipilä et al. (2015a) Sipilä, O., Caselli, P., & Harju, J. 2015a, Astronomy & Astrophysics, 578, A55
  * Sipilä et al. (2015b) Sipilä, O., Harju, J., Caselli, P., & Schlemmer, S. 2015b, Astronomy & Astrophysics, 581, A122
  * Sipilä et al. (2010) Sipilä, O., Hugo, E., Harju, J., et al. 2010, Astronomy & Astrophysics, 509, A98
  * Tafalla et al. (2002) Tafalla, M., Myers, P., Caselli, P., Walmsley, C., & Comito, C. 2002, The Astrophysical Journal, 569, 815
  * Taquet et al. (2016) Taquet, V., Wirström, E. S., & Charnley, S. B. 2016, The Astrophysical Journal, 821, 46
  * Tinacci et al. (2022) Tinacci, L., Germain, A., Pantaleone, S., et al. 2022, ACS Earth and Space Chemistry, 6, 1514
  * Wakelam et al. (2015) Wakelam, V., Loison, J.-C., Herbst, E., et al. 2015, The Astrophysical Journal Supplement Series, 217, 20
  * Wakelam et al. (2017) Wakelam, V., Loison, J.-C., Mereau, R., & Ruaud, M. 2017, Molecular Astrophysics, 6, 22
  * Winnewisser et al. (1979) Winnewisser, G., Churchwell, E., & Walmsley, C. 1979, Astronomy and Astrophysics, vol. 72, no. 1-2, Feb. 1979, p. 215-223., 72, 215

##  Appendix A Relative Ice Abundances of Key Volatiles

In our model, low ice abundances in the central regions of the source,
attributed to the high temperatures (approximately 200K), render these ice
abundance values less meaningful. As a result, in Fig. 7, we present simulated
abundance ratios beyond 2000 AU only, where ice abundances are higher due to
temperatures dropping to around 10-20 K. Our findings for CO and CH3OH align
well with the range observed in low-mass young stellar objects (LYSOs) as
described in Boogert et al. (2015b). Although our estimates for CH4 abundances
are slightly elevated, they remain reasonably close to the values reported in
the Boogert et al. (2015b). In our model, CO2 exhibits levels lower by two
orders of magnitude, while NH3 is overestimated by a factor of three. These
discrepancies may stem from a variety of factors, such as the elemental
abundances used in our model, the impact of background emissions on the
observed line intensities, or missing grain-surface chemistry (in the case of
CO2). A comprehensive discussion of the cause of these variations exceeds the
scope of our current study and is therefore omitted.

Figure 7: Relative ice abundances of key volatiles with respect to water ice.
The solid lines represent values obtained in this work and the shaded zones of
the same colour represent the range of observed values reported in Boogert et
al. (2015b)


# Transdimensional 2D Full-Waveform Inversion and Uncertainty Estimation

Reetam Biswas11footnotemark: 1 and Mrinal K. Sen11footnotemark: 1

###### Abstract

Full-Waveform Inversion (FWI) has now become a widely accepted tool to obtain
high-resolution velocity models from seismic data. It restores the high-
frequency features in a velocity model using both the time and amplitude
information from the seismic data. Typically, the velocity model in its
discrete form is represented on a rectangular grid, which is uniformly
distributed throughout the model space, and we solve for the elastic
properties at these grid points. FWI is mostly solved employing a local
optimization method, where one obtains a velocity update by minimizing the
misfit between the observed and the calculated seismograms. Note also that FWI
is a highly non-linear problem which is known to be prone to non-uniqueness.
The convergence to a globally optimum solution is not guaranteed; it depends
on the choice of the starting model. Thus, a Bayesian formulation of the
inverse problem with subsequent sampling of the posterior distribution is a
preferred choice, since it enables uncertainty quantification. However, with
the increase in the dimension of a model, sampling search space becomes
computationally expensive. We employ a recently developed trans-dimensional
sampling method called Reversible Jump Hamiltonian Monte Carlo (RJHMC), to the
2D full waveform inversion problem. We represent our velocity model using
Voronoi cells, determined from the distribution of certain nuclei points in
the model space. This method offers two advantages. First, it solves for a
variable dimensional velocity updates by using a trans-dimensional reversible
jump Markov Chain Monte Carlo (RJMCMC) step and thus tries to achieve an
optimum number of nuclei to represent the model and minimize the misfit. A
smaller number of parameters helps in an efficient sampling of the model
search space. Second, it applies the gradient-based Hamiltonian Monte Carlo
(HMC) step, which further improves the sampling by allowing the algorithm to
take a large step guided by the gradient. This two-step algorithm proves to be
a useful tool for model exploration and uncertainty quantification in FWI.

Trans-dimensional 2D FWI using RJHMC

## Introduction

Unlike travel-time tomography and Amplitude versus offset (AVO), FWI
(Tarantola,, 1986; Virieux and Operto,, 2009) uses both amplitude, and phase
information to generate subsurface velocity models using an iterative least-
square data-fitting procedure. It minimizes the data misfit by updating the
velocity model with a scaled gradient of the misfit function. Generally, the
adjoint state method (Lailly and Bednar,, 1983; Tarantola,, 1984; Pratt and
Worthington,, 1990) is used to calculate the gradient of the misfit function.
One of the most straightforward formulations of FWI is the constant-density
acoustic FWI, which solves for the P-wave velocity only.

Commonly, we use a pre-determined number of model unknowns in a geophysical
inversion problem, namely, the number of layers or cells to represent a
velocity model. There exist several statistical techniques to understand the
choice of model parametrization required for optimal data misfit like F-tests
(Aster et al.,, 2011). The above formulation requires the FWI problem to have
a fixed dimension of the model parameters. However, the subsurface properties
vary greatly with space, and thus the sensitivity of seismic observation also
varies spatially. In FWI, an over-parameterized model leads to poor
conditioning of the inverse problem, and is mostly tackled by introducing
regularization. The regularization term damps the solution and smooths the
model (Sen and Roy,, 2003; Virieux and Operto,, 2009). In geophysical
problems, fixing the number of parameters may lead to two cases. Starting an
inversion with too few parameters will underfit the data, while on the other
hand, starting with too many parameters will overfit the data (Sen and
Stoffa,, 1991). Both of these conditions lead to erroneous parameter
estimation with extreme uncertainty (Dosso et al.,, 2014). A trans-dimensional
formulation of the problem helps in overcoming these shortcomings.

The problem of determining the model dimension during inversion, i.e.,
treating the number of parameters as a parameter itself, has recently received
some attention in the geophysical community. Solving such an inverse problem
may seem unrealistic. Note that a Bayesian sampling framework has a built-in
parsimonious character; a high-dimensional model is naturally discouraged
(e.g., Malinverno,, 2002). Statisticians have actively developed several
trans-dimensional Markov chain Monte Carlo (MCMC) methods that allow
transition between states having different dimensions, i.e., simultaneously
traversing both parameter and model space (for a review, see Sisson,, 2005).
Out of these, the most commonly used is the reversible jump algorithm (RJMCMC)
of Green, (1995) (see also Green and Mira,, 2001; Sambridge et al.,, 2006) for
exploring variable dimensional models. Betancourt, (2017) showed that the
parameter sub-space becomes more and more complex as the parametrization moves
towards the higher dimensional space, i.e., increasing the number of
parameters, the tension between the density and volume increases exponentially
and the region from where significant contribution could come becomes
narrower. Using a trans-dimensional approach helps to eradicate some of these
drawbacks. It shows a parsimonious character and tries to represent the model
using the optimal number of required parameters, keeping the model simple
(enough such that it is neither under nor over-parameterized); thus, achieving
better exploration and uncertainty quantification of model space in comparison
to the same model with more parameters and using a fixed-dimensional inversion
scheme (Dettmer et al.,, 2016; Olugboji et al.,, 2017).

Malinverno and Leaney, (2000) first applied the reversible jump algorithm to
the problem of inversion of a zero-offset vertical seismic profile (for a
complete treatment see Malinverno and Leaney,, 2005). Subsequently,
Malinverno, (2002) applied it to an electrical-resistivity inversion. Later,
Agostinetti and Malinverno, (2010) extended the work to solve for velocity
models using receiver-function data. Bodin and Sambridge, (2009) applied it to
a seismic tomography problem using Voronoi cells. Ray and Key, (2012) apply
trans-dimensional inversion to marine controlled-source electromagnetic (CSEM)
data, and later Ray et al., (2016) solved a frequency-domain 1D elastic FWI
problem. Dosso et al., (2014) applied RJMCMC to Geo-acoustic profile
estimation; Dadi et al., (2016) for velocity log up-scaling using RJMCMC. Even
though RJMCMC has better model exploration than traditional fixed-dimensional
Bayesian inversion, it faces the problem of having low model acceptance ratio.
Al-Awadhi et al., (2004) proposed a hybrid technique to combine RJMCMC with
some other method to improve the model acceptance probability. Most recently,
Sen and Biswas, (2017) combined RJMCMC with the Hamiltonian Monte-Carlo
(RJHMC) to speed up convergence and applied it to 1D seismic inversion.
Hawkins et al., (2019) used HMC as one of the cases alongside the birth-death
scheme; each gets triggered by some probability one at a time.

Using HMC can boost the convergence as well as the acceptance ratio of RJMCMC;
Sen and Biswas, (2017) show a comparison of using a conventional RJMCMC with
the RJHMC for a 1D seismic inversion problem. Although RJHMC converges in
fewer iterations, it comes at the cost of additional forward problems to
calculate the gradient since it uses gradient information to generate a model
update and minimize the problem of sampling correlated models which most
Bayesian techniques face. Recently, HMC has gained popularity in solving
geophysical problems (for e.g., see Sen and Biswas,, 2017; Fichtner et al.,,
2018; Hawkins et al.,, 2019). However, to achieve a good sampling in HMC, it
needs to be tuned, which is somewhat problem-dependent. Betancourt, (2017),
and Fichtner et al., (2018) discuss different possible ways to tune the Mass
matrix in HMC to have a minimum correlation between the samples and explore
the model space.

One other limitation of FWI is the requirement of a good starting model.
Although this problem has recently been addressed using global optimization
methods (e.g., Datta and Sen,, 2016), no attempt has yet been made to
characterize uncertainty using a fully nonlinear sampling method. Zhu et al.,
(2016) estimate uncertainty assuming the posterior probability distribution
(PPD) to be Gaussian and starting with a model obtained from a local
optimization as the prior. Kotsi et al., (2018) used Bayesian framework to
address uncertainty in estimating velocity changes in 4D FWI. In this paper,
we primarily address these two issues using a trans-dimensional approach. We
treat the model parametrization and velocity values as unknown in the
inversion and estimate the uncertainty. Here, we apply RJHMC to a
computationally intensive 2D FWI problem and demonstrate its feasibility and
usefulness using a noisy synthetic dataset.

## Theory

To solve the FWI problem, we first describe the parametrization used in
describing a given velocity model and then the forward modeling algorithm.
Finally, we describe the inverse formulation for model exploration based on
RJHMC.

###  0.1 Model Parametrization

To solve a trans-dimensional 2D FWI problem using RJHMC, we choose to
represent our model space using an adaptive ensemble of nuclei. The 2D P-wave
velocity model is partitioned using a variable number of non-overlapping
Voronoi tessellation (Okabe,, 1992; Bodin et al.,, 2009). The Voronoi nuclei
constitute parametrization of the model space, and each Voronoi cell encloses
all the points in space which are closer to the nucleus of that Voronoi cell
than any other Voronoi nuclei. Here, we treat the number and position of these
Voronoi nuclei as unknown, thus making the Voronoi cell mesh evolving along
with the nucleus. The Voronoi nuclei determine the elastic property
distribution in the model. The nuclei spatial distribution, as well as their
number, can adaptively change as a function of iteration during sampling to an
efficient number of parameters so as to match the synthetic data with the
observed data and reduce the misfit. These may lead to concentrating nuclei
where high-resolution is required in the model, and sparser mesh elsewhere
(Agostinetti et al.,, 2015; Bodin et al.,, 2012; Sambridge et al.,, 1995;
Agostinetti et al.,, 2015).

To demonstrate an example, LABEL:fig:voronoiexample shows a seismic velocity
model discretized using a set of Voronoi polygons. Let us consider a set C of
$k$ nuclei in the 2D plane ($\textbf{C}={c_{1},...,c_{k}}$ where
$c_{i}\in\mathbb{R}^{2}$ ). The Voronoi tessellation defines $k$ mutually
exclusive regions $R_{1}$, $R_{2}$,…,$R_{k}$ such that the points in the area
which fall within the region $R_{i}$ are closer to $c_{i}$ than other
$c_{j}(j\neq i)$. Note that the cell nuclei may not necessarily represent the
geometric center. The boundaries between the cells represent the perpendicular
bisectors of each pair of neighboring nuclei. Thus in a 2D space, each nucleus
$c_{i}$ is defined by its position ($x$ and $z$ coordinates) and a seismic
wave velocity $v_{i}$. Therefore, the model parameters having $k$ number of
nuclei can be represented by

$\mathbf{m}=(k,\mathbf{C},\mathbf{V}),$ (1)

where $\mathbf{V}$ is the vector of the velocity values assigned to each
nuclei. Thus the effective number of unknowns (dimension of the model) is $3k$
(for each nucleus, it’s X and Z position, and its velocity). During the
inversion, the number of nuclei, their position, and elastic property are
variable and thus the generated Voronoi cell can take different size, shape,
and velocity. This partition modeling approach directly leads to a dynamic
parametrization and makes the data determine the cell geometry and the grid
resolution to achieve a maximum resolution power of data.

To generate the full-scale 2D fine mesh grid velocity model from the nuclei,
we first create the Voronoi map between the nuclei; the grid points which fall
within the Voronoi area are assigned a velocity from that respective nucleus.
In LABEL:fig:voronoiexample, the white cross-hairs (does not represent all the
grid points in the figure for simplicity) show the fine mesh grid laid on top
of the Voronoi map. Using Voronoi cells as a parametrization may have a
discontinuous border for neighboring cells having different velocities,
therefore a smoothing is performed to minimize artifacts due to the edge
reflection.

voronoiexamplewidth=Figure showing the velocity model created from nuclei
cells represented in black dot. The white cross does not represent all the
grid points in the figure for simplicity are the grid points in the velocity
model. The color in the Voronoi cells are the velocity of that respective
nuclei in that area.

###  0.2 Forward Modeling

Once a gridded model is obtained, we perform forward modeling using the
constant-density acoustic wave equation, where the only elastic parameter is
the P-wave velocity. Such type of wave equation is given by

$\frac{1}{c^{2}}\frac{\partial^{2}P}{\partial t^{2}}=\nabla^{2}P+S(x,t),$ (2)

where $P$ is the pressure wavefield, $c$ is the P-wave velocity, $\nabla^{2}$
is the Laplacian given by $\frac{\partial^{2}}{\partial
x^{2}}+\frac{\partial^{2}}{\partial z^{2}}$, and $S(x,t)$ is the source term.

At each iteration, after the velocity model is perturbed, the wave equation is
solved to compute new seismogram ($\mathbf{d}_{cal}$) and the misfit from the
observed seismogram ($\mathbf{d}_{obs}$) is calculated using

$E=\|\mathbf{d}_{obs}-\mathbf{d}_{cal}\|_{2}.$ (3)

The gradient is calculated by the adjoint-state method (Lailly and Bednar,,
1983; Tarantola,, 1984; Pratt and Worthington,, 1990), where the data residual
is backward propagated in reverse time from the receivers to calculate the
adjoint wavefield by

$\frac{1}{c^{2}}\frac{\partial^{2}R}{\partial t^{2}}=\nabla^{2}R+\Delta
S(x,t),$ (4)

where $R$ is the adjoint wavefield, and $\Delta S(x,t)$ is the data residual
between the observed data ($\mathbf{d}_{obs}$) and the calculated data
($\mathbf{d}_{cal}$). Then the gradient of the misfit function is calculated
by summing the zero-lag cross-correlation of the two wavefields over multiple
shots

$\frac{\partial E}{\partial
m}=\frac{1}{c(x,z)^{3}}\sum_{shots}\ddot{P}(x,z,t)R(x,z,t),$ (5)

where P and R are the source and adjoint wavefields in space and time.

###  0.3 Trans-dimensional Bayesian Formulation

We employ a Bayesian approach to full waveform inversion with the goal of
estimating multiple solutions for uncertainty quantification. Choice of prior
probability distribution is crucial to sample the regions near the global
minimum efficiently (Tarantola,, 1984; Buland and Omre,, 2003; Bosch et al.,,
2010; Karimi et al.,, 2010). In a Bayesian framework (Mosegaard and
Tarantola,, 1995), information about the model parameters is expressed by a
probability distribution that describes the associated uncertainties, known as
the posterior probability distribution (PPD). It represents the intersection
of the set of possible solutions, which satisfy our prior knowledge, and our
observations. The PPD can be represented by the Bayes rule as follows:

$\hat{\pi}(\mathbf{m})=\frac{p(\mathbf{m})p(\mathbf{d|m})}{p(\mathbf{d})},$
(6)

where $p()$ is a probability distribution, $A|B$ indicates conditional
dependence between $A$ and $B$, and d denotes the calculated data. Thus,
$\hat{\pi}(\textbf{m})$ is the PPD, $p(\textbf{m})$ represents the prior
information available on the model $\textbf{m}\in\mathbb{R}^{k}$ ($k$ is
generally assumed to be known), and $p(\textbf{d}|\textbf{m})$ is the
likelihood (also represented as $l(\textbf{m})$), which is the probability as
to how well the synthetic data generated by the given model m matches the
observed data d. The denominator $p(\textbf{d})$ is called the “evidence” or
the marginal likelihood (Tarantola and Valette,, 1982; Hirata and Matsu’ura,,
1987; Cary and Chapman,, 1988; Sen and Stoffa,, 1996; Bernardo and Smith,,
2001). It is defined as the numerator on the right-hand side integrated over
the entire model space and thus regarded as a constant. Therefore, we can
write Equation 6 in the form of a proportionality relationship:

$\hat{\pi}(\mathbf{m})\propto p(\mathbf{m})p(\mathbf{d|m}).$ (7)

We use Markov Chain Monte Carlo (MCMC) type of method to sample models from
the PPD using Equation 7. The new model ($\textbf{m}^{\prime}$) is sampled
from an assumed proposal distribution ($q(\textbf{m}^{\prime}|\textbf{m})$)
and is accepted using Metropolis-Hastings acceptance criterion as

$\alpha(\textbf{m}^{\prime}|\textbf{m})=\min\Bigg{[}1,\frac{q(\textbf{m}|\textbf{m}^{\prime})}{q(\textbf{m}^{\prime}|\textbf{m})}\frac{p(\textbf{m}^{\prime})}{p(\textbf{m})}\frac{\ell(\textbf{m}^{\prime})}{\ell(\textbf{m})}\Bigg{]}.$
(8)

The above MCMC algorithm draws samples from a fixed number of model
parameters. The most commonly used trans-dimensional MCMC method is the
reversible jump MCMC algorithm. In RJMCMC, the random walk can jump between
the states having different number of unknowns. Solving a trans-dimensional
problem needs some modification to the Bayes’ rule (Malinverno and Leaney,,
2000) to accommodate the variability in the model parametrization as

$\pi(k,\textbf{m})=\frac{p(k)p(\textbf{m}|k)p(\textbf{d}|\textbf{m},k)}{\sum_{k\in\kappa}\int
p(k)p(\textbf{m}|k)p(\textbf{d}|\textbf{m},k)dm},$ (9)

where $\textbf{m}\in\mathbb{R}^{c\times k}$ represents a model in the $M_{k}$
dimensional parameter search space associated with $k$ layers. In our FWI
problem, $c=3$. Here, $p(k)p(\textbf{m}|k)$ represents the combined prior
probability of the state $(k,\textbf{m})$ (assuming independent prior), which
is the product of prior for the model dimension $p(k)$ and the prior for the
model parameter values $p(\textbf{m}|k)$. The function
$p(\textbf{d}|\textbf{m},k)$ represents the likelihood of the state
$(k,\textbf{m})$, also written as $\ell(\textbf{m})$. The integration in the
denominator represents the normalization factor (also known as evidence) in
case of a trans-dimensional problem. It integrates out the product of prior
and the likelihood for all possible values of $k\in\kappa$, where $\kappa$
represents the different number of parametrization cells possible. Similar to
Equation 7, we can represent it as

$\pi(k,\mathbf{m})\propto p(k)p(\mathbf{m}|k)p(\mathbf{d|m},k).$ (10)

Green, (1995) defines the acceptance ratio for accepting a new sampled model
$(\textbf{m}^{\prime},k^{\prime})$ given the current model $(\textbf{m},k)$ as

$\alpha(k^{\prime},\textbf{m}^{\prime}|k,\textbf{m})=\Bigg{[}1,\frac{\pi(\textbf{m}^{\prime},k^{\prime})q(\textbf{m},k)}{\pi(\textbf{m},k)q(\textbf{m}^{\prime},k^{\prime})}|\textbf{J}|\Bigg{]}.$
(11)

This equation shows the transition from model
$\textbf{m}\in\mathbb{R}^{c\times k}$ (dimension of $c\times k$) to a model
$\textbf{m}^{\prime}\in\mathbb{R}^{c\times k^{\prime}}$ (dimension of $c\times
k^{\prime}$). $q(\textbf{m})$ is the proposal distribution used in sampling a
new model. Here, $|\textbf{J}|$ represents the determinant of the Jacobian
matrix, which takes into account a change in the dimension. The Jacobian term
acts as a normalization factor for the change in volume between two different
dimensional spaces.

##  1 Reversible Jump Hamiltonian Monte Carlo (RJHMC)

Here, we follow the Reversible Jump Hamiltonian Monte Carlo algorithm (Sen and
Biswas,, 2017), which is a stochastic method based on Bayes’ rule. It consists
of two steps (Al-Awadhi et al.,, 2004). The first step is a trans-dimensional
step similar to RJMCMC. In this step, the model $\mathbf{m}$ with $k$ nuclei
may jump towards an increase in dimension to the proposed state
$\mathbf{m}^{\prime}$ having $k^{\prime}=k+n$ nuclei. The prior for the model
dimension is given by $p(\mathbf{m}^{\prime},k^{\prime})$. The second step is
a fixed-dimensional Hamiltonian Monte Carlo (HMC) (Neal,, 2011), where the
model $\mathbf{m}^{\prime}$ moves to a new state $\mathbf{m}^{*}$ having the
same number of nuclei points $k^{\prime}$. Note that in the reverse step of
this move, we first perform the fixed dimensional HMC to change the state of
the model from $\mathbf{m}^{*}$ to $\mathbf{m}^{\prime}$ having $k^{\prime}$
number of nuclei and then perform the trans-dimensional step to move model
$\mathbf{m}^{\prime}$ having $k^{\prime}$ nuclei to model state $\mathbf{m}$
with $k$ nuclei. The RJHMC acceptance probability gives the transition from
model $\mathbf{m}$ $(k)$ to model $\mathbf{m^{\prime}}$ $(k^{\prime})$ to
final model $\mathbf{m}^{*}$ $(k^{\prime})$.

####  1.0.1 First: Trans-dimensional step (RJMCMC)

For a first trans-dimensional step, the random moves can jump between the
states having different model dimensions (number of nuclei here). Here, we
have the transition from model m with $k$ nuclei to model
$\textbf{m}^{\prime}$ having $k^{\prime}$ nuclei. We have adopted the birth-
death type of scheme (Agostinetti et al.,, 2015) to carry out the trans-
dimensional step. This scheme consists primarily of three different moves. In
each of these moves, a candidate model is generated by randomly perturbing the
current model following the given proposal distribution. The three moves are
as follows:

  1. 1.

_No Perturbation_ : In this move, we do not perturb the model space and keep
it the same. This step has a half probability.

  2. 2.

_Birth_ : We add $n$ number of nuclei at a randomly chosen available locations
using a prior probability distribution. The new nucleus is added by keeping in
mind that the distance between any two nearby nuclei is greater than the
minimum distance allowed. This step has a quarter probability.

  3. 3.

_Death_ : We remove $n$ number of nuclei randomly from the current nuclei.
This step also has a quarter probability.

The acceptance ratio for a Bayesian inversion dealing with change in dimension
has a Jacobian matrix as described in Equation 11. Now, for a birth-death type
of approach described above, the determinant of this Jacobian matrix can be
shown to be unity (please refer to Appendix B in Sen and Biswas, (2017) for
full detail), i.e., $|\textbf{J}|=1$. Now we will describe the various
distributions used in the trans-dimensional part of the algorithm.

##### Prior Distribution ($p(\textbf{m})$):

Prior probability in a Bayesian algorithm represents our prior knowledge about
the model. The distribution can be in the form of some complex formulation to
a simple bound applied to the model. Here, we assume that the two parameters
(number of nuclei, and velocity) in the model m are not correlated a priori so
that the complete prior distribution can be represented as a product of the
prior distribution of the two:

$p(\textbf{m})=p(k)p(\textbf{m}|k).$ (12)

The prior on the number of parameters (partitions) $p(k)$ is chosen to be a
bounded uniform distribution between the minimum number of nuclei $k_{min}$
and the maximum number of nuclei possible $k_{max}$, represented as:

$\displaystyle p(k)=\begin{cases}\frac{1}{({k_{max}-k_{min}})}\quad\text{if
}k_{min}\leq k\leq k_{max}\\\
0\quad\quad\quad\quad\quad\quad\quad\quad\text{Otherwise}.\\\ \end{cases}$
(13)

The prior distribution on the parameter values ($p(\mathbf{m}|k)$) is defined
as a Gaussian distribution as

$\displaystyle
p(\mathbf{m}|k)=\frac{1}{\big{[}(2\pi)^{k}\det(\mathbf{\bar{C}}_{p})\big{]}^{1/2}}\exp\Bigg{[}-\frac{1}{2}(\mathbf{m}-\mathbf{\bar{m}})^{T}(\mathbf{\bar{C}}_{p})^{-1}$
$\displaystyle(\mathbf{m}-\mathbf{\bar{m}})\Bigg{]},$ (14)

where $\mathbf{\bar{C}}_{p}$ is the prior covariance matrix and
$\mathbf{\bar{m}}$ represents the prior mean model.

##### Likelihood Function ($p(\textbf{d}|\textbf{m},k)$):

This plays a major role in the acceptance probability of Equation 9. This term
quantifies the measure of error between the predicted data and the observed
data. In general, we assume the error between the observed and the predicted
data to be normally distributed and is quantified by a covariance matrix
$\textbf{C}_{err}$. It takes into account various factors like noise in the
data, approximation error in modeling, etc. Estimating $\textbf{C}_{err}$
requires extensive work and is not considered in our study. Here, we assume a
simple case in which the error between different data points are uncorrelated
and thus can represent the covariance matrix as a diagonal matrix of the same
standard deviation $\sigma_{err}$, i.e.,
$\textbf{C}_{err}=\sigma^{2}_{err}\textbf{I}$, where I is an identity matrix.
Thus, the likelihood function $p(\mathbf{d}|\textbf{m},k)$ in Equation 9 can
be represented as

$p(\textbf{d}|\textbf{m},k)=\frac{1}{\big{[}(2\pi\big{)}\sigma_{err}]^{N/2}}\exp\Big{[}-\frac{1}{2}\Big{(}\frac{E}{\sigma_{err}}\Big{)}^{2}\Big{]}.$
(15)

####  1.0.2 Second: Fixed-dimensional step (HMC)

In this section we describe a fixed Dimensional MCMC step, where the model in
state $\textbf{m}^{\prime}$ with $k^{\prime}$ nuclei moves to a new state
$\textbf{m}^{*}$ with same $k^{\prime}$ nuclei using Hamiltonian Monte Carlo.
Duane et al., (1987) proposed this method based on an approximate simulation
of Hamiltonian dynamical system using the potential energy given by log-target
density. It helps in avoiding issues with high curvature and concentration of
measures. It produces a distant proposal for the Metropolis algorithm and
avoids slow exploration by using gradient information while changing states
(Betancourt,, 2017) . In HMC, we calculate the total potential energy $(U)$
and kinetic energy $(K)$ of the state, i.e., the Hamiltonian Energy. HMC
augments the state space of the target distribution by adding auxiliary sets
of variables termed as momentum p. Each nucleus is augmented with an initial
momentum value derived from a normal distribution $N(0,1)$. The total
Hamiltonian energy of the model m with momentum p can be given by

$H(\textbf{m}^{*},\textbf{p})=U(\textbf{m}^{*})+K(\textbf{p}).$ (16)

The potential energy $U$ is defined as the negative logarithm of the posterior
distribution given as

$U(\textbf{m}^{*})=-\log\big{[}p(\textbf{m}^{*})p(\textbf{d}|\textbf{m}^{*})\big{]},$
(17)

where the prior probability $p(\textbf{m}^{*})$ is taken to be a uniform
bounded distribution and the likelihood $p(\textbf{d}|\textbf{m})$ is same as
the Equation 15. The kinetic energy term $K$ is calculated by the equation
given as

$K(\textbf{p})=\frac{\textbf{p}^{T}\textbf{M}^{-1}\textbf{p}}{2},$ (18)

where M is the mass matrix. Fichtner et al., (2018) have emphasized the
critical role of mass matrix in the tuning of the HMC sampling for a fixed-
dimensional problem. It can be used to put different weights to different
parameters and control their perturbations and achieve sampling with most
uncorrelated samples. In our trans-dimensional case, since the number of
parameters is continuously changing with iterations, we chose an identity
matrix. To perform the Hamiltonian step, the model
$(\textbf{m}^{\prime},k^{\prime})$ at the beginning of the second step is
assigned a momentum $(p)$ value for each nucleus. It has a Hamiltonian energy
$H(\textbf{m}^{\prime},\textbf{p})$. This momentum value is drawn randomly
from a normal distribution $N(0,1)$. Then, we simulate Hamiltonian dynamics
$L$ times. We have used the Leapfrog technique to simulate the Hamiltonian
dynamics (Van Gunsteren and Berendsen,, 1988; Neal,, 2011). In each
Hamiltonian step, we update the momentum $p$ with the help of gradient
information for a $\delta/2$ time step. The gradient is calculated using the
FWI adjoint-state method. We use this half updated momentum to update the
model for the complete $\delta$ time step, and again using the gradient
calculated from this new model, we update the momentum for a $\delta/2$ time
step, such that after this one loop both momentum and model have traversed a
complete $\delta$ time step. At the end of the Hamiltonian simulation, we have
a new state of model $\textbf{m}^{*}$ and momentum $\textbf{p}^{*}$ with
Hamiltonian energy $H(\textbf{m}^{*},\textbf{p}^{*})$. We calculate the
difference in the total Hamiltonian energy of the two states. Note that the
parameters $L$ and $\delta$ have to be carefully tuned for specific
applications. The updates for the position state and the momentum vector in
each iteration of HMC thus can be given by

$\displaystyle p_{i}^{new}$ $\displaystyle=$ $\displaystyle
p_{i}^{old}-\frac{\delta}{2}\times\frac{\partial U}{\partial m_{i}},$
$\displaystyle\textbf{m}_{new}$ $\displaystyle=$
$\displaystyle\textbf{m}_{old}-\delta\times\textbf{p}^{new},$ $\displaystyle
p_{i}^{new}$ $\displaystyle=$ $\displaystyle
p_{i}^{old}-\frac{\delta}{2}\times\frac{\partial U}{\partial m_{i}^{new}}.$
(19)

####  1.0.3 Acceptance Ratio for RJHMC

After completion of the two steps, i.e., a transition from a model
$(\textbf{m},k)\rightarrow(\textbf{m}^{\prime},k^{\prime})\rightarrow(\textbf{m}^{*},k^{\prime})$,
we accept this new model with an RJHMC acceptance probability given in Sen and
Biswas, (2017). Note that in the case of a birth step, we first have a trans-
dimensional step followed by an HMC step. However, in the reverse move, we
first go through the HMC step followed by a trans-dimensional death step. This
is to maintain reversibility of the chain. Now, for the above transition, we
can write the detail balance equation as

$\displaystyle\int_{(\textbf{m},\textbf{m}^{\prime},\textbf{m}^{*})\in A\times
B\times C}$
$\displaystyle\pi(\textbf{m})q(\textbf{m}^{\prime},{k^{\prime}})P(\textbf{m}^{*}|\textbf{m}^{\prime})\alpha(\textbf{m}^{*},k^{\prime}|\textbf{m},k)d\textbf{m}d\textbf{m}^{\prime}d\textbf{m}^{*}$
$\displaystyle=\int_{(\textbf{m},\textbf{m}^{\prime},\textbf{m}^{*})\in
A\times B\times C}$
$\displaystyle\pi(\textbf{m}^{*})P(\textbf{m}^{\prime}|\textbf{m}^{*})q(\textbf{m},{k})\alpha(\textbf{m},k|\textbf{m}^{*},k^{\prime})d\textbf{m}d\textbf{m}^{\prime}d\textbf{m}^{*},$
(20)

for all Borel sets $A\in\mathbb{R}^{c\times k}$, $B\in\mathbb{R}^{c\times
k^{\prime}}$, and $C\in\mathbb{R}^{c\times k^{\prime}}$. Sen and Biswas,
(2017) give the acceptance probability for the above inverse formulation. The
acceptance probability is

$\displaystyle\alpha(\textbf{m}^{*},k^{\prime}|\textbf{m},k)=$
$\displaystyle\min\Bigg{[}1,\frac{p(\textbf{m}^{*}|k^{\prime})p(\textbf{d}|\textbf{m}^{*},k^{\prime})}{p(\textbf{m}|k)p(\textbf{d}|\textbf{m},k)}\times$
(21)
$\displaystyle\frac{q(\textbf{m}|k)}{q(\textbf{m}^{\prime}|k^{\prime})}\times\frac{\exp(-H(\textbf{m}^{\prime}))}{\exp(-H({\textbf{m}^{*}}))}\Bigg{]}.$

where $q(\textbf{m}^{\prime}|k^{\prime})$ is the proposal distribution for the
model parameter values that is carefully chosen to be an approximation of a
posterior distribution, to achieve a high value of acceptance ratio.
LABEL:fig:flowchart shows the flowchart of the algorithm.

flowchartwidth=0.8Flowchart of RJHMC.

##  2 Result

We demonstrate our algorithm using a synthetic dataset generated using the
Marmousi model for P-wave velocity. The model has dimensions of $251$ (grid
points in depth) $\times 767$ (grid points in the horizontal direction), with
a grid spacing of 12 m. A 15 Hz Ricker wavelet was used to generate the
synthetic data set. We used 21 shots at an interval of $36$ grid spaces with
767 receivers per shot. Gaussian noise was added to the traces, which were
then used in the inversion. LABEL:fig:seismogram shows a shot gather for one
of the shots. We started with 1000 nuclei randomly distributed in our model
space, and we performed two steps of the leap-frog Hamiltonian Dynamics for
each iteration of RJHMC. In the trans-dimensional step of RJHMC, the algorithm
can add or delete up to a maximum of 30 nuclei at a time.

seismogramwidth=Figure shows one of the generated shot gathers. Note that the
direct arrivals have been muted.

Before applying the RJHMC algorithm to the Marmousi model, we performed a test
on the prior distribution. A dry run was made with no input data while setting
the likelihood function to be unity to sample the model space. With the
absence of a likelihood function, the algorithm should produce a posterior
identical to the prior distribution. For this test, we set $k_{max}$ to be
20000. LABEL:fig:kPrior represents almost a uniform distribution, i.e., all
possible values of nuclei count within the bound of $k$ are possible, similar
to the defined uniform distribution in Equation 13. Similarly,
LABEL:fig:modelPrior represents prior mean values of the P-wave velocity as
defined in the prior for the model parameter as in Equation 13. These have a
close resemblance with the prior distribution defined and thus we can conclude
that the algorithm is sampling correctly.

2kPrior,modelPriorwidth=(a) Figure showing the distribution of number of
nuclei, and (b) Figure showing the mean velocity model with likelihood
function set to unity.

To throw light on the convergence of the algorithm, we plot the normalized
error curve and the running mean plot of the P-wave velocity at five randomly
chosen locations. The error is generally high at the initial iterations,
since, the starting model can be quite far away from most likely models.
However, with the successive iteration of MCMC, current state reaches closer
to the target distribution (Betancourt,, 2017) and samples from there; at this
point, the error stabilizes. This initial part of exploration is called the
burn-in phase and the models after this phase samples the posterior
distribution and thus is used for calculating the uncertainty in model
parameters. We ran a total of 15000 iterations of RJHMC, considering the
computation complexity of the problem. LABEL:fig:error15k shows the normalized
error progression along with iteration. By examining the error curve, we note
that there is a stability in error at around 5000 iterations. We also
generated a running mean of $V_{P}$ at five randomly chosen locations in
LABEL:fig:runningMean. All five plots show excellent convergence. Based on
these two plots, we infer the completion of burn-in phase to be around 5000
iterations.

2error15k,runningMeanwidth=0.5(a) Plot of error curve along with iteration (b)
Plot of the running mean $V_{P}$ parameters for 5 randomly chosen locations.

The MCMC method performs poorly in a high-dimensional problem. With the
increase in dimension, the sample space narrows down significantly and thus
becomes more difficult to sample (Betancourt,, 2017), and the consecutive
samples become more and more correlated. This makes the posterior biased and
limited to the vicinity of the starting model. To test the capability of
RJHMC, we calculate the autocorrelation of the successive models sampled
during the run. LABEL:fig:corr represents the correlation matrix for 2500
successive models starting from the $5000^{th}$ iteration (after burn-in
phase). The matrix shows a high value near the diagonal and minimal values
close to zero at the far end. This shows that sample independence in RJHMC has
been achieved with successive iterations. LABEL:fig:autocorrelation shows the
auto-correlation coefficient of the starting model with the successive models
up to iteration 2500. This figure shows how rapidly the correlation
coefficient of the model falls to a small value.

2corr,autocorrelationwidth=0.7(a) Plot of the correlation matrix using the
models from iteration 5000 to 7500 (2501 models) (b) Plot of the correlation
coefficient of starting model with the successive models at various iterations
up to iteration 2500.

LABEL:fig:modelswithIter shows the progress of the model during sampling. The
first model is the true velocity model. We used a smooth starting model drawn
from the prior and sampled with 1000 nuclei as shown in
LABEL:fig:modelswithIterb. The models in the other panels show the velocities
at various iterations of RJHMC along with the nuclei count at that step. We
started with 1000 nuclei as a staring model (LABEL:fig:modelswithIterb). Note
that, in general, we observe an increase in the number of nuclei until
convergence. LABEL:fig:modelswithIterd shows the model from the last
iteration.

modelswithIterwidth=Figure showing velocity models at different iterations of
RJHMC: (a) the true velocity model, (b)-(d) models at progressive iterations.
Plots show the iteration number and the number of nuclei count in the model
represented as NC.

LABEL:fig:locationPPD2_1 shows the velocity models at $X=(2,5,8)$ km. It shows
the plot of the marginal PPD (grey cloud), true log (green), mean parameter
estimate (red), and 95% confidence interval (dashed line) for the wells marked
with green location. The well has the target reservoir of the Marmousi model
at a depth of 2.3 km that is fairly well resolved. Comparison of the Posterior
distribution at the three well locations in LABEL:fig:locationPPD2_1 shows
that the true model is almost everywhere within one posterior standard
deviation of the posterior mean model. However, in the deeper part of the
model, the posterior uncertainty is closer to the prior. LABEL:fig:kPPD15k1
shows the marginal posterior density for the number of nuclei counts. There is
a distinct peak at around 15000-16000 number of parameters (nuclei).

locationPPD2_1width=Plot of marginal PPD (grey cloud), true log (green), mean
parameter estimate (red), and 95% confidence interval (dashed line) for well
at location X=(2,5,8)Km marked with Green Location

kPPD15k1width=0.5Plot of marginal PPD for number of nuclei count

We observe that with depth, the uncertainty in the estimation of the P-wave
velocity increases. LABEL:fig:meanModel15k, LABEL:fig:mapModel,
LABEL:fig:medianModel and LABEL:fig:perstd15k1 show the mean model, maximum-a-
posterioiri model, median model and the standard deviation (as a percentage of
the mean model) at each grid. The maximum percentage standard deviation is
nearly 10%, and the standard deviation seems to be increasing with depth as
the data coverage becomes poorer with depth. A Bayesian inversion is an
essential tool for model exploration and uncertainty quantification. The
sampled posterior distribution of the model parameters helps to investigate
the quality and resolution power of the inversion. Comparing the prior
distribution with the posterior distribution gives an idea about the recovered
or improvement of resolution for the seismic model. In our experiment, we set
the standard deviation in the Prior distribution large. LABEL:fig:meanModel15k
and LABEL:fig:perstd15k1 display plots of posterior mean and posterior
standard deviation specifically. Note that, we show posterior standard
deviations relative to the corresponding mean values. Zones of low standard
deviations correspond to regions of well-resolved velocities while velocities
within zones of high standard deviation are poorly resolved. These are
generally dictated by experimental geometry and ray-path coverage in the
subsurface. For example, along the edges, the model resolution is poor, and
the standard deviation is also significant in the deeper parts of the models.
Note also that the standard deviation at the interface has a high relative
standard deviation in the region. It is due to the velocity-depth trade-off.
It represents that position of the nuclei, and the parameter value is not
entirely an independent quantity. To further investigate this observation, we
plot model parameter correlation matrix for a small vertical and horizontal
line containing around 45-50 grid points in LABEL:fig:verticalModelCovar and
LABEL:fig:horizontalModelCovar respectively. It is clearly evident that in the
vertical correlation matrix, as given in LABEL:fig:verticalModelCovar, the
model parameters are spread throughout and there is a very high correlation of
a parameter location with its neighbors. The same positive correlation is not
that prominent in the horizontal direction correlation matrix given in
LABEL:fig:horizontalModelCovar.

4meanModel15k,mapModel,medianModel,perstd15k1width=Plot of the estimated (a)
mean velocity model, (b) Maximum-a-posterior model, and (c) median model. (d)
Plot of the percentage standard deviation of of the estimated velocity model.

2verticalModelCovar,horizontalModelCovarwidth=0.45(a) Plot of the neighbouring
model parameter correlation matrix in (a) vertical direction (b) horizontal
direction.

Finally, LABEL:fig:voronoiComparison shows the comparison of the Voronoi cell
distribution and density at the beginning of the inversion and the end of the
inversion. We started with 1000 nuclei and later at the end of the RJHMC we
have around 15500 nuclei.

voronoiComparisonwidth=0.8Plot showing the comparison of the Voronoi cell
distribution and density at the beginning of the inversion and at the end of
the inversion. We started with 1000 nuclei count and later at the end of the
RJHMC we have around 15500 nuclei count.

##  3 Discussions

In this paper, we have demonstrated the applicability of a trans-dimensional
MCMC method for a realistic 2D full waveform inversion problem. The main
drawback of this method is the high computation and memory cost due to the
cross-correlation of the forward and the adjoint state wavefields at the same
time step. This requires saving the forward wavefields in the memory. Again,
due to the limited GPU memory, data transfers are much slower compared to the
computation, and the GPU waits for the data-copy to finish before continuing
the computation. Therefore, we have followed a two-level checkpointing
technique as given in Datta et al., (2018), where the high overhead of memory
usage is traded with the data movement and computations. For computation, we
used four K40 Nvidia GPUs and NVMe (Non-Volatile Memory Express memory) at
Institute for Geophysics, University of Texas at Austin. This enables us to
have the data movement in the background, while the GPUs are continuously in
use.

As mentioned earlier, the Bayesian inversion has an intrinsic parsimony
character (Malinverno,, 2002). Given the set of earth models having a
different number of model parameters and fitting the data, models with fewer
parameters are favored and have a higher posterior distribution. This feature
can also be seen in our RJHMC result. In the first case, when we set the
likelihood to unity, the sampled model is not influenced by the data. Thus, we
can see in LABEL:fig:kPrior all sampled models have equal probability.
However, when the data influences the likelihood function, parametrization is
dictated by a lower dimensional model due to the parsimonious nature of
Bayesian inversion. This is evident in LABEL:fig:kPPD15k1, where there is a
peak at around 16000 nuclei, though a much higher model dimension is possible.
This is around $8.5\%$ of the total nuclei possible (maintaining minimum
separable distance).

Theoretically, the choice of parametrization used while solving an inverse
problem does not prevent the model exploration to the target distribution if
more than enough samples are drawn. However, for practical purposes, this
sometimes cannot be achieved, and therefore different choice of
parametrization produces different posterior distribution. In our application,
Voronoi cells have been used for representing the velocity structure. However,
there can be several other choices of parametrization which can be exploited
depending on the problem. Hawkins et al., (2019) show a comparison of the
result produced by three different types of parametrization: Voronoi, Delaunay
triangulation with linear interpolants and with Clough-Tocher interpolant.
They showed that for models with a sharp transition or discontinuities,
Voronoi performs better than the other two. However, in case of a smoothly
varying model, Delaunay triangulation with some interpolant would do better.
Several sharp changes in velocity are expected in models like the Marmousi.
Therefore, we choose Voronoi to represent our model over the other
possibilities. Note that, since Voronoi cells can be discontinuous at the
edges, which could produce spurious reflections in FWI. We apply smoothing to
the velocity model before waveform simulations.

##  4 Conclusions

In this paper, we demonstrated the RJHMC algorithm and its application to a 2D
FWI problem using the Marmousi model. We represented our 2D model using
certain nuclei points. RJHMC algorithm solves the 2D FWI problem in a trans-
dimensional way. It solves for the number of nuclei as well as the locations
and velocities of the nuclei. The data dictate entirely the number of nuclei
required, and is mainly determined automatically. The maximum number of nuclei
possible is the dimension of the velocity grid, i.e., $251\times 767=192517$,
but in the example, using only $\approx 15000-16000$ nuclei points, RJHMC is
able to reconstruct the model. Using the gradient information from HMC, the
algorithm can make some large jumps in the model space and efficiently sample
the model space. RJHMC successfully reconstructed the complete model with an
optimal number of model parameters as needed by the data. Multiple models
sampled from the PPD can quantify uncertainty in the estimation of the P-wave
velocity using 2D FWI.

## References

  * Agostinetti et al., (2015) Agostinetti, N. P., G. Giacomuzzi, and A. Malinverno, 2015, Local three-dimensional earthquake tomography by trans-dimensional monte carlo sampling: Geophysical Journal International, 201, 1598–1617.
  * Agostinetti and Malinverno, (2010) Agostinetti, N. P., and A. Malinverno, 2010, Receiver function inversion by trans-dimensional monte carlo sampling: Geophysical Journal International, 181, 858–872.
  * Al-Awadhi et al., (2004) Al-Awadhi, F., M. Hurn, and C. Jennison, 2004, Improving the acceptance rate of reversible jump mcmc proposals: Statistics & probability letters, 69, 189–198.
  * Aster et al., (2011) Aster, R. C., B. Borchers, and C. H. Thurber, 2011, Parameter estimation and inverse problems: Academic Press.
  * Bernardo and Smith, (2001) Bernardo, J. M., and A. F. Smith, 2001, Bayesian theory.
  * Betancourt, (2017) Betancourt, M., 2017, A conceptual introduction to hamiltonian monte carlo: arXiv preprint arXiv:1701.02434.
  * Bodin and Sambridge, (2009) Bodin, T., and M. Sambridge, 2009, Seismic tomography with the reversible jump algorithm: Geophysical Journal International, 178, 1411–1436.
  * Bodin et al., (2009) Bodin, T., M. Sambridge, and K. Gallagher, 2009, A self-parametrizing partition model approach to tomographic inverse problems: Inverse Problems, 25, 055009.
  * Bodin et al., (2012) Bodin, T., M. Sambridge, N. Rawlinson, and P. Arroucau, 2012, Transdimensional tomography with unknown data noise: Geophysical Journal International, 189, 1536.
  * Bosch et al., (2010) Bosch, M., T. Mukerji, and E. F. Gonzalez, 2010, Seismic inversion for reservoir properties combining statistical rock physics and geostatistics: A review: Geophysics.
  * Buland and Omre, (2003) Buland, A., and H. Omre, 2003, Bayesian linearized avo inversion: Geophysics, 68, 185–198.
  * Cary and Chapman, (1988) Cary, P., and C. Chapman, 1988, Automatic 1-d waveform inversion of marine seismic refraction data: Geophysical Journal International, 93, 527–546.
  * Dadi et al., (2016) Dadi, S., R. Gibson, and K. Wang, 2016, Velocity log upscaling based on reversible jump markov chain monte carlo simulated annealingrjmcmc-based sonic log upscaling: Geophysics, 81, R293–R305.
  * Datta et al., (2018) Datta, D., D. Appelhans, C. Evangelinos, and K. Jordan, 2018, An asynchronous two-level checkpointing method to solve adjoint problems on hierarchical memory spaces: Computing in Science & Engineering, 20, 39–55.
  * Datta and Sen, (2016) Datta, D., and M. K. Sen, 2016, Estimating a starting model for full-waveform inversion using a global optimization method: Geophysics, 81, R211–R223.
  * Dettmer et al., (2016) Dettmer, J., R. Hawkins, P. R. Cummins, J. Hossen, M. Sambridge, R. Hino, and D. Inazu, 2016, Tsunami source uncertainty estimation: The 2011 japan tsunami: Journal of Geophysical Research: Solid Earth, 121, 4483–4505.
  * Dosso et al., (2014) Dosso, S. E., J. Dettmer, G. Steininger, and C. W. Holland, 2014, Efficient trans-dimensional bayesian inversion for geoacoustic profile estimation: Inverse Problems, 30, 114018.
  * Duane et al., (1987) Duane, S., A. D. Kennedy, B. J. Pendleton, and D. Roweth, 1987, Hybrid monte carlo: Physics letters B, 195, 216–222.
  * Fichtner et al., (2018) Fichtner, A., A. Zunino, and L. Gebraad, 2018, Hamiltonian monte carlo solution of tomographic inverse problems: Geophysical Journal International, 216, 1344–1363.
  * Green, (1995) Green, P. J., 1995, Reversible jump markov chain monte carlo computation and bayesian model determination: Biometrika, 82, 711–732.
  * Green and Mira, (2001) Green, P. J., and A. Mira, 2001, Delayed rejection in reversible jump metropolis–hastings: Biometrika, 88, 1035–1053.
  * Hawkins et al., (2019) Hawkins, R., T. Bodin, M. Sambridge, G. Choblet, and L. Husson, 2019, Trans-dimensional surface reconstruction with different classes of parameterization: Geochemistry, Geophysics, Geosystems.
  * Hirata and Matsu’ura, (1987) Hirata, N., and M. Matsu’ura, 1987, Maximum-likelihood estimation of hypocenter with origin time eliminated using nonlinear inversion technique: Physics of the earth and planetary interiors, 47, 50–61.
  * Karimi et al., (2010) Karimi, O., H. Omre, and M. Mohammadzadeh, 2010, Bayesian closed-skew Gaussian inversion of seismic AVO data for elastic material properties: GEOPHYSICS, 75, R1–R11.
  * Kotsi et al., (2018) Kotsi, M., A. Malcolm, and G. Ely, 2018, 4d full-waveform metropolis hastings inversion using a local acoustic solver, in SEG Technical Program Expanded Abstracts 2018: Society of Exploration Geophysicists, 5323–5327.
  * Lailly and Bednar, (1983) Lailly, P., and J. Bednar, 1983, The seismic inverse problem as a sequence of before stack migrations: Conference on inverse scattering: theory and application, Siam Philadelphia, PA, 206–220.
  * Malinverno, (2002) Malinverno, A., 2002, Parsimonious bayesian markov chain monte carlo inversion in a nonlinear geophysical problem: Geophysical Journal International, 151, 675–688.
  * Malinverno and Leaney, (2000) Malinverno, A., and S. Leaney, 2000, A monte carlo method to quantify uncertainty in the inversion of zero-offset vsp data: 70th SEG Annual Meeting Expanded Abstracts, 2393–2396.
  * Malinverno and Leaney, (2005) Malinverno, A., and W. S. Leaney, 2005, Monte-carlo bayesian look-ahead inversion of walkaway vertical seismic profiles: Geophysical prospecting, 53, 689–703.
  * Mosegaard and Tarantola, (1995) Mosegaard, K., and A. Tarantola, 1995, Monte carlo sampling of solutions to inverse problems: Journal of Geophysical Research: Solid Earth, 100, 12431–12447.
  * Neal, (2011) Neal, R. M., 2011, Mcmc using hamiltonian dynamics, in Handbook of Markov Chain Monte Carlo: Chapman & Hall/CRC.
  * Okabe, (1992) Okabe, A., 1992, Spatial tessellations: Wiley Online Library.
  * Olugboji et al., (2017) Olugboji, T., V. Lekic, and W. McDonough, 2017, A statistical assessment of seismic models of the us continental crust using bayesian inversion of ambient noise surface wave dispersion data: Tectonics, 36, 1232–1253.
  * Pratt and Worthington, (1990) Pratt, R. G., and M. Worthington, 1990, Inverse theory applied to multi-source cross-hole tomography. part 1: Acoustic wave-equation method: Geophysical prospecting, 38, 287–310.
  * Ray and Key, (2012) Ray, A., and K. Key, 2012, Bayesian inversion of marine csem data with a trans-dimensional self parametrizing algorithm: Geophysical Journal International, 191, 1135–1151.
  * Ray et al., (2016) Ray, A., A. Sekar, G. M. Hoversten, and U. Albertin, 2016, Frequency domain full waveform elastic inversion of marine seismic data from the alba field using a bayesian trans-dimensional algorithm: Geophysical Journal International, 205, 915–937.
  * Sambridge et al., (1995) Sambridge, M., J. Braun, and H. McQueen, 1995, Geophysical parametrization and interpolation of irregular data using natural neighbours: Geophysical Journal International, 122, 837–857.
  * Sambridge et al., (2006) Sambridge, M., K. Gallagher, A. Jackson, and P. Rickwood, 2006, Trans-dimensional inverse problems, model comparison and the evidence: Geophysical Journal International, 167, 528–542.
  * Sen and Biswas, (2017) Sen, M. K., and R. Biswas, 2017, Transdimensional seismic inversion using the reversible jump hamiltonian monte carlo algorithm: GEOPHYSICS, 82, R119–R134.
  * Sen and Roy, (2003) Sen, M. K., and I. G. Roy, 2003, Computation of differential seismograms and iteration adaptive regularization in prestack waveform inversion: Geophysics, 68, 2026–2039.
  * Sen and Stoffa, (1991) Sen, M. K., and P. L. Stoffa, 1991, Nonlinear one-dimensional seismic waveform inversion using simulated annealing: Geophysics, 56, 1624–1638.
  * Sen and Stoffa, (1996) ——–, 1996, Bayesian inference, gibbs’ sampler and uncertainty estimation in geophysical inversion: Geophysical Prospecting, 44, 313–350.
  * Sisson, (2005) Sisson, S. A., 2005, Transdimensional markov chains: A decade of progress and future perspectives: Journal of the American Statistical Association, 100, 1077–1089.
  * Tarantola, (1984) Tarantola, A., 1984, Inversion of seismic reflection data in the acoustic approximation: Geophysics, 49, 1259–1266.
  * Tarantola, (1986) ——–, 1986, A strategy for nonlinear elastic inversion of seismic reflection data: Geophysics, 51, 1893–1903.
  * Tarantola and Valette, (1982) Tarantola, A., and B. Valette, 1982, Generalized nonlinear inverse problems solved using the least squares criterion: Reviews of Geophysics, 20, 219–232.
  * Van Gunsteren and Berendsen, (1988) Van Gunsteren, W. F., and H. J. C. Berendsen, 1988, A Leap-frog Algorithm for Stochastic Dynamics: Molecular Simulation, 1, 173–185.
  * Virieux and Operto, (2009) Virieux, J., and S. Operto, 2009, An overview of full-waveform inversion in exploration geophysics: Geophysics, 74, WCC1–WCC26.
  * Zhu et al., (2016) Zhu, H., S. Li, S. Fomel, G. Stadler, and O. Ghattas, 2016, A bayesian approach to estimate uncertainty for full-waveform inversion using a priori information from depth migration: GEOPHYSICS, 81, R307–R323.