Gravimetric inversion based on model exploration with growing source bodies (Growth) in diverse earth science disciplines

1.   Introduction

Gravimetry is a geophysical discipline that deals with acquisition, processing, and interpretation of gravity data. The data can be collected at the earth surface or in shallow subsurface, in mines, tunnels, boreholes, etc. (terrestrial gravimetry), above the earth surface by means of flying devices (airborne gravimetry), on the sea (shipborne or marine gravimetry), and even at altitudes of artificial satellites (extraterrestrial or satellite gravimetry). The observed gravity field (of the Earth or other planets and moons) is generated by the subsurface density distribution (of the mentioned bodies). Therefore, knowing the gravity field, the subsurface density structure can be studied. In addition, when temporal changes of the gravity field are observed, the temporal changes in the subsurface density distribution can be studied.

The direct gravimetric problem reads: Given a particular subsurface density distribution, determine the gravitational field respective to it. The solution of the direct problem is given by means of the Newton volumetric integral for gravitational potential or any other quantity derived from it, such as gravity anomaly/disturbance, geoidal height, deflection of the vertical, various gradient components, etc. (see e.g., ^[81,82] and references therein). The kernel of the given Newton volumetric integral is then respective to the given gravimetric quantity derived from the potential. The inverse gravimetric problem is formulated: Given the observed gravitational field, determine the density distribution generating it. The inverse problem is known to be non-unique and ill-posed. Gravitational potential itself is not measurable, but other quantities derived from it are. The inverse problem is solved either by iterative forward density modelling seeking the match (fit) between observed and model gravity data, or by direct inversion of the gravity data. The iterative forward modelling can be manual, semi-automatic, or automatic. The inverse problem is nonlinear. If the subsurface is partitioned into a set of bodies or structural elements with fixed geometry, seeking only the constant density of each body or element, then the inverse problem becomes linear. Loosely speaking, we can say the inverse problem is nonlinear in geometry and linear in density estimation.

The applicability of gravimetry across various earth science disciplines is diverse. In physical geodesy, it serves determining heights and height systems as well as the geoid adopted as the reference surface to reckon heights from. Gravimetry, integrated with other geophysical methods and constraints from geology and tectonics, to reduce the ambiguities of the gravimetric inverse problem, can be used in global or regional structural studies aiming at the determination of the structure and tectonic development of the earth crust or lithosphere, reaching as deep as the sub-lithospheric mantle. In regional or local studies, it can assist in prospection for natural resources such as mineral deposits, hydrocarbons, geothermal energy, reservoir monitoring, gas sequestration, radioactive waste storage, etc. In near-surface studies, it can assist archeological prospection or engineering geophysical applications. When observing gravity changes, the subsurface density changes due to various geodynamic phenomena can be studied, such as those due to glacial isostatic adjustment, due to magma migration in volcanic systems, ground water level, or amount changes in aquifers, etc.

In structural crustal or lithospheric studies, the gravimetry operates with signal at the level of several tens or hundreds of mGal (1 mGal = 10^-5 m/s²). In archeological prospection, engineering geophysics, hydrology, or studying volcanoes (volcano gravimetry), we deal with signal at the level of several tens or hundreds of μGal (1 μGal = 10^-8 m/s²). Then, we often speak of (3D or 4D) micro-gravimetry.

The Growth inversion approach (e.g., ^[16]) belongs to inversion methods that make no prior assumptions about the number and geometry of the sought sources, but seek sources in the form of agglomerates of populated cells of a subsurface partition. This can be done either by a linear approach, in which all cells of the subsurface partition are filled with unknown density contrasts by means of a constrained least squares adjustment (e.g., ^[7,37,45,48]), or by a nonlinear approach, seeking the shapes and sizes of subsurface source bodies or shapes of density discontinuity surfaces of prespecified density contrasts in an iterative growth (or shaping) process (e.g., ^{[44,50,53,59,71,72,80,87]}). The Growth inversion explores the model space and lets source bodies grow step-by-step by filling (with positive and negative prescribed incremental density contrasts) and aggregating the filled prismatic cells within the volumetric domain below the topographic surface during an iterative weighted adjustment process (e.g., ^[16]).

Mathematical constraints are usually employed by minimizing a function involving the parameters of the model (e.g., ^[32,46]), like the volume of the causative body (e.g., ^[45]), moment of inertia, e.g., ^[41], depth dependence of the solution ^[48], or distance, flatness, smoothness, and compactness (e.g., ^[3,10,37]). Incorporation of the mathematical constrains into the final solution is performed by means of the regularization approach, where the misfit functional is composed of two parts: the first one describing the fit of the model response to the observed data and the second describing additional mathematical properties of the regularized solution (e.g., ^[61,79]).

Other inversion approaches have also been developed, including the stochastic approach for several kinds of geophysical data (e.g., ^[74]), the spectral approach, fuzzy entropy, graph theory ^[85], and GPU parallel computing ^[31], to name a few. Optimization methods based on the Genetic Algorithm to solve nonlinear inverse problems (e.g., ^[1,47,49]) can also be employed on a variety of geophysical data. Gravity or gravity change data in volcano gravimetry or petroleum exploration have been inverted using the GA optimization (e.g., ^{[52,56,78,86]}). For structural studies in larger regions with deeper reach (lithospheric studies), the JIF3D software offers a framework for joint inversion of different types of geophysical data in 3D, including scalar and tensor gravity data. The LitMod3D package ^[36] serves 3D-integrated geophysical-petrological interactive modelling of the lithosphere and underlying upper mantle using a variety of input datasets: potential fields (gravity and magnetic), surface heat flow, elevation (isostasy), seismic, magnetotelluric, and geochemical. Another approach for crustal studies based on non-linear inversion of gravity and magnetic data was proposed by Prutkin, based on initial line segments adjustment and subsequent local corrections producing structural bodies and interfaces ^[67,69]. This approach was also applied to inverting and interpreting spatiotemporal gravity changes, namely those of the 2004–2005 Tenerife unrest ^[68].

Software tools for direct inversion (among many others) are, for instance: IVIS-3D ^[43], 3DINVER.M ^[38], 3GRAINS ^[75], GICUDA ^[30], DenInv3D ^[88], Grav3CH_inv ^[58], and open-source code SimPEG written in Python, e.g., ^[51]. Heavy software ^[42] is available for modelling gravity changes in hydrological applications. This list is not deemed to be exhaustive.

2.   History of Growth inversion methodology development

The Growth inversion methodology was originally developed for inverting point surface gravity data, the residual (regional trend corrected) complete Bouguer anomalies (CBA), in order to obtain structural density models of the subsurface. It has its roots in least-squares adjustment (LSA) inversion for the subsurface partitioned into prismatic cells (organized into layers of thickness growing with depth), assisted by correlation analysis (for separating signal and noise) and least-squares prediction (LSP), producing input gravity data on a regular grid ^{[11,12,13,54]}. The LSA inversion adopted double (weighted mixed) minimization, minimizing the misfit residuals and maximizing the smoothness of the model density distribution. Alternatively, instead of using a subsurface partition into prismatic cells, it could adopt a set of spheres, aiming at identifying subsurface cavities in karst areas ^[12]. For the gravitational effect of the prismatic cell, to populate the elements of the design matrix in the LSA problem formulation, it adopted the formula of Pick et al. ^[66].

The follow-on work abandoned the linear LSA inversion approach and focused on seeking the geometry of bodies of pre-specified (fixed) density contrasts (e.g., ^[4,60]) using iterative exploration of the model space, partitioned into prismatic cells, and growing bodies ^[14,16]. This approach is somewhat similar to the "bubbling" method ^[90] and the "open-reject-fill" method ^[71]. It was inspired by the latter. However, unlike the two quoted approaches, the one introduced by Camacho ^[14,16] was a 3D method, could work with irregularly spaced data given on the irregular topographic surface, did not need prior seeds for the growing bodies, could grow bodies of both positive and negative density contrasts, could deal with non-uniform partition (bigger prism at deeper placements), and could simultaneously determine (adjust) a regional trend (or additional trend) in the input gravity data that was not removed (or not properly removed) beforehand.

The size and resolution of the subsurface partition was automatically preconditioned by the extent and irregular spacing of the input gravity data, as well as by the assumed capacity (RAM) of the computer. The gravitational attraction of a given cell of the subsurface partition at a given gravity point (station) is given by the formula of Pick et al. ^[66]. The attractions of individual cells at individual gravity stations form the design matrix of the inverse problem formulation. Both positive and negative growing bodies are formed with a prescribed density contrast used in populating new cells in an iterative process. The proportionality between the observed gravity and the gravity respective to the model growing throughout the iterations, by adding more and more populated cells, is controlled by a scaling factor. The iterative procedure adds only one filled cell in each iteration step. It seeks the cell to be filled by a model exploration process aiming to improve the data fit.

Actually, the adjustment in each iteration step is a weighted mixed adjustment, in which the minimization of data misfit is combined with the minimization of the model complexity and scatter-ness (maximization of model compactness or smoothness), mediated by minimizing the total volume of populated cells. The parameter called balance factor weighs the balance between data fit and model compactness. Throughout the iterative process the model expands in positive and negative populated cells forming the growing bodies and the scaling factor decreases. When it effectively reaches the value 1, the process stops, producing the final model and the final trend in gravity anomalies. For a detailed description of the Growth inversion method with due mathematical apparatus, as well as its implementation in the form of the GROWTH software package, we refer the reader to ^[16].

The GROWTH implementation can initiate the iterative inversion process either from a null density model or from a starting density model. It can seek either homogenous bodies (of both positive and negative density contrasts), or heterogenous source bodies. It also allows optional stratification in the sought model. The GROWTH package is composed of three programs: (1) PARAM, which carries out the selection of parameters, options, and unknowns; (2) GROWTH, the main program which performs the gravity inversion; and (3) SECTIONS that allows visualization of the inversion results. Input and output files for and of these three programs, and the values of the parameters they offer as default or require as input, are described in due detail in ^[16].

The Growth approach to CBA data inversion was modified from growing source bodies to growing morphology (i.e., shaping of horizontal density discontinuity surfaces ^[5]). In 2011, the Growth approach progressed from its GROWTH implementation to GROWTH-2 ^[19]. The GROWTH-2 package is composed of three programs: (1) GRID3D is used to generate a 3D partition of the subsurface volume into parallelepiped elements, (2) GROWTH is used to perform the inversion routine and to obtain a 3D anomalous density model, and (3) VIEW is used for visual representation of the input data, of the inversion model, and of data misfit (the residuals). Compared to the predecessor package, GROWTH-2 provides several novelties: (a) a Graphical User Interface (GUI), (b) an optional automated routine for determination of parameter λ (the balance factor), (c) optional determination of values for minimum density contrast, (d) a robust handling of outlier gravity data, and (e) improved automated data reduction for terrain effects based on anticorrelation with topographic data.

In addition, a variation of the inversion approach was developed ^[20], which operates, instead of with isolated anomalous source bodies embedded in a stratified background medium (or a medium with continuously growing density with depth), with sources given as sub-horizontal layers with morphological shapes of the density discontinuity surfaces in between them determined in the growth (shaping) process. This approach can handle both the layered sedimentary basins and their layered basements. For resolving the basin shape and structure more realistically, an additional weighting matrix is introduced in the minimization term of the model, in order to force the negative contrast cells forming the basin to generate a body that is flatter, shallower, and adhering better with its upper boundary to the topographic surface. This can be done for both homogenous and stratified basin fill. Similarly, a weighting matrix was proposed for a structural crustal model with several discontinuity surfaces.

Another upgrade to the Growth approach for CBA gravity data inversion came in 2021 in the form of the GROWTH-3 software package ^[26]. Many things became automated in this version and the inversion tool became much more user-friendly. This package also includes downward density increase and layered source structures. Note that the assumed downward increase of the background medium density (either stratified or continuous) is translated into the structural model bodies or into the shaped model strata. One of the motivations for an upgrade in the package was that it can run on operating system Windows 10 Pro. GROWTH-3 is a more compact package with only two codes (programs GROWTH and VIEW). The new implemented features are: a) optional automatic determination of the prescribed values for the density contrast, b) optional sub-horizontal layered structures in the model (taken from ^[20]); c) optional increase in density with depth (taken from ^[24]), meaning that positive density contrasts increase with depth while negative contrasts decrease with depth, both either linearly or exponentially, d) new graphical presentation throughout the running iterative inversion process, and e) simpler structure of the output file for the resulting 3D model.

The followings are the values and options that the user must set (keep or override in a dialog window) prior to running the inversion process, while default values are automatically determined (based on input gravity data) and offered: a) the (mean) size (side) of the model cells (keeping the default value is recommended); b) the (positive and negative) density contrast value; c) the balance factor (λ) value; d) choosing a homogenous model or a model with downward density increase (default), with three options for the latter: (i) continuous increase, (ii) increase due to a large number of regularly spaced layers (default option), and (iii) increase due to a small number of layers as specified by the user; and e) optional regional offset and trend re-adjustment. The random search coefficient (the default value offered in the dialog box) equal to 1 means systematic yet slower model space exploration and model creation, while higher values mean more random search and a faster iterative inversion process.

A Growth inversion approach was developed as well, which is applicable in volcano-geodetic studies for simultaneous inversion of spatiotemporal (time-lapse) microgravity changes and surface deformation (elevation changes) observed at sparse scattered stations (survey benchmarks) on the topographic surface ^[21]. This inversion procedure seeks growing bodies of pressure changes and of temporal (mass) density changes, based on the aggregation of point sources (the point source of mass change and Mogi point source of pressure), that do not have to necessarily coincide. It was applied to the interpretation of gravity changes, leveling, and interferometric synthetic aperture radar (InSAR) data from the volcanic area of Campi Flegrei (Italy) for the period 1992–2000, and deformation data for the period 1993–2013 ^[23,73]. In addition, this inversion approach is applicable to studying gravity changes and surface deformation in subsidence areas due to ground water exploitation ^[33]. Applied to deformation data, it can track pressure changes associated with magma intrusions in volcanic areas ^[23,29].

A more sophisticated inversion tool based on growing source bodies, dedicated to interpretation of ground deformation data, derived from InSAR and GNSS data, has been developed ^[25]. This approach can handle a combination of pressure source bodies and dislocation sources of different types: dip-slip, strike-slip, and tensile. Such sources represent either magmatic processes or sources due to other processes such as earthquakes, landslides, or groundwater-induced subsidence. It was test applied to deformation data from the Etna volcano ^[25]. It was used to study fluid migration induced by magma injection at depth and hazards associated with flank collapses on La Palma (Canary Islands) based on SAR data and interpreted dislocation sources ^[34], as well as the onset of the unrest on the island. Volcanic unrest on La Palma culminated in the largest historical eruption on the island. The Growth inversion approach was used to study this unrest throughout 2021 using InSAR data ^[35]. The interpretation of the Growth inversion results (pressure and dislocation sources) was constrained by structural controls obtained by Growth inversion of gravity data of the island. The study revealed a shallow magma accumulation 3.5 months before the eruption.

A separate version of the Growth inversion approach was designed for 4D micro-gravimetry, dedicated to the inversion of spatiotemporal (time-lapse) gravity changes observed on the topographic surface, implemented as the GROWTH-dg tool ^[27]. The specifics of such inversion are due to the fact that the input gravity data are typically low in number, sparse, and noisy. Inverting gravity changes has raised new challenges for the Growth approach in terms of severe under-sampling of the information about the sources due to a limited (low) number of observation points. Compared to inverting complete Bouguer anomalies (CBA data) with typically good coverage and high signal-to-noise ratio, in the case of gravity changes, one deals with data of amplitudes at several tens (a few hundred at the most) of μGals with observational error at the level of (at least) 10 to 15 μGal (i.e., data usually with poor signal-to-noise ratio).

The inversion of gravity changes results in models of subsurface temporal density changes. These may represent density changes of static masses, or movement of mobile masses such as magma or its volatiles. Since geodynamic events producing gravity changes usually produce surface deformation, and the vertical displacements (elevation changes) of the surface generate their own component of the gravity change due to the vertical displacement of the observation point in the ambient earth gravity field (a coupling effect), the deformation-induced component of the gravity change must be removed from the input gravity data by means of a correction. Also, the hydrological signal representing the gravitational effects of precipitation, soil moisture changes, and ground water changes must be removed in terms of a hydrological correction. Residual gravity changes, compiled by these two corrections then serve as input gravity data in the Growth inversion.

The point gravity changes usually represent a field with very few gravity highs or lows. Therefore, the Growth solutions usually consist of a low number of isolated source bodies. The bodies can be homogenous or heterogenous. Homogenous bodies are produced in the growth process by allowing one density contrast (both positive and negative) and only one visitation and filling of each subsurface cell by this contrast. Heterogenous bodies are produced by allowing repeated (though limited) visitation and filling of the cells by a pre-specified incremental density contrast. The GROWTH-dg tool provides two ways for stipulating the final density contrast of the model source bodies. Either the density contrast (or the incremental density step and number of density levels) is set as an inversion parameter prior to running the inversion process, or the amount of filled cells (as a percentage out of all subsurface cells) is pre-set, which has impact on the final density contrast of the model when the process stops. The GROWTH-dg package allows the functionality of vertical (downward) weighting of the model cells (upward or downward "pushing" or "forcing" of negative and positive source bodies) and re-weighting of residuals serving to suppress or eliminate outliers in the data and uncorrelated signals such as site effects related to just one observation point ^[27].

The applicability of the GROWTH-dg tool to interpretation of gravity changes in volcano gravimetry was tested by Bódi et al. ^[9] by simulating synthetic sources with focus on thin elongated bodies that pose the greatest challenge in gravimetric inversion based on aggregating filled cells in a growth process.

Let us conclude the historical development overview of the Growth inversion methodology by pointing out several important givens related to all of the GROWTH implementations. The Growth solutions (i.e., the subsurface aggregations of filled cells), the models, are very sensitive to the choice of the balance factor (λ) value and the prescribed density contrast value (or, alternatively, the amount of filled cells). Several diverse models can be obtained, all with a relatively acceptable level of the data misfit in terms of r.m.s., depending on these two key inversion parameters. Although this may seem like a drawback of the methodology, it can actually mean a gain. Several "acceptable" models can be produced and then validated against existing independent geological, tectonic, and volcanological cognition, or other geophysical methods and data with their own interpretation. The user is always advised to run the Growth inversion with parameter values suggested by the software as default values (automatically deduced from the input gravity data), then sensibly vary the values of the inversion parameters in subsequent repeated inversion runs, and carefully observe the behavior of the obtained models. This will enrich the insights into the subsurface structure or the subsurface process. Illustrations of such practices are given and discussed in the following sections.

3.   Current Growth inversion implementation

The latest release of the Growth inversion tool, GROWTH-23 ^[28], merges the functionality of the two preceding tools, the Growth inversion of gravity anomalies (GROWTH-3) and the inversion of spatiotemporal gravity changes (GROWTH-dg). The novelty of this tool dwells in automatic recognition of the input gravity data (gravity anomalies versus gravity changes) and automatic provision of the default values for the inversion parameters tailored to the specific input gravity data. This new implementation integrates the specific aspects of the previous codes for gravity anomalies (CBA data) and gravity changes data. It also offers higher operational simplicity and more graphical information.

The methodology, mathematical apparatus, and its numerical implementation will not be explained here in due detail. Instead, we refer the reader to ^[28] for details. Only some key aspects of the inversion approach and the software package will be explained. The inversion of time-lapse gravity changes has its specifics. The density contrasts of the subsurface partition (see Figure 1) represent temporal density changes. The background density is not only homogenous, it is zero in the entire subsurface. The resulting model consists typically of isolated source bodies. Consequently, not many cells of the subsurface partition are filled with specific density contrasts. On the other hand, the inversion of gravity anomalies (CBA data) typically results in complex structural models, whereby nearly all the subsurface cells (see Figure 1) are filled with density contrasts. GROWTH-23 allows heterogenous models by repeated revisitation of cells and their re-filling with an incremental density contrast (both positive and negative). For inversions with homogenous background (for gravity change data implicitly, for CBA data optionally), the program adopts four levels of density contrast, while for stratified density (of the background and of the model) increasing with depth, it adopts 30 levels of density contrast. The program recognizes CBA data versus gravity change data based on the number of observations, 30 being pre-set as the threshold number. In GROWTH-23, the user can no longer specify the value of the target average density contrast and the number of density levels. This is done automatically.

Figure 1.  The design of the inverse problem based on partitioning the subsurface model space into right-rectangular (parallelepiped) cells and irregular point gravity data given on the topographic surface.

DownLoad: Full-Size Img PowerPoint

The critical inversion parameter is the balance factor (λ) which weighs the misfit residuals minimization and the total mass of the model minimization (both with L2 norm):

Above, $\boldsymbol{v}$ is the column vector of residuals of the data-fit, given as the difference between the observed gravity (gravity anomaly or gravity change, ${g}_{i}$) at the i-th data point (of n points) on the topographic surface and the sum of the gravitational effects of all the cells of the subsurface partition of the model at this i-th point (see Figure 1),

whereby the gravitational effect (vertical attraction) at the i-th point of the j-th cell (${A}_{ij}$) is given by the formula ^[66] for a right-rectangular cell of density contrast $\Delta {\rho }_{j}$. Matrix A, with components ${A}_{ij}$, is the design matrix of the physical configuration of the problem, and Eq (2) represents the observation equations of the problem. Set J represents cells of non-null density contrasts. Optionally, a term representing co-adjusted offset or linear trend enters the formulation of the residuals in Eq (2). Matrix ${\boldsymbol{Q}}_{D}$ is the a priori covariance matrix of the observables populated by the gravity data uncertainties, set as the diagonal matrix of gravity data variances.

Column vector $\boldsymbol{m}$ is formed by the model density contrast values $\Delta {\rho }_{j}$ (positive and negative) of the m cells of the subsurface partition, and ${\boldsymbol{Q}}_{M}$ is the a priori covariance matrix for the uncertainties of the model parameters. This matrix is selected as diagonal and with such elements that cause cells at different depths to have an equal probability of being included in the solution ^[48].

The above weighted mixed adjustment is employed in an iterative model exploration and growth process, in which all the subsurface cells are explored and only one of them is populated by a preset density contrast (Δρ) value (positive or negative). This process starts with a very low number of filled cells and as it progresses, the source bodies represented by aggregated filled cells grow. Therefore, Eq (1) must reflect this growth and adapt to the number of filled cells, aka the strength of the gravitational effect of the model, throughout the growth. To adapt to the growth, it adopts a scaling factor (f):

with

where positive and negative density values $\Delta {\rho }_{j}$ in (3) and (4) have been expressed as ± $f$ and where J+ and J- are the set of indices corresponding to those cells filled with positive and negative density values. The term ${g}_{tr}$ is the optional co-adjusted offset or linear regional trend. The iterative process progresses from a small number of cells filled with high density contrasts to larger amounts of filled cells and smaller values of $f$ and $e$. The iterative growth process ends when the scale factor reaches $f$ = Δρ, or conversely when it is not possible to find a new cell with values $f$ and $e$ smaller than those of the preceding step. In this last case, the last $f$ value would be the resulting model density contrast. Heterogenous Growth models are obtained by allowing the iterative process to visit and fill any cell more than once (yet limited to a prespecified number), with the preset incremental density step Δρ. This was so until the arrival of the GROWTH-dg and GROWTH-23 implementations of the model exploration and cell filling of the growth process. The GROWTH-dg tool optionally and GROWTH-23 implicitly use a different end condition. Instead of pre-specifying the Δρ value at which the scaling factor shall arrive and the process complete, the alternative termination of the growth process is based on pre-specifying the amount of filled cells (in %). The process stops when the prespecified % of filled cells is reached, resulting in the final density contrast of the model (in case of homogenous models) or the final average density contrast (in case of heterogenous models).

The balance factor (λ) is the factor for the selected balance between tight data fit and compactness (smoothness) of the model. Low values produce good fit (low misfit residuals) but overly complex (crowded, messy) models with large total anomalous mass (both positive and negative). High λ values produce simple compact models with low total mass and rounded source bodies but poorer data fit (large residuals). We suggest to pick the optimal λ value through an autocorrelation analysis of the corresponding residuals, obtained as a function of the mutual distance between the stations ^[26,27]. The optimal value would correspond to the null autocorrelation. Default values are set as λ = 50 for CBA data and λ = 30 for spatiotemporal gravity changes. Nevertheless, for different application cases, the user must aim for a better specific λ value by observing the autocorrelation distribution to aim at a null value for zero distance. In addition, the user should observe the misfit r.m.s. and the map of final residuals, as well as the behavior of the adjusted model (aiming at not too noisy, not too compact), in response to changing the value of the balance factor in a trial-and-error process.

Small values of λ produce overfitted models whereby the data noise is translated into model noise, whereas high values of λ produce over-regularized (over-smoothed, over-compacted) models with less source bodies of smaller sizes and rounder shapes also vertically incorrectly displaced. The default values of all inversion parameters can be changed by the user manually prior to running the inversion procedure.

In case of the presence of outliers (blunders), the data reweighting functionality should be enabled, suppressing or eliminating the extreme residuals (outlier data) contribution into the adjustment. This iterative reweighting is based on successive assigning of smaller weights to large residuals throughout the growth process. This reweighting is controlled by the "blunder value" parameter (B), which multiplied by the standard deviation of the gravity data defines the domain of the prevailing normal Gaussian distribution, beyond which the values of residuals are considered outliers to be suppressed or eliminated. In each iteration step of the growth process, the median of the absolute values of the residuals (divided by 0.6745) is used as a robust estimator of the standard deviation. Standardized residuals (i.e., residuals divided by this estimator) are used to compute the weights in each iteration step (see Eq (19) in ^[17]). The scale factor f is used in this equation, too, for eliminating premature outlier exclusion in the early stages of the iterative growth process when the small amount of filled cells cannot guarantee a good enough data fit, yet, tolerating the occurrence of large residuals at early iteration stages. For more details and the mathematics of the iterative residuals reweighting functionality, see Section 5.3 in ^[17]. The cut-off values of the B parameter are set as 3.8 for CBA data and 2.2 for gravity changes. Values of B higher than cut-off values mean weak or no reweighting.

GROWTH-23 uses the pre-set amount of filled cells (in %) out of all subsurface cells as the end condition for the iterative inversion process. Suggested default values are 60 % for CBA data and 3% for gravity changes. Again, these values can be changed by the user prior to running the inversion. This parameter is also a key parameter in the Growth process, as it implies the values of the final model density contrasts. Therefore, the user should sensibly vary this value in a trial-and-error approach (along with sensibly varying the λ value), while observing the character of the resulting model. A small size model (with a small relative number of filled cells) will require larger density contrasts, while a very large model (with a high percentage of filled cells) will require smaller density contrasts. By experimenting with the end value during repeated inversion runs, the user can target the final model to a desired or expected average density contrast solution.

The choice of using a stratified background medium (as well as model) with step-wise linear downward increasing density in structural studies based on CBA data is mediated by the "flattening coefficient". GROWTH-23 includes an option for automatic stepwise linear increase of density of model strata by means of this flattening factor. This factor forces the growth process to alter from growing the source bodies to shaping the strata. Shaping the strata means growing their morphologies by step-by-step growth of the undulations of the density interfaces among the strata. With this feature, the morphology of the resulting models changes from "muffin-type" or "onion-type" to a "bell-shape" pattern. For more details regarding the growth process in stratified background medium, see Section 3.3.7 in ^[28]. The flattening coefficient can attain values from 0 to 100. The higher the flattening factor value, the stronger the flattening. A pre-set zero value of the flattening factor means no stratification (homogenous background). This option is proper for inverting gravity changes. The proposed default value for CBA inversion is 10. Again, trial-and-error repeated inversions with sensibly varying the flattening factor are recommended to get a feeling for the model behavior.

An inherent feature of the inversion of gravity data is the loss of sensitivity with depth, as well as to the periphery of the model domain as related to the horizontal extent of the input gravity data. Another one is the inability of resolving vertically overlapping positive and negative contrasting source bodies. Growth inversion of CBA gravity data generated by shallow layered sedimentary basins tends to produce too deep model solutions. This tendency can be suppressed (counter-acted) by (optionally) including additional depth-weighting in the model covariance matrix (${\boldsymbol{Q}}_{M}$), mediated by the so called "coefficient for upward forcing". The application of this coefficient forces the strata (if applying stratification with flattening) or the isolated source bodies (if applying no stratification—i.e., homogenous background) upwards to shallower depths. Not applying this depth-weighting factor is referred to as (inverse) "modelling with neutral depth". For more details on the upward forcing functionality, see Section 3.3.8 in ^[28]. For inversion of gravity changes the proposed default value for the upward forcing parameter is zero. For inversion of CBA data in complex structural settings and larger areas the proposed default value is pre-set at 50. However, this value should be used with caution or varied sensibly in a trial-and-error process of repeated inversion runs. External constrains must be used wherever available, such as borehole data. If borehole data are available, the stratified Growth model can be adjusted, in terms of the depths of the individual density interfaces among the strata, to the borehole data. The flattening coefficient and the upward forcing coefficient both modify the model covariance matrix (${\boldsymbol{Q}}_{M}$) of the mixed adjustment at the beginning of the iterative growth process. This covariance matrix does not change throughout the iterative growth process.

When the iterations stop, the growth process produces a final model. The 3D model is graphically displayed in terms of several horizontal slices (at automatically selected depths) and several W–E vertical cross-sections (at automatically selected northings). Then, two ASCI output files are produced. File "Mod.txt" contains the values of the inversion parameters applied in the run and information of the resulting Growth model described as a list of filled cells. For each cell, the file contains: UTM coordinates (easting, northing, elevation above sea level) of each cell center, cell sides sx, sy, sz, cell density contrast (kg/m³), and relative sensitivity of the cell. File "Fil.txt" contains a list of observed, modelled, and residual values, including offset/trend. File Mod.txt can be used to input and display the 3D model in third-party software.

Regarding the uncertainty ("precision") of the Growth solution (model), it is possible to compute the posteriori covariance matrix for the adjusted model parameters, at completion of the iterative growth process, including the uncertainties of the density contrast values in each populated cell of the Growth model (cf. Eq (22) in ^[17]). The question arises as how useful would such a quantification of model uncertainty be in the light of the fact that the model density contrast is not a unique product of the Growth inversion, rather a mediator and a free inversion parameter of the Growth inversion. The model density contrast is selected by the user, via selecting the amount of filled cells (%), prior to running the iterative growth process. Several diverse Growth models with diverse density contrasts may be nearly equally acceptable based on their data misfit. Then what value would such an estimated model uncertainty offer? More than with precision, we are concerned with the accuracy of the Growth solution. Therefore, it is more telling to compare the obtained Growth models to simulated synthetic models or known solutions, in order to assess the ability of Growth inversion to reproduce the original source, given the specifics of a particular source and the settings of the environment in which it is located. If seeking an unknown source, all possible available constraining information (geological, tectonic, volcanological, or stemming from other geophysical methods) must be employed in order to reduce the great ambiguity of purely gravimetric solutions.

To conclude the description of the Growth inversion approach, we illustrate the growth process workflow in Figure 2.

Figure 2.  Chart of the Growth inversion workflow.

DownLoad: Full-Size Img PowerPoint

4.   Applications across earth science disciplines

The applicability of the Growth inversion methodology ranges from inversion and interpretation of regional or local gravity, given as complete Bouguer anomalies (CBA data), to assist the determination of the structure of the earth crust, or of sedimentary basins and their basements, through near surface applications for gravimetric detection of cave space in karst areas, void spaces in undermined areas of active or abandoned mines posing threats to infrastructure, cavities such as tombs, crypts, or other buried anthropogenic objects in archeological prospection, to volcano gravimetry for studying volcanic unrest or awakening volcanoes based on observed spatiotemporal gravity changes. We illustrate this applicability here on several selected case studies.

4.1.   Structural geophysics applications–gravity anomalies

The Growth inversion method was applied to obtain gravimetric structural density models of several volcanic islands and sedimentary basins. The GROWTH inversion tool for gravity anomaly inversion was first applied to CBA data on Grand Canaria, Canary Islands ^[14], and to CBA data on Lanzarote ^[15] to analyze the crustal structure of the islands based on the inverted density distribution. Likewise, yet complemented by additional geophysical methods, it was used to determine the structure of Tenerife, Canary Islands ^[2]. A structural gravimetric study based on GROWTH inversion of gravity data was conducted also on Terceira Island ^[55] and Pico Island, in the Azores ^[57]. A 3D structural density model for the volcanoes Merapi and Merbabu in Central Java, Indonesia was constructed based on Growth inversion ^[77]. The intrusive complex under Faial Island of the Azores archipelago was gravimetrically determined using Growth inversion ^[17]. In this application, an enhancement to the method was adopted: the terrain density for topographic correction for compilation of CBA data was simultaneously estimated in the adjustment. The analysis of the structural composition of the Vesuvius volcano and its neighborhood was conducted by applying the Growth inversion of CBA data using the "shaping layers" modification of the GROWTH tool ^[5]. The CBA data set, improved by new data in the central area of Tenerife, was inverted using GROWTH ^[40] to image the high-density core and the remaining structural composition of the island with enhanced resolution and detail. Structural gravimetric study using the Growth approach was carried out also for La Palma of the Canary Islands ^[18].

The enhanced Growth approach implemented as the GROWTH-2 software package was applied to a new, improved CBA dataset on Tenerife ^[19] yielding improved information on the location and shape of the main structural building blocks of the island. A structurally layered model was obtained for the same data, too ^[20]. The shape and stratification of the alluvial basin in the Low Andarax valley (Almeria, Spain) was determined gravimetrically using the Growth inversion approach ^[22]. An improved gravimetric structural model of Lanzarote (Canary Islands, Spain) was obtained ^[24] by inverting an enhanced CBA gravity data set using the GROWTH-2 package ^[19,20]. Program GROWTH-3 was applied to invert the CBA gravity data of El Hierro (Canary Islands) to obtain a 3D density model of the crustal structure of the island ^[26].

4.2.   Volcano geodesy applications–spatiotemporal gravity changes

Due to the limited number of input gravity data, the volcano-gravimetric studies usually opt for inversions based on assumed simple geometric source bodies, their gravitational effect being given by an analytical formula. In such inversions the number of sought source parameters is smaller than the number of input gravity data. The source parameters are then found by optimization methods such as the genetic algorithm. This approach might better suit cases where there is external evidence or justified reasoning for choosing the number of sought sources and their shapes (e.g., spheres, ellipsoids, cylinders, prisms, etc.). On the other hand, inversion methods which do not pre-constrain the number or shapes of sources are better suited in the absence of such externally justified assumptions and may identify the subsurface locations of the most significant density changes, indicate the nature of the sources, and shed light on the subsurface geodynamic process. We compared the two approaches in a case study at the Laguna del Maule volcanic field, Chile ^[83]. Spatiotemporal gravity changes observed on Ischia (Italy), related to the 2017 earthquake, were inverted and interpreted using the Growth approach ^[6]. Also, the 3D displacements of the 2014 Napa Valley, California (USA) earthquake ^[25] were inverted and interpreted using Growth. We have also revisited the gravimetric interpretation of the volcanic unrest of 2004–2005 on Tenerife (Canary Islands, Spain) by applying the Growth inversion to the observed time-lapse gravity changes ^[84].

Here we show the re-inversion of the spatiotemporal gravity changes of the Tenerife unrest using the GROWTH-23 inversion tool. The observed and processed data were provided by ^[39]. Here we aim at illustrating the sensitivity of the Growth model to the opted values of inversion parameters, namely to the balance factor (λ), to the relative amount of filled cells affecting the final average density contrast of the model, and to the strength of residuals reweighting. In Figure 3, we show three Growth solutions obtained with two different values of the balance factor (λ = 20 and 30), two different end conditions (% = 0.2 and 0.4), as well as without and with fairly strong reweighting of residuals (B = 1.2). In the figure, the first row is the top view, the second row is the lateral view from the east, and the third row is an oblique view at the 3D Growth model. The last row lists the values of the applied inversion parameters, the average density contrast of the model, the misfit residuals (r.m.s. and max value) and the amount of injected (intruded) mass. The contour lines represent topographic relief. Green dots are the gravity points of the gravimetric monitoring network. Filled cells with positive density contrast are shown in red.

Figure 3.  Three different Growth solutions due to the diversely selected values of the inversion parameters. For the description, see the text.

DownLoad: Full-Size Img PowerPoint

We do not interpret these three models here. They only serve the purpose of demonstrating the sensitivity of the Growth model to the selected values of the inversion parameters. None of the three presented models represents the final model that served the gravimetric reinterpretation of the unrest. For the choice of the most acceptable (optimum) model, based on external criteria or cognition, and its interpretation, see ^[84]. In any case, even these diverse solutions hint at a stalled magma intrusion and failed eruption. They are useful for detecting possible pathways for the ascent of magma and its volatiles. The filled cells (shown in red) of positive density contrast are not to be understood as injected magma replacing the edifice rock or filling a previously void space defined by the red cells. They represent a volumetric domain (a source body in terms of Growth inversion) of a bulk density change—i.e., a subsurface domain affected, for instance, by a swarm of thin magmatic dykes and sills, or a system of cracks and void spaces filled with magma and/or its rising volatiles (brines).

In the solution presented in the left column, the northern-most source body is considered an artifact, as it is at the border of the area covered by the data and lies in a domain of very low sensitivity of the solution. The shallow SW source body is considered to represent risen volatiles at the outskirts of the caldera, at the depth of the aquifer.

4.3.   Near surface applications–microgravity anomalies in archeological prospection

Unlike in 4D microgravimetry, where the inversion approach faces the challenge of sparse microgravity data low in number, scattered over often an extensive area and with low signal-to-noise ratio, in 3D microgravimetric applications, the gravity data are given typically on a small area, on a dense mesh, and have good signal-to-noise ration. To test how well the Growth inversion performs in near surface microgravimetric studies, we inverted a gravity data set collected for the identification of shallow tombs and crypts in archeological prospection of the Basilica of Saint Nicholas, Trnava, SW Slovakia.

The basilica is relatively large, about 60x30 m² and is known to hide several crypts. Their exact locations were then unknown, except for one. A joint microgravity and ground-penetrating radar (GPR) survey was conducted ^[62,76]. Most of the crypts were detected by the microgravity method and confirmed by GPR. One crypt was not identified by gravimetry, but later was discovered by GPR. Subsequently, the findings were verified by drilling and video inspection, confirming the presence of six new crypts ^[63,65].

The gravity data were collected in a dense regular grid with step of 1 m and corrected for the gravitational effect of the walls of the basilica building. Seven significant anomalies (gravity lows) were identified in the map of the local CBA data, see Figure 4 (left), with amplitudes up to 50 μGal. These anomalies were initially attributed to cavities, which was subsequently confirmed by means of 3D density modelling and Euler deconvolution ^[63]. The GPR survey confirmed all seven crypts and revealed one additional crypt, left undetected by the gravity survey, due to the masking effect of a very close wall.

Figure 4.  A schematic plan of the Basilica of Saint Nicholas with the overlaying residual CBA field (left). Results of the GPR survey with crypt positions indicated (right–reproduced with permission from ^[76]). Inset indicates the geographical location of the study site.

DownLoad: Full-Size Img PowerPoint

In order to understand the applicability and suitability of the Growth inversion in near surface applications such as archeological microgravimetric prospection, we inverted the same gravity data set as used in the previous studies, quoted above, and compared the Growth achievements to the previous results. The Growth inversion was carried out with the GROWTH-23 package. The cell size was selected at 1.1 m. Co-adjustment of offset or trend was disabled. Reweighting of residuals was disabled, as the input gravity data were dense, well correlated, and fairly accurate. No flattening (stratification) and no depth weighting (upward forcing) was applied, since we were exploring the immediate subsurface within several meters. The two key parameters affecting the appearance of the Growth model were the balance factor (λ) and the end condition in terms of the amount (in %) of filled cells. We were experimenting with these two values to find the most suitable model for the identification of the crypts ^[8].

We started with selecting the λ value first, then terminated the model at lower and subsequently higher percentages of filled cells. Models with higher λ allowed us to focus on the most significant sources (crypts) with suppressed model noise, while lower λ values enabled us to explore less influential anomalies, which can also lead to revealing less distinct source objects, unless they are representing artifacts. Rising the end value allows the source bodies to inflate and more sources of small influence to pop up. This is a good example of selecting the inversion parameters based on the situational knowledge of the site. Increasing both λ and the end condition (%) produces more voluminous and less dense anomalies. However, we are seeking cavities (unless filled with debris), so high density contrast in the solutions is preferable.

In Figure 5, we present the Growth solution that best suits the gravimetric detection of the near-surface cavities identified as crypts in the basilica (λ = 150, % = 10, r.m.s. residuals = 5 μGal, average depth = 4 m below the ground floor level, average density contrast = -493.7 kg/m³). For the sensitivity of the Growth solutions to the values of the individual inversion parameters of the GROWTH-23 package, as well as the experience learned in this case study, we refer to ^[8].

Figure 5.  Growth model (λ = 150, % = 10) of immediate subsurface of the Saint Nicholas Basilica. Section (1): Planar (top) view at the 3D model including an overlay of the basilica scheme. Section (2): Vertical cross-sections along the profiles depicted in section (1). Section (3): Horizontal slices at three selected depths.

DownLoad: Full-Size Img PowerPoint

For relations between Figures 4 and 5, negative density (blue) anomalies corresponding to known crypts are marked with red capital letters A–M. In Figure 5, it can be seen that a significant number of crypts were successfully recovered by the inversion, most dominantly those manifested by the strongest gravity low in the residual CBA map (see Figure 4 (left)). The selection of a slightly higher λ in combination with an end value of 10% yielded a reasonably balanced solution in terms of data fitting, with misfit (r.m.s.) of 5 μGal. Horizontal slices of Section (3) reveal the positions of the individual crypts. The volume of the recovered source bodies is reproduced ambiguously, depending on the selection of inversion parameters. The solution includes significant positive density contrast sources around the building's perimeter walls, which could be attributed to the insufficient correction of the gravity data for the gravitational effect of the outer walls and particularly of their underground foundations that remain unknown. The presence of an unknown source marked by "?" below the presbytery was detected throughout multiple Growth solutions, although the previous gravimetric study ^[63] did not suggest such an object, while it was not possible to conduct GPR observation at this spot. Therefore, its true origin remains a mystery for now.

Some crypts, such as objects "A" or "K" in this solution, are obtained quite ambiguously, as it is difficult to distinguish their actual presence from model noise. In ^[8], more solutions are presented, where the appearance of crypt "A" is more evident. This illustrates that a suite of solutions with sensibly varied inversion parameters is more telling than just one solution.

All in all, the experience gained during this case study is that the Growth inversion approach and the GROWTH-23 package are suitable for microgravimetric archeological prospection applications.

4.4.   Near surface applications–microgravimetric prospection for hidden cavities in mining areas

Additionally, we tested the applicability and suitability of the Growth inversion to prospection for hidden cavities in mining areas that may represent sinkhole or surface collapse hazards. We inverted a microgravity data set acquired above an abandoned coal mine in Wolfsberg, Austria. The studied area lies above a defunct coal mine complex, which was abandoned in the 1950s. The mined coal deposit is located at depths of 15–25 m below the surface (e.g., ^[89]). The mining was carried out by the "chamber and pillar" method. Threat is associated with such mining methods, as the collapse of a chamber can cause sinkholes on the surface. This area was micro-gravimetrically surveyed. We used the residual (trend-removed) CBA dataset of ^[89] on a dense regular mesh with step of 2 m with estimated accuracy of 12 μGal. The CBA map exhibited a well-pronounced gravity low with an amplitude of -35 μGal (Figure 6).

Figure 6.  The residual CBA gravity map of the study area. Inset shows the geographical location of the site. The black line indicates the position of the data profile used for 2D modelling in ^[89].

DownLoad: Full-Size Img PowerPoint

To estimate the depth of the source of this anomaly, Zahorec et al. ^[89] used 3D Euler deconvolution ^[70] with regularized derivatives ^[64], resulting in a cluster of depth solutions with mean depth at approximately 3.7 m below the surface. This solution was confirmed by forward gravity modelling obtaining the cavity at an average depth of 4.65 m (spanning from 2.3 to 7 m) using a density of -2200 kg/m³. These solutions were later verified by a borehole exploration campaign of a private company, encountering a substantial cavity with the upper edge at the depth of 3.5 m, with a height of 3 to 4 m, and a diameter of roughly 5 m ^[89].

We inverted the same residual CBA dataset using the Growth approach with the GROWTH-dg software package, seeking a homogenous source body representing a cavity with the density contrast of 2200 kg/m³. The average cell size was set at 1.3 m. We set no offset or trend co-adjustment, no reweighting of residuals (B = 8), and no depth weighting (D = 0). We varied the values of the balance factor in order to find the best-recovery solution, by comparing the obtained Growth solutions with the result quoted above. The best-recovery Growth solution (see Figure 7) was obtained with the following parameters: λ = 75, Δρ = -2200 kg/m³. This is a tight-fit solution with data misfit: r.m.s. = 8 μGal, mean = 1 μGal, max = 20 μGal. Although we ran the inversion on data over the whole data area (cf. Figure 6), in Figure 7, we zoom on the cavity as imaged by the Growth source body.

Figure 7.  The best-recovery homogenous Growth model (density contrast of -2200 kg/m³) of the cavity with a 3D outline (shown as a cylinder) of the cavity verified by the borehole and 2D modelled by Zahorec et al. ^[89]. The two insets on the right show top and lateral (from south) views at the 3D model. Contour lines represent the topographic surface, and the green dots are the gravity data points.

DownLoad: Full-Size Img PowerPoint

The solution in Figure 7 consists of a single source with negative density contrast of the air relative to the rock environment, representing the sought cavity. The position and horizontal size of the cavity was recovered correctly. As for the depth span of the cavity, the Growth solution senses only the upper part of the cavity, just like the Euler deconvolution solutions of ^[89]. Our Growth-reproduced cavity starts at 2.5 m below the surface and reaches 4.5 m below the surface. Although the Growth inversion has not revealed the whole vertical size of the cavity, it caught its upper part— i.e., how close to the surface the cavity developed—which is an essential piece of information for sinkhole hazard assessment.

5.   Discussion and conclusions

The previous and current packages for inversion of gravity anomaly data (CBA data) based on the Growth approach have proven to be successful in structural geophysical studies, as reviewed in this contribution. Likewise, the GROWTH-dg and GROWTH-23 inversion tools have demonstrated their strengths and benefits in inverting spatiotemporal gravity changes, producing subsurface time-lapse density change models that can be successfully interpreted, as reviewed in this work. In addition to these structural and geodynamic applications of the Growth methodology, we tested its applicability in near-surface gravimetric studies. We chose two 3D-microgravimetry case studies: one from the realm of archeological prospection for crypts and tombs, and the other from the field of undermined area hazard detection and evaluation by uncovering cavities that potentially can lead to sinkholes. Both test case studies demonstrated the Growth solutions to be a success. We can conclude that the Growth inversion approach, based on the GROWTH-23 implementation, is applicable across a wide range of earth science disciplines, dealing with a wide range of gravity data, from dense accurate 3D microgravity data on a small area, through sparse inaccurate 4D microgravity data over fairly extensive areas, to regional CBA gravity data. The Growth inversion methodology has manifested, through these individual studies, that it can yield benefits over other standardly used methods.

Use of AI tools declaration

The authors declare that they have not used artificial intelligence (AI) tools in the creation of this article.

Software availability

The GROWTH-23 source codes, user manual, and executable files can be freely downloaded from the repository at https://github.com/josefern/GROWTH-23 under an open-source license. The GROWTH-23 software is written in FORTRAN and is compatible with the Microsoft operating system Windows 10 Pro. It was compiled using the Intel® Visual Fortran Compiler with a Quickwin application project. The complete compilation requires several additional files (codes, auxiliary files for compilation, etc.), which are also freely available on the website.

Acknowledgments

This work was partially supported by the Slovak Research and Development Agency under the contract (project) No. APVV-19-0150 (acronym ALCABA), by the VEGA grant agency under project Nos. 2/0002/23 and 1/0587/24, and under the ERA.net program ERA.MIN-2 project "Deposit-to-Regional Scale Exploration" (acronym D-Rex). A.G.C. and J.F. were supported by a Spanish Agencia Estatal de Investigación (10.13039/501100011033) G2HOTSPOTS (PID2021-122142OB-I00) grant.

Conflict of interest

The authors declare that they have no known conflicts of interest, competing financial interests, or personal relationships that could have appeared to influence the work reported in this paper.