3-D crustal density modeling, Case study: Iran

Various sources of geodetic, geological and geophysical information can improve our understanding about the crustal density. In this paper we propose a methodology for determination of crustal density based on terrestrial gravity observations, geopotential models, Digital Terrain Models (DTM), geological maps and cross-sections, deep geophysical profiles, and geoid models. The method can be algorithmically explained as follows: (i) Extraction of density information of topography (outer part of the crust above the geoid) from the existing geological maps, cross-sections and knowledge about geological formations. (ii) Removal of the global gravity field from terrestrial gravity observations via a geopotential model plus the centrifugal effect. (iii) Removal of the effect of the terrain masses above the geoid using Newton's integral and the computed density model in step (i). (iv) Downward continuation of the resulted residual observations from the surface of the Earth down to the geoid by free-air reduction. (v) Interpolation of the downward continued gravity residuals to develop a regular grid on the geoid. (vi) Setup of the observation equations based on Newton integral for the gridded residual observations. (vii) Solution of the Newton integral equations using the least squares method, stabilized by regularization techniques. As practical capability of the method its patch-wise implementation can be mentioned, which allows changes of resolution and location of the derived density model. This capability enables zooming and applying the method over different geological features such as faults, subduction zones, and other tectonic features.
The efficiency of the mentioned methodology was first approved by a simulation and next was used for the computation of a new residual density model of the crust between the geoid and the Moho discontinuity in the geographical region of Iran. Following is a summary of the practical implementation of the method. From the study of the published regional deep sounding geophysical profiles we came to the conclusion that the case study region is consists of 3 layers and accordingly 3 matrices consist of the depths of the layers with 1 degree resolution was compiled. This part is critical part in our computations since it has been proved that any error in the location of the layers can adversely affect the computed density model. Having derived 42.5 km as the average Moho depth of the region 151 was selected as the maximum related degree and order of spherical harmonic expansion for the removal of the long wavelength signals from the gravity observations. As the global geopotential model to supply the mentioned spherical harmonic coefficients, EGM2008 up to degree and order 151 was implemented. The needed point gravity data of the study were supplied from two sources, namely 6800 point gravity from BGI (Bureau Gravimetrique International) and 6500 point gravity from NCC (National Cartographic Center). Newton integral over topographical masses was used to remove the short wavelength gravitation signals from the gravity data. The boundary of the integration over the topographical masses was provided by SRTM digital terrain model with 30 arc second resolution, and the needed mass densities were extracted from 1/250000 and 1/100000 geological maps of Iran, published papers, and National Iranian Oil Company (NIOC) reports. The regional geoidal heights needed for downward continuation model was taken from the latest geoid model developed by the Surveying and Geomatics department of the University of Tehran. The resulted gravitational effects were then downward continued to the surface of the geoid by free-air reduction.
After application of the band-pass filter to the surface gravity data the resulting residuals were downward continued to the surface of the geoid. Next, the residual gravitation values were introduced to left-hand side of the Newton integral equation to start our constrained inversion process. The final product of the computations, i.e. the density variation models of the 3 crustal layers, provided us with a new 3-D density model of the case study region. According to the results, the range of density variations within the first layer is -120 to 40 (kg/m^3), while within the second, and third layers the range is almost the same, i.e. -40 to 40 (kg/m^3), which is due to the natural characteristic of the density layers, i.e. the deeper the layer the smoother the density variations. Besides, the calculated density models show remarkable correlation with (i) the observed residual gravitations, (ii) main tectonic units of Iran i.e. Zagros, Alborz and at margin of Lut and Dashte Kavir , Kopetdag, and (iii) Moho depth. Those correlations are also depending on the amplitude and wavelength of source crustal anomalies. The minima of both residual gravitation and residual density model are over Zagros, Alborz and at margins of Lut and Dashte Kavir where the maxima of the both quantities are located over Lut, Ddashte Kavir and the west of Tehran province.