3D density imaging using gravity and gravity gradient in the wavenumber domain and its application in the Decorah

Density imaging is a method that uses the inversion of the gravity and gravity gradient spectra in the wavenumber domain to create accurate 3D reconstructions of subsurface density distributions. This approach offers computational efficiency and rapid calculations. This research used preliminary inversions to examine the spectral characteristics of gravity and gravity gradient anomalies, as well as the resulting models, were scrutinized through preliminary inversions. 3D density imaging of gravity and gravity gradient was performed in the wavenumber domain using depth weighting on both noise-added and theoretical data, producing a density model that was consistent with the theoretical one. The technique was then used in the Decorah region of the United States, where 3D density imaging was performed and an examination of the properties of gravity and gravity gradient anomalies was conducted. The results showed where high-density Decorah complexes, low-density siliceous intrusive rocks, and high-density intrusive rock masses, were the distributed within the surrounding rock. Each of these provided comprehensive insights into the intrusive pathways to the rock mass. Thus, the appropriateness and effectiveness of the density imaging method were confirmed, supporting a deeper understanding of the structural division and geological evolution in the region.

However, the journey of geophysical inversion is shadowed by non-uniqueness, requiring the careful application of constraints.To surmount this challenge, a judicious amalgamation of diverse a priori information becomes imperative.The inversion process gains momentum by utilizing physical property ranges, geological understanding, and related geophysical techniques.This combination of sophisticated inversion techniques and contextual constraints emphasizes how important it is to achieve realistic representations of Earth's subterranean complexity.Some experts have improved the inversion technique based on elements such as adaptive pruning L-curve technique and hybrid PCG-bat algorithm, but they are mainly used for two-dimensional inversion.The research on three-dimensional physical property inversion is imminent [9][10][11] .In this context, density imaging in the wavenumber domain becomes a powerful tool that allows one to directly calculate the sub-surface density distribution from gravity anomalies.First of all, the wavenumber domain gravity orthogonal computation can significantly reduce storage space and increase computational speed, by converting complex pleated product computation in the spatial domain to a simple product computation in the wavenumber domain.This effectively increases computational efficiency.Previous studies have been conducted for many models, and wavenumber domain gravity normalization is a better-developed technical method [12][13][14] .While it is not dependent on a priori models, the wavenumber domain density imaging method exploits the advantages of wavenumber domain orthogonal computation.A wide range of approaches has been produced in the research history.For instance, the Cribb imaging technique uses the Fourier transform to reveal density distribution derived from vertical derivatives of observed gravity anomalies 15 , whereas the density equivalent distribution method proposed by Kobrunov and Varfolomeev, operates in the wavenumber domain and inverts density distribution grounded in gravity data spectra 16 .Fedi presented the DEXP approach, which uses extreme points to forecast the depth of anomalous structures 17 .A 3D correlation imaging technique was presented by Guo et al., that provides high spatial resolution for geological body occurrences and residual mass distribution 18 .Priezzhev's method leverages rapid computation in the wavenumber domain, combining a priori information with gravity field data 19,20 .Kobrunov's iterative inversion approach, based on functional representation, expeditiously produces sub-surface density and structural models 21 .Using depth scaling factors, Cui and Guo recently enhanced the wavenumber domain density imaging technique to produce high-resolution density models 22 .Previous studies on gravity and gravity gradient inversion have demonstrated the effectiveness of spatial domain inversion methods, however, they come with a cost in terms of computational time and space requirements.Meanwhile, frequency domain inversion methods struggle with limited depth resolution due to the rapid decay of gravity field spectra.Moreover, there is a clear research vacuum about the spectral characteristics of gravity gradient data.
Our study endeavours to address these voids by examining the spatial anomalies and spectral characteristics of gravity and gravity gradient data in the wavenumber domain of the theoretical model.Drawing on prior knowledge from gravity and gravity gradient spectra, we incorporate a depth weighting factor into the wavenumber domain density imaging method.This innovation aims to improve depth resolution, with the use of wavenumber domain gravity field inversion's inherent computational efficiency.We also extend this novel methodology to gravity gradient data, taking advantage of its unmatched lateral resolution for density imaging.This study presents tangible outcomes rather than staying in the theoretical discourse domain.Our findings highlight the achievement of three-dimensional density models with high precision, that are obtained from the inversion of both gravity and gravity gradient data.We move beyond theory and apply our method to realworld data from the Decorah area in the United States, utilizing measured airborne gravity gradient data.The successful creation of a high-precision three-dimensional density model of subterranean space validates the effectiveness of our approach and increases its utility and applicability.In the pages that follow, we embark on a rigorous journey, delving into the opportunities and challenges that converge as we strive for high-density imaging precision.Through the utilization of high-precision gravity and gravity gradient data inversion, we aim to further explore the geological enigma of the Earth.Our method not only resolves existing questions but also provides the key to unlocking unprecedented depths of comprehension.As we navigate this scientific frontier, our vision encompasses revolutionizing resource management, reshaping exploration paradigms, and deepening our understanding of Earth's hidden complexities.

Gravity and gradient forward modelling theory
The gravity anomaly formula utilizes the rectangular prism models as delineated by Blakely 23 (Fig. 1a), and is confined to the limits of x 1 ≤ x ≤ x 2 , y 1 ≤ y ≤ y 2 and z 1 ≤ z ≤ z 2 , as reported by Plouff 24 : where, g z is the gravity anomaly, γ is 6.67 × 10 -11 N•m 2 /kg 2 , which is the gravitation constant, ρ is the density.
The wavenumber domain spectrum of both gravity and gravity gradient is derived via Fast Fourier Transform (FFT), thereby yielding the gravity anomaly spectrum, as previously demonstrated 19 . (1) where T pq p, q = x, y, z is each component of the gradient whose second-order derivatives satisfy Laplace's equation, T xx + T yy + T zz = 0 , and with T xy = T yx , T xz = T zx , T yz = T zy .Due to these two facts, only five com- ponents of the gradient are independent.Each gravity gradient component responds uniquely to the size, shape, and thickness of the density anomaly, providing a wide range of constraints in the interpretation process.The gravity and its gradient are converted to the wavenumber domain to obtain its spectrum and the spectrum of the gravity anomaly g z is obtained according to the Fourier transform theory: where G k x , k y represents the spectrum of the gravity anomaly.k x and k y are the wavenumbers along the x-axis and y-axis, respectively.The gravity gradient formula in the wavenumber domain calculated is following 25 :

3D density imaging theory of gravity and gravity gradient in the wavenumber domain
According to Priezzhev the density imaging formula in the wavenumber domain reveals the following 19,20,22 : The constant n ( 1 < n < 10 ) reflects the filter resolution, with larger values indicating higher resolution.z is the inversion depth, k is the wavenumber, and ρ is the density of the mass body at a given point.Based on the conversion relationship between the gravity gradient anomaly spectrum and gravity anomaly spectrum 25 , the formula for inverting density distribution according to the gravity gradient anomaly spectrum can be derived through Formula (2), such as

Forward modelling of gravity and gravity gradient of rectangular prism model
We simulated Model I to verify the relationship between gravity, gravity gradient, and the model.The anomalies in gravity and gravity gradient were linked to the edge, boundary, angle, and mass centre of the anomalous mass 26 (Fig. 2).The boundaries were delineated by the zero value of T xx and T yy , which corresponded to the east-west and north-south variations of gravity, respectively.The corner point was determined by the extreme value of T xy .The anomalous axis in the north-south and east-west directions were represented by T xz and T yz , respectively, with their extreme values indicating the boundaries.The abnormal centre was indicated by the extreme value of g z and T zz with, T zz showing higher resolution than g and its zero-value corresponding to the boundary 27 .
(2) The gravity anomaly spectrum displayed rapid attenuation and fluctuation, while the gravity gradient spectrum displayed slower attenuation, pronounced periodicity, and a discernible response in the k x and k y directions.The spectrum characteristics of the gravity gradient were closely linked to the orientation of each component, with T xx and T xy displaying rapid attenuation and unsuitability for further inversion.

Analogous spectral features of gravity anomaly and gradient amplitude spectrum and determination of buried depth
The distributional attributes of the anomaly amplitude are revealed by the amplitude spectrum, which captures the frequency-dependent fluctuations in the gravity anomaly amplitude spectrum.The gravity anomaly and gradient amplitude spectrum exhibit analogous spectral features, which decay rapidly, and the wavenumber corresponding to the zero value varies periodically 28 (Fig. 3).The fluctuation-type amplitude spectrum curve is mainly related to the horizontal width of the geologic body, based on the characteristics of the amplitude spectrum curve that provide a preliminary judgment of the morphology of the anomalous body.The spacing of wavenumbers between the amplitude's zero-value points is closely related to the horizontal width of the body.The rectangular prism model can be inverted and the half-width in the horizontal direction can be derived from the wave value at the first minimum of the spectrum.
where B is the width of the body, π is the constant, which is 3.14, and k is the spacing of wavenumbers between the zero-value points of the amplitude.The gravity anomaly spectra of 3D and 2D bodies are equivalent 29 , allowing for the extraction of profiles at k y = 0 for both the gravity anomaly and gradient amplitude spectra.Due to the small gradient amplitude spectrum, two profiles were required.The first minimum point of the amplitude spectrum corresponds to a wave value of 0.247 (Fig. 4), resulting in a calculated half-width of 12.7 km, which closely matches the actual width of 12 km.The buried depth of a model can be determined by calculating the slope of its average radial logarithmic power spectrum, obtained by logarithmically transforming the average radial power spectrum 30 .where h t is depth, E(r) is the average radial power spectrum, and r is the radial frequency.

Noise tolerance evaluation of density imaging method using gravity gradient data
To test the density imaging method's noise tolerance, we added 10% random noise and calculated the theoretical forward gravity and gravity gradient anomalies of Model II (Fig. 1b).It is commonly known that the gravity anomaly was unable to accurately represent the model boundary; however, the gravity gradient anomalies aligned well with the inflexion point or model boundary, with T zz providing a clear outline of the model position.The presence of noise produced fuzzy characteristic lines for the gravity and gravity gradient data but had no discernible effect on the fundamental properties of the anomaly that corresponded to the model body (Fig. 5).Furthermore, we performed a comparative analysis of the computational speed of Model I and Model II.The spatial domain model orthogonalization is the analytical computation of the rectangular model and the wavenumber domain computation is done using Matlab 2016b compiled software.The computer parameters are as follows: 32 GB of the Random Access Memory (RAM) and a 12th Gen Intel(R) Core (TM) i7-12700 K 3.60 GHz the Central Processing Unit (CPU).The comparison results show that, particularly for large and complex models, the orthorectified modelling in the wavenumber domain is faster than modelling in the spatial domain, improving the inversion efficiency (Table 1).

Iterative inversion based on depth weighting
The proposed method's effectiveness and noise immunity were verified by inverting forward modelling and noiseadded data.The model's depth and width were initially evaluated using the gravity spectrum as a density imaging constraint.To enhance longitudinal resolution, we introduced Commer's depth weighting function 31  www.nature.com/scientificreports/spatial gradient weighting function that improves inversion results by introducing a priori depth information.This function has been widely used since it was proposed [32][33][34] .
where z is the centre depth, dz is the inverse domain, α is the empirical value that determines the weights at the near-surface and is generally taken as α = 0.001 to overcome the skin effect.z c1 , z c2 are the depth of the top and bottom of the model.d 1 and d 2 are the interface constraint factors for the upper and lower interfaces, respectively.The iterative inversion formula with the depth weighting function shows the following: where, ρ i is the density at i iteration, ρ i+1 is at i + 1 iteration, �ρ i is the correction.The subterranean domain to be inverted is partitioned into N horizontal layers, each layer is divided into m × n orthogonal prisms.
As illustrated in Fig. 6, the specific iterative process entails the following steps: Step 1 Begin by creating an initial density file, denoted as ρ 0 , with an initial value of zero.The dimensions of this file are configured to correspond with the range of gravity anomalies.
Step 2 Next, calculate the spectrum of observed gravity anomalies (or gravity gradients), and use Eq. ( 8) to infer depth-related data from the obtained spectrum.
Step 4 Combining the results of Step 2 with other geological information, initiate the first iteration of the density distribution, denoted as ρ 1 .This is accomplished by taking the depths obtained in the second step and combining them with other geologic information, combining the depth constraint function of Eq. ( 9), and then applying Eq. (10).
Step 5 Model the density distribution ρ 1 in the entire spatial domain orthogonal to the wavenumber domain.This facilitates the generation of its gravity anomaly spectrum (or gravity gradient anomaly spectrum), followed by a Fourier inverse transformation to obtain the computed gravity anomalies (or gravity gradient anomalies).
Step 6 After obtaining the computed values, compare them to the observed values.Calculate the root mean square (RMS) difference and evaluate the fitting results.If the difference in gravity anomaly fitting is less than 2.4 × 10 −8 (or gravity gradient fitting difference is less than 1 × 10 −3 ), terminate the iteration.Otherwise, consider the fitting difference as the observed gravity anomaly (or gravity gradient anomaly) and repeat steps two to six in the iteration process.
In this study, the model inversion depth is 20 km, the subsurface is divided into 40 rectangular layers, the grid of each layer is 128 × 128, and the inversion time is 15 s.Compared with traditional inversion methods such as least squares 35,36 , there is a very significant improvement in computational speed while both achieve computational accuracy.The top burial depth of the model is obtained as 0.2 km through the spectrum information, and the bottom burial depth is estimated as 15 km through the different parameters, and the addition of the depth weighting factor can control the depth of the model inversion to avoid the epidermal effect and the depth is too large.Iterating over each prism within the strata is accomplished as follows.This information-independent approach enables direct initialization of the density to zero.The specific iterative process is shown in Fig. 6.We present a novel density imaging methodology utilizing a flowchart approach to iteratively refine density models for gravity anomaly modelling.Notably, amplitude spectrum inversion results serve as key constraints driving the iterative process.We increment density models with each imaging outcome until the residual error between the model gravity forward value and the observed gravity anomaly satisfies predetermined accuracy thresholds, the RMS of g 1 is 2.4 × 10 -8 , while it is 1 × 10 −3 of T zz in the presented inverse problems, halting the iteration.This approach offers an efficient and effective means to refine density models for accurate gravity anomaly modelling.

Comparative analysis of density imaging results for Model I and Model II using gravity and gravity gradient data
The wavenumber domain computation of density imaging for Model I and Model II as tabulated in Table 2 is presented here.The resultant 3D density imaging outcomes based on g z and T zz are depicted in Fig. 7.The inver- sion effects of the data from both models are impeccable, and their corresponding theoretical positions have been ( 9) accurately determined.The horizontal positions of the inversion results of the data sets are delineated in Fig. 7, where the actual positions of the models are represented by black boxes, and the boundaries of the models are clearly illustrated.The vertical profile of the inversion result is evident that the position of the anomaly body in both data sets is in excellent agreement with the positions of the model.We note that the higher density in Model 4 leads to a slight tailing phenomenon in Model 5.

Analysis of the accuracy and noise resistance of 3D density imaging using gravity gradient data
The 3D density imaging results of gravity gradients, apart from g z and T zz exhibit a higher rate of false anomalies (Table 3, Fig. 7).This phenomenon is attributed to the derivative of the data, which leads to distortions in T xx  www.nature.com/scientificreports/and T xz , T yy and T yz , and T xy in the x and y direction.In contrast, the 3D density imaging outcomes of g z and T zz display a remarkable level of accuracy.Table 3 presents the residuals obtained from the comparison between modelled values of g z and T zz , both orthogonal to the theoretical values derived from density imaging through g z and T zz .The table demonstrates the exceptional accuracy of g z and T zz in both Model I and Model II.In Model I, the RMS of g z residual is 0.004, while the RMS for T zz is a mere 2.3 × 10 -4 .Similarly, in Model II, these values further improve to 7.03 × 10 -9 and 1.83 × 10 -5 respectively.When introducing a 10% random noise to the inversion results 37 , a slight enhancement is observed, yielding an RMS of g z residual of 7.09 × 10 -9 and 9.75 × 10 -5 for the two models.Despite this noise, the residuals comfortably adhere to the computational accuracy prerequisites, which mandate g z residual within 1 × 10 -4 and T zz residual within 1 × 10 -5 .We also tested the inversion results introducing a 20% and 50% random noise for Model II and the RMS of g z and T zz residuals were 7.29 × 10 -9 , 1.92 × 10 -4 of the former data, 8.57 × 10 -9 , and 4.84 × 10 -4 of the later, respectively.It follows that the results of 3D density imaging, even with the added noise, remain in accord with those of the original model, with the dense values predominantly situated within the model area.The assessment of the density imaging reveals a nearconsistent forward modelling value of the inverted density model and the residual error of the theoretical value, indicative of the robust noise resistance properties of the method, and, T zz is even better at it.Verification of the method's efficacy was performed in the Decorah region of the United States.

Full tensor gravity and gravity gradient measurements for 3D density imaging in the Decorah region
The Decorah region is rich in mineral resources and offers a promising target for exploration due to its vast sedimentary terrain.Bell Geospace has conducted meticulous gravity and gravity gradient measurements, capturing comprehensive full tensor data, by utilizing cutting-edge technologies.Three-dimensional density imaging has been generated to delineate sub-surface density distributions by using a wavenumber domain method.The efficacy of this approach has been validated through rigorous analyses.

Gravity and gravity gradient anomalies in the Decorah
Anomalies related to both gravity and gravity gradient are caused by variations in density within the Earth's interior.The latter is sensitive to variations in shallow density differences, whereas the former represents information from deep field sources.Comparison of the two allows for a comprehensive understanding of the density distribution within a specific area.Gaussian regional filtering, with a 2 km wavelength cutoff is used to identify the target anomalies and eliminate the impact of sedimentary caprock.The area under investigation, referred to as the Northeastern Iowa intrusive complex, is composed of a lower-density siliceous intrusion in the middle and the comparatively dense Decorah complex in the south (Fig. 8a).Within this area, several small, relatively dense mafic and siliceous rocks have intruded into the surrounding low-density rock 38,39 .

Analysis of gravity anomalies in the Decorah region
Between December 2012 and January 2013, Bell GeoSpace used an aerial gravity full tensor gradient measurement system to survey the subject area.A total of 94 EW survey lines, spaced at 400 m intervals, and 9 connecting lines, spaced at 4000 m intervals, were measured.The resulting gravity and gravity gradient anomalies are displayed in Fig. 9a,b, respectively.Examination of the g z values revealed a broad range of low gravity anomalies in the northeast region, corresponding to siliceous intrusion, as confirmed by geological mapping.Numerous high-value anomaly traps and belts were found in the western portion of the study area.These belts and traps

Physical property analysis of the Decorah region
The 3D density body horizon at a depth of 1250 m in the study region was determined using the Tikhonov regularization inversion method 38 (Fig. 8b).The horizon at the same depth was calculated using a novel technique presented in this study (Fig. 9c,d), with mostly consistent findings.An analysis of the Decorah's geologic map revealed the location of different intrusive rock bodies.On the east and north faces, two anomaly bodies with a relatively small range and high density were identified as mafic intrusion bodies; on the east and north sides, a large-scale SW-trending Decorah complex was found (Fig. 9e,f).The variation in sub-surface 3D density structure was clearly illustrated, revealing the intrusion process of the bodies from bottom to top.Together with the detected entities, the density model's iso-surfaces correspond to distinct geological structural units with minor anomalies.The inversion results were reliable, with residuals conforming to a normal distribution and having a standard deviation of 0.003 mGal for g z and 3.4 E for T zz (Fig. 9g,h).A superior high-precision density model of the research region was generated by utilizing a novel 3D density imaging technique in the wavenumber domain.This allowed for a comprehensive understanding of the geological distribution and mineral resources within this area.

Conclusion
In this study, we conducted an in-depth analysis of gravity and gravity gradient anomalies in the spatial and wavenumber domains.Leveraging the efficient computational capabilities in the wavenumber domain, we obtained a priori information on the model through amplitude spectrum inversion and introduced depth weighting to constrain the wavenumber-domain 3D density imaging method.This enabled us to quickly and effectively create an accurate 3D underground density model, which we then validated using measurements taken in realworld environments.Through a rigorous examination of the model results and spectrum characteristics, we conclude that g z and T zz demonstrate superior 3D density imaging inversion with a significant noise immunity and depth imaging effect.Additionally, the 3D density imaging inversion utilizing g z and T zz in the Decorah area of the United States unveiled the distribution of intrusive rocks and their intrusion path.This work highlights the effectiveness of the 3D density inversion method in the wavenumber domain and provides crucial insights for further research.Our future work will be focused on maximizing the superiority of gravity and gravity gradient through joint inversion and improving the accuracy of the sub-surface 3D density model.

Figure 1 .
Figure 1.Rectangular prism model, (a) is Model I, and (b) is Model II.

( 7 )Figure 2 .Figure 3 .Figure 4 .
Figure 2. Gravity and gravity gradient anomalies of Model I, (a-g) the anomalies of T xx , T xy , T xz , T yy , T yz , g z , T zz , the black line is the model.

2(r 2 −r 1 ),Figure 5 .
Figure 5. Gravity and gravity gradient anomalies of Model II (a-g) the anomalies of T xx , T xy , T xz , T yy , T yz , g z , T zz , the black line is the model.

Figure 6 .
Figure 6.Iterative density imaging with amplitude spectrum inversion constraints for gravity anomaly modelling.

Figure 7 .
Figure 7. g z/T zz 3D density imaging results of Model I and II.(a,b) are inversion results of Model I, g z and T zz results, respectively; (c-f) are results of Model II, (c,e) is the 3D body and slices of g z result, respectively.(d,f) are about T zz result. (g,h) the RMS of g z and T zz iterations.
www.nature.com/scientificreports/extend north-to-east, with a limited high-value anomaly range in the northern sector, and an expansive, highamplitude, high-value range from southwest to south, which is indicative of the Decorah Complex.Further analysis of the T zz map exposed several local anomalies, with the eastern region displaying a prominent high- value area and the southern portion exhibiting a distinct outline of the high-value area.Additional anomalies with minor amplitudes were found in flat regions.Cumulatively, analysis of the gravity and gravity gradient data facilitated the division of the study area into multiple structural units.Density imaging was conducted leveraging g z and T zz , with the sub-surface inversion at the depth of 0-8 km, layer interval of 150 m, 20 iterations, and sub-surface model grid points of (930 × 997 × 53).

Figure 8 .
Figure 8.The geological and density model of the Decorah, (a) is a Geologic map 38 , and (b) is a 3D density model 39 .

Figure 9 .
Figure 9. Measured gravity and gravity gradient anomalies in the Decorah.Blacklines are tectonic lines.(a) g z , (b) T zz , (c,e) are the inversion results of g z , (c) density voxel, (e) isosurfaces, pink marks ρ is 0.11, green marks ρ is − 0.1.(d,f)are the inversion results of T zz , (d) density voxel, (f) isosurfaces, pink marks ρ is 0.2, green marks ρ is − 0.18.(g,h)represent the residuals between the forward and observed values of the model obtained from g z and T zz inversion, respectively. https://doi.org/10.1038/s41598-023-49711-z

Table 1 .
Comparison of forward duration of gravity and gravity gradient in space domain and wavenumber domain.The more complex the model is, the more obvious the velocity advantage in the wavenumber domain is.

Table 3 .
Evaluation of the effect of gravity/gravity gradient 3D density imaging.The residual between the forward gravity anomaly of the theoretical model proposed in this article and the result obtained from model inversion in the wavenumber domain serves as a crucial evaluation metric.