Iterative 2D Gravity Modelling for Basement Depth Estimation

In gravity modelling, iterative refinement of 2D basement geometry is a classical problem where the density contrast has to be determined in advance. We present a relatively simple method to estimate 2D basement depth along with its density contrast. Known basement depths at several points and their associated gravity anomalies determine the gravity - depth relationship that is used to obtain the first estimate of the overall basement depth. The initial model response is calculated with a unitary density. The linear regression between calculated model response and actual gravity data results in the density contrast estimation. The misfit or residuals and the gravity - depth relationship are used to update the model iteratively. The proposed technique is relatively stable and results in a recovered model quite similar to the synthetic model. We also explore the possibility to include estimation of the constant the regional anomaly.


Introduction
The well-known method proposed by Bott in 60s allows 2D basement modeling by iterative algorithm [1]. The Bouguer correction formula is be used to estimate the basement depth from the gravity anomaly at each point, i.e. Bouguer plate or 1D approximation. The horizontal position and depth points constitute a polygon defining the cross-section of an initial 2D model. Talwani's method [2] is used to calculate the gravity response of the initial model. The misfit between calculated and observed gravity anomaly determines the correction factor for basement depth, also by using the Bouguer plate approximation. The process is done iteratively until convergence. As an extension of the iterative 2D modeling, similar method has also been proposed for 3D basement geometry [3]. The initial model and its iterative refinement are similar to the method for 2D modelling, while the 3D gravity forward modelling is calculated by using vertical prism as the building block of the 3D basement model [2].
The classical Bott's method needs a relatively accurate estimate of the density contrast which is critical in gravity modelling. Slightly different density contrast may lead to significant differences in model geometry. In addition, this approach is inherently unstable. However, a regularization or damping can be included to stabilize the iteration. The latter is done by formulating the algorithm in the formal Gauss-Newton method for non-linear inverse problems [4]. Recently, other type of stabilization was implemented for 2D and 3D iterative gravity modelling, i.e. using maximum differences of a factor in reducing the residual [5,6]. In this paper, we adopt the algorithm proposed by Florio [7,8] originally for 3D basement modeling for 2D gravity modeling. We also evaluate the validity of the method in the presence of noise. Furthermore, we propose to include the estimation of the constant regional gravity anomaly in the iterative process. In all cases, we used synthetic gravity data associated with a relatively realistic basement model.

Method
The Bouguer plate approximation is the simplest form of linear relationship between gravity and its anomalous source's geometry, i.e. thickness. In basement modelling, the thickness of the sediment is in fact the depth of the basement topography. A more realistic gravity -depth linear relationship can be deduced from known basement depth at several points with their actual associated gravity anomalies [7][8][9]. Such gravity -depth relationship allows the first estimate of the basement depth without assuming the density contrast value in advance. This avoids one of the difficulties in Bott's and other similar methods.
We consider gravity data distributed along a profile line perpendicular to the strike of a 2D model. The cross-section of the 2D model is represented by a polygon defined by its corners (xk, zk), k = 1, 2, …, N with N is the number of vertices. Overall, the method can be summarized in the following algorithm: 1. A linear regression is done between known basement depths at several points (d, in km) with their associated observed gravity anomalies (g, in mGal). With the observed gravity as horizontal axis (x) and the known basement depths as vertical axis (y), we have, where a is the gradient or slope of the regression line and b is the constant. Equation (1) is used to estimate the basement depth, i.e. the initial model, from the gravity anomaly at every point.

2.
The theoretical gravity response of the initial model is calculated by using Talwani's method [2] with a unit density contrast, i.e. 1 gr/cm 3 or 1000 kg/m 3 in SI.
3. A linear regression is performed between observed gravity (y-axis) and calculated gravity with a unit density contrast (x-axis). The gradient of the straight-line curve gives the estimated density contrast that will be used for all subsequent forward modeling calculations.
4. The forward modeling is done for the initial model with the estimated density contrast to obtain the difference between observed and calculated gravity (E) at each point.
5. The gravity difference (E) associated with the initial model is used to obtain the depth modification at every point by using the following formula, where a is the gradient or slope in Equation (1). 6. The steps number (4) and number (5) are repeated until a convergence is achieved or other stopping criteria are reached.
The criteria for convergence are either of the following: the total misfit falls below the target misfit (e1), or the total misfit difference at two consecutive iterations becomes very small, i.e. less than a threshold value (e2). The maximum number of iterations is set to avoid infinite loop if the convergence cannot be achieved.

Application to Synthetic Data
The iterative 2D gravity modeling algorithm was applied to invert synthetic gravity data associated with a relatively complex basement model. The synthetic model and data without noise (density contrast -0.5 g/cm 3 ) are shown in Figure 1. Since we deal with synthetic data, the linear regression was performed by using depths from the vertices of the synthetic model and calculated gravity anomalies associated with them (Equation 1). The term "observed" gravity will be used for the synthetic gravity data. The results of the regression are: a = -0.0890 and b = -1.8764 ( Figure 2) which are then used for gravity to depth conversion to obtain the initial model. In Figure 2, the slope of the regression line does not represent a negative value since the depth is set to positive downwards. The regression between the initial model response with a unit density (calculated) and actual gravity anomaly (observed) has led to estimated density contrast of -0.4651 g/cm 3 with a constant term of -1.2744 (Figure 3). Note that the regression line reflects the correlation between the negative values of both observed and calculated gravity. Figure 4 shows the initial basement model and its calculated gravity response. The gravity to depth conversion results in a very approximative geometry of the basement model. However, the model response is quite similar to the observed anomaly except for the details in the spatial high frequency variations.    Iterative refinements of the initial model led to the model presented in Figure 5 with overall or average misfit of 0.044 mGal after 31 iterations. Although the synthetic data do not contain noise, the iterations cannot be done further to obtain smaller misfit. In this case, e1 was set to 0.05 mGal. Similar exercise was also performed with synthetic data containing noise of 0.5 mGal with Gaussian distribution. The final results are presented in Figure 6. The misfit of 0.294 mGal was obtained after 15 iterations. The final misfit around or slightly smaller than the noise content in the data is in general considered adequate, i.e. e1 was set to 0.3 mGal. Inverse models from synthetic data with and without noise recovered the synthetic model satisfactorily. Higher noise in the data (e.g. 1 mGal or more) would result in more oscillations of the final model. In such case, a smoothing process, either during or after inversion would be needed to obtain more meaningful results.

Estimation of the Regional Anomaly
In the case of a linear gravity -depth relationship, a constant addition to the gravity anomaly would lead to similar feature of the basement model. Therefore, we investigate the possibility to model gravity data containing the simplest, i.e. constant, regional anomaly. For this preliminary attempt, we consider only the synthetic data without noise to simplify the problem. We add to the synthetic gravity data a +6 mGal constant regional anomaly and proceed with the same iterative algorithm as before.
The regression between depths and associated gravity anomalies gives a = -0.0890 and b = -1.3423 that are used to obtain the initial model. The gradient is the same as before, while the constant term differs due to the presence of the additional regional anomaly. The regression between the initial model response calculated with a unit density and the observed gravity anomaly results in the same density contrast of -0.4651 g/cm 3 with a different constant term of 4.7226. The initial model and its gravity response calculated with the estimated density is shown in Figure 7. With the same initial model and density contrast as before, the discrepancy between calculated and observed gravity is certainly due to the regional anomaly contained in the gravity data.
The constant factor of the second regression (to obtain the density contrast) is proportional to the added regional gravity anomaly. Thus, this constant can be considered as a representative of the regional anomaly, although its value is not exactly the same. Furthermore, an approximation to the regional anomaly led to "correct" basement model. In this case, we add the estimated regional anomaly to the IOP Publishing doi:10.1088/1755-1315/1031/1/012014 5 calculated anomaly associated with the basement model estimates. In such case, the observed gravity anomaly shown in Figure 7 is not changed in the iteration process. The result is presented in Figure 8, with misfit of 0.148 mGal after 50 iterations. The basement model recovers well the synthetic model and it does not differ significantly with the basement model obtained from gravity data without the regional anomaly added (see Figure 5).

Figure 7.
Synthetic data with regional anomaly added and calculated gravity data of the initial model with estimated density contrast (top), synthetic and initial models (bottom).

Figure 8.
Modeling of synthetic data with regional anomaly, observed and calculated data of the final basement model (top), synthetic and inverse models (bottom).

Conclusion
We have presented a relatively simple iterative 2D gravity modelling similar to the iterative method that relies on the gravity -depth or thickness relationship. It does not involve complex matrix multiplications nor inversions. The proposed technique is stable for synthetic data without and with 0.5 mGal Gaussian noise. In general, the final models agree well with the synthetic model. We also improve the algorithm to include the estimation of the constant component of the regional anomaly with satisfactory results. The algorithm can be extended to involve more complex density contrast that varies with depth. In such case, the regression would not be linear but follows the density variation with depth, e.g. exponential.