Remaining non-isostatic effects in isostatic-gravimetric Moho determination—is it needed?

× 1 ◦ . To reach this goal, we present a new model, named MHUU22, formed by the SGGUGM2 g ravitational ﬁeld, Ear th2014 topog raphy, CR UST1.0 and CR UST19 seismic crustal models. Particularly, this study has its main emphasis on the RNIEs on gravity and Moho constituents to ﬁnd out if we can modify the stripping gravity corrections by a speciﬁc correction of the RNIEs. The numerical results illustrate that the RMS differences between MHUU22 MD and the seismic model CRUST1.0 and least-squares combined model MOHV21 are reduced by 33 and 41 per cent by applying the NIEs, and the RMS differences between MHUU22 MDC and the seismic model CRUST1.0 and least-squares combined model MDC21 are reduced by 41 and 23 per cent when the above strategy for removing the RNIEs is applied. Hence, our study demonstrates that the speciﬁc correction for the RNIEs on gravity disturbance is signiﬁcant, resulting in remarkable improvements in MHUU22, which more clearly visualize several crustal structures.


I N T RO D U C T I O N
Isostasy is a key concept in the Earth sciences describing the state of equilibrium (or mass balance) to which the mantle tends to balance the mass of the crust in the absence of external disturbing forces.One of most important application of the isostasy is to map the Moho discontinuity (or Moho), which separates the boundary between the Ear th's cr ust and upper mantle.The Moho constituents (i.e.Moho depth or crustal depth; MD and Moho density contrast; MDC) are generall y computed b y two techniques of seismic and gravimetric.
The seismic observations used in compiling global Moho models are typically sparse, and therefore interpolation of global MDs, particularly over areas without adequate seismic data, frequently yield unrealistic results with large uncertainties.Accordingly, over large parts of the world with a limited coverage of seismic data, the gravimetric or combination of seismic and gravimetric data is fruitfull y of fered (Abrehdar y 2016 ).The g ravimetric-isostatic approach uses the isostatic gravity anomaly (or disturbance) to model the Moho constituents.This anomaly/disturbance is the gravity signal 2066 C The Author(s) 2023.Published by Oxford University Press on behalf of The Royal Astronomical Society.This is an Open Access article distributed under the terms of the Creative Commons Attribution License ( https://cr eativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.corresponding to the isostatic compensation by the Moho constituents.Hence, the anomalous density structures not only within the crust but necessarily within the whole lithosphere (called the non-isostatic effects, NIEs) must be modelled and removed from the obser ved g ravity signal, that is, the obser ved g ravity data must be corrected in two main w ays, namel y for the gravitational contributions of mass density variations in different layers of the Earth's crust such as ice and sediment layers (stripping gravity corrections), as well as for the gravity contribution from deeper mass variations below the crust.Most of the NIEs are applied by removing lower deg ree har monics of the g ravity field (in this study har monics at degrees n < n 0 = 10), which are assumed to be caused by deep Earth density variations, and by directly adding modelled stripping gravity corrections ('additive corrections') for ice, sediment, bathymetry, etc. (see below Section 2).Ho wever , one cannot expect to be able to remove all NIEs in this way.For instance, except for lacking stripping corrections there are also Moho constituents caused by techtonic motions in the Earth.Subsequently, the main purpose of this paper is to investigate if one can reach even better results by applying a technique to remove what we call remaining non-isostatic effects (RNIEs) on gravity.To do so, we will first recover the Moho constituents on global scale to a resolution of 1 • × 1 • by employing the so-called Vening Meinesz-Moritz (VMM) technique (Sj öberg 2009 ; Sj öberg & Bagherbandi 2011 ), and we will then compare the outcomes with and without applying RNIEs using least-squares combined seismic and gravimetric-isostatic MD and MDC models (Sj öberg and Abrehdary 2022 a, b ) as well as the seismic model CRUST1.0(Laske et al. 2013 ).
Up to now, several different isostatic models have been offered for mapping the MD and MDC before and after correcting for NIEs, and it is not specified which model is most suitable, which makes it difficult to judge what is the best method in a certain situation.Also, these studies usually used the same seismic model for determining the NIEs and also for the e v aluation of their effects on the Moho constituents.For instance, the Bagherbandi & Sj öberg ( 2012 ) calculated NIEs on crustal thickness using CRUST2.0 in Fennoscandia and they concluded that the final result for the Moho could be considered as bias-free.See also Bagherbandi et al . ( 2013 ).
Unlike previous studies, we now include additional information such as stripping gravity corrections as well as using the seismic crustal model CRUST19 (Szwillus et al. 2019 ) in estimating the RNIEs, leading to a new Moho model named MHUU22 containing both MD and MDC by using data sets of the SGGUGM2 gravitational field (Liang et al. 2020 ), Earth2014 topography (Hirt & Rexer 2015 ), CRUST1.0 and CRUST19 seismic crustal models.This model will be e v aluated b y comparing it with CRUST1.0,MOHV21 and MDC21 models (Sj öberg and Abrehdary 2022 a, b ).The MHUU22 model includes uncertainty estimates.(It is worth noting that error estimates of MD and MDC are of great importance to any usage of Moho models.Ne vertheless, the e v aluation of such uncertainties are usually overlooked.)

The concept of the VMM method
The VMM inverse problem of isostasy aims at estimating the MD in such way that the isostatic compensating attraction A C ( P ) of the crustal mass totally compensates the Bouguer gravity disturbance on the Earth's surface, δg TBISN B , which is corrected for the gravitational contributions of topograph y/bath ymetry and density contrasts of the oceans, ice and sediments and also RNIEs.Assuming that the MDC ( ρ) is known, Sj öberg ( 2009 ) showed that this constraint leads to a nonlinear Fredholm integral equation in MD (denoted D below), whose second-order solution becomes: where R is sea level radius, σ is the unit sphere.D 1 is a first-order approximation, defined by the spectral form: where f nm is the ne gativ e spherical harmonic coefficient of δg TBISN B di vided b y the gravitational constant times density contrast and Y nm is the fully normalized spherical harmonic of degree n and order m .n 0 is set to 10 under the assumption that the lower degrees are all created by gravity disturbances in the deep mantel and Earth's core.Sj öberg & Bagherbandi ( 2011 ) showed that the above Fredholm equation could also yield similar solutions to the MDC for known MD and to the product of MD and MDC.See also Sj öberg & Bagherbandi ( 2017 , p. 306).

A D D I T I V E C O R R E C T I O N S T O T H E G R AV I T Y D I S T U R B A N C E
As already presented in the Introduction, gravity data contain signals of the whole spectrum of the Earth's structure, that is, signals from the geometries and density distributions of topography and bathymetry, ice and sediment la yers as w ell as variable density layers in the mantle and core/mantle topographic variations.Here we assume that the long-wavelength contribution to the gravity field, say to degree and order 10, is related to the mantle as well as coremantle boundar y topog raphy (Sj öberg 2009 ).To strip the g ravity data caused onl y b y the geometry and density contrast of the Moho interface, all aforementioned signals contributors to the gravity data should be removed by applying the so-called stripping gravity corrections as well as RNIEs.

Stripping gravity corrections
The main input data set used in the VMM method is the the refined Bouguer gravity disturbances, δg R1 B , that is, free-air gravity disturbances corrected for topograph y, bath ymetry, ice thickness and sediment basins (stripping gravity corrections δg TBIS ; see Tenzer et al. 2015 ).
The refined Bouguer gravity disturbance δg R1 B is thus realized from the SGGUGM2 free-air gravity disturbance δg SGGUGM2 FA by: (3)

Remaining NIEs
It needs to be mentioned to the reader that in general the crust is not in complete isostatic equilibrium, and the observed gravity data are not only generated by the topographic/isostatic masses, but also from masses in the deep Earth interior, that lead to the NIEs.Typically, the NIEs are divided into two groups.The first group comprises those effects that, if they were corrected for, make the isostatic gravity anomaly/disturbance consistent with the Moho in isostatic balance.In reality, the Moho is not only formed by isostasy, but there are other causes, which affect the crustal thickness estimation.Hence, in contrast to above, the second group are those gravitational effects that are caused by the deviation of the Moho from its isostatic model.According to Sj öberg ( 2009 ) the major part of the long wavelengths of the geopotential are caused by density variations in the Earth's mantle and core/mantle topography variations.Such NIEs could be the contribution of different factors, such as crustal thickening/thinning, thermal expansion of mass of the mantle (Kaban et al. 2004 ), Glacial Isostatic Adjustment (GIA), plate flexure (Watts 2001 , p. 114) and effect of other phenomena (Tenzer et al. 2009 ).This implies that this contribution to gravity will lead to systematic errors/NIEs of the computed Moho topography and density contrast.Hence, the NIEs should also be corrected on the isostatic gravity disturbance.
It needs a tr ustwor thy seismic Moho model for its implementation.It is worth mentioning that, there are already some global crust models based on seismic methods.CRUST1.0 is the most frequently used crustal model, which is based on a database of crustal thickness data from active source seismic studies as well as from receiver function studies.More recently, Szwillus et al. ( 2019 ) developed a global crustal thickness model and velocity structure from geostatistical analysis of seismic data (i.e.CRUST19).In this section, a short discussion on which seismic model should be used for best estimating the RNIEs is presented among the two selected candidate models CRUST1.0 and CRUST19.To do so, we first compare their power spectra to degree 180 in Fig. 1 , showing that the spectra are close in the low degrees (say, to degrees 10), but then the power of CRUST19 decreases significantly slower than that for CRUST1.0.
In Fig. 2 (a), we plot the difference between CRUST1.0 and CRUST19, w hereas in F ig. 2 (b), the difference between two models di vided b y the standard de viation of CR UST19 is plotted, sho wing that the largest differences occur in South America and Asia, otherwise the two models are overall similar.However, disagreements are observed in the continent-ocean transition which could be explained by the fact that CRUST1.0 uses a smooth transition over the continental margin, while CRUST19 uses a sharp jump from continental to oceanic crust leading to a systematic difference between CRUST19 and CRUST1.0 on the edges of the continents (see Szwillus et al. 2019 ).The assessment of the global Moho models of CRUST1.0 and CRUST19 is not an easy task, but the RMS differences between the two are typically of the order of 4.5 km and reduces to 3.3 km when exempting S. America.
Next, we will rely on CRUST19 model for estimating the RNIEs, as CRUST19 is based on a USGS database of crustal seismic studies including little additional information, following a design philosophy that tries to respect the data as much as possible.This model interpolates MD and average -wave velocity of the crystalline crust using a non-stationary kriging algorithm.A major advantage is that the model is relati vel y straightforw ard in deri ving the uncertainties.
If we assume that CRUST19 is known and correct, the gravity effect of the RNIEs can be estimated from its difference to the obser ved g ravity effect: where The spherical harmonic coefficients c j nm ( j = VMM, CRUST19) in eq. ( 4b) are generated using the following expression (see Bagherbandi et al. 2017 ): Here ρ e ( ≈ 5 .5 g cm −3 ) is the Earth's mean density, R( = 6371 × 10 3 m ) is the Earth's mean radius (which approximates the geocentric radius of the geoid surface), D VMM , D CRUST19 and D 0 are the MDs of VMM, CRUST19 and their global mean value (approximately 23 km based on CRUST19 model), respecti vel y. ρ is the (assumed) constant MDC and the spherical harmonic coefficients ( ρ( nm are defined for the following products of the MDC and MD: . Taking into consideration the gravitational contribution of the crust density heterogeneous, the NIEs are applied to the isostasy gravity disturbance δg I .The refined Bouguer gravity disturbance in eq. ( 3) is then rewritten as: (5)

D ATA
As we already mentioned, the main goals of this investigation are to compute RNIEs on gravity by using the CRUST19 seismic crustal model and to use the VMM model in imaging the Moho constituents on a global scale to a resolution of 1 • × 1 • .To that end, we deliver a new model, named MHUU22, including MD and MDC estimated from the SGGUGM2 gravitational field, Earth2014 topography, CR UST1.0 and CR UST19 seismic crustal models.SGGUGM2 is a new high-resolution Ear th g ravity field model from satellite g ravimetr y, satellite altimetr y and Ear th Gravitational Model 2008 (EGM2008)-derived gravity data based on the theory of the ellipsoidal harmonic analysis and coefficient transformation (Liang et al. 2020 ).
The refined Bouguer gravity disturbances were generated using the SGGUGM2 Earth gravitational model coefficients and DTM2006 topographic model (Pavlis et al. 2007 ), complete to degree and order 180 of spherical harmonics, and the spherical harmonics of the normal gravitational field were computed according to the parameters of the reference system GRS-80 (Moritz 2000 ).Then these gravity disturbances were corrected for the density variation of the oceans, ice sheets and sediment basins (i.e.δg R1 B ; see Tenzer et al. 2015 ) and also the RNIEs (Bagherbandi & Sj öberg 2013 ), (see eqs 3 and 5 ).As already stated, we follow the method presented by Bagherbandi & Sj öberg ( 2012 ) to estimate RNIEs on gravity using the recent seismic crustal thickness in CRUST19 model by comparing the gravimetric and seismic data in the frequency domain.CRUST19 is also used to obtain the nominal MD ( D 0 ).clearly shows significant effects of applying the RNIEs, most obvious along the mid ocean ridges (see also Sj öberg & Abrehdary 2021a ).In the open seas, the RNIEs could also be determined based on empirical formulas related with the inverse crustal ageing (e.g.Seton et al. 2020 ) and/or the lithospheric thermal-pressure variations (e.g.Bagherbandi et al. 2017 ).
We finally calculate the correlation coefficients for the oceanic lithosphere age (M üller et al. 2008 ) versus RNIE, MD and MDC are −0.56,−0.49 and −0.35, respecti vel y.

MD and MDC estimation along with their uncertainties
In this section, we employ the VMM model to estimate the Moho constituents (i.e.MD and MDC) on global scale in a set of 1 • × 1 • blocks.In this w ay, dif ferent hetero geneous data are used, including the global Earth gravity field model (e.g.SGGUGM2), the global topography model (e.g.Earth2014) and the global seismic crustal model (e.g.CRUST1.0 and CRUST19).In the following, we map the free-air gravity disturbances produced by the SGGUGM2 model, RNIEs deducted by CRUST19, and the refined Bouguer gravity disturbances after applying RNIE corrections ( see Table 1 ).
In order to calculate the standard errors of the MHUU22 MD/MDC, we follow the strategy explained by Sj öberg & Abrehdary ( 2021b ) based on an optimal least-squares combination of seismic and gravimetric-isostatic models of MD or MDC at a resolution of 1 • × 1 • .
Finally, the MHUU22 MD and MDC models are geophysically interpreted along with their standard errors in Table 2 and Figs 4  and 5

Assessment of MHUU22 MD/MDC
The main aim of this section is to assess the quality of MHUU22 MD/MDC with respect to some other global models.To do so, we compare MHUU22 before/after correcting for RNIEs (PMHU22/MHUU22) with CRUST1.0,MOHV21 and MDC21.3 demonstrates significant improved agreements in the comparisons when applying the NIEs, in particular in the RMS.In the comparisons with MOHV21 and MDC21, there are reductions of 40 and 23 per cent, respecti vel y.
In the upper part of Table 3 (also in Figs 6 a and b), the MDs estimated by MHUU22 are compared with those in CRUST1.0 and MOHV21, respecti vel y, showing that the RMS dif ferences between MHUU22 and CRUST1.0 and MOHV21 are 3.1/4.6and 1.6/2.7 km before/after considering the RINE corrections.In the lower part of Table 3 (along with Figs 7 a and b), we compare the MDC of MHUU22 with CRUST1.0 and MDC21.In summary, the RMS difference between MHUU22 MDC and CRUST1.0 is 64.2/109.3kg m −3 and the comparison with MDC21 yields 51.6/67.2kg m −3 before/after applying the RINE corrections, respecti vel y.

. C O N C L U S I O N S
The main focus of this paper was to estimate the RNIEs on gravity using a seismic Moho model (CRUST19).For this purpose, the refined Bouguer gravity disturbance was primarily reduced for gravity of topography, density heterogeneities related to bathymetry, ice, sediments and other crustal components using stripping gravity corrections, and it was further corrected for RNIEs of nuisance gravity signals from mass anomalies below the crust due to crustal thickening/thinning, thermal expansion of the mantle, Delayed Glacial Isostatic Adjustment (DGIA) and plate flexure.We then delivered a new gravimetric-isostatic MHUU22 model including the MD and MDC using the VMM gravimetric-isostatic model from the refined Bouguer gravity disturbance and Earth2014 topographic data.
We also validated the RMS of differences between MHUU22 MD and the seismic model CRUST1.0 and MOHV21 when applying the RNIE corrections, yielding reductions of 33 and 41 per cent, respecti vel y.Similarl y, the RMS dif ference of MHUU22 MDC and the seismic model CRUST1.0 and MOHV21 reduced 41 and 23 per cent, respecti vel y.Hence, the specific corrections for the RNIEs on gravity disturbance is significant, resulting in expected noteworthy changes in MHUU22 towards the optimal least-squares models.For instance, this implies that only the RNIE corrected Figs 4 (a) and 5 (a) clearly show small MDs and MDCs that are typical along the mid ocean ridges.Also the regional MD maxima in NW China (about 59 km) and central Finland (about 57 km) are only visible after correcting for the RNIE (Fig. 4 a).
) Here, c RNIEs nm , c VMM nm and c CRUST19 nm are the spherical harmonic coefficients of the gravity disturbances of the RNIEs, VMM and CRUST19, respecti vel y.
Figs 3 (a) and (b) plot the free-air gravity disturbances computed by SGGUGM2 and the RNIEs estimated by CRUST19, whereas Figs 3 (c) and (d) depict the refined Bouguer gravity disturbances stripped by the bathymetry, ice, sediment variations before and after applying the RNIEs.Comparing Figs 3 (a) and (c) show considerable modifications over oceans due to the bathymetric stripping gravity corrections, and also in central Greenland and Antarctica due to the correction for ice density variation and also for sediment stripping over sediment basins.Ho wever , a comparison of Figs 3 (c) and (d)

Figure 6 .
Figure 6.(a) The MD differences between the MHUU22 and CRUST1.0 and (b) the differences between the MHUU22 and MOHV21 (unit: km).
(c).Figs 4 (a) and (b) plot the MD estimated from the MHUU22 model after and before applying the correction for the RNIEs, and similar comparison for MDC are shown in Figs 5 (a) and (b) .As one would already expect, in most areas there are remarkable changes.For example, the shallow MDs with large MDCs along mid-ocean ridges become more pronounced, and significantly deeper MDs with larger MDCs are also seen in the Himalayas and the S. American west coast.Similar effects can also be seen in areas with post-glacial rebound (e.g.Fennoscandia and Hudson Bay in Canada).In areas with existing huge ice masses (Antarctica and Greenland) the MD is also significantly modified.Figs 4 (c) and 5 (c) map the MHUU22 MD and MDC standard er rors, showing er rors smaller than 3 km and 30 kg m −3 in oceanic regions, while for continental regions the errors reach 6 km and 70 kg m −3 .
MOHV21 (Sj öberg & Abrehdary 2022a ) is a-n MD model based on an optimal least-squares combination of five global seismic and gravimetric-isostatic models at a resolution of 1 • × 1 • , while MDC21 (Sj öberg & Abrehdary 2022b ) is a new global MDC model based on a weighted least-squares combination of three available MDC models at a resolution of 1 • × 1 • .Table

Table 1 .
Statistics of the free-air and refined Bouguer gravity disturbance estimated from SGGUGM2 model.STD is the standard deviation of the estimated quantities over the study area, δg SGGUGM2 FA is free-air gravity disturbance, δg TBIS is Bouguer gravity disturbance corrected for topography, bathymetry, ice thickness and sediment basins (i.e.stripping gravity corrections), δg RNIEs is gravity disturbance corrected for RNIEs and δg TBISN R is gravity disturbance corrected for stripping gravity corrections as well as RNIEs.

Table 2 .
Statistics of the MD and MDC estimated from MHUU22 model.