Three-Dimensional Gravity Inversion in the Presence of the Sediment-Basement Interface: A Case Study in Utah, USA

: We introduce a novel approach to three-dimensional gravity inversion in the presence of the sediment-basement interface with a strong density contrast. This approach makes it possible to incorporate the known information about the basement depth in the inversion. It also allows the user to determine the depth-to-basement in the initial inversion phase. One can then use this interface to constrain the ﬁnal inversion phase. First, the inversion generates the depth-to-basement model based on the 3D Cauchy-type integral representation of the gravity ﬁeld. Then, in the second phase, full 3D voxel-type inversion applies the depth-to-basement model determined in the ﬁrst phase as an a priori constraint. We use this approach to the 3D inversion of the Bouguer gravity anomaly data observed in Utah, USA. The results of inversion generated a 3D density model of the top layers of the earth’s crust, including unconsolidated sediments and the top of the crystalline basement.


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.

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, ρ 0 , described by the equation z = h(x, y) − H 0 ( Figure 1). We also assume that S coincides with the horizontal plane P, z = −H 0 , at infinity. Note that in this case, function h(x, y) represents a distance between surface S and plane P. The gravity field, g, 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 d z is a unit vector directed upward along the vertical axis; γ is the universal gravitational constant, and r is the radius-vector of the observation point. Cauchy-type integral, C S r , h(x, y)d z , is a surface integral of ϕ = h(x, y)d z over the density contrast surface, S: where r 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 N m cells, P k , with the constant discretization of ∆x and ∆y in the x and y directions, respectively. We assume that within each cell, P k (k = 1,2, . . . N m ), the corresponding part of the density contrast surface, can be represented by an element of a plane described by the following equation: where (x k , y k ) denotes the center of the cell, P k , and h (k) = h(x k , y k ).
In this case, Equation (1) for vertical component, g z , of the gravity field takes the form: Using the discrete model parameters introduced above and discrete gravity data, g z (r n ), 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 (k) =h(x k , y k ) of the order N m ; d is a vector of observed data, g z , of the order N d , and A h is a rectangular matrix of the size N d × N m , formed by the gravity field kernels, Equation (6).
Note that the operator A h (m) is nonlinear because data kernels (Equation (6)) and the corresponding matrix, A h , depend on m. Therefore, the inverse problem (Equation (7)) is also nonlinear, and the Fréchet derivative, F, of operator A h (m) 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 W d is the data weighting matrix; m apr is the a priori model of the density contrast surface; and W m 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 W e 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].

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; ρ(r) 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 L m small rectangular cells, D l , and assumes that the density is constant within each cell, ρ l . As a result, Equation (11) for gravity field takes the following form: Assuming a relatively small size of rectangular cells, D l , 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 r l =(x l ,y l ,z l ), l = 1, . . . ,L m , and the cell sides are ∆x, ∆y, ∆z. Also, we have a discrete number of observation points r n = (x n , y n , z n ), n = 1, . . . , N d . 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, A ρ nl , is calculated by the following formula (see Equation (12)): and r nl = (x l − x n ) 2 + (y l − y n ) 2 + z l 2 .
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, ρ l ) of the order L m ; d is a vector of observed data; g z , of the order N d ; and A ρ is a rectangular matrix of the size N d × L m , formed by the gravity field kernels, as in Equation (12). Note that the number of discretization cells, L m , in the voxel-type inversion is obviously significantly greater than the number of unknown elevations, N m , in the case of the density contrast surface inversion: N m L m . 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 A ρ is substituted for matrix A h . 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, m apr . 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.

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-tobasement 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 (Figures 2 and 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%.   It is interesting to compare the 3D density models produced by two inversions. For example, in Figures 9 and 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.   Figures 9-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 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 Figures 14-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.

Conclusions
The gravity data are often used for regional geological study. One of the major geological structures affecting the observed data is the surface of the consolidated basement. The goal of gravity inversion is two-fold-to determine the depth-to-basement and find the density distribution within the sedimental layer and the underlying basement. Therefore, it is reasonable to apply the inversion in two stages in this situation. In the first stage, one can run inversion for the depth-to-basement only to determine this dominant geological structure in the survey area. Then, we can use the produced density contrast surface as a soft constraint (a priori model) for the conventional 3D volume inversion in the second stage.
We have applied the developed two-stage approach to the inversion of Bouguer gravity anomaly data in central Utah. The inversion produced one of the first 3D models of the upper crust in central Utah. This model provides detailed information about the density distribution in the top layers of the earth crust under the Great Salt Lake and the surrounding mountains. We plan to extend this project to interpret the gravity data over the entire state of Utah.