We have developed the open-source software Tesseroids, a set of command-line programs to perform forward modeling of gravitational fields in spherical coordinates. The software is implemented in the C programming language and uses tesseroids (spherical prisms) for the discretization of the subsurface mass distribution. The gravitational fields of tesseroids are calculated numerically using the Gauss-Legendre quadrature (GLQ). We have improved upon an adaptive discretization algorithm to guarantee the accuracy of the GLQ integration. Our implementation of adaptive discretization uses a “stack-based” algorithm instead of recursion to achieve more control over execution errors and corner cases. The algorithm is controlled by a scalar value called the distance-size ratio ([Formula: see text]) that determines the accuracy of the integration as well as the computation time. We have determined optimal values of [Formula: see text] for the gravitational potential, gravitational acceleration, and gravity gradient tensor by comparing the computed tesseroids effects with those of a homogenous spherical shell. The values required for a maximum relative error of 0.1% of the shell effects are [Formula: see text] for the gravitational potential, [Formula: see text] for the gravitational acceleration, and [Formula: see text] for the gravity gradients. Contrary to previous assumptions, our results show that the potential and its first and second derivatives require different values of [Formula: see text] to achieve the same accuracy. These values were incorporated as defaults in the software.

Euler deconvolution has been widely used in automatic aeromagnetic interpretations because it requires no prior knowledge of the source magnetization direction and assumes no particular interpretation model, provided the structural index defining the anomaly falloff rate related to the nature of the magnetic source, is determined in advance. Estimating the correct structural index and electing optimum criteria for selecting candidate solutions are two fundamental requirements for a successful application of this method. We present a new criterion for determining the structural index. This criterion is based on the correlation between the total‐field anomaly and the estimates of an unknown base level. These estimates are obtained for each position of a moving data window along the observed profile and for several tentative values for the structural index. The tentative value for the structural index producing the smallest correlation is the best estimate of the correct structural index. We also propose a new criterion to select the best solutions from a set of previously computed candidate solutions, each one associated with a particular position of the moving data window. A current criterion is to select only those candidates producing a standard deviation for the vertical position of the source smaller than a threshold value. We propose that in addition to this criterion, only those candidates producing the best fit to the known quantities (combinations of anomaly and its gradients) be selected. The proposed modifications to Euler deconvolution can be implemented easily in an automated algorithm for locating the source position. The above results are grounded on a theoretical uniqueness and stability analysis, also presented in this paper, for the joint estimation of the source position, the base level, and the structural index in Euler deconvolution. This analysis also reveals that the vertical position and the structural index of the source cannot be estimated simultaneously because they are linearly dependent; the horizontal position and the structural index, on the other hand, are linearly independent. For a known structural index, estimates of both horizontal and vertical positions are unique and stable regardless of the value of the structural index. If this value is not too small, estimates of the base level for the total field are stable as well. The proposed modifications to Euler deconvolution were tested both on synthetic and real magnetic data. In the case of synthetic data, the proposed criterion always detected the correct structural index and good estimates of the source position were obtained, suggesting the present theoretical analysis may lead to a substantial enhancement in practical applications of Euler deconvolution. In the case of practical data (vertical component anomaly over an iron deposit in the Kursk district, Russia), the estimated structural index (corresponding to a vertical prism) was in accordance with the known geology of the deposit, and the estimates of the depth and horizontal position of the source compared favorably with results reported in the literature.

We present a gravity interpretation method for estimating the relief of an arbitrary interface separating two homogeneous media. The upper medium is discretized into rectangular, juxtaposed prisms whose thicknesses represent the depths to the interface and are the parameters to be estimated from the gravity anomaly. The density contrast of each prism is assumed to be constant and known. To stabilize the inversion, we introduce two kinds of constraints on the depths. The first one requires proximity between the observed and computed depths at isolated points such as those obtained from boreholes (absolute equality constraint). The second one requires that groups of depths approximately follow an established linear relationship among the depths (relative equality constraint). Both kinds of constraints are imposed in the least‐squares sense. We illustrate the method performance by applying it to a synthetic anomaly produced by a simulated basement relief consisting of four narrow and adjacent structural lows. Only two structural lows produced isolated gravity lows. Nonetheless, the whole basement topography was successfully reconstructed with an average error of 4% of the maximum relief amplitude. In this example, the relative constraints established that the thicknesses of adjacent prisms should be as close as possible to each other (overall relief smoothness). As absolute constraints we used point information about the basement depth at five points. In addition to the overall relief smoothness, other relevant geologic information, such as localized relief smoothness, occurring at structural terraces, for example, can be incorporated by assigning different weights to the relative equality constraints. The method was applied to the gravity anomaly of Recôncavo Basin, Brazil, leading to a sharper definition of the basement features relative to previous gravity interpretations of this area.

The Python code that produces the results presented here is available under a BSD 3-clause open-source license at github.com/pinga-lab/paper-moho-inversion-tesseroids Moho model and data: The Moho depth model for South America and the data used to generate it are available under a Creative Commons Attribution 4.0

Extending the compact gravity inversion technique by incorporating a priori information about the maximum compactness of the anomalous sources along several axes provides versatility. Thus, the method may also incorporate information about limits in the axes lengths or greater concentration of mass along one or more directions. The judicious combination of different constraints on the anomalous mass distribution allows the introduction of several kinds of a priori information about the (arbitrary) shape of the sources. This method is particularly applicable to constant, linear density sources such as mineralizations along faults and intruded sills, dikes, and laccoliths in a sedimentary basin. The correct source density must be known with a maximum uncertainty of 40 percent; otherwise, the inversion produces thicker bodies for densities smaller than the true value and vice‐versa. Because of the limitations of the inverse gravity problem, the proposed technique requires an empirical technique to analyze the sensitivity of solutions to uncertainties in the a priori information. The proposed technique is based on a finite number of acceptable solutions, presumably representative of the ambiguity region. By using standard statistical techniques, each parameter is assigned a coefficient measuring its uncertainty. The known hematite and magnetite ore body shape, in the vicinity of Iron Mountain, MO, was reproduced quite well using this inversion technique.

To produce a unique and stable solution in potential‐field interpretation, an inversion method must introduce particular constraints. These constraints will inevitably restrict the type of geological setting where the method may be applied. We present a nonmathematical overview of most stabilizing constraints used in inversion methods. Our purpose is to demonstrate that the inversion results are valuable only if the mathematical stabilizing constraints are translated from the geological setting. We identify five basic types of constraints: (1) lower and upper bounds of parameter estimates; (2) proximity of a parameter estimate to a specified value; (3) proximity between pairs of parameter estimates; (4) concentration of the anomalous source about a geometrical element such as an axis; and (5) source compactness. In practice, if used in isolation, constraints (1), (2), (4), and (5) will not produce geologically meaningful results, regardless of the geological setting of the interpretation area. Constraint (3) may produce geologically meaningful results if the anomalous source has a spatially smooth attribute such as the physical property. We illustrate that constraints 1–4, if used in isolation, cannot delineate the geometry of a simulated sill intruded into a sedimentary basin. The basic constraints may (and should) be combined in inversion to produce geologically meaningful results. We present two examples of effective constraint combination: (1) proximity to a specific value and mass concentration about an axis (used to delineate the thickness variation of a sill intruded in a sedimentary basin) and (2) inequality, proximity of a parameter estimate to a specified value, and proximity between pairs of parameter estimates (used to map a discontinuous basement relief). Usually, the stabilizing constraints are too restrictive to hold at all points of a given geological environment. In this case, we use different constraints in different sub‐areas. Each constraint is based on its compatibility with the actual geology of the subarea.

Geophysics is the science of using physical observations of the Earth to infer its inner structure. Generally, this is done with a variety of numerical modeling techniques and inverse problems. The development of new algorithms usually involves copy and pasting of code, which leads to errors and poor code reuse. Fatiando a Terra is a Python library that aims to automate common tasks and unify the modeling pipeline inside of the Python language. This allows users to replace the traditional shell scripting with more versatile and powerful Python scripting. The library can also be used as an API for developing stand-alone programs. Algorithms implemented in Fatiando a Terra can be combined to build upon existing functionality. This flexibility facilitates prototyping of new algorithms and quickly building interactive teaching exercises. In the future, we plan to continuously implement sample problems to help teach geophysics as well as classic and state-of-the-art algorithms.

We present a method to estimate the basement relief as well as the density contrast at the surface and the hyperbolic decaying factor of the density contrast with depth, assuming that the gravity anomaly and the depth to the basement at a few points are known. In both cases, the interpretation model is a set of vertical rectangular 2D prisms whose thicknesses are parameters to be estimated and that represent the depth to the interface separating sediments and basement. The solutions to both problems are stable because of the incorporation of additional prior information about the smoothness of the estimated relief and the depth to the basement at a few locations, presumably provided by boreholes. The method was tested with synthetic gravity anomalies produced by simulated sedimentary basins with smooth relief, providing not only well-resolved estimated relief, but also good estimates for the density contrasts at the surface and for the decaying factors of the density contrast with depth. The method was applied to the Bouguer anomaly from Recôncavo Basin, estimating the surface density contrast and the decaying factor of the density contrast with depth as [Formula: see text] and [Formula: see text], respectively.