Modelling the Gravitational Effects of Random Underground Density Variations

A mathematical model for small-scale spatial variations in gravity above the Earth’s surface is presented. Gravity variations are treated as a Gaussian random process arising from underground density variations which are assumed to be a Gaussian random process. Expressions for two-point spatial statistics are calculated for both the vertical component of gravity and the vertical gradient of the vertical component. Results are given for two models of density variations: a delta-correlated model and a fractal model. The effect of an outer scale in the fractal model is investigated. It is shown how the results can be used to numerically generate realisations of gravity variations with fractal properties. Such numerical modelling could be useful for investigating the feasibility of using gravity surveys to locate and characterise underground structures; this is explored through the simple example of a tunnel detection scenario.


Introduction
This paper is concerned with gravitational clutter in gravity surveys. Gravity surveys measure small variations in the Earth's gravity, and these variations can provide information about the nature and structure of the Earth and the underground environment. Surveys are carried out in a wide variety of scenarios: on land, in bore holes, on or under the sea surface and from space (Hinze et al. 2013). They cover a large range of scales, from small-scale surveys used to detect sinkholes (Styles et al. 2005) (Panisova et al. 2016), to large-scale geophysical prospecting for minerals (Hinze et al. 2013).
The work reported here was motivated by the desire to develop a simple model of clutter that can be used in computer simulations of small-scale gravity surveys. Of particular interest are surveys used to detect underground objects such as sinkholes, or pipes and other artificial structures. The ground in which the anomaly of interest is embedded may have a density which varies from place to place. These density variations will produce their own gravitational signal, which can be treated as random clutter (Boddice et al. 2019). The characteristics of the clutter may determine whether particular underground structures can be detected and identified. In this paper, results are given for both the vertical component of gravity and the vertical gradient of the vertical component, with particular emphasis on the latter because of its relevance to the recent development of new, highly sensitive gradiometers based on cold-atom technology (Sorrentino et al. 2014). These new instruments may find application in near-surface geophysical monitoring. The gradiometer is particularly attractive for near-surface studies. Insensitivity to Earth tides and seismic noise, combined with high sensitivity, allow for more rapid and accurate surveys (Boddice et al. 2019;Kennedy et al. 2014). Hydrological studies, pertaining to groundwater or reservoirs, are one potential application area (Piccolroaz et al. 2015).
In Sect. 2, theoretical results for relevant statistical properties are derived. These results relate the statistics of the gravitational variations to the statistics of the underground density variations. In Sect. 3, a method is presented which enables computer generation of data having the correct statistics; this method can be used to generate either random density variations, which can then be used in a forward model to calculate the gravitational field, or to directly generate gravitational field data.
The approach taken is to model the clutter as a stationary Gaussian random process. Such a process is fully described by its autocorrelation function (or, equivalently, its power spectral density) (Middleton 1996). The task is then twofold: first, to relate the autocorrelation function of the gravity variations to the autocorrelation function of the underground density variations, and second, to produce a computer model which can generate Gaussian-distributed gravity data having the prescribed autocorrelation function. These tasks are covered in Sects. 2 and 3, respectively.
In Sect. 2.1, general expressions are given for the autocorrelation functions and power spectral densities of gravity and gravity gradients produced by homogeneous underground density variations; the detailed derivations are contained in an "Appendix". Then, in Sects. 2.2 and 2.3, specific results are produced for a number of different models of density correlation. Some of these results are in the existing literature, but some are new to the best of the author's knowledge; the latter include the gravity gradient results and the results for a power-law model with an outer scale.
One of the models of density variations is the power-law spectrum. Crustal density variations in borehole samples have been shown to have spectra that display a powerlaw dependence on spatial frequency (Pilkington and Todoeschuck 2004;Bansal and Dimri 2010). Such power-law behaviour implies fractal, multi-scale spatial properties. There is also evidence for power law behaviour of density variations in surface samples and gravity surveys (Miranda et al. 2015). Section 4 presents an example of the use of computer-generated data to investigate the problem of the detection of an underground tunnel in clutter, when the clutter has a power-law spectrum.

Gravity Correlation Functions
Random processes can be characterised by average quantities, for example, mean and variance. For spatially varying quantities, the degree of correlation between values at two points is often an important metric. Indeed, for a Gaussian process, it provides a complete description (Middleton 1996).
Consider a gravity measurement carried out at a point {x,y,z}. The z (vertical) component of the gravitational field produced by random density variations σ r , z within a volume V is given by Hinze et al. (2013) where r is a vector in the horizontal plane, that is, r {x, y} and G is the gravitational constant. It will be assumed that the data are, if necessary, corrected for the Earth's curvature, so that z can be considered as the vertical component of a Cartesian coordinate system. The geometry is shown in Fig. 1. A gravity sensor is placed at a height z above the ground, the ground surface being at z 0 and lying in the x-y plane. A point below the ground of density σ contributes a gravitational force at the sensor, and the integral in Eq. (1) is over all of these contributions within the volume of integration. A similar expression applies to the gravity gradient, F zz , but with the kernel of Eq. (1) being replaced by its derivative with respect to z. In practice, a vertical gravity gradient is measured as a difference in gravity between two sensors at slightly different heights. For simplicity, this is approximated as a true gradient, that is, the limit of infinitesimal height separation. Now consider the spatial autocorrelation function of the gravity signal. This is given by a double volume integral. The integral for the autocorrelation function of the random gravity component is where the two outer integrals are over the infinite horizontal (x,y) plane and the random density variations descend to a depth of z a ; for some results the limit z a → ∞ will be taken. Here, the angled brackets denote ensemble averages, so the left-hand side of Eq. (2) is an average over the quantity inside the brackets. The average on the lefthand side is shown explicitly by the angled brackets, the corresponding average on the right-hand side is implicit in ρ σ r 1 , r 2 , z 1 , z 2 , the correlation function of the density variations σ . This correlation function is assumed to depend on the difference coordinates only, that is, the density variations are assumed to be a stationary process. Note that the convention used here is that the correlation function, ρ, has a subscript which shows which quantity it refers to; if there is no subscript, the expression is a generic one (the same convention is applied to the power spectral density S and the structure function D). Thus, the correlation function of gravity is written as It is clear that the correlation function of the z derivative of gravity can be found by differentiating Eq. (3) twice. That is Now, the correlation function of the density variations can be written as the inverse Fourier transform of a power spectral density S σ where the spatial wave vector is k k x , k y , k z . The approach that will be taken is to substitute Eq. (5) into Eq. (2), which yields a ninefold integral, and to evaluate the spatial integrals first, which can be done without including an explicit functional form for the power spectral density. The details of this calculation can be found in the "Appendix". For the case of z a → ∞, the resulting expression for the correlation function of F z is where x x 1 − x 2 and y y 1 − y 2 . It can be seen from Eq. (6) that the correlations in the horizontal plane depend on the difference coordinates, x and y, only, whereas the vertical correlations depend on the sum coordinate, z 1 + z 2 only. In principle, any desired power spectral density S σ can be used in Eq. (6). In practice, the resulting integrals may or may not prove tractable. It can be useful to have results for the two-dimensional power spectra of the gravity variations. These can, for example, be used to generate numerical realisations of random gravity variations, as in the simulations in Sect. 3. The two-dimensional spectra are Fourier transforms of the correlation functions with respect to the x and y coordinates only: the z coordinates remain untransformed. If the density variations are isotropic in the horizontal plane, and the limit of z a → ∞ is taken, then the twodimensional power spectral densities of the gravity and gravity gradient are given by [see "Appendix" Eqs. (A6) and (A7)] and where ω is the (radial) angular frequency. The derivation of Eqs. (6) and (7) is essentially the same as that of Naidu (1968); however, that paper has a multiplying factor of 4/π 2 rather than 1/(2π ) [Eqs. 32 and 33 in Naidu (1968), where 2d 1 z 1 + z 2 ]. The reason for this is twofold: First, Naidu (1968) uses a definition of the Fourier transform in which the factor of 2π is in the forward transform rather than inverse transform [compare with Eq. (5)], so to compare the expressions, one needs to multiply Eq. (6) by (2π ) 3 , which gives a multiplying factor of 4π 2 . Second, there is a minor error in the result given for one of the integrals, where a value of 2/π is given instead of the correct 2π [this error affects Eqs. 24 and 26 in Naidu (1968)]; this results in an erroneous division by π 4 , and thus a multiplying factor of 4/π 2 rather than the correct 4π 2 . This error also appears in the results of Maus (1999) (Eq. 7 in that paper), who referred to Naidu (1968), but the correct version is given by Naidu and Mathew (1998).
In the present paper, the following functional forms of S σ will be considered: a delta-correlated (spectrally flat) model and a power-law model. The former is the simplest possible model and is sometimes used in gravity modelling (Zhdanov and Cox 2013;Brown et al. 2016); the latter is based on the experimental data discussed earlier (Pilkington and Todoeschuck 2004;Bansal and Dimri 2010;Miranda et al. 2015).

Delta-Correlated Model
In forward modelling of gravity, it is often the case that the volume of interest is divided into a grid of small regions each of which is assigned a particular density. The simplest way to introduce random density variations into such a model is to assign an independent random density component to each region. In the limit as the grid size tends to zero this becomes a delta-correlated model. The power spectral density is given by the Fourier transform of the delta-correlated autocorrelation function, which results in a flat power spectral density From this it can be inferred that d 0 has units of kg/m 3/2 , that is, density times square root of volume. Substituting Eq. (9) into Eq. (6) leads to the following results for the gravity and gravity gradient autocorrelation functions (see Appendix, "Delta-Correlated Model" section, note the "Appendix" also has results for the autocorrelation function from a layer of finite thickness)

Power-Law Spectrum
A power-law spectral model has a power spectral density that is proportional to a given power of the spatial frequency over some range of that variable. The mathematically simplest form of power law spectrum is where ν is the exponent of the power-law and A is a multiplying factor that determines the strength of the density fluctuations. Here, ν is a dimensionless number, S has units of kg 2 m −3 and A has units of kg 2 m −(ν+3) . This model is spatially isotropic. One thing that is apparent with this model is that it has a singularity when the spatial frequency is zero. For some values of ν, the singularity causes integrals containing S σ to diverge. One way of avoiding this unphysical divergence is to modify the functional form of Eq. (12) to remove the singularity; this is equivalent to having an outer scale in the spatial domain that removes large-scale correlations. This can be done by using the following spectrum where k 0 is a constant: an outer scale parameter. Another way of avoiding the singularity is to use the structure function of the gravity variations rather than the correlation function. This latter approach has the advantage of being more mathematically tractable, but the disadvantage of only working for a certain range of values of the exponent ν. The structure function is the variance of the difference in the density between two points separated by distance r, and is thus zero for r = 0. Since one is usually only interested in data taken over a finite region, the divergence of the structure function for large spatial separations can be ignored. When the correlation function, here denoted ρ, is not divergent at r = 0, it is easy to show that the structure function is given by In situations where the integral for ρ diverges as r → 0 [as Eq. (5) may], the structure function can be written by taking the difference in Eq. (14) inside the integral sign, as long as the kernel exists in the limit r → 0. Note that, in the case of an isotropic spectrum, the kernel in Eq. (5) can be re-written in spherical polar coordinates and reduced to a single integral This is a useful result that links the variance of the density variations within a given region to the spectral model. For the model of Eq. (12) the integral can be evaluated to give where G is the Euler gamma function. Note that Eq. (16) is only valid for 3 < ν < 5. The upper limit is due to the divergence as k goes to zero, and the lower limit is due to another divergence as k goes to infinity [i.e. the lack of an inner scale in Eq. (12)]. According to Bansal and Dimri (2010), experimental data suggest that models with values of ν between 2.86 and 3.97 are appropriate for rocks in continental crusts. Note that the lower value of this range is less than the lower limit for convergence of Eq. (15). However, this divergence of Eq. (15) does not mean that the pure power-law model cannot be used for values of ν that are less than 3. As shall be shown later, the correlation functions for F z and F zz only have divergences associated with the lack of an outer scale (i.e. at small k); the lack of an inner scale does not produce divergent integrals. This is because gravity acts like a low-pass spatial filter which suppresses high spatial frequencies, which correspond to small-scale density variations. It does mean, however, that care is needed in relating the parameters of the models of Eqs. (12) and (13) to experimental rock density data: one cannot always simply fit to the measured variance of rock density. In fact, the same issue arises with the delta-correlated model of Eq. (9), which also suffers from an infinite density variance. Equations (7) and (8) can be used to find the two-dimensional power spectral densities of the gravity variations with the power-law models. Making use of the integral applied to Eq. (A23) in the "Appendix", one finds that the model of Eq. (13) gives and Here, B is the beta function and 2 F 1 is a hypergeometric function. For the case where there is no outer scale [model of Eq. (12)], these reduce to and The result of Eq. (19) is the same as that given by Maus (1999), apart from the erroneous multiplying factor discussed previously. It is shown in the "Appendix" that when Eq. (12) is substituted into Eq. (6), the resulting integral for the autocorrelation function diverges when ν > 1. In the corresponding case for the gravity gradient, the resulting integral diverges for ν > 3. However, the structure function can be used for larger values of the power-law exponent. The structure function for the vertical component of gravity can be expressed in terms of a hypergeometric function [see "Appendix", Eq. (A19)] where the reflection formula for the gamma function has also been employed (Abramowitz and Stegun 1970). This is valid for ν < 3: for larger ν the integral does not converge. By using the asymptotic expansion for the hypergeometric function (Abramowitz and Stegun 1970), one finds that for large r, the structure function goes as r ν−1 . The corresponding result for the gravity gradient is (see "Appendix"; note, the "Appendix" also contains results for a finite layer) which is valid for ν < 5. The asymptotic result for large r goes as r ν−3 . In Appendix "Power-Law Spectrum with Outer Scale" section it is shown that the model with an outer scale, given by Eq. (13), produces a result for the z component of gravity which can be reduced to a single integral The function B with three arguments is the incomplete beta function. It does not seem possible to express this integral in terms of known functions, but it can be evaluated numerically and has the advantage of being finite for ν > 3. By making k 0 sufficiently small, the outer scale can, if necessary, be made large enough to have no influence in the range of r values that are of interest. The corresponding result for the gravity gradient simply has an extra factor of k 2 multiplying the kernel of the integral. The structure function of the density variations, Eq. (16), is plotted in Fig. 2 on a log-log scale. In fact, the square root of the structure function is plotted, which gives a root mean square density variation in kg/m 3 . Lines are shown with three different values of the power-law exponent ν; the values of the magnitude parameter A were chosen to give the same density variation, of 36.4 kg/m 3 , at a separation of 0.1 m; note that this is merely a choice of convenience, because A is a free parameter in this model. The results in Fig. 2 show the standard deviation of the differences in density between a pair of points separated by the distance given on the horizontal axis. Densities become increasingly dissimilar as the separation between points increases. It can be seen that the power-law model gives a structure function which is a straight line when plotted on a log-log scale. This shows that variations in the density occur across all scale sizes, which is a characteristic feature of power-law models. The slopes of the lines in Fig. 2 are simply one half of the power-law exponent. Obviously, there is an unphysical aspect to this pure-law model: at large separations, the density variations become unrealistically large. The larger the power-law exponent, the more rapidly the unphysical densities are reached. This issue could be overcome by modifying the model to include an outer scale. However, in many circumstances, this modification Eo of the power-law model is not necessary, and the simplest form can be retained. If one is only considering a small-scale survey (say, of order 100 m), the most important density variations will be those on the scale size equal to or less than the dimensions of the region being surveyed, and the predictions of the model outside of this region can be ignored. Result for the square root of the structure function of the gravity gradient are plotted in Fig. 3 [using Eq. (22), and Eq. (14) applied to Eq. (11)]. Here, the separation distance r is confined to a horizontal plane at a distance of z 1 z 2 1 m above the ground surface. The units on the vertical axis are Eötvös units, that is, 10 −9 s −2 .
The solid curves are for the power-law model and the values of the parameters are those corresponding to the upper two lines in Fig. 2. The dashed line is for the delta-correlated model. It can be seen that the delta-correlated model gives a flat line at large separations, that is, separations much greater than the height of the measurement point above ground. This is a result of the lack of correlation in the underlying density variations. The power-law models, however, produce a structure function which continues to increase as separation distance increases, which is a result of the multi-scale nature of these models and the fact that they do not have a limiting outer scale. As mentioned previously, one would expect this increase to end at a sufficiently large separation, so in reality there will be an outer scale; however, for small-scale surveys, this outer scale may not be encountered, in which case the pure power-law model will be adequate.
Unlike the density structure function of Fig. 2, which is characterised by just a single slope, the power-law gravity gradient structure function has two distinct regions, each with a different slope. At large separations, the slope approaches the asymptotic limit of ν-3. At shorter separations, however, the structure function decreases more rapidly with decreasing separation. This corresponds to a more rapid reduction in small-scale spatial structure in the gravity gradient, which is a consequence of the spatial filtering effect of gravity: the fluctuations in gravity are coming primarily from structures near the surface, and deeper structures do not contribute.

Numerical Simulation
Gaussian random processes can be simulated by using the method of filtering of Gaussian noise (Jakeman and Ridley 2006). One approach to modelling gravitational clutter would be to generate realisations of underground density variations having the desired correlation properties and calculate the gravitational effects of these on the sensor using a forward gravity model. However, this has the disadvantage that a large three-dimensional underground volume may need to be modelled, which requires a large amount of computer memory and processing power. In the case of deltacorrelated density fluctuations, the requirement on computer memory may be reduced by treating the underground environment as a series of separate layers, the gravitational effects of which can be calculated separately and then summed. This works because in the delta-correlated model there are no correlations between different layers, so they can be treated independently of one another. For models that are not delta-correlated, such as the power-law models of Eqs. (12) and (13), this approach cannot be used. It would be possible to have a model that is delta-correlated in the z-direction, but correlated in the horizontal plane, and this could then be treated as an independent set of layers. However, a different approach, which will be used here, is to directly generate the random gravity, or gravity gradient, signals by the filtering of Gaussian noise. The one potential drawback to this is that it assumes the underground environment contains nothing apart from the random density variations. Of course, an underground structure of interest can be introduced to the model by adding its gravitational signal to that of the random density variations, but this results in a spurious clutter component coming from the region occupied by the structure of interest. So if, for example, one were to add the signal from an underground tunnel to the model, the modelled clutter would include spurious density variations coming from that part of the ground occupied by the tunnel, and so the magnitude of the clutter would be overestimated. If, however, the volume of the tunnel is only a small fraction of the total volume that contributes to the clutter, the errors introduced may be acceptably small.
Two different methods were used. The basic method of filtering of Gaussian random noise was used to generate results for the delta-correlated model, and the method of sub-harmonics (Frehlich 2000) was used for the power-law model. Both methods are based on the use of the discrete Fourier transform (DFT). In the basic method, a realisation of the gravity gradient is given by Here, the real part of the right-hand side is taken; the imaginary part yields a second, independent realisation. The gravity gradient is produced on an N x by N y sized grid in the horizontal plane, with spacing of x and y in the x and y directions, respectively. a and b are arrays of zero-mean Gaussian random numbers with variances given by The step sizes in spatial frequency space being given by ω x 2π /(N x x) and ω y 2π /(N y y). The method of sub-harmonics modifies the basic technique to more accurately include large-scale spatial variations, which are negligible in the deltacorrelated model but significant in the power-law model. A set of sub-harmonics of the following form are added to result (24) The sub-harmonics are powers of 3, so for the results given here, which use N p = 2, the 3rd and the 9th sub-harmonics are added. a and b, are Gaussian random numbers with Results showing an example of the spatial variation of the gravity gradient over an area 100 m by 100 m are given in Fig. 4 for both the delta-correlated model (using Eq. A9 from the "Appendix") and the power-law model (using Eq. 20). Both are sampled on a square grid with 0.25 m spacing. The number of points on the original grid was N x = N y = 1024, so the plots in Fig. 4 show a reduced section of the original grid. The size of the regions in Fig. 4 is that used for the simulations in Sect. 4, and the reason for the larger original grid size was to improve the representation of low spatial frequencies in the power-law model. For both cases the height of the sensor above ground was 1 m. The delta-correlated model used d 0 5 kg/m −3/2 and the powerlaw model used ν 3.5 and A = 100 kg 2 m −6.5 . The result from the delta-correlated model has the unstructured appearance that is characteristic of white noise. In fact, the only correlated structure is on a small scale, commensurate with the height of the sensor above ground. This small-scale structure is a consequence of the filtering effect of gravity, mentioned earlier. In contrast, the power-law result has structures on multiple scales; this includes scales larger than the region shown, which is why all the data values plotted are significantly less than the theoretical mean of zero: this region comprises part of a larger structure. Note that the models describe density variations of zero mean, so a constant underground density needs to be added in order to obtain realistic, non-negative density values.

Tunnel Detection in Clutter
In this section, the problem of detection of the presence of a tunnel in random clutter is presented as a simple example of how computer-generated clutter can be used in a model of a gravity survey. The scenario that is modelled is one in which a gravity gradiometer is scanned in one dimension along a 25 m-long straight line, below which there may or may not be an underground tunnel. For simplicity, it is assumed that the following things are known a priori: (i) the tunnel (if present) is at right angles to the line of data; (ii) the tunnel is sufficiently long and uniform that its gravitational response matches that of an infinite cylinder; (iii) the tunnel centre is a known distance below ground; and (iv) the density contrast between the tunnel and its surroundings is 2,000 kg/m 3 (corresponding to an air-filled tunnel). The gravitational gradient from an infinite cylinder lying along the y direction, of radius r 0 , distance z 1 below the plane containing the sensor and density contrast η, is (Pan et al. 2017) where u x/z 1 and the centre of the tunnel is at x = 0. It can be seen that the shape of the gravity gradient produced by the tunnel, in the x direction, F zz (x), is a function of the depth only. Three example profiles are shown in Fig. 5 for different tunnel depths. The value of the peak signal is inversely proportional to the square of z 1 , and the tunnel profile becomes wider when the tunnel is deeper.
In the absence of knowledge of the tunnel depth, one approach would be to simply identify the largest peak in the line of data as a potential tunnel signal. The model for random clutter can be used to find the probability that such a peak was produced by the random density variations, i.e. the probability of false alarm. However, if there is prior knowledge of depth, a better approach is to apply a matched filter: if the depth is known, then the shape of the tunnel signal is also known and can be used as a matched filter. The matched filter processing involves taking the measured line of data and cross-correlating it with the known profile function (for simplicity, DFT-based circular correlation was used for this numerical calculation). This approach will not only give better rejection of random clutter; it will also give a value for the radius of the candidate tunnel. It is not difficult to show that if the profile of Eq. (28) is normalised to the form then the peak of the cross-correlation of Eq. (28) with Eq. (29) is simply ηr 2 0 (this can be shown by Fourier-transforming Eqs. (28) and (29) with respect to x, multiplying them and then inverse-transforming to obtain the cross-correlation function). Thus, the only unknown quantity, the tunnel radius r 0 , can be obtained from the value of the cross-correlation peak. Of course, this requires prior knowledge of the tunnel depth. Since the width of the signal peak increases with depth, errors will be introduced if the true tunnel depth differs from the assumed depth. More sophisticated approaches are possible, using a range of different matched filters each for a tunnel at a different depth, for example. However, here it will be assumed that the tunnel depth is known sufficiently well to obtain good results with a single matched filter.
Assuming there is no more than one tunnel along the 25 m line of data, the largest cross-correlation peak is identified as a candidate tunnel, its magnitude giving a corresponding radius for the tunnel. The question that will be addressed by the model is, what is the probability that a particular tunnel radius will be identified if only clutter is present? This will result in a "false alarm" probability as a function of tunnel size, that is, it will give some information about whether a tunnel of a given size could be readily detected in the presence of clutter. Therefore, in the simulations here, no tunnel is actually present. Each simulation run identifies a candidate tunnel with given radius and represents a false alarm. The probability of false alarm will depend on both the magnitude of the random density variations and the magnitude of intrinsic sensor noise. The results plotted in Fig. 6, show the fraction of candidate tunnels with radius above a given value. Here, it was assumed that the sensor is subject to white, additive Gaussian noise. The three different results correspond to no sensor noise and sensor noise with standard deviations of 2Eo and 5Eo. Obviously, the probability goes to unity as the tunnel radius goes to zero: a very small tunnel would give a very small response, and there is always enough clutter present to produce a cross-correlation peak which exceeds the expected response. The utility of the probability of false alarm plot is that it shows what size of tunnel could be reliably detected in the presence of this amount of clutter: by only accepting a candidate tunnel as genuine if its radius exceeds a certain threshold value, one can reduce false alarms to an acceptable level.
The method described in Sect. 3 was used to generate 100 realisations of power-law clutter with the same parameters that were used for the right-hand plot in Fig. 4. Each realisation was generated on a grid with 1024 by 1024 points and a grid spacing of 0.25 m. The central 25 m by 25 m region of the resulting 256 m by 256 m area was used for the simulated surveys. Although there are only 100 truly independent realisations, each realisation of clutter does contain multiple structures; thus, by using every row of data along the y-axis as a distinct survey, the effective number of realisations used was somewhat greater than 100. Sensor noise was included by adding random white Gaussian noise to the data. Then, for each row of data, the mean was subtracted and a straight line was fitted, the resulting linear trend also being subtracted; this was done to remove the effect of structures that are much larger than the tunnel signal, and thus allow a consistent cross-correlation to be obtained for each line of data. The largest cross-correlation peak was then identified and used to calculate a candidate tunnel radius. The values from all the rows were used to calculate a cumulative distribution P f . This is shown in Fig. 6 and is the probability that the candidate tunnel radius is less than a given value. Figure 6 gives results for the three different tunnel depths used in  Fig. 5, with 2Eo of sensor noise. As one would expect, the deeper the tunnel the larger it has to be in order to be detected. A tunnel of, for example, 20 cm radius could be detected at 1 m depth with a very low probability of false alarm; however, the same tunnel would be subject to nearly 10% false alarms at 2 m depth and 70% at 4 m depth. Figure 7 gives results for a single tunnel depth with three different levels of sensor noise, starting with the case of negligible noise. Even in the absence of sensor noise, a tunnel of radius, say, 0.1 m could not be readily detectable in this level of random clutter because setting a detection threshold at 0.1 m would result in around 80% of runs producing a false alarm. However, setting a detection threshold at a tunnel radius of 0.25 m would mean that tunnels of this size and larger would be detectable with a very small probability of false alarm, as long as the sensor noise was 2Eo or lower.

Conclusions
The methods presented in this paper allow one to model the effects of random underground density variations on gravity measurements. The theoretical results predict the statistical properties of the gravity signals from the statistical properties of the underground density variations. The numerical modelling methods allow one to generate realisations of random gravity fields, which can be used in a gravity survey model. The major underlying assumption is that the underground density variations, and therefore the gravity variations, constitute a Gaussian random process, that is, the two-point spatial statistics are jointly Gaussian. As with any modelling assumption, this needs to be verified by experimental measurement. This is a task that may be suited to a new generation of gravity gradiometers based on cold atoms, which are, in principle, capable of carrying out rapid surveys with high sensitivity. These instruments are still in the early stages of development, and significant future improvements in sensitivity are anticipated. In Sorrentino et al. (2014), a sensitivity of order 2Eo was reported, albeit with laboratory-based equipment and with a 2-h integration time. Gradiometers also have the advantage of requiring fewer corrections than gravimeters, although careful measurement and correction of surface topography will still be necessary.
The same need for experimental validation applies to the spectral models considered here. Although there is some existing evidence for power-law fractal behaviour, it should be noted that the surface-sample and gravity survey results in the literature (Pilkington and Todoeschuck 2004;Miranda et al. 2015) tend to consider variations over large scale sizes only, of order kilometres, so their relevance to the kind of smallscale survey considered here is uncertain. Borehole data does give statistics for smallscale density variations (sub-metre vertical separations), and there is evidence that power-law behaviour extends down to these small scales, see Bansal and Dimri (2010). However, this data uses vertical drill samples, most of which come from many hundreds of metres below the surface, so their relevance to near-surface features is not clear. Another feature of the spectral models presented here is the assumption of homogeneity in the spatial statistics of the density variations. This is a condition that may not be met in practice, particularly in the vertical direction, where a certain amount of stratification might be expected. However, for a gravity survey taking place in the horizontal plane, a lack of homogeneity in the vertical direction would be expected to be of minor importance: as long as density variations in horizontal layers are homogeneous, the assumption of overall homogeneity should not lead to large errors in the predicted horizontal correlations. Another relevant point regarding vertical stratification is that it means that data from vertical boreholes cannot be relied upon to infer horizontal statistics of density variations. If a borehole is going through varying strata, then the density variations seen in the vertical sample may be completely different in magnitude and nature from those within the individual, horizontal layers.
In conclusion, the results and techniques presented in this paper may be useful for modelling of future gravity surveys where clutter from random underground density variations is a significant factor. However, the assumptions of Gaussian statistics and the use of particular spectral models (such as the fractal model) would need to be verified by measurement campaigns before the results can be used with confidence.
The z integrals are just elementary exponential integrals and can be easily evaluated to give Eq. (6) in the limit of z a → ∞, that is, for a layer of infinite depth. Results can also be obtained for a horizontal layer of finite depth, in which case the z 1 and z 2 integrals evaluate to give where k k 2 x + k 2 y 1 2 .Assuming the density variations are isotropic in the horizontal plane, one can convert to a radial spatial coordinate r and change k x and k y in Eq. (6) to the polar coordinates k and θ , θ being the angle between k and r. The angular integral can then be carried out to give a J 0 Bessel function (Abramowitz and Stegun 1970), resulting in where r 2 x 2 + y 2 . The two-dimensional power spectral density of the gravity variations can be written as the Fourier transform of the correlation function Converting to polar coordinates and again integrating over the angular coordinate yields Substituting Eq. (A3) into Eq. (A5) and changing the order of integration gives The integrals with respect to r and k constitute a Fourier-Bessel transform, which obeys for any function F (Watson 1944). Comparison of Eqs. (A7) and (A6) leads to Eq. (7) in the main text. Differentiating Eq. (7), that is, using Eq. (4), via Eq. (A5), gives Eq. (8) for the gravity gradient.
For the correlation function, the k z integral in Eq. (6) is an elementary one and can be easily evaluated. Converting to polar coordinates and integrating over the angular coordinate gives ρ Fz (x, y) π d 2 0 G 2 ∞ 0 exp (−2z 0 k) J 0 (k r) dk.
infinite summation that results when Eq. (A16) is substituted into Eq. (A15) can be written as ∞ n 1 − r 2 (z 1 + z 2 ) 2 n Γ n + 1 2 − ν 2 Γ n + 1 − ν 2 (n!) 2 . (A18) Comparing this with the series expansion of a hypergeometric function (Abramowitz and Stegun 1970) leads to the following result A slight further simplification can be obtained using the reflection formula for the gamma function (Abramowitz and Stegun 1970) which gives Eq. (21) in the main text. Note that in the special case of ν 2, Eq. (21) can be expressed in terms of elementary functions Equation 22 for the gravity gradient can be derived by differentiating Eq. (A19) twice, using Eqs. (4) and (14), and then using the contiguity relations for the hypergeometric functions (Sneddon 1961). Alternatively, one can obtain the same result by differentiating Eq. (6).
A result can also be derived for the case of a horizontal layer of finite thickness, using Eq. (A2). The exponential terms in Eq. (A2) evaluate in the same way as before; however, there is the additional cosine term which contains k z . When the k z integral is performed for the power-law model, a further cosine transform results, which can be found in Erdélyi (1953). The structure functions can then be written in the following integral forms, which can be evaluated numerically.
Here K is a modified Bessel function of the second kind.

Power-Law Spectrum with Outer Scale
After substituting Eq. (13) into Eq. (6), one can convert to polar coordinates in the {k x ,k y } plane and carry out the angular integral to give where, as before, k 2 k 2 x + k 2 y . The k z integral can be found in Gradshteyn and Ryzhik (1980) in terms of a hypergeometric function, which can be written as an incomplete beta function (Abramowitz and Stegun 1970) and gives Eq. (23) in the main text.