Next Article in Journal
Development and Calibration of a 3D Micromodel for Evaluation of Masonry Infilled RC Frame Structural Vulnerability to Earthquakes
Previous Article in Journal
Defining the Geotourism Potential of the CILENTO, Vallo di Diano and Alburni UNESCO Global Geopark (Southern Italy)
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Seismic Constrained Gravity Inversion: A Reliable Tool to Improve Geophysical Models Away from Seismic Information

by
Daniele Sampietro
*,† and
Martina Capponi
Geomatics Research & Development s.r.l., c/o ComoNExT, Via Cavour 2, 22074 Lomazzo, CO, Italy
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Geosciences 2021, 11(11), 467; https://doi.org/10.3390/geosciences11110467
Submission received: 15 October 2021 / Revised: 9 November 2021 / Accepted: 10 November 2021 / Published: 12 November 2021
(This article belongs to the Section Geophysics)

Abstract

:
The exploitation of gravity fields in order to retrieve information about subsurface geological structures is sometimes considered a second rank method, in favour of other geophysical methods, such as seismic, able to provide a high resolution detailed picture of the main geological horizons. Within the current work we prove, through a realistic synthetic case study, that the gravity field, thanks to the availability of freely of charge high resolution global models and to the improvements in the gravity inversion methods, can represent a valid and cheap tool to complete and enhance geophysical modelling of the Earth’s crust. Three tests were carried out: In the first one a simple two-layer problem was considered, while in tests two and three we considered two more realistic scenarios in which the availability on the study area of constraints derived from 3D or 2D seismic surveys were simulated. In all the considered test cases, in which we try to simulate real-life scenarios, the gravity field, inverted by means of an advanced Bayesian technique, was able to obtain a final solution closer to the (simulated) real model than the assumed a priori information, typically halving the uncertainties in the geometries of the main geological horizons with respect to the initial model.

Graphical Abstract

1. Introduction

The study of the Earth’s crust by means of indirect geophysical methods, the so called potential field methods, is crucial for a great variety of scientific and industrial activities. In particular, focusing on the gravity field, nowadays there is a broad availability of worldwide datasets and models thanks to the observations of several dedicated satellite missions, to the exploitation of satellite altimetry (offshore) and to the advancements in airborne and shipborne gravimetry. These datasets and models can be profitably used to determine the Earth crust major characteristics from a global to a very local scale (see, for instance, [1,2,3]). The recently developed and freely available Global Gravity field Models (GGMs), such as the XGM2019e model [4] with a spatial resolution of about 3600 m and an accuracy of the order of 3 mGal (for the sake of readability, within the current work we choose to express gravity related quantities in mGal instead of in m/s 2 , i.e., in the International System of Units, recalling that 1 mGal = 1 5 m/s 2 ), are paving the way to geophysical applications based on the processing and interpretation of gravity field data both for scientific as well as industrial purposes without the need of further gravity surveys. As for the industrial applications, which are mainly related to natural resources exploration and in particular hydrocarbons, the gravitational field is usually exploited in two scenario: It can be used in the preliminary stage of the exploration activities or in complex areas where seismic surveys fail. In the former, gravity is exploited to perform a first study at large scale (e.g., about 5 · 10 5 m × 5 · 10 5 m or even more) in order to optimise and further focus other, and usually more expensive, geophysical surveys. In the latter scenario, where the presence of peculiar geological formations (e.g., salt diapirism) lead the seismic survey to fail, gravity can be profitably used to map the deepest layers. The main characteristics of both these scenarios are summarised in the following:
Scenario 1, advanced modelling of deepest layers:
  • a local (smaller than 1 · 10 5 m × 1 · 10 5 m) area to be investigated;
  • a good geological knowledge of the area available from previous studies;
  • several seismic profiles or even 3D seismic available, properly constraining the shallowest geological structures and units;
  • well logs to constrain density distribution in the shallowest layers;
  • Moho and basement depths known with poor spatial resolution and accuracy.
Scenario 2, exploration in frontier areas:
  • a large (greater than 1 · 10 5 m × 1 · 10 5 m) area to be investigated;
  • a general geological knowledge of the area available from literature;
  • few seismic profiles available, giving sparse information on the shallowest geological structures and units;
  • few, or even no, constraints for the density distribution;
  • Moho and basement depths known only from global models (such as the CRUST1.0 model [5]) but with poor spatial resolution and accuracy.
In this context, the latest GGMs eventually combined with more detailed gravity data collected at a ground level, can be exploited to improve the geophysical modelling of the study area without the need for further dedicated gravity surveys. This is especially true in offshore regions where the contribution of satellite altimetry has significantly improved GGMs and since accurate gravity is usually collected together with seismic acquisitions. In this respect, the bottleneck for both the above scenarios is now represented by the development of ad-hoc techniques to fully exploit the gravity observations together with the available a priori knowledge. It is in fact well known that gravity inversion, in order to be effective, should rely on a set of external constraints (such as well logs or information from seismic profiles) to face the inherent problem related to the non-uniqueness of its solution (see e.g., [6]). In the current work, after recalling the main idea behind the Bayesian inversion algorithm presented in [2,7,8], we will assess the accuracy and reliability of the inversion on a realistic synthetic case study similar to those usually faced in industrial resource exploration activities trying to simulate the above Scenarios 1 and 2. In doing this, we are also motivated by the dual objective of generating a complete realistic case study that can be in general used to assess gravity inversion solutions and of finalising a methodology to improve geophysical modelling based on seismic data by means of gravity inversion.
In the next section (Section 2), the main ideas behind the Bayesian inversion algorithm used in the current study are recalled, while in Section 3 the developed synthetic realistic case study is described. The different tests performed are reported in Section 4, and finally, in Section 5, some conclusions are drawn.

2. Methods

The inversion of the gravitational field to retrieve 3D density variations within the Earth’s crust used in this work is based on the Bayesian algorithm presented in [2,7,8]. We just recall here its main concept and characteristics, leaving the interested reader to the above references and in particular to [8] for a detailed description of the theory behind the inversion. The most important feature of the algorithm, which differentiates it from the majority of other classical inversion schemes, is that it allows to invert at the same time for the 3D density distribution and for the geometry of the main geological units in a harmonised way. This is obtained by discretizing the study volume V in a set of small volumetric elements (made by right rectangular prisms), and introducing as unknowns for each element V i , an ancillary categorical variable defining the geological unit (here called label, L i ) and a continuous variable ρ i , describing the density of each element V i . In other words, given a volumetric element i, the label L i defines whether it is made by water, Plio-Quaternary sediments, Messinian salt, etc. In order to reduce as much as possible the space of the possible solutions the proposed inversion algorithm relies on the definition of a probabilistic a priori model. This a priori model is based on a set of maps of the depth of the main geological horizons with the corresponding predicted accuracy (used to define the a priori model in terms of label) and on information (expressed in terms of average value and range of variability) of the density within each layer. Given these data, the algorithm builds a complete 3D model of the considered volume equipped with a posterior probability distribution P ρ , L | Δ g o able to evaluate any perturbation of the initial model. Here Δ g o is the vector of the observed gravity anomalies, ρ is the vector of the densities ρ i and L is the vector of the labels L i . In particular, the posterior probability distribution is P ρ , L | Δ g o e ϕ with ϕ given by:
ϕ = 1 2 Δ g o A ρ T C Δ g 1 Δ g o A ρ + 1 2 ρ ρ ¯ T C 1 ρ ρ ¯ + 1 2 γ i = 1 N s 2 L i + + 1 2 λ i = 1 N j Δ i q 2 L i , L j i = 1 N χ ρ ¯ | 3 σ 2 ρ i
where A is the forward operator, ρ ¯ is the vector composed by the set of ρ ¯ = ρ ¯ L i , i.e., the average density of a given label , C is a diagonal matrix containing the set of σ 2 = σ 2 L i , i.e., the variance of the density of a given label , s 2 L i is a weight measuring the distance of the given model from the a priori one in terms of label, q 2 L i , L j is a weight measuring the clustering of two neighbours elements i , j and χ ρ ¯ | 3 σ 2 ρ i is the characteristic function, truncating the distribution between ρ ¯ 3 σ and ρ ¯ + 3 σ .
This function is highly non-linear, non-convex and contains a huge number of unknowns, so it is particularly difficult to be sampled or even plotted. For these reasons, the chosen solution is to find the set of parameters ρ i and L i that maximises the posterior probability. Wishing to estimate a maximum a posteriori of our posterior distribution, we face an optimization problem with the part of the variables that are discrete. The solution resorts to an application of a Gibbs sampler combined with a simulated annealing [9], as can be found in a large part of the literature; here the application of the method to the image analysis, with the seminal paper by [10], is worth mentioning. Basically, the final solution is computed as the perturbation of the initial model which maximises the a priori probability and which fits the gravity observations. Note that the algorithm will search the solution only in the range of the admissible models given the accuracy defined by a priori information. Here γ and λ are empirically set trying to equalise the order of magnitude of the different terms of Equation (1).
This Bayesian inversion scheme is particularly suitable for the two scenarios identified in Section 1 since it would allow for a proper modelling of the reliability and accuracy of different available information: Seismic profiles can locally constrain the geometries along a given line, 3D seismic can fix the model on a specific area and well logs data can locally reduce the uncertainty in the density distribution of each geological unit. Moreover, the main geological horizons’ depth together with the predicted accuracy, required to build the a priori model, can be easily computed, e.g., by means of a kriging procedure [11] or by an analogous Least Squares Collocation [12], directly applied to the interpreted seismic profiles. Both these interpolation algorithms will give similar results predicting maps coherent with the interpreted seismic lines and with accuracy worsening when moving away from seismic. As for the deepest layers, as mentioned before, in the absence of better information, we can suppose to use global crustal models (or even flat surfaces) but assigning them a very high uncertainty.

3. The Synthetic Case Study

In order to evaluate the quality of a gravity inversion algorithm the easiest way is to make use of synthetic models which, however, should be complex enough to simulate real-world problems. This is in fact the only possibility to quantitatively evaluate the closeness of the retrieved solution to the true model. For this reason, we have simulated, starting from a real set of interpreted seismic profiles in the Mediterranean Sea, a realistic synthetic model made of nine geological layers, namely bathymetry, Plio-Quaternary sediments, Messinian salt, Tortonian, Oligocene, Creataceous and Jurassic sediments, continental crust and mantle. In order to build a synthetic model, we start from analysing the available interpreted seismic profiles as follows:
1.
we compute a thickness of each geological unit;
2.
we estimate an empirical covariance of the thickness of each layer and interpolate it by means of Gaussian theoretical covariance function;
3.
for each layer, we compute a new random thickness with the same stochastic characteristics of the original dataset by classical triangular decomposition of the covariance matrix [13];
4.
we apply the same procedure to the first layer, namely the bathymetry;
5.
finally, starting from the new bathymetry and thicknesses, we build a set of simulated seismic profiles with the same stochastic characteristics as the initial one.
The main features of each layer in terms of mean, standard deviation (STD) and correlation length are reported in Table 1.
The modelled layers up to the base of the Jurassic sediments are characterised by a general increase of the STD with depth and by correlation lengths of the order of 25,000 m; the deepest layers, namely the basement and the Moho, instead show larger correlation lengths and STD. The synthetic seismic profiles have been interpolated by means of kriging procedure [11] on a set of regular grids with spatial resolution of 1500 m, also obtaining maps of the accuracy for each layer. Considering the densities, we assumed a constant density distribution for each layer, apart from the pre-Jurassic sediments and continental crust, where a linear trend with depth was added in the model. Densities of the model, together with their gradient in the z direction, are reported in Table 2.
Note that the constant density hypothesis is probably the largest approximation made to create the synthetic model. However, this simplification is not exploited in the inversion where an unknown density for each volumetric element has been considered. Starting from the horizons and densities, we have been able to build the 3D volume of our synthetic model up to a depth of 35,000 m and with a spatial resolution of 1500 m in the x and y directions and of 100 m in the z direction. As an example, we show the depth-to-basement and depth-to-Moho maps and two profiles of the 3D synthetic model (Figure 1).
The figure shows the complexity of the model which simulates a real-world scenario; in particular, black lines represent the discontinuities between the different geological units. The linear gradients added to the last sedimentary layer and to the continental crust are also recognisable from the figure. In the rest of the paper, we will refer to this model as the “true model” to distinguish it from the “a priori model”, which is the one used as input to the Bayesian inversion, and the “a posteriori model”, which is the result of the inversion itself.
To complete the synthetic dataset gravity observations at zero level have been simulated by applying a forward operator [14] to the modelled volume and adding a 3 mGal coloured noise with a correlation length of 4000 m (both the gravitational effect of the 3D volume and the observation error are shown in Figure 2).
In order to remove border effects in computing the forward model a border of 3 · 10 6 m was added around the volume. The border was created by extending in a constant way the limits of the volume according to what is described in [15]. This simulated dataset will be used as observation during the different inversions performed.

4. Results

We will present now, the results obtained by applying the inversion algorithm, briefly summarised in Section 2, to a set of tests derived by exploiting the synthetic model presented in Section 3. In particular, apart from an initial extremely simplified case that mimics a classical two-layer problem [16] (i.e., the Moho depth estimate supposing all the other information perfectly known), the two scenarios identified in Section 1 are treated. These last two cases show how the proposed gravity inversion can be profitably used to improve the modelling of a given volume also in the presence of seismic observations. In all the tests presented (except for the first one), we used the same densities as in the a priori models, obtained by adding a constant (random) value to the true density of each layer of the true model. The obtained average values (i.e., the elements of the vector ρ ¯ in Equation (1)) are reported in Table 3 together with the density accuracy in terms of STD used in the inversion (i.e., the square root of elements of the C matrix in Equation (1)).

4.1. Test 1: The Two-Layer Problem (Moho Estimate)

In the first investigated case study, we consider a simple situation in which everything is perfectly known apart from the Moho depth. Even if this is a not realistic scenario, it is useful to test the algorithm performance in a case in which the uniqueness of the solution is theoretically guaranteed, see, for instance, [6]. Moreover, this hypothesis, which is the basis of the two-layer model, is usually assumed in many scientific works [17,18]. In the absence of better information to build the a priori model, required by the inversion algorithm, we used a flat Moho at a depth of 28,000 m (i.e., the actual average depth of the Moho in the area). As for its uncertainty, we suppose a range of ± 3000 m. The a priori model for the two same profiles of Figure 1 is reported in Figure 3, where the dashed lines represent the accuracy of the a priori Moho.
The gravitational effect of the difference between the true and the a priori model is shown in Figure 4 on the left, together with the residual field after the inversion on the right, while in Figure 5 on the left the final estimated Moho is reported.
It can be seen that the algorithm has been able to almost recover the complete signal of the Moho, with differences between the true and the estimated depths of only 300 m (STD). For the sake of comparison, we also reported the Moho obtained by applying a classical Parker–Oldenburg inversion [16]. It can be seen how the Bayesian inversion is able to recover in a better way the true shape of the Moho. This is mainly due to approximations used in the forward and inverse operators of the Parker–Oldenburg inversion, which cut the high frequencies of the signal and by the fact that the Parker–Oldenburg inversion requires a constant density contrast, which is not the case for the considered model. The crustal density is in fact about 2900 kg/m 3 at a depth of 24,000 m (the minimum Moho depth) and more than 3000 kg/m 3 at a depth of 35,000 m (the maximum Moho depth) with a corresponding reduction in density contrast of about 25%. Looking to the gravity residuals of Figure 4, we can see that the Bayesian inversion over-fits the signal, showing a final residual of only 0.2 mGal (STD). This is possible because the algorithm is allowed to change the 3D density distribution inside each layer: The final model shows a difference in terms of density of about 4 kg/m 3 (STD). Note that the resultant density model of this first Bayesian inversion test is similar to the true one even if, in the a priori model, a not perfectly correct density is used. We have in fact performed a set of tests in which the a priori crustal density has been increased up to 30 kg/m 3 and the inversion is still able to recover the correct 3D density distribution (error STD smaller than 8 kg/m 3 ) and the proper Moho depth (with differences with respect to the true model of 500 m in terms of STD).

4.2. Scenario 1: Advanced Modelling of Deepest Layers (Basement and Moho Estimate)

In a second case study, a more realistic scenario has been tested in which only the geometries of the shallowest part of the model, i.e., up to the crystalline basement, are considered perfectly known, while densities and the deepest layers, namely the Moho and the basement, are unknown. This case has an important counterpart in real-life problems and can be considered as an example of Scenario 1 mentioned in Section 1, in which the shallowest layers are constrained by a 3D seismic survey and only the lower part of the volume is unknown. Again in the a priori model, we use as starting points for Moho and basement flat surfaces at a depth equal to the average value of the considered layers. As for their accuracy, we keep 3000 m for the Moho and we suppose 2000 m for the basement. Since the algorithm allows to move each layer in the range defined by the a priori accuracy, where the distance between the true surface and the a priori one is larger than the accuracy, the accuracy was increased up to the real surface, thus allowing the inversion to obtain the correct result. The new a priori model with the accuracy of each layer is shown in Figure 6 while the differences, in terms of gravity, between the true and the a priori models, and between the true and the a posteriori models, are shown in Figure 7.
As can be seen comparing Figure 7 with Figure 4, the most important gravimetric effect is still the one due to the Moho variations which dominates the two maps at the low-frequencies with an STD of about 11 mGal, while the effect due to the basement acts mainly at medium frequencies and is smaller (STD of about 7 mGal). Looking at the gravity residuals after the inversion, it can be seen that they are higher and more spatially correlated than in the previous scenario. This is basically due to the fact that in the current test the system has a greater degree of freedom and the most probable solution is found on the basis of the gravity residuals but also of the constraints imposed on the inversion, such as the closeness to the a priori model (third term in Equation (1)). Considering the retrieved layers, we can see from Figure 8 that the inversion has been able to qualitatively retrieve all the major elements of both the Moho and the basements. The retrieved Moho in particular is more regularised than the solution found in the previous scenario but still shows a good accuracy with differences with respect to the true model of only 900 m (STD). As for the basement the result shows, on the one hand an over-regularised low frequency and, on the other hand, the presence of a small high frequency signal, used mainly to compensate the observation error. In any case, also for the basement, the result contains all the main features and shows a difference with respect to the true model of 1400 m (STD).
We can conclude from this simulation that the Bayesian inversion algorithm is able to deal with the problem of estimating the deepest layers once the shallowest are fixed (e.g., by accurate 3D seismic surveys), qualitatively retrieving all the main features of the simulated basement and Moho and quantitatively reducing the model error from about 3000 m and 1700 m of the a priori model, to 1400 m and 900 m for the basement and Moho, respectively.

4.3. Scenario 2: Exploration in Frontier Areas (Shallowest Layers Fixed in Few Points)

In the last case study, we simulate an example of Scenario 2 of Section 1, in which the shallowest layers are constrained by a set of seismic profiles (here we have simulated a profile approximately every 30,000 m). In order to simulate this scenario, we computed a set of new geological horizons by perturbing the discontinuity surfaces used to build the true model. In particular, we applied a deformation to each surface with the same correlation length of the surface itself and with an amplitude varying point-wise and depending on the distance to the simulated seismic profile (the procedure used to simulate deformations is reported in Appendix A). The first layer, namely the bathymetry, is supposed known and fixed and the last layer, namely the Moho, is supposed not to be observed by the seismic profiles. In this way we are able to build a new model to be used as the a priori model in the inversion, which is equal to the true one only along the simulated seismic lines. As an example, we report in Figure 9 the deformation applied to the base of the Messinian salt layer and Moho, while in Figure 10 two sections of the a priori model together with the estimated uncertainty are shown. As can be observed, the resulting a priori model is coherent with the simulated seismic profiles but different from the true model. In Figure 11, the difference in terms of gravity between the true and the a priori models together with the residuals after the inversion are reported. It can be seen that, similarly to the other tests, the signal is dominated by the low frequencies coming from the deepest layers; however, in contrast to the other scenarios, the differences with respect to the true model also show medium and high frequencies. Two profiles of the estimated model are shown in Figure 12 together with the geological horizons of the true model (cyan lines) and of the a priori model (dashed lines). The inversion tends to move the layers of the a priori model toward the true ones. This is particularly evident for the Moho and basement but it is in general true also for the shallowest layers. This improvement is confirmed by Figure 13 where, as an example, the depth of the base of the Messinian salt for the true, a priori and a posteriori model is reported. It can be observed that the inversion is able to properly reduce the depth of the layer outside the domain covered by seismic profiles, whereas close to the profiles only minor corrections are introduced to the model. The improvement is confirmed as well by the statistics reported in Table 4, where the mean and the STD of the differences for each layer between the true and the a priori model is compared with the difference between the true and the a posteriori one. Here, we can see how the inversion is able to catch the main shape of each layer (decreasing the STD of the differences for all the layers). As for the average difference the situation is more complex and the inversion seems to be able to improve only the shallowest layers.
In any case, also for this scenario, we can conclude that the algorithm is able to improve, away from the seismic profiles, the geophysical modelling by halving the distances from the true model. It can also be observed that the smallest relative improvement is found for the base of the Plio-Quaternary sediments, i.e., the shallowest layer (we recall the fact that the bathymetry is considered known and fixed). This is due to the fact that this layer is used by the inversion algorithm to compensate the observation error. Looking at the gravity residuals, we see that the model overfits the data, with a final residual of only 0.6 mGal (STD), which is much lower than the observation error. Finally, for the sake of completeness, we report in Figure 14 a horizontal section of the density distribution at a depth of 5000 m of the true, a priori and a posteriori models. It can be seen that, as expected, the main density variations are due to the presence of the different layers. However, in the a posteriori model, density variations inside each horizon can also be found. For instance, looking to Cretaceous sediments, which is one of the most important layers of the considered section, densities between 2470 kg/m 3 and 2484 kg/m 3 (coherent with the a priori assumption in Table 3) are found.

5. Conclusions

Within the presented work, we developed a realistic synthetic case study, which can be in general exploited in order to test gravity inversion algorithms. In detail, the developed dataset simulates a 3D volume in the Mediterranean area over a region of 3 · 10 5 m × 2 · 10 5 m up to a depth of 35,000 m and with a spatial resolution of 1500 m in the x and y directions. The model has been used to assess the accuracy of a Bayesian inversion algorithm on a set of typical scenarios ranging from the classical two-layer, used in general for scientific investigations, where the only unknown is the discontinuity between two units of known density, to a complex case study where all the layers and densities are considered unknown and only partially constrained by the presence of simulated seismic profiles. Actually, the simulated case study proved that, when the simple two-layer problem is tackled, with the aim of estimating the Moho, the proposed inversion is able to retrieve almost perfect results, with differences between the simulated and resulting models of only 300 m (STD), also in the presence of an important observation error and when the density contrast between the two layers is not perfectly known. We also assess the quality of the Bayesian inversion on two more realistic scenarios: In the former, the shallowest layers’ geometries are known and only the deepest layers should be retrieved, while in the latter also the shallowest layers are known only concerning a few profiles. The results obtained show that in both these situations, the Bayesian inversion improves the 3D modelling of the considered volume. For the first scenario, the algorithm was able to catch all the main features of the basement and Moho, retrieving the surfaces with an accuracy of 1400 m and 900 m, respectively. Similar results are obtained for the last simulated scenario with a general improvement of geometries of the initial 3D model. Here, the accuracy of the final solution ranges, in terms of geometry, between 400 m and 800 m. In terms of gravity field, it is observed that the Bayesian inversion tends to overfit the gravity data with residuals after the inversion always smaller than 1 mGal. This is due to a large value degree of freedom of the performed inversion. We can conclude that the proposed algorithm can safely be used also in a complex scenario to improve 3D modelling of the subsurface even when seismic profiles or even 3D seismic data are available.

Author Contributions

D.S. supervised the study; M.C. and D.S. conceived and performs the study; D.S. and M.C. performs the bibliographic research; M.C. and D.S. wrote the paper. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Simulation of Non-Stationary Random Fields

In order to simulate the starting point for the case study investigated in Section 4.3, we have to compute, starting from the maps of depth to the geological horizons of the true model, a set of alternative maps, which respect the available information, namely which are equal to the true model along the simulated seismic lines and are coherent from the stochastic point of view, with the initial maps. The problem can be ultimately brought back to the simulation of non-stationary random fields. In this work, we will apply the following scheme:
1
we compute the empirical variogram, using the available interpreted seismic profiles data as observations;
2
we fit the empirical variogram by means of a proper theoretic variogram function;
3
knowing the theoretical variogram, we can estimate a map of the depth of the considered geological unit and its predicted accuracy by means of a kriging solution;
4
we randomly select on the obtained map a set of points;
5
we extract a random value on each of the points obtained at step 4 from a Gaussian probability function with a STD equal to the one predicted by the kriging map on that specific location;
6
we spatially correlate the extracted sample by applying a kriging procedure on the random points with the same variogram function estimated at step 2.
We report here, for the sake of clarity an example of the application of the above scheme applied to the first layer (namely the bathymetry) of the synthetic data presented in this work.
The starting point is a set of interpreted seismic profiles, observing the depth of the considered layer (Figure A1a).
From these observations the empirical variogram is computed and fitted by means of a stable variogram model (Figure A1b). At this point, by means of an ordinary kriging, we compute the depth of the considered layer on a regular grid over the whole study area together with its predicted accuracy (see Figure A1c,d). Note that the available profiles only partially covered the area of interest, as is often the case with this kind of problem. Now, we randomly select n points with n = Δ x Δ y d 0 , where Δ x and Δ y are the total size of the considered area in the x and y directions, respectively, while d 0 is the correlation length as obtained by the variogram itself. At these selected points, we generate a white noise with STD accordingly to the accuracy predicted by the kriging interpolation. The variogram function, used to interpolate seismic profiles, is than used to correlate the noise thus obtaining a possible perturbation of the initial considered model (see Figure A1e). A new bathymetry model consistent with the seismic data (accordingly to the point-wise predicted accuracy) is obtained by summing up the initial map with the random perturbation (Figure A1f). As it can be observed from Figure A1 the specific realisation is coherent with available information: it has in fact a correlation length similar to the one of initial seismic information and a variance locally consistent with the accuracy predicted by the kriging.
Figure A1. Available seismic profiles observing the bathymetry over the considered area (a). Empirical and theoretical variograms (b). Estimated bathymetry from the seismic profiles (c) and predicted accuracy in terms of STD (d). Random realisation of possible deformation (e). Possible realisation of bathymetry coherent with seismic data (f).
Figure A1. Available seismic profiles observing the bathymetry over the considered area (a). Empirical and theoretical variograms (b). Estimated bathymetry from the seismic profiles (c) and predicted accuracy in terms of STD (d). Random realisation of possible deformation (e). Possible realisation of bathymetry coherent with seismic data (f).
Geosciences 11 00467 g0a1

References

  1. Reguzzoni, M.; Sampietro, D. GEMMA: An Earth crustal model based on GOCE satellite data. Int. J. Appl. Earth Obs. Geoinf. 2015, 35, 31–43. [Google Scholar] [CrossRef]
  2. Sampietro, D.; Capponi, M. Gravity for lIthosphere Architecture Determination and Analysis: The Central Eastern Mediterranean Case Study. Geophys. Prospect. 2021. [Google Scholar] [CrossRef]
  3. Ebbing, J.; Braitenberg, C.; Götze, H.J. Forward and inverse modelling of gravity revealing insight into crustal structures of the Eastern Alps. Tectonophysics 2001, 337, 191–208. [Google Scholar] [CrossRef]
  4. Zingerle, P.; Pail, R.; Gruber, T.; Oikonomidou, X. The combined global gravity field model XGM2019e. J. Geod. 2020, 94, 66. [Google Scholar] [CrossRef]
  5. Laske, G.; Masters, G.; Ma, Z.; Pasyanos, M. Update on CRUST1. 0—A 1-degree global model of Earth’s crust. Geophys. Res. Abstr. 2013, 15, 20132658abstrEGU. [Google Scholar]
  6. Sampietro, D.; Sanso, F. Uniqueness theorems for inverse gravimetric problems. In VII Hotine-Marussi Symposium on Mathematical Geodesy; Springer: Berlin/Heidelberg, Germany, 2012; pp. 111–115. [Google Scholar]
  7. Marchetti, P.; Sampietro, D.; Capponi, M.; Rossi, L.; Reguzzoni, M.; Porzio, F.; Sansò, F. Lithological constrained gravity inversion. A Bayesian approach. In Proceedings of the 81st EAGE Conference and Exhibition 2019, London, UK, 3–6 June 2019; pp. 1–5. [Google Scholar]
  8. Sansó, F.; Sampietro, D. Analysis of the Gravity Field: Direct and Inverse Problems; Birkhäuser: Basel, Switzerland, 2021. [Google Scholar]
  9. Robert, C.; Casella, G. Monte Carlo Statistical Methods; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2013. [Google Scholar]
  10. Geman, S.; Geman, D. Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images. IEEE Trans. Pattern Anal. Mach. Intell. 1984, PAMI-6, 721–741. [Google Scholar] [CrossRef] [PubMed]
  11. Wackernagel, H. Multivariate Geostatistics: An Introduction with Applications; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2003. [Google Scholar]
  12. Sampietro, D.; Capponi, M.; Mansi, A.; Gatti, A.; Marchetti, P.; Sansò, F. Space-Wise approach for airborne gravity data modelling. J. Geod. 2017, 91, 535–545. [Google Scholar] [CrossRef] [Green Version]
  13. Davis, M.W. Production of conditional simulations via the LU triangular decomposition of the covariance matrix. Math. Geol. 1987, 19, 91–98. [Google Scholar]
  14. Nagy, D. The gravitational attraction of a right rectangular prism. Geophysics 1966, 31, 362–371. [Google Scholar] [CrossRef]
  15. Sampietro, D.; Capponi, M. Practical tips for 3D regional gravity inversion. Geosciences 2019, 9, 351. [Google Scholar] [CrossRef] [Green Version]
  16. Oldenburg, D.W. The inversion and interpretation of gravity anomalies. Geophysics 1974, 39, 526–536. [Google Scholar] [CrossRef]
  17. Ebadi, S.; Barzaghi, R.; Safari, A.; Bahroudi, A. Evaluation of different gravimetric methods to Moho recovery in Iran. Ann. Geophys. 2019, 62, GD561. [Google Scholar] [CrossRef]
  18. Reguzzoni, M.; Sampietro, D.; Rossi, L. The gravimetric contribution to the Moho estimation in the presence of vertical density variations. Rend. Lincei Sci. Fis. E Nat. 2020, 31, 69–81. [Google Scholar] [CrossRef]
Figure 1. Depth to basement and Moho depth of the synthetic model (upper images). Example of two profiles (AA’ and BB’) on the synthetic 3D model. Solid black lines represent discontinuity between different geological units, colours represent density in kg/m 3 . W stands for oceanic water, P Plio-Quaternary sediments, M Messinian salt, T Tortonian sediments, O Oligocene sediments, C Cretaceous sediments, J Jurassic sediemnts, PJ Pre-Jurassic sediments.
Figure 1. Depth to basement and Moho depth of the synthetic model (upper images). Example of two profiles (AA’ and BB’) on the synthetic 3D model. Solid black lines represent discontinuity between different geological units, colours represent density in kg/m 3 . W stands for oceanic water, P Plio-Quaternary sediments, M Messinian salt, T Tortonian sediments, O Oligocene sediments, C Cretaceous sediments, J Jurassic sediemnts, PJ Pre-Jurassic sediments.
Geosciences 11 00467 g001
Figure 2. Gravitational effect of the synthetic model (left), and simulated observation error (right).
Figure 2. Gravitational effect of the synthetic model (left), and simulated observation error (right).
Geosciences 11 00467 g002
Figure 3. Example of two profiles (AA’ and BB’ in Figure 1) for the two-layer problem inversion. Solid black lines represent discontinuity between different geological units, dashed lines represent uncertainty, colours represent density in kg/m 3 .
Figure 3. Example of two profiles (AA’ and BB’ in Figure 1) for the two-layer problem inversion. Solid black lines represent discontinuity between different geological units, dashed lines represent uncertainty, colours represent density in kg/m 3 .
Geosciences 11 00467 g003
Figure 4. Two-layer problem inversion: Residuals between the gravitational effect of the true and a priori model (left) and of the true and the a posteriori model (right).
Figure 4. Two-layer problem inversion: Residuals between the gravitational effect of the true and a priori model (left) and of the true and the a posteriori model (right).
Geosciences 11 00467 g004
Figure 5. Two-layer problem inversion: Moho estimated by the Bayesian inversion (left) and from classical Parker–Oldenburg inversion (right).
Figure 5. Two-layer problem inversion: Moho estimated by the Bayesian inversion (left) and from classical Parker–Oldenburg inversion (right).
Geosciences 11 00467 g005
Figure 6. Example of two profiles (AA’ and BB’ in Figure 1) for Scenario 1 inversion:. Solid black lines represent discontinuity between different geological units, dashed lines represent uncertainty, colours represent density in kg/m 3 .
Figure 6. Example of two profiles (AA’ and BB’ in Figure 1) for Scenario 1 inversion:. Solid black lines represent discontinuity between different geological units, dashed lines represent uncertainty, colours represent density in kg/m 3 .
Geosciences 11 00467 g006
Figure 7. Scenario 1 inversion: Difference between the gravitational effect of the true and a priori model (left) and of the true and the a posteriori model (right).
Figure 7. Scenario 1 inversion: Difference between the gravitational effect of the true and a priori model (left) and of the true and the a posteriori model (right).
Geosciences 11 00467 g007
Figure 8. Scenario 1 inversion: Estimated depth to basement (left) and Moho depth (right).
Figure 8. Scenario 1 inversion: Estimated depth to basement (left) and Moho depth (right).
Geosciences 11 00467 g008
Figure 9. Differences between true and a priori model of Scenario 2 inversion for the base of the Messinian layer (left) and for the Moho (right). Black lines represent the simulated seismic profiles.
Figure 9. Differences between true and a priori model of Scenario 2 inversion for the base of the Messinian layer (left) and for the Moho (right). Black lines represent the simulated seismic profiles.
Geosciences 11 00467 g009
Figure 10. Example of two profiles (AA’ and BB’ in Figure 1) for Scenario 2 inversion. Solid black lines represent discontinuity between different geological units, dashed lines represent uncertainty, colours represent density in kg/m 3 .
Figure 10. Example of two profiles (AA’ and BB’ in Figure 1) for Scenario 2 inversion. Solid black lines represent discontinuity between different geological units, dashed lines represent uncertainty, colours represent density in kg/m 3 .
Geosciences 11 00467 g010
Figure 11. Scenario 2 inversion: Difference between the gravitational effect of the true and a priori model (left) and of the true and the a posteriori model (right).
Figure 11. Scenario 2 inversion: Difference between the gravitational effect of the true and a priori model (left) and of the true and the a posteriori model (right).
Geosciences 11 00467 g011
Figure 12. Example of two profiles (AA’ and BB’ in Figure 1) of the estimated model for Scenario 2 inversion. Cyan lines represent the geological horizons of the true model, black lines represent the horizons estimated by the inversion, dashed black lines represent the horizons of the a priori model and colours represent density in kg/m 3 .
Figure 12. Example of two profiles (AA’ and BB’ in Figure 1) of the estimated model for Scenario 2 inversion. Cyan lines represent the geological horizons of the true model, black lines represent the horizons estimated by the inversion, dashed black lines represent the horizons of the a priori model and colours represent density in kg/m 3 .
Geosciences 11 00467 g012
Figure 13. Base of the Messinian salt from the true (a), a priori (b) and a posteriori model (c).
Figure 13. Base of the Messinian salt from the true (a), a priori (b) and a posteriori model (c).
Geosciences 11 00467 g013
Figure 14. Density distribution at a depth of 5000 m from the true (a), a priori (b) and a posteriori models (c).
Figure 14. Density distribution at a depth of 5000 m from the true (a), a priori (b) and a posteriori models (c).
Geosciences 11 00467 g014
Table 1. Statistics of the simulated geological horizons.
Table 1. Statistics of the simulated geological horizons.
LayerAverage Depth [m]Depth STD [m]Correlation Length [m]
Bathymetry201073024,000
Base Plio-Quaternary222071023,000
Base Mess. Salt355097026,000
Base Tortonian4210146029,000
Base Oligocene4990155026,000
Base Creataceous5800148028,000
Base Jurassic6920168026,000
Basement11,610296046,000
Moho28,360173031,000
Table 2. Density values used to generate the synthetic model.
Table 2. Density values used to generate the synthetic model.
LayerAverage Density [kg/m 3 ]Density Gradient [kg/m 4 ]
Water10300
Plio-Quaternary22200
Mess. Salt21600
Tortonian22600
Oligocene24000
Creataceous24800
Jurassic25500
Pre-Jurassic26200.005
Continental Crust26700.012
Mantle33000
Table 3. Density values and accuracy of the a priori models for Scenario 1 and 2 inversions (Section 4.2 and Section 4.3).
Table 3. Density values and accuracy of the a priori models for Scenario 1 and 2 inversions (Section 4.2 and Section 4.3).
LayerAverage Density [kg/m 3 ]Density STD [kg/m 3 ]
Water10300
Plio-Quaternary220610
Mess. Salt217210
Tortonian225110
Oligocene241110
Creataceous247310
Jurassic254910
Pre-Jurassic261422
Continental Crust268031
Mantle332131
Table 4. Scenario 2 inversion: Statistics of the differences between true, a priori and a posteriori models.
Table 4. Scenario 2 inversion: Statistics of the differences between true, a priori and a posteriori models.
LayerA-PRIORI vs. True ModelA-POSTERIORI vs. True Model
Mean [m]STD [m]Mean [m]STD [m]
Water0000
Base Plio-Quaternary0100085
Base Mess. Salt15065947378
Base Tortonian32374294393
Base Oligocene316883168409
Base Cretaceous370861309418
Base Jurassic250930707517
Basement801718267818
Moho3611736837825
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Sampietro, D.; Capponi, M. Seismic Constrained Gravity Inversion: A Reliable Tool to Improve Geophysical Models Away from Seismic Information. Geosciences 2021, 11, 467. https://doi.org/10.3390/geosciences11110467

AMA Style

Sampietro D, Capponi M. Seismic Constrained Gravity Inversion: A Reliable Tool to Improve Geophysical Models Away from Seismic Information. Geosciences. 2021; 11(11):467. https://doi.org/10.3390/geosciences11110467

Chicago/Turabian Style

Sampietro, Daniele, and Martina Capponi. 2021. "Seismic Constrained Gravity Inversion: A Reliable Tool to Improve Geophysical Models Away from Seismic Information" Geosciences 11, no. 11: 467. https://doi.org/10.3390/geosciences11110467

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop