Gravitational field calculation in spherical coordinates using variable densities in depth

We present a new methodology to compute the gravitational fields generated by tesseroids (spherical prisms) whose density varies with depth according to an arbitrary continuous function. It approximates the gravitational fields through the Gauss-Legendre Quadrature along with two discretization algorithms that automatically control its accuracy by adaptively dividing the tesseroid into smaller ones. The first one is a preexisting two dimensional adaptive discretization algorithm that reduces the errors due to the distance between the tesseroid and the computation point. The second is a new density-based discretization algorithm that decreases the errors introduced by the variation of the density function with depth. The amount of divisions made by each algorithm is indirectly controlled by two parameters: the distance-size ratio and the delta ratio. We have obtained analytical solutions for a spherical shell with radially variable density and compared them to the results of the numerical model for linear, exponential, and sinusoidal density functions. The heavily oscillating density functions are intended only to test the algorithm to its limits and not to emulate a real world case. These comparisons allowed us to obtain optimal values for the distance-size and delta ratios that yield an accuracy of 0.1% of the analytical solutions. The resulting optimal values of distance-size ratio for the gravitational potential and its gradient are 1 and 2.5, respectively. The density-based discretization algorithm produces no discretizations in the linear density case, but a delta ratio of 0.1 is needed for the exponential and most sinusoidal density functions. These values can be extrapolated to cover most common use cases, which are simpler than oscillating density profiles. However, the distance-size and delta ratios can be configured by the user to increase the accuracy of the results at the expense of computational speed. Lastly, we apply this new methodology to model the Neuquén Basin, a foreland basin in Argentina with a maximum depth of over 5000 m, using an exponential density function.


Introduction
The lithosphere's density variation with depth has been studied for close to a century. Over this time period, several density-depth relations have been proposed for different rock types (e.g. Maxant, 1980;Rao, 1986;Rao et al., 1993Rao et al., , 1994. Furthermore, depth-variable densities have been used in forward and inverse gravity modelling, mostly applied to sedimentary basins (Cordell, 1973;Rao, 1986;Cowie & Karner, 1990;Rao et al., 1993Rao et al., , 1994Zhang et al., 2001;Welford et al., 2010). These forward gravity models have been developed for two or three dimensional bodies in Cartesian coordinates, which limits the applications to local scales. The advent of satellite gravimetry has provided gravity field measurements with global coverage, enabling modelling and interpretation on regional and global scales. Hence, designing forward modelling methods that reproduce the gravity anomalies for such scales is of high importance.
To take into account the curvature of the Earth, many global forward modelling methods are defined in geocentric spherical coordinates. A common approach is to discretize the Earth Figure 1: A tesseroid (spherical prism) in a geocentric spherical coordinate system, with a computation point P and its local north oriented Cartesian coordinate system. After Uieda (2015).
into tesseroids (Anderson, 1976), i.e. spherical prisms, which are defined by pairs of latitude, longitude, and radial boundaries (Fig. 1). The gravitational fields generated by a tesseroid on any external point are given by volume integrals that must be numerically approximated. The literature offers two main approaches: one involves Taylor series expansion (Heck & Seitz, 2007;Grombein et al., 2013) while the other makes use of Gauss-Legendre Quadrature (GLQ) (Asgharzadeh et al., 2007;Wild-Pfeiffer, 2008;Li et al., 2011;Uieda et al., 2016;Lin & Denker, 2019). The Taylor series expansion is not well suited to develop an algorithm for a density varying with depth according to an arbitrary continuous function. Different series expansion terms would have to obtained for each density function desired. Conversely, an arbitrary continuous density function can be included in the GLQ without any change to the integration method. Thus, we will focus only on GLQ-based methods henceforth.
The main challenge of the GLQ integration is the loss of accuracy that occurs when the computation point approaches the tesseroid (Ku, 1977). Uieda et al. (2016) built on the three dimensional adaptive discretization algorithm of Li et al. (2011) to automatically obtain integration results with 0.1% accuracy. The algorithm recursively divides the tesseroid into smaller ones when a threshold is exceeded, namely when the normalised distance to the computation point is greater than a "distance-size ratio" parameter (D). Uieda et al. (2016) have also obtained standard values of D for the gravitational potential, its gradient and Marussi tensor components by comparing the numerical integration results with the fields generated by a spherical shell.
Two recent publications present alternative approaches for calculating the gravitational fields of homogeneous tesseroids and incorporate methodologies for tesseroids with variable densities in depth. Fukushima (2018) analytically integrates the volumetric integral for the gravitational potential in the radial direction, obtaining a surface integral, which is then numerically solved by conditionally splitting the tesseroid and applying the double exponential quadrature rule. The gradient of the potential and the Marussi tensor components are computed by finite differences. Fukushima (2018) also generalized their method to tesseroids with a radial polynomial density function of arbitrary degree. Lin & Denker (2019) compared the different integration and discretization methodologies for homogeneous tesseroids. From this analysis they developed a combined method: for computation points near the tesseroid, they use a GLQ integration with an adaptive discretization based on Uieda et al. (2016) but only applied to the horizontal dimensions. If the computation point is farther than a certain truncation distance, a second order Taylor series approximation is applied instead along with the regular subdivision developed by Grombein et al. (2013). Lin & Denker (2019) also introduced a variation of their combined method to compute the gravitational fields generated by tesseroids with a linearly varying density in the radial dimension.
Both the Lin & Denker (2019) and the Fukushima (2018) studies limit the radial density variation to polynomial functions. While most continuous and smooth functions can be approximated by piecewise linear functions, the choice of a discretization interval is neither straight forward nor automatic for the general case. Although many piecewise linear approximation algorithms that automate this process exist (Ketkov, 1969;Vandewalle, 1975;Imamoto & Tang, 2008;Ahmadi et al., 2013), they either require a fixed number of discretization intervals or are meant to be used only on convex functions. Using such algorithms would limit the range of density functions that can be assigned to the tesseroids and would not guarantee a fully automatic process. Furthermore, it is well known that the use of high-degree polynomials to approximate a highly variable function produces unstable results when extrapolating beyond the data domain. These shortcomings could make piecewise linear or high-degree polynomial density functions unwieldy for non-linear gravity inversions (e.g. Uieda & Barbosa, 2017) if the density function in depth is highly variable.
We present a new algorithm for computing the gravitational potential and its gradient generated by a tesseroid with an arbitrary continuous density function on an external point. It is based on the three dimensional GLQ integration, a two dimensional version of the adaptive discretization of Uieda et al. (2016) (following Lin & Denker (2019), and a new densitybased radial discretization algorithm. To ensure the accuracy of the numerical approximation, we empirically determine optimal values for the controlling parameters by comparing the numerical results with analytical solutions for spherical shells. Finally, we apply the methodology to model the Neuquén basin, Argentina, using tesseroids with linear and exponentially increasing density with depth.

Methodology
Consider a tesseroid in a geocentric spherical coordinate system defined by pairs of geocentric latitudinal (φ 1 , φ 2 ), longitudinal (λ 1 , λ 2 ), and radial (r 1 , r 2 ) boundaries. We define an external computation point P (r, φ, λ) located at radius r, geocentric latitude φ, and longitude λ. Grombein et al. (2013) provide efficient formulations for the volume integrals of the gravitational potential and its first-and second-order derivatives of a tesseroid with homogeneous density. The derivatives of the gravitational potential are taken with respect to the local north-oriented Cartesian coordinate system with origin at P (Fig. 1). Here, we will assume that the tesseroid has a density varying with r according to an arbitrary continuous function ρ(r). Thus, the integrals for the gravitational potential and its first-order derivatives are slightly modified to and in which α ∈ {x, y, z}, G = 6.674 × 10 −11 m 3 kg −1 s −1 is the gravitational constant and = r 2 + r 2 − 2rr cos ψ, cos ψ = sin φ sin φ + cos φ cos φ cos(λ − λ).

Gauss-Legendre Quadrature integration
Applying a N th order GLQ, we can approximate each integral in equations 1 and 2 by a weighted sum of the integration kernel evaluated on the roots of an N th order Legendre polynomial (Hildebrand, 1987, p. 390). Unlike the homogeneous density case, the radial density function ρ(r) must also be included in the integration and evaluated on the Legendre polynomial roots (i.e. quadrature nodes). where f (r , φ , λ ) is an integral kernel for a homogeneous tesseroid (Grombein et al., 2013), (r i , φ j , λ k ) are the coordinates of the quadrature nodes, N r , N φ , M λ are the quadrature orders and W r i , W φ j , W λ k are the quadrature weights in the radial, latitudinal, and longitudinal directions, respectively. Hildebrand (1987, p. 391) provides the formulas for calculating the GLQ weights and pre-computed values of the Legendre polynomial roots for low orders. See Uieda et al. (2016) for a more detailed description using a similar notation to the one used here. It is worth noting that a GLQ is equivalent to approximating the tesseroid by N r × N φ × N λ point masses located on the quadrature nodes (Ku, 1977;Asgharzadeh et al., 2007). Ku (1977) noticed that the GLQ integration becomes less accurate when the computation point is closer to the mass element. One way to prevent this from happening would be to increase the GLQ order. Doing so would uniformly increase the number of point masses inside the tesseroid volume. However, an increase in the point mass concentration is only required close to the computation point (Uieda et al., 2016). Alternatively, Li et al. (2011) proposed an adaptive discretization algorithm which keeps the GLQ order fixed and divides the tesseroid based on a ratio between the distance to the computation point and its dimensions. This algorithm produces a more efficient computation because an increased concentration of point masses is produced where it is needed more. Uieda et al. (2016) developed a modified version of this algorithm along with an efficient computational implementation. Both algorithms perform tesseroid subdivisions in the latitudinal, longitudinal and radial directions, thus we can define them as three-dimensional adaptive discretization algorithms. On the other hand, Lin & Denker (2019) proposed a two dimensional discretization algorithm that subdivides the tesseroid only on the latitudinal and longitudinal directions. Removing a dimension from the discretization makes the computation more efficient by reducing the number of tesseroids in the model, while retaining an acceptable accuracy (Lin & Denker, 2019).

Two Dimensional Adaptive Discretization
Here we will follow Lin & Denker (2019) and use a two dimensional version of the adaptive discretization of Uieda et al. (2016). What follows is a summary of the algorithm and the reader is referred to Uieda et al. (2016) for a detailed description.
Step 1: Check that the tesseroid satisfies the following inequality for its latitudinal and longitudinal dimensions L i (i ∈ {λ, φ}): in which D is a positive scalar called the distance-size ratio, d is the distance between the computation point and the geometric centre of the tesseroid and the dimensions of the tesseroid are defined as L λ = r 2 arccos(sin 2 φ t + cos 2 φ t cos(λ 2 − λ 1 )), and L φ = r 2 arccos(sin φ 2 sin φ 1 + cos φ 2 cos φ 1 ).
Step 2: If none of the dimensions of the tesseroid fail inequality 11, then compute the gravitational effect of the tesseroid using a second-order GLQ (Eq. 9). Add the computed effect to a running total.
Step 3: If inequality 11 does not hold for one of the dimensions (longitudinal or latitudinal), split the tesseroid in half along that dimension. Repeat steps 1-3 for all smaller tesseroids until none are left that violate inequality 11.
Final step: By the end of the algorithm, a second-order GLQ will have been applied to each smaller tesseroid and the running total gathered in step 2 will be the gravitational effect of the original tesseroid.
The distance-size ratio D determines how many times the tesseroids will be divided. Therefore, it effectively regulates both the accuracy of the algorithm and its computation time. An optimal value for D cannot be directly calculated from the desired accuracy level. Instead, it is empirically determined by comparing the numerical results with the analytical solution for a spherical shell. Uieda et al. (2016) used a shell with homogeneous density to determine optimal values of D. Here, we will repeat the numerical experiment using analytical expressions for shells with density varying according to functions of r. This experiment will also test whether the same values of D determined by Uieda et al. (2016) can be used for the two dimensional adaptive discretization.

Density-based Discretization Algorithm
The numerical integration of an arbitrary continuous density function introduces a new type of problem: the integration error from using only a few quadrature nodes to account for the variation of the density function. The three dimensional adaptive discretization may help to reduce this kind of error by adding more point masses in the radial direction. However, it does not take into account the density function itself when dividing the tesseroid and hence it is not well suited to fully perform this task.
We have developed a complementary discretization algorithm in the radial direction that takes into account the variations of the density function. This density-based discretization is applied prior to the two dimensional adaptive discretization described in the previous section. In short, the algorithm divides the tesseroid along the radial dimension at the depths at which the maximum density variations take place.
Consider an original tesseroid with a density given by the continuous function ρ(r ). Before the density-based discretization starts, we normalise the density function to the range [0, 1] as follows in which ρ min and ρ max are the minimum and maximum density values inside the tesseroid boundaries. We emphasize that this normalised density function will not be modified throughout the algorithm. In case the density function is constant, both maximum and minimum densities will be equal and the density-based discretization algorithm will not be applied. The algorithm is comprised of the following steps ( Fig. 2): Step 1: Define a straight line ρ l (r ) that assumes the same values as the normalised density function ρ n (r ) at the boundaries of the tesseroid (r 1 and r 2 ): Step 2: Evaluate the normalised density function and the straight line on a range of N radii between r 1 and r 2 . We have opted for N = 101 but the specific value of N is not critical to the algorithm.
Step 3: Compute the absolute difference between the values of the straight line and the normalised density function: Step 4: If the following inequality holds, the tesseroid will not be divided: in which L r is the radial dimension of the tesseroid being considered for division, L orig r is the radial dimension of the original tesseroid, and δ is a positive constant henceforth called the delta ratio.
Step 5: If inequality 20 is not satisfied, then the tesseroid is split in two parts at the radius r max at which the maximum absolute difference (Eq. 19) takes place. Repeat steps 1-5 for each smaller tesseroid produced in this step.
Final step: Once all smaller tesseroids satisfy inequality. 20, each one is subjected to the two dimensional adaptive discretization algorithm described earlier to calculate their gravitational effects.
On the first iteration, the ratio L r /L orig r equals one because the tesseroid being divided is the original one. For further iterations, the ratio will be progressively smaller than one as the tesseroids get smaller. This is intended to limit the number of divisions to the ones that will significantly reduce the numerical error: dividing a large tesseroid with a small max{∆ρ(r )} would improve the integration accuracy more than dividing a small tesseroid with a higher max{∆ρ(r )}.
The higher δ is, the fewer divisions will be made, and viceversa. Thus, it controls how many times the tesseroids will be divided based on the density function and, indirectly, determines the accuracy and computation time of numerical integration. This raises the need to determine a maximum value of δ that ensures an acceptable accuracy while minimising the computation time.

Algorithm summary
In summary, given a tesseroid with density varying in depth according to an arbitrary continuous function, we propose the following steps to numerically compute its gravitational fields on an external point: Step 1: Apply the density-based discretization algorithm to discretize the tesseroid in the radial dimension, producing a set of tesseroids with the same longitudinal and latitudinal dimensions as the original one but with different radial boundaries.
Step 2: Apply the two dimensional adaptive discretization algorithm for each tesseroid obtained in the previous step. If needed, the algorithm will divide each tesseroid in the latitudinal and longitudinal direction, generating a set of smaller tesseroids.
Step 3: Apply a second-order GLQ to numerically compute the gravitational fields (Eq. 9) generated by each tesseroid obtained in the previous step. The numerical integration includes the density function and can be applied without modification to any continuous function. The sum of these results is the gravitational field of the original tesseroid.

Software implementation
We have implemented the algorithms described in the previous sections in the Python programming language. The software is based on the Python implementation of Uieda et al. (2016) that is included in the Fatiando a Terra library v0.5 (Uieda et al., 2013). The more time consuming parts of the algorithm are written in the Cython language to achieve higher performance. We leverage the dynamic nature of the Python language to allow user-defined radial density functions as inputs to the software. Thus, our code can evaluate linear, exponential, polynomial, sinusoidal, cubic splines, or any other continuous density function without modification. This new code is freely available under the BSD 3-clause opensource license. All source code, Python scripts, data, and model results are made available through an online repository doi.org/10.6084/m9.figshare.8239622 (Soler et al., 2019) or github.com/pinga-lab/tesseroid-variable-density. The repository also contains instructions for replicating all results presented here.
3 Determination of the distance-size and delta ratios The distance-size ratio D of the adaptive discretization and the delta ratio δ of the density-based discretization determine how many times each tesseroid will be divided and thus indirectly control the numerical error of the integration. Optimal values for D and δ must be determined in order to ensure both acceptable numerical accuracy and computation efficiency for the algorithm. Uieda et al. (2016) compared the numerical integration of homogeneous density tesseroids with the analytical solution of a spherical shell (Mikuška et al., 2006;Grombein et al., 2013) in order to obtain default values for the distance-size ratio D. We will follow this idea but for our needs the spherical shell must have the same density function of radius as our tesseroid model. Lin & Denker (2019) show the analytical solution of the gravitational potential generated by a spherical shell with linear density in the radial coordinate. Applying the Newton's Shell Theorem (Chandrasekhar, 1995;Binney & Tremaine, 2008), we derive expressions for the gravitational potential of a spherical shell with linear, exponential, and sinusoidal density functions (see Appendix A).
In order to compare the numerical results with the analytical solution we must build spherical shell models out of tesseroids. We divide the spherical shell along the latitudinal and longitudinal directions to obtain a shell model made out of 6 × 12 = 72 tesseroids of size 30 • × 30 • . In order to asses the density-based discretization in the radial dimension, we use spherical shell models with different thicknesses (Table 1). The thickness values were chosen to represent a range of applications, from topographic to lithospheric scale models. Because the amount of tesseroid divisions in the adaptive discretization will be proportional to the tesseroid size (Eq. 11), some of these configurations represent worse-case scenarios. Most practical applications will use tesseroids smaller than 30 • × 30 • × 1000 km.
The differences between the analytical and numerical solutions could be calculated on a single computation point due to the rotational symmetry of the spherical shell. However, the Table 1: Description of the tesseroid models used to build spherical shells and characterize the accuracy of the numerical integration. The outer radius (R 2 ) of every shell model is equal to the mean Earth radius (6378.137 km), while the inner radius (R 1 ) is determined by its thickness. The horizontal dimensions of the tesseroids and the total number of tesseroids in the shell model are given in the latitudinal and longitudinal dimensions, respectively.
Thickness Tesseroid size Number of tesseroids 0.1 km numerical results depend on the relative location between the computation point and the tesseroid (Ku, 1977;Asgharzadeh et al., 2007;Uieda et al., 2016). We account for this effect by calculating the differences on regular grids and keeping only the maximum absolute difference. All calculations are repeated for four different computation grids (Table 2): a local grid at the pole, a local grid at the equator, a global grid at zero height above the shell, and a global grid at a height of 260 km (representing the nominal height of the GOCE satellite). These grids span a broad scenario of use cases and ensure an acceptable accuracy on each one of them. We perform comparisons between the analytical solutions for the spherical shell and the numerical integration results for linear, exponential, and sinusoidal density functions. The sinusoidal density case is included in order to test the numerical approximation to its limits. The comparisons are repeated for all combinations of tesseroid models in Table 1 and computation grids in Table 2. From these results, we generalize optimal values for D and δ that ensure a numerical error lower than 0.1% of the spherical shell values for most use cases.

Linear Density
A spherical shell with a linear density function given by has an analytical solutions for the gravitational potential and its vertical derivative given by Eqs. 34 and 32. The values of the angular and linear coefficients (a and b) can be chosen so that the density assumes the values of ρ in = 3300 kg/m 3 and ρ out = 2670 kg/m 3 on the inner (R in ) and outer (R out ) radii of the shell, respectively, The absolute density difference defined on equation 19 will always be zero for the linear density case. As a result, the inequality 20 will always be satisfied and no divisions will ever be performed during the density-based discretization. Therefore, the two dimensional adaptive discretization algorithm is the only mechanism that controls the accuracy of the numerical integration for the linear density case. For this reason, we will ignore the values of δ and only determine the minimum value of D needed in order to guarantee an acceptable accuracy.
We compute the gravitational potential (V ) and its vertical derivative (g z ) for each shell model in Table 1 on each computation grid in Table 2. The horizontal derivatives of the potential are equal to zero outside of the shell due to the rotational symmetry and are thus omitted from the analysis. The computations are repeated for values of the distance-size ratio D ranging from 0.5 to 5 in increments of 0.5. We then calculate the absolute difference between the numerical results and the analytical solution for the shell. Fig. 3 shows the maximum absolute difference for each shell model and computation grid as a function of D. The differences are relative to the shell value. Finally, we set the optimal value of D as the smallest value for which the difference is below 0.1%.
We observe from Fig. 3 that the relative errors for the potential and g z fall below the 0.1% threshold at D = 1 and D = 2.5, respectively. Notably, a value of D = 2 would suffice for g z in the case of the satellite height grid. For all other configurations, these values are consistent and independent of the shell model thickness or geographic location.

Exponential Density
For an exponential density function, the density-based discretization will be applied before the adaptive discretization algorithm. This means that optimal values for both the distance-size ratio D and the delta ratio δ (Eq. 20) must be determined. We perform an error analysis similar to what was done for the linear density case. Now the spherical shell will have an exponential density function that assumes the values of ρ out = 2670 kg/m 3 and ρ in = 3300 kg/m 3 on the outer and inner surfaces, respectively, defined as follows: where and b is a dimensionless constant that determines the variability of the function. A higher value of b increases the maximum slope of the density function (Fig. 4).

D-δ space exploration
We aim to find a combination of the D and δ that produces a numerical error lower than the 0.1% threshold while minimizing computation time. We use a grid search method and compute the numerical error for every (D, δ) pair belonging to a grid on the D-δ space (Fig. 5). For optimum algorithm performance, we search for the (D, δ) pair that minimizes the number of tesseroid division while keeping a numerical error under the 0.1% threshold. This requirement translates into the smallest possible value of D and the highest possible value of δ. Because this is a time consuming computation, we limit the analysis to a high value of b = 30 (Fig.4) and the global grid defined in Table 2. We compute the relative difference between the numerical and analytical results of the gravitational potential (V ) and its vertical derivative (g z ) for all shell models defined in Table 1 The smallest values of D that are within the 0.1% threshold coincide with D linear for both V and g z . These results indicate that D linear can be safely extrapolated to non-linear density cases. This is not surprising considering that the two dimensional adaptive discretization and the density-based discretization are independent of each other. The former divides the tesseroid in the horizontal dimensions, while the latter only divides in the radial dimension. Hence, the optimal value of δ is likely to be independent of the optimal value of D. Because the grid search was limited to a specific computation grid and value of b, we perform a more detailed analysis in the  Figure 3: Differences between the gravitational fields generated by each tesseroid shell model and the analytical solution as a function of the distance-size ratio D. Each model has a linear density function (Eq. 22). The computations were performed on the four grids described in Table 2 and using the shell models detailed in Table 1. Each line represents the maximum absolute difference between the numerical results and the analytical solution for a given shell model. Due to the linearity of the density function, the density-based discretization algorithm is not applied. Differences are reported as a percentage of the analytical solutions. The horizontal dashed black line represents a target difference of 0.1%. following section to determine an optimal value of δ for the exponential density case.

Delta ratio determination
Having chosen values of D equal to the ones obtained for the linear density case, we are free to explore the integration error as a function of δ in more detail. We will compute the difference between the numerical and analytical results for all combinations of computation grid (Table 2) and spherical shell model (Table 1), varying δ from 10 −3 to 10 1 . The calculations are repeated for each b ∈ {1, 2, 5, 10, 30, 100} (Fig. 4) to examine the accuracy of the method for density functions of different sharpness. Because larger δ values result in fewer tesseroid divisions, our intention is to find the highest value of δ that produces a relative difference bellow the 0.1% threshold. Fig. 6 shows the resulting relative differences for V and g z as a function of δ. For the sake of brevity, each curve corresponds to the maximum difference between all shell models. The curves for b = 1 and b = 2 are below the 0.1% threshold for all values of δ and do not change for δ > 0.8, indicating that the density functions are sufficiently smooth and no density discretizations are necessary. For all other values of b, the difference falls below the 0.1% threshold if δ = 0.1. These results indicate that there is no significant relationship between the sharpness of the exponential density function and the numerical error.

Sinusoidal Density
So far we have tested the density-based discretization algorithm against linear and exponential density functions. Nev- The percentage difference values were obtained from the comparison between the analytical solution and the numerical approximation of the gravitational fields (V and g z ) generated by a spherical shell with an exponential density function (Eq. 25). These comparisons were carried out on the "Global" grid (Table 2), with the spherical shell models detailed in Table 1, and an exponential density function with b = 30. The percentage difference values are obtained as the maximum difference between every shell model. The points inside the dashed line are the ones that present an error lower than 0.1%. The D value obtained for the linear density case for each gravitational field is also shown (D linear ).  Figure 6: Difference between the numerical and analytical solutions as a function of δ ratio for different exponential density functions. These comparisons were carried out for every shell model (Table 1) and computation grid (Table 2) using a fixed value of D. Each curve corresponds to the maximum difference out of all shell models for a particular value of b (Eq. 25). Differences are reported as a percentage of the analytical solutions. The horizontal dashed black line represents a target difference of 0.1%. ertheless, this new algorithm is well suited to more complex continuous functions, for example non-monotonic functions or those with multiple inflection points. Although such density functions are rarely, if ever, seen in geological structures, we want to subject our algorithm to such cases in order to show that it can solve more complex situations than linear and exponential functions. Furthermore, acquiring an optimal value for the δ ratio in case of unrealistic density functions allow to extrapolate it to simpler realistic cases. Thus this tests are only intended to submit the algorithm to the extreme and not to emulate a real case scenario.
We will consider spherical shells with a sinusoidal density function defined as follows: in which A is a constant that controls the amplitude and vertical shift of the sine function, R is the mean Earth radius, and b is a dimensionless constant that regulates how many periods of the trigonometric function are included inside the inner and outer radii. The analytical solutions for V and g z of a spherical shell with a sinusoidal density function can be found on Appendix A. We compute the relative difference between the numerical and analytical results for V and g z for all combinations of spherical shell models (Table 1) and computation grids (Table 2). We fixed the distance-size ratio D to the ones obtained for the linear density case and explored values of δ ranging from 10 −4 to 1. The calculations are repeated for values of b equal to 1, 2, 5 and 10 (Fig. 7). For all computations, the value of A is fixed at 1650 kg/m 3 . Fig. 8 shows the relative differences between the analytical and the numerical solutions for sinusoidal density case. Once again, each curve corresponds to the maximum difference between all shell models. For all values of b except for b = 10, the differences fall below the 0.1% threshold for δ = 0.1. In the case of b = 10, a lower value of δ = 0.01 is required to achieve 0.1% difference. We note that, even for the case of b = 10, the differences for δ = 0.1 are below 1%.

Algorithm performance
Because the density-based discretization algorithm introduces more divisions along the radial dimension, it is reasonable that the computation time for variable densities will be higher than the homogeneous density case. Directly comparing the computation time of both cases cannot be done in a meaningful way because it is highly dependent on the implementation of the algorithm, the choice of programming language, and the particular density function used. In order to obtain an indicator of the increase in the computation time, we chose the number of tesseroid divisions in the density-based discretization as a proxy measure.  Figure 8: Difference between the numerical and analytical solutions as a function of δ ratio for different sinusoidal density functions. These comparisons were carried out for every shell model (Table 1) and computation grid (Table 2) using a fixed value of D. Each curve corresponds to the maximum difference out of all shell models for a particular value of b (Eq. 28). Differences are reported as a percentage of the analytical solutions. The horizontal dashed black line represents a target difference of 0.1%.
We analyse exponential density functions (Eq. 25) with values of b equal to 1, 2, 5, 10, 30 and 100, and sinusoidal density functions (Eq. 28) with values of b equal to 1, 2, 5, and 10. We then apply the density-based discretization algorithm on each function and record the number divisions carried out. Figs. 9a and 9b show the exponential and sinusoidal density functions and the resulting discretization points as orange dots.
For the exponential density functions, the algorithm performs a single division independently of the value of b. Therefore, the computation time would likely be at least twice that of the homogeneous density case (assuming equal implementations). On the other hand, Fig. 9c shows an almost linear relation between the number of divisions in the sinusoidal density case and the value of b. Thus, the computation time is likely to be dependent on the number of wavelengths of the sinusoidal function that are contained within the tesseroid.

Application to the Neuquén Basin
We applied the new algorithms and optimal values of D and δ determined previously to calculate the gravitational effects of the Neuquén Basin, a sedimentary basin located to the east of the Andes, between 32 • S and 40 • S latitude (Fig. 10a). The basin includes continental and marine siliciclastics, carbonates, and evaporites accumulated over the Jurassic and the Cretaceous constituting a stratigraphic record up to 5000m of depth (Howell et al., 2005).
The thickness of the sediment pack was digitized from Heine (2007) on a regular grid with a resolution of 0.05 • on both longitude and latitude directions (Fig. 10b). We created a tesseroid model of the sediment pack by placing a 0.05 • × 0.05 • tesseroid on each node of the grid. The top of each tesseroid was fixed at the median of the topography of the basin (845 m above mean Earth radius) and the bottom at corresponding thickness of the basin.
We must also define a density function for the tesseroid model. Sigismondi (2012) measured a minimum and maximum density contrast for the Neuquén basin of -412kg/m 3 and -275kg/m 3 , respectively. We have chosen an exponential density variation (Eq. 25) that assumes the minimum value on the top surface and the maximum at 5014m depth (the thickest part of the basin), with a value of b equal to 3. This density variation is in the order of magnitude of the ones used by Cowie & Karner (1990) and Cordell (1973). This density function can be seen on Fig. 11.
Finally, we computed the gravitational potential V and the vertical component of the gradient (g z ) on a computation grid of 159 × 163 nodes (0.05 • spacing on both longitude and latitude) at a 10 km height over the mean Earth radius. The resulting fields can be seen in Fig. 10c-d. We computed the differences between the results for the exponential density function and those generated by the same model but now with a constant density (the mean of -412kg/m 3 and -275kg/m 3 ). We also computed the differences between the exponential density and a linear density that assumes these two values on the highest and the deepest point of the basin (Fig. 11). Fig. 12a-b and Fig. 12c-d show the differences with the constant and the linear density, respectively. The maximum absolute difference in the computed g z is approximately 8 mGal for the homogeneous case and approximately 6 mGal for the linear density case. Both are well above the accuracy of most available data products.

Discussion
By including the density function into the GLQ integration, the method described here can be applied without modification to tesseroids with density following any continuous function of radius. A density-based discretization algorithm divides the tesseroid in the radial dimension to ensure accurate integration of the density function. This algorithm is independent of the GLQ integration and could potentially be used to determine an optimal discretization when approximating a density function by piecewise linear (Lin & Denker, 2019) or piecewise polynomial (Fukushima, 2018) functions.
Our numerical experiments show that the two dimensional adaptive discretization is enough to achieve 0.1% accuracy with a second-order GLQ in the case of a linear density function (Fig. 3). Values of the distance-size ratio determined here for the gravitational potential (D = 1) and its vertical derivative (D = 2.5) are compatible with Uieda et al. (2016). These results also show that there is no significant relationship between the accuracy of the method and the thickness of the tesseroid model.
The D-δ space exploration for the exponential density function (Fig. 5) showed that values of D determined for the linear density case are equally applicable to the exponential density case. Furthermore, the difference with respect to the analytical solution only falls below 0.1% for values of δ lower than 0.1. It follows that the density-based discretization is required to achieve the desired accuracy level for non-linear density functions.
More detailed error analyses showed that δ = 0.1 is sufficient to guarantee an accuracy of 0.1% for all exponential density functions tested (Fig. 6) and most of the sinusoidal functions tested (Fig. 8). The exception is the case of a sinusoidal function with b = 10 (Eq. 28), for which δ = 0.01 is required to achieve 0.1% accuracy. Nevertheless, using δ = 0.1 in this case would still produce results with an error of less than 1%.
While the algorithms proposed here can be applied to any continuous density function, the optimal values for D and δ are empirically determined only for linear, exponential, and sinusoidal density functions. Thus, we can only say with certainty that these values will produce results with 0.1% accuracy for these functions. However, all of our numerical tests include worse-case scenarios (zero computation heights, large tesseroids, highly variable density functions, etc). For this reason, these optimal values of D and δ can plausibly be extrapolated to any realistic continuous density function. Notwithstanding, we encourage users of the algorithms and software to perform similar tests in order evaluate the accuracy when using density functions more complex than the ones tested here.
The algorithm performance analysis shows that the computation time using variable densities is likely to be at least two times larger than the homogeneous density case. It is also likely to increase proportionally to the number of inflection points of the density function. As is the case for most numerical methods, there is a trade-off between computation time and accuracy. Nevertheless, density profiles will have few inflection points for most geophysical applications. Therefore, the computation time would be in the same order of magnitude as the one for the homogeneous density in most real world applications.

Conclusions
We have developed a new methodology to compute the gravitational potential and its gradient for a tesseroid with a density given by a continuous function of radius. It numerically solves the volume integrals through Gauss-Legendre Quadrature (GLQ) by including the density function in the numerical integration. By implementing the algorithm in the dynamic Python programming language, users can define their own density function to supply to the software. This allows the use of any arbitrary continuous function without modification to the method or software. The accuracy of the numerical integration is automatically controlled by a two dimensional adaptive discretization algorithm and a new density-based discretization algorithm. The former divides the tesseroid in half if the ratios of the distance to the computation point and the latitudinal or the longitudinal size of the tesseroid are smaller than a predefined distance-size ratio D. This algorithm minimises the integration error when the computation point is close to the tesseroid. Nevertheless, the adaptive discretization alone is not sufficient to guarantee the accuracy of the method in case of tesseroids with non-linear density functions. To overcome this challenge, the density-based discretization algorithm divides the tesseroid along the radial dimension on the points at which the maximum variations of the density function take place. The number of radial divisions performed, and thus the accuracy of the computation, is controlled by the δ parameter. This new algorithm is intended to minimise the error due to the inability of the GLQ to produce precise approximations of density functions with sharp variations.
We empirically determined optimal values for D and δ by comparing the numerical results with analytical solutions. Our analysis included a range of tesseroid models and computation grids as well as common density functions. These values minimize the computational load while maintaining the numerical error below 0.1% of an analytical solution. The density functions used to establish the optimal values were a linear, an exponential, and a sinusoidal function. The linear function represents the smoothest variation of density and does : Linear and exponential densities used to compute the gravitational fields generated by a tesseroid model of the Neuquén sedimentary basin. The height is defined above the mean Earth radius, and its axis is spanned between the deepest and the highest point of the basin. not require density-based discretization at all. Through this analysis, we have obtained optimal values for the distancesize ratio D of 1 and 2.5 for the potential and its gradient, respectively. We analysed the error for exponential functions ranging from smooth to sharp and sinusoidal function of different wavelengths to test the accuracy of the density-based discretization. It worth noting that the sinusoidal density tests were carried out in order to test the algorithm to the extreme and not to emulate a real world scenario. Values of δ equal to 0.1 are sufficient to guarantee a 0.1% accuracy for most realist cases tested. These results could be extrapolated to other continuous density functions that are sufficiently smooth, as is the case for most common geophysical applications. Smaller values of δ can be used for more highly variable density functions to increase the accuracy of the integration accordingly.
Some computational performance is sacrificed in order to obtain an accurate and general purpose method. The densitybased discretization adds to the computation time by increasing the number of GLQ integrations. Likewise, allowing users to supply custom density functions incurs an overhead cost because of the added function evaluations. It is also impossible to optimize the source code and formulas without knowledge of the specific density function. Nonetheless, the densitybased discretization is independent of the adaptive discretization and GLQ integration. It can be viewed as a type of preprocessing step that can be combined with other more specialized integration methods.
An application to modelling the Neuquén Basin, Argentina, demonstrates that the effects of sediment compaction cannot be ignored. When compared to an exponential density function, a homogeneous density with depth would result in an error of up to -8 mGal, while a linear density function would result in up to -6 mGal error. Accurate and robust forward modelling is a key component to any gravity inversion method. Errors of this magnitude could result in significant overestimation of sedimentary basin thickness, for example.  Figure 12: Differences between the gravitational fields generated by the tesseroid model of the Neuquén basin with an exponential density contrast and with a homogeneous and linear density variation. (a)-(b) Differences on V and g z between the exponential density model and the homogeneous density one, (c)-(f) differences on V and g z between the exponential density model and the linear density one, calculated at 10km of height over the mean Earth radius. Zhang, J., Zhong, B., Zhou, X., & Dai, Y., 2001. Gravity anomalies of 2D bodies with variable density contrast, Geophysics, 66 (3), 809-813.

A Analytical Solutions for Spherical Shell
Consider a spherical shell with inner radius R 1 and outer radius R 2 , whose density is a function ρ(r ) of the radial coordinate. We want to get an analytical expression for the gravitational fields the shell generates on any external point located at a distance r from the centre of the shell (r > R 2 ). According to Newton's Shell Theorem (Theorem XXXI) (Chandrasekhar, 1995;Binney & Tremaine, 2008), the gravitational potential generated by the shell on any external point is the same as it would be if its mass were concentrated on a point located at its centre: where M is the total mass of the shell, which can be easily computed as follows: where Ω symbolizes the volume of the shell. Combining the two previous equations, we get the following expression for the potential: which is in agreement with the one obtained by Binney & Tremaine (2008, p.62). The gradient of potentials that depends solely on r have only one non zero components: the vertical component of the gradient (g z ). Following Grombein et al. (2013): From Eq. 31 we can obtain expressions of the gravitational potential for different density functions. The integration of the following density functions have been carried out by SymPy (Meurer et al., 2017), a Python library for symbolic mathematics.

A.1 Linear density
For a linear density function ρ(r ) = ar + b , the gravitational potential at any external point is The first term on this equation reproduces the potential generated by a spherical shell with variable density ρ(r ) = ar , while the second term constitutes the potential generated by a spherical shell with homogeneous density ρ = b (Mikuška et al., 2006;Grombein et al., 2013). Eq 34 is in agreement with the one obtained by Lin & Denker (2019).