Integration of Residual Terrain Modelling and the Equivalent Source Layer Method in Gravity Field Synthesis for Airborne Gravity Gradiometer Test Site Determination

Round 1

Reviewer 1 Report

Comments and Suggestions for Authors

I am satisfied with the authors' efforts to improve the manuscript, and I have no further comments. I think the current version can be accepted.

Comments on the Quality of English Language

Minor editing of English language required.

Author Response

I would like to extend my heartfelt thanks for your time, expertise, and careful evaluation of our manuscript. We are genuinely appreciative of your positive assessment of the manuscript.

Reviewer 2 Report

Comments and Suggestions for Authors

Overview and general recommendation:

Airborne gravity measurement is a commonly used remote sensing method to obtain the underground density distribution. The establishment of an airborne gravity gradiometer test site in China is important to geodesy and geophysics. 

A new gravity field modeling method in the article is introduced, which combines the GGM long-wavelength gravity field, the gravity field implied by the RTM terrain of constant density, and the gravity field generated by density anomalies. This results in a better figure of the gravity field. The author further recommends that an area covering 5 km × 5 km area around BJS and WHS volcanos for building the airborne gravity gradiometer test site. Because of this, the current study is on a topic of relevance and general interest to the readers of the journal. 

On the one hand, I found the paper to be overall professionally written and much of it to be well described. I felt confident that the authors conducted rigorous experiments and reasonable analysis. On the other hand, I found some questions of the paper. I have little confidence in one analysis and came away with some questions to be able to recommend this paper for publication as it stands. Therefore, I recommend that a minor revision is warranted. I explain my concerns in more detail below.

Major comments:

1. Page6, lines 242. The distance between the gravity observation points in the article is approximately 500m. There are 213 measurement points that arranged evenly within an area of 15km×15km. In this case, the distance between measuring points is approximately 1000m, please check.

2. Page6, lines 240. There are 213 measurement points that arranged in the measurement area. The amount of gravity measurement data is small and the distances between the measurement points are slightly larger. What effect does this have on the subsequent calculations? Can you add a numerical simulation study to minimize its impact?

Minor comments:

3. Page1, lines 6. Please specify what kind of terrain leads to larger errors when calculating gravity using RTM with constant density conditions.

4. Page6, lines 234. On page 6, line 234, the gravity anomaly range is -39.55 mGal ∼ 28.54 mGal, while the color scale range in Fig. 1 (b) is about -30 mGal ∼ 30 mGal. And the statistical gravity anomaly signal minimum is 68.53 in Table 2 on page 10, while Fig. 3 shows that the minimum value is about -30 mGal, so please check.

5. Page7, lines 248. "GGM gGGM", check the syntax here.

6. Page7, lines 288. The density anomaly covers an area of 2278 - 297 km, which does not correspond to the previous statement, please check. 

 

Comments on the Quality of English Language

The paper has a clear topic, smooth flow of language, clear hierarchy, and accurate expression of ideas. Minor English touch-ups are recommended.

Author Response

Airborne gravity measurement is a commonly used remote sensing method to obtain the underground density distribution. The establishment of an airborne gravity gradiometer test site in China is important to geodesy and geophysics. 

A new gravity field modeling method in the article is introduced, which combines the GGM long-wavelength gravity field, the gravity field implied by the RTM terrain of constant density, and the gravity field generated by density anomalies. This results in a better figure of the gravity field. The author further recommends that an area covering 5 km × 5 km area around BJS and WHS volcanos for building the airborne gravity gradiometer test site. Because of this, the current study is on a topic of relevance and general interest to the readers of the journal. 

On the one hand, I found the paper to be overall professionally written and much of it to be well described. I felt confident that the authors conducted rigorous experiments and reasonable analysis. On the other hand, I found some questions of the paper. I have little confidence in one analysis and came away with some questions to be able to recommend this paper for publication as it stands. Therefore, I recommend that a minor revision is warranted. I explain my concerns in more detail below.

Thank you for your valuable feedback and kind suggestions. We appreciate your time and effort in thoroughly reviewing our work. We have carefully considered your comments and have made the necessary revisions and corrections to address the mentioned issues.

Major comments:

1. Page6, lines 242. The distance between the gravity observation points in the article is approximately 500m. There are 213 measurement points that arranged evenly within an area of 15km×15km. In this case, the distance between measuring points is approximately 1000m, please check.

Many thanks. Here it should be 1000m. We have made corrections in the revised manuscript. Please check.

2. Page6, lines 240. There are 213 measurement points that arranged in the measurement area. The amount of gravity measurement data is small and the distances between the measurement points are slightly larger. What effect does this have on the subsequent calculations? Can you add a numerical simulation study to minimize its impact?

Many thanks for your suggestions. As mentioned in our manuscript, the normal approach of gravity field recovery is GGM+RTM (constant density 2670 kg/m³). In order to overcome the large errors due to the constant density assumption, we involved the density anomalies which are obtained from the equivalent source layer method and sparse gravity measurements. By means of the Tikhonov regularization, sparse gravity data could transform into dense gridding density anomalies. The only problem is, density distribution details are dependent on the resolution of gravity data. With sparse data, near-surface density details may not be fully recovered.

As the density in this study is the equivalent surface density distribution rather than the real distribution of masses, it is difficult to define a criterion to estimate the effect of gravity inversion on density anomalies. The purpose of this study is for gravity field determination. Therefore, a numerical simulation is implemented to estimate its effect on gravity field determination. The simulated tests show that the equivalent surface density would achieve high accuracy in gravity field modeling. In the simulation test, a model composed of five masses is assumed.

Number	Density (g/cm3)	Size (km)	Position of Center (km)
1	1.0	2ⅹ2ⅹ2	(-5,-5,-3)
2	-1.0	2ⅹ2ⅹ2	(5,-5,-3)
3	-1.0	2ⅹ2ⅹ2	(-5,5,-3)
4	1.0	2ⅹ2ⅹ2	(5,5,-3)
5	0.5	10ⅹ10ⅹ8	(0,0,-12)

The gravity and gradients on the surface at a height of 0 are computed analytically with real density and provide the true values which denote the observations for inversion and reference values for comparison. The simulated gravity and gradients are computed with densities inversed using the equivalent source method. The comparison between simulated values and reference values gives insight into the quality of the gravity field estimated from derived densities. Our results show that the root mean square (RMS) value of differences between observed and simulated gravity is 0.029 mGal, which is ~1.02% of the mean value of gravity. The RMS value of differences between the observed and simulated gradient is 0.074 E, which is ~0.95% of the mean value of the gradient. This shows the errors could be acceptable in the mGal-level gravity determination and E-level gradient determination.

Still, in our numerical experiments, 193 of 213 gravity measurements, which is more than 90% of all available datasets, are used for the density anomaly inversion in this study. The left 10% of the measurements are used for validation. The recovered gravity field would be significantly improved from the recovery rate of 54.62% to 86.22%. This proved the method in this study is effective in areas with large density anomalies and sparse gravity measurements.

 

Minor comments:

3. Page1, lines 6. Please specify what kind of terrain leads to larger errors when calculating gravity using RTM with constant density conditions.

Thanks. This has been revised. Large errors would occur over rugged areas.

4. Page6, lines 234. On page 6, line 234, the gravity anomaly range is -39.55 mGal ∼ 28.54 mGal, while the color scale range in Fig. 1 (b) is about -30 mGal ∼ 30 mGal. And the statistical gravity anomaly signal minimum is 68.53 in Table 2 on page 10, while Fig. 3 shows that the minimum value is about -30 mGal, so please check.

The gravity anomalies on page 6, line 234 are observations, they range from -39.55 mGal to 28.54 mGal. In table 2 and Fig 3, they show the statistical information and features of computed gravity anomalies. The computed gravity anomalies are in a grid with a resolution of 50 m while the resolution of measured gravity anomalies in Fig. 1(b) is 1000 m. Therefore, the computed gravity anomalies in Table 2 and Fig 3 give more details. It shows the gravity signals in the crater while this was not measured. Therefore, it shows that there is a great difference in the minimum values of Fig 3 and Fig 1(b). According to your suggestion, the color scale ranges are revised in Fig.1 (b), Fig.3, and Fig.6. Please check.

5. Page7, lines 248. "GGM gGGM", check the syntax here.

Thanks. Here, indicates the long-wavelength signals provided by the GGM. In order to distinguish GGM and gGGM, we revised the manuscript as follows

“To study the detailed features of gravity field over WVF areas, the gravity and gradient signals are divided into three parts, the long-wavelength signals provided by GGM (), the high-frequency signals generated by residual terrain with constant density assumption of 2670 kg/m³ (), and the residual gravity field signals provided by subsurface density anomaly ().”

6. Page7, lines 288. The density anomaly covers an area of 2278 - 297 km, which does not correspond to the previous statement, please check. 

Thanks. This has been revised. It should be “2278 - 2297” km.

Author Response File: Author Response.docx

Reviewer 3 Report

Comments and Suggestions for Authors

Comments for author File: Comments.pdf

Author Response

The manuscript presents the results of interesting analyses regarding the structure of the local gravity field model. I believe this is valuable research material. However, I have some comments that I hope can help the authors improve this manuscript.

Thank you for your overall assessment of the manuscript. We have takenyou’re your suggestions into consideration and have made the necessary revisions and corrections accordingly. The revised manuscript now includes changes indicated in red for easy identification. We believe that your input has further improved the clarity and quality of the manuscript.

1. The introduction should include clear information that the gravity field model determined by the proposed method will be verified with measured data only in terms of the determined gravity. The values of the second derivatives of the potential determined from the model will not be verified, and the proposed solution only helps in selecting an area where a test field (indicated in the title) can be located. Although this is indicated in the summary, it should also be clearly written in the final parts of the introduction.

Thanks for your suggestion. In the introduction, the following text is included “Need to mention that, the promoted method in this study is validated through comparison with gravity measurements, which have affirmed its effectiveness in gravity determination over the Wudalianchi Volcano Field. When it comes to the determination of gradients, only the performance of ESL is validated through a numerical simulation test due to the unavailability of actual gradient measurements. Considering the fundamental theory of gravity field determination and the confirmed performance of ESL, this validates the method's capability in determining the gravity field over the study area and its validity in determining the test site.”

2. The lack of evaluation of the determined second derivatives slightly reduces the value of this study.

Thanks. We agree that it would be better if the promoted method could be validated through comparison with accurate gradient measurements. However, limited by technology, the gradiometer is under development. Therefore, it is difficult to obtain reference values for such validation. Based on the fundamental theory of gravity field, the gravity field components, including the gravity and gradients, are due to the features of the Earth’s mass distributions. The indicated long-wavelength signals could be provided with global gravity field models (GGM) and the short-wavelength signals could be provided with residual terrain modelling (RTM). Therefore, the normal approach of gravity field recovery is GGM+RTM (constant density 2670 kg/m³) over areas absent of actual measurements. The performance of GGM+RTM has been confirmed by various studies, such as Yildiz (2012), and Liu (2020). However, the density assumption of 2670 kg/m³ would involve large errors over our study area. In order to overcome such a problem, we involved the density anomalies which are obtained from the equivalent source layer (ESL) method and sparse gravity measurement. In our method, the gradient would be GGM+RTM (constant density 2670 kg/m³)+RTM(density of ESL). The performance of ESL in the gradient determination is validated in a simulation numerical test. In the simulation test, a model composed of five masses is assumed.

Number	Density (g/cm3)	Size (km)	Position of Center (km)
1	1.0	2ⅹ2ⅹ2	(-5,-5,-3)
2	-1.0	2ⅹ2ⅹ2	(5,-5,-3)
3	-1.0	2ⅹ2ⅹ2	(-5,5,-3)
4	1.0	2ⅹ2ⅹ2	(5,5,-3)
5	0.5	10ⅹ10ⅹ8	(0,0,-12)

The gravity and gradients on the surface at a height of 0 are computed analytically with real density and provide the true values which denote the observations for inversion and reference values for comparison. The simulated gravity and gradients are computed with densities inversed using the ESL method. The comparison between simulated values and reference values gives insight into the quality of the gravity field estimated from derived densities. Our results show that the root mean square (RMS) value of differences between the observed and simulated gradient is 0.074 E, which is ~0.95% of the mean value of the gradient. This shows the errors could be acceptable in the E-level gradient determination. Besides, the external validation experiments using gravity measurements showed that, GGM+RTM (constant density 2670 kg/m³)+RTM(density of ESL) recovered gravity field would be significantly improved from the recovery rate of 54.62% to 86.22%. Both experiments give insights into the better performance of the GGM+RTM (constant density 2670 kg/m³)+RTM(density of ESL) than GGM+RTM (constant density 2670 kg/m³) in gradient determination.

Liu J. Combining GGM and RTM to model the gravity gradient tensor[C]//EGU General Assembly Conference Abstracts. 2020: 4407.

Yildiz H. A study of regional gravity field recovery from GOCE vertical gravity gradient data in the Auvergne test area using collocation[J]. Studia Geophysica et Geodaetica, 2012, 56: 171-184.

3. In formula (2) there is a value ???? (harmonic correction), which does not appear in formula (1), but should appear there if formula (2) is a derivative of formula (1). The use of the ???? value is not explained clearly enough and I think that this part of the description requires improvement, especially since the problem of this correction is still valid.

Thanks. The HC for RTM gravitational potential is added in formula (1) and formula (3). The reason for requiring HC and its expressions are detailed described in the revised manuscript. Please check it.

4. It would be good to cite the sources of equations (1), (2) and (3).

Many thanks. The related references are added.

5. Similarly, equations (9) - (11) - please provide sources.

Many thanks. The related references are added.

6. L-curve method - should be described in more detail and cite appropriate source.

The text introduce the L-curve was added in the revised manuscript.

“The regularization parameter is determined by L-curve method. This method was initially introduced by Lawson and Hanson (1974) and later advocated by Hansen (1999) for addressing the linear inverse problem. This method leverages the characteristic behavior of the Tikhonov curve, which involves a series of inversions to construct the complete Tikhonov curve. It identifies the point on the curve that corresponds to the maximum curvature. In this study, we determine the suitable regularization parameter by analyzing the trade-off features between the regularized solution norm from the provided data as the regularization parameter varies, and  is adopted following Fig.2 (pointed by red arrow).”

Fig. 2 The L-curve between the regularized solution norm from the provided data with changing of the regularization parameter varies.

7. Line 229-233 – “not calibrated gravity anomalies” should be better explained. What point is the reference point? If it is the difference in gravity between points, it may introduce another, less misleading name instead of gravity anomalies.

Here, “not calibrated gravity anomalies” means that absolute gravity measurements were not implemented. The values are measured by relative gravimeters. In order to avoid confusion, the text is reformulated as

“It is important to mention that, the gravity anomalies in Zhang et al. (2017) and in Fig. 1(b) are not defined as general gravity anomalies in geodesy because absolute gravity values are not provided, instead, they are relative gravity measurements. This implies that a consistent bias would exist between measured and actual gravity anomalies.”

To keep consistency with previous studies, especially following the definition in Zhang et al. (2017), they are still named gravity anomalies in the manuscript. This will not affect the following numerical experiments and the conclusions. I hope you will be satisfied with such corrections.  

8. Line 255-256 – shouldn't the RMS value be greater than or equal to the average value? Please check this.

Many Thanks. This has been corrected. Besides, STD and RMS values are both provided in the revised manuscript.

“The differences between EGM2008 and terrestrial measurements have a mean value of 7.28 mGal, a standard deviation (STD) of 6.39 mGal and a root-mean-square (RMS) value of 9.98 mGal with respect to the mean value of 7.66 mGal, TD value of 6.75 mGal and RMS value of 10.20 mGal for the differences between EIGEN-6C4 and terrestrial measurements.”

9. Line 288 – The authors write: "The maximum depth of inversion is set to 15 km." From the previous description (e.g. equation (12)) it follows that the inversion will concern topography (according to RTM). Please explain in detail how defined is the layer for which the inversion was performed and the densities were determined.

Many thanks for the comments. We have add some more details about the equivalent source layer (ESL) method in the revised manuscript. As discussed in the first paragraph of Section 2.2: “we assume a finite ESL consisting of uniformly sized prism grids, which represent compact mass anomalies in the shallow crust. These grids are aligned with DEM units, and their equivalent densities can be utilized in the RTM process. Regardless of the distribution of terrestrial gravity data, densities can be obtained through linear regularization given the defined ESL range and size.”

And we discuss how to determine the ESL parameters in part of Section 2.2 and 3.2. The pre-determined parameters of ESL, like horizontal range, layer thickness, and buried depth, could influence the equivalent density results and RTM effects. The ESL should cover a larger area than the ground gravity data to cope with the boundary effects, the output density anomalies cover 1 km larger areas in this study. The choice of depth interval and regularization parameter could have an effect on inversion results. We define a 1 km thick ESL between 1-2 km depth by trail and test, and determine  by L-curve method.

10. I recommend marking the test areas, e.g. with letters that can be used in the text instead of coordinates.

Many thanks. The coordinates of test areas is deleted and test area is marked as “SW region” in the revised manuscript. Please check.

11. Below Table 1 is a group of formulas. I think that at least the formula for the value of ε (if not all formulas) should be defined in the text along with cite the source.

Thanks. That’s a good suggestion. All equations are referred to Yang et al., (2018). This has been given in the revised manuscript.

12. It would be good to include also standard deviations in the statistics provided. RMS depends on the average, so at large average values it does not best reflect the variability of a given quantity.

Thanks. That’s a good suggestion. The STD values are provided in Tables 1-3.

13. The manuscript provides detailed formulas for determining the RTM components. However, there is no appropriate description (including formulas) of how the global components (from the GGM model) and the heights from the global MERIT2160 model were determined.

Thanks. As you know, the GGM is generally provided with spherical harmonic coefficients. The generated gravity field components could be directly computed from them with specific equations. There are no more that could be discussed further. Displaying the equations will make the manuscript extremely long. Therefore, the formulas to compute GGMs components are not displayed in the manuscript. Instead, we added the related equations through reference and involved the software used for such computation in the revised manuscript, as following “In the following computations, long-wavelength signals generated by GGM are computed with Graflab software with Eqs. (56) and (65-71) in Buhca and Janak (2021) for gravity anomaly and gradients, respectively.” Besides, the information “the performance of EGM2008 and EIGEN-6C4 expanding to a degree and order of 2159” indicates the truncating degree for computation.

       For the MERIT2160, it was provided and generated by Hirt et al. (2019). They are in grid format. Briefly, the spherical harmonic analysis was made for MERIT DEM and obtained the spherical harmonic coefficients of heights, then the spherical harmonic coefficients of heights were transformed into grids based on the spherical harmonic synthesis. These procedures were completed by Hirt et al. (2019). In this study, we use the MERIT2160 product which is in DEM format. In order to avoid confusion, the related text is revised as “The smoothed model MERIT2160 DEM was created by Hirt et al. (2019). It was obtained based on the spherical harmonic analysis of MERIT DEM, the resulting spherical harmonic coefficients were truncated to a degree and order of 2159, and then was transformed into MERIT2160 DEM with spherical harmonic synthesis.”

14. Please also provide the coordinate system in which the analyses were carried out and the results are presented. It would also be good to present the location of the analysed area on a smaller scale map.

Thanks. The coordinates are defined in the Gauss-Kruger Projection coordinate system. This information is added in the revised manuscript.

 

Author Response File: Author Response.docx

Round 2

Reviewer 3 Report

Comments and Suggestions for Authors

I am satisfied with the corrections made and the answers/comments provided by the Authors. I have no further comments.