To reduce uncertainties in reconstructed images, geologic information must be introduced in a numerically robust and stable way during the geophysical data inversion procedure. In the context of potential (gravity) data inversion, it is important to bound the physical properties by providing probabilistic information on the number of lithologies and ranges of values of possibly existing related rock properties (densities). For this purpose, we have introduced a generalization of bounding constraints for geophysical inversion based on the alternating direction method of multipliers. The flexibility of the proposed technique enables us to take into account petrophysical information as well as probabilistic geologic modeling, when it is available. The algorithm introduces a priori knowledge in terms of physically acceptable bounds of model parameters based on the nature of the modeled lithofacies in the region under study. Instead of introducing only one interval of geologically acceptable values for each parameter representing a set of rock properties, we define sets of disjoint intervals using the available geologic information. Different sets of intervals are tested, such as quasidiscrete (or narrow) intervals as well as wider intervals provided by geologic information obtained from probabilistic geologic modeling. Narrower intervals can be used as soft constraints encouraging quasidiscrete inversions. The algorithm is first applied to a synthetic 2D case for proof-of-concept validation and then to the 3D inversion of gravity data collected in the Yerrida Basin (Western Australia). Numerical convergence tests show the robustness and stability of the bound constraints that we apply, which is not always trivial for constrained inversions. This technique can be a more reliable uncertainty reduction method as well as an alternative to other petrophysically or geologically constrained inversions based on the more classic “clustering” or Gaussian-mixture approaches.

SUMMARY: Gravity inversion methods are able to recover density distributions in the Earth but they require strong constraints based on various prior information. In the present study, our aim is to invert gravity anomalies using existing geological models and density information in orogenic areas like the extensively studied Pyrenees region. Geophysical exploration for resource assessment using gravity and seismic data already accessible in this area has been performed to identify potential fluid/gas reservoirs of economic significance. For a given pre-existing model defined by different geological units and summarizing our knowledge of the area, we aim at retrieving the density distribution within each unit. For this, we use an Alternate Direction Multiplier Method to perform gravity data inversion by constrained by Interval Bound Constraints (IBC-ADMM) defined as bounding distribution intervals of possible density values. To estimate these bounding intervals, we first use a prior density model geologically compatible and obtained using geological information, gravity modelling, seismological data and seismic models. In a second step, we apply a Taguchi statistical analysis on representative density variations inside each unit of the prior density model to estimate their impact on data residuals, and reduce the parameter space prior to inversion. We perform the gravity data inversion constrained by these Taguchi-derived density intervals and we estimate model uncertainties for these density intervals. We apply the technique to the entire Pyrenees range at 2 km resolution. The sensitivity analysis shows that the strongest variations on gravity data are obtained essentially in the Axial Zone, in the middle and lower crust, and in the Upper mantle. Inverted densities are compared to the prior model. They are increased in the central part of the Pyrenees Axial Zone while, in the deeper structures (from middle crust down to upper mantle), they are decreased in the Southern (Spanish) Iberian upper mantle and increased in the Northern (French) Eurasian Upper mantle. A possible shortening of the Iberian lower crust–upper mantle system can be assumed below the western part of the Axial Zone.

In previous studies, the Auxiliary Differential Equation (ADE) method has been applied to 2-dimensional (2D) seismic-wave propagation modelling in viscoelastic media. This method is based on the separation of the wave propagation equations derived from the constitutive law defining the stress-strain relation. We make here a 3D extension of a finite-difference (FD) scheme to solve a system of separated equations consisting in the stress-strain rheological relation, the strain-velocity and the velocity-stress equations. The current 3D FD scheme consists in the discretization of the second order formulation of a nonlinear viscoelastic wave equation with a time actualization of the velocity and displacement fields. Compared to the usual memory variable formalism, the ADE method allows flexible implementation of complex expressions of the desired rheological model such as attenuation/viscoelastic models or even nonlinear behaviours, with physical parameters that can be provided from dispersion analysis. The method can also be associated with optimized perfectly matched layers-based boundary conditions that can be seen as additional attenuation (viscoelastic) terms. We present the results obtained for a nonlinear viscoelastic model made of a Zener viscoelastic body associated with a nonlinear quadratic strain term. Such nonlinearity is relevant to define unconsolidated granular model behaviour. Thanks to a simple model, but without loss of generality, we demonstrate the accuracy of the proposed numerical approach.

We image the internal density structure of the Pyrenees by inverting gravity data using an a priori density model derived by scaling a V_p model obtained by full waveform inversion of teleseismic P-waves. Gravity anomalies are computed via a 3-D high-order finite-element integration in the same high-order spectral-element grid as the one used to solve the wave equation and thus to obtain the velocity model. The curvature of the Earth and surface topography are taken into account in order to obtain a density model as accurate as possible. The method is validated through comparisons with exact semi-analytical solutions. We show that the spectral-element method drastically accelerates the computations when compared to other more classical methods. Different scaling relations between compressional velocity and density are tested, and the Nafe–Drake relation is the one that leads to the best agreement between computed and observed gravity anomalies. Gravity data inversion is then performed and the results allow us to put more constraints on the density structure of the shallow crust and on the deep architecture of the mountain range.

Greenstone belts in the West African Craton (WAC) are separated by several generations of granitoids intruded at ca. 2.18–1.98 Ga. Simultaneous folding and exhumation play an important role in the formation of greenstone-granitoid belts. However, the overall tectonic regime and origin of granitoids remain controversial. In this study, we present the estimates of the mantle potential temperature (Tp) for the WAC, which yields values of about 1500–1600 °C, pressure estimates of initial and final melting yield values of about 3.7–5.2 GPa and 1–1.3 GPa, respectively. Subsequently, 2D thermo-mechanical models have been constructed to explore the width of volcano-sedimentary basin on spatial-temporal evolution of diapirs that emplaced in the lower-middle crust during compression. The models show that the width of the volcano-sediment layer plays an important role in the formation mechanisms of greenstone-granitoid belts. The lower crust beneath sedimentary sequences is deformed into a buckle fold during the first compressional stage, through which relief uplifts slowly. Subsequently, the buckle fold is further deformed into several individual folds. Diapirs made of lower crust rocks ascend and emplace in the middle-upper crust resulting from instability. Benefitting from the mantle temperature, the pressure estimates and the numerical modelling results, a new geodynamic model was constructed. This model indicates that a series of sheet-like granitoids possibly derived from either subducted mélanges, lower crust and/or mantle melting that are accumulated at depths of the subcontinental mantle would channel along diapirs before feeding the upper crust. When the granitoids arrive at the solidified lids of the diapirs, they would favour migrating horizontally and intrude into the upper crust through weakening zones between the diapirs. Our geodynamic model also suggests an asymmetry of structures between the upper and middle-lower crust, with the dome-like granitoids overlying high-grade sedimentary synforms and high-grade diapirs underlying low-grade greenstone belts.