Inversion of the density structure of the lithosphere in the North China Craton from GOCE satellite gravity gradient data

Lithospheric thinning in the North China Craton (NCC), one of the oldest cratons on Earth, has been a focus of geoscientists across the globe for a long time. In this study, seismic tomographic P-wave velocity variations are converted to density perturbations as the initial constraint values, and the four components of the GOCE Level-2 gravity gradient product are subjected to topographic effect correction, relief correction on the Moho and sedimentary layer boundary, and long-wavelength correction. In addition, the anomalous gravity gradient effect resulting from the uncertainty of the Moho and sedimentary layer depth is considered. Considering the effects of temperature variations in the lithosphere, which leads to the non-uniform density distribution, the anomalous gravity gradient effect resulting from the temperature variations is corrected at the first time in this area. Inverse computation is performed based on the preconditioned conjugate gradient algorithm. Because the Lagrange empirical parameter in the algorithm exhibits uncertainty in data inversion, the regularization parameter, which is typically an empirical parameter, is replaced with the value of the inflection point of the L-curve. The results show the followings. The non-uniform density distribution in the lithosphere is related not only to the composition but also to the internal temperature variations in the lithosphere. The density of the lithosphere in the NCC is significantly non-uniform in both the horizontal and vertical directions and exhibits a notable segment-wise spatial distribution pattern. Two tectonic units, namely the Taihangshan tectonic zone and the Linfen–Weihe graben, constitute a central gravity gradient transition zone. There is a significant difference in the density distribution on the two sides of this central transition zone.


Introduction
The lithosphere in the North China Craton (NCC) has suffered destruction since the Phanerozoic. Since the concept of lithospheric thinning in North China was proposed in the early 1990s, research on the destruction of the NCC has engendered increasingly heated debates. Particularly, in recent years, geophysicists, geochemists and geologists in China and elsewhere have conducted extensive research on the destruction of the NCC and achieved a clearer understanding of the deep structure, formation and evolution of the NCC (e.g. Zhao et al. 2001;Kusky and Li 2003;Zhai and Santosh 2011;Sun and Kennett 2017). However, currently, there is still a debate on a number of issues related to the NCC, such as the destruction mechanism, time of occurrence and geodynamic mechanism, based on the researchers' varying backgrounds and approaches (Wu et al. 2008). As one of the oldest cratons on Earth, the NCC is a notable example of a craton that has undergone modification and destruction. Studying the destruction of the NCC has provided a window for understanding the formation, evolution, stabilization and destruction of paleocontinents.
Currently, there are two mainstream views on what caused the destruction of the NCC, namely delamination (Gao et al. 2009) and thermal erosion (Zhu et al. 2012), both of which involve tectonic deformation and material distribution inside the lithosphere in the NCC. Tian and Wang Earth, Planets and Space (2018) 70:173 Therefore, it is of vital importance to acquire high-resolution structural data of the lithosphere. As seismic tomography technology advances, numerous researchers have extensively studied the crustal and upper mantle velocity structures in the NCC and its surrounding area and achieved some remarkable results (e.g. Hearn et al. 2004;Zheng et al. 2006;Chen 2009;Huang and Zhao 2009;Zhao et al. 2009;Tian and Zhao 2011). In contrast, there are relatively few studies on the deep underground density structure. The internal density structure of the Earth is the basis for studying the composition and internal structural deformation of the lithospheric crust. Fang (1996) obtained the density distribution in the lithosphere in North China through inversion of Bouguer gravity anomaly data using a constrained least-squares method. Wang et al. (2014) obtained a high-resolution three-dimensional density structure of the lithosphere in the NCC through inversion of Bouguer gravity anomaly data using an algebraic reconstruction inversion method with the seismic travel time as the constraint. Xu et al. (2016) obtained the density of the lithospheric crust as well as the depth of the lithosphere-asthenosphere boundary in the NCC through rapid joint inversion of gravity, geoid and topography data.
In recent years, techniques for measuring gravitational potential fields have evolved from static to dynamic methods, and measurement instruments have also gradually advanced from ground measurement systems to aerial and satellite measurement systems. As the second derivative of the gravitational potential, the gravity gradient can satisfactorily reflect a high-frequency density change and help increase the resolution of the density structure. Full tensor of the gravity gradient can effectively help increase the amount of observation data. Gravity gradient data can reflect the shape and size of an underground anomalous body from various perspectives. While there are many components of the gravity gradient, the inversion solution is non-unique. To constrain the non-unique problem of gravity gradient data inversion, seismic data are introduced, which, together with the multiple components of the gravity gradient, are subjected to a joint inversion to reduce the non-uniqueness of the inversion solution. In this study, based on a high-resolution P-wave velocity structure obtained from seismic tomography, a 3D initial density model is constructed for the lithosphere in the NCC based on the empirical velocity-density conversion relation for North China. The GOCE satellite gravity gradient data are processed. Considering that the anomalous gravity gradient effect is a comprehensive reflection of a myriad of factors and non-unique problem, the anomalous gravity gradient effect caused by the topographic relief at each underground layer interface is eliminated. In addition, by fully considering that the density distribution in the lithosphere is primarily collectively affected by its internal temperature and composition, the anomalous gravity gradient effect caused by the internal temperature of the lithosphere is also eliminated. Separation of the temperature effects from the residual gravity gradient anomalies provides constraints of the density changes due to variations in lithosphere and mantle composition. Thus, the anomalous gravity gradient effect caused by the nonuniform matter density distribution is obtained. Based on the anomalous gravity gradient components obtained, the preconditioned conjugate gradient inversion algorithm was adopted to calculate the horizontal density structure in the NCC. Because the Lagrange empirical parameter in the preconditioned conjugate gradient (PCG) inversion algorithm exhibits uncertainty in data inversion, this algorithm is improved. The distribution of the density structure at various depths in the NCC is obtained through inversion and is analysed and interpreted based on the available geological and geophysical data.

Initial density model data processing
The area between 33°N and 43°N and between 100°E and 120°E is selected as the study area. As shown in Fig. 1a, the NCC has a relatively diverse geological structure, with the North China Basin and the Tancheng-Lujiang fault zone as the primary geological units in the east, the Taihangshan tectonic zone and the Linfen-Weihe graben in the centre and the Ordos block and basins surrounding it to the west. In addition, the Yinshan-Yanshan orogenic belt, the Qilianshan orogenic belt, the Qinling-Dabie Shan Orogen and the Inner Mongolia-Northern Hebei Paleoproterozoic Orogen are also distributed in the surrounding area. Figure 1a shows the topography and the main structure in the area, which was modified from Kusky and Li (2003) and .
The gravitational field is the most basic and direct physical quantity reflecting the density and dynamics of the Earth's internal structure. Due to the non-uniqueness (See figure on next page.) Fig. 1 Topography and initial density model. a Topography and main tectonics in the area NCC, NCC: North China Craton, NSGL: North-South Gravity Line. b Absolute density at the depth of 10 km. c Absolute initial density at the depth of 25 km. d Absolute density at the depth of 42 km. e Absolute density at the depth of 60 km. f Absolute density at the depth of 80 km. g Absolute density at the depth of 120 km Tian and Wang Earth, Planets and Space (2018) 70:173 issue of the inversion solution and the fact that the kernel function decays with distance from the source, seismic data can be used to constrain the inversion calculation. Therefore, we acquired 3D 0.5° × 0.5° P-wave tomography data of the study area (Tian et al. 2009). By studying a large amount of manually collected seismic data of North China, Feng et al. (1989) proposed a wave velocity-density relation suitable for the study area: The collected P-wave velocity values were converted to density values using Eq. (1). Figure 1b-g shows the density values converted from the P-wave velocity values. According to the P-wave tomography data of the study area (Tian et al. 2009), values of 2840 kg/m 3 (10 km), 2980 kg/m 3 (25 km), 3320 kg/m 3 (42 km), 3320 kg/m 3 (60 km), 3320 kg/m 3 (80 km), and 3330 kg/m 3 (120 km) were used as the absolute density values of the layers at different depths. The differences between converted density and absolute density represent the density perturbations. We adopted the density perturbations as the initial values for density inversion.
The high-frequency signals induced by the geological characteristics of the shallow and middle layers are more significant in the gravity gradient data. In the previous research (Wang et al. 2014), there has been uncertainty to some extent at the depth of 150 km for the gravity data inversion. Accordingly, the bottom layer of the inversion model is determined to be at the depth of 120 km. Combined with the initial value for density inversion velocity perturbations, an initial inversion model was established for the inverse calculation based on the initial value. The inversion model is divided into six layers in the depth direction with the lower boundary at depths of 10 km, 25 km, 42 km, 60 km, 80 km and 120 km. In the horizontal direction, a 0.25° × 0.25° mesh with equidistantly spaced nodes is generated on each of the top four layers, and a 0.5° × 0.5° mesh with equidistantly spaced nodes is generated on each of the bottom two layers. This model partition method not only increases the abnormal effect of the deep blocks at the surface but also effectively reduces the number of parameters of the inversion model.

Preprocess of the GOCE data
As GOCE satellite gravity gradiometer consists of threeaxis accelerometers, two of the axis are high sensitive and the other one is not so sensitive. In the GOCE data, the precisions of the original gravity gradient observation components V xy and V yz are lower than those of the (1) ρ =    2.78 + 0.56(υ p − 6.0) (6.0 ≥ υ p ≥ 5.5) 3.07 + 0.29(υ p − 7.0) (7.5 ≥ υ p ≥ 6.0) 3.22 + 0.20(υ p − 7.5) (8.5 ≥ υ p ≥ 7.5) high-precision tensors V xx , V xz , V yy and V zz (Rummel et al. 2011;Yi et al. 2013). To use the gravity gradient data to detect anomalies, four high-accuracy anomalous gravity gradient components ( T xx , T xz , T yy , T zz ) provided by the GOCE satellite can be used to study regional internal structures and obtain more realistic and reliable structural models of the Earth (e.g. Bouman et al. 2015;Fecher et al. 2015;Li et al. 2017).
We collected the anomalous gravity gradient component (Fig. 2a) on the radius of the sphere (http:// goce.kma.zcu.cz/) by subtracting the orbital reductioncorrected observed gravity gradient component by the normal gravity gradient component. The preprocessed gravity gradient component with the spatial resolution of 10 arc-min at the average height of the orbit (250 km) was obtained, by correcting the 36th month (November, 2009 through October, 2013) GOCE GO_CONS_EGG_TRF_2 along-orbit data product for average orbital correction (Sebera et al. 2014). And with the Geodetic Reference System 1980 (GRS 80) reference ellipsoid as the reference, the normal gravity gradient component generated by the normal ellipsoid under the same average orbital radius (250 km) condition was calculated using the spherical harmonic equation for the normal gravity gradient in a local north-oriented coordinate system (with the z-axis pointing towards the Earth's centre, the x-axis pointing towards north and the y-axis pointing towards west) proposed by Šprlák (2012).
Obtaining more accurate mass and density anomalies especially helps with the geometrical interpretation, as they allow all of the local extrema to be better localized. Although the data signal becomes stronger and the noise is amplified, this problem can be resolved by improving the signal-to-noise ratio. We did not select spherical harmonic data because they are often regularized with respect to potential but not measured gradients.
The anomalous gravity gradient components were then extended downward by 240 km to obtain the anomalous gravity gradient component at a height of 10 km with the iterative spherical downward continuation method by Sebera et al. (2014) in Fig. 2b. As the highest point of topography is about 5 km in this area (Fig. 1a), for the space that mass-free outside the Earth surface, the Poisson integral equations can be valid for the downward continuation. The iterative downward continuation method (Sebera et al. 2014) was adopted. The method is based on an iterative scheme and the Poisson integral equations applied to gravitational data. The iterative scheme was theoretically justified by Sansò and Sideris (2016) and Mansi et al. (2018). All gravity gradient anomaly components were obtained at the same height of 10 km. Figure 2c represents the gravity gradient T zz obtained with the 220th-order model of the EGM2008 at an altitude of 10 km, which coincides with the trend and value of the measured gravity gradient after downward continuation (Fig. 2b). Compared with those in Fig. 2c, the local extremes in Fig. 2b are obvious and not conterminous. The statistical difference between the used GOCE data and a GOCE-only global model harmonic synthesis is shown in Table 1.

The topographic correction of the gravity gradient data
The gravity gradient tensors are not sensitive to the density variation vertically, and the lithospheric horizontal density variation is what we really concerned about. Therefore, we have to reduce the other known gravity gradient anomalies that can be determined by the other reliable measures. As the layer interface undulation data can be determined by the seismic wave velocity data, the topographical data determined by the other methods can meet the higher resolution, which are more reliable. Therefore, we adopted these established interface models (CRUST 1.0, ETOPO1) and reduced the corresponding gravity gradient anomaly effect from the corrected gravity gradient data.
The spatial resolution of the four preprocessed GOCE gravitational gradient grid tensors ( T xx , T xz , T yy , T zz ) at a noise-free satellite altitude is 10 arc-min (http:// goce.kma.zcu.cz/), the spatial resolution of the ETOPO1 is 1′ × 1′, and the spatial resolution of the CRUST1 is 1° × 1°. As the spatial resolutions of the three models The gravity gradient effect T zz after downward continuation to the altitude of 10 km. c The gravity gradient T zz with the 220th-order model of the EGM2008 at the altitude of 10 km. d The gravity gradient effect T zz caused by the topographic masses Table 1 The statistic difference value between gravity gradient after downward continuation and the 220thorder model of the EGM2008 dT represents the difference value of gravity gradient between preprocessed GOCE data and EGM2008 are different, the downloaded data are used with their original spatial resolution for the downward continuation and forward-modelling calculation, which makes the calculated results more reliable. We homogenized the same spatial resolution of 0.5° × 0.5° before calculating the remaining gravity gradient. For the higher spatial resolution of the data (preprocessed GOCE gravitational gradient grid data, topography data ETOPO1), we extracted data from the calculated results (Fig. 2b, d).
For the gravity gradient effect caused by the underground layer interface (CRUST1.0), we used the common kriging interpolation method to obtain the data needed for this calculation.
As the North China Craton is such a large area, it is more suitable to use the tesseroids instead of the planar prisms for forward modelling. We used the tesseroids as the modelling units under the spherical coordinates. We adopted the software named "Tesseroids" (Uieda et al. 2016), which helps us to implement the forward-modelling calculations (https ://tesse roids .readt hedoc s.io/en/ lates t/). The software "Tesseroids" can be used to convert the tesseroids in the spherical coordinates into the planar prisms in the Cartesian coordinates, with the converted planar prisms can be used, then the software can implement the forward-modelling calculations under the Cartesian coordinates which corresponds into the forward modelling calculations.
In addition, 1′ × 1′ ETOPO1 mesh topography data (Amante and Eakins 2009) were also acquired. The whole NCC was segmented to 1′ × 1′ mesh points. The longitude and latitude coordinates of the mesh points were used as the length and width boundaries of the upright cuboid unit, the centre of the mesh as the cuboid centre, and the topographic height and geoid as the cuboid upper and lower boundaries, respectively. A density of 2 670 kg/m 3 was used for the part of the upright cuboid unit above sea level, and a density of 1 030 kg/m 3 was used for the part below sea level. Using the prism unitbased forward-modelling method (Nagy et al. 2000), the reduction in the gravity gradient at a height of 10 km generated by the topographic mass was calculated. To reduce the effects of high-frequency topographic errors, the result from forward modelling of the topographic mass was filtered based on the previous research results (Jekeli and Rapp 1980). The anomalous gravity gradient effect in this area caused by the topographic masses was thus obtained (Fig. 2c). As there are too many figures of the four components, we only display the figures of one component T zz in the initial data processing.

The underground layer interface correction and uncertainty analysis
The topographically corrected gravity gradient anomalies are a comprehensive reflection of the relief at each underground layer interface, the internal temperature variations and the non-uniform distribution of material density in the lithosphere. To examine the non-uniform density of the lithosphere, the gravity gradient effects caused by factors other than density should be first eliminated. First, we collected 1° × 1° relief data of the sedimentary and Moho layer interfaces (CRUST 1.0) in the study area (Laske et al. 2013), as shown in Fig. 3a, e. Forward calculations were performed using the prism unit-based forward-modelling method (Nagy et al. 2000). When forward calculating the gravity gradient effect resulting from the relief at the sedimentary layer interfaces, the average depth of the sedimentary layer interfaces and their density difference were set to 4 km and − 200 kg/m 3 (Laske et al. 2013), respectively; when calculating the gravity gradient effect resulting from the relief at the Moho interface, the average depth of the Moho was set to 35 km, and the density difference between the crust and mantle was set to − 420 kg/m 3 (Laske et al. 2013). The values of the two mean depths are the average values of the sediment and Moho depths from CRUST 1.0. We calculated and adopted the mean density difference value of every point between the mean depth and the corresponding layer (sediment and Moho) from CRUST 1.0. Figure 3b, f shows the anomalous gravity gradient component T zz caused by the relief at the sedimentary and Moho layer interfaces at the height of 10 km. The crustal thickness data in CRUST1.0 are derived from active source seismic studies and receiver function studies. In areas where such constraints are missing, crustal thicknesses are estimated using gravity constraints (Laske et al. 2013). There are uncertainties in the sedimentary depth and Moho depth (Zheng et al. 2011;Laske et al. 2013;Shen et al. 2016;Abrehdary et al. 2017). As shown in Fig. 3c, g, we adopted the previously determined sedimentary (Shen et al. 2016) and Moho depth (Zheng et al. 2011;Shen et al. 2016) values of the NCC area. The published results were obtained using combined seismic and satellite gravity observations in this area, which have a 95% confidence interval. These data ( Fig. 3c, g) represent the difference between the CRUST 1.0 and regional results (Zheng et al. 2011;Shen et al. 2016) in this area. The gravity gradient effect caused by the uncertainty of the sedimentary and Moho depth was calculated (Fig. 3d, h) by the same forward-modelling method (Nagy et al. 2000), the figures show little variations in most areas and most gravity gradient values are around the numeric value 0 besides some small areas. And the detailed values of the effect caused by the different layer uncertainty depths were presented and are analysed in Tables 2 and 3 where the AGSU is the average gravity gradient effect caused by the uncertainty of the sedimentary layer depth, AGRV is the average gravity gradient residual value (Fig. 7), the R s represents the fluctuation of each average gravity gradient caused by the uncertainty of the sedimentary layer depth compared with the total remaining gravity gradient.
where the AGMU is the average gravity gradient effect caused by the uncertainty of the Moho layer depth, the R m represents the fluctuation of average gravity gradient caused by the uncertainty of the Moho layer depth compared with the total remaining gravity gradient. The percentage of R s and R m are within the reasonable range compared with the effect of uncertainties on residual lithosphere and upper mantle gravity by the parameter tests (Herceg et al. 2015). The effect of uncertainty layer depth on the inversion results is analysed in Section Antinoise ability test of PCG algorithm.

The gravity gradient effect correction of the mantle
To highlight the gravity gradient anomalies induced by the lithosphere, a long-wavelength correction was performed to correct the gravity gradient effect resulting from the non-uniform matter density of the lithosphere below 120 km. The relation between the buried depth of the field source and the order of the spherical harmonic function of the gravitational potential presented by Bowin et al. (1986) is used: where R is the Earth's radius; Z is the buried depth of the field source; and n is the order of the spherical harmonic function. The long-wavelength gravity gradient anomalies obtained from calculations using the second-to 44thorder spherical harmonic coefficients (Li et al. 2017) were deducted. The interfacial relief-corrected anomalous gravity gradient component at the height of 10 km was (4) Z = R/(n − 1)  obtained by subtracting the topographic relief-corrected anomalous gravity gradient component by the gravity gradient anomalies resulting from the relief at the Moho and sedimentary layer interfaces as well as the longwavelength gravity gradient effect at the same height, as shown in Fig. 4.

Gravity gradient effect caused by the temperature variations in the lithosphere
Previous research results (Kaban et al. 2003Mooney and Kaban 2010) show that the cause of density variations in the lithosphere is multifaceted and that the density distribution in the lithosphere is primarily affected by its internal temperature and composition. And the results by Herceg et al. (2015) show that the residual mantle gravity anomalies are mainly caused by a heterogeneous density distribution in the lithosphere and uppermost mantle, which depends on variations in lithosphere mantle temperature and composition. Separation of the temperature effects from the residual gravity anomalies provides constraints of the density changes due to variations in lithosphere and mantle composition. Lithospheric thinning in North China not only reflects the variation of the thickness, but also accompanies the change of the substance properties and the transformation of the thermal state. Both mainstream views, namely delamination (Gao et al. 2009) and thermal erosion (Zhu et al. 2012), would induce inhomogeneous density distribution of the lithosphere, but the characteristics of the density distribution induced by different mechanisms are various. For further discussions of the inhomogeneous density distribution caused by the component in the lithosphere, the gravity gradient effect caused by the temperature variations in the lithosphere should be separated.
While in previous studies on the study area, the nonuniform density of the lithosphere resulting from temperature variations and composition was not distinguished, and the gravity gradient effect caused by the temperature variations in the lithosphere was not differentiated; consequently, the density obtained through inversion contained the gravity gradient effect caused by the temperature variations in the lithosphere. In this study, we first collected 0.5° × 0.5° temperature data of the study area. Two groups of temperature data were used. The temperature data of the area above a depth of 70 km were obtained using the steady-state heat conduction equation with the temperature at the depth of 70 km as the boundary condition (An and Shi 2007). The temperature data of the area between the depths of 70 km and 120 km were obtained using the Euler state equation based on S-wave tomography data (Yang et al. 2013).
A temperature data model was established based on the temperature data of the area between the surface and the depth of 120 km. The temperature model was divided into eight layers along the 120 km depth at intervals of 15 km. Based on the average temperature T 0 in each layer, as well as the coefficient of thermal expansion α (Tesauro et al. 2009), the temperature variation-density variation relation for the lithosphere provided by Humphreys and Dueker (1994) and Kaban and Trubitsyn (2012) was used: where ρ 0 is the estimated density at a certain depth provided by Crust1.0 ( ρ 0 is set to the same value for model elements in the same layer); T = T 0 − T is the temperature variation in each model element, the normal symbol for the temperature variation is the Kelvin; and �ρ is the corresponding density variation in each model element. The global average coefficient of thermal expansion is α = 3.5 × 10 −5 K −1 (Kaban et al. 2003). The value of the coefficient of thermal expansion is primarily determined by the temperature and composition and generally varies within the range of (3.5 ± 0.25) × 10 −5 K −1 (Poudjom Djomani et al. 2001;Sobolev et al. 1996). Due to the limits of the study area and to achieve stable results, the same coefficient of thermal expansion the same as global average coefficient was used for the whole NCC. The temperature variation in each model element was converted to the corresponding density variation.
In addition, the corresponding kernel function was also calculated based on the temperature model. The gravity gradient effect caused by the temperature variations at a height of 10 km was determined using the prism unitbased forward-modelling method (Nagy et al. 2000). Figure 5 shows the flow chart of gravity gradient effect caused by the temperature variations in the lithosphere. Figure 6a, b shows the density variations at the depth of 70 km and 120 km in terms of the temperature variations. Figure 6c-f shows the four gravity gradient components caused by the temperature variations in the lithosphere. Here, Fig. 6f is used as an example for analysis. The maximum value of the temperature effect-induced gravity gradient component T zz can reach as high as 5.82E, which is on the same order of magnitude as the maximum value of the interfacial relief effect-induced gravity gradient component T zz . In addition, there is a certain correlation between the distribution of anomalous values of the gravity gradient caused by the temperature variations in the lithosphere and the surface structure. Overall, the anomalous values of the gravity gradient in the Qilian block are negative, whereas the anomalous values of the gravity gradient in the nearby Alxa block and Ordos block are both positive. The anomalous values of the gravity gradient in the Tancheng-Luyan fault zone and the Yinshan-Yanshan orogenic belt are both significantly negative. The positive anomalous values of the gravity gradient in the central Taihangshan Mountains are surrounded by the negative anomalous values of the gravity gradient in the eastern NCC and the northern Taihangshan Mountains. The values of the gravity gradient along the east and west sides of the gravity gradient zone undergo significant changes. Overall, the anomalous values of the gravity gradient on the east side of the gravity gradient zone are negative, and the distribution of anomalous values of the gravity gradient on the west side of the gravity gradient zone exhibits regional characteristics.
The components of the remaining gravity gradient anomalies (Fig. 7) were obtained by subtracting the gravity gradient effect caused by the temperature variations in the lithosphere from the previously calculated gravity gradient anomalies. Of the four components of the remaining gravity gradient anomalies, two components, namely T xz and T zz , undergo relatively significant changes. Figure 7d shows the cross section at the latitude of 33.5°N, 37.5°N and 41.5°N; Fig. 8a-c compares the gravity effect T zz of the residuals without temperature consideration (Fig. 4d), temperature variations gravity gradient effect T zz (Fig. 6f ) and the The density variations at the 120 km in terms of the temperature variations. c T xx . d T xz . e T yy . f T zz remaining gravity gradient effect T zz (Fig. 7d). The average percentage of the gravity gradient value T zz regarding the every point between the temperature variations ( Fig. 6f ) and the residuals without temperature consideration (Fig. 4d) is 26.31%. This indicates that the gravity gradient effect caused by the temperature variations in the lithosphere is a factor that should not be overlooked and that the non-uniform density distribution in the lithosphere is related not only to its composition but also to its temperature variations. The remaining gravity gradient anomalies are mainly caused by the non-uniform density of the lithosphere.

Inversion method
The inversion of the gravity gradient tensor is used to solve linear equations. Because the unknown variable m is much larger than the sampled data vector, the inversion is an underdetermined problem and the solution is not unique. Proper restrictions are needed for the objective function in order to narrow the solution area. We seek an objective function that enables us to account for prior information; the constructed objective function aligns with both the observations and additional geophysical constraints. Therefore, in the geophysical inversion theories, objective functions often consist of a data fitting function and a model objective function (Constable et al. 1987). The inversion problem is thus transformed to seek for a model m, which satisfies the data fitting function and minimizes the objective function.
Occam's inversion method, which was proposed by Constable et al. (1987) and DeGroot-Hedlin and Constable (1990), was employed in this study. The objective function can be written as: where µ represents the regularization parameter; the trade-off factors balance the data fitting function φ d and the model objective function φ m .
(6) minimize: φ = φ d + µφ m Fig. 7 The remaining gravity gradient anomalies components. a T xx . b T xz . c T yy . d T zz , the cross sections at the latitude of 33.5°N, 37.5°N, and 41.5°N are displayed Tian and Wang Earth, Planets and Space (2018) 70:173 The data fitting function can be written as: where m = m − m 0 is the modification of model parameter vector m from the initial model m 0 . The kernel matrix G is the linear projection operator from the model space to the observation space. d is the corresponding modification of the observation. σ i is the standard deviation of the error associated with the ith observation. The Lagrange multiplier is used as the regularization parameter for the conjugate gradient inversion algorithm. A large amount of computation space is needed for the Lagrange multiplier to solve a large matrix, and therefore, empirical values are normally selected as the optimal choice for the regularization parameter. Because the regularization parameter balances the data fitting function and the model objective function, a large value will (7) 1 σ 1 , 1 σ 2 , . . . , 1 σ N lead to a significant difference between model outputs and observations, whereas a small value will lead to unreliable inversion results due to the inconsistency between models and reality. Selecting empirical parameters for diverse observed data sets introduces uncertainty into the results, which weakens the applicability of Occam's inversion method. Considering the above-mentioned issues, we select the radius of the curvature inflection point of the L-curve as the regularization parameter. The L-curve is a trade-off standard based on the data fitting and objective function and thus is suitable to solve largescale problems (Hansen 1992). This method has been validated in a previous study (Tian et al. 2018a, b).
where ⌢ ρ = log(φ d ) , ⌢ η = log(φ m ) , and the symbols ′ and ″ represent the first and second derivatives of the function, respectively. Therefore, identifying the maximum curvature of the L-curve is considered to be the value of the regularization parameter in the inversion process as a whole. To construct the model objective function, the model roughness of Occam's inversion method helps us to identify a smooth model that agrees reasonably well with the observations. The roughness matrix (DeGroot-Hedlin and Constable 1990) is included in the model objective function in order to restrict the spatial structure of the model and to guarantee the continuity of image along the three directions. The 3-D model vector R is the sum of the squared first-order derivatives of the three components of the model vector m along the x, y and z directions. By meshing the model units into grids and replacing derivatives with finite differentials, the above formula can be written in the matrix form as: where R x , R y , and R z are the roughness matrices along the x, y and z directions, respectively.
Because the gravity data do not have fixed depth resolutions and kernel functions decay fast with depth, the inversion results close to ground surface cannot adequately represent the accurate location of the anomaly (Li and Oldenburg 1998). To overcome the tendency of density to concentrate at the surface, a depth weighting function was added to the model objective function to counteract the natural decay of the kernel function. Therefore, the depth weighting function (Eq. 12) is included in the model objective function so that the kernel function can accurately present the weights of anomaly at various depths.
We adopt the depth weighting function for ground surface gravity that was proposed by Li and Oldenburg (1998): where Z is the block unit depth, β and Z 0 are constant, Z 0 depends on the observation height and the cell size of the model discretization. When the exponent β = 3 and the form of the function is (Z + Z 0 ) −3 / 2 , which is consistent with the attenuation of the kernel function. (10) The inverse problem is formulated as an optimization problem where an objective function of the model is minimized subject to the constraints in Eq. (17) (Li and Oldenburg 1996). To find a particular model, Li and Oldenburg (1998) defined an objective function of the density and minimized that quantity subject to adequately fitting the data, the details of the objective function are problem dependent, but generally the model is required close to the reference model m 0 ( m = m − m 0 ), and the model is smooth in three spatial directions. Therefore, the model objective function (Li and Oldenburg 1998) that includes the roughness function and depth weighting function can be written as: replacing derivatives with finite differentials, the model objective function in the matrix form is (Li and Oldenburg 1998): where R i is the differential operator along each of the directions, and D represents the discretized matrix of the depth weighting function, α i is the weighting coefficients of the terms in the objective function, as the purpose of the objective function is to make the model required closed to the reference model m 0 . The detailed and adequate value of α i can be referred to Tian et al. (2018a, b).
Adding the model objective function to the objective function gives the following formula (Li and Oldenburg 1998): which can be written as the equation set: Page 15 of 22 Tian and Wang Earth, Planets and Space (2018) 70:173 Let A = G √ kW i represent the Jacobian matrix and b = d 0 , and Eq. (17) can be simplified to: We selected the four well-observed components with the same spatial resolution T xx ,T xz ,T yy and T zz as the observed data, which implies the Jacobian matrix A = Because the condition number of the Jacobian matrix is normally very large, to increase the convergence speed and the solution stability, Eq. (18) is rewritten as: where S is the preconditioned factor and is usually approximated as A T A −1 . By doing this, the eigenvalues of Eq. (19) will be concentrated along the diagonal and the condition number will be improved so that the iteration efficiency is improved (Pilkington 1997).
The PCG method is an iterative method that provides linear or nonlinear solutions. It is mainly used to solve for equations such as φ = 1 2 m T Am − m T b whose solution is m = A −1 b , equivalent to Eq. (18). This algorithm takes the gradient of initial point as the initial conjugate direction and seeks for the minimum of the objective function along this direction to improve stability. Following the computing process of the preconditioned conjugate gradient inversion, model m is solved. The fitting error is used as a criterion for the iteration. If the model m satisfies the predetermined difference β , the iteration ceases. Details regarding the PCG algorithm can be found in Pilkington (1997Pilkington ( , 2009) and Li and Yang (2011).
The fitting error in the inversion algorithm is where N is the number of observations, d obs i is the measured data, and d fwd i is the theoretical response from inversion models. Figure 9 shows the flow chart of the procedure of the PCG inversion algorithm.

Anti-noise ability test of PCG algorithm
The uncertainty of the gravity gradient effect due to the total uncertainty sedimentary and Moho interface depth is the − 0.07% ( T xx ), − 2.03% ( T xz ), − 3.92% ( T yy ), − 4.51% ( T zz ), respectively. To ensure the reliability of the inversion results, we designed two sets Fig. 9 Flow chart of the procedure of the PCG inversion algorithm of complex models (M1 and M2), and the 8% Gaussian noise was added to the four gravity gradient components in the synthetic model tests. To test the feasibility of the inversion method for several anomalies, model M1 consists of four anomaly blocks with different sizes, depths and densities. To demonstrate the effectiveness of the inversion method for large blocks, we designed two large blocks with density anomalies in Model M2.
Compared with Fig. 10b, e, Fig. 10c, f shows that there are only some same anomalies appear around the anomaly, and the density at the upper layer is a little higher than the anomalies in Fig. 10b, e, the inversion result reveals that the algorithm has robust anti-noise ability in Table 4.

Results
For the inversion of the density structure in the NCC, the maximum outer iteration number was set to 100; the absolute value of the defined RMS misfit β was set to 1 × 10 −2 . As shown in Fig. 11, after 19 outer iterations, the curve of defined RMS misfit tends to be horizontal, the variation is more and more slight between the values of defined RMS misfit, an defined RMS misfit of  3.5 × 10 −3 was obtained, at which point the inverse iterative computation was complete.

Discussion
As shown in Fig. 12, the inversion results show that the density of the lithosphere in the NCC exhibits notable non-uniformity in both the horizontal and depth directions and a significant segment-wise spatial distribution pattern. A gravity gradient zone composed of two tectonic units, namely the Taihangshan tectonic zone and the Linfen-Weihe graben, functions as a central transition zone. There is a significant difference in the distribution of the density anomalies between the two sides of the central transition zone. Hence, the NCC is divided into three areas, namely the eastern NCC, which consists of the North China Basin and Bohai Bay, the central NCC, which consists of the central transition zone, and the western NCC, which consists of the Ordos Basin and its surrounding area. The density distribution at various depths in each area is discussed separately.

Characteristics of the distribution of density anomalies in the eastern NCC
A notable spatial distribution pattern can be found for density anomalies in the whole eastern NCC, which adjoins to Bohai Bay. This is in agreement with the theory (Tian and Zhao 2011) that this area has suffered significant intense destruction.
There is a notable high-density anomalous body A1 in the North China Basin at a depth of 10 km (Fig. 12a) in the eastern NCC. The amplitude of the anomalies is significantly smaller at the top of the lower crust and the upper mantle (at depths of 25 km and 42 km, respectively). This anomalous body (Fig. 12b, c) is mainly distributed along the Beijing-Baoding line. This finding matches the Tangshan-Xingtai seismic zone and is in good agreement with the results concerning the P-wave velocity derived from a study conducted by Huang and Zhao (2004).
Density anomalies vary significantly between the two sides of the north-north-east-trending Tancheng-Lujiang fault zone in the eastern NCC. The anomalous body A2 on the west side of the fault zone exhibits notable characteristics, whereas the block on the east side of the fault zone exhibits no significant density anomalies. At the depths of 10-42 km (Fig. 12a-c), the anomalous body A2 on the west side exhibits significant low-density anomalies. At larger depths of 42-60 km, the maximum amplitude of the density anomalies in the anomalous body A2 on the west side decreases, while the body A3 exhibits no density variations. At the depth of 60 km, the density change of the anomalous body A2 and A3 is relatively gentle, and both the anomalous body A2 on the west side and the anomalous body A3 on the east side exhibit low-density anomalies, but the density characteristic of the body A3 is not evident as the A2 obviously. At depths of 80-120 km, the anomalous body A2 on the west side transitions from a low-density anomalous body to a high-density anomalous body, whereas the A3 still exhibits no obvious density variations. These results show that there is always a significant difference in the density distribution between the blocks on the two sides of the Tancheng-Lujiang fault zone and that the Tancheng-Lujiang fault zone has already penetrated the lithosphere.
In Bohai Bay (area A4), at depths of 42-120 km, particularly in the upper mantle, density anomalies are distributed along the Tancheng-Lujiang fault zone, low-density anomalies are more common, and the maximum amplitude of the density anomaly variations is large. Previous research results (Teng et al. 1997;Su et al. 2009) show that Bohai Bay experiences the combined action of the extension of the main fault at the eastern boundary of the Tancheng-Lujiang fault zone and potential plumes forming in the upper mantle. The gravity gradient effect caused by the upper mantle temperature has already been eliminated; however, density anomalies are still notably distributed in this area, and the maximum value of the low-density anomalies can reach − 50 kg/m 3 at a depth of 120 km (Fig. 12e). These results show that density anomalies in Bohai Bay are primarily caused by an extension of the main fault at the eastern boundary of the Tancheng-Lujiang fault zone.
Notable low-density anomalies can be observed at a depth of 10 km (Fig. 12a) in the transition area A5 between the Yanshan block and the Songliao Basin. As   Fig. 11 With the iterative calculation, residual mean square between the forward calculated theoretical gravity gradient and the gravity gradient measurements, the red line represents the given absolute value of 1 × 10 −2 the depth increases, continuous high-density anomalies can be found in the upper mantle. These density anomalies may be a result of the western extension of the highdensity mantle lithosphere beneath the Songliao Basin towards beneath the Yanshan orogenic belt (Tian and Zhao 2011).
Compared with the central and western NCC, there are no wide areas of the density anomalies but a notable spatial distribution pattern in the whole eastern NCC, which indicates that this area has suffered different degrees destruction. Furthermore, comparing Fig. 12e, f with Fig. 6a, b which has been marked on the figures, we find that the notable density distribution exists with the hightemperature distribution. Combined with our results and previous research results, Teng et al. (1997); Su et al. (2009) show that there are obvious plumes forming in the upper mantle in the eastern NCC, which means that the high temperature of the lithosphere and mantle in the eastern NCC is accompanied with the non-uniform density of the lithosphere resulting from the composition.

Characteristics of the distribution of density anomalies in the central NCC
The central NCC, composed of the central transition zone, mainly exhibits notable low-density anomalies. The anomalous bodies in this area exhibit significant segment-wise characteristics. The Taihangshan tectonic zone, with an overall north-east-south-west-trending strike direction, can be divided into three blocks, namely the southern block B1, the central block B2 and the northern block B3. At the same depth, the magnitude and distribution of anomalies both vary between these three blocks. At depths of 25-120 km, the left side of the southern block B1 adjoins to the positive-density anomalies in the western NCC, the central block B2 adjoins to the high-density anomalies in the eastern NCC and the northern block B3 adjoins to the high-density anomaly zone in the North China Basin in the east and adjoins to the Yinshan-Yanshan block in the west, forming a large low-density anomaly zone. The maximum values of the low-density anomalies in the southern block B1 and the central block B2 are more significant.
Overall, the Linfen-Weihe graben block is characterized by high-velocity anomalies alternating with low-velocity anomalies. Block B4 in the northern Linfen-Weihe Graben is mainly distributed in the Datong volcanic zone. Density anomalies in block B4 undergo significant changes with increasing depth. At a depth of 25 km (Fig. 12b), block B4 exhibits high-density anomalies. However, as the depth increases from 42 to 120 km, notable low-density anomalies can be found in block B4. This suggests that the Datong volcano is no longer active (Tian et al. 2009), but low-density anomalies are still distributed beneath this area.
At depths of 60-120 km, block B5 in the central NCC is surrounded by the low-density anomalies in the central Taihangshan block B2 and the northern Linfen-Weihe graben block B4. However, the high-density anomalies in block B5 become increasingly prominent as the depth increases and adjoin to the high-density anomaly zone in the North China Basin in the eastern NCC at the depth of 80 km (Fig. 12e). This indicates that the destruction of the lithosphere in the NCC may have already extended from the North China Basin to the central NCC.
At depths of 10-25 km, notable low-density anomalies can be found in the Yinshan-Yanshan block B6, and they are surrounded by high-density anomalies. These lowdensity anomalies are distributed in a relatively shallow zone along an orogenic belt and are mainly affected by the granitic bodies distributed along the orogenic belt (Hou et al. 2014). In addition, as the depth increases to 60 km (Fig. 12d), the aforementioned low-density anomalies adjoin to block B3 in the northern section of the Taihangshan tectonic zone and block B4 at the junction of the Linfen graben and the Weihe graben, forming a large low-density anomaly zone. Compared to previous findings (Wang et al. 2014), in this study, low-density anomalous blocks are more widely distributed in this area. The whole central NCC block is located in an orogenic belt, where the mantle temperature is relatively high (An and Shi 2007). This indicates that this area is significantly affected by both temperature and composition. A high heat-flow environment is present in this area, resulting in the upwelling of hot material in the deep asthenosphere, which modifies the lithospheric mantle.

Characteristics of the distribution of density anomalies in the western NCC
Tectonic deformation and intense magmatic activity have occurred in the Ordos area since the Mesozoic-Cenozoic, resulting in the formation of a special tectonic framework in Ordos and its surrounding area . The Ordos area is surrounded by two rift valleys, namely the Yinchuan-Hetao rift valley and the Shaanxi-Shanxi rift valley. Overall, the Ordos area displays as a stable and solid tectonic unit.
At a depth of 42 km (Fig. 12c), the central Ordos block C1 exhibits high-density anomalies, whereas the northern Ordos block exhibits no notable anomalies. Based on the receiver function , the central Ordos block is 41 km in thickness, and the northern Ordos block is 45 km in thickness. Based on this information, it is inferred that the central and southern areas of the Ordos are in the lithospheric mantle at a depth of 42 km. At depths of 60-120 km, high-density anomalies in the central Ordos block C1 transition to low-density anomalies, which are surrounded by high-density anomalies in the southern and northern areas of Ordos. The central Ordos block adjoins to the Yinshan orogenic belt in the north and the Qinling orogenic belt in the south. Because the characteristics of the distribution of density anomalies are consistent with those of the surface structure, it is inferred that the characteristics of the distribution of density anomalies on the south and north sides of Ordos are affected by the deep dynamic processes in the orogenic belts on the two sides.
The Qilian block is composed of block C2 in the southeast and block C3 in the north-west. Overall, a part of the Qilian block in the crust exhibits notable low-density anomalies. The upper boundary of this area corresponds to the Haiyuan fault zone. The direction of density anomalies in this area is basically consistent with the strike direction of the surface fault zone. At the junction of the Qilian block and the Ordos block, the distribution of density anomalies exhibits a significant abrupt change-from low-density anomalies to high-density anomalies. Generally, the anomalies in the Qilian block are north-westsouth-east-trending, which is a behaviour that is closely related to the regional tectonic framework-the movement of the Tibetan Plateau changes from northward to eastward upon being blocked by the stable Ordos block and Alxa block. At a depth of 60 km (Fig. 12d), block C2 in the north-east exhibits continuous high-density anomalies; in addition, as the depth increases, the distribution of high-density anomalies becomes increasingly prominent, and the amplitude of the high-density anomalies increases. In comparison, block C3 in the north-west still exhibits notable low-density anomalies, in which the amplitude decreases as the depth increases. The seismic sounding results show that the Moho depth in the Qilian block is less than 60 km (Tian and Zhao 2011). In addition, the tomography results show that intense lowvelocity anomalies at the bottom of the crust transition to intense high-velocity anomalies at the top of the Moho (Li et al. 2002). High-density anomalies found in block C2 in the south-east suggest that the block enters the lithospheric mantle at this depth.
The Alxa block C4 exhibits continuous notable highdensity anomalies from the crust to the mantle, indicating that this area is a stable tectonic unit. Part of the Qaidam block C5 in the crust exhibits low-density anomalies. At a depth of 42 km (Fig. 12c), the lower boundary of the density anomaly zone corresponds to the Riyueshan fault, and the direction of the density anomalies is basically consistent with the strike direction of the surface fault zone. At a depth of 60 km (Fig. 12d), this zone is characterized by low-density anomalies alternating with high-density anomalies that adjoin to the high-and low-density anomalies in the Qilian block. The northern Bayan Har block C6, which adjoins to the Qaidam block C5, exhibits intense low-density anomalies. At depths of 80-120 km, the direction of the density anomalies is no longer consistent with the strike direction of the fault zone. The Qaidam block C5 mainly exhibits high-density anomalies, whereas the low-density anomalies found in the northern Bayan Har block C6 become increasingly prominent as the depth increases.
Considered the whole central and western NCC, the destruction of the lithosphere in the NCC has extended from the North China Basin only to the regional area of the central NCC. Meanwhile, referred to the temperature variation in this area (Fig. 6a, b), the temperature of the lithosphere and mantle in western NCC is relatively low. These shallow responses to the deep characteristics show that the whole main region of the central and western NCC has not been destroyed yet the overall NCC still maintains stable attributes.

Conclusions
The GOCE data were processed with several corrections to obtain the remaining gravity gradient components, and the non-uniform density distribution in the lithosphere related to the internal temperature variations was assessed. The effect of the gravity gradient on the density variation in the lithospheric mantle due to temperature was separated from that due to composition. The data corrections were completed by analysing the uncertainty related to the models used.
The eastern NCC exhibits notable regional characteristics, the distribution of density anomalies varies significantly within this area and the results show that there is always a significant difference in the density distribution between the blocks on the two sides of the Tancheng-Lujiang fault zone. The Tancheng-Luyan fault zone has already penetrated the lithosphere and has extended to Bohai Bay.
The central NCC, composed of the central transition zone, mainly exhibits notable low-density anomalies. The anomalous bodies in this area exhibit significant segment-wise characteristics. The whole central NCC block is located in an orogenic belt, where the mantle temperature is relatively high; the results show that this area is significantly affected by both temperature and composition. A high heat-flow environment exists in this area, causing the upwelling of hot material in the deep asthenosphere, which modifies the lithospheric mantle.
In the western NCC area, the Ordos block and the Alxa block appear collectively as a stable and solid tectonic unit. The directions of the density anomalies in this part of the Qilian block and the Qaidam block in the crust are basically consistent with the strike directions of the surface fault zones corresponding to the upper and lower boundaries, respectively. It is inferred that the characteristics of the distribution of density anomalies on the south and north sides of Ordos are affected by the deep dynamic processes in the orogenic belts on the two sides. The anomalous density body in the Qilian block is northwest-south-east trending.