Lithospheric elastic thickness estimates in central Eurasia

We estimate the elastic thickness of a continental lithosphere by using two approaches that combine the Vening Meinesz-Moritz (VMM) regional isostatic principle with isostatic flexure models formulated based on solving flexural differential equations for a thin elastic shell with and without considering a shell curvature. To model the response of the lithosphere on a load more realistically, we also consider lithospheric density heterogeneities. Resulting expressions describe a functional relation between gravity field quantities and mechanical properties of the lithosphere, namely Young’s modulus and Poisson’s ratio that are computed from seismic velocity models in prior of estimating the lithospheric elastic thickness. Our numerical study in central Eurasia reveals that both results have a similar spatial pattern, despite exhibiting also some large localized differences due to disregarding the shell curvature. Results show that cratonic formations of North China and Tarim Cratons, Turan Platform as well as parts of Siberian Craton are characterized by the maximum lithospheric elastic thickness. Indian Craton, on the other hand, is not clearly manifested. Minima of the elastic thickness typically correspond with locations of active continental tectonic margins, major orogens (Tibet, Himalaya and parts of Central Asian Orogenic Belt) and an extended continental crust. These findings generally support the hypothesis that tectonically active zones and orogens have a relatively small lithospheric strength, resulting in a significant respond of the lithosphere on various tectonic loads, compared to a large lithospheric strength of cratonic formations. Article history: Received 26 February 2018 Revised 14 June 2018 Accepted 28 September 2018


IntroductIon
The (effective) elastic thickness of the lithosphere defines its integrated strength in response to various tectonic loads, while depending also on the lithospheric density structure and its rheological properties.According to theoretical models, large values of the continental lithospheric elastic thickness apply typically to old cratonic formations, while the elastic thickness of active tectonic margins and young orogens is relatively small.
Various methods were developed and applied to estimate the lithospheric elastic thickness.The Fourier transform (e.g., Watts 2001, p. 195) is probably the most commonly used method for the elastic thickness estimation; see studies, for instance, by Calmant et al. (1990) for the oceanic lithosphere, Filmer et al. (1993) for Marquesas and Society Islands, Audet and Mareschal (2004) for Canada, or Gómez-Ortiz et al. (2005) for Iberian Peninsula.Turcotte et al. (1981) presented the approach that combines Jeffrey's (1976) method for a Moho depth recovery with a flexural theory.Forsyth (1985) proposed the coherence method based on comparing gravity and topographic models.McGovern et al. (2002) applied the admittance analysis.Ojeda and Whitman (2002) used the coherence analysis to estimate the elastic thickness in the northern part of South America and Swain and Kirby (2003a) applied the same approach for Australia.Later, Swain and Kirby (2006) repeated the study for Australia, but applying the method developed by Forsyth (1985).McKenzie (2003) modelled the elastic thickness, while considering the effect of internal loads.Jordan and Watts (2005) used both, the forward and Terr.Atmos.Ocean. Sci., Vol. 30, No. 1, 73-84, February 2019 non-spectral inverse gravity modelling techniques to determine the elastic thickness along Indo-Eurasian continental collision zone.Tassara (2005) applied the flexural analysis along Andean margin and Tassara et al. (2007) used the wavelet form of a classical spectral isostatic analysis for South American and surrounding tectonic plates.Pérez-Gussinyé et al. (2004, 2007, 2009) conducted similar studies for Fennoscandia, South America and Africa.Galán and Casallas (2010) applied the admittance analysis to estimate the elastic thickness of Colombian Andes.McKenzie (2010) compared the coherence and admittance methods and demonstrated that their results differ significantly.Tesauro et al. (2013) presented the global model of the lithospheric elastic thickness, while considering variations of Young's modulus within the lithosphere, and Tesauro et al. (2018) took into consideration also temperature, composition and strain rates of the lithosphere.
Theoretical foundations describing a flexural response of the lithosphere on the loads were given by Vening Meinesz (1931).He assumed that the lithosphere is an elastic shell and formulated a regional isostatic principle based on a flexural model.Later, Moritz (1990) used this isostatic principle to define the method for a gravimetric Moho recovery.Sjöberg (2009) developed this method further and called it the Vening Meinesz-Moritz (VMM) theory of isostasy.The original idea of Vening Meinesz (1931), however, was based on a loading theory.In this way, mechanical properties of the lithosphere (such as rigidity, elastic thickness, Young's modulus, and Poisson's ratio) are taken into the consideration instead of gravity information.The VMM and flexural models thus describe two different methods for computing the Moho depth, but their solutions should theoretically be similar if the lithospheric shell is elastic, while an isostatic equilibrium exits in all regions globally.Eshagh (2016a) confirmed this numerically at the study area in Tibet.Eshagh (2018) facilitated this finding in deriving a numerical model for estimating the lithospheric elastic thickness based on combining the VMM regional isostatic principle and the solution to a flexural differential equation for a thin elastic layer, but without taking into the consideration a spherical curvature of the shell.In this study, we extend this definition by taking into consideration the shell-curvature term in deriving the corresponding solution to a flexural differential equation and compare the results from both flexural models in central Eurasia.Moreover, the elastic thickness of the lithosphere depends not only on mechanical properties of the lithosphere, but also on the lithospheric density structure that was not taken into the consideration in theoretical definitions by Eshagh (2018).To address this aspect, we further adopt a generalized formulation of the VMM isostatic theory that accounts for the gravitational contribution of lithospheric density heterogeneities, and utilize this generalization in both methods for estimating the elastic thickness.

MEthod
In this section, we first briefly recapitulate the VMM and flexural isostatic theories and then combine them in order to determine the elastic thickness of a continental lithosphere based on solving flexural differential equations for a thin elastic shell with and without considering the shell curvature.

VMM model
Eshagh (2016b) formulated the VMM inverse problem of isostasy for finding a Moho depth T VMM from gravity disturbances g d based on assuming only a uniform lithospheric density.Since density variations especially within a continental lithosphere could be significant, we extend his definition by taking into the consideration also the lithospheric density heterogeneities.The VMM isostatic model is then given in the following spectral form where T 0 is the mean Moho depth, R is the Earth's mean radius, G is Newton's gravitational constant, t D is the Moho density contrast (between the crust and the uppermost mantle), Y n,m are the surface spherical harmonic functions of degree n and order m, i, and m denote the spherical colatitude and longitude respectively, n is the upper summation index of spherical harmonics, and the degree-dependent parameters n b read (cf.Eshagh 2017) The gravity disturbances g d used to compute the Moho depth in Eq. ( 1) are defined in the spectral domain by the spherical harmonic coefficients g n d of the external gravity field.The gravity disturbances are corrected for the gravitational contributions of topography and density contrasts of bathymetry (i.e., the ocean density contrast), sediments and consolidated crust.These gravity corrections are described by the corresponding spherical harmonics of topography/ bathymetry g nm TB d , sediments g nm

d
, and density heterogeneities within the consolidated (crystalline) crust g nm

d
. We note here that the gravitational contribution of continental glaciers that reaches maxima of about 300 mGal (cf.Tenzer et al. 2010) should be taken into the consideration in polar regions, while the gravitational contribution of atmosphere is completely negligible (cf.Tenzer et al. 2009a).It is worth mentioning also that the application of the bathymetricstripping gravity correction is required not only in studies of the oceanic but also continental lithospheric structure (such as a Moho depth) or processes (such as a flexural responds of the lithosphere on the topographic load) because the application of this correction modifies the gravity field globally (not only offshore).
By disregarding the first term on the right-hand side of Eq. ( 1), the expression for computing the Moho depth flexure ΔT (with respect to the mean Moho depth T 0 ) is simplified into the following form , If we further consider only a continental lithosphere, the substitution for the parameters n b from Eqs. ( 2) to ( 3) yields The VMM model in Eq. ( 4) is formulated based on assuming that the crust (corrected for anomalous density structures) is in an isostatic equilibrium.However, the original idea of Vening Meinesz was that the load of topographic masses is bending the lithosphere.Below we present this concept in detail.

Flexural Models
Here, we assume two models for a flexural isostasy with and without considering the curvature of the spherical shell.

Flexural Model Without considering shell's curvature
Let us assume that a thin elastic layer approximates the lithosphere.Such a thin layer deforms under the load.In this case, the flexural deformations are described by the following partial differential equation (see Watts 2001;Artemieva 2011, p. 554) where g is the mean gravity at sea level, and d denotes the gradient operator.
The flexural rigidity D is defined as a function of Young's modulus E, the lithospheric elastic thickness T e and Poisson's ratio v in the following form The parameters E, v, T e , and D in Eq. ( 6) describe mechanical properties of the lithosphere.We further defined the density function K by where Sed t and H Sed are the density and thickness parameters of sediment layers, and the corresponding parameters Crys t and H Crys are specified for the consolidated crustal layers.The density distribution function t within the topography and bathymetry in Eq. ( 7) reads ) where H is the topographic height on land and the (negative) bathymetric depth offshore, c t = 2670 kg m -3 is the topographic density, and w t = 1027 kg m -3 is the seawater density.We note that the bathymetric gravity correction could be evaluated more accurately by adopting a depthdependent seawater density model (Gladkikh and Tenzer 2012; see also Tenzer et al. 2011Tenzer et al. , 2012)).Eshagh (2016b) presented the flexural differential equation [from Eq. ( 5)] in the following spectral form where The expression in Eq. ( 9) defines a functional relation between the spherical harmonics (ΔT) nm and Knm that describe the Moho flexure and the loading masses respectively.According to Vening Meinesz (1931), the flexural respond of the lithosphere on the load can be expressed in terms of the Moho flexure T Flex D by rearranging the expression in Eq. ( 9) into the following form (cf. Eshagh 2016b, 2018) where the degree-dependent compensation coefficients C n describe the compensation according to the rigidity of deformed shell by the following formula (cf.Eshagh 2016b) It is noted that for D = 0, the expression in Eq. ( 11) describes a Moho depth according to Airy's (1855) theory, meaning that T e = 0.In other words, Airy's isostatic principle does not hold for a flexural deformation, because it assumes a local compensation mechanism.

Flexural Model considering shell's curvature
To describe mathematically a flexural deformation due to the load, the elastic shell theory could also be adopted (cf.Turcotte and Schubert 1982).The loading theory for an elastic lithospheric shell is described in the following form (Eshagh 2016b) By analogy with Eqs. ( 5) and ( 9), the flexural differential equation [from Eq. ( 13)] can be described in terms of the spherical harmonic coefficients of the load and the flexure as follows ( ) The solution of Eq. ( 14) for (ΔT) nm and the substitution of obtained result into the spherical harmonic expansions yields the same formula as presented in Eq. ( 11), but in this case with the degree-dependent compensation coefficients (cf.Turcotte et al. 1981) If we disregard the term related to a shell curvature, the differential equation in Eq. ( 13) becomes identical with Eq. ( 5).In this way, we could see that differences in definitions of the degree-dependent compensation coefficients according to Eqs. ( 12) and ( 15) are related only to the shellcurvature term.It is important to clarify here that inverse solutions to the flexural differential equations described in Eqs. ( 5) and ( 13) are ill-posed because they comprise gradient operators of 4 th or even 6 th order depending on whether a shell's curvature is included Eq. ( 13) or not Eq.( 5).Nevertheless, the reformulation of these flexural differential equations in the spectral domain [given in Eqs. ( 9) and ( 14)] allows us to stabilize their solutions by limiting a spherical harmonic series to a specific degree and order.

Elastic thickness
Let us now assume that both Moho solutions obtained by solving the VMM and flexural isostatic models in Eqs. ( 4) and ( 11) or ( 15) are the same (or at least similar).A direct mathematical model can then be established in order to estimate the elastic thickness parameter T e , but the important issue is that in the classical flexural isostatic theory it is assumed that the lithospheric density heterogeneities are absent (cf.Artemieva 2011, p. 556).Therefore, we generalized the functional model in Eq. ( 4) for crustal density heterogeneities, while sub-crustal lithospheric density heterogeneities are also incorporated into the model, but implicitly by assuming a variable Moho density contrast.
For ΔT VMM = ΔT Flex , we then postulate that The right-hand side of Eq. ( 16) contains the gravity disturbances g d corrected for the gravitational contributions of crustal density heterogeneities, while the left-hand side comprises the same components, but described (via C n ) in terms mechanical properties of the lithosphere including the elastic thickness T e .The expression in Eq. ( 16) thus represents a basic functional relation for deriving the values of T e that could further be rearranged into the following form ( ) The expression in Eq. ( 17) defines a mathematical model that functionally relates the lithospheric mechanical properties (including the elastic thickness T e ) with the gravity disturbances (corrected for the gravitational contributions of crustal density heterogeneities).Mechanical properties of the lithosphere (including the elastic thickness T e ) are involved in the degree-dependent compensation coefficients C n .However, by solving the flexural differential equations with and without considering the curvature of a spherical shell we obtained two different definitions of the coefficients C n that are given in Eqs. ( 12) and ( 15).Now the question is which definition provides a result that is consistent with geophysical and geological evidences.We investigated this aspect numerically at the study area covering central Eurasia.

study ArEA And dAtA AcquIsItIon
The study area of central Eurasia (Fig. 1) is bounded by the parallels 0 and 90 arc-deg of the northern latitudes and the meridians 40 and 140 arc-deg of the eastern longitudes.A tectonic configuration of this study area comprises parts of Indian, Eurasian, Arabian, African, and Sunda plates including Tibetan, Iranian, and Yangtze tectonic blocks.The most prominent geological feature there is an active continent-to-continent collision between Indian and Eurasian plates that is responsible for the formation of Himalaya and Tibet Orogens.Ural Orogen, on the other hand, represents the oldest known orogenic structure that is mostly deeply buried into younger sediments.Ural Orogen separates Russian Craton from Central Asian Orogenic Belt, while Baikal Rift Zone represents a geological boundary between Central Asian Orogenic Belt and Siberian Craton.Most of the southern and central parts of Indian plate are formed by Indian Craton.Large cratonic provinces in China are belong to North China, Yangtze, Sichuan, and Tarim Cratons.
Regional maps of the CRUST1.0(Laske et al. 2013) Moho parameters are shown in Fig. 2. The Moho depth in central Eurasia is typically 30 -55 km, while deepens to more than 60 km under Himalayan and Tibet Orogens.The Moho density contrast in central Eurasia varies typically from about 600 -770 kg m -3 , with maxima under Tibet.
We further computed Poisson's ratio and Young's modulus from the CRUST1.0P and S wave velocity model (using the Matlab code provided by Bevis).The results are shown in Fig. 3. Maximum values of Young's modulus and Poisson's ratio apply for Siberian Craton, and minima are irregularly distributed along the western part of Tethyan Orogenic Belt and whole Indian Craton.
To take into the consideration the lithospheric density heterogeneities, we compute gravitational contributions individually for the topography-bathymetry, sediments and consolidated crust, while the gravitational contribution of the sub-crustal lithosphere is included implicitly by assuming a variable Moho density contrast (cf.Tenzer et al. 2015).The results are shown in Fig. 4. All computations of gravity and gravity corrections are realized on a 1 × 1 arc-deg grid within the study area, with a spectral resolution complete to the spherical harmonic degree 180.The combined topographic-bathymetric gravitational contribution reaches the largest negative values (about -500 mGal) over Himalaya and Tibet, while negative values also prevail over surrounding regions in central Asia.Elsewhere, this contribution is typically positive.The gravitational contribution of sediment density contrast is mostly within a relatively small interval ±50 mGal, except for large positive values exceeding regionally even 100 mGal on both sides of Ural Orogen (i.e., Russian Craton to the west and the West Siberian Mesozoic-Cainozoic sedimentary basin to the east).The gravitational contribution of consolidated crust density variations is everywhere negative, with small negative values over eastern China and Siberia and the largest negative values over the western part of Asia.

rEsuLts
The methodology of estimating the effective elastic thickness T e from gravity field quantities and mechanical properties of the lithosphere comprises in principle three numerical steps.Firstly, we compute the gravity disturbances from the coefficients of the EIGEN-6C4 (Förste et al. 2014) global gravitational model corrected for the GRS80 (Moritz 2000) normal gravity component.We then remove from the EIGEN-6C4 gravity disturbances the gravitational contributions of crustal density heterogeneities as seen on the left-hand side of Eq. ( 17).In particular, we apply the gravity corrections due to the combined effect of topographybathymetry (Fig. 4a), sediments (Fig. 4b), and consolidated crust (Fig. 4c).The intermediate results of this procedure are shown in Fig. 5.As seen in Fig. 5a, the free-air gravity disturbances are typically distributed within a small interval of values between ±50 mGal, except for gravity highs coupled by gravity lows along Himalaya and the northern margins of Tibetan Plateau.The refined gravity disturbances corrected for the gravitational contributions of topography, bathymetry, sediments, and consolidated crust (Fig. 5d) are mostly positive with maxima over Himalaya and Tibet.It is worth mentioning that the computation of refined gravity disturbances shown in Fig. 5d differs from computing the Bouguer gravity disturbances by means of subtracting the gravitational contributions instead of adding them to the free-air gravity disturbances.Even so, both gravity data types a highly spatially correlated with the Moho geometry (cf.Tenzer et al. 2009b).
In the second step, we compute the corresponding refined gravity disturbances as a function of mechanical parameters of the lithosphere (i.e., Poisson's ratio and Young's modulus) according to the expression on the right-hand side of Eq. ( 17), while using two different definitions of the coefficients C n from Eqs. ( 12) and ( 15) obtained by solving the flexural differential equations for a thin elastic layer with and without considering a shell curvature.This computation is realised on a 1 × 1 arc-deg grid within the study area, while for each computation surface point we compute sets of the gravity disturbances for different values of the elastic thickness T e from 0 -140 km with a step of 1 km.The value of T 0 was set to be 25 km.
Finally, we adopt the principle of minimizing differences between the gravity disturbances computed according to the left-hand side and the right-hand side of Eq. ( 17) as a criterion for selecting values of the lithospheric elastic thickness T e .To decrease the level of noise from results of the lithospheric elastics thickness, we further apply a median filter with the length of 2 arc-deg, which corresponds to approximately 200 km.The results after filtering are plotted in Fig. 6, and their differences are shown in Fig. 7.
As seen in Fig. 6, both results have a similar prevailing spatial pattern, except for some large but localized differences (see Fig. 7).Overall, we could see that maxima of the lithospheric elastic thickness generally coincide with locations of cratonic formations of North China and Tarim Cratons, Turan Platform and large parts of Siberian Craton.In contrast, large values of the lithospheric thickness of Indian and Russian Cratons are missing.Minima are found under active tectonic margins, major orogens (Tibet, Himalaya) and extended crust.Central Asian Orogenic Belt is characterized by small values of the elastic thickness, except for its south-west part.Small lithospheric thickness is also observed along Baikal Rift Zone.

dIscussIon
To better understand differences between two results plotted in Fig. 7, we first analyse them from a theoretical perspective and then numerically.The elastic thickness determination is in principle an inverse problem even if we have computed it in a forward modelling manner.Moreover, we use the same type of input data to estimate it.Hence, differences between the elastic thickness estimates are only caused by disregarding the shell-curvature term in the definition of degree-dependent compensation coefficients C n given in Eq. ( 12), meaning that if we disregard the shell curvature the load is entirely compensated by bending stresses.Alternatively, when taking into the consideration a shell curvature a part of the load is compensated also by membrane stresses (cf.Turcotte et al. 1981).In other words, differences in Fig. 7 occur at locations where the load is compensated not solely by bending stresses [Eq.( 12)], but also membrane stresses [Eq.( 15)] inside the shell.Turcotte et al. (1981) mentioned also that in order for membrane stresses to support topographic loads on the Earth, the elastic thickness of the lithosphere should be about 1000 km.However, as observed in Fig. 7, differences in the estimated elastic thickness reach up to about 100 km when disregarding the shell curvature and membrane stresses, which is much less than a theoretical estimate.Moreover, the contribution of membrane stresses on the elastic thickness T e is degree dependent.This is illustrated in Fig. 8, where we investigate the degree-dependent compensation parameters C n that are computed differently according to Eqs. ( 12) and (15).For this numerical example we used the following values: the Moho density contrast t D = 600 kg m -3 , Young's modulus E = 100 GPa, and Poisson's ratio v = 0.25 and plotted the parameters C n for the elastic thickness T e of 10 and 100 km.
As seen in Fig. 8, for T e = 10 km, the contribution of membrane stresses is not significant.The differences between values of C n computed with and without considering the shell-curvature term are small.For T e = 100 km, however, these differences become considerably larger, especially at low degrees, while again negligible at higher degrees.From these findings, we can conclude that the contribution of membrane stresses for small values of the elastic thickness is not important, because the compensation parameters in both cases are more or less equivalent.However, for large values of the elastic thickness, membrane stresses influence considerably low degrees of the compensation.This numerical example clearly illustrates the significance of considering or ignoring membrane stresses in describing a mechanism of the flexural isostasy.Hence, we could say that membrane stresses in general compensate some part of the load, thus weakness the compensation mechanism.As also seen in Fig. 8, the lithospheric structures which are controlled to a large extent by compensation mechanisms (e.g., orogens) are characterized by a small elastic thickness.Moreover, the compensation mechanism weakens much faster with increasing degrees in regions characterized by a small elastic thickness (e.g., cratons).
Let us now investigate the following two dimensionless parameters, mentioned by Turcotte et al. (1981), in order to demonstrate the importance of membrane and bending stresses to support topographic loads, namely When the parameter x is large membrane stresses can fully support topographic loads without flexure, while membrane stresses are negligible if this parameter is small.The parameter v is a measure of resistance of the lithospheric shell to bending or flexure.As seen in Fig. 9a, values of x are typically very small, meaning that membrane stresses do not support significantly topographic loads, except for some large values of this parameter that are closely correlated with areas of large elastic thickness (see Fig. 6a) where membrane stresses support loads more significantly.Similarly, we could see in Fig. 9b that places with large values of Fig. 7. Differences between values of the lithospheric elastic thickness estimated based on solving the flexural differential equations for a thin elastic layer (Fig. 6a) and for a membrane stress inside the spherical shell (Fig. 6b).
v correspond to a continental lithosphere that could better resist bending stresses.In both cases, we could see that these parameters are controlled to a large extent by variations in the elastic thickness, while their dependence on Poisson's ratio and Young's modulus is much less significant.This finding has important consequences, confirming that the strength of the lithosphere is mainly controlled by its elastic thickness, while other mechanical properties of the lithosphere are less significant.

suMMAry And concLudIng rEMArks
We apply two methods for modelling the elastic thickness of the continental lithosphere and compare their nu-merical performance at the study area of central Eurasia.Both methods are formulated based on combining the VMM isostatic theory with the solutions to flexural differential equations given for a thin elastic layer with and without considering membrane stress of the spherical shell.In this way, the elastic thickness could be estimated from gravity information and mechanical properties of the lithosphere.Moreover, we incorporate into both solutions also the effect of lithospheric density structure by subtracting the gravitational contributions of crustal density heterogeneities and taking into consideration the sub-crustal lithospheric density heterogeneities by assuming a variable Moho density contrast.
Our results reveal that both solutions are similar, except for some localized differences.A finding of a minimum  elastic thickness typically along tectonically active zones and orogens agrees with a significant respond of the lithosphere to tectonic loading deformations, while a significant lithospheric strength of cratonic formations propagates into large values of the elastic thickness.

Fig. 1 .Fig. 2 .
Fig. 1.Topography of the study area.Tectonic plate boundaries are indicated by red lines.

Fig. 8 .
Fig. 8.The degree-dependent compensation parameters C n computed from Eqs. (12) and (15) for the elastic thickness T e of 10 and 100 km.