1. Introduction
Gravity data are the source of important information about the structure and composition of the earth’s crust, which is critical in mineral resource exploration. Three-dimensional inversion of the regional gravity data represents a challenging task due to the large size of the inversion domain and intrinsic nonuniqueness and instability of the potential field inversion. One can reduce the ambiguity of inversion by considering some a priori information about the internal structure of the earth’s crust. For example, the presence of the sediment-to-basement interface can be used as an a priori constraint on the inversion. When the depth to the basement is known from seismic survey data, the basement surface can be incorporated into the inversion algorithm. Otherwise, we can subdivide the gravity inverse problem into two phases. First, we apply the depth-to-basement inversion to generate the initial model of the density interface. Then, we use this model as a priori information for the full voxel-based gravity inversion in the second phase.
In this paper, we discuss first the method of the inversion to the depth of the density interface based on the 3D Cauchy-type integral representation of the gravity field. We also briefly describe the principles of the voxel-based gravity inversion. Finally, we apply the developed method to interpret the Bouguer gravity anomaly observed over the central part of Utah. The preliminary results were presented as a conference paper at the international conference IMAGE 2021 [
1].
The U.S. Geological Survey (USGS) has been developing a National Crustal Model of the continental United States [
2]. As a part of this project, regional maps of the depth to the basement were produced. The depth of the basement in those publications was defined as “the depth to the top surface of igneous or metamorphic rocks; rocks above this depth are either sediments or sedimentary rocks.” However, this definition did not provide a unique characteristic of the basement because sedimentary rocks could experience some metamorphism, resulting in a significant increase in density. In this paper, we use the gravity data to determine the depth of the strong density contrast in the subsurface, which in some cases can be associated with the surface of the crystalline basement. We demonstrate that our inversion results correlate well with the depth to the basement determined by USGS using the analysis of available geological and geophysical data [
3], without running the inversion of the entire gravity dataset.
2. Application of 3D Cauchy-Type Integrals in Depth-to-Basement Gravity Inversion
Inversion of gravity data for the depth to the basement requires an application of practical modeling algorithms designed for this specific model with strong density contrast. In this situation, it is not suitable to use the voxel-type discretization of the density models. Indeed, one would require applying a very fine discretization with a large number of cells to represent a density contrast accurately. Therefore, the natural approach is based on only discretization of the sediment-to-basement interface. This approach was initially introduced in [
4,
5] and practically implemented in [
6,
7]. For completeness, we briefly discuss this approach below, following [
7]. The interested reader can find the complete mathematical description of the application of 3D Cauchy-type integrals in depth-to-basement gravity inversion in [
7]. In the framework of this approach, the gravity response of the density contrast surface is represented in the form of the surface integral over the density contrast interface. Remarkably, this surface integral has a special form of the 3D Cauchy-type integral, which extends to 3D cases all the properties of the classical Cauchy integral of the theory of functions of a complex variable [
6]. This allows applying the reach properties of these integrals to the effective numerical representation of the gravity field of the density contrast surface.
We consider a model of the sediment-to-basement interface represented by the surface,
S, with a density contrast,
described by the equation
(
Figure 1). We also assume that
S coincides with the horizontal plane
P,
, at infinity. Note that in this case, function
represents a distance between surface
S and plane
P. The gravity field,
, generated by the density interface,
S, can be represented as a 3D Cauchy-type surface integral over this surface [
5,
6,
7] as follows:
where
is a unit vector directed upward along the vertical axis;
γ is the universal gravitational constant, and
is the radius-vector of the observation point.
Cauchy-type integral,
is a surface integral of
over the density contrast surface,
S:
where
is the radius-vector of the integration point on surface
S.
We should note that the numerical analog of Equation (1) is based only on a discretization of the density contrast surface. At the same time, the standard forward modeling and inversion methods require the volume discretization of the modeling domain. Thus, using this approach allows reducing the gravity forward modeling from a 3D to a 2D integration domain. This significantly reduces the computing resources necessary for 3D inversion for the density contrast surface.
Following [
6], we can discretize the Cauchy-type integral (Equation (2)) by dividing the horizontal plane
P into a rectangular grid of
cells,
, with the constant discretization of
Δx and
Δy in the
x and
y directions, respectively. We assume that within each cell,
(
k = 1,2,… ), the corresponding part of the density contrast surface, can be represented by an element of a plane described by the following equation:
where
denotes the center of the cell,
, and
In this case, Equation (1) for vertical component,
, of the gravity field takes the form:
where
.
Using the discrete model parameters introduced above and discrete gravity data,
, we can approximate the forward modeling operator for the gravity field produced by the density contrast surface
S as follows:
where
and
Note that the accuracy of numerical approximation (Equation (6)) of integral Equation (4) was examined in detail in [
6]. One can find the proper recommendations for determining the discretization parameters in [
6] as well. Therefore, we can write the forward modeling problem for gravity field using the matrix representation as follows:
Here m is a vector of model parameters (the elevations, h(, ) of the order ; d is a vector of observed data, , of the order , and is a rectangular matrix of the size ×, formed by the gravity field kernels, Equation (6).
Note that the operator is nonlinear because data kernels (Equation (6)) and the corresponding matrix, , depend on m. Therefore, the inverse problem (Equation (7)) is also nonlinear, and the Fréchet derivative, F, of operator depends on the model parameters.
Following the standard framework of the regularization method, we reduce the solution of the nonlinear inverse problem (Equation (7)) to minimization of the corresponding parametric functional:
where
is the data weighting matrix;
is the a priori model of the density contrast surface; and
is a diagonal matrix of the model parameter weights based on integrated sensitivity:
This selection of model weights provides an equal sensitivity of the observed data to the elements of the density contrast surface located at different depths and horizontal positions [
8,
9].
Matrix
is also a diagonal matrix of the minimum support stabilizing functional providing focusing inversion [
7,
8]:
We solve the minimization problem (Equation (8)) using the re-weighted regularized conjugate gradient (RRCG) method. A detailed description of this method can be found in [
8,
9].
3. Voxel-Type 3D Gravity Inversion
The conventional voxel-type 3D inversion is based on discretization of the subsurface in rectangular prisms. The following formula describes the gravity forward modeling problem:
where
r is the source location;
r′ is the receiver location;
is the density distribution within some domain
D; and
γ is the universal gravitational constant.
The discrete form of forward modeling operator (Equation (11)) is obtained by dividing domain
D into
small rectangular cells,
, and assumes that the density is constant within each cell,
. As a result, Equation (11) for gravity field takes the following form:
Assuming a relatively small size of rectangular cells,
, we can use the point-mass approximation, which dramatically speeds up the processing time while yielding very accurate results [
10].
Let us assume that the coordinates of the cell centers are
=(
,
,
),
l = 1,…,, and the cell sides are Δ
x, Δ
y, Δ
z. Also, we have a discrete number of observation points
,
n = 1, …, . Using discrete model parameters and discrete data, we can present the forward modeling operator for gravity field (Equation (12)) as follows:
where the gravity field kernel,
is calculated by the following formula (see Equation (12)):
and
One can find a detailed study of the accuracy of the point mass approximation (Equation (14)) in [
10]. Using the discrete model parameters introduced above, we can approximate the forward modeling operator for the gravity field produced by the volume distribution of density as follows:
Here m is a vector of model parameters (densities, ) of the order ; d is a vector of observed data; of the order ; and is a rectangular matrix of the size ×, formed by the gravity field kernels, as in Equation (12). Note that the number of discretization cells, , in the voxel-type inversion is obviously significantly greater than the number of unknown elevations, in the case of the density contrast surface inversion: ≫.
The inverse problem (Equation (15)) can be solved by minimizing the same Tikhonov parametric functional as in Equation (8). The only difference is that vector m of the model parameters is formed by the density values within prismatic cells of the volume discretization grid, m = ρ, and matrix is substituted for matrix . In addition, in the case of the voxel-type inversion, we select the two-layered models with the density contrast surface determined on the first step of the inversion as an a priori density model, However, in this case, we consider a relatively weak density contrast. Thus, the 3D voxel-type inversion is guided by the inversion results for the density contrast surface. At the same time, the volume distribution of the density is adjusted to fit the observed data better.
4. Inversion of Bouguer Gravity Anomaly in the Utah State
We have applied the developed method of two-step inversion to the complete Bouguer gravity anomaly dataset over Utah (
Figure 2). The United States Geological Survey (USGS) extracted these data from the gravity database maintained by the National Geophysical Data Center and combined them with the USGS data and several university theses and dissertations. The observed gravity data were reduced to the Bouguer anomaly using the 1967 gravity formula and a reference density of 2.67 g/cc. USGS also calculated the terrain corrections radially outward from each station to a distance of 167 km, using a method developed by Plouff (USGS Open-file Report 77-535). The minimum curvature technique was used to convert the data to a 1 km grid (
https://pubs.usgs.gov/of/1998/ofr-98-0761/utah_boug.html (accessed on 1 March 2022)).
In Utah, USGS has reported a presence of a relatively strong density contrast surface associated with the Mesozoic basement [
3]. However, it is challenging to recover this density contrast surface using a conventional voxel-type 3D inversion. The problem is that traditional gravity field inversion usually generates smooth images of the density distribution, which cannot resolve a sharp density contrast associated with the depth to the basement. As discussed above, to overcome this problem, we use a novel approach that involves a two-step inversion. In the first step, we invert the Bouguer gravity anomaly to the depth-to-basement interface using Cauchy-type integral representation discussed above. In the second step, we apply the voxel-type 3D gravity inversion using the depth-to-basement model as a soft constraint. This approach has produced one of the first 3D density distributions of the crustal model in Utah based on the Bouguer gravity anomaly data.
This paper focuses on inverting the Bouguer gravity anomaly, selected in the central part of the State of Utah (
Figure 2 and
Figure 3). A blue rectangle shown in the center of the map (
Figure 2) outlines the selected area of interest (AOI).
First, we apply depth-to-basement inversion. In this case, the USGS depth to Mesozoic basement (
Figure 4) serves an a priori model [
3]. This map was constructed using two parameters characterizing the thickness of key geologic layers: (1) the vertical size of unconsolidated sediments (e.g., depth to bedrock) and (2) the depth to the basement. These parameters were estimated based on different methods, e.g., seismic reflection, well data, gravity, magnetic surveys, etc.
We present in
Figure 5 the density contrast surface produced by the depth-to-basement inversion with the density contrast equal to 0.3 g/cc. This value was determined based on known data and the best convergence of the inversion. The iterative inversion ran until the L2 misfit between the observed and predicted data reached a level of approximately 7%.
We should note that we do not know if we have reconstructed the actual depth of the crystalline basement. The reason is that tectonic processes have created highly complex geology in the western United States (e.g., [
3]), which makes it difficult to identify a specific geologic layer associated with the basement. Therefore, the surface we have recovered should be considered a density contrast surface in the crust that underlies the unconsolidated sediment thickness. At the same time, this surface may coincide with the surface of the crystalline basement in some areas, while it is difficult to draw a specific conclusion in other areas.
Comparison between the USGS model (
Figure 4) and gravity contrast inverse model (
Figure 5) shows an amazing similarity between these two maps, confirming the USGS map’s validity.
Next, we use a conventional 3D volume integral representation to invert the Bouguer gravity anomaly, filtered with a simple plane removal (
Figure 6). The 3D inversion domain was discretized into 99 × 79 × 20 = 156,420 rectangular cells. We run inversion both without and with the a priori model in this step. The two-layered model with the density contrast surface found in the previous step (
Figure 5) and a relatively weak contrast of 0.1 g/cc was used as an a priori model. At the same time, the voxel-type density inversion applied on the second phase updated the density everywhere in the inversion domain, including the volumes both over and under the density contrast surface determined in the first stage. The iterative inversions in both cases (without and with the a priori model) reached the L2 misfit between the observed and predicted data by less than 5%.
Figure 7 and
Figure 8 demonstrate that both inversions reproduce the observed data well.
It is interesting to compare the 3D density models produced by two inversions. For example, in
Figure 9 and
Figure 10, we show the vertical sections along profile BB′ (see
Figure 6) of 3D density models produced by the inversion without and with an a priori density contrast model.
Figure 11 and
Figure 12 represent the corresponding vertical sections along profile DD′, respectively. The solid line in
Figure 9,
Figure 10,
Figure 11 and
Figure 12 indicates the outline of the a priori density contrast surface used in the inversion. These figures show that by including the a priori density contrast model in the inverse problem solution, we can better delineate the bottom of unconsolidated sediments and the top of the basement. In addition, the produced inverse density model corresponds well to the depth of Mesozoic basement constructed by USGS (the dashed lines in
Figure 9,
Figure 10,
Figure 11 and
Figure 12). This indicates the ability of the guided inversion to adjust the a priori model to fit the observed data better. In short, the guided inversion is the data-driven approach.
Figure 13 presents combined vertical sections along all five profiles shown in
Figure 6 of a 3D density model produced by the inversion with an a priori density contrast model. One can see in this figure that the density inverse model provides a fascinating picture of the deep geological structure in the area of investigation. First of all, we observe in
Figure 13 a significant density deficit under the Great Salt Lake, indicating a large depression formed by the steeply deeping faults serving as the flanks of surrounding mountains. The images also show more dense roots under the Lakeside, Cedar, Stansbury, and Oquirrh mountains and the Wasatch Range, and under the Stansbury and Antilope islands. In
Figure 14,
Figure 15 and
Figure 16, we have also plotted the horizontal sections of the density inverse model at elevations of 0.7 km, 0.2 km, and −0.3 km above sea level.
Thus, the produced 3D density model of the area of interest in central Utah provides important information about the complex geology in the area and the thickness of unconsolidated sediments.