1 Introduction

The study of Moho structures can provide crucial information for applications in geology, geophysics, and geodesy. For example, such information is useful in geodynamic studies of the Earth’s interior, in unraveling the gravitational signals of anomalous subsurface distributions and in understanding geoid anomalies (e.g., Abrehdary 2016). Models of the Moho depth (MD) are resolved either by seismic, isostatic compensation or gravimetric-isostatic techniques. Typically, seismic data are generated by the travel-time of seismic waves from ground-based man-made explosions reflected at the denser Moho boundary and returned to the Earth’s surface. As this method is expensive, particularly at sea, such data become sparse, inhomogeneous and/or interpolated, which leads to low resolution and accuracy in the process. Hence, over large parts of the world, particularly at high seas, with little or no coverage of seismic data, methods based on gravimetric-isostatic techniques, and in combinations with seismic data (combined approaches), offer attractive alternatives. In this context, the global and mostly homogeneous data available free of charge from the dedicated satellite gravity missions GRACE (Tapley et al. 2004) and GOCE (Drinkwater et al. 2003; Floberghagen et al. 2011) are of particular interest.

Different attempts to estimate a global crustal model have appeared during the years with different and often unclear accuracies. A pioneer model, named CRUST5.1, was determined at a resolution of 5° × 5° by Mooney et al. (1998) based on seismic refraction data, succeeded by CRUST2.0 (Bassin 2000), which was compiled with a resolution of 2° × 2°. Another model was created by Zhou et al. (2006)-based constraining Rayleigh and Love wave data in the crust and upper mantel, while Meier et al. (2007) directly used surface wave data in their global MD model. Today, CRUST1.0 (Laske et al. 2013) is the most frequently cited and applied model, compiled from active seismic source and receiver function data. A recent global model of Szwillus et al. (2019), here named CRUST19, uses geostatistical Kriging technique (closely related with the method of collocation) for interpolation of discrete crustal depths derived from seismic data to a global model.

The isostatic compensation of the gravity field originates with Airy (1855) and Pratt (1855, 1859) [local isostatic compensation] and Vening Meinesz (1931) [regional compensation] and, finally, Moritz (1990, Section 8.3.2) [global compensation], whose technique was solved to 2nd order by Sjöberg (2009) under the name Vening Meinesz-Moritz (VMM) method. The gravimetric-isostatic compensation is dependent on both variables MD and Moho Density Contrast (MDC), but the methods above are only solved for the MD by assuming a fixed MDC. Sjöberg and Bagherbandi (2011) modified the method to solve for the MDC alone or in combination with the MD using both seismic and gravity data. Later, Tenzer and Bagherbandi (2012) demonstrated that replacing the gravity anomaly by gravity disturbance improves the VMM solution, and Sjöberg (2013) confirmed that the gravity anomaly and kernel function of the original VMM solution are inconsistent, so that either the gravity anomaly should be replaced by the gravity disturbance, or the kernel function needs a slight change.

The uncertainties in the gravimetric-isostatic MD models are typically attributed to stochastic and systematic errors (NIEs) in the global gravity, topographic and ice and sediment model data sets, but also to systematic biases (so-called Non-Isostatic Effects; NIEs) due to unmodeled mantle (and crustal) density structures and geophysical processes. As the observed gravity data include signals generated by sources in the Earth’s interior that may cause significant biases in the computational results, such NIEs must be reduced below noise level. The method used here to remove those effects was introduced by Bagherbandi and Sjöberg (2012).

Errors in the seismic crustal models arise from several factors such as the used survey method, the spatial resolution of the survey (for example the spacing of the shot points and the recording stations), and the analytical techniques utilized to process the data (Christensen and Mooney 1995; Chulick et al. 2013). A major problem in validating the accuracies in available seismic crustal models of MD and MDC is that their qualities, which vary from place to place, are typically not specified. The uncertainties are different in various seismic methods and can differ even for the same method in different surveys. The largest uncertainties of the MDs are strongly attributed to inaccuracies of crustal models currently available (Grad and Tiira 2012). Čadek and Martinec (1991) estimated MD uncertainties in their global Moho model of the order 20% (5 km) for the oceanic crust and of the order 10% (3 km) for the continental crust. The results of more recent seismic and gravity studies, however, demonstrate that these error estimates are too optimistic. Grad et al. (2009), for example, revealed that the uncertainties in the estimated MDs based on the seismic data under the European continent are about 10 km with an average error of more than 4 km. Much larger Moho uncertainties are expected over large portions of the world, particularly at sea, where the seismic data are sparse.

In the present article, some new features are used to determine a global MD model called MOHV21 from five available models, namely (a) as the standard errors (STEs) of the models and correlations among them are frequently missing, these statistics are estimated (cf. Sjöberg and Abrehdary 2021) as a preparation to (b) a proper adjustment of the individual models in a weighted least squares approach, c) including modeling the STEs of MOHV21. The following MD models are used in the merging process: MDN07 (Meier et al. 2007), CRUST1.0 (Laske et al. 2013), GEMMA1.0 (Reguzzoni and Sampietro 2015), KTH15C (Abrehdary et al. 2017) and CRUST19 (Szwillus et al. 2019). The new model is determined at a resolution of 1° × 1°.

Admittedly, as the input data were not specified for continental and oceanic crusts, we could not with reasonable efforts specify the resulting MDs accordingly, but we use the terms land and ocean region MDs.

2 Methodology

We aim at combining five global Moho models, pixel by pixel, by a weighted least squares procedure. Only very few Moho models have been determined along with detailed STEs, and many of them may need to be independently checked. Another problem is that many models partly use the same data, which means that they are correlated. Hence, such correlations should be determined and considered in the combination. Once these problems are solved, one may start from the general observation equation for each pixel:

$$ {\mathbf{e}}x = {\mathbf{L}} - {{\boldsymbol{\upvarepsilon}}};\quad E\left\{ {{\boldsymbol{\upvarepsilon}}\, {\boldsymbol{\upvarepsilon }}}^{{\text{T}}} \right\} = {\mathbf{Q}} = \sigma_{0}^{2} {\mathbf{P}}^{ - 1} $$
(1)

where eT = (1, 1, …, 1), x is the unknown MD of the combination, L and ε are the vectors of observations and their errors, Q is the error variance–covariance matrix, P is the weight matrix and \(\sigma_{0}^{2}\) is the variance of unit weight. Superscript T denotes the transpose of a vector or matrix.

Assuming that the errors are stochastic with vanishing expectations, the least squares solution for x becomes

$$ \hat{x} = {\mathbf{e}}^{{\text{T}}} {\mathbf{PL}}/\left( {{\mathbf{e}}^{{\text{T}}} {\mathbf{Pe}}} \right) $$
(2)

with the error estimate

$$ {\hat{\varvec{\varepsilon }}}_{x} = {\mathbf{e}}^{{\text{T}}} {\mathbf{P}}{\hat{\varvec{\varepsilon }}}/\left( {{\mathbf{e}}^{{\text{T}}} {\mathbf{Pe}}} \right) $$
(3)

and standard error (STE)

$$ \sigma_{x} = \sigma_{0} \left( {{\mathbf{e}}^{{\text{T}}} {\mathbf{Pe}}} \right)^{ - 1/2} . $$
(4)

The variance of unit weight (\(\sigma_{0}^{2}\)) can be unbiasedly estimated by

$$ s^{2} = \left( {{\mathbf{L}} - {\mathbf{e}}\hat{x}} \right)^{{\text{T}}} {\mathbf{PL}}/\left( {n - 1} \right). $$
(5)

where n is the number of observations in the adjustment.

To be able to determine the weight matrix P, we consider the following assumptions and solutions (see also Sjöberg and Abrehdary 2021). Firstly, we assume that y and z are two stochastic variables (read observed MDs) with the same expectation, i.e.,

$$ E\left\{ y \right\} = E\left\{ z \right\}. $$
(6)

Then, the variance of z follows from that of y as

$$ \sigma_{z}^{2} = \sigma_{y}^{2} + E\left\{ {z^{2} - y^{2} } \right\}. $$
(7)

Proof

Let \(y = u + \varepsilon_{y}\) and \(z = u + \varepsilon_{z}\), where \(u\) is the expected value of y and z. \(\varepsilon_{y}\) and \(\varepsilon_{z}\) are the stochastic errors with vanishing expectations. Then, Eq. (7) is verified by inserting y and z into Eq. (7).

Secondly, the covariance between z and y is given by

$$ q_{zy} = \left( {\sigma_{z}^{2} - G} \right) $$
(8a)

where

$$ G = E\left\{ {z\left( {z - y} \right)} \right\}. $$
(8b)

Proof

Upon substituting y and z in Eq. (8b) as in the proof of Eq. (7), it follows that

$$ G = \sigma_{z}^{2} - q_{zy} , $$
(9)

which verifies Eq. (8a).

In the practical applications below of Eqs. (7) and (8b) for the MD at a specific pixel, the expectation operator will be replaced by a weighted mean value over a local set of pixels centred at the pixel of interest. More precisely, only 5 pixels are used in the average, where the central pixel is given the weight 1/3, and each of the closest north, south, east, and west pixels is given weight 1/6. Also, to avoid errors from biases in z and y, we will remove the (weighted) mean values of these quantities in the region around the computation point are removed prior to computations. (These refinements were proposed but not utilized in Sjöberg and Abrehdary 2021). Finally, note that, we did not use latitude weighting in the gridding, which to some extent biases the determined mean values toward data from high latitudes.

3 Description of the Moho models in the adjustment

Here follows a short description of the five Moho models used in the study.

3.1 MDN07

Meier et al. (2007) presented the MDN07 model, which also comprise STEs, was produced by phase and group velocities of Rayleigh and Love waves at a resolution of 2° × 2°. The model, calculated by a nonlinearized neural network method, is limited by the lateral resolution of the input phase and group velocity maps, ranging from 500 to 1000 km.

3.2 CRUST1.0

Laske et al. 2013 released the CRUST1.0 Moho model, which is tailored from active source seismic and receiver function surveys with a resolution of 1° × 1°. It is composed of 35 key crustal features, incorporating its thickness, density, and compressional wave velocities (Vp) and shear waves (Vs) for eight layers of ice, water, upper, middle, and lower sediments, upper, middle, and lower crust. The Vp values are estimated by field measurements, whereas Vs and density are computed using empirical Vp–Vs and Vp-density relationships, respectively. As large parts of Africa, South America, Antarctica, and Greenland lack field measurements, the seismic velocity structure of the crust is extrapolated from the average crustal structure for areas with similar crustal age and tectonic setting. The topography and bathymetry are achieved from ETOPO1 (Amante and Eakins 2009), while the sediment basins are covered by the sediment model of Laske (1997) with some near-coastal updates. The CRUST1.0 marine sediment structures are compiled by Divins (2003).

3.3 The GEMMA1.0

Reguzzoni and Sampietro (2015) presented the Earth crustal model GEMMA1.0 by applying an inversion algorithm to satellite GOCE-derived gravity data and some seismic data at a resolution of 0.5° × 0.5°. The model incorporates different additional external information such as topography, bathymetry, and ice sheet models from ETOPO1, a recent 1° × 1° sediment global model and some prior hypotheses on crustal density. The model also considers lateral density variations in the upper mantle.

3.4 The KTH15C

Abrehdary et al. (2017) produced the KTH15C model using a preliminary combined VMM solution for the MD at a resolution of 1° × 1° by applying the data sets of GOCO05 global gravitational field (Mayer-Güerr 2015), the Earth14 solid Earth topographic model (Hirt and Rexer 2015) and CRUST1.0 crustal model in a least squares method.

3.5 The CRUST19

Szwillus et al. (2019) determined a global Moho model (here named CRUST19) and velocity structure based on the analysis of seismic data at a resolution of 1° × 1° utilizing a nonstationary kriging algorithm for data interpolation of an MD model with STEs. The discrete data were obtained from the U.S. Geological Survey Global Seismic Catalog database. The only additional data are a model of the age of the ocean floor, which is needed to separate continental and oceanic domains. As CRUST19 does not contain any data sets of gravity anomalies, the result is predicted to be vague in areas where the data coverage is missed or sparse. However, CRUST19 claims that this is not drawback of the method but consider it as an illustration of the limitations of the current database. The method allows lateral variation in covariance parameters expected on a global scale (Risser and Calder 2015).

4 Variance–covariance models among the Moho models

Equations (7) and (8a) are applied (pixel by pixel) to find the elements of the Q matrix for the adjustment. We start from the MD and variance data of the KTH15C model to determine its variances of and correlations with the other four models (The reason for starting from KTH15C is that it generally has the smallest STEs among all included models). Although variance estimates are available also for some other models, we believe that it is more appropriate to recompute also these parameters to provide better weights for the joint adjustment.

Some statistics of the estimated MDs of the five models, as well as their STEs computed by Eq. (7), are shown in Table 1. Generally, the smallest STEs are given by KTH15C, while those of GEMMA1.0 and CRUST19 are larger for other models on land and ocean, respectively.

Table 1 Statistics of the MDs given by MDN07, CRUST1.0, GEMMA1.0, KTH15C and CRUST19 to a resolution of 1° × 1°
Full size table

In Table 2, the correlations among the five models determined by Eqs. (8a) and (8b) are demonstrated. Both rows and columns represent in order MDN07, CRUST1.0, GEMMA1.0, KTH15C, and CRUST19.

Table 2 The four matrices represent the ID of maximum positive, average, and maximum negative correlations
Full size table

The first matrix in Table 2 is given for identification of the correlations that follow numerically in the three matrices to the right. The three matrices represent the maximum positive, average, and maximum negative correlations. Both rows and columns are ordered as the MD models in Table 1. The table shows that the major maximum positive correlations are obtained between CRUST1.0 and KTH15C, slightly followed by CRUST1.0 versus CRUST19. This result is reasonable, as both KTH15C and CRUST19 rely on CHRUST1.0. The smallest maximum and mean correlations and largest negative correlations are seen between MDN07 and all other models.

5 The adjustment and results

Here, Eqs. (2), (3) and (5) are applied pixel by pixel, and the result is summarized in Table 3.

Table 3 Statistics of model MOHV21
Full size table

Comparing Tables 1 and 3, one can see that, in general, the STEs in MOHV21 are significantly smaller than those for other models.

Figure 1a–b visualizes the new MD model MOHV21 and its STE. The largest MDs, of order 70–75 km, with STEs of order ± 6–6.5 km, are found on land in the Himalayas and along the central westcoast of S. America and in the Arequipa Massif between the Andes and the Pacific Ocean. A notable feature is the regional MD maxima in Tarim basin in NW China of 59 ± 6.5 km and in Central Finland of 57 ± 4.7 km. At sea, the MD is generally within 20 ± 2.5 km and well within 10 km along mid-ocean ridges.

Fig. 1
figure 1

a The MD model MOHV21 and b its standard error. Unit: km

Full size image

For a comparison, we compiled data also from some regional MD models in Tables 4, 5 and 6 and Fig. 2, they are compared with MOHV21. As one can see, there are large differences, ranging to ± 25 km in Africa, Antarctic, and part of S. America, while the differences are generally gentle on the parts of oceans covered by the regional models. Note that, the regional models do not include gravimetric-isostatic data, which may be the main reason for their large differences to MOHV21. Another advantage of MOHV21 should be the least squares combination of different models.

Table 4 Sources of regional MD models used in this study
Full size table
Table 5 Statistics of the differences between MDs from MOHV21 and regional seismic models
Full size table
Table 6 Statistics of the differences between MOHV21 and all regional seismic models (WORLD)
Full size table
Fig. 2
figure 2

MD differences in km of MOHV21 minus a mosaic of regional seismic MD models shown in Table 4

Full size image

6 Conclusions and final remarks

Model MOHV21 combines MD estimates from five global models, based on seismic and gravimetric-isostatic data and models, at a resolution of 1° × 1°. The combination is performed pixel by pixel by highlighting good and toning down bad data in a least squares sense. By including correlations among the data in the weighting procedure, no information is repeatedly used. The result is a MD model that statistically performs better than any of its five ingredients. Its standard error is mostly within 3 km and ranges to 6.5 km only in the central west coast of S. America and in the Himalayas, where it reaches its maximum of 75 ± 6.8 km. The shallow MDs along mid-ocean ridges are well exposed in MOHV21, and much of the outliers seen in other models are eliminated, but some special regional features deserve mentioning. For instance, in the Tarim basin in NW China and in Central Finland, the MDs reach 59 ± 6.5 km and 57 ± 4.7 km. In Antarctica, the MOHV21 MD is rather constantly around 35 km, and smaller MDs, typically 20 km or less, are seen in the oceanic areas, particularly along mid-ocean ridges, where the MD is mostly limited to 10 km or less.

Further studies are needed to explain the large differences (ranging to ± 25 km) between MOHV21 and some regional models in Africa, Antarctica, and parts of S. America. However, in this context, one important factor is that the regional models are merely generated from seismic data, while MOHV21 demonstrates the combined result of seismic and gravimetric-isostatic modeling.