A local lithospheric structure model for Vietnam derived from a high-resolution gravimetric geoid

High-resolution Moho and lithosphere–asthenosphere boundary depth models for Vietnam and its surrounding areas are determined based on a recently released geoid model constructed from surface and satellite gravity data (GEOID_LSC_C model) and on 3ʹʹ resolution topography data (mixed SRTM model). A linear density gradient for the crust and a temperature-dependent density for the lithospheric mantle were used to determine the lithospheric structure under the assumption of local isostasy. In a first step, the impact of correcting elevation data from sedimentary basins to estimate Moho depth has been evaluated using CRUST1.0 model. Results obtained from a test area where seismic data are available, which demonstrated that the sedimentary effect should be considered before the inversion process. The geoid height and elevation-corrected sedimentary layer were filtered to remove signals originating below the lithosphere. The resulting Moho and lithosphere–asthenosphere boundary depth models computed at 1ʹ resolution were evaluated against seismic data as well as global and local lithospheric models available in the study region. These comparisons indicate a consistency of our Moho depth estimation with the seismic data within 1.5 km in standard deviation for the whole Vietnam. This new Moho depth model for the study region represents a significant improvement over the global models CRUST1.0 and GEMMA, which have standard deviations of 3.2 and 3.3 km, respectively, when compared to the seismic data. Even if a detailed geological interpretation of the results is out of scope of this paper, a joint analysis of the obtained models with the high-resolution Bouguer gravity anomaly is finally discussed in terms of the main geological patterns of the study region. The high resolution of our Moho and lithosphere–asthenosphere boundary depth models contribute to better constrain the lithospheric structure as well as tectonic and geodynamic processes of this region. The differences in Moho depth visible in the northeast and southwest sides of the Red River Fault Zone confirmed that the Red River Fault Zone may be considered the boundary between two continental blocks: South China and Indochina blocks. However, no remarkable differences in lithosphere–asthenosphere boundary depth were obtained from our results. This suggests that the Red River Fault Zone developed within the crust and remained a crustal fault.


Introduction
The Vietnam region, lying on the southeast part of the Eurasian tectonic plate, has a complex tectonic history inherited from the tectonic evolution of Southeast Asia. Molnar and Tapponnier (1975) and Jolivet et al. (1990) have proposed that the tectonics of Southeast Asia can be entirely related to the India-Eurasia collision initiated 80 My ago, which has influenced the structure of the Earth far beyond the Himalaya-Tibet orogeny. The study region is considered to be mainly on two tectonic blocks (Fig. 1a, b): the Indochina (or Sundaland) terrane and the South China terrane. The Indochina terrane includes east Thailand, southwest part of Vietnam, Cambodia, a large portion of Laos, and small parts of Malaysia and Indonesia. These tectonic features have played an important role on the geological evolution of Vietnam and on the formation of economic resources of the country. The boundary between South China and Indochina blocks is located along the Ma River Fault Zone (MRFZ) (Helmcke 1985;Hutchison 1989;Findlay 1997;Findlay and Trinh 1997;Lepvrier et al. 1997Lepvrier et al. , 2011Faure et al. 2014;Metcalfe 2013). However, some authors suggested that the Red River Fault Zone (RRFZ) is the major geological discontinuity, which separates South China from Indochina block (Leloup et al. 1995(Leloup et al. , 2001Gilley et al. 2003;Trần et al. 2013;Dinh et al. 2018). The northeastern part of Vietnam, belonging to the South China terrane, is situated immediately southeast of the eastern Himalayan syntaxis. The Hoang Lien Son Mountain Range is the highest relief of this region, with its Fansipan peak at the height of about 3,140 m above sea level. This northern part of Vietnam has a complicated tectonic setting, dominated by active faults, such as the Red River Fault (RRF), Chay River Fault (CRF), Lo River Fault (LRF)-all within the Red River Fault System (RRFS), Dien Bien Phu Fault (DBPF), Da River Fault (DRF), Ma River Fault (MRF), and Son La Fault (SLF) (Fig. 1a, b). Most of these active faults are strike-slip faults and extend into the northwestsoutheast (NW-SE) direction (Huang et al. 2009). The nature of displacements of the Da River terrane, located between the RRF and MRF, provided key constraints on the accommodation of the southward extrusion of the Indochina Peninsula. The Da River Zone represents a key area for the understanding of the geodynamic evolution of the Indochina and South China blocks during the Indonesian orogeny (Lepvrier et al. 2004). The RRFZ, located in the center of the northern part of the study region, extends 1000 km from the eastern Himalayas through southernmost China to northern Vietnam and is considered to be mechanically connected with Fig. 1 a Tectonic situation in the study region and panel in the bottom right corner is a sketch regional tectonic map, modified from Gilley et al. (2003), and b topography and major active faults of the Vietnam region the Indian-Eurasian collision. This fault zone is considered as a pure crustal structure (Jolivet et al. 2001;Searle 2006). However, some authors suggested that the strikeslip RRFZ is a lithospheric structure (Gilley et al. 2003;Leloup et al. 2007). According to the earthquake catalog compiled by the Institute of Geophysics (IGP)-Vietnam Academy of Science and Technology (VAST) (Nguyen 2005), 90% of earthquakes in Vietnam are located in the northwest part (Huang et al. 2009;Tran 2012). Even if Vietnam is considered as a low-seismicity region when compared to the major tectonic activity regions in the world, it remains tectonically active, as indicated by the occurrence of moderate earthquakes (Huang et al. 2009). The earthquake catalog also indicated 30 earthquake active areas, with magnitudes of 5.5-6.8. Such earthquakes have caused heavy damage in the northwest as well as offshore from central Vietnam.
Classically, one of the fundamental constraints to be considered for improving the knowledge of the lithospheric structure is the determination of the Moho discontinuity. In the study region, several isostatic-gravimetric Moho depth surfaces were developed for the northern part and the sea of Vietnam (Bui 1983;Dang 2003;Braitenberg et al. 2006;Dinh 2010;Nguyen et al. 2018). However, the results of these studies differ considerably due to the lack of data and independent constraint conditions . Moreover, these results are largely outdated now due to the gravity model used, which was based on relatively sparse and irregular ground measurements. Recently, other studies were carried out using seismic data for determining the Moho depth in Southeast Asia, e.g., , Li et al. (2008), Noisagool et al. (2014, Yu et al. (2017) and Su et al. (2018). However, these studies did not provide a detailed Moho depth model for this region. The lithosphere-asthenosphere boundary (LAB) is another fundamental piece of information in plate tectonics. Several types of geophysical data are classically used in LAB determinations (e.g., seismic, thermal, electromagnetic, gravity, and rheological) making the resultant LAB depth models significantly different (Eaton et al. 2009;Artemieva 2011). The structure of the lithosphere beneath Vietnam and its surrounding areas is even less well understood than the Moho topography. The isostatic-gravimetric approach is commonly used to constrain such models.
The availability of accurate Global Gravity field Models (GGM) and very high-resolution Digital Terrain Models (DTM) gives easy access to gravity and height data, making the gravity inversion method an expedient tool to study lithospheric structure. Recently, through the local isostatic assumption, Fullea et al. (2006) developed a new method to model Moho and LAB depth using geoid heights derived from a GGM, in combination with elevation data from a DTM. A refinement of this approach was presented in Fullea et al. (2007) by considering the lithospheric mantle density variation with temperature instead of using a constant one. This method was successfully used in the determination of the Moho and LAB depth of the Arabia-Eurasia collision (Iran) , Central Asia (Robert et al. 2012), southern India (Kumar et al. 2014), the Iberian Peninsula (Torne et al. 2015), Africa (Globig et al. 2016), and Central Eurasia (Robert et al. 2017). All of these studies used the geoid height and elevation data derived from the GGMs EGM96 at 30ʹ resolution (Lemoine et al. 1998) or EGM2008 (Pavlis et al. 2012) at 5ʹ resolution and the DTMs ETOPO2 at 2ʹ resolution (Smith and Sandwell 1997) or ETOPO1 (Amante and Eakins 2009) at 1ʹ resolution, respectively. Presently, for the purposes of height system modernization, most countries are improving their local or regional geoid model. Thanks to the use of the land gravity data collected at national level, the current accuracy and resolution of local geoid models are then significantly improved. These local geoid models are likely to provide more accurate determination of the lithospheric structure than with a GGM only at least for the regions where the fill-in data are used as this study region. In 2019, a gravimetric-only quasigeoid model (GEOID_LSC) at 5′ resolution was determined for Vietnam and its surrounding areas based on new gravity data (Vu et al. 2019). This model has a standard deviation (STD) of 8.3 cm when compared to GNSS/leveling data in Vietnam, while with EGM2008, it is 29 cm . The use of the GEOID_LSC model is thus expected to significantly improve the accuracy of the Moho and LAB depth in Vietnam. This model is the basis for this research.
The advantage of the isostatic hypothesis is that it removes a significant part of the influence of deep density heterogeneity to obtain a much more balanced anomaly (Kaban et al. 2016). Nevertheless, the reality is not so simple, because the isostatic anomalies may contain significant contributions due to other effects, e.g., the effect of sedimentary layers as indicated by Evans and Crompton (1946). Therefore, sedimentary corrections should be considered before hypotheses of isostatic equilibrium. Presently, models of sedimentary basins are sufficiently reliable to calculate their effects, e.g., CRUST1.0 (Laske et al. 2012). The approach has been developed to refine the isostatic anomalies. In Bouguer gravity anomalies, the sedimentary basins, as lying on the upper part of the crust, are well known to produce a significant gravity signal at all wavelengths. Such effect must be removed before the gravity inversion process to estimate Moho depth (Wienecke 2006). Sediment corrections were thus applied in many studies (Kaban et al. 2004(Kaban et al. , 2016Braitenberg and Ebbing 2009;Tenzer and Hamayun 2010;Tenzer et al. 2012;Reguzzoni et al. 2019). Similarly, the effect of sedimentary basins in the elevation data should be also removed before applying the local isostatic assumption in combination with the geoid height to visualize the lithosphere structure. Sandwell and Renkin (1988) used the method presented by Crough (1983) to calculate the sediment correction for bathymetry. This study indicated the importance of the sediment correction on bathymetry in improving its relationship with geoid height. Despite the importance of eliminating the sedimentary effect, it was not considered in the method proposed by Fullea et al. (2007). In Globig et al. (2016), the sedimentary thickness was used to account for lateral changes in the crustal density, and then, the method of Fullea et al. (2007) was used to determine the lithospheric structure. However, no improvement in the accuracy of the Moho depth was obtained when compared with seismic data. In this study, the topography-corrected sedimentary basin model is proposed to refine the model presented by Fullea et al. (2007). Unfortunately, seismic data are not available on the sedimentary basins of the study region, e.g., Khorat Plateau or Gulf of Tonkin, to validate the proposed method. Therefore, a calculation is performed on a test area in Central Eurasia where reliable seismic data are available to assess the efficiency of applying the sedimentary corrections to topography in computing the Moho and LAB depths. Then, the Moho and LAB depths will be determined for the study area including Vietnam using the high-resolution local gravimetric geoid model based on the proposed method. The validation of our models is performed using the seismic data as well as the global and local models available in the study region. A comparison with the Moho depth calculated using the classical approach from the inversion of Bouguer gravity anomalies is also presented. Finally, we further discuss the implications of this Moho and LAB depth model with respect to the regional geological structure of the study area.

Methodology
The Moho and LAB depths are calculated by combining elevation above mean sea level and geoid data together with thermal analysis in a 1D approach proposed by Fullea et al. (2007). This method assumes local isostasy with a compensation depth below the LAB, and a fourlayer model: water, crust, lithospheric mantle, and asthenosphere. According to Fullea et al. (2006), the relation between the Moho and LAB depths is then as follows: (1) where z c is the depth of the crust-mantle boundary (Moho), z L is the depth of the lithosphere-asthenosphere boundary (LAB), ρ c is the mean density of the crust, ρ m is the mean density of the lithospheric mantle, ρ w is the density of the sea water, ρ a is the density of the asthenosphere, L 0 is the depth of the free asthenospheric level, and E is the elevation. The geoid height (N) can be determined using the following relation: where z max is the depth of the compensation level, G is the universal gravitational constant (6.674 10 -11 m 3 kg −1 s −2 ), and g is the Earth's surface gravitational acceleration (9.81 m s −2 ).
As absolute densities are used in this approach, an integration constant N 0 is needed to adjust the zero level of the geoid heights (Fullea et al. 2006). N 0 is determined by applying Eq. (2) to a selected lithospheric reference column (z cREF , z LREF ) that eliminates geoid height.
The expressions (1) and (2) form a system of equations for the depths of the Moho (z c ) and LAB (z L ) and for the densities of the crust (ρ c ) and lithospheric mantle (ρ m ). Thus, there are four variables and only two constraints. If densities in the crust and lithospheric mantle are known for the whole model, then variable depths for the Moho and LAB can be retrieved from Eqs. (1) and (2). Fullea et al. (2007) proposed a refinement of this approach by considering a lithospheric mantle density variation with temperature and a linear density gradient from upper and lower crustal densities for the crustal layer instead of using a constant one. According to Parsons and Sclater (1977), the lithospheric mantle density depends on the temperature as follows: where α is the linear coefficient of thermal expansion (K −1 ), T a is the temperature at the LAB, and T m (z) is the temperature at depth z in the lithosphere. This equation is solved with boundary conditions considering fixed temperature at the surface (T s ) as well as at the LAB. The vertically average value of the lithospheric mantle density can be determined by integrating Eq. (3) between z c and z L : where T mh is the temperature at the Moho. The geoid and DTM data, which are used for deriving geoid height and (2) elevation data, respectively, contain signals representing all masses within the Earth. Therefore, to study the large scale lithospheric structure only (i.e., Moho and LAB interfaces), a filtering procedure (low-pass filter) must be applied to the data. Thus, a MATLAB program was developed to filter geoid height and elevation data, and subsequently estimate the Moho and LAB depth.

Gravity data and anomalies
The terrestrial gravity data used in our geoid and gravity anomaly computation were collected from the IGP-VAST, the Vietnam Institute of Geodesy and Cartography (VIGAC), and the Bureau Gravimétrique International [BGI; (Bonvalot 2020)]. The total number of gravity points is 31,102 (Vu et al. 2019) and their distribution is shown in Fig. 2a. A set of 19,267 land gravity points was collected from the IGP, but these surveys have been done between 1961 and 1984 for the purpose of geological survey, exploration geophysics, and mineral prospecting when positioning was of poor quality. In this study, the cleaned-up IGP data are used (see Vu et al. (2019) for details). Most of the country has been resurveyed from 2003 to 2011, through the project "Measurement and Improvement of Vietnam National Gravity Data". This dataset comprises 10,940 points, including an absolute gravity network of 11 base reference stations with an accuracy of better than ± 5 μGal and their tie points, 29 first-order gravity stations with an accuracy of better than ± 15 μGal, and their tie points and 92 third-order gravity points and more than 10,000 detailed points measured from 2005 to 2009 (green points in Fig. 2a). Finally, the land gravity data set was also complemented by a set of gravity data provided by BGI for Vietnam (895 points) and surrounding areas (229 points in Cambodia, and 351 points in Thailand). The land gravity measurements were used in combination with fill-in data where no gravity data existed, e.g., mountainous or coastal areas. A mixed model [GOCE DIR-R5 (Bruinsma et al. 2014) plus EGM2008 (Pavlis et al. 2012)], called the mixed DIR/EGM model, up to degree/order 2159, plus RTM effects, and gravity field derived from altimetry satellites, were used to provide the fill-in information over land and marine areas, respectively. The land gravity measurements were interpolated to grid nodes using GRAVSOFT GEOGRID (Forsberg and Tscherning 2008) based on the Least-Squares Collocation (LSC) method. To avoid step-like discontinuities when merging the fill-in and measurement grids, a transition procedure was applied (see Vu et al. (2019) for details). The grid Fig. 2 a Distribution of land gravity data used in this study: red dots are the VIGAC points, absolute points are indicated by triangles, green dots are the IGP points, and blue dots are obtained from the BGI database; b, c corresponding maps of free-air and complete Bouguer gravity anomalies, respectively of free-air gravity anomalies is shown in Fig. 2b. They vary from − 213.0 to 175.8 mGal with mean and STD of − 12.8 and 26.9 mGal, respectively. Large negative values, − 150 to − 200 mGal, are seen in the northwestern region. The short wavelengths of the free-air gravity anomalies correlate strongly with topography. This map of free-air anomalies shows the same small details as seen in the topographic map (Fig. 1b).
The Bouguer gravity anomalies are obtained by removing the effect of the topographic and bathymetric masses from the free-air anomalies. The planar approach is used to compute topographic reductions. The mixed SRTM model, described in the section "Topography and sedimentary effect", is used to compute the topographic reductions. Terrain corrections were calculated within a radius of 20 km for the inner zone using the detailed DTM grid at 3ʹʹ and 200 km for outer zone using the coarse DTM grid at 3ʹ. Constant values for crustal density (ρ c ) of 2670 kg/m 3 and for sea water (ρ w ) of 1030 kg/m 3 were used. The estimated terrain corrections vary from 0 to 48.6 mGal. Large terrain corrections, at the 20-30 mGal level, are in the northwest region. The Bouguer plate corrections were also computed using the elevation derived from the mixed SRTM model. Finally, the Bouguer gravity anomalies, shown in Fig. 2c, are obtained by subtracting the Bouguer plate and adding the terrain corrections to the free-air anomalies. They vary from − 377.6 to 166.7 mGal with mean and STD of − 38.0 and 49.5 mGal, respectively. A northwest-southeast (NW-SE) trend is visible in Fig. 2c with large negative anomalies up to 150-200 mGal in the northwest. This is due to the Hoang Lien Son Mountain Range, with its Fansipan peak at the height of about 3140 m above the sea level. Generally, the most negative Bouguer anomalies correspond to the highest summits, which are located in the northwest. The large negative anomalies suggest a large root beneath this area. According to the Airy-Heiskanen model, the crust will thicken there. In contrast, there is an antiroot in the southeast corner of the study area, which is over sea with depths down to 2500 m. Figure 2c shows that, even if the standard deviation of the free-air anomalies is smaller (22.6 mGal), the Bouguer gravity anomalies are smoother, i.e., they vary less over small distances. The advantage of the Bouguer gravity anomaly is that the topographical masses outside the geoid are almost completely removed. This topographic effect dominates the free-air gravity anomalies. After its removal, the Bouguer gravity anomalies represent mainly the effects of lateral density variations and the variations of the Moho topography. Therefore, they are commonly used to determine and interpret the structure and depth of the crust.
In this study, the grid and map of Bouguer gravity anomalies will be used to (1) convert height anomalies to geoid height and (2) interpret in relation with Moho depth.

Geoid height and filtering
An accurate local gravimetric quasigeoid, called GEOID_LSC, was determined using the available gravity data in Vietnam by Vu et al. (2019). Geoid undulations are used to determine Moho and LAB depth, so first the quasigeoid must be converted to the geoid. The grid of Bouguer gravity anomalies and the mixed SRTM model were used to compute the differences between height anomaly and geoid undulation for converting the quasigeoid to the geoid (Hofmann-Wellenhof and Moritz 2006), called GEOID_LSC_C. This conversion procedure was not applied in recent studies using the geoid height inversion approach despite the fact that GGMs provide the height above the quasigeoid (Hofmann-Wellenhof and Moritz 2006). The geoid is an equipotential surface of the Earth's gravity field, while the quasigeoid is not. According to Hofmann-Wellenhof and Moritz (2006), the quasigeoid-to-geoid correction can reach − 0.1 m if the Bouguer gravity anomaly is − 100 mGal and elevation is 1000 m. This correction varies from − 0.771 to 0.020 m in the study region.
The Moho depth is often estimated by modeling gravity and/or geoid under the assumptions of local or regional isostasy. These assumptions need to be verified and this can be done using seismic data, even when sparsely distributed. First, correlations between Moho depth derived from seismic data, geoid height, and elevation need to be determined. Figure 3a shows for the study area the correlations between Moho depth from 24 receiver functions analysis ) and the geoid height derived from GEOID_LSC_C and EGM2008. Thicker crusts correspond to low geoid heights. The geoid height derived from GEOID_LSC_C (R 2 = 0.44) agrees better than EGM2008 (R 2 = 0.37) with the seismic-derived Moho depths. Hence, a more accurate Moho depth beneath Vietnam is likely to be determined from GEOID_LSC_C model, even if the correlations between geoid height and Moho depth remain quite low. The explanation for such low correlation is that all masses within the Earth contribute to all wavelengths of the observed geoid. According to Bowin (1983), the degree 10 geoid, i.e., wavelengths of 4000 km and longer, provides an estimation of the combined contributions from mass anomalies at the core-mantle boundary region and from the deep mantle and deeper parts, i.e., 400 km depth, of plate convergent zones. Removing the long wavelengths of the geoid height eliminates the effect of sub-lithospheric sources from geoid undulations. However, different maximum degree values have been used in past studies to filter the geoid signal for different regions. Table 1 shows the degree/order (d/o) used in several studies.
To determine the optimum d/o for Vietnam, the removal procedure was carried out on the GEOID_ LSC_C model from degree 8-11 of GOCE DIR-R5 in steps of 1 degree. The residual geoid undulations were then compared to Moho depth derived from 24 receiver function analysis in Vietnam (Fig. 3b). The correlation increases significantly after removing the effect of sub-lithospheric sources from long wavelengths of geoid undulations, and Fig. 3b indicates that the optimum degree for removing the low wavelengths is 9 (R 2 = 0.64). After this removal, higher geoid heights correspond to thicker crust. The correlation between the Bouguer gravity disturbance and Moho depth (R = − 0.54) was given by Sjöberg and Bagherbandi (2017) for South America to verify that the Bouguer gravity can be used to estimate the Moho depth. Here, the high correlation (R = 0.80) Fig. 3 Relation between Moho depths beneath seismic stations and corresponding geoid undulation retrieved from a GEOID_LSC_C and EGM2008 and b GEOID_LSC_C after removing low degrees between the geoid undulation and Moho depth suggests that the Moho geometry can also be estimated with geoid undulations.
Traditionally, we can remove the long wavelength components by subtracting the geoid undulations calculated from low d/o spherical harmonics. However, according to Sandwell and Renkin (1988), a too large amplitude of residual geoid undulation between northeast and southeast Hawaii was caused by the sharp cut-off of the spherical harmonic coefficients at d/o 10. This rigorous elimination leads to side lobes in the spatial structures of the truncated fields. To minimize this effect, high-pass filtering with a 1D-Gaussian function, or a gentler cutoff was used. In the latter case, coefficients are rolled off smoothly. According to 1D-Gaussian truncation, prior to summing the series, coefficients of degree n were multiplied by: when n reaches σ n , the coefficients are reduced by 0.6. In practice, σ n was set to 9 and the coefficients for degrees ranged from 2 to 25 have been multiplied by w n calculated using Eq. (5).
In gentle truncation, coefficients of degree n were multiplied by: and the rigorous truncation of coefficients at n = 9 is replaced by gently cutting the spherical harmonic series from n min = 2 to n max = 16 (i.e., the coefficients for degrees ranged from 2 to 16 have been multiplied by the function w n , Eq. (6)).
The resulting geoid heights from the GEOID_LSC_C model and their residuals are shown in Fig. 4. The geoid (5) w n = exp −(n − 2) 2 2(σ n − 2) 2 , w n = n − n min n max − n min 4 − 2 n − n min n max − n min 2 + 1, height varies from − 35.3 to 16.7 m. The low harmonic coefficients (d/o = 9) from the GOCE DIR-R5 were used to remove the long wavelengths of GEOID_LSC_C. The residuals vary from − 11.9 to 0.3 m, from − 7.9 to 1.5 m, and from − 8.8 to 1.0 m when using rigorous truncation, 1D-Gaussian function, and gentle truncation, respectively. The amplitude (maximum minus minimum) after rigorous truncation (12.2 m) is much larger than after Gaussian (9.4 m) and gentle (9.8 m) truncations. The amplitudes of 1D-Gaussian and gentle truncation are at the same level with slightly better results with the Gaussian function, which is selected for computing Moho and LAB depth.

Topography and sedimentary effect
In Vu et al. (2019), a mixed Shuttle Radar Topography Mission (SRTM) model was constructed as follows: over land areas, the 90 m resolution SRTM3arc_v4.1 (Farr et al. 2007) was used as the detailed DTM. The 15ʹʹ resolution Digital Bathymetry Model (DBM) SRTM15arc_ plus (Becker et al. 2009) was used over sea, and after re-gridding to 3ʹʹ, it was then merged with SRTM3arc_ v4.1 using the full-resolution coastline in Generic Mapping Tools (GMT) (Wessel and Smith 1998). This mixed SRTM model is used in the following computations.
In this study, the topography-corrected sedimentary basins model is used to refine the model presented by Fullea et al. (2007) in determination of the lithospheric structure. To compute the sedimentary effect on the topography, a density-depth relationship that incorporates a sediment compaction model (Ebbing et al. 2007) is used. The function defining the density of the sedimentary basin (ρ s ) according to the depth (d) was given as follows:  Globig et al. (2016) where 0 is the effect of porosity, ρ f is the fluid density, ρ g is the grain density, and b 1 , b 2 are the depth-decay parameters. According to Kaban et al. (2016), the sedimentary effect was corrected on the elevation (E) as follows: where E c is the elevation-corrected sedimentary layer and ρ topo is the density of topography. ρ topo = 2670kg/m 3 is used for the topography above the sea level, and for the sea areas, the density is replaced by The sedimentary effect is evaluated using the global sediment thickness CRUST1.0 model, which is shown in Fig. 5a. The thickness of the sediment layer ranges from 0 to 8 km for this region. The parameters used to determine the effect of sediment on topography are listed in Table 2.
According to Fullea et al. (2006), Fullea et al. (2007), and Globig et al. (2016), the high-frequency components (wavelength of 100 km and shorter) need to be removed from the elevation data set to avoid mapping of unrealistic signals, which are related to flexural support of topographic loads, into the modeled crustal and lithospheric topography. To that end, a low-pass filter is applied to eliminate these short wavelengths from the elevation data. A simple filter was tested by averaging values of all grids being within a cut-off wavelength of the filtered grid (100 km). The disadvantage of this simple approach is that all grid nodes contribute equally. It is more realistic to use a weighting system where grid nodes close to the filtered grid node contribute more than distant nodes. Applying such a 2D-Gaussian low-pass filter instead of simple filter will lead to a smoother grid. The procedure to correct for the influence of the sedimentary layer and filtering of elevation data is as follows: • First, the 1 × 1 0 sediment thickness derived from CRUST1.0 model was re-gridded to 1ʹ to compute sediment corrections for the topography using Eq. (8). The 3 arc-second mixed SRTM model was also re-sampled to 1ʹ and added the sediment corrections, called the mixed SRTM-corrected sediment, for determining the Moho and LAB depth model; • Then, using simple or 2D-Gaussian low-pass filtering, the short wavelengths (i.e., < 100 km) were removed from the mixed SRTM-corrected sediment grid, as shown in Fig. 5b, c, respectively. The 2D-Gaussian low-pass filtered elevation data are much smoother, especially in mountainous regions with rapid elevation changes.
The correlations between Moho depth derived from seismic data with elevation, elevation-corrected sediment removing short wavelengths with a 2D-Gaussian or simple filter are displayed in Fig. 6a. A high correlation with elevation data after removing short-wavelength signals (R = 0.81 with 2D-Gaussian filter and R = 0.80 with simple filter), identical to between Moho depth with geoid height after removing long wavelengths, can be seen. The two filtering methods give nearly identical results in terms of correlation, but the 2D-Gaussian low-pass filtered elevation data are much smoother and will be used to compute Moho and LAB depth.
Then, the geoid height, removing long wavelengths with a 1D-Gaussian high-pass filter for low spherical harmonic coefficients (d/o 9), and elevation-corrected sediment layer removing short-wavelength components using a 2D Gaussian low-pass filter (i.e., 100 km), will be used in the section "Results" to determine Moho and LAB depth. However, in this study region, no seismic data are available on the sedimentary basins, e.g., Khorat Plateau or Gulf of Tonkin, to validate the proposed method. Therefore, a calculation is first performed on a test area where reliable seismic data allow to assess the efficiency of applying the sedimentary correction to topography in computation of the Moho depth. Consequently, the Moho and LAB depth will be determined for Vietnam using the proposed method.

Moho and LAB depth estimation for the test area
The selected test area is located in Central Eurasia where a Moho depth model was determined by Robert et al. (2017) using geoid height and elevation data derived from EGM2008 and ETOPO1. In this area, bounded by 35-50° N and 70-100° E, lies the Tarim basin, one of the largest  basins in the world. A set of relatively complete seismic data was collected to validate the model proposed by Robert et al. 2017. From these seismic data, 30 stations located on the sedimentary basins are selected to validate our approach (Fig. 7a). Similar to Robert et al. (2017), the geoid height derived from EGM2008 after removing signals to d/o 11 and elevation data from ETOPO1 after filtering wavelength shorter than 100 km are used to calculate the Moho depth. The sedimentary thickness derived from the CRUST1.0 model is used to calculate the effect on the topography. The procedure described in the section "Topography and sedimentary effect" is used to correct for the effect of sedimentary layer and the filtering of elevation data, and the parameters are listed in Table 2.
The thickness of sediment layer, topography, topography after short-wavelength filtering, and topography-corrected sedimentary effect after short-wavelength filtering are shown in Fig. 7a-d, respectively. The correlations between Moho depth and elevation after removing short wavelengths, and elevation-corrected sediment after removing short-wavelength were computed and are displayed in Fig. 6b. A higher correlation coefficient is obtained when the sediment effects are taken into account (R = 0.64 versus R = 0.59).
We then computed the Moho depth for the test region using the parameters given in Table 3. The absolute discrepancies beneath seismic stations, between our different estimations of Moho depth (using elevationcorrected sediment after removing short wavelengths and elevation after removing short wavelengths only), as well as those obtained from Robert et al. (2017) are shown in Fig. 8. The average discrepancies are 4.2, 5.0, and 6.3 km, respectively. The comparison indicates that  there is significant improvement in the accuracy of the Moho depth when the elevation correction for the sedimentary effect is used, and it will be applied in the determination of the lithospheric structure for Vietnam in the section "Moho and LAB depth estimation for Vietnam".

Moho and LAB depth estimation for Vietnam
The densities of the crust ( ρ c ) and lithospheric mantle ( ρ m ), and reference Moho depth (z cREF ) play an important role in the computation model of Moho and LAB using the geoid height and elevation data. The crustal density is determined by assuming a linear increase with depth between constant values at the upper ( ρ T c = 2670) and lower ( ρ B c = 2820) crust, corresponding to an average crustal density of 2745 kg/m 3 . The lithospheric mantle density is considered to vary with temperature. To calculate the lithospheric mantle density, the temperatures at the surface and at the base of the lithosphere are fixed. For the reference Moho depth value, an average depth of 30.6 km was estimated using the 50 seismic receiver functions in the study region. One has to take into account that these stations are mainly in the north of Vietnam. There are no stations in the Mekong Delta where the thickness is expected to be much thinner than in the north. The Moho depth was computed with z cREF ranging from 28 to 31 km in steps of 0.5 km. The Moho depths were then compared to those derived from receiver function analysis. Finally, the best Moho depth model was obtained for z cREF = 28.5 km. The remaining parameters are used as listed in Table 3. After filtering of geoid height and elevation data, the maps of Moho and LAB depth, called Moho_GEOID and LAB_GEOID, respectively, were obtained. The results are listed in Table 4 and shown in Fig. 9.
The Moho depths vary from 15.3 to 37.8 km with mean and STD of 28.0 and 4.0 km, respectively. The lithospheric depths vary from 82.3 to 144.7 km with mean and STD of 123.1 and 12.3 km, respectively. Small Moho depths are generally seen over sea ranging from 15 to 25 km, increasing with topographic elevation in mainland Fig. 8 Absolute differences of Moho depth beneath seismic stations and those from the new maps: a using elevation-corrected sediment after removing short wavelengths (100 km), b using elevation after removing short wavelength (100 km), and c Moho depth from Robert et al. (2017) areas with maxima in the northwest mountainous region, ranging from 32 to 35 km (up to 37 km in the south of China). The topography-corrected sediment layer after removing short-wavelength components (Fig. 5c) is in good agreement with the trend in estimated Moho depth model. High elevations (up to 3100 m) in the northwest region are linked to a thick crust corresponding to the high loading there. Moreover, the geological and geodynamic settings of the area are fairly complex there. On the mainland, local minima are observed of about 27 km in the Mekong Delta (southern Vietnam). Similar to Moho depth, LAB depth also tends to be small over sea and larger in the north of Vietnam.

Evaluation and validation
The assessment of our models of Moho and LAB depth is carried out here in two steps: first, from a comparison with available global or regional models (statistics summarized in Table 4) and then from a confrontation with the results derived from seismic data.

Comparison with global and regional models
The CRUST1.0 (Laske et al. 2013) global model, based on seismic data, has a resolution of 1 arc-degree. The model includes various types of information for the Earth's crust, e.g., density, sedimentary thickness, as well as Moho depth values. Figure 10a shows the Moho depth of CRUST1.0 for the study region, which varies from 20.0 to 40.9 km. The comparison with our geoid-derived model (Moho_GEOID) reveals large differences, ranging from − 14.4 to 6.1 km with mean bias and STD of − 4.6 and 3.1 km, respectively (Table 4). This large bias is probably due to inaccuracies of the CRUST1.0 model for this area, as seismic data are very sparse. Moreover, it should be stressed that both models have different resolutions (1′ and 1° for our geoid-derived model and CRUST1.0, respectively).
The global gravity Moho model (GEMMA; Reguzzoni and Sampietro 2015) was computed by inverting the global grid of second radial derivatives of the gravitational potential observed by the ESA mission GOCE and by considering the density model from (Christensen and Mooney 1995) for the continental crust and from (Carlson and Raskin 1984) for the oceanic crust. GEMMA is available with a resolution of 0.5 0 . Figure 10b shows the Moho depth derived from GEMMA for the study region with a range of variation between 17.6 and 54.4 km ( Table 4). Figure 9a, b shows that variable trends of Moho_GEOID and the GEMMA model are close. This is reasonable, because both models are based on the isostatic gravity hypothesis. Nevertheless, their differences range from − 17.3 to 3.8 km with mean bias and STD of − 5.8 and 4.1 km, respectively. A possible reason for the bias is that the average density contrast of the Moho interface used in GEMMA is significantly different from the one adopted in this study. We can see also small discrepancies over sea and flat areas with differences less than 5 km, and the largest discrepancies are over rugged areas. This might be due to the resolution of GOCE of about 100 km allowing to resolve long wavelengths of Moho undulations only. A Moho depth model was previously calculated for north Vietnam (above 16 0 N), by Nguyen et al. (2018) using Bouguer gravity anomalies derived from a map at 1:500,000 in 2011 for the onshore part and satellite altimetry inferred gravity [UCSD V23.1; (Sandwell et al. 2014)] model for the offshore part. The comparison of Moho depths in North Vietnam derived from this study and from Nguyen et al. (2018) show differences ranged from − 9.4 to 3.3 km with mean bias and STD of − 1.2 and 1.8 km, respectively (Table 4). The differences are mainly within ± 1 km and confirm the good agreement between both models. The larger differences are over sea and surrounding areas. This might be attributed to different sediment corrections and gravity data used in the calculations. Note also that the source of the gravity data used for the surrounding areas is not disclosed in Nguyen et al. (2018).
To better evaluate the impact of different approaches of Moho modeling using gravity data, we have also estimated the Moho depth from the inversion of our Bouguer anomaly map following the classical method proposed by Parker-Oldenburg (1974). The 3DINVER.M (Gómez-Ortiz and Agarwal 2005) program was used for this purpose. For a better comparison with our Moho_GEOID model, a consistent mean density of crust and mantle, and averaged Moho depth had to be determined. Using Eq. (4) and our depth models of Moho and LAB, we have first determined the variation of the averaged lithospheric mantle density for the study area (see Fig. 11a). From these variations (3247 to 3263 kg/m 3 ), a mean value of 3252 kg/m 3 was determined. A mean crustal density of 2745 kg/m 3 and an average Moho depth of 28.5 km were finally used to compute the Moho depth surface through the Parker-Oldenburg's gravity inversion method. The Fig. 11 a Average value of the lithospheric mantle density and b differences between Moho_GEOID and Moho_GRAVITY models resulting model (Fig. 10c), called Moho_GRAVITY, varies from 12.2 to 38.5 km with mean and STD of 28.5 and 4.3 km, respectively. As expected, this model is in good agreement (differences mostly within ± 1 km) with the trend in Moho_GEOID model (see Table 4 and Fig. 11b). The largest differences (ranged from − 4.3 to 7.5 km) observed in the northern and southern edges of the study region are due to edge effects in the inversion procedure. A notable difference is visible between the sedimentary basin in the Gulf of Tonkin and and a small sediment area with thickness up to 3 km (Fig. 5a). This difference is related to the sedimentary corrections and low-pass filtering procedures in two models. Then, the consistency of the two approaches Fullea et al. (2007) and Parker-Oldenburg (1974) can be validated here for Vietnam.

Validation from seismic data
Tele-seismic results obtained from a total of 50 seismic stations were also compiled to validate our model. They include 24 stations from , 2 stations from Li et al. (2008), 1 station from Noisagool et al. (2014), 13 stations from Yu et al. (2017), and 10 stations of MRF, Vietnam from Su et al. (2018). Four of the stations are in Thailand, 2 stations in southern China, and 1 station on Hainan Island, China. To evaluate the possible improvement of the models developed in this study, these seismic data were also used to compare with existing Moho depth models in the study region. Moho depths derived from the receiver function are calculated beneath the topographic elevation, instead of referenced to sea level, so the elevation of seismic stations had to be subtracted before comparing the results.
Statistics of these comparisons are listed in Table 5 and shown in Fig. 12. Generally, our model reveals a significant improvement in terms of mean and STD when compared with CRUST1.0 and GEMMA. There is almost no bias (0.2 km) between Moho depths derived from the Moho_GEOID model and all 50 tele-seismic stations, and STD is only 2.2 km; with CRUST1.0, mean bias and STD are − 3.9 and 3.2 km, respectively, and with GEMMA − 5.9 and 3.3 km, respectively. From Fig. 12a, c, a significant improvement can be seen in the northern part, especially in the mountainous northwest region. This is thanks to the land gravity data included in our new estimations, whereas GEMMA was constructed with lower resolution GOCE data only, making it less accurate in mountainous areas. Figure 12c shows that the GEMMA model has good quality in plain areas (northeastern Vietnam), it is even better than CRUST1.0 there. Also, thanks to the ground gravity data, the model has higher accuracy in Vietnam than for surrounding regions, where fill-in data have been used. This is likely to be demonstrated in Fig. 12a, which shows that the three seismic stations on the Khorat Plateau in Thailand have larger differences than those in Vietnam. The elevation of this plateau is just 100-250 m corresponding to Bouguer gravity anomalies about − 60 mGal (Fig. 2c). The Moho depth, derived from seismic data, is up to 40 km, even if it lies in the same Indochina terrane as the northwest of Vietnam. The Bouguer anomalies in the Khorat are not sufficiently negative to explain the thick crust. According to Noisagool et al. (2014), a simple possible explanation for the difference in Moho depth within the same terrane is the fact that the Khorat Plateau is a forefront for many past collision activities. One more reason to explain the large differences between our model and results from seismic data is that the thickness of the sedimentary layer is up to 8 km for the Khorat Plateau. Another sedimentary density may be needed in calculating corrections for topography as well as Bouguer gravity anomalies. However, there are not enough seismic data on this plateau to determine an optimal density. This is another reason for the large difference between our model and results from seismic data there.
The improvement with respect to CRUST1.0 in the northeast plain region (Red River Delta) is clearly seen in Fig. 12a, b. This may be due to seismic data of this area not being assimilated in CRUST1.0. In contrast, in the northwestern part and central regions of Vietnam, the CRUST1.0 model has good quality and differences are only a few km. For the purpose of evaluating the accuracy of the new model in Vietnam, a set of 43 stations was selected to validate the Moho_GEOID model. The results are listed in Table 5. A STD of 1.5 km is found for the territory of Vietnam only. A set of 40 seismic stations in northern part (> 16 0 in latitude) was also used to compare the new model with Nguyen et al. (2018). The results are listed in Table 5. The comparison indicates that our estimation and Nguyen et al. (2018) are consistent in terms of mean value and STD in the northern part of the study area with slight improvement of our model: STD of 1.2  For the assessment of the LAB depth model, the recent global tomography model SL2013sv at 0.5° resolution (Schaeffer and Lebedev 2013) is used. S-velocity derived from this model corresponding to lithospheric depth is shown in Fig. 13. Generally, high S-velocity is located in south China and the Khorat Plateau corresponding to thick lithosphere. Major differences between our estimation and S-wave velocity-derived model are observed in northern Vietnam where the geological setting is fairly complicated. Moreover, the different approaches used for the lithospheric modeling can also lead to this difference. Unlike the Moho depth, the validation of the LAB depth model cannot be more investigated due to the poor knowledge of the LAB depth in the region. This is also a more general problem of uncertainties in determining LAB depth in other regions of the world. The global lithospheric model is still poorly constrained, as stated by Artemieva (2011): "none of the geophysical or petrologic techniques can resolve the base of the lithosphere (regardless of its definition) with a resolution better than ~ 50 km".

Discussion
The results obtained from the present study can be highlighted in the context of the geodynamic evolution of the region. The lithospheric structure of the study area is not only linked to the collision between the Indochina and South China blocks but also to the extrusion process caused by the collision between the Indian Plate and Eurasian Plate, which began at about 50 Ma (Huchon et al. 1994). The rapid increase of the Moho depth in the northwest region is considered due to the southeastern extension of the eastern Tibetan plateau along the eastern Himalayan syntaxis (Dinh et al. 2018). The thinner crust in the Hanoi plain and its eastern coastal plain region is suggestive of a recent rifting process of the opening of the East Sea (Dinh 2010). This is called the Red River rift, and it is located between the LRF and RRF (Mazur et al. 2012). Hence, an uplifted structure characterizes the Moho surface in this section where the Moho depth is about 28 km. The obvious differences in the trends of structures in the northeast and southwest sides of RRFZ are signatures of crustal scale deformation along this fault system in northern Vietnam. However, no remarkable differences in LAB depth were obtained from our results. This suggests that the RRFZ developed within the crust and remained a crustal fault. These differences in Moho depth also confirm that the RRFZ may be considered to be the boundary between two continental blocks: South China and Indochina blocks. The NW-SE trend in the northern part corresponds to the extension direction of major active fault systems, as shown in Fig. 1b. The Moho depth models agree that the RRFS extended southward in the sea of Vietnam and joined the 109° meridian fault zone. The 109° meridian fault, originating in the south Hainan island area, is the major fault line trends south, passing through the whole of the East Vietnam continental shelf, then the west coast of Kalimantan island until it reaches the Sunda bay of Indonesia (Fig. 1a). This confirms the important role this fault plays in the tectonic movement across the sea as well as the southeastern continental shelf of Vietnam (Hong Nguyen et al. 2012). Another hypothesis by Tapponnier et al. (1986) suggested that it is affected by the southeastward extrusion and large counter-clockwise motions of the Indochina block along the RRF and Wang Chao fault which in turn were reactivated by the collision of the Indian subcontinent with Eurasia. Obviously, deeper geological and geophysical investigations would be required to refine the complex geodynamic evolution of this region.

Conclusion
We have determined new models of Moho and LAB depth for Vietnam and its surrounding areas at 1′ resolution based on the high-accuracy and -resolution GEOID_LSC_C recently developed (Vu et al. 2019) and mixed SRTM models. Sedimentary thickness data derived from CRUST1.0 model were used to correct the elevation data. The topography-corrected sedimentary basin model is employed to refine the inverse model presented by Fullea et al. (2007). The results obtained for a test area demonstrated that a better consistency with seismic data is reached when elevation corrected for the sedimentary effect is used. The average discrepancies between the Moho depth derived from seismic data with those inferred from this study using elevation corrected for sediment and elevation data only are 4.2 and 5.0 km, respectively. This elevation-corrected model was used to provide new insights in the lithospheric structure beneath Vietnam. It was determined using a linear density increase with depth in the crust and a thermal analysis for the lithospheric mantle density under the assumption of local isostasy. High-pass filtering with a 1D-Gaussian function was used to remove the long wavelengths in the geoid height using GOCE DIR-R5 up to d/o = 9. 2D-Gaussian low-pass filtering was used to remove the short wavelengths (i.e., < 100 km) from the mixed SRTM-corrected sedimentary thickness. These filtered data were used to determine the Moho and LAB depth models.
The Moho depths derived from our model vary across the study region from 15.3 to 37.8 km with mean and STD of 28.0 and 4.0 km, respectively. The lithospheric depths vary from 82.3 to 144.7 km with mean and STD of 123.1 and 12.3 km, respectively. These results were compared to the seismic data as well as to global and local models of the study region. A set of 50 seismic stations distributed over Vietnam was used to validate locally the proposed model of Moho depth. This comparison indicates an agreement within 1.5 km in STD for the entire country. This new Moho depth model thus represents a significant improvement over the global models CRUST1.0 and GEMMA, which have STDs of 3.2 and 3.3 km, respectively, when compared to the same seismic data. Regarding the geoid-derived LAB depth model also inferred from this study, it is consistent with a recent global model of S-wave velocity model from this region. However, the current uncertainties on the structure of the lithosphere beneath Vietnam and its surrounding areas do not allow a more thorough evaluation. The differences in Moho depth visible on the northeast and southwest sides of RRFZ confirmed that the RRFZ may be considered the boundary between the South China and Indochina blocks. However, no remarkable differences in LAB depth were revealed, which suggests that the RRFZ developed within the crust and remained a crustal fault.