Salt structures from inversion of residual gravity anomalies: application in Santos Basin, Brazil

The Santos Basin, with an area of about 350.000 km2, is the largest salt basin of the South Atlantic, and due to its high economic hydrocarbon potential, it is a recurrent theme in scientific studies. The salt structures over the region present great importance for hydrocarbon accumulation and the geological/geophysical studies are performed from seismic reflection data, which requires time and efforts for acquisition and data processing. We identify salt structures using a new workflow based on inversion of residual gravity anomalies, where we use the Moho and basement depths obtained from gravity inversion, followed by the calculation of the gravity residual anomaly, assumed to be representative of the salt structures. This workflow is tested for a geological profile in the Santos Basin, and the results are evaluated along a 2D seismic section tied to well markers. The geometry of the stratified salt obtained from gravity inversion correlates with the seismic interpretation, with the advantage of estimating the entire salt package, including halite and stratified salt. With only seismic data, sometimes the stratified salt can be misinterpreted as sediments. The procedure can be applied to identify salt in sedimentary basins where seismic data is unavailable or of low quality.


INTRODUCTION
The gravity modeling of salt structures in sedimentary basins has been subject of several scientific studies. Santos Basin, located in the southeastern Brazilian margin, is the salt basin in the South Atlantic (Mohriak et al. 2012). The geological knowledge related to these salt structures presents high interest for the petroleum industry in the region, due to their properties in sealing the reservoirs, guiding the migration of hydrocarbons and also facilitating the deformation responsible for many structural and stratigraphic hydrocarbon traps (Demercian et al. 1993).
Evaporites are characterized by physical properties, such as density and velocity, that differentiate them from other sedimentary rocks; hence, geophysical methods can be used to help the geological characterization of salt structures. In the seismic method, chrono-stratigraphic reflectors are identified representing mainly stratigraphic layer boundaries and unconformities reflecting changes in acoustic impedance. In this sense, salt structures composed mainly by halite show high seismic velocity (around 4,000-5,500 m/s) and low density (about 2.17 kg/m 3 ) when compared with surrounding rocks.
Furthermore, seismic interpretation can be easily accomplished where the structures are basically salt diapirs. In the gravimetric study, the changes in the gravity field occur due to density changes and, therefore, a gravity anomaly can be associated to buried features such as diapirs (Dobrin 1976, Blakely 1995. This method can be useful in regions of exploratory potential with little available geological information, providing an initial geological framework with important information to guide future studies involving other methods, such as seismic acquisition. Most of the available information in literature aimed to study salt structures in sedimentary basins used the seismic method as the main one to define their geometry. While the reflection seismic method presents a subsurface image that shows geological structures and horizons for some depths in detail, the gravimetric method presents the contribution of all geological bodies from all depths. The combined analysis of seismic and gravimetric data may help geoscientists to identify geological features, such as diapiric structures, in a more reliable way than using only one of the methods (Dobrin & Savit 1988).
In the Santos Basin, only few studies of crustal modeling that use seismic and gravity data are available. Mio (2005), and Lima and Mohriak (2013) modeled crustal structures in order to create geological models for the basin, focusing primarily on the interpretation of the seismic data, adjusted to the gravity field. When interpreting 2D seismic reflection data for areas with the presence of salt bodies, the basement is not so clearly identifiable like the top and base of the salt layer, due to the loss of data resolution below salt packages. The Moho interface, even deeper, scarcely appears in 2D seismic reflection data. In such cases, obtaining this interface by a distinct method may be a good solution.
This paper describes a technique to obtain the depth of the basement and the Moho interfaces using only gravity data, acquired for this purpose or available from global models. After the calculation of these surfaces, the forward modeling provides a residual gravity anomaly that can be used to determine a good approximation of the geometry of the salt deposits, allowing even the detection of layered salt structures, in an independent methodology that does not require seismic data. This paper uses the available seismic data in the Santos Basin, Brazil, to check the results obtained by this method and to ensure that the procedure is feasible to provide a good estimate of these structures with no significant cost, allowing a better planning of the seismic data acquisition and a refined geophysical interpretation.

Geological outline of the Santos Basin
Located in the southeastern region of the Brazilian continental margin, between latitudes 23º00'S and 28º00'S, the Santos Basin extends to water depths of 3,000 m ( Fig. 1). It has as area of about 350,000 km 2 , including the coasts of Rio de Janeiro, São Paulo, Paraná and Santa Catarina. The basin is bounded to the north by Cabo Frio High and to the south by the Florianópolis lineament (Moreira et al. 2007).
The formation of the basin is related to the breakup of the Supercontinent Gondwana, which resulted in the separation of South America and Africa. A deposition of thick evaporite packages, composed mainly of halite (82%), anhydrite (11%), and bittern salts (7%), ( Jackson et al. 2015, Rodriguez et al. 2018, occurred during the Upper Aptian, due to an arid environment and to periodic marine transgressions in a shallow gulf extending from the Santos Basin to the Sergipe-Alagoas Basin (Gamboa et al. 2008).
At the Upper Cretaceous (Santonian-Campanian), an uplift in the continental area was registered due to the reactivation of the previous formed basement faults, resulting in a remarkable erosion and progradation of clastic wedges into the basin, which was subject to a subsidence process at the same period (Macedo 1987, 1989, Almeida & Carneiro 1998. This siliciclastic progradation was partially responsible for deforming and moving the evaporite packages through the basin depocenter in the region known as São Paulo Plateau. In addition, an oceanward tilt occurred while the basin evolved into one dominated by thermal subsidence, which resulted in the translation of the salt basinward (Modica & Brush 2004, Guerra & Underhill 2012. The lower halite layer started its movement immediately after its deposition, furthering the formation of mini basins, where variable thickness of stratified salt was deposited (Gamboa et al. 2008, Davidson et al. 2012. Both halite and anhydrite show the same seismic velocity with depth and thus, the seismic amplitude contrast between these salt layers and adjacent sediments depend only on their acoustic properties. The halite Figure 1. Location of Santos Basin. The red line shows the profile used for forward and inverse modeling. The black dotted line represents the continent-ocean boundary according to Cainelli and Mohriak (1999).
shows high seismic velocity (about 4,500-5,500 m/s), but low density (around 2.17 g/cm 3 ). The anhydrite is also represented by high velocity (about 4,000-5,500 m/s), and high density value (approximately 2.98 g/cm 3 ) (Bassiouni 1994, Mavko 2005. Due to these differences in density, the gravity modeling integrated with seismic data presents a good solution to identify types of evaporates within the interbedded salt.

METHODOLOGY
The diagram displayed in Figure 2 summarizes the main steps of the proposed methodology, discussed next.
Step 1 The first step consists of gravity inversion to model Moho and basement depth variations. In this paper, we use the Moho and the basement depths from Constantino et al. (2016), where the method is comprehensively explained. The main points are briefly described in the next sections.

Moho depth
To estimate the Moho depth from the gravity inversion, the free-air anomaly ( Fig. 3) is corrected for the bathymetry (Bouguer correction) and for the sediment thickness variations. The resulting gravity field is then inverted by applying an iterative constrained inverse modeling (Braitenberg et al. 1997).
The gravity effect due to a sedimentary package can be calculated by the Parker Algorithm (Parker 1973) with a constant density contrast along the discontinuity; however, to more realistic results, the sediment compaction with depth should be considered. In this study a compaction model described by Sclater and Christie (1980), based on an exponential reduction of porosity with depth, is used. The density in dependence of the depth below the seafloor is calculated with Equation 1: where: = initial porosity of sediments; d = decay parameter.
The values of density, porosity and decay parameters were calibrated from well data information. The gravity effect of sediments is subtracted from the Bouguer anomaly and the resulting filed is inverted in order to obtain Moho undulations.
For this, only the long wavelength part of the observed gravity field is used. The cutoff wavelength is estimated from the decay of the amplitude spectrum of the gravity field (Russo & Speed 1994). The method requires some parameters, including the reference depth of the interface, and also the density contrast across the interface (between crust and mantle). To constrain these parameters, the reference depth (d) and the density contrast are varied for the different combinations of parameters and the root mean square (RMS) difference with punctual values obtained from another source is used. The best fit between gravity-derived and punctual values is the one considered as the best values for reference depth and density contrast.

Basement depth
To find the basement topography, the gravity effect of sediments and Moho are subtracted from observed gravity anomaly (g obs ) . The gravity effect of Moho g( moho ) is calculated using a constant density contrast along the Moho interface by applying the Parker algorithm (Parker 1973). The gravity effect of the sediments, already calculated for inversion of the Moho, should be estimated again (g seed ). When inverting for the ocean bottom from the gravity signal, we must use the density contrast between crust and water. In order to be able to also use this same density difference for the inversion of the basement, the sediments-gravity effect must be calculated against the reference density of water. Sediments can be regarded as a surplus mass relative to water (Braitenberg et al. 2006).
The obtained residual field (Equation 2) is inverted by applying the iterative constrained inverse modeling. Density contrast between basement and water is used and the reference is set to level zero. The procedure results in the basement topography (Equation 2).
Step 2 In this step, a two-dimensional density model is built, and the gravity response of the model is calculated based on Talwani et al. (1959) and compared with the observed anomaly. The goal of this step is to model the crust using the interfaces with clearly defined and known contrasts, that is, the Moho and the basement. The density model is created with the interfaces obtained by gravity inversion (step 1), resulting in a modeled crust with the main discontinuities: Moho, basaltic basement and seabed. Due to the intentional absence of the salt package in this model, residual anomaly (calculated -observed) is obtained. The residual anomaly is inverted in the last step to highlight the salt bodies.
Step 3 In this step, the approach developed by Silva and Barbosa (2006) for gravity inversion is applied. The method locates and delimits one or more anomalous bodies inside a known area where gravimetric data is available. For this method, the main assumptions are: • the anomalous body and the embedding rock must have an infinite dimension along the axis perpendicular to the profile; • the anomalous body must be compact (without holes in its interior) and should contain the smallest possible volume around the geometric elements (points and lines) defined by the user; • the anomalous body must present a uniform density distribution, and the contrast with the embedding rock must be known.
Considering all these elements and assuming that salt bodies are consistent with the assumptions required by the method, the residual anomaly obtained in step 2 can be inverted in order to obtain the total thickness of a possible salt body for a selected profile in the study area.
The interpretation model is a grid composed of 2D juxtaposed prisms, of which density contrasts are to be determined through inversion. Outlines of the gravity sources are specified in terms of geometric elements (line segments and points) and the density contrast, associated with the geometric elements, defines each gravity source framework. The density-contrast distribution that fits the observed anomaly is then estimated and the gravity sources closest to the specified geometric elements are represented. Each source is admitted as a homogeneous body with a known density Figure 3. Free-air anomaly for the study area. The free-air anomaly is taken from Molina (2009), a satellite-derived data with 2 × 2 arc min original resolution, resampled on a regular grid spacing with of 6 km. For reference, bathymetric isolines for water depths 1,000, 2,000 m and 3,000 m are shown. The red line shows the profile used for forward and inverse modeling.
contrast. The method from Silva and Barbosa (2006) is briefly described next.
The gravity anomaly ( ) In which: g = an N x 1 vector whose ith element is A = an N ×M matrix whose generic element is given by Aij.
The method searches for a solution to this linear system that satisfies the gravity anomaly, presenting most of its mass excess (or deficiency) around the geometric elements defined interactively by the user. This is achieved by the method of Guillen and Menichetti (1984), following some steps described below.
First, a standard minimum-norm solution is obtained by Equation 5: In which: μ = a nonnegative scalar; T = a transposition operator; I = the identity matrix.
The larger the value of μ, the smaller the Euclidean norm of . Thus, parameter is updated iteratively by Equation 6: In which (Equation 7): W (k) = a diagonal matrix whose nonzero elements are given by With W (k) a nonlinear matrix, there is a nonlinear inverse problem that must be solved iteratively. The iteration is initialized by setting wjj = 1wjj = 1 and it stops when (Equation 10): In which: τ = a positive scalar controlling the degree of homogeneity assumed for the estimated gravity sources; v j = the target contrast.
During the inversion process, three parameters must be defined: Parameter μ favors solutions in which the elementary prisms lying in the neighborhood of the geometric elements receive the largest density-contrast estimates in absolute values. It also controls the stability of the solutions; parameter f is the freezing parameter. A large value assigned to f tends to maintain the frozen density-contrast estimates along iterations; and tolerance (τ), is a positive scalar controlling the degree of the density contrast homogeneity, avoiding unrealistic density-contrasts for the estimated gravity sources. The larger the value of τ, the greater its homogeneity. A typical value is 0.01, corresponding to a source with density-contrast variation less than or equal to 1% of the target density contrast.

RESULTS
In order to perform the forward model, a seismic profile from the database of the Brazilian National Agency of Petroleum, Natural Gas and Biofuels (ANP) is selected and used to constrain the results (Fig. 1).
The model is built based on the uppermost density discontinuities of the crust, which are: • ocean floor: transition from water to sediments; • basement: transition from oceanic sediments to the basaltic rock; • Moho: transition from the crust and to the mantle.
The bathymetry (ocean bottom) was obtained from the General Bathymetric Chart of the Oceans (GEBCO) and the data estimated from gravity inversion were used for Moho and basement depths (Figs. 4 and 5). The density values assigned to each layer are average values, 2.5 g/cm 3 for the sediments, 2.67 g/cm 3 for the crust and 3.3 g/cm 3 for the mantle.
Usually, in order to apply a forward modeling based on gravity field data, it is necessary to create a model with the maximum geological information available for a specific area, fitting observed and calculated gravity anomalies (Fig. 6A). However, the objective of this paper was to create a simple model, disregarding a priori information about the presence of the salt package, and then to calculate the inversion of the residual gravity anomaly to obtain the dimension of this package.
The residual gravity anomaly (Fig. 6B) is shifted in order to make it negative all along the profile. The shift is used for two main reasons: the inversion procedure used in this paper is not Figure 4. Moho depth obtained from gravity inversion. For reference, bathymetric isolines for water depths 1,000, 2,000 and 3,000 m are shown. Data provided by Constantino et al. (2016). The red line shows the profile used for forward and inverse modeling. suitable for both positive and negative anomalies. Beyond that, the negative anomaly referring to the body of interest (salt) is a priori information, so it is reasonable to adopt this assumption. The geometric element is placed at a 5 km of depth (average depth for the salt layer) and extending along the entire profile. The anomaly is inverted using an interpretation model consisting of a grid of 126 × 64 cells with dimensions of 1 and 0.125 km in the x-and z-directions, respectively, and setting μ = 0.25, f = 500, and τ = 0.01, with a maximum of 30 iterations.
The choice of these parameters is done looking for the stability and unicity of the solution. An unstable solution is one that presents high variability of parameters due to small variations in data. Lack of uniqueness is characterized by the existence of several sets of parameters that produce the same predicted data.
The parameter μ is commonly adopted for inverse problems in geophysics to overcome problems of instability and no unicity. It controls the stability of the solution and the choice of its value is of great importance. High values favor the mismatch of the data and very low values of μ make the solution unstable. An optimum value for μ is the smallest positive value still producing stable solutions .
The stability of a solution is verified by contaminating the data with different sequences of pseudorandom Gaussian noise with mean 0 and standard deviation of 0.1 mGal. The contaminated data for each series of pseudorandom noise values was inverted to values between 0.05 and 1. Standard deviation was calculated for each value between the parameters obtained for each cell. The solution with a standard deviation less than 0.1 mGal was considered stable.
The tolerance (τ) is maintained at 0.01, which corresponds to a variation of 1% or less between the solution and the target density contrast. The frozen parameter (f) is established as 500. This parameter is the value assigned to the weight to "freeze" the solution over a certain interval, and to do this, it must be high enough. If, during the inversion, the last iterations present approximately constant values and maximum density contrast values greater than the target density contrast, it means that the value of f is not high enough to "freeze" the solution and, according to Barbosa and Silva (2006), should be increased in multiples of 5 or 10. During the inversion, the maximum density contrast value converged to the target density contrast for f = 500.
The inversion of the residual gravity field is shown in Figure 7. Cells filled with dark red show density contrast (Δρ) closest to target contrast (-0.42 g/cm 3 ).

DISCUSSION
From the inversion of the residual gravity anomaly field applied over a profile in Santos Basin, a salt body located at the average depth of 5 km with density contrast of -0.42 g/cm 3 is suggested (see Fig. 7). To check the reliability of the inversion results, the salt body is compared with the interpretation of a seismic line located along the chosen profile (Fig. 8).
The seismic horizons were converted to depth using the methodology described by Vincentelli (2008), with velocities  of 1,500, 2,500 and 4,000 m/s for seabed, sediments and salt, respectively. All of these values were obtained from interval velocities measured on sonic logs. Some uncertainties are pointed out, representing potential packages of stratified salt (ES), which, according to Gamboa et al. (2008), consists of interbedded salt layers such as halite, anhydrite and complex salts deposited in shallow water within mini basins. The question marks (?) correspond to possible depositions of ES.
The comparison of the salt body from the inversion procedure and the salt interpreted in the seismic line is shown in Figure 9, assuming that the suggested uncertainties are actually ES. To check the consistency of this affirmation, a complete density model of the crust is formulated expecting a good adjustment of the observed and fitted anomaly. For this model, all the important density contrast layers from ocean bottom to mantle are considered: bathymetry from GEBCO and seismic data, marking the transition from water to sandy, muddy or solid material, generally sediments; top of the salt package interpreted from seismic, marking the transition from upper sediments to halite; base of the salt (also interpreted from seismic), marking the transition from halite to lower sediments; basement and Moho from Constantino et al. (2016), marking the transition from lower sediments to basaltic rock and between crust and mantle, respectively (Fig. 10).
Default density values of 3.3 and 2.67 g/cm 3 are used for the mantle and the crust, respectively. For the sedimentary layer, considering that it is composed mainly of sandstones, shales and carbonates (Mio 2005) with average densities of 2.37, 2.6 and 2.65 g/cm 3 (Dobrin & Savit 1988), Figure 9. Comparison between the inverse model, obtained with gravimetric data only (blocks) and the outline of the salt obtained from seismic interpretation (pink line), assuming that the suggested uncertainties are actually ES. The pink contour is the same presented in figure 10 (base and top of the entire salt layer interpreted from seismic, with the ES adjusted from gravity modeling). a density of 2.58 g/cm 3 is used for the upper layer. Considering the compaction with depth and consequent density increase, a slightly higher density of 2.6 g/cm 3 is used for the lower layer, which contains a similar lithology. For salt package and ES, the average density of halite 2.17 g/cm 3 is chosen. Even with the presence of denser salts in the ES, this density can be justified by two arguments: first, the good adjustment of the observed and calculated anomalies, with an RMS of only 2.6 mGal, and second, because the gravimetric method cannot discern this type of horizontal stratification when present inside a larger sedimentary package, therefore, an average density value is reasonable.

CONCLUSION
We draw the following conclusions from this study: • The proposed methodology does not require seismic data and estimates salt structures from inversion of residual gravity anomalies. Although the method does not discern between halite-rich intervals and stratified salt, it gives an estimative of the entire salt package present in the basin. When using only seismic data to interpret salt structures, the stratified salt can be interpreted as sediments other than salt. A good solution to resolve this problem would be analyzing the seismic data together with the salt estimative obtained from the gravity inversion proposed in this paper.
1. The procedure can be fast applied and with very low costs to help identifying the location of salt packages before designing for a seismic acquisition; 2. For future studies, the following steps can be applied to estimate the salt structures in Santos Basin: 1. obtain Moho and basement interfaces from gravity inversion and use them to calculate a gravity forward model containing three layers: sediments, crust and mantle; 2. estimate the residual anomaly of the forward model; 3. invert the gravity residual anomaly using a flat geometry (line) located in the average depth of the package as the target geometric element and a contrast density of -0.42 g/cm 3 , eventually working with a shift in the residual anomaly so that it becomes entirely negative; 3. The method can be applied to other salt basins. In this case, a prior study should be performed to adjust mainly the parameters of density contrast and average depth of the salt layer.