Grav3CH_inv: A GUI-based MATLAB code for estimating the 3-D basement depth structure of sedimentary basins with vertical and horizontal density variation

Highlights

• The code estimates the 3D depth structure of a sedimentary basin from gravity data.

• The routine combines advantages of both the wavenumber and space domain techniques.

• It allows consideration of density variations in vertical and horizontal directions.

• Synthetic and real data applications demonstrate the efficiency of the algorithm.

Abstract

This study presents a software tool Grav3CH_inv developed to estimate the three-dimensional depth structure of sedimentary basins from their gravity data through an iterative process. The algorithm linked to the developed code operates recursively both in the wavenumber domain and in the space domain based on a triple method combination. The modelling strategy allows considerations of the model space with an exponential increase in density with depth and also with density variations in the horizontal direction. The accuracy of computation of gravity anomalies in the FFT-based forward procedure is also increased by using the shift-sampling technique which minimizes discretization effects during transformations between the space and wavenumber domains. Given the observed gravity anomalies and the density design of the basin, the iterative procedure performs automatically until the goodness of the fit between the observed and modeled anomaly, either the root-mean-square error or the largest error is below its pre-assigned value. As an advantage, the computing time is acceptably short for such a kind of an modelling problem. The GUI-enabled interactive control functions of the Grav3CH_inv code allow users to set optional settings, style of outputs and export formats, and facilitate operations without requiring coding expertise to perform the relevant procedures of the algorithm. The feasibility and accuracy of the proposed software is demonstrated by evaluating the synthetically produced gravity anomalies of various 3D basin models and also by analyzing an actual gravity data from the Los Angeles basin, California.

Introduction

Investigations on sedimentary basins are generally associated with their development and economic significance. Economically, sedimentary basins are the most important sources of energy-related products such as oil, gas, coal, uranium and geothermal fluids, and the largest deposit reservoirs for many minerals. Determining the presence and spatial position of smooth depressions and uplifts in the basement relief of a sedimentary basin may be guiding in locating stratigraphic and structural oil traps (Chakravarthi et al., 2007; Silva et al., 2010). Besides, basement relief estimates can also provide important contributions in understanding aquifer structures in hydrogeological studies (Bohidar et al., 2001; Himi et al., 2017; Lekula et al., 2018), and in understanding the flow rate of the discharge in glaciological investigations (Venteris and Miller, 1993), or can have important implication as a density determination tool in landslide studies (Mantlík et al., 2009).

The gravity method is an effective way to approximate the depth structure of a sedimentary basin because variations in gravity can be observable due to the presence of substantial contrast in density between the sedimentary infill and the mass below the basement interface (Silva et al., 2006, Chakravarthi and Ramamma, 2015). The commonly used mathematical geometries in sedimentary basin modelling are the polygonal model (Talwani et al., 1959) or the stacked prism model (Bott, 1960). A number of algorithms available apply a uniform density in their inversion schemes for the units above the basement interface (Cordell and Henderson, 1968; Murthy and Rao, 1989; Barbosa et al., 1997, 1999; Mendonca, 2004; Pallero et al., 2015, 2017; Ekinci et al., 2021). However, for the thick basins, the density of sedimentary rocks varies with depth, so this assumption is often unrealistic (Chakravarthi et al., 2013, 2017; Mallesh et al., 2019; Florio, 2020). Several studies have shown that the densities of the sedimentary fills generally increase with depths due to increasing pressure (Athy, 1930; Hedberg, 1936;Hughes and Cooke, 1953; Cordell, 1973; Crosby et al., 2006; Tenzer and Gladkikh, 2014; Cai and Zhdanov, 2015). Therefore, the practical benefits of algorithms considering constant density for the interpretation would be insufficient, particularly when the densities in the mass above the interface are non-uniform with depth. Hence, other methods that incorporate varying density models with depth for more accurate modellings have been proposed, for example, the linear density function (Murthy and Rao, 1979), the quadratic density depth relation (Bhaskara Rao, 1986; Feng et al., 2016), the hyperbolic law (Litinsky, 1989), the cubic function (Garcia-Abdeslem, 2005), the parabolic density depth relation (Chakravarthi et al., 2002; Chakravarthi and Sundararajan, 2004; Silva and Santos, 2017) and the exponential density function (Granser, 1987).

In general, the sediment density increases rapidly at shallow depths and progressively less rapidly towards deeper. Although each density function has its own benefit for a reliable interpretation, density information from wells, lithological records, and sample density measurements are always needed to define an appropriate function that could give the most effective solution. The density-depth dependency of a sedimentary basin may often be simulated by an exponential relation for the case the high contrast in density exists at the shallower depths, while the density contrast attains its minimum at the lower stratigraphic unit of a thick sedimentary sequence (Athy, 1930; Cordell, 1973; Maxant, 1980; Garcia-Abdeslem, 1992; Chappell and Kusznir, 2008; Chakravarthi et al., 2013; Gu et al., 2014; Mallesh et al., 2019). If the appropriate function is exponential, it is not possible to derive an analytical expression in the space domain for the gravity effect of anomalous sources (Chai and Hinze 1988; Chakravarthi, 2009). To achieve this, a number of frequency domain algorithms have been developed to simulate the gravity effect of source bodies varying in density exponentially with depth. For instance, Cordell (1973) divided the model space of a profile section into prismatic compartments and presented an iterative routine which incorporated both the gravity effect and its vertical gradient. Granser (1987) developed a frequency-domain algorithm based on Parker's (1973)method for calculations of the model gravity anomalies of structures including exponential change in density with depth. Later, Chai and Hinze (1988) derived a forward formula for a 3D vertical rectangular prism realized in the wavenumber domain in which an exponential change in density contrast was considered, followed by a shift-sampling technique to improve the precision of transformation from the wave number domain to the space domain. In another attempt by Bhaskara Rao et al. (1993), graphical methods deduced from characteristic values of frequency-domain transforms derived for a number of simple shaped bodies with exponential density contrast were used for the analysis of the gravity anomalies of sedimentary basins. Feng et al. (2014) extended the Parker-Oldenburg algorithm for a 3D case and introduced a gravity data inversion compatible with lateral and vertical density consideration.

It is evident that frequency domain techniques in calculating gravity anomalies are faster than space domain techniques, thus, FFT-based methods result in computer time-efficient practical applications (Granser, 1987). On the other hand, the utility of many methods depends on the precision of the transformations between the space and wavenumber domains (Chai and Hinze, 1988). Several researchers have built computer programs based on FFT algorithms, including inversion schemes to reconstruct the density interface geometry from the gravity data. Although those based on Parker (1973) and Oldenburg (1974) have a short estimation time and are capable for large data sets (Nagendra et al., 1996; Gomez–Ortiz and Agarwal, 2005; Shin et al., 2006), a fixed datum plane and/or filtering of the data is required for their procedural convergence of inversion. Another computer program for the gravity modelling of sedimentary basins by Pham et al. (2018) performs a recursive procedure by computing gravity anomalies using the forward algorithm of Granser (1987) and modifying the depth estimates in the spatial domain using the modelling scheme of Cordell and Henderson (1968). Although their algorithm by this combination does not require filtering and the average depth plane of the interface, the density-depth dependence of the model space is limited to the assumption that it varies only in the vertical direction.

In this paper, it is aimed at presenting an alternative software tool for simulating the 3D depth structure of a sedimentary basin from its gravity field. The proposed algorithm is based on utilizing the inverted Bouguer slab relation of Cordell (1973) in a recursive procedure which combines both the frequency domain method of Chai and Hinze (1988) and the space domain technique of Cordell and Henderson (1968). By this combination, the modelling scheme is capable to allow considerations of the density variety of the model space both in the vertical and horizontal directions simultaneously. The precision of the transformations between the space and wavenumber domains during computation of gravity anomalies is also improved by adopting the shift-sampling method suggested by Chai and Hinze (1988). The presented algorithm linked to the code Grav3CH_inv.m is built in MATLAB (R2013b) with a graphical interface which guides the user to manage the operation simply without any requirement of coding experience. The software is tested for its practical application and accuracy on both synthetic and actual gravity data.