Monte Carlo Inversion Research Articles

Ice-sheet burial and erosion inferred from cosmogenic nuclide bedrock depth profiles: Implications for the glaciation history of northeastern Fennoscandia

In this study, we evaluate the application of shallow (<2.5 m) cosmogenic depth profiles in bedrock to constrain long-term ice-burial and erosion histories. Using Markov Chain Monte Carlo inversion modelling on a series of synthetic scenarios, we demonstrate that cosmogenic 10Be and 26Al profiles provide more robust constraints on ice-burial duration and erosion histories than surface samples alone, particularly when erosion rates are low (<5–10 m Myr−1) and/or non-steady. We apply this method to new depth profiles of 10Be and 26Al measurements from two tors in the Parkajoki region in northeastern Sweden. Our results indicate erosion depths of ∼2–10 m and ice burial for ∼20–35% of the time since 500 ka. These estimates imply more erosion and less ice burial than previously inferred from the same tors. However, by re-assessing the extent of ice cover during the Weichselian from existing records, we show that some cosmogenic nuclide inheritance predates the penultimate glacial maximum (Late Saalian), implying limited glacial erosion in the Parkajoki region during the last glacial cycle.

An early giant planet instability recorded in asteroidal meteorites

Giant planet migration appears widespread among planetary systems in our Galaxy. However, the timescales of this process, which reflect the underlying dynamical mechanisms, are not well constrained, even within the Solar System. As planetary migration scatters smaller bodies onto intersecting orbits, it would have resulted in an epoch of enhanced bombardment in the Solar System’s asteroid belt. Here, to accurately and precisely quantify the timescales of migration, we interrogate thermochronologic data from asteroidal meteorites, which record the thermal imprint of energetic collisions. We present a database of 40K–40Ar system ages from chondrite meteorites and evaluate it with an asteroid-scale thermal code coupled to a Markov chain Monte Carlo inversion. Simulations require bombardment to reproduce the observed age distribution and identify a bombardment event beginning 11.3−6.6+9.5Myr\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$11.{3}_{-6.6}^{+9.5}\\, {\\mathrm{Myr}}$$\\end{document} after the Sun formed (50% credible interval). Our results associate a giant planet instability in our Solar System with the dissipation of the gaseous protoplanetary disk.

Markov chain Monte Carlo inversions of the internal rotation of Kepler subgiants

The measurement of the internal rotation of post-main-sequence stars using data from space-based photometry missions has demonstrated the need for an efficient angular momentum transport in stellar interiors. No clear solution has emerged so far, and it remains a challenge for stellar modellers to explain the observed trends. We constrained the shape of the internal rotation profile of six Kepler subgiants that were studied in details in 2014 and also the properties of the missing angular momentum transport process that acts in stellar interiors from Markov chain Monte Carlo (MCMC) inversions of the internal rotation. We applied a new MCMC inversion technique to existing Kepler subgiant targets and tested various shapes of the internal rotation profile of the six subgiants that were observed in 2014. We also constrained the limitations on the number of free parameters that can be used in the MCMC inversion, showing the limitations in the amount of information in the seismic data. First, we show that large-scale fossil magnetic fields are not able to explain the internal rotation of subgiants, similarly to what was determined from detailed studies of Kepler red giants. We are also able to constrain the location of the transition in the internal rotation profile for the most evolved stars in the available set of subgiants. We find that some of them exhibit a transition that is located close to the border of the helium core, but one object exhibit a transition located much higher in radius. We conclude that various processes might be at play that would explain our observations, but a consistent detailed modelling of all available subgiants is required to reveal the physical nature of the angular momentum process, in particular, for the least evolved objects. In addition, it is paramount to increase the number of stars for which these inferences are possible (e.g. with the future PLATO mission) because they play a key role in validating candidates for the transport process.

Some practical applications of shear wave velocity measurements in dense sand

Shear wave velocity (Vs) profiles were obtained using two techniques at a dense sand site in Blessington near Dublin, Ireland. It was shown that both the multichannel analysis of surface waves (MASW) and seismic CPT (SCPTU) techniques provided reliable data. Use was made of a novel approach involving a Monte Carlo inversion of the MASW output to assess the fit between the measured and theoretical surface wave data. Differences between the SCPTU and MASW data are likely to be due to a combination of natural anisotropy in the sand and to uncertainty in the two methods involved. Previous published correlations between CPTU data and Vs also work well for this site. The Vs data was also used to accurately predict (Class 3 prediction) measured settlement data from shallow footing tests at the site.

Mapping structures on the core–mantle boundary using Sdiff postcursors: Part I. Method and Validation

SUMMARY Ultra-low velocity zones (ULVZs) are patches of extremely slow seismic velocities on the core–mantle boundary (CMB). Here, we target them using the postcursors to S core-diffracted phases (Sdiff) caused by ULVZs. We use traveltimes of these postcursors to make probabilistic maps using a reversible-jump Markov chain Monte Carlo inversion setup. For the forward model, we extend 2-D wave front tracking (2DWT) software, previously developed for surface wave multipathing studies, to the CMB. The 2DWT is able to model the full multipathing behaviour of Sdiff postcursors and compute arrival times for a given ULVZ input velocity structure on the order of a few CPU seconds, as opposed to 100s of CPU hours required for 3-D full waveform synthetics. We validate the method using synthetic data sets produced by the 2DWT, as well as 3-D full waveform synthetics, using a parametrisation formed from a collection of ellipses. We also test idealistic data coverage versus a case of more realistic coverage. We show ULVZ size and velocity reduction can typically be well recovered, and our maps show the inherent trade-off between these parameters around the edge of the ULVZ. Our method cannot directly constrain ULVZ height; tests show that it underestimates ULVZ velocity reductions and overestimates ellipticity for thinner ULVZs due to neglecting mantle effects.

Bayesian Monte Carlo Inversion of InSAR Time Series Deformation Induced by Wastewater Injection: A Case Study in West Texas

AbstractWastewater disposal can induce detectable surface uplift, which may cause ground instability and threaten infrastructure. The distributions of local hydro‐geomechanical parameters, especially Young's modulus and hydraulic conductivity, play an essential role in these geohazards. To constrain these parameters, we have inverted spatio‐temporal deformation measured by Interferometric Synthetic Aperture Radar (InSAR) and injection information using a Bayesian Monte Carlo approach with a poroelastic finite element model. Sentinel‐1A/B imagery from 2014 to 2020 is processed to track the spatio‐temporal deformation in Winkler county, West Texas, USA. The posterior distribution of subsurface effective volumes reveals under‐reported volumes in the well near the deformation center. In addition, the inversion results provide better constraints for the parameters than those solely obtained based on the cumulative spatial deformation or temporal development of the deformation center.

Efficient Bayesian inversion of borehole geophysical measurements with a gradient‐based Markov chain Monte Carlo method

AbstractWell logs are geophysical measurements of rock properties acquired continuously along a borehole. Because of the physics underlying their operation, borehole logging instruments perform local spatial averages of rock properties in the vicinity of the borehole, yielding well logs that are often affected by tool design, layer boundaries, borehole environmental conditions and mud‐filtrate invasion. Such environmental effects are ubiquitous and can lead to errors in petrophysical–property estimations from well logs if not accounted for in the interpretation. Separate well‐log inversion mitigates environmental effects by matching well logs with their numerical simulations. The latter simulations rely on specific assumptions about the relative geometry of the borehole and the rocks penetrated by the well. For vertical wells penetrating horizontal layers, it is common to assume piecewise‐constant layer properties that we collectively refer to as the earth model; some of these properties have a direct relationship to well logs (e.g. resistivity, gamma ray, density and acoustic slowness). Well‐log inversion is traditionally approached with deterministic methods such as the Levenberg–Marquardt algorithm to minimize the differences between well logs and their numerical simulations, often without accounting for the relationship between measurement noise and model uncertainty. Bayesian inversion, on the other hand, yields the posterior probability distribution of estimated properties and intrinsically quantifies their uncertainty. However, Bayesian Well‐log inversion usually involves the implementation of Markov chain Monte Carlo sampling, which requires a prohibitive number of forward simulations and is therefore not suitable for rapid petrophysical and/or elastic and mechanical evaluations of rocks. We introduce an efficient Bayesian Well‐log inversion method using a gradient‐based Markov chain Monte Carlo method. Gradient‐based Markov chain Monte Carlo is a group of relatively new Markov chain Monte Carlo algorithms that combine gradient updates with Hessian‐based sampling. gradient‐based Markov chain Monte Carlo draws samples efficiently from the posterior probability distribution using the gradient information to guide samples towards high‐probability regions and the Hessian to approximate the posterior probability distribution locally. We verify the gradient‐based Markov chain Monte Carlo inversion method with the inversion of synthetic well logs, including gamma ray, resistivity, density, neutron porosity, photoelectric factor and compressional/shear‐wave slowness. Results show that gradient‐based Markov chain Monte Carlo decreases the computational cost by more than 90% compared to conventional Markov chain Monte Carlo sampling. Next, gradient‐based Markov chain Monte Carlo inversion is applied to a field example from the Central North Sea, where conventional interpretation methods yield shale concentration, sandstone porosity and water saturation with errors up to 19.2%, 14.9% and 48.8%, respectively. Gradient‐based Markov chain Monte Carlo inversion products satisfactorily reproduce the available measurements with only a modest increase in computational cost compared to deterministic inversion approaches.

Sulfuric acid weathering counteracts CO2 drawdown from silicate weathering in mountainous catchments from southwest China

Sulfide oxidation, coupled with carbonate weathering is thought to be a source of CO2 over geological timescales to offset CO2 drawdown from silicate weathering. Its impact on the carbon budget of weathering at the catchment scale is generally constrained by utilizing riverine water chemistry and the isotope value (δ34S) of sulfate. However, the multiple origins of riverine sulfate and the high variability of sulfide-derived δ34S certainly limit the determination of sulfate sources and thus our understanding of the role of sulfuric acid-induced weathering in the carbon budget. In this study, to quantify the carbon budget between weathering by carbonic and sulfuric acids in mountainous catchments (the Yinjiang River, the Shiqian River, and the Yuqing River) in southwest China, a Monte Carlo inversion model is employed for major element concentrations and sulfate δ34S (δ34SSO4). This inversion quantitatively partitions the cations and sulfate between different sources, and further identifies the pyrite-derived δ34S and the flux ratios of alkalinity to dissolved inorganic carbon (ΔALK/ΔDIC). Inversion results reveal that the majority of major cations (>90 %) is derived from carbonate (calcite and dolomite) dissolution, and that riverine sulfate budget is predominately contributed by pyrite oxidation and precipitation. In the Yinjiang River, the average δ34SSO4 is 8.4 ‰ in the upper and middle reaches, but declines significantly to 0.4 ‰ in the lower reaches. In the Shiqian River and Yuqing River, δ34SSO4 mainly cluster around 3.6–5.7 ‰ and 5.5–8.5 ‰, respectively. Inversion results suggest that the primary controls on δ34SSO4 are the variability of lithology and the variation in the proportional contribution of pyrite oxidation to riverine sulfate. In the Yinjiang River, the fraction of H2SO4-driven weathering is highly variable, and higher in the lower reaches than in the upper reaches, attributing to distinct underlying lithology in the catchment (a widespread distribution of coal-bearing strata). The weathering fluxes with 1 < ΔALK/ΔDIC < 2 evidenced that the negative alkalinity flux from carbonate weathering driven by sulfuric acid from both pyrite oxidation and acid deposition counteracts much of the alkalinity flux from silicate weathering. This is implicated on the timescale between marine carbonate precipitation and sulfate reduction (>5–10 kyr and <10 Myr), indicating that weathering in these catchments acts as short-term sinks of atmospheric CO2 but long-term sources of it. Moreover, reanalysis of prior datasets of karst-dominated catchments in China also highlight the important role of pyrite oxidation and acid precipitation in the carbon budget.

Geothermal heat flux is the dominant source of uncertainty in englacial-temperature-based dating of ice rise formation

Abstract. Ice rises are areas of locally grounded, slow-moving ice adjacent to floating ice shelves. Temperature profiles measured through ice rises contain information regarding changes to their dynamic evolution and external forcings, such as past surface temperatures, past accumulation rates and geothermal heat flux. While previous work has used borehole temperature–depth measurements to infer one or two such parameters, there has been no systematic investigation of parameter sensitivity to the interplay of multiple external forcings and dynamic changes. A one-dimensional vertical heat flow forward model developed here examines how changing forcings affect temperature profiles. Further, using both synthetic data and previous measurements from the Crary Ice Rise in Antarctica, we use our model in a Markov chain Monte Carlo inversion to demonstrate that this method has potential as a useful dating technique that can be implemented at ice rises across Antarctica. However, we also highlight the non-uniqueness of previous ice rise formation dating based on temperature profiles, showing that using nominal values for forcing parameters, without taking into account their realistic uncertainties, can lead to underestimation of dating uncertainty. In particular, geothermal heat flux represents the dominant source of uncertainty in ice rise age estimation. For instance, in Crary Ice Rise higher heat flux values (i.e. about 90 mW m−2) yield grounding timing of 1400 ± 800 years, whereas lower heat flux of around 60 mW m−2 implies earlier ice rise formation and lower uncertainties in the ice rise age estimations (500 ± 250 years). We discuss the utility of this method in choosing future ice drilling sites and conclude that integrating this technique with other indirect dating methods can provide useful constraints on past forcings and changing boundary conditions from in situ temperature–depth measurements.

Joint estimation ofS-wave velocity and damping ratio of the near-surface from active Rayleigh wave surveys processed with a wavefield decomposition approach

SUMMARYThe use of surface wave measurements to derive an S-wave velocity profile of the subsurface has become a widely applied procedure; however, their potential use to reconstruct the S-wave material damping properties of the subsoil is generally overlooked, due to the difficulties in obtaining consistent surface wave amplitude information from field data and translating them into robust estimates of the dissipative properties of the near-surface. In this work, we adapt a wavefield decomposition technique for the processing of elastic surface wave data to the extraction of the complete set of properties of Rayleigh waves generated by a controlled source and propagating in dissipative geomaterials. Retrieved information includes multimodal phase velocity and ellipticity as well as the frequency-dependent attenuation coefficient. We exploit the key advantages of wavefield decomposition processing (joint interpretation of multicomponent recordings, coupled estimation of wave propagation parameters, modelling of multiple superimposing modes) to maximize the robustness of the retrieval of Rayleigh wave properties, especially of the dissipative ones. For the subsequent interpretation of Rayleigh wave dispersion, ellipticity and attenuation data we implement a joint Monte Carlo inversion yielding a coupled estimate of S-wave velocity and damping ratio profile for the subsurface; we incorporate a series of geophysical constraints to narrow down the searched parameter space to realistic soil models. We apply this processing and inversion scheme to a bespoke synthetic data set and to a field survey for the characterization of a strong motion station; in both cases, we succeed in retrieving Rayleigh wave multimodal dispersion, ellipticity and attenuation curves. From the inversion of data from the simulated seismogram we are able to reconstruct the properties of the synthetic model. As for the real case, we determine an S-wave velocity and damping ratio model for the soil column below the station, through which we are able to model the inelastic earthquake local response observed at the site. Basing on the results obtained for the real case, we argue that one of the advantages brought by our processing method—the possibility to process active Rayleigh wave data acquired by a 2-D array illuminated by different source positions—may play a key role in allowing to retrieve dissipative properties of the near-surface closer to the material damping of the soil materials, and less influenced by scattering determined by possible discontinuities in the subsurface.

A 3D crustal shear speed model of the Auckland volcanic field, New Zealand, from multi-component ambient noise tomography

We present a 3D shear-wave speed model of the crust in the intra-plate Auckland Volcanic Field (AVF), that provides new depth information on the basement geology and is a step towards determining lithospheric controls on volcanism of the AVF, within the Zealandia continent. Beneath the AVF, our model finds a mid-crustal contact of the upturning ophiolite-rich Dun Mountain-Maitai Terrane (DMT), and the western end of the accretionary Waipapa (Caples) Terranes which are dominated by terrigenous sedimentary rocks. Westward and at depth, the DMT presents as a high-speed, thick, and horizontal part of the Gondwana-ward verging syncline seen in the DMT elsewhere in New Zealand. Along the western side of the Junction Magnetic Anomaly (which indicates the presence of the DMT), a major vertical discontinuity probably marks the contact between the horizontal DMT and the Waipapa (Caples) Terranes. The Waipapa (Caples) Terrane contain a large wedge-shaped low-speed zone tapering eastward and extending NNW-SSE alongside the DMT. In addition, the model identifies a sharp vertical discontinuity at the eastern edge of the AVF that extends from the near-surface through the low-speed wedge. We interpret this discontinuity as a relatively young subvertical fault with a significant downthrow to the west. It appears to be a reactivation of a fault zone inherited from the break-up of Gondwana. The model also imaged the largest massif of crystalline ultramafic rocks in the DMT that causes the Takapuna Gravity Anomaly and the NE striking Port Waikato-Mangatangi cross-fault system. The strike directions of the faults mirror a rhombic-rectangular fault grid inherited from the breakaway of Zealandia from Gondwana.The shear speed model was produced using frequency-time analysis of multi-component Rayleigh wave Green's functions, reversible-jump Markov chain Monte Carlo inversion of group arrival times, and the inversion of dispersion relations for shear speed as a function of depth. Plain language summaryA 3D model of shear wave speeds in the crust provides new information on the geological structure beneath the Auckland Volcanic Field (AVF). The field is positioned above the contact between ocean floor rocks of the Dun Mountain-Maitai Terrane (DMT) in the west with fast wave speeds and a prominent wedge-shaped low-speed zone in the sandstone-dominated Waipapa/Caples Terrane in the east. The model images a fold structure of the DMT similar to that seen throughout New Zealand and at depth delineates the contact as a major sub-vertical discontinuity. A second vertical discontinuity just east of the AVF is very abrupt, is also oriented NNW-SSE and is probably caused by faulting. The model has also imaged the NE-oriented Port Waikato-Mangatangi cross-fault system in the south of the region and the large massif of coarse-grained intrusive rocks in the DMT that causes the previously recognised Takapuna gravity anomaly.

Joint inversion of temperature, vitrinite reflectance and fission tracks in apatite with examples from the eastern North Sea area

AbstractAs sediment accumulation indicates basin subsidence, erosion often is understood as tectonic uplift, but the amplitude and timing may be difficult to determine because the sedimentary record is missing. Quantification of erosion therefore requires indirect evidence, for example thermal indicators such as temperature, vitrinite reflectance and fission tracks in apatite. However, as always, the types and quality of data and the choice of models are important to the results. For example, considering only the thermal evolution of the sedimentary section discards the thermal time constant of the lithosphere and essentially ignores the temporal continuity of the thermal structure. Furthermore, the types and density of thermal indicators determine the solution space of deposition and erosion, the quantification of which calls for the use of inverse methods, which can only be successful when all models are mutually consistent. Here, we use integrated basin modelling and Markov Chain Monte Carlo inversion of four deep boreholes to show that the erosional pattern along the Sorgenfrei–Tornquist Zone (STZ) in the eastern North Sea is consistent with a tectonic model of tectonic inversion based on compression and relaxation of an elastic plate. Three wells in close proximity SW of the STZ have different data and exhibit characteristic differences in erosion estimates but are consistent with the formation of a thick chalk sequence, followed by minor Cenozoic erosion during relaxation inversion. The well on the inversion ridge requires ca. 1.7 km Jurassic‐Early Cretaceous sedimentation followed by Late Cretaceous–Palaeocene erosion during inversion. No well demands thick Cenozoic sedimentation followed by equivalent significant Neogene exhumation. When data are of high quality and models are consistent, the thermal indicator method yields significant results with important tectonic and geodynamic implications.

Metamorphic bedrock geometry of Santorini using HVSR information and geophysical modeling of ambient noise and active-source surface-wave data

Santorini is the most active volcano of the Southern Aegean Volcanic Arc, with volcanic material that has accumulated during the last 600+ Kyrs on top of the pre-volcanic Santorini island. The geometry of the pre-volcanic basement not only provides constraints on the volcanic history, but is also fundamental for seismic hazard assessment, as it partly controls site-effects on strong ground motions. We investigate the geometry of the metamorphic (Cycladic) basement of Santorini using information from both passive (noise) and active source data. We performed a large number (∼280) of single-station ambient noise measurements and collected additional data from 42 previous studies to compute the Horizontal-to-Vertical Spectral Ratio (HVSR) curves. The majority of the recovered HVSR curves show prominent spectral ratio peaks (locally larger than 7–8) indicating a large impedance contrast between pyroclastic volcanic strata and bedrock metamorphic formations. We also performed active (MASW) and passive (Noise Array) surface wave measurements at 4 selected sites to determine local dispersion curves. To find the 1D Vs structure from the surface down to the metamorphic bedrock we performed a joint Monte Carlo inversion of the dispersion curves and the fundamental frequency (f0HVSR) values from the HVSR curves. These models, as well as deep borehole information and HVSR data along the caldera cliff, where the thickness and stratigraphy of the volcanic formations above the bedrock can be directly observed, allowed us to determine a linear relationship between the thickness of the pyroclastic formations and their fundamental resonance period for central-southern Santorini. The resulting geometry of pre-volcanic island shows deep basins around the pre-Alpine bedrock outcrop (Profitis Ilias) and generally agrees with recent results from larger-scale tomographic and gravity data but with some significant local differences. The basins are 200+ meters thick in the Kamari and Perissa area (southeastern Santorini) and >400 m thick in the central (Fira-Imerovigli) area, while there is a thick (∼70-130 m) pyroclastic layer on top of a metamorphic basement “plateau” in the Megalochori area.

Autocorrelation R2 on Mars

AbstractA purpose of the Interior Exploration Using Seismic Investigations, Geodesy and Heat Transport (InSight) mission is to reveal the Martian interior structure with seismic data. In this work, ambient noise autocorrelation of the continuously recorded vertical‐component seismic signals has extracted the Rayleigh waves that propagate around Mars for one cycle, R2. The Mars orbiting surface waves are observed at a lag time of ∼6,000 s in the stacked autocorrelation series filtered between 0.005 and 0.01 Hz. Synthetic seismograms from a set of radially concentric velocity models were computed to find the best‐fitting one as the starting model for a Monte Carlo inversion. The starting model was randomly perturbed iteratively to increase the correlation coefficients and reduce the absolute time shifts between the synthetic and observed R2. An S‐wave low‐velocity layer in the inverted velocity model extends to ∼400 km depth, consistent with Marsquake observations, geophysical inversion, and high‐pressure experiments.

Along-strike variation in fault structural maturity and seismic moment deficits on the Yushu-Ganzi-Xianshuihe fault system revealed by strain accumulation and regional seismicity

The 900-km-long Yushu-Ganzi-Xianshuihe fault system that separates the Bayan Har and Qiangtang/Chuandian blocks is a primary tectonic structure in central eastern Tibet. Being one of the world's most seismically active fault systems, it hosted 12 M>7 earthquakes since the 18th century, and its future seismic risk has been the subject of extensive research over the past two decades. Previous studies have suggested that the western segment of the fault system differs largely from its eastern segment in strain accumulation and release, but the origin of such difference remains unknown. Here we use Sentinel-1 data to produce a complete map of present-day deformation for the entire fault system, in order to explore the along-strike variation in fault slip behaviour. Based on the InSAR deformation map, we solve for fault kinematic parameters including slip rate, locking depth and dip angle using a Bayesian Markov Chain Monte Carlo inversion approach. Slip rate on the Dangjiang-Yushu-Ganzi fault is consistently ∼5–6 mm/yr, only half of that on the Xianshuihe fault. This rate difference originates from nonuniform deformation mechanism within the Bayan Har block, with the western part of the block deforming more rigidly than the eastern part. Both historical and modern earthquake catalogues show that the western fault segment has only recorded M>7 earthquakes whereas the eastern segment has seen many small-to-moderate ruptures in addition to large earthquakes. These observations imply that the western fault segment may be more structurally mature than the eastern segment, which agrees with the long-term fault propagation as indicated by geochronological data and fault geometrical complexities. We calculate the seismic moment deficits along the fault system and find that several segments, e.g. the Dangjiang fault and the Manigange segment of the Ganzi fault, have accumulated a large seismic moment (equates to an Mw>7 earthquake) since the last major event.

Bayesian Markov Chain Monte Carlo inversion of surface-based transient electromagnetic data

Conventional linearized deterministic inversions of transient electromagnetic (TEM) data inherently simplify the non-uniqueness and ill-posed nature of the problem. While Monte-Carlo-type approaches allow for a comprehensive search of the solution space, gaining the ensemble of inferred solutions as comprehensive as possible may be limited utility in high-dimensional problems. To overcome these limitations, we utilize a Markov Chain Monte Carlo (MCMC) inversion approach for surface-based TEM data, which incorporates Bayesian concepts into Monte-Carlo-type global search strategies and can infer the posterior distribution of the models satisfying the observed data. The proposed methodology is first tested on synthetic data for a range of canonical earth models and then applied to a pertinent field dataset. The results are consistent with those obtained by standard linearized inversion approaches, but, as opposed to the latter, allow us to estimate the associated non-linear, non-Gaussian uncertainty.

Evaluation of the 3D Near-Surface Velocity Structure in an Urban Environment from Ambient Noise Array Tomography: The Case of the City of Thessaloniki (Northern Greece)

Abstract We examine the implementation of ambient noise array tomography in an urban environment to assess the 3D near-surface shear wave velocity (VS) structure at an intermediate spatial scale (∼1 km2, depth range 200–300 m). The application employs cross correlation traces of vertical component ambient noise recordings from a local network installed in Thessaloniki city (Northern Greece), allowing the determination of Rayleigh wave travel times for the frequency range of 1.5–14 Hz. The results confirm the presence of a complex subsurface with strong lateral variations in the geology, with travel times varying up to almost one order of magnitude. A surface wave travel time tomography approach was applied for each frequency to determine the spatial variability of the group velocity, involving the use of approximate Fresnel volumes, as well as damping and spatial smoothing constraints to stabilize the results. We also employed an interfrequency smoothing scheme to obtain smooth but data-compatible dispersion curves at the cost of inverting all travel time data simultaneously. Following the application of several quality cutoff criteria, we reconstructed local group slowness dispersion curves for a predefined tomographic grid in the study area. The final 3D velocity model was determined by a modified Monte Carlo inversion of these dispersion curves and the spatial integration of the obtained 1D VS profiles. Different model parameterizations were tested for the inversion to determine the optimal datafit. The final 3D velocity model is in a very good agreement with the local geology, previous larger scale studies, and other geophysical surveys, providing additional structural constraints (such as hidden fault identification) for the complex sedimentary deposits and bedrock formation in Thessaloniki, up to the depth of ∼250–300 m. The introduction of the aforementioned modifications to the ambient noise array tomography suggests that it can be efficiently adjusted and employed as a reliable tool for imaging the 3D seismic structure in urban environments with complex geology.

Bayesian tomography with prior-knowledge-based parametrization and surrogate modelling

SUMMARYWe present a Bayesian tomography framework operating with prior-knowledge-based parametrization that is accelerated by surrogate models. Standard high-fidelity forward solvers (e.g. finite-difference time-domain schemes) solve wave equations with natural spatial parametrizations based on fine discretization. Similar parametrizations, typically involving tens of thousand of variables, are usually employed to parametrize the subsurface in tomography applications. When the data do not allow to resolve details at such finely parametrized scales, it is often beneficial to instead rely on a prior-knowledge-based parametrization defined on a lower dimension domain (or manifold). Due to the increased identifiability in the reduced domain, the concomitant inversion is better constrained and generally faster. We illustrate the potential of a prior-knowledge-based approach by considering ground penetrating radar (GPR) traveltime tomography in a crosshole configuration with synthetic data. An effective parametrization of the input (i.e. the permittivity distributions determining the slowness field) and compression in the output (i.e. the traveltime gathers) spaces are achieved via data-driven principal component decomposition based on random realizations of the prior Gaussian-process model with a truncation determined by the performances of the standard solver on the full and reduced model domains. To accelerate the inversion process, we employ a high-fidelity polynomial chaos expansion (PCE) surrogate model. We investigate the impact of the size of the training set on the performance of the PCE and show that a few hundreds design data sets is sufficient to provide reliable Markov chain Monte Carlo inversion at a fraction of the cost associated with a standard approach involving a fine discretization and physics-based forward solvers. Appropriate uncertainty quantification is achieved by reintroducing the truncated higher order principle components in the original model space after inversion on the manifold and by adapting a likelihood function that accounts for the fact that the truncated higher order components are not completely located in the null space.

Machine learning‐accelerated gradient‐based Markov chain Monte Carlo inversion applied to electrical resistivity tomography

ABSTRACTExpensive forward model evaluations and the curse of dimensionality usually hinder applications of Markov chain Monte Carlo algorithms to geophysical inverse problems. Another challenge of these methods is related to the definition of an appropriate proposal distribution that simultaneously should be inexpensive to manipulate and a good approximation of the posterior density. Here we present a gradient‐based Markov chain Monte Carlo inversion algorithm that is applied to cast the electrical resistivity tomography into a probabilistic framework. The sampling is accelerated by exploiting the Hessian and gradient information of the negative log‐posterior to define a proposal that is a local, Gaussian approximation of the target posterior probability. On the one hand, the computing time to run the many forward evaluations needed for both the data likelihood evaluation and the Hessian and gradient computation is decreased by training a residual neural network to predict the forward mapping between the resistivity model and the apparent resistivity value. On the other hand, the curse of dimensionality issue and the computational effort related to the Hessian and gradient manipulation are decreased by compressing data and model spaces through a discrete cosine transform. A non‐parametric distribution is assumed as the prior probability density function. The method is first demonstrated on synthetic data and then applied to field measurements. The outcomes provided by the presented approach are also benchmarked against those obtained when a computationally expensive finite‐element code is employed for forward modelling , with the results of a gradient‐free Markov chain Monte Carlo inversion, and also compared with the predictions of a deterministic inversion. The implemented approach not only guarantees uncertainty assessments and model predictions comparable with those achieved by more standard inversion strategies, but also drastically decreases the computational cost of the probabilistic inversion, making it similar to that of a deterministic inversion.

Upscaling and downscaling Monte Carlo ensembles with generative models

SUMMARYMonte Carlo methods are widespread in geophysics and have proved to be powerful in non-linear inverse problems. However, they are associated with significant practical challenges, including long calculation times, large output ensembles of Earth models, and difficulties in the appraisal of the results. This paper addresses some of these challenges using generative models, a family of tools that have recently attracted much attention in the machine learning literature. Generative models can, in principle, learn a probability distribution from a set of given samples and also provide a means for rapid generation of new samples which follow that approximated distribution. These two features make them well suited for application to the outputs of Monte Carlo algorithms. In particular, training a generative model on the posterior distribution of a Bayesian inference problem provides two main possibilities. First, the number of parameters in the generative model is much smaller than the number of values stored in the ensemble, leading to large compression rates. Secondly, once trained, the generative model can be used to draw any number of samples, thereby eliminating the dependence on an often large and unwieldy ensemble. These advantages pave new pathways for the use of Monte Carlo ensembles, including improved storage and communication of the results, enhanced calculation of numerical integrals, and the potential for convergence assessment of the Monte Carlo procedure. Here, these concepts are initially demonstrated using a simple synthetic example that scales into higher dimensions. They are then applied to a large ensemble of shear wave velocity models of the core–mantle boundary, recently produced in a Monte Carlo study. These examples demonstrate the effectiveness of using generative models to approximate posterior ensembles, and indicate directions to address various challenges in Monte Carlo inversion.