# Comprehensive gravitational modeling of the vertical cylindrical prism by Gauss–Legendre quadrature integration

Comprehensive gravitational modeling of the vertical cylindrical prism by Gauss–Legendre... SUMMARY Forward modeling is the basis of gravitational anomaly inversion that is widely applied to map subsurface mass variations. This study uses numerical least-squares Gauss–Legendre quadrature (GLQ) integration to evaluate the gravitational potential, anomaly and gradient components of the vertical cylindrical prism element. These results, in turn, may be integrated to accurately model the complete gravitational effects of fluid bearing rock formations and other vertical cylinder-like geological bodies with arbitrary variations in shape and density. Comparing the GLQ gravitational effects of uniform density, vertical circular cylinders against the effects calculated by a number of other methods illustrates the veracity of the GLQ modeling method and the accuracy limitations of the other methods. Geological examples include modeling the gravitational effects of a formation washout to help map azimuthal variations of the formation’s bulk densities around the borehole wall. As another application, the gravitational effects of a seismically and gravimetrically imaged salt dome within the Laurentian Basin are evaluated for the velocity, density and geometric properties of the Basin’s sedimentary formations. Geopotential theory, Gravity anomalies and Earth structure, Downhole methods, Numerical modelling, Statistical seismology 1 INTRODUCTION Calculating the gravitational effects (i.e. the gravitational potential and its higher order spatial gradients) of an extended mass body requires integrating at the observation point the differential gravitational effects of the unit volume elements making up the body. Any volume symmetries of the source considerably simplify the related scalar, vector and tensor integrals. Thus, for example, the gravitational effects of the prism’s mass are readily accessible in Cartesian (Nagy 1966), spherical (Asgharzadeh et al. 2007) and elliptical (Roussel et al. 2015) coordinates. The considerable interest in computing the gravitational effects of vertical cylinder-like bodies has spawned a relatively large number of quantitative modeling approaches, which generally consider the case of either the finite or semi-infinite vertical axis. Nettleton (1942), for example, presented a closed formula for the vertical-component gravitational attraction ($$g_{o,Z_o}$$) of the vertical semi-infinite cylinder (VL) where the cylinder’s radius is considerably smaller than its thickness. Nettleton (1942) also presented the closed formula for $$g_{o,Z_o}$$ based on the solid angle (SA) subtended at the observation point O by the vertical circular cylinder with thickness that is less than or equal to about a half of its average depth. Van Baak (1989) calculated $$g_{o,Z_o}$$ of the vertical cylinder using line integrals (LIs) around the boundaries of the volume-filling discs. The cross-sectional integral of each disc is obtained from the bounding curve either in closed analytical form if it is a complete circle or numerically as the sum of the n line segment integrals if an approximating polygon with n line vertices is applicable. The line segment integrals for all discs are evaluated and summed at each observation point with an accuracy that increases as the number of line segments increases in the bounding curve approximations of the discs. Additional closed $$g_{o,Z_o}$$ formulae are available for the circular vertical cylinder (P) using Cartesian coordinate integrals (Parasnis 1961), as well as elliptic integrals (EIs) of the first and second kind and Heuman’s Lambda function (Nabighian 1962) and the Fourier–Hankel (FH) transform (Singh 1977). Rim & Li (2016) presented closed expressions of $$g_{o,Z_o}$$ and its gradients for the vertical and tilted cylinder. Numerical approximations of $$g_{o,Z_o}$$ are also available involving the first few terms of the Legendre polynomial (LP) expansion (e.g. Parasnis 1961; Reilly 1969; Jenkins et al. 1983; Telford et al. 1990) and the conjugate complex (CC) variables formulation of Green’s Theorem (Kwok 1991). In general, the above gravitational effect estimates typically assume a constant density throughout the cylinder, although Jenkins et al. (1983) present a method to accommodate density variations as annuli and discs along the respective radial and axial dimensions. However, they do not accommodate azimuthal density variations within the cylinder, which may appear, for example, along borehole walls through physically or chemically altered fluid-bearing rock formations. These alterations occur in various sizes to complicate borehole gravimetry applications for interpreting surface gravity, seismic and borehole geology observations (e.g. Smith 1950). Borehole gravity data also constrain bedding dip estimates assuming the formation densities are available (Brown & Lautzenhiser 1982). For dipping beds, Nekut (1989) showed that borehole gravity gradiometry also can be used to estimate the bedding dips, strikes and densities, as well as the average density of the hosting rock formation. Modern applications for borehole gravity gradiometry also facilitate CO2 sequestration, enhanced oil recovery and mining and exploration applications (Rim & Li 2012). In this paper, the complete gravitational effects of the vertical cylinder with discrete variations of density along the radial, azimuthal and axial dimensions are developed in the context of the cylindrical unit volume element or prism shown in Fig. 1. These effects are numerically calculated from the Gauss–Legendre quadrature (GLQ) distribution of point gravity poles within the unit prism elements that fill up the cylinder (Stroud & Secrest 1966). These equations are applicable for calculating the complete gravitational effects of salt domes, granitic and other igneous plume-like intrusions, carbonate reefs, gas-hydrate columns, kimberlite pipes, crustal meteorite impact structures, mantle plumes, formation washouts and dipping rock formations along wellbores, and other homogeneous or heterogeneous vertical partial or full cylinder-like geological features. The accuracy of the GLQ-based gravitational modeling is tested by comparing the $$g_{o,Z_o}$$ estimates against those of the other methods described above that vary with changes in the cylinder’s radius, thickness and burial depth. Figure 1. View largeDownload slide The cylindrical ρϕZ coordinate system for modeling at the observation point O(ρo, ϕo, Zo), the gravitational effects of a vertical cylindrical prism by integrating the differential source mass dms through the volume of the prism with bottom and top surfaces that include the corner points A(ρ1, ϕ1, Z1) and B(ρ2, ϕ2, Z2). The Cartesian perspective is given by the superposed, co-registered Cartesian XYZ coordinate system axes. Figure 1. View largeDownload slide The cylindrical ρϕZ coordinate system for modeling at the observation point O(ρo, ϕo, Zo), the gravitational effects of a vertical cylindrical prism by integrating the differential source mass dms through the volume of the prism with bottom and top surfaces that include the corner points A(ρ1, ϕ1, Z1) and B(ρ2, ϕ2, Z2). The Cartesian perspective is given by the superposed, co-registered Cartesian XYZ coordinate system axes. To help demonstrate the geological utility of the GLQ method, the gravitational effects of a formation washout are estimated in the context of the related borehole gravity measurements. The gravitational effects of a seismically constrained salt dome also are modeled within the Laurentian Basin off the eastern coast of Canada. Insights on the Basin’s average sediment density, velocity, and geometric properties are obtained from comparing the modeled gravitational anomaly to shipborne-constrained free-air gravitational anomaly estimates from the Eigen-6c4 spherical harmonic coefficient expansion (Förste et al. 2014). 2 GRAVITATIONAL EFFECTS FOR THE CYLINDRICAL PRISM Assuming (XYZ)-Cartesian and (ρϕZ)-cylindrical coordinate systems, with the vertical Z-axis upwards positive from origins on the orthogonal (XY)- and (ρϕ)-planes (Fig. 1), the far-field gravitational effects of a differential vertical cylindrical prism element in both coordinate systems is inversely proportional to the distance Ros between the observation (O) and source (S) points. This distance in cylindrical coordinates is   $${R_{os} = \left[{\rho _{o}^2 + \rho _{s}^2 - 2 \rho _{o} \rho _{s} {\rm cos}(\phi _{o} - \phi _{s}) + (Z_o - Z_{s})^2}\right]^\frac{1}{2},}$$ (1)where ρo and ρs are the radial distances of the respective observation and source points from the Z-axis, with respective ϕo- and ϕs-longitude coordinates, at the distances Zo and Zs from the origin’s XY- and ρϕ-planes. In addition, the differential volume dvs of the cylindrical prism’s differential mass dms is   $${dv_{s} = \rho _{s} d\rho _{s} d\phi _{s} dZ_{s}.}$$ (2) Thus, the total gravitational potential Po of the cylindrical prism at the observation point is the sum of the differential mass potentials given there by   $${P_{o} = \int _{V} dP_{os} = \int _{Z_{s} = Z_{1}}^{Z_{2}} \int _{\phi _{s} = \phi _{1}}^{\phi _{2}} \int _{\rho _{s} = \rho _{1}}^{\rho _{2}} \left[\frac{k}{R_{os}}\right] d\rho _{s} d\phi _{s} dZ_{s},}$$ (3)where k = Gσsρs, G = 6.67408 × 10−11 m3kg−1s−2 is the universal gravitational constant and σs is the mass density of the cylindrical prism. In addition, the spatial gradients of Po at the observation point in the ρo-radial, ϕo-longitudinal, and Zo-axial directions yield the vector gravitational field components ($$g_{{o,\rho _{o}}}$$, $$g_{{o,\phi _{o}}}$$, $$g_{{o,Z_o}}$$) given by   \begin{eqnarray} g_{{o,\rho _{o}}} &=& \int _{V} dg_{os,\rho _{o}} = \int _{Z_{s} = Z_{1}}^{Z_{2}} \int _{\phi _{s} = \phi _{1}}^{\phi _{2}} \int _{\rho _{s} = \rho _{1}}^{\rho _{2}} \nonumber\\ &&\left[-k \frac{(\rho _{o} - \rho _{s} {\rm cos} (\phi _{o} - \phi _{s}))}{R_{os}^3}\right] d\rho _{s} d\phi _{s} dZ_{s}, \end{eqnarray} (4a)  \begin{eqnarray} g_{{o,\phi _{o}}} &=& \int _{V} dg_{os,\phi _{o}} = \int _{Z_{s} = Z_{1}}^{Z_{2}} \int _{\phi _{s} = \phi _{1}}^{\phi _{2}} \int _{\rho _{s} = \rho _{1}}^{\rho _{2}} \nonumber\\ &&\left[-k \frac{\rho _{s} {\rm sin}(\phi _{o} - \phi _{s})}{R_{os}^3}\right] d\rho _{s} d\phi _{s} dZ_{s}, \end{eqnarray} (4b)and   \begin{eqnarray} g_{{o,Z_o}} &=& \int _{V} dg_{os,Z_o} = \int _{Z_{s} = Z_{1}}^{Z_{2}} \int _{\phi _{s} = \phi _{1}}^{\phi _{2}} \int _{\rho _{s} = \rho _{1}}^{\rho _{2}} \nonumber\\ &&\left[-k \frac{(Z_o - Z_{s})}{R_{os}^3}\right] d\rho _{s} d\phi _{s} dZ_{s}. \end{eqnarray} (4c) Furthermore, the spatial gradients of the gravitational field components of the cylindrical prism involve the elements of the 3 × 3 matrix given by   $$\nabla dg_{os} = \left( \begin{array}{c{@}{\quad}c{@}{\quad}c}dg_{os,\rho _{o},\rho _{o}} & {dg_{os,\rho _{o},\phi _{o}}} & {dg_{os,\rho _{o},Z_o}} \\ {dg_{os,\phi _{o},\rho _{o}}} & {dg_{os,\phi _{o},\phi _{o}}} & {dg_{os,\phi _{o},Z_o}} \\ {dg_{os,Z_o,\rho _{o}}} & {dg_{os,Z_o,\phi _{o}}} & {dg_{os,Z_o,Z_o}} \\ \end{array} \right)\!,$$ (5)where $$dg_{os,\rho _{o},\rho _{o}}$$ is the spatial gradient of $$dg_{os,\rho _{o}}$$ in the ρo-radial direction, $$dg_{os,\rho _{o},\phi _{o}}$$ is the spatial gradient of $$dg_{os,\rho _{o}}$$ in the ϕo-longitudinal direction and so on. These elements are formulated in Table 1 with the related partial differentials listed in Table 2. In general, the gravitational gradient tensor components at the observation point may be expressed over the source region by   $${g_{o,m,u} = \int _{V} dg_{os,m,u} = \int _{Z_{s} = Z_{1}}^{Z_{2}} \int _{\phi _{s} = \phi _{1}}^{\phi _{2}} \int _{\rho _{s} = \rho _{1}}^{\rho _{2}} dg_{os,m,u},}$$ (6)where m is the ρo-, ϕo-, or Zo-components of go and dgos, and u is the ρo-, ϕo-, or Zo-directions along which the gradient is taken. Table 1. Formulae for the gravitational tensor gradient components of the matrix in eq. (5). Component  [Formula]· dρsdϕsdZs  $$dg_{os,\rho _{o},\rho _{o}}$$  $$\frac{\partial \left(dg_{os,\rho _{o}}\right)}{\partial \rho _{o}}$$  $$dg_{os,\rho _{o},\phi _{o}}$$  $$\frac{1}{\rho _{o}} \left(\frac{\partial (dg_{os,\rho _{o}})}{\partial \phi _{o}} - dg_{os,\phi _{o}}\right)$$  $$dg_{os,\rho _{o},Z_o}$$  $$\frac{\partial (dg_{os,\rho _{o}})}{\partial Z_o}$$  $$dg_{os,\phi _{o},\rho _{o}}$$  $$\frac{\partial (dg_{os, \phi _{o}})}{\partial \rho _{o}}$$  $$dg_{os,\phi _{o},\phi _{o}}$$  $$\frac{1}{\rho _{o}} \left(\frac{\partial (dg_{os, \phi _{o}})}{\partial \phi _{o}} + dg_{os, \rho _{o}}\right)$$  $$dg_{os,\phi _{o},Z_o}$$  $$\frac{\partial (dg_{os,\phi _{o}})}{\partial Z_o}$$  $$dg_{os,Z_o,\rho _{o}}$$  $$\frac{\partial (dg_{os, Z_o})}{\partial \rho _{o}}$$  $$dg_{os,Z_o,\phi _{o}}$$  $$\frac{1}{\rho _{o}} \frac{\partial (dg_{os, Z_o})}{\partial \phi _{o}}$$  $$dg_{os,Z_o,Z_o}$$  $$\frac{\partial (dg_{os,Z_o})}{\partial Z_o}$$  Component  [Formula]· dρsdϕsdZs  $$dg_{os,\rho _{o},\rho _{o}}$$  $$\frac{\partial \left(dg_{os,\rho _{o}}\right)}{\partial \rho _{o}}$$  $$dg_{os,\rho _{o},\phi _{o}}$$  $$\frac{1}{\rho _{o}} \left(\frac{\partial (dg_{os,\rho _{o}})}{\partial \phi _{o}} - dg_{os,\phi _{o}}\right)$$  $$dg_{os,\rho _{o},Z_o}$$  $$\frac{\partial (dg_{os,\rho _{o}})}{\partial Z_o}$$  $$dg_{os,\phi _{o},\rho _{o}}$$  $$\frac{\partial (dg_{os, \phi _{o}})}{\partial \rho _{o}}$$  $$dg_{os,\phi _{o},\phi _{o}}$$  $$\frac{1}{\rho _{o}} \left(\frac{\partial (dg_{os, \phi _{o}})}{\partial \phi _{o}} + dg_{os, \rho _{o}}\right)$$  $$dg_{os,\phi _{o},Z_o}$$  $$\frac{\partial (dg_{os,\phi _{o}})}{\partial Z_o}$$  $$dg_{os,Z_o,\rho _{o}}$$  $$\frac{\partial (dg_{os, Z_o})}{\partial \rho _{o}}$$  $$dg_{os,Z_o,\phi _{o}}$$  $$\frac{1}{\rho _{o}} \frac{\partial (dg_{os, Z_o})}{\partial \phi _{o}}$$  $$dg_{os,Z_o,Z_o}$$  $$\frac{\partial (dg_{os,Z_o})}{\partial Z_o}$$  View Large Table 2. Formulae for partial differentials in Table 1 with k = Gσsρs. Partial differential  Formula  $$\frac{\partial (dg_{os,\rho _{o}})}{\partial \rho _{o}}$$  $$-k(R_{os}^2 - 3(\rho _{o} - \rho _{s}{\rm cos}(\phi _{o} - \phi _{s}))^2)/R_{os}^5$$  $$\frac{\partial (dg_{os,\rho _{o}})}{\partial \phi _{o}}$$  $$-k(\rho _{s}R_{os}^2{\rm sin}(\phi _{o} - \phi _{s}) - 3\rho _{o}\rho _{s}{\rm sin}(\phi _{o} - \phi _{s})(\rho _{o} - \rho _{s}{\rm cos}(\phi _{o} - \phi _{s})))/R_{os}^5$$  $$\frac{\partial (dg_{os,\rho _{o}})}{\partial Z_o}$$  $$3k(Z_o - Z_{s})(\rho _{o} - \rho _{s}{\rm cos}(\phi _{o} - \phi _{s}))/R_{os}^5$$  $$\frac{\partial (dg_{os,\phi _{o}})}{\partial \rho _{o}}$$  $$3k\rho _{s}{\rm sin}(\phi _{o} - \phi _{s})(\rho _{o} - \rho _{s}{\rm cos}(\phi _{o} - \phi _{s}))/R_{os}^5$$  $$\frac{\partial (dg_{os,\phi _{o}})}{\partial \phi _{o}}$$  $$-k\rho _{s}(R_{os}^2{\rm cos}(\phi _{o} - \phi _{s}) -3\rho _{o}\rho _{s}({\rm sin}(\phi _{o} - \phi _{s}))^2)/R_{os}^5$$  $$\frac{\partial (dg_{os,\phi _{o}})}{\partial Z_o}$$  $$3k\rho _{s}(Z_o - Z_{s}){\rm sin}(\phi _{o} - \phi _{s})/R_{os}^5$$  $$\frac{\partial (dg_{os,Z_o})}{\partial \rho _{o}}$$  $$3k(Z_o - Z_{s})(\rho _{o} - \rho _{s}{\rm cos}(\phi _{o} - \phi _{s}))/R_{os}^5$$  $$\frac{\partial (dg_{os,Z_o})}{\partial \phi _{o}}$$  $$3k\rho _{o}\rho _{s}(Z_o - Z_{s}){\rm sin}(\phi _{o} - \phi _{s})/R_{os}^5$$  $$\frac{\partial (dg_{os,Z_o})}{\partial Z_o}$$  $$-k(R_{os}^2 - 3(Z_o - Z_{s})^2)/R_{os}^5$$  Partial differential  Formula  $$\frac{\partial (dg_{os,\rho _{o}})}{\partial \rho _{o}}$$  $$-k(R_{os}^2 - 3(\rho _{o} - \rho _{s}{\rm cos}(\phi _{o} - \phi _{s}))^2)/R_{os}^5$$  $$\frac{\partial (dg_{os,\rho _{o}})}{\partial \phi _{o}}$$  $$-k(\rho _{s}R_{os}^2{\rm sin}(\phi _{o} - \phi _{s}) - 3\rho _{o}\rho _{s}{\rm sin}(\phi _{o} - \phi _{s})(\rho _{o} - \rho _{s}{\rm cos}(\phi _{o} - \phi _{s})))/R_{os}^5$$  $$\frac{\partial (dg_{os,\rho _{o}})}{\partial Z_o}$$  $$3k(Z_o - Z_{s})(\rho _{o} - \rho _{s}{\rm cos}(\phi _{o} - \phi _{s}))/R_{os}^5$$  $$\frac{\partial (dg_{os,\phi _{o}})}{\partial \rho _{o}}$$  $$3k\rho _{s}{\rm sin}(\phi _{o} - \phi _{s})(\rho _{o} - \rho _{s}{\rm cos}(\phi _{o} - \phi _{s}))/R_{os}^5$$  $$\frac{\partial (dg_{os,\phi _{o}})}{\partial \phi _{o}}$$  $$-k\rho _{s}(R_{os}^2{\rm cos}(\phi _{o} - \phi _{s}) -3\rho _{o}\rho _{s}({\rm sin}(\phi _{o} - \phi _{s}))^2)/R_{os}^5$$  $$\frac{\partial (dg_{os,\phi _{o}})}{\partial Z_o}$$  $$3k\rho _{s}(Z_o - Z_{s}){\rm sin}(\phi _{o} - \phi _{s})/R_{os}^5$$  $$\frac{\partial (dg_{os,Z_o})}{\partial \rho _{o}}$$  $$3k(Z_o - Z_{s})(\rho _{o} - \rho _{s}{\rm cos}(\phi _{o} - \phi _{s}))/R_{os}^5$$  $$\frac{\partial (dg_{os,Z_o})}{\partial \phi _{o}}$$  $$3k\rho _{o}\rho _{s}(Z_o - Z_{s}){\rm sin}(\phi _{o} - \phi _{s})/R_{os}^5$$  $$\frac{\partial (dg_{os,Z_o})}{\partial Z_o}$$  $$-k(R_{os}^2 - 3(Z_o - Z_{s})^2)/R_{os}^5$$  View Large Analytical solutions for the triple integrals in eqs (3), (4) and (6) clearly are complicated relative to their least-squares numerical approximations by GLQ integration (e.g. Stroud & Secrest 1966; Ku 1977; von Frese et al. 1981) as given by   \begin{eqnarray} &&{ \int _{Z_{s} = Z_{1}}^{Z_{2}} \int _{\phi _{s} = \phi _{1}}^{\phi _{2}} \int _{\rho _{s} = \rho _{1}}^{\rho _{2}} [f(\rho _{o}, \phi _{o}, Z_o, \rho _{s}, \phi _{s}, Z_{s})] d\rho _{s} d\phi _{s} dZ_{s} }\nonumber\\ && \cong\,\frac{(\rho _{2} - \rho _{1}) (\phi _{2} - \phi _{1}) (Z_{2} - Z_{1})}{8} \sum _{nk = 1}^{K} \sum _{nj = 1}^{J} \sum _{ni = 1}^{I}\nonumber\\ && [f(\rho _{o}, \phi _{o}, Z_o, {\hat{\rho }_{ni}}, {\hat{\phi }_{nj}}, {\hat{Z}_{nk}})] A_{ni} A_{nj} A_{nk} \end{eqnarray} (7)Here, the square bracketed f-functions are generalizations of the appropriate integrands in observation and source point coordinates, and the Gauss–Legendre coefficients Ani, Anj and Ank correspond to the ρni, ϕnj and Znk coordinates of the nth Gaussian node within the interval (−1, 1). To transform the integration from the unit interval into the $${\hat{\rho }_{ni}}$$, $${\hat{\phi }_{nj}}$$and $${\hat{Z}_{nk}}$$ coordinates within the source body requires scaling the body point coordinates as   $${{\hat{\rho }_{ni}} = \frac{\rho _{ni} (\rho _{2} - \rho _{1}) + (\rho _{2} + \rho _{1})}{2},}$$ (8a)  $${{\hat{\phi }_{nj}} = \frac{\phi _{nj} (\phi _{2} - {\phi }_{1}) + ({\phi }_{2} + {\phi }_{1})}{2},}$$ (8b)and   $${{\hat{Z}_{nk}} = \frac{Z_{nk} (Z_{2} - Z_{1}) + (Z_{2} + Z_{1})}{2}.}$$ (8c) Thus, eq. (7) computes the least-squares gravitational effects of the cylindrical prism by summing at each observation point the GLQ-coefficient weighted gravitational effects of I × J × K equivalent point poles. The accuracies of the GLQ estimates in eq. (7) are directly controlled by the selected number of (I, J, K) nodes that, in turn, directly controls the computational time. In practice, effective accuracy is achieved by selecting the node spacing to be smaller than the distance between the source nodes and the observation point (e.g. Ku 1977; von Frese et al. 1981; Li et al. 2011). The cylindrical coordinate system is clearly advantageous for evaluating volume integrals in determining the gravitational effects of vertical tubular sources. However, to compare these effects with observed anomaly values, they must be projected into the coordinate system of the observations (e.g. Rim & Li 2016). Thus, the respective cylindrical-to-Cartesian projections of the gravitational vector and tensor components at the observation point are   $$\left( \begin{array} {c}g_{o,X_{o}} \\ {g_{o,Y_{o}}} \\ {g_{o,Z_o}} \\ \end{array} \right) = P(\phi _o) \left( \begin{array}{c}g_{o,\rho _{o}} \\ {g_{o,\phi _{o}}} \\ {g_{o,Z_o}} \\ \end{array} \right),$$ (9)and   \begin{eqnarray} &&{\left( \begin{array}{c{@}{\quad}c{@}{\quad}c}g_{o,X_{o},X_{o}} & g_{o,X_{o},Y_{o}} & g_{o,X_{o},Z_o} \\ g_{o,Y_{o},X_{o}} & g_{o,Y_{o},Y_{o}} & g_{o,Y_{o},Z_o} \\ g_{o,Z_o,X_{o}} & g_{o,Z_o,Y_{o}} & g_{o,Z_o,Z_o} \\ \end{array} \right)}\nonumber\\ &&= P(\phi _o) \left( \begin{array}{c{@}{\quad}c{@}{\quad}c}g_{o,\rho _{o},\rho _{o}} & g_{o,\rho _{o},\phi _{o}} & g_{o,\rho _{o},Z_o} \\ g_{o,\phi _{o},\rho _{o}} & g_{o,\phi _{o},\phi _{o}} & g_{o,\phi _{o},Z_o} \\ g_{o,Z_o,\rho _{o}} & g_{o,Z_o,\phi _{o}} & g_{o,Z_o,Z_o} \\ \end{array} \right) P^T(\phi _o), \end{eqnarray} (10)where PT(ϕo) is the transpose of:   $$P(\phi _o) = \left( \begin{array}{c{@}{\quad}c{@}{\quad}c}{\rm cos} (\phi _o) & -{\rm sin} (\phi _o) & 0 \\ {}{\rm sin} (\phi _o) & {\rm cos} (\phi _o) & 0 \\ 0 & 0 & 1 \\ \end{array} \right).$$ (11) 3 SENSITIVITY ANALYSIS AND FORWARD MODELING This section tests the GLQ-based $$g_{o,Z_o}$$ estimates against the uniform density, vertical cylinder modeling from the other methods described above, illustrates and compares the GLQ-based gravitational estimates of a uniform density, vertical cylinder against the gravitational estimates of the same cylinder from Rim & Li (2016), and illustrates the comprehensive GLQ gravitational effects of the cylindrical prism element. Fig. 2 shows how the tested gravitational effects were modeled using the observation-centred convention for a vertical circular cylinder of constant density σs, length L and radius R with respective depths Z1 and Z2 to its bottom and top surfaces along the XZ-plane passing through the cylinder’s axis. Figs 3–6 give axis-centred XZ-profile comparisons of the $$g_{o,Z_o}$$ anomaly estimates for a vertical cylinder from the GLQ method, CC method (Kwok 1991), EI method (Nabighian 1962), FH method (Singh 1977), LI method (Van Baak 1989), LP method (Telford et al. 1990), P method (Parasnis 1961) plus SA and VL methods (Nettleton 1942). Here, the vertical cylinder was assumed with density σs = 2800 kg m−3, length L = 500 m, top surfaces Z2 at -50, -500, -2500 and -5000 m below the observation plane, and radii R of 50 m (Fig. 3), 500 m (Fig. 4), 2500 m (Fig. 5) and 5000 m (Fig. 6). Figure 2. View largeDownload slide Vertical circular cylinder of constant density σs, length L and radius R with the respective bottom and top surfaces at depths Z1 and Z2 and the XZ-plane passing through the cylinder’s axis. Figure 2. View largeDownload slide Vertical circular cylinder of constant density σs, length L and radius R with the respective bottom and top surfaces at depths Z1 and Z2 and the XZ-plane passing through the cylinder’s axis. Figure 3. View largeDownload slide Gravitational anomaly estimates of a vertical circular cylinder of density σs = 2800 kg m−3, length L = 500 m, radius R = 50 m and top surface at Z2 from the GLQ, CC, EI, FH, LI, LP, P, SA and VL methods. Vertical axis is in logarithmic units of $$g_{o,Z_o}$$ (mGal). Figure 3. View largeDownload slide Gravitational anomaly estimates of a vertical circular cylinder of density σs = 2800 kg m−3, length L = 500 m, radius R = 50 m and top surface at Z2 from the GLQ, CC, EI, FH, LI, LP, P, SA and VL methods. Vertical axis is in logarithmic units of $$g_{o,Z_o}$$ (mGal). Figure 4. View largeDownload slide Gravitational anomaly profile comparisons modeled using the methods listed in Fig. 3 for a vertical circular cylinder of density σs = 2800 kg m−3, length L = 500 m, radius R = 500 m and top surface at Z2. Vertical axis is in logarithmic units of $$g_{o,Z_o}$$ (mGal). Figure 4. View largeDownload slide Gravitational anomaly profile comparisons modeled using the methods listed in Fig. 3 for a vertical circular cylinder of density σs = 2800 kg m−3, length L = 500 m, radius R = 500 m and top surface at Z2. Vertical axis is in logarithmic units of $$g_{o,Z_o}$$ (mGal). Figure 5. View largeDownload slide Gravitational anomaly profile comparisons modeled using the methods listed in Fig. 3 for a vertical circular cylinder of density σs = 2800 kg m−3, length L = 500 m, radius R = 2500 m and top surface at Z2. Vertical axis is in logarithmic units of $$g_{o,Z_o}$$ (mGal). Figure 5. View largeDownload slide Gravitational anomaly profile comparisons modeled using the methods listed in Fig. 3 for a vertical circular cylinder of density σs = 2800 kg m−3, length L = 500 m, radius R = 2500 m and top surface at Z2. Vertical axis is in logarithmic units of $$g_{o,Z_o}$$ (mGal). Figure 6. View largeDownload slide Gravitational anomaly profile comparisons modeled using the methods listed in Fig. 3 for a vertical circular cylinder of density σs = 2800 kg m−3, length L = 500 m, radius R = 5000 m and top surface at Z2. Vertical axis is in logarithmic units of $$g_{o,Z_o}$$ (mGal). Figure 6. View largeDownload slide Gravitational anomaly profile comparisons modeled using the methods listed in Fig. 3 for a vertical circular cylinder of density σs = 2800 kg m−3, length L = 500 m, radius R = 5000 m and top surface at Z2. Vertical axis is in logarithmic units of $$g_{o,Z_o}$$ (mGal). To better visualize the modeled anomaly ranges and details near the cylinder’s axis, the $$g_{o,Z_o}$$ estimates in mGal are displayed in base-10 logarithmic values. For the GLQ modeling, the cylinder’s radial dimensions were filled using 2 (R = 50 m), 11 (R = 500 m), 51 (R = 2500 m) and 101 (R = 5000 m) nodes, whereas the azimuthal dimensions were filled using 360 (R = 50 m), 360 (R = 500 m), 360 (R = 2500 m) and 720 (R = 5000 m) nodes. The length dimensions of all 16 cylinders were filled using 2 nodes. Study of Figs 3(e), 4(e) and (f), 5(e)–(g), 6(e)–(h) reveals that the $$g_{o,Z_o}$$ estimates from the LP, SA and VL methods increasingly differ from the GLQ, CC, EI, FH, LI and P estimates towards the central axis. These differences however, decrease as the R (radii) and Z2 (depths) of the cylinder become smaller and deeper. Fig. 7 compares the rms differences between the $$g_{o,Z_o}$$ estimates from the GLQ modeling and those of the other eight methods described above. To remove the effects of singularities, the gravitational anomaly estimates at X = 0 were deleted from the rms difference calculations. Figs 7(a)–(d) suggest that for constant L (lengths) and R (radii), the rms differences become smaller as the cylinder’s Z2 (depth) increases. Likewise, increasing rms differences result for constant L (lengths) and Z2 (depths) as the R (radii) of the cylinder increases as shown in Figs 7(e)–(h). Fig. 7 also suggests that for all R and Z2 values, the CC, EI and FH estimates are very similar. Figure 7. View largeDownload slide Comparison of $$g_{o,Z_o}$$ rms differences between the GLQ method and the other eight methods described in Section 3. Subfigures (a)–(d) illustrate for a constant radius of R, how changes in the depth to the top surface Z2 affects the rms differences, whereas subfigures (e)–(h) show for a constant depth to the top surface of Z2, how changes in the radius R affects the rms differences. Figure 7. View largeDownload slide Comparison of $$g_{o,Z_o}$$ rms differences between the GLQ method and the other eight methods described in Section 3. Subfigures (a)–(d) illustrate for a constant radius of R, how changes in the depth to the top surface Z2 affects the rms differences, whereas subfigures (e)–(h) show for a constant depth to the top surface of Z2, how changes in the radius R affects the rms differences. Figs 8 and 9 illustrate at 0 m altitude, the cylindrical coordinate gravitational anomaly components for a uniformly dense, vertical cylinder with σs = 1000 kg m−3 and dimensions ρ1 = 0 m to ρ2 = 100 m and Z1 = −350 m to Z2 = −50 m. These effects were GLQ-modeled using 4 × 360 × 2 = 2880 nodes spanning the radial, longitudinal and vertical dimensions of the cylinder’s volume. Figs 8(a)–(d), respectively, show the modeled gravitational Po-potential and vector gravitational $$g_{o,\rho _{o}}$$-, $$g_{o,\phi _{o}}$$- and $$g_{o,Z_o}$$- anomaly components, whereas Fig. 9 gives the related gravitational gradient anomaly components from the three trace and three unique off-diagonal elements of the symmetric matrix in eq. (5). Note that four of the nine elements of this matrix are redundant because $$g_{o,\rho _{o},\phi _{o}} = g_{o,\phi _{o},\rho _{o}}$$, $$g_{o,\rho _{o},Z_o} = g_{o,Z_o,\rho _{o}}$$, and $$g_{o,\phi _{o},Z_o} = g_{o,Z_o,\phi _{o}}$$ and any two components of $$g_{o,\rho _{o},\rho _{o}}, g_{o,\phi _{o},\phi _{o}}$$ and $$g_{o,Z_o,Z_o}$$ may be related to the third one via Laplace’s equation. Within working precision, these results match the estimates from the closed-form gravitational expressions based on complete Lipschitz–Hankel type integral solutions involving elliptic integrals of the first and second kinds and Heumann’s lambda function (Rim & Li 2016). Figure 8. View largeDownload slide Maps (a)–(d) give the respective gravitational potential (Po) and vector radial ($$g_{o,\rho _o}$$), longitudinal ($$g_{o,\phi _o}$$) and vertical ($$g_{o,Z_o}$$) anomaly components of a cylinder with ρ1 = 0 m, ρ2 = 100 m, Z1 = −350 m, and Z2 = −50 m with density σs = 1000 kg m−3 in cylindrical coordinates. The gravitational effects were modeled at 0 m altitude for the cylinder using I = 4, J = 360 and K = 2 nodes. See Table 3 for statistical attributes on these maps. Figure 8. View largeDownload slide Maps (a)–(d) give the respective gravitational potential (Po) and vector radial ($$g_{o,\rho _o}$$), longitudinal ($$g_{o,\phi _o}$$) and vertical ($$g_{o,Z_o}$$) anomaly components of a cylinder with ρ1 = 0 m, ρ2 = 100 m, Z1 = −350 m, and Z2 = −50 m with density σs = 1000 kg m−3 in cylindrical coordinates. The gravitational effects were modeled at 0 m altitude for the cylinder using I = 4, J = 360 and K = 2 nodes. See Table 3 for statistical attributes on these maps. Figure 9. View largeDownload slide Cylindrical coordinate gravitational gradient components for the cylinder in Fig. 8. See Table 3 for statistical attributes on these maps. Figure 9. View largeDownload slide Cylindrical coordinate gravitational gradient components for the cylinder in Fig. 8. See Table 3 for statistical attributes on these maps. Figs 10 and 11 illustrate at 500 m altitude, the cylindrical coordinate gravitational effects of a vertical cylindrical volume element or prism with σs = 2800 kg m−3 and dimensions ρ1 = 3000 m to ρ2 = 4000 m, ϕ1 = 30° to ϕ2 = 60° and Z1 = −10000 m to Z2 = −5000 m. These effects were GLQ modeled using 2 × 2 × 2 = 8 nodes with each nodal pair spanning the radial, longitudinal and vertical dimensions of the prism’s volume. Figs 10(a)–(d), respectively, show the modeled gravitational Po-potential and vector gravitational $$g_{o,\rho _{o}}$$-, $$g_{o,\phi _{o}}$$-, and $$g_{o,Z_o}$$- anomaly components, whereas Fig. 11 gives the gravitational gradient anomaly components from the three trace and three unique off-diagonal elements of the symmetric matrix in eq. (5). The gravitational effects of the cylindrical prism in Figs 10 and 11 are, respectively, mapped in Cartesian coordinates in Figs 12 and 13 via the transformation eqs (9)–(11). In general, the Cartesian mapping is more typical of exploration geophysical applications even though the data in cylindrical and Cartesian projections provide equivalent spatial anomaly detail. Figure 10. View largeDownload slide Maps (a)–(d) give the respective gravitational potential (Po) and vector gravitational radial ($$g_{o,\rho _o}$$), longitudinal ($$g_{o,\phi _o}$$) and vertical ($$g_{o,Z_o}$$) anomaly components of a cylindrical prism with ρ1 = 3000 m, ρ2 = 4000 m, ϕ1 = 30°, ϕ2 = 60°, Z1 = −10  000 m and Z2 = −5000 m with density σs = 2800 kg m−3 in cylindrical coordinates. The gravitational effects were modeled at 500 m altitude for the cylindrical prism using I = 2, J = 2 and K = 2 nodes. See Table 3 for statistical attributes on these maps. Figure 10. View largeDownload slide Maps (a)–(d) give the respective gravitational potential (Po) and vector gravitational radial ($$g_{o,\rho _o}$$), longitudinal ($$g_{o,\phi _o}$$) and vertical ($$g_{o,Z_o}$$) anomaly components of a cylindrical prism with ρ1 = 3000 m, ρ2 = 4000 m, ϕ1 = 30°, ϕ2 = 60°, Z1 = −10  000 m and Z2 = −5000 m with density σs = 2800 kg m−3 in cylindrical coordinates. The gravitational effects were modeled at 500 m altitude for the cylindrical prism using I = 2, J = 2 and K = 2 nodes. See Table 3 for statistical attributes on these maps. Figure 11. View largeDownload slide Cylindrical coordinate gravitational gradient components for the cylindrical prism in Fig. 10. See Table 3 for statistical attributes on these maps. Figure 11. View largeDownload slide Cylindrical coordinate gravitational gradient components for the cylindrical prism in Fig. 10. See Table 3 for statistical attributes on these maps. Figure 12. View largeDownload slide Maps (a)–(d) give the respective gravitational potential (Po) and vector left–right ($$g_{o,X_o}$$), top–bottom ($$g_{o,Y_o}$$) and vertical ($$g_{o,Z_o}$$) anomaly components for the cylindrical prism in Fig. 10 in Cartesian coordinates. Note that anomaly maps of Figs 10(a) and 12(a), as well as Figs 10(d) and 12(d) are the same for both cylindrical and Cartesian coordinate systems. See Table 3 for statistical attributes on these maps. Figure 12. View largeDownload slide Maps (a)–(d) give the respective gravitational potential (Po) and vector left–right ($$g_{o,X_o}$$), top–bottom ($$g_{o,Y_o}$$) and vertical ($$g_{o,Z_o}$$) anomaly components for the cylindrical prism in Fig. 10 in Cartesian coordinates. Note that anomaly maps of Figs 10(a) and 12(a), as well as Figs 10(d) and 12(d) are the same for both cylindrical and Cartesian coordinate systems. See Table 3 for statistical attributes on these maps. Figure 13. View largeDownload slide Gravitational gradient components for the cylindrical prism in Fig. 10 in Cartesian coordinates. Note that anomaly maps of Figs 11(f) and 13(f) are the same for both cylindrical and Cartesian coordinate systems. See Table 3 for statistical attributes on these maps. Figure 13. View largeDownload slide Gravitational gradient components for the cylindrical prism in Fig. 10 in Cartesian coordinates. Note that anomaly maps of Figs 11(f) and 13(f) are the same for both cylindrical and Cartesian coordinate systems. See Table 3 for statistical attributes on these maps. 4 GLQ GRAVITATIONAL MODELING OF A FORMATION WASHOUT Formation washouts are sources of uncertainty in borehole logging. They may form whenever the drill bit penetrates a weak rock formation (e.g. a loosely consolidated sandstone layer), and may extend partially or completely around the borehole wall. In this section, procedures for locating formation washouts and modeling their gravitational effects by GLQ integration are considered. The well-established method based on Compton scattering of gamma rays uses a continuous beam of photons typically emitted from a 0.66 MeV Cs137 source into the formation (Dewan 1983). These photons interact with the electrons of formation material to scatter in all directions with intensities proportional to the electron densities that can be related to the bulk densities of the formation. Modern implementation of this method allows for determination of azimuthal mass density variations (Bargach et al. 2000), which can ultimately be used to locate any formation washouts (density anomalies) around the borehole wall. Another approach for mapping formation washouts around a borehole may be based on borehole gravimetry and related gradiometry (Nekut 1989). Here, a grid of point gravity poles or masses is setup around the circumference of the wellbore. Iterative gravitational modeling is then adopted to update the masses of the gridded point sources until their gravitational effects fit the gravitational observations to within prescribed error limits. The resulting distribution of gridded point masses may finally be used to help map out washouts around the borehole wall. In theory, eqs (3), (4) and (6) can be used for estimating the washout’s gravitational effects at any point within the borehole on and outside the washout’s surface envelop. As examples, Figs 14 and 15 present a washout’s Cartesian coordinates gravitational effects along the respective central XZ- and YZ-planes of a borehole with radius ρBorehole = 0.5 m. The washout extends from ρ1 = 0.5 m to ρ2 = 1 m, ϕ1 = 150° to ϕ2 = 210° and Z1 = −1000.25 m to Z2 = −999.75 m along the respective radial, azimuthal and axial dimensions of the borehole. The assumed density contrast between the washout and its confining formation is σs = −1500 kg m−3. These effects were initially GLQ modeled in cylindrical coordinates via eqs (3), (4) and (6) using 5 × 21 × 20 = 2100 nodes, and subsequently projected into Cartesian coordinates via eqs (9)–(11). These figures indicate that $$g_{o,Y_{o}}$$, $$g_{o,X_{o},Y_{o}} = g_{o,Y_{o},X_{o}}$$ and $$g_{o,Y_{o},Z_{o}} = g_{o,Z_{o},Y_{o}}$$ show little-to-no variation across the XZ-plane in contrast to the other components that vary strongly across this plane, whereas all components vary strongly across the YZ-plane. Clearly, to use the inversion of borehole gravity data to map a formation washout, the related anomaly variations must be observed on a plane within the sensitivity ranges of modern gravimeter and gravity gradiometer sensors (Chen & Cook 2005). Figure 14. View largeDownload slide Maps (a)–(d) give in Cartesian coordinates the respective gravitational Po potential and $$g_{o,X_o}$$, $$g_{o,Y_o}$$ and $$g_{o,Z_o}$$ vector anomaly components in the respective X-, Y-, and Z-directions of a formation washout with ρ1 = 0.5 m, ρ2 = 1 m, ϕ1 = 150°, ϕ2 = 210°, Z1 = −1000.25 m and Z2 = −999.75 m with density σs = −1500 kg m−3. The gravitational effects were modeled along the wellbore’s central XZ- and YZ-planes using I = 5, J = 21 and K = 20 nodes. See Table 3 for statistical attributes on these maps. Figure 14. View largeDownload slide Maps (a)–(d) give in Cartesian coordinates the respective gravitational Po potential and $$g_{o,X_o}$$, $$g_{o,Y_o}$$ and $$g_{o,Z_o}$$ vector anomaly components in the respective X-, Y-, and Z-directions of a formation washout with ρ1 = 0.5 m, ρ2 = 1 m, ϕ1 = 150°, ϕ2 = 210°, Z1 = −1000.25 m and Z2 = −999.75 m with density σs = −1500 kg m−3. The gravitational effects were modeled along the wellbore’s central XZ- and YZ-planes using I = 5, J = 21 and K = 20 nodes. See Table 3 for statistical attributes on these maps. Figure 15. View largeDownload slide Gravitational gradient components for the formation washout in Fig. 14 in Cartesian coordinates. See Table 3 for statistical attributes on these maps. Figure 15. View largeDownload slide Gravitational gradient components for the formation washout in Fig. 14 in Cartesian coordinates. See Table 3 for statistical attributes on these maps. 5 GLQ GRAVITATIONAL MODELING OF A SALT DOME As another application of the GLQ gravitational modeling method, the salt intrusion in Fig. 16 is considered that was mapped within the sedimentary units of the Laurentian Basin off Canada’s eastern coast by gravity and seismic data. The Late Triassic breakup of the North American and African continents (Uchupi & Austin Jr 1979) triggered the formation of a number of basins along the Scotian passive margin including the Laurentian Basin (Adam & Krezsek 2012). This passive margin is known to be associated with evaporate deposition, as well as halokinesis in some areas (Uchupi & Austin Jr 1979; Adam & Krezsek 2012). Figure 16. View largeDownload slide (a) Geographic location of the Laurentian Basin study area involving seismic lines STP-3, STP-4 and STP-21. (b) Processed seismic reflection section for line STP-21. (c) Seismic meta-attribute section for line STP-21, where the dashed white polygon outlines the inferred salt body. Figure 16. View largeDownload slide (a) Geographic location of the Laurentian Basin study area involving seismic lines STP-3, STP-4 and STP-21. (b) Processed seismic reflection section for line STP-21. (c) Seismic meta-attribute section for line STP-21, where the dashed white polygon outlines the inferred salt body. The study region includes the horseshoe-shaped negative gravitational anomaly at 500 m above mean sea level shown in Fig. 17(a) from shipborne survey-constrained free-air gravitational anomaly (FAGA) estimates of the Eigen-6c4 model (Förste et al. 2014). Elements of this anomaly were imaged in 1984–1985 by 29 lines of 2-D multichannel seismic reflection surveying that collected some 3100 km of data. The survey obtained roughly 7000 ms of raw data at a 4-ms sampling rate using source and receiver intervals of 25 m each. The data were reprocessed in 2006 and are now accessible at the Department of Natural Resources Canada. The south-central section of the horseshoe-shaped gravitational minimum includes a circular feature that sits north of the NE-SW trending seismic line STP-21, and between the NW-SE trending seismic lines STP-3 and STP-4 (Fig. 17a). Figure 17. View largeDownload slide (a) Free-air gravitational anomaly (FAGA) estimates from the Eigen-6c4 model, (b) digital elevation model (DEM), (c) GLQ-modeled terrain gravitational effects (TGE) and (d) the complete Bouguer anomalies (CBA) for the study area. The FAGA, TGE and CBA estimates were evaluated at 500 m above sea level. See Table 3 for statistical attributes on these maps. Figure 17. View largeDownload slide (a) Free-air gravitational anomaly (FAGA) estimates from the Eigen-6c4 model, (b) digital elevation model (DEM), (c) GLQ-modeled terrain gravitational effects (TGE) and (d) the complete Bouguer anomalies (CBA) for the study area. The FAGA, TGE and CBA estimates were evaluated at 500 m above sea level. See Table 3 for statistical attributes on these maps. The seismic data from line STP-21 along the southern margin of the circular negative gravitational anomaly reveal a possible salt body for the anomaly’s source (Fig. 16b). To enhance the intrusive body’s edges, seismic pattern recognition analysis (Hashemi et al. 2008) was carried out on the line STP-21 data. A set of seismic data samples within the salt body was selected to compare against another data set sampled from the remaining sedimentary formations within the Basin. The comparison facilitated calculating the relevant seismic attributes in a non-redundant feature space (Hashemi & Javaherian 2009) that, in turn, were subjected to principal component analysis for a supervised multilayer perception (MLP) neural network classifier function (Hashemi et al. 2008). The MLP mapping function was applied to all traces of the STP-21 line to infer the spatial extent of the salt body shown in Fig. 16(c). This analysis estimated the body’s average width, average two-way-time (TWT) thickness and the average TWT between its top surface and the seabed at about 7700 m (R = 3850 m), 5283 ms and 173 ms, respectively. The thickness of the seawater column above the salt body amounts to roughly TWT = 480 ms, so that for assumed seawater, sediment and salt intrusion compressional velocities of VP,Seawater = 1484 m s−1, VP,Sediment = 1600 m s−1 and VP,Salt = 4450 m s−1, the depth below sea level of the seabed, and the body’s top surface below the seabed and the salt body’s thickness amount to roughly 356, 138 and 11755 m, respectively. Assuming an average salt density of σSalt = 2170 kg m−3, the cylindrical body’s vertical-component gravitational anomaly ($$g_{o,Z_o}$$) at 500 m above sea level was GLQ modeled as shown by the red curve in Fig. 18(a). The GLQ gravitational modeling used 8 × 360 × 2 = 5760 nodes in the respective radial, azimuthal and depth dimensions. Using GLQ-spherical prism gravitational modeling (Asgharzadeh et al. 2007), the terrain gravitational effects (TGE) of the digital elevation model (DEM, Smith & Sandwell 1997) in Fig. 17(b) were calculated at 500 m above sea level (Fig. 17c) assuming average sea water, sediment, and basement rock densities of σSeawater = 1030 kg m−3, σSediment = 2250 kg m−3 and σBasement = 3030 kg m−3, respectively. The TGE, in turn, were subtracted from the FAGA to obtain the region’s complete Bouguer gravitational anomalies (CBA) in Fig. 17(d). In Fig. 18(a), the magenta and green curves illustrate the CBA profile and its linear regional trend along line STP-21. The linear trend was subtracted from the CBA as possible effects due to lack of isostatic compensation and other regional modeling errors within the basement topography and intrusive body. The removal of the linear regional yields a close fit between the residual CBA (blue curve) and the modeled $$g_{o,Z_o}$$ effects of the salt intrusion (Fig. 18b). Accordingly, the choices for sea water, sediment, salt and basement density, velocity and depth variations yield gravity signals that are consistent with the region’s DEM, seismic reflection and FAGA observations. Figure 18. View largeDownload slide (a) CBA (magenta curve) with its linear regional trend (green curve), and the superimposed GLQ modeled $$g_{o,Z_o}$$ gravitational anomaly of the seismically constrained intrusive salt body (red curve) along seismic line STP-21. (b) Residual CBA (blue curve) with the superimposed GLQ modeled $$g_{o,Z_o}$$ gravitational anomaly of the intrusive salt body (red curve) along the STP-21 seismic line. Figure 18. View largeDownload slide (a) CBA (magenta curve) with its linear regional trend (green curve), and the superimposed GLQ modeled $$g_{o,Z_o}$$ gravitational anomaly of the seismically constrained intrusive salt body (red curve) along seismic line STP-21. (b) Residual CBA (blue curve) with the superimposed GLQ modeled $$g_{o,Z_o}$$ gravitational anomaly of the intrusive salt body (red curve) along the STP-21 seismic line. 6 CONCLUSIONS AND RECOMMENDATIONS The complete gravitational effects (i.e. potential, anomaly and gradient components) of the vertical cylinder with radial-, azimuthal- and axial-density variations were developed in the context of the least-squares numerical gravitational effects of the cylindrical unit volume element or prism estimated by GLQ integration. The GLQ modeling was developed to accommodate gravity studies of rock formations with heterogeneous densities due to physical or chemical alterations. Comparisons against the results of modeling simple uniform density vertical cylinders by a number of conventional methods illustrated the veracity of the GLQ approach. As an application in borehole gravity, the gravitational effects of a formation washout were estimated to assist in determining azimuthal variations of formation bulk densities around the borehole wall. A further application considered the gravitational effects of a seismically and gravimetrically imaged salt body within the Laurentian Basin for constraining the velocity, density and geometric properties of the Basin’s sedimentary formations. Table 3 lists the statistical attributes for Figs 8–15 and 17 to facilitate the efforts of readers to replicate these results. Table 3. Statistical attributes for the maps in Figs 8–15 and 17 include the map altitude (Zo), as well as amplitude minimum (Min), maximum (Max), average (Mean) and standard deviation (Std). The amplitude units (AU) and the power (n) of the amplitude multiplier (× 10n) are also listed. Figure  Map  Zo(m)  Min  Max  Mean  Std  n  AU  8a  Po  0  43.9899  343.8151  92.7064  48.1505  0  mGal × m  8b  $$g_{o,\rho _o}$$  0  −0.7800  −0.0304  −0.1467  0.1478  0  mGal  8c  $$g_{o,\phi _o}$$  0  0  0  0  0  0  mGal  8d  $$g_{o,Z_o}$$  0  −1.9163  −0.0043  −0.0806  0.1977  0  mGal  9a  $$g_{o,\rho _o,\rho _o}$$  0  −10.1395  2.6059  0.3666  0.8223  −3  mGal m−1  9b  $$g_{o,\rho _o,\phi _o} = g_{o,\phi _o,\rho _o}$$  0  0  0  0  0  −3  mGal m−1  9c  $$g_{o,\rho _o,Z_o} = g_{o,Z_o,\rho _o}$$  0  0.0088  9.8842  0.5130  1.4259  −3  mGal m−1  9d  $$g_{o,\phi _o,\phi _o}$$  0  −10.1970  −0.0215  −0.4492  1.1290  −3  mGal m−1  9e  $$g_{o,\phi _o,Z_o} = g_{o,Z_o,\phi _o}$$  0  0  0  0  0  −3  mGal m−1  9f  $$g_{o,Z_o,Z_o}$$  0  −0.3725  20.3365  0.0826  1.5095  −3  mGal m−1  10a  Po  500  6563.1  22051  11623  3453.9  0  mGal × m  10b  $$g_{o,\rho _o}$$  500  −1.1071  0.9918  −0.5909  0.2327  0  mGal  10c  $$g_{o,\phi _o}$$  500  −1.0101  1.0101  0  0.2282  0  mGal  10d  $$g_{o,Z_o}$$  500  −2.9206  −0.0765  −0.5441  0.5504  0  mGal  11a  $$g_{o,\rho _o,\rho _o}$$  500  −0.3990  0.0752  0.0260  0.0688  −3  mGal m−1  11b  $$g_{o,\rho _o,\phi _o} = g_{o,\phi _o,\rho _o}$$  500  −0.1059  0.1059  0  0.0335  −3  mGal m−1  11c  $$g_{o,\rho _o,Z_o} = g_{o,Z_o,\rho _o}$$  500  −0.3353  0.3363  0.0795  0.0764  −3  mGal m−1  11d  $$g_{o,\phi _o,\phi _o}$$  500  −0.3941  −0.0092  −0.0612  0.0613  −3  mGal m−1  11e  $$g_{o,\phi _o,Z_o} = g_{o,Z_o,\phi _o}$$  500  −0.3360  0.3360  0  0.0600  −3  mGal m−1  11f  $$g_{o,Z_o,Z_o}$$  500  −0.0125  0.7930  0.0352  0.1208  −3  mGal m−1  12a  Po  500  6563.1  22051  11623  3453.9  0  mGal × m  12b  $$g_{o,X_o}$$  500  −1.1065  1.1066  0  0.4772  0  mGal  12c  $$g_{o,Y_o}$$  500  −1.1065  1.1066  0  0.4772  0  mGal  12d  $$g_{o,Z_o}$$  500  −2.9206  −0.0765  −0.5441  0.5504  0  mGal  13a  $$g_{o,X_o,X_o}$$  500  −0.3965  0.0752  −0.0176  0.0736  −3  mGal m−1  13b  $$g_{o,X_o,Y_o} = g_{o,Y_o,X_o}$$  500  −0.1084  0.1078  0  0.0430  −3  mGal m−1  13c  $$g_{o,X_o,Z_o} = g_{o,Z_o,X_o}$$  500  −0.3356  0.3355  0  0.0888  −3  mGal m−1  13d  $$g_{o,Y_o,Y_o}$$  500  −0.3965  0.0752  −0.0176  0.0736  −3  mGal m−1  13e  $$g_{o,Y_o,Z_o} = g_{o,Z_o,Y_o}$$  500  −0.3356  0.3355  0  0.0888  −3  mGal m−1  13f  $$g_{o,Z_o,Z_o}$$  500  −0.0125  0.7930  0.0352  0.1208  −3  mGal m−1  14a  Po  −  −0.0050  −0.0009  −0.0018  0.0006  0  mGal × m  14b  $$g_{o,X_o}$$  −  0.0001  0.0107  0.0013  0.0013  0  mGal  14c  $$g_{o,Y_o}$$  −  −0.0012  0.0012  0  0.0003  0  mGal  14d  $$g_{o,Z_o}$$  −  −0.0051  0.0051  0  0.0011  0  mGal  15a  $$g_{o,X_o,X_o}$$  −  −34.8223  1.9380  −1.8250  3.8013  −3  mGal m−1  15b  $$g_{o,X_o,Y_o} = g_{o,Y_o,X_o}$$  −  −3.4642  3.4642  0  0.7612  −3  mGal m−1  15c  $$g_{o,X_o,Z_o} = g_{o,Z_o,X_o}$$  −  −25.9132  25.9132  0  3.0304  −3  mGal m−1  15d  $$g_{o,Y_o,Y_o}$$  −  0.2082  10.7654  1.5485  1.6250  −3  mGal m−1  15e  $$g_{o,Y_o,Z_o} = g_{o,Z_o,Y_o}$$  −  −1.1407  1.1407  0  0.3227  −3  mGal m−1  15f  $$g_{o,Z_o,Z_o}$$  −  −7.5209  24.0788  0.2765  2.5012  −3  mGal m−1  17a  FAGA  500  −36.2590  32.6437  −3.4772  13.2437  0  mGal  17b  DEM  −  −3250  −3  −681  787  0  m  17c  TGE  500  −76.7469  −0.6997  −17.4067  19.6542  0  mGal  17d  CBA  500  −27.1880  61.7900  13.7376  19.6185  0  mGal  Figure  Map  Zo(m)  Min  Max  Mean  Std  n  AU  8a  Po  0  43.9899  343.8151  92.7064  48.1505  0  mGal × m  8b  $$g_{o,\rho _o}$$  0  −0.7800  −0.0304  −0.1467  0.1478  0  mGal  8c  $$g_{o,\phi _o}$$  0  0  0  0  0  0  mGal  8d  $$g_{o,Z_o}$$  0  −1.9163  −0.0043  −0.0806  0.1977  0  mGal  9a  $$g_{o,\rho _o,\rho _o}$$  0  −10.1395  2.6059  0.3666  0.8223  −3  mGal m−1  9b  $$g_{o,\rho _o,\phi _o} = g_{o,\phi _o,\rho _o}$$  0  0  0  0  0  −3  mGal m−1  9c  $$g_{o,\rho _o,Z_o} = g_{o,Z_o,\rho _o}$$  0  0.0088  9.8842  0.5130  1.4259  −3  mGal m−1  9d  $$g_{o,\phi _o,\phi _o}$$  0  −10.1970  −0.0215  −0.4492  1.1290  −3  mGal m−1  9e  $$g_{o,\phi _o,Z_o} = g_{o,Z_o,\phi _o}$$  0  0  0  0  0  −3  mGal m−1  9f  $$g_{o,Z_o,Z_o}$$  0  −0.3725  20.3365  0.0826  1.5095  −3  mGal m−1  10a  Po  500  6563.1  22051  11623  3453.9  0  mGal × m  10b  $$g_{o,\rho _o}$$  500  −1.1071  0.9918  −0.5909  0.2327  0  mGal  10c  $$g_{o,\phi _o}$$  500  −1.0101  1.0101  0  0.2282  0  mGal  10d  $$g_{o,Z_o}$$  500  −2.9206  −0.0765  −0.5441  0.5504  0  mGal  11a  $$g_{o,\rho _o,\rho _o}$$  500  −0.3990  0.0752  0.0260  0.0688  −3  mGal m−1  11b  $$g_{o,\rho _o,\phi _o} = g_{o,\phi _o,\rho _o}$$  500  −0.1059  0.1059  0  0.0335  −3  mGal m−1  11c  $$g_{o,\rho _o,Z_o} = g_{o,Z_o,\rho _o}$$  500  −0.3353  0.3363  0.0795  0.0764  −3  mGal m−1  11d  $$g_{o,\phi _o,\phi _o}$$  500  −0.3941  −0.0092  −0.0612  0.0613  −3  mGal m−1  11e  $$g_{o,\phi _o,Z_o} = g_{o,Z_o,\phi _o}$$  500  −0.3360  0.3360  0  0.0600  −3  mGal m−1  11f  $$g_{o,Z_o,Z_o}$$  500  −0.0125  0.7930  0.0352  0.1208  −3  mGal m−1  12a  Po  500  6563.1  22051  11623  3453.9  0  mGal × m  12b  $$g_{o,X_o}$$  500  −1.1065  1.1066  0  0.4772  0  mGal  12c  $$g_{o,Y_o}$$  500  −1.1065  1.1066  0  0.4772  0  mGal  12d  $$g_{o,Z_o}$$  500  −2.9206  −0.0765  −0.5441  0.5504  0  mGal  13a  $$g_{o,X_o,X_o}$$  500  −0.3965  0.0752  −0.0176  0.0736  −3  mGal m−1  13b  $$g_{o,X_o,Y_o} = g_{o,Y_o,X_o}$$  500  −0.1084  0.1078  0  0.0430  −3  mGal m−1  13c  $$g_{o,X_o,Z_o} = g_{o,Z_o,X_o}$$  500  −0.3356  0.3355  0  0.0888  −3  mGal m−1  13d  $$g_{o,Y_o,Y_o}$$  500  −0.3965  0.0752  −0.0176  0.0736  −3  mGal m−1  13e  $$g_{o,Y_o,Z_o} = g_{o,Z_o,Y_o}$$  500  −0.3356  0.3355  0  0.0888  −3  mGal m−1  13f  $$g_{o,Z_o,Z_o}$$  500  −0.0125  0.7930  0.0352  0.1208  −3  mGal m−1  14a  Po  −  −0.0050  −0.0009  −0.0018  0.0006  0  mGal × m  14b  $$g_{o,X_o}$$  −  0.0001  0.0107  0.0013  0.0013  0  mGal  14c  $$g_{o,Y_o}$$  −  −0.0012  0.0012  0  0.0003  0  mGal  14d  $$g_{o,Z_o}$$  −  −0.0051  0.0051  0  0.0011  0  mGal  15a  $$g_{o,X_o,X_o}$$  −  −34.8223  1.9380  −1.8250  3.8013  −3  mGal m−1  15b  $$g_{o,X_o,Y_o} = g_{o,Y_o,X_o}$$  −  −3.4642  3.4642  0  0.7612  −3  mGal m−1  15c  $$g_{o,X_o,Z_o} = g_{o,Z_o,X_o}$$  −  −25.9132  25.9132  0  3.0304  −3  mGal m−1  15d  $$g_{o,Y_o,Y_o}$$  −  0.2082  10.7654  1.5485  1.6250  −3  mGal m−1  15e  $$g_{o,Y_o,Z_o} = g_{o,Z_o,Y_o}$$  −  −1.1407  1.1407  0  0.3227  −3  mGal m−1  15f  $$g_{o,Z_o,Z_o}$$  −  −7.5209  24.0788  0.2765  2.5012  −3  mGal m−1  17a  FAGA  500  −36.2590  32.6437  −3.4772  13.2437  0  mGal  17b  DEM  −  −3250  −3  −681  787  0  m  17c  TGE  500  −76.7469  −0.6997  −17.4067  19.6542  0  mGal  17d  CBA  500  −27.1880  61.7900  13.7376  19.6185  0  mGal  View Large Further exploitation of the GLQ modeling results is possible by adapting them for modeling the gravitational effects of arbitrarily tilted cylinders (Rim & Li 2016). By Poisson’s theorem (Hinze et al. 2013), these results may also be extended to model the magnetic and thermal effects of arbitrarily oriented tubular geological bodies. ACKNOWLEDGEMENTS We thank the National Elite Foundation of Iran and the Institute of Geophysics at the University of Tehran for supporting this work. We also thank the anonymous reviewers for their constructive comments that helped to improve the manuscript. REFERENCES Adam J., Krezsek C., 2012. Basin-scale salt tectonic processes of the Laurentian Basin, Eastern Canada: insights from integrated regional 2D seismic interpretation an 4D physical experiments, Geol. Soc. Lond. Spec. Publ. , 363, 331– 360. Google Scholar CrossRef Search ADS   Asgharzadeh M.F., von Frese R.R.B., Kim H., Leftwich T.E., Kim J.W., 2007. Spherical prism gravity effects by Gauss-Legendre quadrature integration, Geophys. J. Int. , 169( 1), 1– 11. Google Scholar CrossRef Search ADS   Bargach S., Bornemann T., Codazzi D., Ford G., Grether B., Rohler H., 2000. Real-time LWD: logging for drilling, Oilfield Rev. , 12( 3), 58– 78. Brown A.R., Lautzenhiser T.V., 1982. The effect of dipping beds on a borehole gravimeter survey, Geophysics , 47( 1), 25– 30. Google Scholar CrossRef Search ADS   Chen Y.T., Cook A., 2005. Gravitational Experiments in the Laboratory , Cambridge University Press, Cambridge. Dewan J.T., 1983. Essentials of Modern Open-hole Log Interpretation , PennWell Publishing, Tusla, Oklahoma. Förste C. et al.  , 2014. EIGEN-6C4 The latest combined global gravity field model including GOCE data up to degree and order 2190 of GFZ Potsdam and GRGS Toulouse, GFZ Data Services, doi:10.5880/icgem.2015.1. Hashemi H., Javaherian A., 2009. Seismic attribute redundancy reduction using statistical feature extraction technique, in 1st International Petroleum Conference and Exhibition, EAGE, Shiraz, Iran , doi: 10.3997/2214-4609.20145883. Hashemi H., Tax D.M.J., Duin R.P.W., Javaherian A., de Groot P., 2008. Gas chimney detection based on improving the performance of combined multilayer perceptron and support vector classifier, Nonlinear Process. Geophys. , 15, 863– 871. Google Scholar CrossRef Search ADS   Hinze W.J., von Frese R. R.B., Saad A.H., 2013. Gravity and Magnetic Exploration: Principles, Practices, and Applications , Cambridge University Press, Cambridge. Google Scholar CrossRef Search ADS   Jenkins A.J.O., Messfin D., Moon W., 1983. Gravity modelling of salt domes and pinnacle reefs, J. Can. Soc. Explor. Geophys. , 19( 1), 51– 56. Ku C.C., 1977. A direct computation of gravity and magnetic anomalies caused by 2- and 3-dimensional bodies of arbitrary shape and arbitrary magnetic polarization by equivalent point method and a simplified cubic spline, Geophysics , 42, 610– 622. Google Scholar CrossRef Search ADS   Kwok Y.-K., 1991. Singularities in gravity computation for vertical cylinders and prisms, Geophys. J. Int. , 104, 1– 10. Google Scholar CrossRef Search ADS   Li Z., Hao T., Xu Y., Xu U., 2011. An efficient and adaptive approach for modeling gravity effects in spherical coordinates, J. appl. Geophys. , 73, 221– 231. Google Scholar CrossRef Search ADS   Nabighian M.N., 1962. The gravitational attraction of a right vertical circular cylinder at points external to it, Geofis. pura appl. , 53( 1), 45– 51. Google Scholar CrossRef Search ADS   Nagy D., 1966. The gravitational attraction of a right rectangular prism, Geophysics , 31( 2), 362– 371. Google Scholar CrossRef Search ADS   Nekut A.G., 1989. Borehole gravity gradiometry, Geophysics , 54( 2), 225– 234. Google Scholar CrossRef Search ADS   Nettleton L.L., 1942. Gravity and magnetic calculations, Geophysics , 7( 3), 293– 310. Google Scholar CrossRef Search ADS   Parasnis D.S., 1961. Exact expressions for the gravitational attraction of a circular lamina at all points of space and of a right circular vertical cylinder at points external to it, Geophys. Prospect. , 9( 3), 382– 398. Google Scholar CrossRef Search ADS   Reilly W.I., 1969. Gravitational and magnetic effects of a right circular cylinder, N. Z. J. Geol. Geophys. , 12( 2–3), 497– 506. Google Scholar CrossRef Search ADS   Rim H., Li Y., 2012. Single-hole imaging using borehole gravity gradiometry, Geophysics , 77( 5), G67– G76. Google Scholar CrossRef Search ADS   Rim H., Li Y., 2016. Gravity gradient tensor due to a cylinder, Geophysics , 81( 4), 59– 66. Google Scholar CrossRef Search ADS   Roussel C., Verdun J., Cali J., Masson F., 2015. Complete gravity field of an ellipsoidal prism by Gauss-Legendre quadrature, Geophys. J. Int. , 203( 3), 2220– 2236. Google Scholar CrossRef Search ADS   Singh S.K., 1977. Gravitational attraction of a vertical right circular cylinder, Geophys. J. R. astr. Soc , 50, 243– 246. Google Scholar CrossRef Search ADS   Smith N.J., 1950. The case for gravity data from boreholes, Geophysics , 15( 4), 605– 636. Google Scholar CrossRef Search ADS   Smith W. H.F., Sandwell D.T., 1997. Global seafloor topography from satellite altimetry and ship depth soundings, Science , 277, 1957– 1962. Stroud A.H., Secrest D., 1966. Gaussian Quadrature Formulas , Prentice-Hall, Upper Saddle River, NJ. Telford W.M., Geldart L.P., Sheriff R.E., 1990. Applied Geophysics , Cambridge University Press, Cambridge. Google Scholar CrossRef Search ADS   Uchupi E., Austin Jr J.A., 1979. The geologic history of the passive margin off New England and the Canadian maritime provinces, Tectonophysics , 59, 53– 69. Google Scholar CrossRef Search ADS   Van Baak D.A., 1989. Efficient computation of the gravitational field of a generalized cylindrical prism, Geophysics , 54( 3), 402– 405. Google Scholar CrossRef Search ADS   von Frese R. R.B., Hinze W.J., Braile L.W., Luca A.J., 1981. Spherical earth gravity and magnetic anomaly modeling by Gauss-Legendre quadrature integration, J. Geophys. , 49, 234– 242. © The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# Comprehensive gravitational modeling of the vertical cylindrical prism by Gauss–Legendre quadrature integration

, Volume 212 (1) – Jan 1, 2018
21 pages

/lp/ou_press/comprehensive-gravitational-modeling-of-the-vertical-cylindrical-prism-6BSx2pBqEy
Publisher
Oxford University Press
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx413
Publisher site
See Article on Publisher Site

### Abstract

SUMMARY Forward modeling is the basis of gravitational anomaly inversion that is widely applied to map subsurface mass variations. This study uses numerical least-squares Gauss–Legendre quadrature (GLQ) integration to evaluate the gravitational potential, anomaly and gradient components of the vertical cylindrical prism element. These results, in turn, may be integrated to accurately model the complete gravitational effects of fluid bearing rock formations and other vertical cylinder-like geological bodies with arbitrary variations in shape and density. Comparing the GLQ gravitational effects of uniform density, vertical circular cylinders against the effects calculated by a number of other methods illustrates the veracity of the GLQ modeling method and the accuracy limitations of the other methods. Geological examples include modeling the gravitational effects of a formation washout to help map azimuthal variations of the formation’s bulk densities around the borehole wall. As another application, the gravitational effects of a seismically and gravimetrically imaged salt dome within the Laurentian Basin are evaluated for the velocity, density and geometric properties of the Basin’s sedimentary formations. Geopotential theory, Gravity anomalies and Earth structure, Downhole methods, Numerical modelling, Statistical seismology 1 INTRODUCTION Calculating the gravitational effects (i.e. the gravitational potential and its higher order spatial gradients) of an extended mass body requires integrating at the observation point the differential gravitational effects of the unit volume elements making up the body. Any volume symmetries of the source considerably simplify the related scalar, vector and tensor integrals. Thus, for example, the gravitational effects of the prism’s mass are readily accessible in Cartesian (Nagy 1966), spherical (Asgharzadeh et al. 2007) and elliptical (Roussel et al. 2015) coordinates. The considerable interest in computing the gravitational effects of vertical cylinder-like bodies has spawned a relatively large number of quantitative modeling approaches, which generally consider the case of either the finite or semi-infinite vertical axis. Nettleton (1942), for example, presented a closed formula for the vertical-component gravitational attraction ($$g_{o,Z_o}$$) of the vertical semi-infinite cylinder (VL) where the cylinder’s radius is considerably smaller than its thickness. Nettleton (1942) also presented the closed formula for $$g_{o,Z_o}$$ based on the solid angle (SA) subtended at the observation point O by the vertical circular cylinder with thickness that is less than or equal to about a half of its average depth. Van Baak (1989) calculated $$g_{o,Z_o}$$ of the vertical cylinder using line integrals (LIs) around the boundaries of the volume-filling discs. The cross-sectional integral of each disc is obtained from the bounding curve either in closed analytical form if it is a complete circle or numerically as the sum of the n line segment integrals if an approximating polygon with n line vertices is applicable. The line segment integrals for all discs are evaluated and summed at each observation point with an accuracy that increases as the number of line segments increases in the bounding curve approximations of the discs. Additional closed $$g_{o,Z_o}$$ formulae are available for the circular vertical cylinder (P) using Cartesian coordinate integrals (Parasnis 1961), as well as elliptic integrals (EIs) of the first and second kind and Heuman’s Lambda function (Nabighian 1962) and the Fourier–Hankel (FH) transform (Singh 1977). Rim & Li (2016) presented closed expressions of $$g_{o,Z_o}$$ and its gradients for the vertical and tilted cylinder. Numerical approximations of $$g_{o,Z_o}$$ are also available involving the first few terms of the Legendre polynomial (LP) expansion (e.g. Parasnis 1961; Reilly 1969; Jenkins et al. 1983; Telford et al. 1990) and the conjugate complex (CC) variables formulation of Green’s Theorem (Kwok 1991). In general, the above gravitational effect estimates typically assume a constant density throughout the cylinder, although Jenkins et al. (1983) present a method to accommodate density variations as annuli and discs along the respective radial and axial dimensions. However, they do not accommodate azimuthal density variations within the cylinder, which may appear, for example, along borehole walls through physically or chemically altered fluid-bearing rock formations. These alterations occur in various sizes to complicate borehole gravimetry applications for interpreting surface gravity, seismic and borehole geology observations (e.g. Smith 1950). Borehole gravity data also constrain bedding dip estimates assuming the formation densities are available (Brown & Lautzenhiser 1982). For dipping beds, Nekut (1989) showed that borehole gravity gradiometry also can be used to estimate the bedding dips, strikes and densities, as well as the average density of the hosting rock formation. Modern applications for borehole gravity gradiometry also facilitate CO2 sequestration, enhanced oil recovery and mining and exploration applications (Rim & Li 2012). In this paper, the complete gravitational effects of the vertical cylinder with discrete variations of density along the radial, azimuthal and axial dimensions are developed in the context of the cylindrical unit volume element or prism shown in Fig. 1. These effects are numerically calculated from the Gauss–Legendre quadrature (GLQ) distribution of point gravity poles within the unit prism elements that fill up the cylinder (Stroud & Secrest 1966). These equations are applicable for calculating the complete gravitational effects of salt domes, granitic and other igneous plume-like intrusions, carbonate reefs, gas-hydrate columns, kimberlite pipes, crustal meteorite impact structures, mantle plumes, formation washouts and dipping rock formations along wellbores, and other homogeneous or heterogeneous vertical partial or full cylinder-like geological features. The accuracy of the GLQ-based gravitational modeling is tested by comparing the $$g_{o,Z_o}$$ estimates against those of the other methods described above that vary with changes in the cylinder’s radius, thickness and burial depth. Figure 1. View largeDownload slide The cylindrical ρϕZ coordinate system for modeling at the observation point O(ρo, ϕo, Zo), the gravitational effects of a vertical cylindrical prism by integrating the differential source mass dms through the volume of the prism with bottom and top surfaces that include the corner points A(ρ1, ϕ1, Z1) and B(ρ2, ϕ2, Z2). The Cartesian perspective is given by the superposed, co-registered Cartesian XYZ coordinate system axes. Figure 1. View largeDownload slide The cylindrical ρϕZ coordinate system for modeling at the observation point O(ρo, ϕo, Zo), the gravitational effects of a vertical cylindrical prism by integrating the differential source mass dms through the volume of the prism with bottom and top surfaces that include the corner points A(ρ1, ϕ1, Z1) and B(ρ2, ϕ2, Z2). The Cartesian perspective is given by the superposed, co-registered Cartesian XYZ coordinate system axes. To help demonstrate the geological utility of the GLQ method, the gravitational effects of a formation washout are estimated in the context of the related borehole gravity measurements. The gravitational effects of a seismically constrained salt dome also are modeled within the Laurentian Basin off the eastern coast of Canada. Insights on the Basin’s average sediment density, velocity, and geometric properties are obtained from comparing the modeled gravitational anomaly to shipborne-constrained free-air gravitational anomaly estimates from the Eigen-6c4 spherical harmonic coefficient expansion (Förste et al. 2014). 2 GRAVITATIONAL EFFECTS FOR THE CYLINDRICAL PRISM Assuming (XYZ)-Cartesian and (ρϕZ)-cylindrical coordinate systems, with the vertical Z-axis upwards positive from origins on the orthogonal (XY)- and (ρϕ)-planes (Fig. 1), the far-field gravitational effects of a differential vertical cylindrical prism element in both coordinate systems is inversely proportional to the distance Ros between the observation (O) and source (S) points. This distance in cylindrical coordinates is   $${R_{os} = \left[{\rho _{o}^2 + \rho _{s}^2 - 2 \rho _{o} \rho _{s} {\rm cos}(\phi _{o} - \phi _{s}) + (Z_o - Z_{s})^2}\right]^\frac{1}{2},}$$ (1)where ρo and ρs are the radial distances of the respective observation and source points from the Z-axis, with respective ϕo- and ϕs-longitude coordinates, at the distances Zo and Zs from the origin’s XY- and ρϕ-planes. In addition, the differential volume dvs of the cylindrical prism’s differential mass dms is   $${dv_{s} = \rho _{s} d\rho _{s} d\phi _{s} dZ_{s}.}$$ (2) Thus, the total gravitational potential Po of the cylindrical prism at the observation point is the sum of the differential mass potentials given there by   $${P_{o} = \int _{V} dP_{os} = \int _{Z_{s} = Z_{1}}^{Z_{2}} \int _{\phi _{s} = \phi _{1}}^{\phi _{2}} \int _{\rho _{s} = \rho _{1}}^{\rho _{2}} \left[\frac{k}{R_{os}}\right] d\rho _{s} d\phi _{s} dZ_{s},}$$ (3)where k = Gσsρs, G = 6.67408 × 10−11 m3kg−1s−2 is the universal gravitational constant and σs is the mass density of the cylindrical prism. In addition, the spatial gradients of Po at the observation point in the ρo-radial, ϕo-longitudinal, and Zo-axial directions yield the vector gravitational field components ($$g_{{o,\rho _{o}}}$$, $$g_{{o,\phi _{o}}}$$, $$g_{{o,Z_o}}$$) given by   \begin{eqnarray} g_{{o,\rho _{o}}} &=& \int _{V} dg_{os,\rho _{o}} = \int _{Z_{s} = Z_{1}}^{Z_{2}} \int _{\phi _{s} = \phi _{1}}^{\phi _{2}} \int _{\rho _{s} = \rho _{1}}^{\rho _{2}} \nonumber\\ &&\left[-k \frac{(\rho _{o} - \rho _{s} {\rm cos} (\phi _{o} - \phi _{s}))}{R_{os}^3}\right] d\rho _{s} d\phi _{s} dZ_{s}, \end{eqnarray} (4a)  \begin{eqnarray} g_{{o,\phi _{o}}} &=& \int _{V} dg_{os,\phi _{o}} = \int _{Z_{s} = Z_{1}}^{Z_{2}} \int _{\phi _{s} = \phi _{1}}^{\phi _{2}} \int _{\rho _{s} = \rho _{1}}^{\rho _{2}} \nonumber\\ &&\left[-k \frac{\rho _{s} {\rm sin}(\phi _{o} - \phi _{s})}{R_{os}^3}\right] d\rho _{s} d\phi _{s} dZ_{s}, \end{eqnarray} (4b)and   \begin{eqnarray} g_{{o,Z_o}} &=& \int _{V} dg_{os,Z_o} = \int _{Z_{s} = Z_{1}}^{Z_{2}} \int _{\phi _{s} = \phi _{1}}^{\phi _{2}} \int _{\rho _{s} = \rho _{1}}^{\rho _{2}} \nonumber\\ &&\left[-k \frac{(Z_o - Z_{s})}{R_{os}^3}\right] d\rho _{s} d\phi _{s} dZ_{s}. \end{eqnarray} (4c) Furthermore, the spatial gradients of the gravitational field components of the cylindrical prism involve the elements of the 3 × 3 matrix given by   $$\nabla dg_{os} = \left( \begin{array}{c{@}{\quad}c{@}{\quad}c}dg_{os,\rho _{o},\rho _{o}} & {dg_{os,\rho _{o},\phi _{o}}} & {dg_{os,\rho _{o},Z_o}} \\ {dg_{os,\phi _{o},\rho _{o}}} & {dg_{os,\phi _{o},\phi _{o}}} & {dg_{os,\phi _{o},Z_o}} \\ {dg_{os,Z_o,\rho _{o}}} & {dg_{os,Z_o,\phi _{o}}} & {dg_{os,Z_o,Z_o}} \\ \end{array} \right)\!,$$ (5)where $$dg_{os,\rho _{o},\rho _{o}}$$ is the spatial gradient of $$dg_{os,\rho _{o}}$$ in the ρo-radial direction, $$dg_{os,\rho _{o},\phi _{o}}$$ is the spatial gradient of $$dg_{os,\rho _{o}}$$ in the ϕo-longitudinal direction and so on. These elements are formulated in Table 1 with the related partial differentials listed in Table 2. In general, the gravitational gradient tensor components at the observation point may be expressed over the source region by   $${g_{o,m,u} = \int _{V} dg_{os,m,u} = \int _{Z_{s} = Z_{1}}^{Z_{2}} \int _{\phi _{s} = \phi _{1}}^{\phi _{2}} \int _{\rho _{s} = \rho _{1}}^{\rho _{2}} dg_{os,m,u},}$$ (6)where m is the ρo-, ϕo-, or Zo-components of go and dgos, and u is the ρo-, ϕo-, or Zo-directions along which the gradient is taken. Table 1. Formulae for the gravitational tensor gradient components of the matrix in eq. (5). Component  [Formula]· dρsdϕsdZs  $$dg_{os,\rho _{o},\rho _{o}}$$  $$\frac{\partial \left(dg_{os,\rho _{o}}\right)}{\partial \rho _{o}}$$  $$dg_{os,\rho _{o},\phi _{o}}$$  $$\frac{1}{\rho _{o}} \left(\frac{\partial (dg_{os,\rho _{o}})}{\partial \phi _{o}} - dg_{os,\phi _{o}}\right)$$  $$dg_{os,\rho _{o},Z_o}$$  $$\frac{\partial (dg_{os,\rho _{o}})}{\partial Z_o}$$  $$dg_{os,\phi _{o},\rho _{o}}$$  $$\frac{\partial (dg_{os, \phi _{o}})}{\partial \rho _{o}}$$  $$dg_{os,\phi _{o},\phi _{o}}$$  $$\frac{1}{\rho _{o}} \left(\frac{\partial (dg_{os, \phi _{o}})}{\partial \phi _{o}} + dg_{os, \rho _{o}}\right)$$  $$dg_{os,\phi _{o},Z_o}$$  $$\frac{\partial (dg_{os,\phi _{o}})}{\partial Z_o}$$  $$dg_{os,Z_o,\rho _{o}}$$  $$\frac{\partial (dg_{os, Z_o})}{\partial \rho _{o}}$$  $$dg_{os,Z_o,\phi _{o}}$$  $$\frac{1}{\rho _{o}} \frac{\partial (dg_{os, Z_o})}{\partial \phi _{o}}$$  $$dg_{os,Z_o,Z_o}$$  $$\frac{\partial (dg_{os,Z_o})}{\partial Z_o}$$  Component  [Formula]· dρsdϕsdZs  $$dg_{os,\rho _{o},\rho _{o}}$$  $$\frac{\partial \left(dg_{os,\rho _{o}}\right)}{\partial \rho _{o}}$$  $$dg_{os,\rho _{o},\phi _{o}}$$  $$\frac{1}{\rho _{o}} \left(\frac{\partial (dg_{os,\rho _{o}})}{\partial \phi _{o}} - dg_{os,\phi _{o}}\right)$$  $$dg_{os,\rho _{o},Z_o}$$  $$\frac{\partial (dg_{os,\rho _{o}})}{\partial Z_o}$$  $$dg_{os,\phi _{o},\rho _{o}}$$  $$\frac{\partial (dg_{os, \phi _{o}})}{\partial \rho _{o}}$$  $$dg_{os,\phi _{o},\phi _{o}}$$  $$\frac{1}{\rho _{o}} \left(\frac{\partial (dg_{os, \phi _{o}})}{\partial \phi _{o}} + dg_{os, \rho _{o}}\right)$$  $$dg_{os,\phi _{o},Z_o}$$  $$\frac{\partial (dg_{os,\phi _{o}})}{\partial Z_o}$$  $$dg_{os,Z_o,\rho _{o}}$$  $$\frac{\partial (dg_{os, Z_o})}{\partial \rho _{o}}$$  $$dg_{os,Z_o,\phi _{o}}$$  $$\frac{1}{\rho _{o}} \frac{\partial (dg_{os, Z_o})}{\partial \phi _{o}}$$  $$dg_{os,Z_o,Z_o}$$  $$\frac{\partial (dg_{os,Z_o})}{\partial Z_o}$$  View Large Table 2. Formulae for partial differentials in Table 1 with k = Gσsρs. Partial differential  Formula  $$\frac{\partial (dg_{os,\rho _{o}})}{\partial \rho _{o}}$$  $$-k(R_{os}^2 - 3(\rho _{o} - \rho _{s}{\rm cos}(\phi _{o} - \phi _{s}))^2)/R_{os}^5$$  $$\frac{\partial (dg_{os,\rho _{o}})}{\partial \phi _{o}}$$  $$-k(\rho _{s}R_{os}^2{\rm sin}(\phi _{o} - \phi _{s}) - 3\rho _{o}\rho _{s}{\rm sin}(\phi _{o} - \phi _{s})(\rho _{o} - \rho _{s}{\rm cos}(\phi _{o} - \phi _{s})))/R_{os}^5$$  $$\frac{\partial (dg_{os,\rho _{o}})}{\partial Z_o}$$  $$3k(Z_o - Z_{s})(\rho _{o} - \rho _{s}{\rm cos}(\phi _{o} - \phi _{s}))/R_{os}^5$$  $$\frac{\partial (dg_{os,\phi _{o}})}{\partial \rho _{o}}$$  $$3k\rho _{s}{\rm sin}(\phi _{o} - \phi _{s})(\rho _{o} - \rho _{s}{\rm cos}(\phi _{o} - \phi _{s}))/R_{os}^5$$  $$\frac{\partial (dg_{os,\phi _{o}})}{\partial \phi _{o}}$$  $$-k\rho _{s}(R_{os}^2{\rm cos}(\phi _{o} - \phi _{s}) -3\rho _{o}\rho _{s}({\rm sin}(\phi _{o} - \phi _{s}))^2)/R_{os}^5$$  $$\frac{\partial (dg_{os,\phi _{o}})}{\partial Z_o}$$  $$3k\rho _{s}(Z_o - Z_{s}){\rm sin}(\phi _{o} - \phi _{s})/R_{os}^5$$  $$\frac{\partial (dg_{os,Z_o})}{\partial \rho _{o}}$$  $$3k(Z_o - Z_{s})(\rho _{o} - \rho _{s}{\rm cos}(\phi _{o} - \phi _{s}))/R_{os}^5$$  $$\frac{\partial (dg_{os,Z_o})}{\partial \phi _{o}}$$  $$3k\rho _{o}\rho _{s}(Z_o - Z_{s}){\rm sin}(\phi _{o} - \phi _{s})/R_{os}^5$$  $$\frac{\partial (dg_{os,Z_o})}{\partial Z_o}$$  $$-k(R_{os}^2 - 3(Z_o - Z_{s})^2)/R_{os}^5$$  Partial differential  Formula  $$\frac{\partial (dg_{os,\rho _{o}})}{\partial \rho _{o}}$$  $$-k(R_{os}^2 - 3(\rho _{o} - \rho _{s}{\rm cos}(\phi _{o} - \phi _{s}))^2)/R_{os}^5$$  $$\frac{\partial (dg_{os,\rho _{o}})}{\partial \phi _{o}}$$  $$-k(\rho _{s}R_{os}^2{\rm sin}(\phi _{o} - \phi _{s}) - 3\rho _{o}\rho _{s}{\rm sin}(\phi _{o} - \phi _{s})(\rho _{o} - \rho _{s}{\rm cos}(\phi _{o} - \phi _{s})))/R_{os}^5$$  $$\frac{\partial (dg_{os,\rho _{o}})}{\partial Z_o}$$  $$3k(Z_o - Z_{s})(\rho _{o} - \rho _{s}{\rm cos}(\phi _{o} - \phi _{s}))/R_{os}^5$$  $$\frac{\partial (dg_{os,\phi _{o}})}{\partial \rho _{o}}$$  $$3k\rho _{s}{\rm sin}(\phi _{o} - \phi _{s})(\rho _{o} - \rho _{s}{\rm cos}(\phi _{o} - \phi _{s}))/R_{os}^5$$  $$\frac{\partial (dg_{os,\phi _{o}})}{\partial \phi _{o}}$$  $$-k\rho _{s}(R_{os}^2{\rm cos}(\phi _{o} - \phi _{s}) -3\rho _{o}\rho _{s}({\rm sin}(\phi _{o} - \phi _{s}))^2)/R_{os}^5$$  $$\frac{\partial (dg_{os,\phi _{o}})}{\partial Z_o}$$  $$3k\rho _{s}(Z_o - Z_{s}){\rm sin}(\phi _{o} - \phi _{s})/R_{os}^5$$  $$\frac{\partial (dg_{os,Z_o})}{\partial \rho _{o}}$$  $$3k(Z_o - Z_{s})(\rho _{o} - \rho _{s}{\rm cos}(\phi _{o} - \phi _{s}))/R_{os}^5$$  $$\frac{\partial (dg_{os,Z_o})}{\partial \phi _{o}}$$  $$3k\rho _{o}\rho _{s}(Z_o - Z_{s}){\rm sin}(\phi _{o} - \phi _{s})/R_{os}^5$$  $$\frac{\partial (dg_{os,Z_o})}{\partial Z_o}$$  $$-k(R_{os}^2 - 3(Z_o - Z_{s})^2)/R_{os}^5$$  View Large Analytical solutions for the triple integrals in eqs (3), (4) and (6) clearly are complicated relative to their least-squares numerical approximations by GLQ integration (e.g. Stroud & Secrest 1966; Ku 1977; von Frese et al. 1981) as given by   \begin{eqnarray} &&{ \int _{Z_{s} = Z_{1}}^{Z_{2}} \int _{\phi _{s} = \phi _{1}}^{\phi _{2}} \int _{\rho _{s} = \rho _{1}}^{\rho _{2}} [f(\rho _{o}, \phi _{o}, Z_o, \rho _{s}, \phi _{s}, Z_{s})] d\rho _{s} d\phi _{s} dZ_{s} }\nonumber\\ && \cong\,\frac{(\rho _{2} - \rho _{1}) (\phi _{2} - \phi _{1}) (Z_{2} - Z_{1})}{8} \sum _{nk = 1}^{K} \sum _{nj = 1}^{J} \sum _{ni = 1}^{I}\nonumber\\ && [f(\rho _{o}, \phi _{o}, Z_o, {\hat{\rho }_{ni}}, {\hat{\phi }_{nj}}, {\hat{Z}_{nk}})] A_{ni} A_{nj} A_{nk} \end{eqnarray} (7)Here, the square bracketed f-functions are generalizations of the appropriate integrands in observation and source point coordinates, and the Gauss–Legendre coefficients Ani, Anj and Ank correspond to the ρni, ϕnj and Znk coordinates of the nth Gaussian node within the interval (−1, 1). To transform the integration from the unit interval into the $${\hat{\rho }_{ni}}$$, $${\hat{\phi }_{nj}}$$and $${\hat{Z}_{nk}}$$ coordinates within the source body requires scaling the body point coordinates as   $${{\hat{\rho }_{ni}} = \frac{\rho _{ni} (\rho _{2} - \rho _{1}) + (\rho _{2} + \rho _{1})}{2},}$$ (8a)  $${{\hat{\phi }_{nj}} = \frac{\phi _{nj} (\phi _{2} - {\phi }_{1}) + ({\phi }_{2} + {\phi }_{1})}{2},}$$ (8b)and   $${{\hat{Z}_{nk}} = \frac{Z_{nk} (Z_{2} - Z_{1}) + (Z_{2} + Z_{1})}{2}.}$$ (8c) Thus, eq. (7) computes the least-squares gravitational effects of the cylindrical prism by summing at each observation point the GLQ-coefficient weighted gravitational effects of I × J × K equivalent point poles. The accuracies of the GLQ estimates in eq. (7) are directly controlled by the selected number of (I, J, K) nodes that, in turn, directly controls the computational time. In practice, effective accuracy is achieved by selecting the node spacing to be smaller than the distance between the source nodes and the observation point (e.g. Ku 1977; von Frese et al. 1981; Li et al. 2011). The cylindrical coordinate system is clearly advantageous for evaluating volume integrals in determining the gravitational effects of vertical tubular sources. However, to compare these effects with observed anomaly values, they must be projected into the coordinate system of the observations (e.g. Rim & Li 2016). Thus, the respective cylindrical-to-Cartesian projections of the gravitational vector and tensor components at the observation point are   $$\left( \begin{array} {c}g_{o,X_{o}} \\ {g_{o,Y_{o}}} \\ {g_{o,Z_o}} \\ \end{array} \right) = P(\phi _o) \left( \begin{array}{c}g_{o,\rho _{o}} \\ {g_{o,\phi _{o}}} \\ {g_{o,Z_o}} \\ \end{array} \right),$$ (9)and   \begin{eqnarray} &&{\left( \begin{array}{c{@}{\quad}c{@}{\quad}c}g_{o,X_{o},X_{o}} & g_{o,X_{o},Y_{o}} & g_{o,X_{o},Z_o} \\ g_{o,Y_{o},X_{o}} & g_{o,Y_{o},Y_{o}} & g_{o,Y_{o},Z_o} \\ g_{o,Z_o,X_{o}} & g_{o,Z_o,Y_{o}} & g_{o,Z_o,Z_o} \\ \end{array} \right)}\nonumber\\ &&= P(\phi _o) \left( \begin{array}{c{@}{\quad}c{@}{\quad}c}g_{o,\rho _{o},\rho _{o}} & g_{o,\rho _{o},\phi _{o}} & g_{o,\rho _{o},Z_o} \\ g_{o,\phi _{o},\rho _{o}} & g_{o,\phi _{o},\phi _{o}} & g_{o,\phi _{o},Z_o} \\ g_{o,Z_o,\rho _{o}} & g_{o,Z_o,\phi _{o}} & g_{o,Z_o,Z_o} \\ \end{array} \right) P^T(\phi _o), \end{eqnarray} (10)where PT(ϕo) is the transpose of:   $$P(\phi _o) = \left( \begin{array}{c{@}{\quad}c{@}{\quad}c}{\rm cos} (\phi _o) & -{\rm sin} (\phi _o) & 0 \\ {}{\rm sin} (\phi _o) & {\rm cos} (\phi _o) & 0 \\ 0 & 0 & 1 \\ \end{array} \right).$$ (11) 3 SENSITIVITY ANALYSIS AND FORWARD MODELING This section tests the GLQ-based $$g_{o,Z_o}$$ estimates against the uniform density, vertical cylinder modeling from the other methods described above, illustrates and compares the GLQ-based gravitational estimates of a uniform density, vertical cylinder against the gravitational estimates of the same cylinder from Rim & Li (2016), and illustrates the comprehensive GLQ gravitational effects of the cylindrical prism element. Fig. 2 shows how the tested gravitational effects were modeled using the observation-centred convention for a vertical circular cylinder of constant density σs, length L and radius R with respective depths Z1 and Z2 to its bottom and top surfaces along the XZ-plane passing through the cylinder’s axis. Figs 3–6 give axis-centred XZ-profile comparisons of the $$g_{o,Z_o}$$ anomaly estimates for a vertical cylinder from the GLQ method, CC method (Kwok 1991), EI method (Nabighian 1962), FH method (Singh 1977), LI method (Van Baak 1989), LP method (Telford et al. 1990), P method (Parasnis 1961) plus SA and VL methods (Nettleton 1942). Here, the vertical cylinder was assumed with density σs = 2800 kg m−3, length L = 500 m, top surfaces Z2 at -50, -500, -2500 and -5000 m below the observation plane, and radii R of 50 m (Fig. 3), 500 m (Fig. 4), 2500 m (Fig. 5) and 5000 m (Fig. 6). Figure 2. View largeDownload slide Vertical circular cylinder of constant density σs, length L and radius R with the respective bottom and top surfaces at depths Z1 and Z2 and the XZ-plane passing through the cylinder’s axis. Figure 2. View largeDownload slide Vertical circular cylinder of constant density σs, length L and radius R with the respective bottom and top surfaces at depths Z1 and Z2 and the XZ-plane passing through the cylinder’s axis. Figure 3. View largeDownload slide Gravitational anomaly estimates of a vertical circular cylinder of density σs = 2800 kg m−3, length L = 500 m, radius R = 50 m and top surface at Z2 from the GLQ, CC, EI, FH, LI, LP, P, SA and VL methods. Vertical axis is in logarithmic units of $$g_{o,Z_o}$$ (mGal). Figure 3. View largeDownload slide Gravitational anomaly estimates of a vertical circular cylinder of density σs = 2800 kg m−3, length L = 500 m, radius R = 50 m and top surface at Z2 from the GLQ, CC, EI, FH, LI, LP, P, SA and VL methods. Vertical axis is in logarithmic units of $$g_{o,Z_o}$$ (mGal). Figure 4. View largeDownload slide Gravitational anomaly profile comparisons modeled using the methods listed in Fig. 3 for a vertical circular cylinder of density σs = 2800 kg m−3, length L = 500 m, radius R = 500 m and top surface at Z2. Vertical axis is in logarithmic units of $$g_{o,Z_o}$$ (mGal). Figure 4. View largeDownload slide Gravitational anomaly profile comparisons modeled using the methods listed in Fig. 3 for a vertical circular cylinder of density σs = 2800 kg m−3, length L = 500 m, radius R = 500 m and top surface at Z2. Vertical axis is in logarithmic units of $$g_{o,Z_o}$$ (mGal). Figure 5. View largeDownload slide Gravitational anomaly profile comparisons modeled using the methods listed in Fig. 3 for a vertical circular cylinder of density σs = 2800 kg m−3, length L = 500 m, radius R = 2500 m and top surface at Z2. Vertical axis is in logarithmic units of $$g_{o,Z_o}$$ (mGal). Figure 5. View largeDownload slide Gravitational anomaly profile comparisons modeled using the methods listed in Fig. 3 for a vertical circular cylinder of density σs = 2800 kg m−3, length L = 500 m, radius R = 2500 m and top surface at Z2. Vertical axis is in logarithmic units of $$g_{o,Z_o}$$ (mGal). Figure 6. View largeDownload slide Gravitational anomaly profile comparisons modeled using the methods listed in Fig. 3 for a vertical circular cylinder of density σs = 2800 kg m−3, length L = 500 m, radius R = 5000 m and top surface at Z2. Vertical axis is in logarithmic units of $$g_{o,Z_o}$$ (mGal). Figure 6. View largeDownload slide Gravitational anomaly profile comparisons modeled using the methods listed in Fig. 3 for a vertical circular cylinder of density σs = 2800 kg m−3, length L = 500 m, radius R = 5000 m and top surface at Z2. Vertical axis is in logarithmic units of $$g_{o,Z_o}$$ (mGal). To better visualize the modeled anomaly ranges and details near the cylinder’s axis, the $$g_{o,Z_o}$$ estimates in mGal are displayed in base-10 logarithmic values. For the GLQ modeling, the cylinder’s radial dimensions were filled using 2 (R = 50 m), 11 (R = 500 m), 51 (R = 2500 m) and 101 (R = 5000 m) nodes, whereas the azimuthal dimensions were filled using 360 (R = 50 m), 360 (R = 500 m), 360 (R = 2500 m) and 720 (R = 5000 m) nodes. The length dimensions of all 16 cylinders were filled using 2 nodes. Study of Figs 3(e), 4(e) and (f), 5(e)–(g), 6(e)–(h) reveals that the $$g_{o,Z_o}$$ estimates from the LP, SA and VL methods increasingly differ from the GLQ, CC, EI, FH, LI and P estimates towards the central axis. These differences however, decrease as the R (radii) and Z2 (depths) of the cylinder become smaller and deeper. Fig. 7 compares the rms differences between the $$g_{o,Z_o}$$ estimates from the GLQ modeling and those of the other eight methods described above. To remove the effects of singularities, the gravitational anomaly estimates at X = 0 were deleted from the rms difference calculations. Figs 7(a)–(d) suggest that for constant L (lengths) and R (radii), the rms differences become smaller as the cylinder’s Z2 (depth) increases. Likewise, increasing rms differences result for constant L (lengths) and Z2 (depths) as the R (radii) of the cylinder increases as shown in Figs 7(e)–(h). Fig. 7 also suggests that for all R and Z2 values, the CC, EI and FH estimates are very similar. Figure 7. View largeDownload slide Comparison of $$g_{o,Z_o}$$ rms differences between the GLQ method and the other eight methods described in Section 3. Subfigures (a)–(d) illustrate for a constant radius of R, how changes in the depth to the top surface Z2 affects the rms differences, whereas subfigures (e)–(h) show for a constant depth to the top surface of Z2, how changes in the radius R affects the rms differences. Figure 7. View largeDownload slide Comparison of $$g_{o,Z_o}$$ rms differences between the GLQ method and the other eight methods described in Section 3. Subfigures (a)–(d) illustrate for a constant radius of R, how changes in the depth to the top surface Z2 affects the rms differences, whereas subfigures (e)–(h) show for a constant depth to the top surface of Z2, how changes in the radius R affects the rms differences. Figs 8 and 9 illustrate at 0 m altitude, the cylindrical coordinate gravitational anomaly components for a uniformly dense, vertical cylinder with σs = 1000 kg m−3 and dimensions ρ1 = 0 m to ρ2 = 100 m and Z1 = −350 m to Z2 = −50 m. These effects were GLQ-modeled using 4 × 360 × 2 = 2880 nodes spanning the radial, longitudinal and vertical dimensions of the cylinder’s volume. Figs 8(a)–(d), respectively, show the modeled gravitational Po-potential and vector gravitational $$g_{o,\rho _{o}}$$-, $$g_{o,\phi _{o}}$$- and $$g_{o,Z_o}$$- anomaly components, whereas Fig. 9 gives the related gravitational gradient anomaly components from the three trace and three unique off-diagonal elements of the symmetric matrix in eq. (5). Note that four of the nine elements of this matrix are redundant because $$g_{o,\rho _{o},\phi _{o}} = g_{o,\phi _{o},\rho _{o}}$$, $$g_{o,\rho _{o},Z_o} = g_{o,Z_o,\rho _{o}}$$, and $$g_{o,\phi _{o},Z_o} = g_{o,Z_o,\phi _{o}}$$ and any two components of $$g_{o,\rho _{o},\rho _{o}}, g_{o,\phi _{o},\phi _{o}}$$ and $$g_{o,Z_o,Z_o}$$ may be related to the third one via Laplace’s equation. Within working precision, these results match the estimates from the closed-form gravitational expressions based on complete Lipschitz–Hankel type integral solutions involving elliptic integrals of the first and second kinds and Heumann’s lambda function (Rim & Li 2016). Figure 8. View largeDownload slide Maps (a)–(d) give the respective gravitational potential (Po) and vector radial ($$g_{o,\rho _o}$$), longitudinal ($$g_{o,\phi _o}$$) and vertical ($$g_{o,Z_o}$$) anomaly components of a cylinder with ρ1 = 0 m, ρ2 = 100 m, Z1 = −350 m, and Z2 = −50 m with density σs = 1000 kg m−3 in cylindrical coordinates. The gravitational effects were modeled at 0 m altitude for the cylinder using I = 4, J = 360 and K = 2 nodes. See Table 3 for statistical attributes on these maps. Figure 8. View largeDownload slide Maps (a)–(d) give the respective gravitational potential (Po) and vector radial ($$g_{o,\rho _o}$$), longitudinal ($$g_{o,\phi _o}$$) and vertical ($$g_{o,Z_o}$$) anomaly components of a cylinder with ρ1 = 0 m, ρ2 = 100 m, Z1 = −350 m, and Z2 = −50 m with density σs = 1000 kg m−3 in cylindrical coordinates. The gravitational effects were modeled at 0 m altitude for the cylinder using I = 4, J = 360 and K = 2 nodes. See Table 3 for statistical attributes on these maps. Figure 9. View largeDownload slide Cylindrical coordinate gravitational gradient components for the cylinder in Fig. 8. See Table 3 for statistical attributes on these maps. Figure 9. View largeDownload slide Cylindrical coordinate gravitational gradient components for the cylinder in Fig. 8. See Table 3 for statistical attributes on these maps. Figs 10 and 11 illustrate at 500 m altitude, the cylindrical coordinate gravitational effects of a vertical cylindrical volume element or prism with σs = 2800 kg m−3 and dimensions ρ1 = 3000 m to ρ2 = 4000 m, ϕ1 = 30° to ϕ2 = 60° and Z1 = −10000 m to Z2 = −5000 m. These effects were GLQ modeled using 2 × 2 × 2 = 8 nodes with each nodal pair spanning the radial, longitudinal and vertical dimensions of the prism’s volume. Figs 10(a)–(d), respectively, show the modeled gravitational Po-potential and vector gravitational $$g_{o,\rho _{o}}$$-, $$g_{o,\phi _{o}}$$-, and $$g_{o,Z_o}$$- anomaly components, whereas Fig. 11 gives the gravitational gradient anomaly components from the three trace and three unique off-diagonal elements of the symmetric matrix in eq. (5). The gravitational effects of the cylindrical prism in Figs 10 and 11 are, respectively, mapped in Cartesian coordinates in Figs 12 and 13 via the transformation eqs (9)–(11). In general, the Cartesian mapping is more typical of exploration geophysical applications even though the data in cylindrical and Cartesian projections provide equivalent spatial anomaly detail. Figure 10. View largeDownload slide Maps (a)–(d) give the respective gravitational potential (Po) and vector gravitational radial ($$g_{o,\rho _o}$$), longitudinal ($$g_{o,\phi _o}$$) and vertical ($$g_{o,Z_o}$$) anomaly components of a cylindrical prism with ρ1 = 3000 m, ρ2 = 4000 m, ϕ1 = 30°, ϕ2 = 60°, Z1 = −10  000 m and Z2 = −5000 m with density σs = 2800 kg m−3 in cylindrical coordinates. The gravitational effects were modeled at 500 m altitude for the cylindrical prism using I = 2, J = 2 and K = 2 nodes. See Table 3 for statistical attributes on these maps. Figure 10. View largeDownload slide Maps (a)–(d) give the respective gravitational potential (Po) and vector gravitational radial ($$g_{o,\rho _o}$$), longitudinal ($$g_{o,\phi _o}$$) and vertical ($$g_{o,Z_o}$$) anomaly components of a cylindrical prism with ρ1 = 3000 m, ρ2 = 4000 m, ϕ1 = 30°, ϕ2 = 60°, Z1 = −10  000 m and Z2 = −5000 m with density σs = 2800 kg m−3 in cylindrical coordinates. The gravitational effects were modeled at 500 m altitude for the cylindrical prism using I = 2, J = 2 and K = 2 nodes. See Table 3 for statistical attributes on these maps. Figure 11. View largeDownload slide Cylindrical coordinate gravitational gradient components for the cylindrical prism in Fig. 10. See Table 3 for statistical attributes on these maps. Figure 11. View largeDownload slide Cylindrical coordinate gravitational gradient components for the cylindrical prism in Fig. 10. See Table 3 for statistical attributes on these maps. Figure 12. View largeDownload slide Maps (a)–(d) give the respective gravitational potential (Po) and vector left–right ($$g_{o,X_o}$$), top–bottom ($$g_{o,Y_o}$$) and vertical ($$g_{o,Z_o}$$) anomaly components for the cylindrical prism in Fig. 10 in Cartesian coordinates. Note that anomaly maps of Figs 10(a) and 12(a), as well as Figs 10(d) and 12(d) are the same for both cylindrical and Cartesian coordinate systems. See Table 3 for statistical attributes on these maps. Figure 12. View largeDownload slide Maps (a)–(d) give the respective gravitational potential (Po) and vector left–right ($$g_{o,X_o}$$), top–bottom ($$g_{o,Y_o}$$) and vertical ($$g_{o,Z_o}$$) anomaly components for the cylindrical prism in Fig. 10 in Cartesian coordinates. Note that anomaly maps of Figs 10(a) and 12(a), as well as Figs 10(d) and 12(d) are the same for both cylindrical and Cartesian coordinate systems. See Table 3 for statistical attributes on these maps. Figure 13. View largeDownload slide Gravitational gradient components for the cylindrical prism in Fig. 10 in Cartesian coordinates. Note that anomaly maps of Figs 11(f) and 13(f) are the same for both cylindrical and Cartesian coordinate systems. See Table 3 for statistical attributes on these maps. Figure 13. View largeDownload slide Gravitational gradient components for the cylindrical prism in Fig. 10 in Cartesian coordinates. Note that anomaly maps of Figs 11(f) and 13(f) are the same for both cylindrical and Cartesian coordinate systems. See Table 3 for statistical attributes on these maps. 4 GLQ GRAVITATIONAL MODELING OF A FORMATION WASHOUT Formation washouts are sources of uncertainty in borehole logging. They may form whenever the drill bit penetrates a weak rock formation (e.g. a loosely consolidated sandstone layer), and may extend partially or completely around the borehole wall. In this section, procedures for locating formation washouts and modeling their gravitational effects by GLQ integration are considered. The well-established method based on Compton scattering of gamma rays uses a continuous beam of photons typically emitted from a 0.66 MeV Cs137 source into the formation (Dewan 1983). These photons interact with the electrons of formation material to scatter in all directions with intensities proportional to the electron densities that can be related to the bulk densities of the formation. Modern implementation of this method allows for determination of azimuthal mass density variations (Bargach et al. 2000), which can ultimately be used to locate any formation washouts (density anomalies) around the borehole wall. Another approach for mapping formation washouts around a borehole may be based on borehole gravimetry and related gradiometry (Nekut 1989). Here, a grid of point gravity poles or masses is setup around the circumference of the wellbore. Iterative gravitational modeling is then adopted to update the masses of the gridded point sources until their gravitational effects fit the gravitational observations to within prescribed error limits. The resulting distribution of gridded point masses may finally be used to help map out washouts around the borehole wall. In theory, eqs (3), (4) and (6) can be used for estimating the washout’s gravitational effects at any point within the borehole on and outside the washout’s surface envelop. As examples, Figs 14 and 15 present a washout’s Cartesian coordinates gravitational effects along the respective central XZ- and YZ-planes of a borehole with radius ρBorehole = 0.5 m. The washout extends from ρ1 = 0.5 m to ρ2 = 1 m, ϕ1 = 150° to ϕ2 = 210° and Z1 = −1000.25 m to Z2 = −999.75 m along the respective radial, azimuthal and axial dimensions of the borehole. The assumed density contrast between the washout and its confining formation is σs = −1500 kg m−3. These effects were initially GLQ modeled in cylindrical coordinates via eqs (3), (4) and (6) using 5 × 21 × 20 = 2100 nodes, and subsequently projected into Cartesian coordinates via eqs (9)–(11). These figures indicate that $$g_{o,Y_{o}}$$, $$g_{o,X_{o},Y_{o}} = g_{o,Y_{o},X_{o}}$$ and $$g_{o,Y_{o},Z_{o}} = g_{o,Z_{o},Y_{o}}$$ show little-to-no variation across the XZ-plane in contrast to the other components that vary strongly across this plane, whereas all components vary strongly across the YZ-plane. Clearly, to use the inversion of borehole gravity data to map a formation washout, the related anomaly variations must be observed on a plane within the sensitivity ranges of modern gravimeter and gravity gradiometer sensors (Chen & Cook 2005). Figure 14. View largeDownload slide Maps (a)–(d) give in Cartesian coordinates the respective gravitational Po potential and $$g_{o,X_o}$$, $$g_{o,Y_o}$$ and $$g_{o,Z_o}$$ vector anomaly components in the respective X-, Y-, and Z-directions of a formation washout with ρ1 = 0.5 m, ρ2 = 1 m, ϕ1 = 150°, ϕ2 = 210°, Z1 = −1000.25 m and Z2 = −999.75 m with density σs = −1500 kg m−3. The gravitational effects were modeled along the wellbore’s central XZ- and YZ-planes using I = 5, J = 21 and K = 20 nodes. See Table 3 for statistical attributes on these maps. Figure 14. View largeDownload slide Maps (a)–(d) give in Cartesian coordinates the respective gravitational Po potential and $$g_{o,X_o}$$, $$g_{o,Y_o}$$ and $$g_{o,Z_o}$$ vector anomaly components in the respective X-, Y-, and Z-directions of a formation washout with ρ1 = 0.5 m, ρ2 = 1 m, ϕ1 = 150°, ϕ2 = 210°, Z1 = −1000.25 m and Z2 = −999.75 m with density σs = −1500 kg m−3. The gravitational effects were modeled along the wellbore’s central XZ- and YZ-planes using I = 5, J = 21 and K = 20 nodes. See Table 3 for statistical attributes on these maps. Figure 15. View largeDownload slide Gravitational gradient components for the formation washout in Fig. 14 in Cartesian coordinates. See Table 3 for statistical attributes on these maps. Figure 15. View largeDownload slide Gravitational gradient components for the formation washout in Fig. 14 in Cartesian coordinates. See Table 3 for statistical attributes on these maps. 5 GLQ GRAVITATIONAL MODELING OF A SALT DOME As another application of the GLQ gravitational modeling method, the salt intrusion in Fig. 16 is considered that was mapped within the sedimentary units of the Laurentian Basin off Canada’s eastern coast by gravity and seismic data. The Late Triassic breakup of the North American and African continents (Uchupi & Austin Jr 1979) triggered the formation of a number of basins along the Scotian passive margin including the Laurentian Basin (Adam & Krezsek 2012). This passive margin is known to be associated with evaporate deposition, as well as halokinesis in some areas (Uchupi & Austin Jr 1979; Adam & Krezsek 2012). Figure 16. View largeDownload slide (a) Geographic location of the Laurentian Basin study area involving seismic lines STP-3, STP-4 and STP-21. (b) Processed seismic reflection section for line STP-21. (c) Seismic meta-attribute section for line STP-21, where the dashed white polygon outlines the inferred salt body. Figure 16. View largeDownload slide (a) Geographic location of the Laurentian Basin study area involving seismic lines STP-3, STP-4 and STP-21. (b) Processed seismic reflection section for line STP-21. (c) Seismic meta-attribute section for line STP-21, where the dashed white polygon outlines the inferred salt body. The study region includes the horseshoe-shaped negative gravitational anomaly at 500 m above mean sea level shown in Fig. 17(a) from shipborne survey-constrained free-air gravitational anomaly (FAGA) estimates of the Eigen-6c4 model (Förste et al. 2014). Elements of this anomaly were imaged in 1984–1985 by 29 lines of 2-D multichannel seismic reflection surveying that collected some 3100 km of data. The survey obtained roughly 7000 ms of raw data at a 4-ms sampling rate using source and receiver intervals of 25 m each. The data were reprocessed in 2006 and are now accessible at the Department of Natural Resources Canada. The south-central section of the horseshoe-shaped gravitational minimum includes a circular feature that sits north of the NE-SW trending seismic line STP-21, and between the NW-SE trending seismic lines STP-3 and STP-4 (Fig. 17a). Figure 17. View largeDownload slide (a) Free-air gravitational anomaly (FAGA) estimates from the Eigen-6c4 model, (b) digital elevation model (DEM), (c) GLQ-modeled terrain gravitational effects (TGE) and (d) the complete Bouguer anomalies (CBA) for the study area. The FAGA, TGE and CBA estimates were evaluated at 500 m above sea level. See Table 3 for statistical attributes on these maps. Figure 17. View largeDownload slide (a) Free-air gravitational anomaly (FAGA) estimates from the Eigen-6c4 model, (b) digital elevation model (DEM), (c) GLQ-modeled terrain gravitational effects (TGE) and (d) the complete Bouguer anomalies (CBA) for the study area. The FAGA, TGE and CBA estimates were evaluated at 500 m above sea level. See Table 3 for statistical attributes on these maps. The seismic data from line STP-21 along the southern margin of the circular negative gravitational anomaly reveal a possible salt body for the anomaly’s source (Fig. 16b). To enhance the intrusive body’s edges, seismic pattern recognition analysis (Hashemi et al. 2008) was carried out on the line STP-21 data. A set of seismic data samples within the salt body was selected to compare against another data set sampled from the remaining sedimentary formations within the Basin. The comparison facilitated calculating the relevant seismic attributes in a non-redundant feature space (Hashemi & Javaherian 2009) that, in turn, were subjected to principal component analysis for a supervised multilayer perception (MLP) neural network classifier function (Hashemi et al. 2008). The MLP mapping function was applied to all traces of the STP-21 line to infer the spatial extent of the salt body shown in Fig. 16(c). This analysis estimated the body’s average width, average two-way-time (TWT) thickness and the average TWT between its top surface and the seabed at about 7700 m (R = 3850 m), 5283 ms and 173 ms, respectively. The thickness of the seawater column above the salt body amounts to roughly TWT = 480 ms, so that for assumed seawater, sediment and salt intrusion compressional velocities of VP,Seawater = 1484 m s−1, VP,Sediment = 1600 m s−1 and VP,Salt = 4450 m s−1, the depth below sea level of the seabed, and the body’s top surface below the seabed and the salt body’s thickness amount to roughly 356, 138 and 11755 m, respectively. Assuming an average salt density of σSalt = 2170 kg m−3, the cylindrical body’s vertical-component gravitational anomaly ($$g_{o,Z_o}$$) at 500 m above sea level was GLQ modeled as shown by the red curve in Fig. 18(a). The GLQ gravitational modeling used 8 × 360 × 2 = 5760 nodes in the respective radial, azimuthal and depth dimensions. Using GLQ-spherical prism gravitational modeling (Asgharzadeh et al. 2007), the terrain gravitational effects (TGE) of the digital elevation model (DEM, Smith & Sandwell 1997) in Fig. 17(b) were calculated at 500 m above sea level (Fig. 17c) assuming average sea water, sediment, and basement rock densities of σSeawater = 1030 kg m−3, σSediment = 2250 kg m−3 and σBasement = 3030 kg m−3, respectively. The TGE, in turn, were subtracted from the FAGA to obtain the region’s complete Bouguer gravitational anomalies (CBA) in Fig. 17(d). In Fig. 18(a), the magenta and green curves illustrate the CBA profile and its linear regional trend along line STP-21. The linear trend was subtracted from the CBA as possible effects due to lack of isostatic compensation and other regional modeling errors within the basement topography and intrusive body. The removal of the linear regional yields a close fit between the residual CBA (blue curve) and the modeled $$g_{o,Z_o}$$ effects of the salt intrusion (Fig. 18b). Accordingly, the choices for sea water, sediment, salt and basement density, velocity and depth variations yield gravity signals that are consistent with the region’s DEM, seismic reflection and FAGA observations. Figure 18. View largeDownload slide (a) CBA (magenta curve) with its linear regional trend (green curve), and the superimposed GLQ modeled $$g_{o,Z_o}$$ gravitational anomaly of the seismically constrained intrusive salt body (red curve) along seismic line STP-21. (b) Residual CBA (blue curve) with the superimposed GLQ modeled $$g_{o,Z_o}$$ gravitational anomaly of the intrusive salt body (red curve) along the STP-21 seismic line. Figure 18. View largeDownload slide (a) CBA (magenta curve) with its linear regional trend (green curve), and the superimposed GLQ modeled $$g_{o,Z_o}$$ gravitational anomaly of the seismically constrained intrusive salt body (red curve) along seismic line STP-21. (b) Residual CBA (blue curve) with the superimposed GLQ modeled $$g_{o,Z_o}$$ gravitational anomaly of the intrusive salt body (red curve) along the STP-21 seismic line. 6 CONCLUSIONS AND RECOMMENDATIONS The complete gravitational effects (i.e. potential, anomaly and gradient components) of the vertical cylinder with radial-, azimuthal- and axial-density variations were developed in the context of the least-squares numerical gravitational effects of the cylindrical unit volume element or prism estimated by GLQ integration. The GLQ modeling was developed to accommodate gravity studies of rock formations with heterogeneous densities due to physical or chemical alterations. Comparisons against the results of modeling simple uniform density vertical cylinders by a number of conventional methods illustrated the veracity of the GLQ approach. As an application in borehole gravity, the gravitational effects of a formation washout were estimated to assist in determining azimuthal variations of formation bulk densities around the borehole wall. A further application considered the gravitational effects of a seismically and gravimetrically imaged salt body within the Laurentian Basin for constraining the velocity, density and geometric properties of the Basin’s sedimentary formations. Table 3 lists the statistical attributes for Figs 8–15 and 17 to facilitate the efforts of readers to replicate these results. Table 3. Statistical attributes for the maps in Figs 8–15 and 17 include the map altitude (Zo), as well as amplitude minimum (Min), maximum (Max), average (Mean) and standard deviation (Std). The amplitude units (AU) and the power (n) of the amplitude multiplier (× 10n) are also listed. Figure  Map  Zo(m)  Min  Max  Mean  Std  n  AU  8a  Po  0  43.9899  343.8151  92.7064  48.1505  0  mGal × m  8b  $$g_{o,\rho _o}$$  0  −0.7800  −0.0304  −0.1467  0.1478  0  mGal  8c  $$g_{o,\phi _o}$$  0  0  0  0  0  0  mGal  8d  $$g_{o,Z_o}$$  0  −1.9163  −0.0043  −0.0806  0.1977  0  mGal  9a  $$g_{o,\rho _o,\rho _o}$$  0  −10.1395  2.6059  0.3666  0.8223  −3  mGal m−1  9b  $$g_{o,\rho _o,\phi _o} = g_{o,\phi _o,\rho _o}$$  0  0  0  0  0  −3  mGal m−1  9c  $$g_{o,\rho _o,Z_o} = g_{o,Z_o,\rho _o}$$  0  0.0088  9.8842  0.5130  1.4259  −3  mGal m−1  9d  $$g_{o,\phi _o,\phi _o}$$  0  −10.1970  −0.0215  −0.4492  1.1290  −3  mGal m−1  9e  $$g_{o,\phi _o,Z_o} = g_{o,Z_o,\phi _o}$$  0  0  0  0  0  −3  mGal m−1  9f  $$g_{o,Z_o,Z_o}$$  0  −0.3725  20.3365  0.0826  1.5095  −3  mGal m−1  10a  Po  500  6563.1  22051  11623  3453.9  0  mGal × m  10b  $$g_{o,\rho _o}$$  500  −1.1071  0.9918  −0.5909  0.2327  0  mGal  10c  $$g_{o,\phi _o}$$  500  −1.0101  1.0101  0  0.2282  0  mGal  10d  $$g_{o,Z_o}$$  500  −2.9206  −0.0765  −0.5441  0.5504  0  mGal  11a  $$g_{o,\rho _o,\rho _o}$$  500  −0.3990  0.0752  0.0260  0.0688  −3  mGal m−1  11b  $$g_{o,\rho _o,\phi _o} = g_{o,\phi _o,\rho _o}$$  500  −0.1059  0.1059  0  0.0335  −3  mGal m−1  11c  $$g_{o,\rho _o,Z_o} = g_{o,Z_o,\rho _o}$$  500  −0.3353  0.3363  0.0795  0.0764  −3  mGal m−1  11d  $$g_{o,\phi _o,\phi _o}$$  500  −0.3941  −0.0092  −0.0612  0.0613  −3  mGal m−1  11e  $$g_{o,\phi _o,Z_o} = g_{o,Z_o,\phi _o}$$  500  −0.3360  0.3360  0  0.0600  −3  mGal m−1  11f  $$g_{o,Z_o,Z_o}$$  500  −0.0125  0.7930  0.0352  0.1208  −3  mGal m−1  12a  Po  500  6563.1  22051  11623  3453.9  0  mGal × m  12b  $$g_{o,X_o}$$  500  −1.1065  1.1066  0  0.4772  0  mGal  12c  $$g_{o,Y_o}$$  500  −1.1065  1.1066  0  0.4772  0  mGal  12d  $$g_{o,Z_o}$$  500  −2.9206  −0.0765  −0.5441  0.5504  0  mGal  13a  $$g_{o,X_o,X_o}$$  500  −0.3965  0.0752  −0.0176  0.0736  −3  mGal m−1  13b  $$g_{o,X_o,Y_o} = g_{o,Y_o,X_o}$$  500  −0.1084  0.1078  0  0.0430  −3  mGal m−1  13c  $$g_{o,X_o,Z_o} = g_{o,Z_o,X_o}$$  500  −0.3356  0.3355  0  0.0888  −3  mGal m−1  13d  $$g_{o,Y_o,Y_o}$$  500  −0.3965  0.0752  −0.0176  0.0736  −3  mGal m−1  13e  $$g_{o,Y_o,Z_o} = g_{o,Z_o,Y_o}$$  500  −0.3356  0.3355  0  0.0888  −3  mGal m−1  13f  $$g_{o,Z_o,Z_o}$$  500  −0.0125  0.7930  0.0352  0.1208  −3  mGal m−1  14a  Po  −  −0.0050  −0.0009  −0.0018  0.0006  0  mGal × m  14b  $$g_{o,X_o}$$  −  0.0001  0.0107  0.0013  0.0013  0  mGal  14c  $$g_{o,Y_o}$$  −  −0.0012  0.0012  0  0.0003  0  mGal  14d  $$g_{o,Z_o}$$  −  −0.0051  0.0051  0  0.0011  0  mGal  15a  $$g_{o,X_o,X_o}$$  −  −34.8223  1.9380  −1.8250  3.8013  −3  mGal m−1  15b  $$g_{o,X_o,Y_o} = g_{o,Y_o,X_o}$$  −  −3.4642  3.4642  0  0.7612  −3  mGal m−1  15c  $$g_{o,X_o,Z_o} = g_{o,Z_o,X_o}$$  −  −25.9132  25.9132  0  3.0304  −3  mGal m−1  15d  $$g_{o,Y_o,Y_o}$$  −  0.2082  10.7654  1.5485  1.6250  −3  mGal m−1  15e  $$g_{o,Y_o,Z_o} = g_{o,Z_o,Y_o}$$  −  −1.1407  1.1407  0  0.3227  −3  mGal m−1  15f  $$g_{o,Z_o,Z_o}$$  −  −7.5209  24.0788  0.2765  2.5012  −3  mGal m−1  17a  FAGA  500  −36.2590  32.6437  −3.4772  13.2437  0  mGal  17b  DEM  −  −3250  −3  −681  787  0  m  17c  TGE  500  −76.7469  −0.6997  −17.4067  19.6542  0  mGal  17d  CBA  500  −27.1880  61.7900  13.7376  19.6185  0  mGal  Figure  Map  Zo(m)  Min  Max  Mean  Std  n  AU  8a  Po  0  43.9899  343.8151  92.7064  48.1505  0  mGal × m  8b  $$g_{o,\rho _o}$$  0  −0.7800  −0.0304  −0.1467  0.1478  0  mGal  8c  $$g_{o,\phi _o}$$  0  0  0  0  0  0  mGal  8d  $$g_{o,Z_o}$$  0  −1.9163  −0.0043  −0.0806  0.1977  0  mGal  9a  $$g_{o,\rho _o,\rho _o}$$  0  −10.1395  2.6059  0.3666  0.8223  −3  mGal m−1  9b  $$g_{o,\rho _o,\phi _o} = g_{o,\phi _o,\rho _o}$$  0  0  0  0  0  −3  mGal m−1  9c  $$g_{o,\rho _o,Z_o} = g_{o,Z_o,\rho _o}$$  0  0.0088  9.8842  0.5130  1.4259  −3  mGal m−1  9d  $$g_{o,\phi _o,\phi _o}$$  0  −10.1970  −0.0215  −0.4492  1.1290  −3  mGal m−1  9e  $$g_{o,\phi _o,Z_o} = g_{o,Z_o,\phi _o}$$  0  0  0  0  0  −3  mGal m−1  9f  $$g_{o,Z_o,Z_o}$$  0  −0.3725  20.3365  0.0826  1.5095  −3  mGal m−1  10a  Po  500  6563.1  22051  11623  3453.9  0  mGal × m  10b  $$g_{o,\rho _o}$$  500  −1.1071  0.9918  −0.5909  0.2327  0  mGal  10c  $$g_{o,\phi _o}$$  500  −1.0101  1.0101  0  0.2282  0  mGal  10d  $$g_{o,Z_o}$$  500  −2.9206  −0.0765  −0.5441  0.5504  0  mGal  11a  $$g_{o,\rho _o,\rho _o}$$  500  −0.3990  0.0752  0.0260  0.0688  −3  mGal m−1  11b  $$g_{o,\rho _o,\phi _o} = g_{o,\phi _o,\rho _o}$$  500  −0.1059  0.1059  0  0.0335  −3  mGal m−1  11c  $$g_{o,\rho _o,Z_o} = g_{o,Z_o,\rho _o}$$  500  −0.3353  0.3363  0.0795  0.0764  −3  mGal m−1  11d  $$g_{o,\phi _o,\phi _o}$$  500  −0.3941  −0.0092  −0.0612  0.0613  −3  mGal m−1  11e  $$g_{o,\phi _o,Z_o} = g_{o,Z_o,\phi _o}$$  500  −0.3360  0.3360  0  0.0600  −3  mGal m−1  11f  $$g_{o,Z_o,Z_o}$$  500  −0.0125  0.7930  0.0352  0.1208  −3  mGal m−1  12a  Po  500  6563.1  22051  11623  3453.9  0  mGal × m  12b  $$g_{o,X_o}$$  500  −1.1065  1.1066  0  0.4772  0  mGal  12c  $$g_{o,Y_o}$$  500  −1.1065  1.1066  0  0.4772  0  mGal  12d  $$g_{o,Z_o}$$  500  −2.9206  −0.0765  −0.5441  0.5504  0  mGal  13a  $$g_{o,X_o,X_o}$$  500  −0.3965  0.0752  −0.0176  0.0736  −3  mGal m−1  13b  $$g_{o,X_o,Y_o} = g_{o,Y_o,X_o}$$  500  −0.1084  0.1078  0  0.0430  −3  mGal m−1  13c  $$g_{o,X_o,Z_o} = g_{o,Z_o,X_o}$$  500  −0.3356  0.3355  0  0.0888  −3  mGal m−1  13d  $$g_{o,Y_o,Y_o}$$  500  −0.3965  0.0752  −0.0176  0.0736  −3  mGal m−1  13e  $$g_{o,Y_o,Z_o} = g_{o,Z_o,Y_o}$$  500  −0.3356  0.3355  0  0.0888  −3  mGal m−1  13f  $$g_{o,Z_o,Z_o}$$  500  −0.0125  0.7930  0.0352  0.1208  −3  mGal m−1  14a  Po  −  −0.0050  −0.0009  −0.0018  0.0006  0  mGal × m  14b  $$g_{o,X_o}$$  −  0.0001  0.0107  0.0013  0.0013  0  mGal  14c  $$g_{o,Y_o}$$  −  −0.0012  0.0012  0  0.0003  0  mGal  14d  $$g_{o,Z_o}$$  −  −0.0051  0.0051  0  0.0011  0  mGal  15a  $$g_{o,X_o,X_o}$$  −  −34.8223  1.9380  −1.8250  3.8013  −3  mGal m−1  15b  $$g_{o,X_o,Y_o} = g_{o,Y_o,X_o}$$  −  −3.4642  3.4642  0  0.7612  −3  mGal m−1  15c  $$g_{o,X_o,Z_o} = g_{o,Z_o,X_o}$$  −  −25.9132  25.9132  0  3.0304  −3  mGal m−1  15d  $$g_{o,Y_o,Y_o}$$  −  0.2082  10.7654  1.5485  1.6250  −3  mGal m−1  15e  $$g_{o,Y_o,Z_o} = g_{o,Z_o,Y_o}$$  −  −1.1407  1.1407  0  0.3227  −3  mGal m−1  15f  $$g_{o,Z_o,Z_o}$$  −  −7.5209  24.0788  0.2765  2.5012  −3  mGal m−1  17a  FAGA  500  −36.2590  32.6437  −3.4772  13.2437  0  mGal  17b  DEM  −  −3250  −3  −681  787  0  m  17c  TGE  500  −76.7469  −0.6997  −17.4067  19.6542  0  mGal  17d  CBA  500  −27.1880  61.7900  13.7376  19.6185  0  mGal  View Large Further exploitation of the GLQ modeling results is possible by adapting them for modeling the gravitational effects of arbitrarily tilted cylinders (Rim & Li 2016). By Poisson’s theorem (Hinze et al. 2013), these results may also be extended to model the magnetic and thermal effects of arbitrarily oriented tubular geological bodies. ACKNOWLEDGEMENTS We thank the National Elite Foundation of Iran and the Institute of Geophysics at the University of Tehran for supporting this work. We also thank the anonymous reviewers for their constructive comments that helped to improve the manuscript. REFERENCES Adam J., Krezsek C., 2012. Basin-scale salt tectonic processes of the Laurentian Basin, Eastern Canada: insights from integrated regional 2D seismic interpretation an 4D physical experiments, Geol. Soc. Lond. Spec. Publ. , 363, 331– 360. Google Scholar CrossRef Search ADS   Asgharzadeh M.F., von Frese R.R.B., Kim H., Leftwich T.E., Kim J.W., 2007. Spherical prism gravity effects by Gauss-Legendre quadrature integration, Geophys. J. Int. , 169( 1), 1– 11. Google Scholar CrossRef Search ADS   Bargach S., Bornemann T., Codazzi D., Ford G., Grether B., Rohler H., 2000. Real-time LWD: logging for drilling, Oilfield Rev. , 12( 3), 58– 78. Brown A.R., Lautzenhiser T.V., 1982. The effect of dipping beds on a borehole gravimeter survey, Geophysics , 47( 1), 25– 30. Google Scholar CrossRef Search ADS   Chen Y.T., Cook A., 2005. Gravitational Experiments in the Laboratory , Cambridge University Press, Cambridge. Dewan J.T., 1983. Essentials of Modern Open-hole Log Interpretation , PennWell Publishing, Tusla, Oklahoma. Förste C. et al.  , 2014. EIGEN-6C4 The latest combined global gravity field model including GOCE data up to degree and order 2190 of GFZ Potsdam and GRGS Toulouse, GFZ Data Services, doi:10.5880/icgem.2015.1. Hashemi H., Javaherian A., 2009. Seismic attribute redundancy reduction using statistical feature extraction technique, in 1st International Petroleum Conference and Exhibition, EAGE, Shiraz, Iran , doi: 10.3997/2214-4609.20145883. Hashemi H., Tax D.M.J., Duin R.P.W., Javaherian A., de Groot P., 2008. Gas chimney detection based on improving the performance of combined multilayer perceptron and support vector classifier, Nonlinear Process. Geophys. , 15, 863– 871. Google Scholar CrossRef Search ADS   Hinze W.J., von Frese R. R.B., Saad A.H., 2013. Gravity and Magnetic Exploration: Principles, Practices, and Applications , Cambridge University Press, Cambridge. Google Scholar CrossRef Search ADS   Jenkins A.J.O., Messfin D., Moon W., 1983. Gravity modelling of salt domes and pinnacle reefs, J. Can. Soc. Explor. Geophys. , 19( 1), 51– 56. Ku C.C., 1977. A direct computation of gravity and magnetic anomalies caused by 2- and 3-dimensional bodies of arbitrary shape and arbitrary magnetic polarization by equivalent point method and a simplified cubic spline, Geophysics , 42, 610– 622. Google Scholar CrossRef Search ADS   Kwok Y.-K., 1991. Singularities in gravity computation for vertical cylinders and prisms, Geophys. J. Int. , 104, 1– 10. Google Scholar CrossRef Search ADS   Li Z., Hao T., Xu Y., Xu U., 2011. An efficient and adaptive approach for modeling gravity effects in spherical coordinates, J. appl. Geophys. , 73, 221– 231. Google Scholar CrossRef Search ADS   Nabighian M.N., 1962. The gravitational attraction of a right vertical circular cylinder at points external to it, Geofis. pura appl. , 53( 1), 45– 51. Google Scholar CrossRef Search ADS   Nagy D., 1966. The gravitational attraction of a right rectangular prism, Geophysics , 31( 2), 362– 371. Google Scholar CrossRef Search ADS   Nekut A.G., 1989. Borehole gravity gradiometry, Geophysics , 54( 2), 225– 234. Google Scholar CrossRef Search ADS   Nettleton L.L., 1942. Gravity and magnetic calculations, Geophysics , 7( 3), 293– 310. Google Scholar CrossRef Search ADS   Parasnis D.S., 1961. Exact expressions for the gravitational attraction of a circular lamina at all points of space and of a right circular vertical cylinder at points external to it, Geophys. Prospect. , 9( 3), 382– 398. Google Scholar CrossRef Search ADS   Reilly W.I., 1969. Gravitational and magnetic effects of a right circular cylinder, N. Z. J. Geol. Geophys. , 12( 2–3), 497– 506. Google Scholar CrossRef Search ADS   Rim H., Li Y., 2012. Single-hole imaging using borehole gravity gradiometry, Geophysics , 77( 5), G67– G76. Google Scholar CrossRef Search ADS   Rim H., Li Y., 2016. Gravity gradient tensor due to a cylinder, Geophysics , 81( 4), 59– 66. Google Scholar CrossRef Search ADS   Roussel C., Verdun J., Cali J., Masson F., 2015. Complete gravity field of an ellipsoidal prism by Gauss-Legendre quadrature, Geophys. J. Int. , 203( 3), 2220– 2236. Google Scholar CrossRef Search ADS   Singh S.K., 1977. Gravitational attraction of a vertical right circular cylinder, Geophys. J. R. astr. Soc , 50, 243– 246. Google Scholar CrossRef Search ADS   Smith N.J., 1950. The case for gravity data from boreholes, Geophysics , 15( 4), 605– 636. Google Scholar CrossRef Search ADS   Smith W. H.F., Sandwell D.T., 1997. Global seafloor topography from satellite altimetry and ship depth soundings, Science , 277, 1957– 1962. Stroud A.H., Secrest D., 1966. Gaussian Quadrature Formulas , Prentice-Hall, Upper Saddle River, NJ. Telford W.M., Geldart L.P., Sheriff R.E., 1990. Applied Geophysics , Cambridge University Press, Cambridge. Google Scholar CrossRef Search ADS   Uchupi E., Austin Jr J.A., 1979. The geologic history of the passive margin off New England and the Canadian maritime provinces, Tectonophysics , 59, 53– 69. Google Scholar CrossRef Search ADS   Van Baak D.A., 1989. Efficient computation of the gravitational field of a generalized cylindrical prism, Geophysics , 54( 3), 402– 405. Google Scholar CrossRef Search ADS   von Frese R. R.B., Hinze W.J., Braile L.W., Luca A.J., 1981. Spherical earth gravity and magnetic anomaly modeling by Gauss-Legendre quadrature integration, J. Geophys. , 49, 234– 242. © The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Jan 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations