Trans-dimensional Bayesian Inversion Research Articles

Uncertainty Quantification in Radial Anisotropy Models Based on Transdimensional Bayesian Inversion of Receiver Functions and Surface-Wave Dispersion: Case Study Sri Lanka

ABSTRACT In geophysical inference problems, quantification of data uncertainties is required to balance the data-fitting ability of the model and its complexity. The transdimensional hierarchical Bayesian approach is a powerful tool to evaluate the level of uncertainty and determine the complexity of the model by treating data errors and model dimensions as unknown. In this article, we take account of the uncertainty through the whole procedure, thus developing a two-step fully Bayesian approach with coupled uncertainty propagation to estimate the crustal isotropic and radial anisotropy (RA) model based on Rayleigh and Love dispersion as well as receiver functions (RFs). First, 2D surface-wave tomography is applied to determine period-wise ambient noise phase velocity maps and their uncertainty for Rayleigh and Love waves. Probabilistic profiles of the isotropic average VS and RA as a function of depth are then derived at station sites by inverting the local surface-wave dispersion and model errors and RFs jointly. The workflow is applied to a temporary seismic broadband array covering all of Sri Lanka. The probabilistic results enable us to effectively quantify the uncertainty of the final RA model and provide robust inferences. The shear-wave velocity results show that the range of Moho depths is between 30 and 40 km, with the thickest crust (38–40 km) beneath the central Highland Complex. Positive RA (VSH>VSV) observed in the upper crust is attributed to subhorizontal alignment of metamorphic foliation and stretched layers resulting from deformation. Negative RA (VSV>VSH) in the midcrust of central Sri Lanka may indicate the existence of melt inclusions and could result from the uplift and folding process. The positive RA in the lower crust could be caused by crustal channel flow in a collision orogeny.

Shallow crustal imaging beneath NW Indian terrane from teleseismic P-wave coda using a Bayesian approach

Deciphering the crustal architecture of the Cretaceous Barmer-Sanchore Basin (BSB) and Neoproterozoic provinces in the NW Indian Shield presents significant challenges due to extensive sediment deposition and widespread felsic igneous rocks. To address such intricacies, we employ a transdimensional Bayesian joint inversion of teleseismic P-wave polarizations and receiver functions to resolve the detailed structure of the upper 10 km of the crust beneath 29 seismic stations deployed in the study area. The new data is analyzed separately within each of the three different backazimuthal groups to identify any strong variations in the shallow subsurface structure and derive a more robust, averaged model for every station site. The resulting 1-D models reveal heterogeneous sediment cover in BSB and Marwar Basin, with respective thicknesses ranging from ∼1.3–5.5 km and 1–1.5 km. Remarkably low shear-wave velocities (Vs < 1.0 km/s) of sediments with increased Vp/Vs (>2.3) are observed in BSB and adjacent areas, which are related to highly unconsolidated Quaternary-Tertiary deposits. The Marwar Basin is imaged to have 0.3–0.5 km of uncompacted sediments that overlie 0.5–1 km thick Marwar Supergroup sequence. The Erinpura Granites (EG) and Malani Igneous Province (MIP) are characterized by a variably thick sediment blanket, typically ranging within 1 km in thickness. The Vs of the sediments in these regions exceeds 2.4 km/s that possibly corresponds to Mesozoic and older sediments. The basement beneath EG exhibits Vs faster than 3.4 km/s, consistent with moderately metamorphosed granites. In MIP, the basement is characterized by non-uniform shear velocities, suggesting the influence of intrusions linked to magmatic episodes in NW India.

Trans-dimensional joint inversion of modal amplitude and arrival- time data for geoacoustic properties in layered muddy sediments

Long-range sound fields in a shallow-water waveguide can be expressed as a summation of several dispersive propagating modes. Modal arrival times as a function of frequency, extracted via warping time-frequency analysis, have been used by many researchers to estimate sound speeds and densities for layered models of the seabed sediments. This paper explores joint inversion of modal amplitude spectra together with arrival times to also estimate sediment-layer attenuations. To consider the common case in underwater-acoustic applications where the source amplitude spectrum is unknown, we include this spectrum as unknown parameters in the inversion, formulated implicitly within the likelihood function. Trans-dimensional Bayesian inversion is applied to estimate marginal probability profiles of the geoacoustic properties, providing quantitative uncertainties. The data considered here were collected on the New England Mud Patch during the 2017 Seabed Characterization Experiment. [Work supported by the Office of Naval Research.]

Parallel tempering in trans-dimensional Bayesian inversion for seabed geoacoustic models with many parameters per layer

Trans-dimensional (trans-D) Bayesian inversion is a powerful tool for estimating seabed geoacoustic models from ocean-acoustic data, combining quantitative model selection with parameter/uncertainty estimation. The approach applies reversible-jump Markov-chain Monte Carlo methods to sample probabilistically over the number of seabed layers and the corresponding geoacoustic parameters for each layer. Layers are added and removed during sampling, referred to as birth and death moves, respectively, changing the dimension of the model. However, the probability of accepting birth and death moves can approach zero for formulations that include many parameters per layer. This paper considers the use of parallel tempering to mitigate this degradation in efficiency. Parallel tempering employs a series of interacting Markov chains with successfully-relaxed acceptance criteria, achieved by raising the likelihood to powers of 1/T, with T greater than or equal to 1 referred to as the sampling temperature. While only the T = 1 chain provides unbiased sampling, probabilistic interchange between chains provides a robust ensemble sampler that mixes more readily over the trans-D model space. The approach is illustrated for wide-angle reflection-coefficient inversion including compressional and shear parameters in the seabed model, resulting in a total of 5 unknown parameters per layer.

Integrating Seismic Methods for Characterizing and Monitoring Landslides: A Case Study of the Heinzenberg Deep-Seated Gravitational Slope Deformation (Switzerland)

While geodetic measurements have long been used to assess landslides, seismic methods are increasingly recognized as valuable tools for providing additional insights into subsurface structures and mechanisms. This work aims to characterize the subsurface structures of the deep-seated gravitational slope deformation (DSGSD) at Heinzenberg through the integration of active and passive seismic measurements. Seismic techniques can hereby deliver additional information on the subsurface structure and mechanisms involved, e.g., the degree of rock mass degradation, the resonant frequencies of the potentially unstable compartments, and the local fracture network orientations that are influenced by wavefield polarization. By employing advanced methods such as H/V analysis, site-to-reference spectral ratios, polarization analysis, surface wave analysis, and the joint multizonal transdimensional Bayesian inversion of velocity structures, we establish a comprehensive baseline model of the landslide at five selected sites. This baseline model shall help identify potential changes after the refilling of Lake Lüsch, which started in 2021. Our results reveal the rupture surface of the DSGSD at various depths ranging from 30 m at the top to over 90 m in the middle of the slope. Additionally, we estimate key parameters including the shear wave velocities of the different rock masses. The 2D geophysical profiles and rock mass properties contribute to the understanding of the subsurface geometry, geomechanical properties, and potential water pathways. This study demonstrates the significance of integrating seismic methods with traditional geodetic measurements and geomorphologic analysis techniques for a comprehensive assessment of landslides, enhancing our ability to monitor and mitigate hazardous events.

Reconstruction of Multiple Target Bodies Using Trans-Dimensional Bayesian Inversion With Different Constraints

Reconstruction of Multiple Target Bodies Using Trans-Dimensional Bayesian Inversion With Different Constraints

Transdimensional Bayesian Inversion and Uncertainty Estimation of DC Resistivity Data for Undercover Cavern Characterization

The presence of subsurface karst caverns could pose significant threat to highway construction safety, and the knowledge about the locations and volumes of the caverns is crucial for evaluating potentially hazardous areas. In this work, we present a transdimensional Bayesian approach for characterizing undercover caverns by inverting direct current resistivity data., This approach provides a robust measure of uncertainty associated with the resulting resistivity model. We evaluate the effectiveness of the algorithm on a synthetic 2D resistivity model. The results demonstrate that quantitative uncertainty information on the derived model parameters is vital for reliable data interpretation and facilitates decision making.

Trans-Dimensional Bayesian Inversion of Rayleigh Wave Dispersion Curves for Shallow Surface Two-Dimensional Imaging

In this study, we achieve 2-D spatial structure directly from dispersion measurements using the reversible jump Markov chain Monte Carlo (rj-McMC) algorithm with a fully 2-D model parameterization. New parsimonious parameterization with low dimensional representation is introduced, which enable us to accurately reproduce the wave velocity field and the associated uncertainties. Field data test demonstrates that the proposed method estimates directly the shear-wave velocity profile with better lateral resolution. The associated uncertainties provide directly interpretable uncertainty on structures of interest.

Depth and frequency dependence of geoacoustic properties on the New England Mud Patch from reflection coefficient inversiona).

Muddy sediments cover significant portions of continental shelves, but their physical properties remain poorly understood compared to sandy sediments. This paper presents a generally applicable model for sediment-column structure and variability on the New England Mud Patch (NEMP), based on trans-dimensional Bayesian inversion of wide-angle, broadband reflection-coefficient data in this work and in two previously published reflection-coefficient inversions at different sites on the NEMP. The data considered here include higher frequencies and larger bandwidth and cover lower reflection grazing angles than the previous studies, hence, resulting in geoacoustic profiles with significantly better structural resolution and smaller uncertainties. The general sediment-column structure model includes an upper mud layer in which sediment properties change slightly with depth due to near-surface processes, an intermediate mud layer with nearly uniform properties, and a geoacoustic transition layer where properties change rapidly with depth (porosity decreases and sound speed, density, and attenuation increase) due to increasing sand content in the mud above a sand layer. Over the full frequency band considered in the new and two previous data sets (400-3125 Hz), there is no significant sound-speed dispersion in the mud, and attenuation follows an approximately linear frequency dependence.

A general model for sediment-column structure on the New England Mud Patch from Bayesian geoacoustic inversion of seabed reflection data

Muddy sediments cover significant portions of continental shelves, but their physical properties remain poorly understood compared to sandy sediments. To explore the spatial and frequency dependencies of mud properties, wide-angle seabed reflection coefficients versus grazing angle and frequency were measured on the New England Mud Patch (NEMP) during the 2017 Seabed Characterization Experiment. This paper presents trans-dimensional Bayesian inversion of reflection coefficients within a frequency band of 1–3 kHz and an angular range of ∼15–25° to obtain geoacoustic profiles and associated uncertainties, as well as frequency dependencies of sound speed and attenuation. The estimated geoacoustic profiles are similar to those from previous inversions of reflection-coefficient data at lower frequencies (0.4–1.3 kHz) collected at two different sites on the NEMP. Based on the inversion results at all three sites, a general interpretive model for sediment-column structure and variability is synthesized for the NEMP. This model includes an upper mud layer in which sediment properties change slightly with depth due to near-surface processes, an intermediate mud layer with uniform properties, and a transition layer where properties change rapidly with depth due to increasing sand content in the mud above a sand layer. [Work supported by the Office of Naval Research.]

Relation Between Oceanic Plate Structure, Patterns of Interplate Locking and Microseismicity in the 1922 Atacama Seismic Gap

AbstractWe deployed a dense geodetic and seismological network in the Atacama seismic gap in Chile. We derive a microseismicity catalog of &gt;30,000 events, time series from 70 GNSS stations, and utilize a transdimensional Bayesian inversion to estimate interplate locking. We identify two highly locked regions of different sizes whose geometries appear to control seismicity patterns. Interface seismicity concentrates beneath the coastline, just downdip of the highest locking. A region with lower locking (27.5°S–27.7°S) coincides with higher seismicity levels, a high number of repeating earthquakes and events extending toward the trench. This area is situated where the Copiapó Ridge is subducted and has shown previous indications of both seismic and aseismic slip, including an earthquake sequence in 2020. While these findings suggest that the structure of the downgoing oceanic plate prescribes patterns of interplate locking and seismicity, we note that the Taltal Ridge further north lacks a similar signature.

Generating 1D depth profile of shear-wave velocity from ambient noise cross-correlation: application to Jakarta array

In the past decade, cross-correlations of ambient seismic noise have been exploited in various applications to model the shallow-to-deep structure of Earth’s interior through tomographic inversions. The stack of cross-correlations between a 2-station pair represents empirical Green’s function and comprises the information of the subsurface structure between those stations. In practice, noise correlation function (NCF) is analyzed to reconstruct surface wave group or phase velocity dispersion; then, the dispersion data is used to model shear-wave velocity (Vs). This study presents a case for temporary seismic networks deployed in the Jakarta Basin; we applied a two-step routine to obtain a representative 1D Vs profile beneath an array. First, we extracted our array’s average phase velocity dispersion based on the relationship between NCF’s spectra and the Bessel function. Then, we invert for the 1D depth profile of Vs using a transdimensional Bayesian inversion to allow for exploring a number of layers in parameterizations. We successfully generate a 1D Vs profile up to 5 km depth reflecting the regional stratigraphy of the Jakarta Basin. In general, a sedimentary basin fill covers the area reaching a depth of 650 m. We suggest that this simple routine can be undertaken for other ambient noise cross-correlation cases; such a 1D depth profile would be beneficial to be used as a reference model.

A New Shear‐Velocity Model of Continental Australia Based on Multi‐Scale Surface‐Wave Tomography

AbstractThe Australian Seismological Reference Model (AuSREM) represents a state‐of‐the‐art geophysical synthesis of the Australian continent. To date, its shear‐wave component has limited resolution at lower‐crustal to uppermost‐mantle depths, where it is mainly constrained by sparse measurements collected at the local scale. In this study, we compile a large data set of surface‐wave phase velocities based on seismic ambient noise and teleseismic earthquakes, to produce Rayleigh and Love phase‐velocity maps of continental Australia in a broad period range (4–200 s). Via transdimensional Bayesian inversion, we translate the phase‐velocity maps into a 3‐D shear‐wave velocity model extending to 300 km depth. Owing to the unprecedented seismic coverage and to the joint use of ambient‐noise and teleseismic data, the retrieved model fills a tomographic gap in the known shear‐wave velocity structure of the continent, comprising lower‐crustal to uppermost mantle depths. Consistent with AuSREM, strong velocity heterogeneities in our model highlight the (faster) cratonic blocks and the (slower) sedimentary basins at upper‐crustal depths. At mantle depths, the most prominent feature of the continent is a large‐scale eastward decrease in shear‐wave velocity. We interpret our observations in light of the relevant literature, and produce depth maps of the Moho and lithosphere‐asthenosphere boundary (LAB). Notably, our LAB proxy features a stripe of thicker lithosphere extending to the east coast, which is not visible in the AuSREM LAB model. This observation supports the idea that lateral variations in lithospheric thickness control both the composition and volume of surface volcanism in eastern Australia.

The crust-mantle transition beneath Northeast China from P–wave receiver functions

Introduction: The continental lithosphere of the Northeast China Plain (NECP) is probed using P–wave receiver function analysis of passive seismic data. The NECP is much–discussed as it includes different geological provinces of varying tecto–magmatic origin and dates back from Archean to Holocene. Quantifying the tectonic and magmatic influences on the structure and composition of the lithosphere puts important constraints on evolution of NECP. For this, we explore 75 sites across NECP using receiver function analysis.Methods: A recently developed technique of inverting for 1–D S–wave velocity profiles beneath seismic stations that is based on the principles of Bayesian statistics (hierarchical transdimensional Bayesian Inversion; HTBI) is applied to receiver functions from the NECP. In addition, an improved crustal thickness–compressional to shear wave velocity ratio (H–κ) analysis was conducted to retrieve the crustal thickness and Vp/Vs ratio of the region. These estimated point measurements are integrated and systematically studied for a regional view of the current crustal architecture.Results and Discussion: We observe a laterally varying and highly complex lithosphere beneath the NECP. A shallower crust–mantle transition (≤32 km) characterises the Precambrian North China craton and Late Mesozoic–Cenozoic Songliao Basin from the adjacent Central Asian Orogenic Belt and the Changbaishan Volcanic field (35–40 km). Beneath the latter, low Vp/Vs ratios (~1.65) are obtained, whereas all other regions feature ratios in excess of 1.75. Multiple velocity gradients are observed at crustal depths within the craton, in contrast to the adjacent orogen, which indicates a higher degree of crustal complexity of the former. The width of the crust–mantle transition across the NECP is found to be mainly intermediate (2–7 km) and occasionally sharp (≤2 km). From our observations, we infer that there is a substantial difference between the eastern North China Craton’s lithospheric architecture and the rest of the NECP, with most of the NECP exhibiting more complexity than previously reported.

3-D shallow shear velocity structure of the Jakarta Basin from transdimensional ambient noise tomography.

Situated on the northern coast of the Indonesian island of Java, Jakarta and its metropolitan area (Greater Jakarta) are subject to significant earthquake hazards from a subduction zone south of Java and nearby active crustal faults. The seismic risk may be even higher because Greater Jakarta resides on a sedimentary basin filled with thick Pliocene-Pleistocene sediments. A comprehensive study of Jakarta Basin's properties and geometry is important for creating robust seismic hazard and risk assessments. The main objective of this study is to develop a 3-D model of Jakarta Basin's shallow shear-wave velocity (VS ) structure and improve on previous models that did not cover the basin edge due to the extent of data coverage. Between April and October 2018, we deployed a new temporary seismic network to extend the spatial coverage beyond that of a previous deployment in 2013, and sampled 143 locations through sequential deployments of 30 broad-band sensors covering Jakarta and its adjacent areas. We conducted a 2-stage transdimensional Bayesian inversion of Rayleigh wave phase velocity dispersion curves derived from seismic noise. To begin, we applied tomography and constructed 2-D phase velocity maps for periods 1-5s. Then, at each point in a regular grid defined on these maps, we invert each dispersion curve into 1-D depth profiles of VS . Finally, these profiles at gridpoints with ∼2km spacing are interpolated to form a pseudo-3-D VS model. Our results reveal the edge of the Pliocene-Pleistocene sediments along the south. Also, we resolve a basement offset across south Jakarta that we suggest may be related to the western extension of the Baribis Fault (alternatively, the West Java Backarc Thrust). We recommend using this 3-D model of the Jakarta Basin for scenario earthquake ground motion simulations. Such simulations would help establish how important it might be to re-assess seismic hazard and risk in Greater Jakarta so that basin resonance and amplification are considered.

Trans-dimensional Bayesian joint inversion of magnetotelluric and geomagnetic depth sounding responses to constrain mantle electrical discontinuities

SUMMARYJoint inversion of magnetotelluric (MT) and geomagnetic depth sounding (GDS) responses can constrain the crustal and mantle conductivity structures. Previous studies typically use either deterministic inversion algorithms that provide limited information on model uncertainties or using stochastic inversion algorithms with a predetermined number of layers that is generally not known a priori. Here, we present a new open-source Bayesian framework for the joint inversion of MT and GDS responses to probe 1D layered Earth’s conductivity structures. Within this framework, model uncertainties can be accurately estimated by generating numerous models that fit the observed data. A trans-dimensional Markov Chain Monte Carlo (MCMC) method is employed to self-parametrize the model parameters, where the number of layers is treated as an inversion parameter that is determined automatically by the data. This adaptability can overcome the under or over-parametrization problem and may be able to automatically detect the conductivity discontinuities in the Earth’s interior. To accelerate the computations, a large number of Markov chains with different initial states can be run simultaneously using the MPI parallel technique. Synthetic data sets are used to validate the feasibility of our method and illustrate how separate and joint inversions, as well as various priors affect the posterior model distributions. The trans-dimensional MCMC algorithm is then applied to jointly invert the MT and GDS responses estimated at the Tucson geomagnetic observatory, North America. Our results not only contain model uncertainty estimates but also indicate two distinct conductivity discontinuities at around 85 and 440 km, which are likely related to the lithosphere-asthenosphere boundary and the upper interface of the mantle transition zone, respectively.

Resolving minute temporal seismic velocity changes induced by earthquake damage: the more stations, the merrier?

SUMMARYGround shaking induced by earthquakes often introduces transient changes in seismic velocity monitored with ambient noise. These changes are usually attributed to relaxation behaviour following the coseismic damage in the subsurface and are of relevance for post-seismic hazard mitigation. However, the velocity evolution associated with this phenomenon can occur at very small timescales and amplitudes that are not resolved with seismic interferometry and are therefore challenging to link to laboratory experiments. A way to improve the temporal resolution of the velocity time-series is to test whether the estimation of the relative seismic velocity changes dv/v obeys the ergodic hypothesis in which the joint use of colocated stations would lead to better resolved measurements. In this study, we present results from a dense seismic array that was deployed for 2 weeks at the remarkable Patache site in Chile. Thanks to high temporal averaging capabilities, we are able to resolve seismic velocity changes in the 3–6 Hz frequency band at a 10-min resolution around the occurrence of a moderate earthquake (PGV ∼1 cm s–1). We report a velocity drop of ∼0.4 per cent in the first 10 min after ground shaking. Half of this initial drop was recovered within the 2 following days. The shape of the recovery follows a log-linear shape over the whole observed recovery phase, analogous to slow dynamics experiments. When normalized by the total amount of processed data, we show that the ergodic hypothesis almost perfectly holds in our network: the dv/v signal-to-noise ratio (SNR) obtained when averaging a few observation with large stacking durations for the correlation functions is almost equal to the SNR when using a large number of observations with small stacking durations. To understand if the ergodicity is linked to a particular site property, we use the array capabilities to identify the surf at the shoreline as the source of the noise and to derive a 1-D shear velocity profile with the focal spot imaging technique and a transdimensional Bayesian inversion framework. The inversion shows that hard rocks lie close to the surface indicating that this material hosts the observed shallow velocity changes. We discuss our high-resolution measurements and attribute them to a stable noise source excited by the shore, the ergodicity property and an ideal subsurface structure. Finally, we discuss the effect of moderate earthquakes on subsurface damage and the potential relaxation processes in hard rocks.

Inversion of ocean acoustic modal-dispersion and amplitude data for seabed geoacoustic models including attenuation

The long-range acoustic field due to an initial sound source in a shallow-water waveguide can be expressed as the sum of a number of dispersive, propagating modes. Modal-dispersion data (arrival time as a function of frequency), as extracted with warping time-frequency analysis, have been inverted to estimate seabed sediment sound-speed and density profiles. In this work, measurements of modal amplitudes as a function of frequency are combined with modal-dispersion data in a joint (simultaneous) inversion which also provides sensitivity to frequency-dependent sound attenuation coefficients in the sediment layers. Since the source signature (spectrum) is often unknown in practice, the relative amplitudes between modes can be considered as data, or the source spectrum can itself be included as additional unknowns in the inversion. Here, trans-dimensional Bayesian inversion is applied to modal-dispersion and amplitude data to estimate marginal probability profiles for sediment geoacoustic properties, including attenuation. The signal processing and inversion methods are demonstrated using modal-dispersion and amplitude data extracted from acoustic measurements made by a hydrophone-equipped underwater glider during the 2017 Seabed Characterization Experiment (Work supported by the Office of Naval Research).

Bayesian geoacoustic inversion of 1-2 kHz seabed reflection data for layered muddy sediments

This paper presents trans-dimensional Bayesian geoacoustic inversion of seabed reflection data for sub-bottom geoacoustic profiles and associated uncertainty estimation at the New England Mud Patch. The data considered here are wide-angle seabed reflection coefficients as a function of grazing angle and frequency measured during the 2017 Seabed Characterization Experiment. Chirp pulses over a frequency band of 1–6 kHz were transmitted by an omnidirectional acoustic source towed by a research vessel at a speed of about 4 knots and recorded at a bottom-moored hydrophone. High signal-to-noise-ratio reflection coefficients from 1–2 kHz and angular coverage of ∼15–25° are considered here for geoacoustic inversion. This frequency range is higher than for previous reflection-coefficient data sets on the Mud Patch. The angular range, although relatively narrow, includes strong Bragg resonances which provide information on the sediment layering properties. The inversion applies the viscous grain-shearing sediment acoustics model, which provides dispersive (frequency-dependent) results for sound speeds and attenuations. The seabed structure estimated here is compared to previous inversion results and to core measurements in the vicinity. [Work supported by the Office of Naval Research.]

A Julia software package for transdimensional Bayesian inversion of electromagnetic data over horizontally stratified media

A quantitative assessment of model parameter uncertainty is vital for a reliable interpretation of electromagnetic (EM) data due to the nonuniqueness inherent to EM inverse problems. Conventional gradient-based inversion approaches typically produce a single preferred model with limited information about parameter uncertainty. The inverse problems can be alternatively postulated into a sampling-based Bayesian inference framework where the solution is represented by a posterior probability distribution of model parameters, which can provide an effective way to rigorously estimate parameter uncertainty related to the recovered solution. We have implemented a Bayesian inversion framework for probabilistic inversion of EM data, which is an open-source software package implemented in the Julia programming language. The key feature of the framework is that it allows the model complexity to be adaptively adjusted to an appropriate level compatible with the data by implementing a reversible jump Markov chain Monte Carlo algorithm, thus allowing the data to infer the appropriate level of model complexity and associated parameter uncertainty. The authors have elaborated the structure of the package with a focus on code modularity and extensibility. Finally, the authors determine the capacity and versatility of the software package through three synthetic examples that simulate different EM scenarios. The inversion results demonstrate the performance and inherent resolving abilities of different EM surveys for conductive and/or resistive structures. In addition, the model ensemble produced by the transdimensional Bayesian inversion conveys a wealth of information. Many important quantities for data interpretation such as parameter uncertainty and correlation between model parameters of interest can be measured by exploring the model ensemble.