University of Birmingham Development of an infinite element boundary to model gravity for subsurface civil engineering applications

A comprehensive of to reproduce the boundary effect using layer of The developed The model validated of the and its results show excellent performance. The proposed method can significantly improve the postprocessing and interpretation of data analysis relevant to micro-gravity sensors. The new method is applied to subsurface civil engineering although its applicability is manifold.

FEM have been successfully used by different researchers for forward modelling of gravity data (e.g., Aitken et al, 16 Schaa et al, 17 and Zhang et al 18 ). There are a number of advantages of FEM over analytical methods, such as the ability to model complex geometries and complicated density distributions. 19 Additionally, numerical methods offer the advantage of obtaining gravity values over the whole numerical domain without requiring additional computational time. A summary of the approaches for modelling gravity in numerical models is discussed below, highlighting gaps in the current knowledge.
Numerical formulations of gravity include the solution of the Poisson equation (Equation (1)) for regions where mass density values exist, and Laplace equation (Equation (2)) for zones where the density contrast vanishes. Furthermore, gravity is an unbounded problem, and when modelled numerically, it requires a boundary condition, which has to reproduce a case where gravity becomes zero at infinity (ie, Dirichlet boundary condition). 20 where is the gravitational potential, Ω ∞ is the entire domain, G is the gravitational constant, ( , , ) = 0 outside of the mass subdomain (Ω M ), i.e., in the zone where there is no density contrast, and , , are components of a local coordinate system. As shown in Equation (1), the solution of the Poisson's equation gives the scalar gravitational potential ( ). The gravity field, g(x, y, z), can be obtained by taking the vertical derivative of (Equation (3)).
Since it is impossible to numerically model infinity due to computational limitations, efforts have been made by researchers (e.g., 18,19 ) to approximate the "infinity" condition of the unbounded domain of the gravitational potential field with various degrees of accuracy. A simple technique is the truncation of the solution domain at a large but finite distance in a way that the effect of the boundary on the results is significantly reduced (i.e., as much as possible depending on the problem, computational limitations, and required degree of accuracy). Zhang et al 18 and Jahandari and Farquharson 11 among others have applied this approach by adding many empty layers (with zero density contrast) around the mass source and assigning = 0 to the outer boundary of the empty layers. This method is called the zero potential boundary condition for the rest of this paper.
Cai and Wang 19 assumed that, at a specific distance from the mass subdomain, the value of the gravitational potential can be computed and applied to the domain boundary. This can result in a reasonable reduction in the size of the numerical model. At a large distance, r, from the mass, they approximated gravity at the far-field domain boundary using Equation (4).
where S is the domain boundary (surface). May and Knepley 21 presented a discussion about the above two numerical boundary conditions (zero potential and the Cai and Wang 19 methods). By making comparisons with other methods, for example, the summation based and fast multipole methods, they concluded that, although the boundary condition of Cai and Wang 19 improved the performance of the numerical analyses (compared to the zero potential boundary case), the convergence rates observed from the finite element analyses were lower than that of the other methods. In addition to the outcomes of May and Knepley, 21 both inaccuracy and computational expense are frequently (separately or together) associated with the truncation of infinite domains. 22 An effective solution for unbounded problems is to add a layer of infinite elements as the boundary condition for the numerical model. Infinite elements are used to reproduce the effect of the far-field boundary on an unbounded finite element domain. As explained by Kumar,23 there are two approaches to implement infinite elements in a finite element analysis: displacement descent and coordinate ascent formulations. In the former method, an infinite element is formed by stretching a finite element, which requires integration over a semi-infinite range in order to compute element properties. The latter method offers more convenience since an infinite element is compressed to be a regular shape finite element where the Gauss-Legendre numerical integration can be applied conveniently to calculate the properties of the infinite element. 23 Butler and Sinha 24 used the Comsol Multiphysics Package for gravity forward modelling and compared different boundary conditions including infinite elements. They simulated a simple model consisting of a cylinder with a radius of 1000 m and a density contrast of 100 kg m −3 with the surrounding material. The cylinder was inside a cube of 10-km side length, of which 2 km were boundary infinite elements from each side. Butler and Sinha 24 reported good predictions of the gravity data using infinite elements. It should be mentioned that the focus of Butler and Sinha 24 was on showing the capability of the Comsol Multiphysics Package for gravity forward modelling using different boundary conditions. The paper mainly showed the application of Comsol rather than formulation details of the boundary conditions used in their work (e.g., infinite elements).
Gharti and Tromp 25 and Gharti et al 26 used an infinite element boundary in a spectral finite element analysis to solve the Poisson's equation and to model gravity anomalies. They called the method "spectral-infinite-element method (SIEM)." A spectral FEM, as proposed by Patera,27 is an integration of the spectral methods (in which a high-order orthogonal expansion is involved in the solution of a differential equation) and finite element analysis. Combining spectral techniques with FEM adds accuracy to the numerical solution of the differential equations. 27,28 Gharti and Tromp 25 used their developed SIEM to solve the unbounded Poisson's equation for a sphere with a radius of 1000 m and a mass density of 1.92 kg m −3 . They implemented parallel iterative Krylov solvers using PETSc (portable and extensible toolkit for scientific computation). Gharti et al 26 showed the application of SIEM, developed by Gharti and Tromp, 25 for computing gravity anomalies using MeshAssist 29 and Trelis/CUBIT 30 to prepare their models. While this approach resolved the features accurately, it should be mentioned that the details of the infinite element shape functions used by Gharti and Tromp 25 and Gharti et al 26 were not presented in their work. Additionally, the implementation of their approach needs a high computational power and access to commercial packages.
Despite several useful features, such as efficient computational time with a high degree of accuracy, there is very limited research in the literature of using infinite elements as the boundary condition for gravity forward modelling, and they all lack detail and an in-depth study resulting in a significant reduction of their widespread use and adaptation. Furthermore, the existing developments in modelling gravity data using numerical methods mainly focus on applications with large density contrasts, such as soil stratification, deep ore mines, big cavities, and underground water/oil reservoirs (e.g., Gharti et al, 26 Farquharson, 31 and Dai et al 32 ). This paper presents a detailed formulation of boundary infinite elements in a 3D classical finite element analysis for forward modelling of gravity data. Mapping functions for elements with three directional infinities are derived to reproduce infinity at the boundaries of the domain. This formulation can be recreated using any programming language, e.g., MATLAB, without the need to access commercial software packages. The outcomes of the work are then compared with the results obtained from accurate closed-form solutions derived from the exact solution for prisms. Additionally, this work is the first attempt to introduce the application of gravity modelling, using numerical methods, to the field of civil engineering, in which small objects such as utility tunnels, large sewer pipes, shallow mineshafts, and culverts are considered. To assess the suitability of the developed approach, it is tested on a subsurface civil engineering-related application of a shallow multi-utility tunnel with input from field gravity observations as a practical example.

FINITE ELEMENT FORMULATION
In finite element analysis, the problem domain is first discretised into small elements, and Equation (1) is numerically approximated over each element. The domain in this work is discretised into brick-shaped elements having eight nodes at their corners. These nodes govern interconnection between elements. The problem is solved for each finite element in a local coordinate system (−1 ≤ ≤ 1, −1 ≤ ≤ 1, −1 ≤ ≤ 1), and then the global solution is obtained. Equation (1) can be expressed numerically as Equation (?) for each finite element.
where K e is the gravitational stiffness matrix of an element, Φ e = [ 1 , 2 , 3 , ..., n p ] T is the gravitational potential vector, n p is the total number of nodes in an element, T stands for the transpose of the matrix/vector, and F e is a vector containing the density distribution of the domain. The unknown in Equation (?) is Φ e , which can be solved after evaluating K e and F e . It is worth noting that FEM solves the problem at discrete nodes. To obtain solutions over an element, shape functions have to be introduced in order to interpolate the solution between the nodal values. For each element, n p shape functions are required in order to consider the effect of all nodes. The vector of shape functions can be expressed as where h i is the shape function of node i. Shape functions are subjected to the following conditions: h i = 1 at node i and 0 at all other nodes, and ∑n p i=1 h i = 1. Global coordinates of each element, (x, y, z), can be expressed as where x i , y i , z i are global coordinates of node i. The conversion between the global and local coordinate systems can be done using the Jacobian matrix (J, Equation (7)), which is formed from the derivatives of x, y, and z (Equation (6)). It should be noted that the shape functions used to interpolate the solution between the nodes are the same as those used to form the Jacobian (transformation) matrix, which makes the analysis isoparametric for the finite element domain.
K e and F e can be calculated from the following equations: where detJ is the determinant of the Jacobian matrix, and B is calculated using Equations (10) and (11).
If the density distribution is constant over an element, which is the case in this paper, the density can be taken out of the integration term in Equation (9). The integration of partial differential equations in this work is performed using an eight-point Gauss-Legendre numerical integration. After completing the above calculations, the assembly of the elemental matrices and vectors leads to the formulation of the global linear equation system: The next section shows how the infinite elements are formulated.

FORMULATION OF INFINITE ELEMENTS
The boundary of the problem domain is modelled using a single layer of infinite elements. The coordinate ascent formulation of infinite elements is used for modelling the infinite element boundary. In this section, the concept of infinite elements is explained in detail for a one dimensional problem, after that additional detail is added for the 3D case considered in this research. Figure 1A shows a one-dimensional domain extending from a known point to infinity. If the known point is the edge of a finite domain modelled with conventional finite elements, the infinite part of Figure 1A can be mapped to the finite domain using infinite elements. Figure 1B shows the physical and the mapped infinite elements. A pole (point 0 with Figure 1B) is arbitrarily specified for the infinite element, which is then used to determine the length of the infinite element. It should be mentioned that node 2 in Figure 1B represents the end of the infinite element, and its distance (x 2 ) is computed from Equation (13). 33 This shows that the infinity point is not modelled explicitly.
The position of the pole has an impact on the results and should be chosen carefully. An incorrect pole position will lead to a disturbance or distortion to the mapping. 34 More information about the effect of the pole can be found in Marques and Owen 33 and Zienkiewicz et al. 22 The geometry and the physical characteristics of the problem have to be considered when choosing the position of the pole. 33 Moreover, fitting the results to the expected decay pattern could also be helpful to position the pole. 22 Mapping functions for infinite elements contain a singular term, which tends to infinity when tends to +1 ( Figure 1B). Therefore, mainly nodes 1 and 2 are considered in the calculations. It should be noted that using Gauss-Legendre numerical integration eliminates the problem of singularity since integration points are taken within the element rather than at its corners. The mapping functions of nodes 1 and 2 in Figure 1B are expressed by Equation (14). 33 The summation of the mapping functions (M) is done over the finite nodes (1 and 2 in Figure 1B). Note that ∑ 2 i=1 M i = 1 is satisfied for any not equal to +1.
Mapping functions of the infinite elements can also be formed using nodes 0 and 2 ( Figure 1B). 22

Formulation of 3D infinite elements
The one dimensional mapping function (Equation (14)) can be extended to 2D and 3D. For the 3D case considered in this paper, there are three types of infinite elements: one-direction, two-direction (edge), and three-direction (corner) infinite elements. It should be mentioned that for the interpolation of the field variables (i.e., Φ), interpolation functions of the finite elements are used. Furthermore, the nodes of the infinite elements that are exposed to infinity are excluded from the analysis of computing gravitational potential since in theory, the value of the gravitational potential becomes zero at infinity. 33,35 The calculation procedures of infinite elements are the same as that of finite elements except in the computation of the Jacobian matrix.
Mapping functions in the infinite elements are used to compute the Jacobian (transformation) matrix. Note that the nodes exposed to infinity, which are excluded in the calculation of field variables, have to be added to the calculation of the Jacobian matrix in order to achieve analysis convergence. 33 In this paper, mapping functions for one-direction elements are obtained from Marques and Owen 33 ; for two-direction elements, they are obtained from Xia and Zhang 36 ; and for three-direction elements, they are derived as part of the contributions of this work and presented in Table 1. Figure 2 displays a typical 3D model showing the position of one-, two-, and three-direction infinite elements. Global coordinates of the infinite elements (x inf , y inf , z inf ) can be computed from Equation (15) (using mapping functions), and the Jacobian matrix (J inf ) to transform results from local to global coordinates is computed from Equation (16).
Note that in the computations related to the infinite elements, the interpolation functions of field variables are of lower degree than the mapping functions of the element coordinates. This makes the analysis superparametric (versus isoparametric for finite elements).
The next section describes how this information can be combined into a 3D numerical model for gravity forward modelling.

MODEL DESCRIPTION
An illustrative numerical model of this work is shown in Figure 3A. The zone of density contrast in the model was a cube of L M = 80 m (side length). One layer of infinite elements was added to the boundary of the domain-end of the anomalous zone-to model the boundary condition. Figure 3B shows a 2D cross section of the 3D model analysed in this work. The density contrast in the anomalous zone was 1800 kg/m 3 and zero in the infinite element boundary layer. The performance of the adopted infinite elements was tested for different scenarios, such as investigating the effects of the mesh fineness, and the position of the pole. In order to validate the results of the developed infinite element boundary, synthetic gravity data were generated using the closed-form solution of the right rectangular prism (Equation (17)) presented in Li and Chouteau. 37 These closed-form solutions were derived from the Newtonian gravity equation. In addition, the whole density cube was taken as one right rectangular prism for the synthetic data. It is worth noting that Equation (17) can be used to compute gravity values in the inside and outside of the prism.
where ijk = (−1) i (−1) j (−1) k , y, z) are coordinates of the observation points, and ( i , j , k ) are corner coordinates of the density cube. In addition, the error between the numerical and synthetic data was calculated using the L 2 norm, Equation (18).
where n m is the number of measurement points, g num , and g synth are the gravity values computed numerically (FEM) and synthetically (Equation (17)), respectively. The results from using the 3D numerical model are presented and discussed in the next section.

RESULTS AND DISCUSSION
This section presents the results obtained from the numerical formulation of Section 4. The section contains the discussion of the effects of mesh fineness and the position of the pole.  Figure 4 shows the numerical results for different element sizes and three positions in the density zone (vertical and horizontal profiles, Figure 3A). Figure 4A,B shows gravity results for the surface of the numerical model at y = 16 m and y = 32 m, respectively. The figures illustrate that there is a significant effect of the element size on the numerical results. The performance of the numerical model improves considerably as the mesh is refined. Figure 4C displays the numerical results of the vertical profile. It shows an influence of the element size on the results but to a smaller extent than that of the surface horizontal profiles. Figure 4D shows the degree of agreement between the analytical results and the formulated numerical model outcomes for the vertical and horizontal profiles. The values of the L 2 error are calculated from Equation (18). The results show a significant convergence between the numerical and analytical solutions as the mesh is refined. It is also shown that the results of the surface horizontal profiles, in comparison with that of the vertical profile, are more sensitive to the element size. For the vertical profile, Figure 4 displays that an element size smaller than 4 m does not add an improvement to the accuracy of the results for the vertical profile, while for the horizontal profiles, there is an improvement when the element size decreases to 2.3 m. The gravity results shown in Figure 4A to C and the L 2 errors shown in Figure 4D indicate that the performance of the formulated numerical model is very good as compared with the analytical results (e.g., an average difference of 8 to 10 Gal between the numerical and analytical results for the element length of 2.3 m). In order to obtain a better insight into the distribution of the error over different areas of the model, the difference between the analytical and numerical gravity values at each node was calculated and plotted in Figure 5 for all three profiles (vertical and horizontal). Figure 5 shows that the error is mainly associated with areas close to the model boundaries. The reason could be due to the order of the shape functions used in this work. Since the type of the decay of gravity is 1∕r 2 , adding additional nodes to the boundary infinite elements could improve the results.

Effect of the pole position
The effect of the pole position (i.e., length of the infinite element) was studied by simulating models with variable lengths of the boundary infinite elements. The numerical results are presented in Figure 6 for the vertical and horizontal profiles in the model.  Figure 6C shows the E L 2 error computed from Equation (18) for different lengths of infinite elements, ranging from 5.0 to 45 m. The results are for three different profiles (vertical and horizontal), and there is an interesting change in the trends between the results for different positions. For small values of infinite element length (e.g., 5.0 m), the magnitude of the error between the numerical and synthetic data is significant. The error decreases with increasing infinite element length until a length of 17.5 m at which an optimum behaviour of the numerical model is achieved. This is due to the fact that practically, the error at different positions of the model converges to the same value. At a length of 20 m, the error in the vertical profile bounces slightly; however, there is a little improvement to the results for the surface horizontal profiles. After that, the error at all positions starts to increase.
It is worth noting that the choice of the pole position (i.e., infinite element length) needs to be done accurately to avoid unrealistic results. For the numerical domain of Figure 3A, the length of 17.5 m was chosen for the infinite elements since there is an overall convergence of the results in different positions. However, for a case where the surface of the model is the main focus, a length of 20 m could be chosen to obtain slightly better results at the desired position.
The next section demonstrates the application of this approach to a practical subsurface civil engineering application.

SUBSURFACE CIVIL ENGINEERING APPLICATION OF THE FORMULATED INFINITE/FINITE ELEMENT MODEL
This section presents the application of the formulated method to a subsurface civil engineering related application using a practical example. The formulated model in this example reproduces a set of field gravity data shown in Figure 7A after some simplifications of the ground density, which is explained in more detail below. The test area from which the gravity data were observed ( Figure 7A) is located at the University of Birmingham UK (Edgbaston) campus. There is a buried multi-utility tunnel with dimensions of 2×2m in cross-section, passing through the survey area. The tunnel is shown by a solid straight line in Figure 7A. Figure 7A also highlights two lines of gravity observations running perpendicular to the tunnel axis (lines 1 and 2). The survey was conducted using a Scintrex CG5, which is based on a mass on a spring instrument. It has a practical sensitivity of approximately 8 to 10 Gal dependent on how well the unwanted signals such as environmental noise (see Boddice et al 38 ) can be removed from the measurements leaving the signal of interest. The survey was carried out on a grid, taking at least three measurements per observation point subject to the variations between each measurement. Each measurement took 30 seconds in order to average out the high-frequency noise. Figure 7B presents relative values of gravity along the survey lines 1 and 2. Each survey line is 12 m long with zero positioned at the left. The points located at 6 m are directly above the multi-utility tunnel. It should be mentioned that accurate information about the distribution of the ground density in the survey area is not available. However, gravity values along lines 1 and 2 illustrate that the ground density is not distributed uniformly. The middle zone is generally less dense than the zones located at the ends of the survey lines, which is as expected given that the multi-utility tunnel is at the centre of the survey line representing an air-filled void surrounded by soil.  To simplify the problem for this example, an average set of field gravity data, symmetric about the tunnel centreline, was used for the comparison with the numerical prediction of the gravity. Furthermore, the field gravity data, representing relative values, were normalised using Equation (19) in order to be comparable to the numerical results. Figure 7C shows the normalised gravity data along lines 1 and 2, as well as the curve of the average values.
where g i and g i,norm are the gravity and normalised gravity values at observation point i, and g min and g max are the minimum and maximum values of gravity in each observation line.
In order to obtain the symmetric average field gravity data, the numerical domain was divided into the middle and the two end zones, as shown in Figure 8A. The two end zones had the same values of density (1900 kg/m 3 ) but different from that of the middle zone (1800 kg/m 3 ). Moreover, the density within each zone was considered uniform. The correct density for these zones to produce the average field gravity of Figure 7C was found in a trial and error process using an analytical model based on Equation (17). The model was identical to that of Figure 8A, which had dimensions of 32m×25m×20m (x, y, and z directions).
To reproduce the average field gravity data, a density of 1900 kg/m 3 for the soil end zones, 1800 kg/m 3 for the soil middle zone, 0 for the tunnel opening, and 2500 kg/m 3 for the tunnel lining were used. It is worth mentioning that the tunnel lining was assigned a thickness of 0.25 m. Furthermore, density values were given to the model in terms of contrast. In this way, a density reference of 1900 kg/m 3 was subtracted from each density value, which gave a 0 density contrast to the soil end zones, a −100 density contrast to the soil in the middle zone, a 600 density contrast to the tunnel lining, and a −1900 density contrast to the void created by the tunnel, as shown in Figure 8A.
After obtaining the density distribution, a numerical model was formulated for the domain shown in Figure 8A. The model had a total of 30 240 elements (33 292 nodes). A layer of infinite elements was added to the domain for modelling the boundary. The lengths were L Mx and L My in the x and y directions, respectively, and 0.1L Mz in the z direction, where subscripts Mx, My, and Mz denote the model length in the directions x, y, and z (i.e., 32 m in the x-direction, 25 m in the y-direction, and 2 m in the z-direction). The normalised gravity values predicted by the formulated numerical model as well as the average field values are plotted in Figure 8B. The results show a very good agreement, which demonstrates the capability of infinite elements for modelling the boundary of a finite element domain and reproducing accurate gravity values (compared with the average field data in Figure 8B). Figure 8C shows the residual (in Gal) between the predicted and average field gravity values.
It is worth noting that the modelled domain consists of zones of different density contrast with considerable changes from a high density material (concrete) to a void (zero density). There is also a significant difference between the dimensions of the anomalous zones (e.g. a 0.25 m wide tunnel wall in a 32×25×20 m domain). Such a domain needs relatively fine meshing within the zones of density contrast in order to obtain accurate results. The utilisation of a truncation method for modelling the boundaries could either lead to a high cost of computation or a low quality of results. The high cost is because of producing an extremely big number of elements due to having a relatively fine mesh within the zone of density contrast and adding many layers of zero density around the domain in order to obtain satisfactory results. The low quality of the results is when fewer elements are used for meshing the zone of density contrast, or boundaries are placed insufficiently far from the contrast zone. The use of the infinite element boundary eliminates the high cost of the computations with good levels of accuracy. Furthermore, more complicated situations can be modelled due to focusing on the zone of the density contrast without spending much storage and computational capacity with the added layers of zero density to model the boundary. It is worth noting that Haji et al 39 simulated 2D numerical models of truncation and infinite element boundaries for a simple gravity modelling problem. The zone of the density contrast in their model was identical (in terms of size and meshing) for both boundary methods. For an L 2 error of approximately 0.20, they showed that the infinite element model needed a total of 2704 elements, whereas the truncation methods had a total of 40 000 elements.
As explained earlier in this section, average symmetric values of field gravity produced by a set of uniform densities were used in this example. This is obviously a considerable simplification of the problem; however, the aim of this example is to show the ability and application of the numerical infinite element boundary in buried infrastructure problems rather than finding the actual density distribution.

CONCLUSIONS
This paper presented a detailed and comprehensive procedure to develop the formulation of an innovative infinite element boundary in 3D to model gravity data using finite element analysis. Current work in this field provides limited knowledge and is somewhat less developed particularly for the application of forward modelling of gravity. The significant advantage of this type of boundary condition over the existing methods (e.g. truncation) is that the domain size of the problem can be reduced to a great extent due to excluding the zones of zero density around the anomalous density domain. This contributes significantly to reducing the computational time and effort of the numerical modelling. Additionally, more complicated cases can be modelled using the infinite element boundary since the elements that would be used in the truncation methods for modelling the boundary (the zone of zero density around that of the density contrast) will be available for modelling complex zones within the density contrast.
Interpolation functions were derived in the paper for the reproduction of infinity of the gravitational potential at the unbounded boundaries of the domain. The effects of the mesh fineness and the position of the pole of the infinite elements (i.e., infinite element length) on the accuracy of the developed formulation were examined.
The results showed a significant effect of both the mesh fineness and the pole position on the performance of the infinite element boundary in the numerical model. Having a reference to accurately choose a length for the infinite elements was shown to be crucial. The reference could be engineering knowledge, background information about the problem, a known decay pattern, or an analytical solution.
Furthermore, a subsurface civil engineering-related problem was analysed using the formulated infinite element boundary. The results showed a high performance of the model for predicting gravity. This model could be implemented in an inversion process to reconstruct density contrast within the ground. Once used in an inversion process, the proposed infinite element formulation can overcome the problems associated with the currently used boundary conditions in numerical analyses of gravity by significantly minimising the computational time and reducing storage space whilst maintaining a very good level of accuracy.