Porosity prediction from gravity inversion with constraints from resistivity data

: This work describes a method to carry out 2-D inversion of gravity data in terms of porosity and matrix density distribution using previous DC resistivity inversion results to constraint the fractional pore-water content in the rocks. The inversion is carried out using a controlled random search (CRS) algorithm for global optimization. The method was tested on synthetic data generated from a model representing a graben, and the results show that it can estimate accurate values of contrast-density and porosity. The method was also applied to gravity and dc experimental data collected in NE Portugal, showing results that agree quite well with the known geological information.

seeks a solution for the inverse problem in terms of the rock-matrix density, porosity, and water-content of the geological formation. The method, which is not a joint inversion procedure, can be considered a cooperative inversion of two coincident geophysical data sets.

Materials and Methods
Bulk density e and fractional porosity , for porous media partially filled with water, are related by the volume weighted density expressed by: where m is the rock matrix density and w the pore-water density, and Sw is the fractional amount of saturation of pore water. The effect of the air in the pores is not considered in this equation.
Archie's law has been currently used for porosity assessment from resistivity (geoelectromagnetic) surveys carried out either on earth surface or in boreholes. For claycleaned rocks and low-frequency or dc methods, the Archie's law is expressed by [5], e = a w Sw -n  -m (2) where a, n and m are empirical constants, w and e are the water and the bulk resistivity, respectively, and Sw is the fractional amount of saturation of pore space. Values between 0.6 to 1.0 for a and, of 2.0 for m and n are commonly used for clay-free sandstones.
The gravity data was inverted using the controlled random search (CRS) algorithm for global optimization (see, e.g., [6]). This algorithm produces a random sample of the objective function hypersurface in the parameter space [6]). All sampled points of the parameter space that produce objective functions smaller than a threshold value, which is a function mainly of the data misfit, are saved and used in a statistical analysis. The mean of all sampled models is then taken as the best parameter estimation.
The CRS algorithm is fully described in the literature (e.g., [6]) and, for this reason, only the relevant points of our implementation will be described. Let g(r) be the Ndimensional vector of gravity observations (Bouguer anomaly) and the theoretical 2-D gravity model responses be represented by f (r, m). The vector r defines the position of the gravity measurements and m is the M-dimensional vector of model parameters, i.e. the porosity  and the matrix density m. To calculate theoretical responses, the region containing the gravitational sources is divided into M cells (Figure 1), consisting of rectangular blocks in the 2-D approach (see [7]). The computed gravity anomaly is then, where A is the gravitational attraction matrix, b is the background density, and (e -b) is the vector containing the model density-contrast values. The element Aij of this matrix A represents the gravitational attraction at the point ri due to cell j having unit density.

Figure 1.
Vertical cross-section of the mesh used in 2-D gravity modelling. ej represents the bulk density of the block j.
The inversion problem consists of finding an estimate of the model parameters m * such that the objective function is defined by where  is a quantity depending of the observational errors, mu and ml representing upper and lower a priori bounds of parameters. These bounds can be estimated from previous geological and/or geophysical results.
The algorithm starts by generating L models, each made up of M cells (or blocks).
The matrix density and the porosity are randomly generated to construct each of these models. For each model, the Sw is calculated from the dc resistivity model and porosity

Ambiguity and sensitivity analysis
The density-contrast () of a block in the model can be expressed by, According to (3), the contribution of each block to the calculated gravity anomaly depends on the respective density-contrast. It is then possible to have several models with different combinations of porosity and matrix density that fit the data equally well. This assertion is the expression of the known problem of nonunique solution in gravity inversion. If n = m, a frequently adopted condition, equation (6) simplifies, This equation shows that, in these conditions, it is possible to decrease the ambiguity limiting the values of m.

Synthetic data example
The described method was firstly applied to gravity and dc-resistivity synthetic data generated from the model shown in Figure 3. Basically, the model simulates a 2-D graben (600 m wide and 437 m depth) filled with sedimentary rocks. The value of the different parameters (geometry, density, porosity, water content, electrical resistivity) used in the data generation are shown in Figure 3 and Table 1. The synthetic gravity data and dipoledipole apparent resistivity pseudosection are also presented in Figure 3.   Table 1. The solid line represents the response of the model shown in Figure 5a). c) Apparent resistivity pseudosection calculated from the same model (b) (dipole distance is 100 m). Contours are in ohm-m. Table 1. Parameterisation of the synthetic model shown in Figure 3.  Contours are in ohm-m.
regional density, to calculate the contrast density. The matrix density was allowed to vary between 2000 and 3100 kg/m 3 . It must be noted that the values of the water resistivity and  Between 1990 and 1993, geophysical, geological, and geochemical research was performed in the Chaves region (NE Portugal) to construct the hydrothermal model associated with hot springs located close to Chaves city [8, 9 and 10]). The hot waters reach temperatures of 78ºC and are known since Roman times. Morphologically the Chaves sedimentary basin (Quaternary), a graben striking in NNE-SSW direction that is bounded by granite (Hercynian) and metamorphic schistose formations (Figure 6a and b), dominates the region. The main fault controlling the structure is the NNE-SSW sinistral fault known as the Chaves-Verin fault. Various geophysical methods, gravity, resistivity, self-potential audiomagnetotellurics and magnetotellurics, have been used to investigate the shallow and deep structures of the Chaves graben. The resistivity method was used to detect and define the geometry of the shallow water circulation within the graben. Figure 6d) shows a dipole-dipole resistivity model corresponding to a profile crossing the graben in the WNW-ESE direction (Figure 6a). This model was obtained from 2-D inversion of the field apparent resistivity data (not shown) using a regularization algorithm [11]). The main characteristic of the dipole-dipole model is the low resistivity zone (<30 ohm-m) depicted in its central part at depths greater than 300 m. As a complement of the regional gravity survey, a gravity profile along the dipole-dipole profile was also carried out ( Figure 6c).
The residual gravity anomaly is characterized by a wide asymmetrical minimum over the graben, corresponding to its sedimentary filling.
The inversion was performed with porosity values constrained to the 0.0 to 0.7 range.
The matrix density was allowed to vary between 2000 and 3100 kg/m 3 . A value of 3.4 ohmm for the pore-water was used [12]) in the calculation of saturation degree. The constant parameters in the Archie's law were a=1.0, n=m=2. In the calculation of the contrast-density a value of 2650 kg/m 3 was assumed as regional density.  The saturation degree distribution estimated during the inversion process is shown in Figure 7b). Except in the uppermost part of the graben the saturation degree is greater than 90%. This is in accordance with the available information that suggests that the hydrothermal aquifer lies at a depth of 300 m. The estimated porosity distribution shown in Figure 7c) is in general lesser than 6%, except in the graben area where values greater than 10% have been estimated. Exceptionally, values greater than 45% appear in some parts of the graben. It is probable that such values are originated by a certain degree of overfitting of the gravity data together with inappropriate values of n and m parameters in Archie law. The most relevant feature in Figure 7c) is the area in the centre of the graben with porosity values of 10-12% that reach higher depths. The location of this relatively high porosity zone shows a strong correlation with the location of the Chaves-Verin fault and reveals the importance of this tectonic-structure in the hydrothermal circulation.
Similar porosity results were obtained from 2D joint inversion of dc and scalar AMT data acquired along the same profile [13]) and applying the Waxman and Smits [14] and the Sen's [15] models. The results suggested that the porosity of the reservoir (in the central part of the graben) is not uniform and might be in the range from 12% to 24%.

Discussion
Prediction of porosity distribution from gravity inversion is feasible if resistivity data (from dc, transient or magnetotelluric surveys) are available to constraint the pore-water saturation. The use of Archie law is limited to free clay areas. However, clay is present in several environments. This can be taken into account using a modified Archie law. This study shows the importance of the a, m and n parameters. It is recommended to obtain experimental values for these parameters. The CRS method is easy to code and is fast enough to be used with advantages in the combined inversion here presented. The results obtained from a particular geological structure (a graben) suggest that the joint inversion of gravity and resistivity data might be a useful strategy in regional hydrological.
However, further work is required to consider more general structures.
Author Contributions: Conceptualization, F.A.M.S. and P.R.; methodology, F.A.M.S. and P.R.; writing-review and editing, F.A.M.S. and P.R. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.