Imaging the Earth's thermochemical structure is crucial for understanding its dynamics and evolution. Moreover, the increased demand for critical minerals and geothermal energy driven by the energy transition has intensified the need for reliable subsurface models. Multi-Observable Thermochemical Tomography (MTT) is a simulation-based, probabilistic inversion platform designed to harness the combined sensitivities of multiple geophysical data sets and thermodynamic modeling. It produces internally consistent estimates of the Earth's interior as probability distributions, offering a powerful means for uncertainty quantification. Here, we present an updated MTT formalism and assess its benefits and limitations to image the thermochemical structure of the lithosphere-asthenosphere system. Individual and combined sensitivities of different observables to parameters of interest (e.g., temperature, composition, crustal architecture) are explored using challenging synthetic models. Our findings demonstrate that a judicious combination of observables can retrieve complex thermochemical structures relevant to greenfields exploration. We then apply MTT to study two cratonic regions of geological and economic significance. In the Superior Craton, we jointly invert receiver functions, gravity anomalies, gravity gradients, geoid anomalies, Rayleigh-wave dispersion curves, absolute elevation and surface heat flow. In the North Australian Craton, we incorporate new data from the AusArray and add teleseismic P- and S-phase travel times to the data sets. The imaged lithospheric architectures provide new insights into the tectonic evolution of these two regions and the physical meaning of geophysical signatures. Additionally, these models offer unique proxies to guide exploration efforts for clean energy and critical minerals and serve as reference models for future high-resolution studies.

The solid earth structure beneath Greenland, meaning the rocky part of Earth from the ice-bed interface to depth, has gained increased interest in recent years as it provides a critical boundary condition for the dynamic evolution of the Greenland ice sheet (GrIS), one of the largest sources of sea-level rise contributions since the early 2000s. However, no consensus has been reached regarding the key internal or surface earth properties influencing this boundary condition and thus GrIS behaviour. One important surface property is the subglacial heat flow, which affects sliding conditions of the ice sheet including the onset of major ice streams and is related to subglacial geology. Lithospheric architecture and mantle viscosity structure are internal properties that influence ice sheet evolution through changes in the height and slope of the ice-bed interface caused by glacial isostatic adjustment. Because there is no general agreement regarding crustal and lithospheric structures, some glaciological studies use an ensemble of solid earth models to incorporate uncertainties into their GrIS predictions, but it is unclear how these variations ultimately affect estimates of future sea-level rise. Here we describe the main solid earth properties that are important for GrIS evolution (heat flow, temperature, viscosity), from the base of the ice sheet to the upper mantle, and we provide some perspectives on how future collaborative efforts and integrated studies could lead to better agreement regarding these key characteristics.

The North Atlantic region is a complex geodynamic setting that comprises multiple continental blocks, sedimentary basins, mid-ocean ridge systems and prominent hotspots. Recent geophysical surveys of the near-surface have enhanced our understanding of crustal elements and the shallow lithosphere. However, our knowledge of the deep lithospheric structure and the physical state and dynamics of the upper mantle is still limited. Here, we exploit the combined sensitivity of surface-wave data, geoid anomalies, absolute topography and surface heat flow to obtain full thermochemical models of the region from the surface down to 350 km. We jointly invert these data sets using a simulation-based, multi-observable probabilistic framework. We validate our results with independent thermobarometric and chemical information from mantle xenoliths and test the effects of using different seismic models on the inversion results. Our model reveals an intricate sublithospheric flow system, driven by the interaction of deep upwellings with the highly irregular lithospheric structure. We corroborate that the main thermal anomaly in the sublithospheric mantle shows a tilted geometry, moving toward Greenland with depth. We reveal that this large-scale anomaly transition into a more complex pattern once it reaches depths of similar to 150 km beneath the North Atlantic. Small-scale downwellings originate from the margins of continental domains, resulting in a complex circulation pattern that limits the radial spread of the deep upwellings and preferentially focuses them within regions of thin lithosphere along a N-S direction. Distinct compositional anomalies in the Greenland lithosphere delineate the North Atlantic Craton, the Nagssugtoqidian mobile belt, and the covered remnants of the Disko Craton. In continental Europe, the East European Craton shows clear indications of depletion in incompatible elements, with the Kola-Karelian cratonic region showing the highest levels of depletion. Our model serves as a base to make interpretations on the enigmatic paleotectonic history of the North-Atlantic region.

The receiver function technique is widely used to image crustal structure using P-to-S converted phases at the Moho discontinuity. However, the presence of sedimentary layer generates additional P-to-S conversions and reverberations, which can overprint the Moho phases and pose problems in imaging crustal structure with standard receiver function techniques. We introduce a robust two-step method that uses H-kappa stacking to determine average thickness and Vp/Vs of the sedimentary layer, followed by waveform-fitting of the observed receiver function to constrain the average crustal thickness and sub-sediment Vp/Vs. We tested the method using both synthetic data and real-data from stations located on sedimentary layers in the Netherlands and USA. We show that the new method outperforms other common approaches in retrieving accurate Moho depth and sub-sediment Vp/Vs estimates, even in cases where the Moho phases are completely overprinted by large-amplitude phases related to sedimentary layers. We propose a two-step method for modeling crustal structure using receiver functions from stations overlying sedimentary layers. In the first step, the sediment layer thickness and velocity ratio (Vp/Vs) are derived from depth versus Vp/Vs stacking of receiver functions. In the second step, the crust-mantle interface (Moho) depth and sub-sediment Vp/Vs are derived from the waveform-fitting grid search of the observed receiver functions with synthetic ones. We tested this method using synthetic receiver functions and real-data from stations in the Netherlands and USA. The new method works well in retrieving Moho depth and sub-sediment Vp/Vs, even in cases where the Moho signals are completely obscured by signals from the sediment layer. New approach for modeling receiver functions for stations overlying a sedimentary layer The approach is effective even when Moho phases are completely overprinted by positive or negative amplitudes from the sedimentary layer Synthetic tests and real-data applications show the advantages of the developed approach over common methods for addressing sediment impact

Abstract The formation of large igneous provinces (LIPs) has been widely believed to be linked to mantle plume activity. However, how the plume modifies the overlying lithosphere, particularly its compositional structure, remains uncertain. Here, we characterize the deep thermochemical structure beneath the Emeishan LIP (ELIP), which is a well‐known Permian plume‐related LIP in China, by taking a multi‐observable probabilistic inversion. Our results find a clear correlation between the lithospheric composition with the ELIP's concentric zones. We infer that the fertile feature of the lithospheric mantle in the ELIP's inner zone was caused by the plume‐derived fertile magmas which infiltrated into and chemically refertilized the ambient depleted lithosphere. This plume‐modified lithospheric compositional structure is likely to be preserved after the plume event, while the present lithospheric thermal structure has been mainly influenced by the subsequent thermal‐tectonic activity. Our results improve our understanding of the physicochemical interactions between the lithosphere and ancient plume.

The Archean Superior craton was formed by the assemblage of continental and oceanic terranes at ∼2.6 Ga. The craton is surrounded by multiple Proterozoic mobile belts, including the Paleoproterozoic Trans-Hudson Orogen which brought together the Superior and Rae/Hearne cratons at ∼1.9-1.8 Ga. Despite numerous studies on Precambrian lithospheric formation and evolution, the deep thermochemical structure of the Superior craton and its surroundings remains poorly understood. Here we investigate the upper mantle beneath the region from the surface to 400 km depth by jointly inverting Rayleigh wave phase velocity dispersion data, elevation, geoid height and surface heat flow, using a probabilistic inversion to obtain a (pseudo-)3D model of composition, density and temperature. The lithospheric structure is dominated by thick cratonic roots (>300 km) beneath the eastern and western arms of the Superior craton, with a chemically depleted signature (Mg# > 92.5), consistent with independent results from mantle xenoliths. Beneath the surrounding Proterozoic and Phanerozoic orogens, the Mid-continent Rift and Hudson Strait, we observe a relatively thinner lithosphere and more fertile composition, indicating that these regions have undergone lithospheric modification and erosion. Our model supports the hypothesis that the core of the Superior craton is well-preserved and has evaded lithospheric destruction and refertilization. We propose three factors playing a critical role in the craton’s stability: (i) the presence of a mid-lithospheric discontinuity, (ii) the correct isopycnic conditions to sustain a strength contrast between the craton and the surrounding mantle, and (iii) the presence of weaker mobile belts around the craton.

The thermochemical structure of the lithosphere exerts control on melting mechanisms in the mantle as well as the location of volcanism and ore deposits. Imaging the complex interactions between the lithosphere and asthenospheric mantle requires the joint inversion of multiple data sets and their uncertainties.In particular, the combination of seismic velocity and electrical conductivity with data proxies for bulk composition and elusive minor phases is a crucial step towards fully understanding large-scale lithospheric structure and melting.We apply a novel probabilistic approach for joint inversions of 3D magnetotelluric and seismic data to image the lithosphere beneath southeast Australia. Results show a highly heterogeneous lithospheric structure with deep conductivity anomalies that correlate with the location of Cenozoic volcanism. In regions where the conductivities have been at odds with sub-lithospheric temperatures and seismic velocities, we observe that the joint inversion provides conductivity values consistent with other observations. The results reveal a strong relationship between metasomatized regions in the mantle and i) the limits of geological provinces in the crust, which elucidates the subduction-accretion process in the region; ii) distribution of leucitite and basaltic magmatism; iii) independent geochemical data, and iv) a series of lithospheric steps which constitute areas prone to generating small-scale instabilities in the asthenosphere. This scenario suggests that shear-driven upwelling and edge-driven convection are the dominant melting mechanisms in eastern Australia rather than mantle plume activity, as conventionally conceived. Our study offers an integrated lithospheric model for southeastern Australia and provides insights into the feedback mechanism driving surface processes.

We present a particular derivation of the face-centred finite volume (FCFV) method and study its performance in non-linear, coupled transport problems commonly encountered in geoscientific and geotechnical applications. The FCFV method is derived from the hybridisable discontinuous Galerkin formulation, using a constant degree of approximation for the discretization of the unknowns defined on the mesh faces (edges in two dimensions). The piecewise constant degrees of freedom are determined in a global problem over the mesh skeleton. Then, the solution and its gradient are recovered at the cells centroid in a set of element-by-element independent postprocesses, both exhibiting linear convergence. The formulation of the transient advection-diffusion-reaction equation is presented in detail and the numerical analysis under challenging advective/diffusive regimes is studied. Finally, we use several numerical examples to illustrate the advantages and limitations of the FCFV method to solve problems of geoscientific and geotechnical relevance governed by the non-linear coupling between advection-diffusion-reactive transport and Stokes flow. Our results show that the FCFV method is an attractive and highly competitive alternative to other commonly used methods.

The thermochemical structure of the subcontinental mantle holds information on its origin and evolution that can inform energy and mineral exploration strategies, natural hazard mitigation and evolutionary models of Earth. However, imaging the fine-scale thermochemical structure of continental lithosphere remains a major challenge. Here we combine multiple land and satellite datasets via thermodynamically constrained inversions to obtain a high-resolution thermochemical model of central and southern Africa. Results reveal diverse structures and compositions for cratons, indicating distinct evolutions and responses to geodynamic processes. While much of the Kaapvaal lithosphere retained its cratonic features, the western Angolan–Kasai Shield and the Rehoboth Block have lost their cratonic keels. The lithosphere of the Congo Craton has been affected by metasomatism, increasing its density and inducing its conspicuous low-topography, geoid and magnetic anomalies. Our results reconcile mantle structure with the causes and location of volcanism within and around the Tanzanian Craton, whereas the absence of volcanism towards the north is due to local asthenospheric downwellings, not to a previously proposed lithospheric root connecting with the Congo Craton. Our study offers improved integration of mantle structure, magmatism and the evolution and destruction of cratonic lithosphere, and lays the groundwork for future lithospheric evolutionary models and exploration frameworks for Earth and other terrestrial planets. Cratons in central and southern Africa exhibit diverse structures, compositions and responses to geodynamic settings, according to a high-resolution thermochemical regional model constructed from land- and satellite-based geophysical observations.

<p>Seismic data from several long-running broadband sensors around the globe will be used to investigate the statistical significance and geologic interpretation of negative velocity discontinuities in the upper-most mantle. Several previous studies have identified negative polarity arrivals in S-wave receiver function data which are variously interpreted as lithosphere-asthenosphere and/or mid-lithosphere boundaries.</p> <p>One-dimensional joint-inversion is applied using the LitMod framework, which is a Bayesian statistical method driven by a Markov Chain Mote Carlo algorithm.</p> <p>LitMod uses a thermodynamically consistent physical model of the mantle and thus provides important constraints for the interpretation of the receiver function results.</p> <p>Joint inversions combine Rayleigh wave phase velocity measurements, both P and S-wave receiver functions, absolute elevation, and geoid height.</p> <p>Particular attention is given to the calculation and inversion of S-wave receiver function data, which represents a new addition to the LitMod framework. </p>

Geoid anomalies offer crucial information on the internal density structure of the Earth, and thus, on its constitution and dynamic state. In order to interpret geoid undulations in terms of depth, magnitude and lateral extension of density anomalies in the lithosphere and upper mantle, the effects of lower mantle density anomalies need to be removed from the full geoid (thus obtaining a residual signal known as the 'upper mantle geoid'). However, how to achieve this seemingly simple filtering exercise has eluded consensus for decades in the solid Earth community. While there is wide agreement regarding the causative masses of degrees > 10 in spherical harmonic expansions of the upper mantle geoid, those contributing to degrees < 7-8 remain ambiguous. Here we use spherical harmonic analysis and recent tomography and density models from joint seismic-geodynamic inversions to derive a representative upper mantle geoid, including the contributions from low harmonic degrees. We show that the upper mantle geoid contains important contributions from degrees 5 and 6 and interpret the causative masses as arising from the coupling between the long-wavelength lithospheric structure and the sublithospheric upper mantle convection pattern, including subducted slabs. Importantly, the contributions from degrees 3 < l < 8 do not show a simple power-law behaviour (e.g. Kaula's rule), which precludes the use of standard filtering techniques in the spectral domain. Our model of the upper mantle geoid will be useful in a wide range of geodynamic and geophysical applications, including the study of i) the thermochemical structure of the lithosphere, ii) dynamic topography and mantle viscosity, iii) the nature of the mechanical coupling of the lithosphere-asthenosphere system and iv) the global state of stress within the lithosphere and its associated hazards.

The recycling of oceanic lithosphere into the deep mantle at subduction zones is one of the most fundamental geodynamic processes on Earth. During the closure of an ocean, ancient oceanic slabs are thought to be consumed entirely in subduction zones due to their negative buoyancy. Yet, it is recently suggested that small pieces of oceanic slabs could be trapped along paleo-subduction zones. What remains far more enigmatic is whether significant portions of paleo-oceanic lithosphere could eventually avoid the fate of subduction and be accreted to continental lithosphere, thus contributing to continental growth through time. We present seismic evidence for a preserved paleo-oceanic lithosphere beneath the Junggar region in northwestern China. We show that unsubducted oceanic lithosphere in the West Junggar has been preserved beneath the Junggar Basin, becoming a piece of the Eurasian continent. This scenario is likely to have occurred in other continents throughout Earth’s history, providing an additional and commonly underestimated contribution to the growth of continental lithosphere.

The Middle and Lower Reaches of the Yangtze River metallogenic belt (MLYMB) is one of the most important Fe‐Cu polymetallic belts in China. However, the mechanism and deep geodynamical process for the formation of this belt are still controversial. Here, we obtain the crustal and the uppermost mantle structures using ambient noise data from a dense seismic profile. A low velocity zone is revealed beneath the Moho of MLYMB, interpreted as the source of the deep mineralization materials. In addition, a low velocity layer (LVL) and a high velocity layer (HVL) are observed in the crust of the southern segment of the profile. The LVL is interpreted as a tectonic detachment layer between the upper and the lower crust, and the HVL is interpreted as the aggregation zone for mineralizing melts or crystallized magma chambers. Based on the observed velocity features, we propose a three‐stage model for the formation of ore deposits in MLYMB. Our model suggests that an upwelling of asthenosphere triggered by the delamination of a previously thickened lithosphere leads to the partial melting of upper mantle rocks, which eventually ponders under the Moho. The magma then infiltrates through the ductile lower crust and reaches a depth of ∼7–13 km, forming a minerals‐enriched magma chamber. Minerals‐rich hot fluids originating from the magma chamber continue to move upward along the pre‐existent faults and the minerals finally precipitate in dense veinlets when reaching shallow depths, forming the ore deposits in and around the MLYMB.

Geoid anomalies offer crucial information on the internal density structure of the Earth, and thus, on its constitution and dynamic state. In order to interpret geoid undulations in terms of depth, magnitude and lateral extension of density anomalies in the lithosphere and upper mantle, the effects of lower mantle density anomalies need to be removed from the full geoid (thus obtaining a residual signal known as the 'upper mantle geoid'). However, how to achieve this seemingly simple filtering exercise has eluded consensus for decades in the solid Earth community. While there is wide agreement regarding the causative masses of degrees > 10 in spherical harmonic expansions of the upper mantle geoid, those contributing to degrees < 7-8 remain ambiguous. Here we use spherical harmonic analysis and recent tomography and density models from joint seismic-geodynamic inversions to derive a representative upper mantle geoid, including the contributions from low harmonic degrees. We show that the upper mantle geoid contains important contributions from degrees 5 and 6 and interpret the causative masses as arising from the coupling between the long-wavelength lithospheric structure and the sublithospheric upper mantle convection pattern, including subducted slabs. Importantly, the contributions from degrees 3 < l < 8 do not show a simple power-law behaviour (e.g. Kaula's rule), which precludes the use of standard filtering techniques in the spectral domain. Our model of the upper mantle geoid will be useful in a wide range of geodynamic and geophysical applications, including the study of i) the thermochemical structure of the lithosphere, ii) dynamic topography and mantle viscosity, iii) the nature of the mechanical coupling of the lithosphere-asthenosphere system and iv) the global state of stress within the lithosphere and its associated hazards.

&lt;p&gt;The present-day thermochemical structure of the subcontinental mantle holds key information on its origin and evolution and informs exploration strategies, natural hazard management and evolutionary model of the Earth system. As such, unravelling the nature of the continental lithosphere, its modification through time and its interactions with the sublithospheric mantle and the atmosphere/hydrosphere constitute some of the main goals of modern geoscience. Despite its fundamental importance, imaging the fine-scale thermochemical structure of the lithosphere using indirect (remote) data is plagued with difficulties, which has traditionally left the analysis of xenoliths and xenocrysts as the only reliable approach.&lt;/p&gt; &lt;p&gt;In recent years, however, &amp;#8216;simulation-based&amp;#8217; inverse methods that integrate multiple geophysical and geochemical datasets within an internally- and thermodynamically-consistent platform have opened new and promising ways to address this &amp;#8216;grand challenge&amp;#8217;. In this presentation, I will discuss i) some recent progress, case studies and future directions on the mapping of the thermochemical structure of the continental lithosphere, and ii) their predictive power for the energy and critical minerals sectors and possible implications for planetary exploration in general.&lt;/p&gt;

The formation and preservation of compositional heterogeneities inside the Earth affect mantle convection patterns globally and control the long-term evolution of geochemical reservoirs. However, the distribution, nature, and size of reservoirs in the Earth's mantle are poorly constrained. Here, we invert measurements of travel times and amplitudes of seismic waves interacting with mineralogical phase transitions at 400-700-km depth to obtain global probabilistic maps of temperature and bulk composition. We find large basalt-rich pools (up to 60% basalt fraction) surrounding the Pacific Ocean, which we relate to the segregation of oceanic crust from slabs that have been subducted since the Mesozoic. Segregation of oceanic crust from initially cold and stiff slabs may be facilitated by the presence of a weak hydrated layer in the slab or by weakening upon mineralogical transition due to grain-size reduction.

Joint probabilistic inversions of magnetotelluric (MT) and seismic data has great potential for imaging the thermochemical structure of the lithosphere as well as mapping fluid/melt pathways and regions of mantle metasomatism. In this contribution we present a novel probabilistic (Bayesian) joint inversion scheme for 3D MT and surface-wave dispersion data particularly designed for large-scale lithospheric studies. The approach makes use of a recently developed strategy for fast solutions of the 3D MT forward problem (Manassero et al.,2020) and combines it with adaptive Markov chain Monte Carlo (MCMC) algorithms and parallel-in-parallel strategies to achieve extremely efficient simulations. To demonstrate the feasibility, benefits and performance of our joint inversion method to image the temperature and conductivity structures of the lithosphere, we apply it to two numerical examples of increasing complexity. The inversion approach presented here is timely and will be useful in the joint analysis of MT and surface wave data that are being collected in many parts of the world. This approach also opens up new avenues for the study of translithospheric and transcrustal magmatic systems, the detection of metasomatised mantle and the incorporation of MT into multi-observable inversions for the physical state of the Earth's interior.

Methodology, supplemental Figures S1–S11, and Tables S1 and S2.

In the context of whole-lithosphere structure, the joint inversion of magnetotelluric (MT) with seismic data is particularly interesting as they provide complementary information on the thermal structure, fluid pathways and water content. Both data sets can put tight constrains on the first-order thermal structure and mineralogical structure of the lithosphere, but only MT is strongly sensitive to anomalous features such as hydrogen content, minor conductive phases and/or small volumes of fluid or melt. This makes joint inversions of MT with other observables a powerful means to detect fluid pathways in the lithosphere including the locus of partial melting, ore deposits and hydrated (or metasomatized) lithologies. This unique potential of joint inversions of MT with other datasets has given impetus to the acquisition of collocated MT and seismic data over large regions. Concrete examples are the US Array, Sinoprobe in China, and the AusLAMP/AusArray in Australia. These multi-disciplinary programs are providing high-quality seismic and MT data with unprecedented resolution and coverage, allowing the pursuit of large-scale 3D joint inversions to image the structure, dynamics and evolution of the whole lithosphere and upper mantle. Within probabilistic approaches the solution to the inverse problem is given by the so-called posterior probability density function which provides complete information about the unknown parameters and their uncertainties conditioned on the data and modelling assumptions. Joint probabilistic inversions of MT and seismic data have been successfully implemented in the context of 1D MT data only. For the cases of 2D and 3D MT data, however, the large computational cost of the MT forward problem has been the main impediment for pursuing probabilistic inversions, as the number of forward solutions required are typically on the order of 105 – 107. To overcome this limitation, we have recently presented a novel strategy [2,3], called RB+MCMC, that computes 3D MT surrogate models and uses complementary parameterizations to couple different data sets. This strategy reduces the computational cost of the 3D MT forward solver and allow us to perform full joint probabilistic inversions of MT and other datasets for the 3D imaging of deep thermochemical anomalies. In this contribution, we first illustrate the benefits and general capabilities of our method for 3D joint probabilistic inversions of MT with other datasets using whole-lithosphere synthetic models. Last, as part of the Exploring for the Future program, we present results of the first joint probabilistic inversion of 3D MT in southeast Australia using the AusLAMP data and a seismic velocity model derived from teleseismic tomography [4]. These results demonstrate the capabilities of our conceptual and numerical framework for 3D joint probabilistic inversions of MT with other geophysical data sets and open up exciting opportunities for elucidating the Earth’s interior in other regions. References [1] Afonso, J.C. et al., (2016), Journal of Geophysical Reseach, 121, doi:10.1002/2016JB013049 [2] Manassero, M. C., et al., (2020), Geophysical Journal International, 223(3), doi: 10.1093/gji/ggaa415 [3] Manassero, M. C., et al., (2021), doi: 10.1029/2021JB021962 [4] Rawlinson, N., et al., (2016), Tectonophysics, doi: 10.1016/j.tecto.2015.11.034