邵翠法，王世民

(中國科學(xué)院大學(xué)地球與行星科學(xué)學(xué)院 中國科學(xué)院計(jì)算地球動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室， 北京 100049)

地殼和地幔中的固體巖石具有顯著的黏彈性性質(zhì)，即在低溫和快速變化的載荷作用下變形由可逆的彈性應(yīng)變占主導(dǎo)，而在高溫和緩慢變化的載荷作用下變形以不可逆的黏性流動(dòng)為主。

彈性巖石圈的定義以及彈性巖石圈厚度的確定是地球動(dòng)力學(xué)研究的中心課題。事實(shí)上，彈性巖石圈并沒有一個(gè)截然的下界面。彈性巖石圈的一種定義方法是借助于巖石的黏彈性松弛時(shí)間(黏度與彈性模量之比)。如定義彈性巖石圈內(nèi)松弛時(shí)間大于一千萬年[1]，相當(dāng)于干橄欖石溫度低于675 ℃、濕橄欖石溫度低于475 ℃。考慮到一般以1 300 ℃等溫線定義熱巖石圈與軟流圈界面[1]，彈性巖石圈應(yīng)看作是熱巖石圈的一個(gè)組成部分。為論述方便起見，下文中一律將彈性巖石圈以下的熱巖石圈部分稱作黏性巖石圈。顯然，整個(gè)熱巖石圈都具有黏彈性性質(zhì)(雖然在彈性巖石圈中以彈性變形為主，而在黏性巖石圈中以黏性變形為主)，且彈性巖石圈向黏性巖石圈的過渡是主要由溫度連續(xù)變化導(dǎo)致的漸變而非突變。根據(jù)松弛時(shí)間定義的彈性巖石圈厚度在原理上無懈可擊，但由于地下巖石的黏度無法精確測量，并不能實(shí)際應(yīng)用于具體區(qū)域的彈性巖石圈厚度確定。

在地球動(dòng)力學(xué)研究中，彈性巖石圈厚度通常由彈性板彎曲分析反推，即通過以彈性薄板在靜巖恢復(fù)力(浮力)作用下彎曲的解析解[1]擬合實(shí)測巖石圈彎曲導(dǎo)致的地形(如島鏈、海溝的外緣隆起位置)來反演巖石圈的抗彎剛度進(jìn)而確定彈性巖石圈等效厚度。這樣得到的彈性巖石圈厚度與根據(jù)松弛時(shí)間定義的彈性巖石圈厚度是否一致并沒有經(jīng)過嚴(yán)格的檢驗(yàn)。

嚴(yán)格來說，彈性巖石圈的厚度應(yīng)當(dāng)是載荷作用時(shí)間的函數(shù)。在地震波傳播的時(shí)間尺度上，整個(gè)地殼和地幔都表現(xiàn)出彈性性質(zhì)，而在地質(zhì)時(shí)間尺度上，彈性變形只在近地表的彈性巖石圈內(nèi)占主導(dǎo)，黏性流動(dòng)則主導(dǎo)彈性巖石圈以下地幔巖石的變形(如地幔對(duì)流)。巖石圈對(duì)各種不同時(shí)間尺度載荷的響應(yīng)可以反映相應(yīng)時(shí)間尺度下的彈性巖石圈厚度。如對(duì)震后變形[2-7]、冰后回彈[8-10]、俯沖過程[11-14]和巖石圈彎曲[15-19]等問題的黏彈性模擬可以約束相應(yīng)載荷作用時(shí)間下的彈性巖石圈厚度。

描述巖石圈彎曲的物理模型經(jīng)歷了從簡單到復(fù)雜的過程。單層彈性薄板加無黏流體地基的模型[20-28]是描述巖石圈彎曲變形最簡單且應(yīng)用最廣泛的模型。通過對(duì)比模型預(yù)測和火山、沉積物等載荷作用下地表形態(tài)的實(shí)際觀測數(shù)據(jù)，可以研究彈性巖石圈厚度與板塊年齡及載荷作用時(shí)間的依賴關(guān)系[29-31]。基于這種模型，Watts等[32]曾通過海底地形反演計(jì)算全球范圍內(nèi)9 758處彈性巖石圈的厚度。

為克服彈性板彎曲模型結(jié)果與時(shí)間無關(guān)的缺陷，之后又發(fā)展出了均勻黏度黏彈性薄板模型[33-36]。 Walcott[33]采用均勻黏度黏彈性薄板模型研究大洋和大陸巖石圈的彈性厚度; Lambeck 和 Nakiboglu[37]給出多層彈性板、均勻黏度黏彈性板以及彈性層疊加黏彈性層模型的解析解，并計(jì)算海洋巖石圈的有效彎曲剛度。為了更好地反映黏度隨深度的變化，文獻(xiàn)[18, 38]采用多層黏彈性模型，并解釋彈性厚度大范圍變化的觀測數(shù)據(jù)。

除解析方法外，數(shù)值模擬方法也越來越多地應(yīng)用于巖石圈彎曲問題研究。如，Zhong等[10]采用三維球坐標(biāo)有限元模型研究不可壓縮地球介質(zhì)在自重作用下的冰后回彈問題；Zhong 和 Watts[17]采用三維直角坐標(biāo)有限元模型計(jì)算巖石圈在表面載荷作用下的黏彈性變形和應(yīng)力，其中考慮了復(fù)雜的非線性黏彈性流變效應(yīng)，但同樣假設(shè)材料為不可壓縮的，與地震波速給出的巖石圈泊松比遠(yuǎn)低于0.5的觀測事實(shí)明顯不符。

本文采用黏度隨深度連續(xù)變化的黏彈性本構(gòu)模型統(tǒng)一描述彈性巖石圈與黏性巖石圈，并建立有限元模型求解在表面垂向載荷(對(duì)海山或島鏈載荷等的一般化處理)作用下的巖石圈彎曲問題。利用經(jīng)典彈性板彎曲解析解，我們將本研究黏彈性模擬計(jì)算得到的外緣隆起位置隨載荷作用時(shí)間的變化轉(zhuǎn)換成彈性巖石圈厚度隨時(shí)間的變化，并與由松弛時(shí)間定義的彈性巖石圈厚度隨時(shí)間的變化進(jìn)行定量對(duì)比。因此，本研究不但實(shí)現(xiàn)了彈性巖石圈向黏性巖石圈的逐漸自然過渡，而且能夠?qū)椥詭r石圈兩種傳統(tǒng)定義(即松弛時(shí)間定義和彈性板彎曲定義)的一致性加以嚴(yán)格檢驗(yàn)。

1 彈性巖石圈的兩種傳統(tǒng)定義及其聯(lián)系

彈性巖石圈的傳統(tǒng)定義包括彈性板彎曲定義和松弛時(shí)間定義兩種[1]。

在彈性板彎曲定義中，將彈性巖石圈之下的地幔介質(zhì)近似為應(yīng)力完全松弛因而處于靜水壓狀態(tài)的流體，而彈性巖石圈由彈性薄板近似，如圖1所示[1]。

圖1 大洋巖石圈受外載荷作用而發(fā)生彎曲Fig.1 The oceanic lithosphere deflected by an applied load

依據(jù)彈性薄板彎曲理論，利用實(shí)測巖石圈在表面載荷作用下發(fā)生彎曲導(dǎo)致的地形(如島鏈、海溝的外緣隆起位置)能夠反演巖石圈的抗彎剛度進(jìn)而確定彈性巖石圈等效厚度[39]。有效彈性厚度Te與薄板抗彎剛度D之間的定量關(guān)系為

(1)

式中：彈性參數(shù)E和ν分別為楊氏模量和泊松比。

根據(jù)松弛時(shí)間定義[1]的彈性巖石圈實(shí)質(zhì)上是要求彈性巖石圈內(nèi)的黏彈性松弛時(shí)間大于載荷作用時(shí)間，即載荷作用在彈性巖石圈中尚未松弛，而在下伏黏性巖石圈中已經(jīng)松弛。根據(jù)Turcotte和Schubert[1]的分析，初始應(yīng)力σ0松弛一半需要的時(shí)間為

(2)

假定楊氏模量E=70 GPa，松弛時(shí)間τr=107a，初始應(yīng)力σ0=100 MPa，對(duì)于干的橄欖石Ea=523 kJ·mol-1，C1=4.2×105MPa-3·s-1，由公式(2)給出彈性巖石圈底界相當(dāng)于干橄欖石溫度675 ℃，而對(duì)于濕的橄欖石Ea=398 kJ·mol-1，C1=5.5×105MPa-3·s-1，由公式(2)給出彈性巖石圈底界相當(dāng)于475 ℃等溫線。

以上兩種彈性巖石圈的定義之間沒有直接聯(lián)系?；谒沙跁r(shí)間定義的彈性巖石圈理論上合理，但因黏度難以直接測量無法廣泛應(yīng)用，而基于彈性板彎曲定義的彈性巖石圈則廣泛應(yīng)用于實(shí)際彈性巖石圈厚度的確定。

2 有限元模型

2.1 物理模型

本研究采用的物理模型如圖2所示。在線狀分布的地表垂向載荷(海山或島鏈載等)作用下的巖石圈彎曲問題近似模擬為平面應(yīng)變問題。巖石圈底部采用彈簧型邊界條件模擬靜水恢復(fù)力。上表面加載區(qū)域指定為垂向載荷線性變化，其余區(qū)域指定為自由邊界。左側(cè)邊界指定為對(duì)稱邊界，而右側(cè)邊界位移固定。

圖2 物理模型Fig.2 Physical model

我們用Maxwell模型黏彈性體代表巖石圈。假設(shè)彈性變形體積可壓縮，但彈性模量和泊松比不隨深度變化；黏度則指定為隨深度連續(xù)變化的情形[1]

(3)

式中：C為指前系數(shù)、Ea為活化能、T為絕對(duì)溫度、R為普適氣體常數(shù)。因?yàn)閹r石圈所受壓力較低，忽略壓力對(duì)黏度的影響。本文考慮的黏度變化是與Karato 和 Wu[40]的結(jié)論相符的。他們指出位錯(cuò)蠕變?cè)诟邷?例如大洋中脊)時(shí)是重要的，而低溫時(shí)擴(kuò)散蠕變占主導(dǎo)。由于海山或島鏈載荷作用的時(shí)間尺度較大(>1 Ma)，在這個(gè)時(shí)間尺度內(nèi)巖石圈有相當(dāng)部分的淺部低溫區(qū)域經(jīng)歷了應(yīng)力松弛，因此擴(kuò)散蠕變是一個(gè)合理的假設(shè)[18]。設(shè)軟流圈溫度為Tm=1 300 ℃，指前系數(shù)C可以用參考軟流圈黏度ηref表示

(4)

因而式(3)可以重新寫為

(5)

只要給定溫度隨深度的變化以及活化能Ea與參考軟流圈黏度ηref就可以用式(5)確定黏度隨深度的變化。Watts和Zhong[18]以及Watts等[16]曾用這樣的擴(kuò)散蠕變黏度公式結(jié)合板塊冷卻模型[41]和分層均勻黏彈性結(jié)構(gòu)研究巖石圈彎曲和等效彈性厚度問題。他們推薦的參數(shù)組合是Ea=120 kJ/mol，ηref=1020Pa·s。本文采用同樣的參數(shù)設(shè)定，以方便模型及結(jié)果比較。

從式(5)可以看出溫度是控制黏度的關(guān)鍵因素。我們采用板塊冷卻模型[41]計(jì)算巖石圈內(nèi)溫度分布

T(y)=T0+(T1-T0)

(6)

式中：y是深度，ts為大洋巖石圈年齡，yL0是熱巖石圈最大厚度，κ是熱擴(kuò)散系數(shù)，Tm和T0分別是軟流圈和海底溫度。模型參數(shù)取值見表1。

表1 模型參數(shù)Table 1 Modeling parameters

對(duì)年齡給定為80 Ma的大洋巖石圈，得到的溫度和黏度隨深度變化分別由圖3(a)和3(b)給出。

圖3 模擬采用的巖石圈溫度和黏度剖面Fig.3 Adopted profiles of temperature and viscosity in the modeled lithosphere

假設(shè)巖石圈上表面受到不隨時(shí)間變化的垂向載荷作用是關(guān)于y軸對(duì)稱的，為減小計(jì)算量，只考慮x>0部分的巖石圈彎曲即可。垂向載荷函數(shù)指定為

(7)

由于模擬問題是線性的，計(jì)算位移與應(yīng)力均與施加的載荷成正比，而黏彈性松弛過程及外緣隆起的時(shí)間演化均與施加載荷的絕對(duì)大小無關(guān)。

巖石圈下表面受到軟流圈靜水恢復(fù)力與巖石圈下表面垂向位移成正比，即

q(x)=kv(x),

(8)

這里k=(ρm-ρw)g，ρm和ρw分別是巖石圈之下軟流圈和巖石圈之上海水的密度，g是重力加速度，v(x) 是巖石圈下表面的垂向位移。

2.2 數(shù)學(xué)模型

巖石圈彎曲過程中，應(yīng)力張量σij滿足準(zhǔn)靜態(tài)平衡方程

σij,j=0，

(9)

位移與小應(yīng)變之間滿足幾何方程

(10)

式中：ui為i方向上的位移，而應(yīng)力與應(yīng)變由體積可壓縮的三維Maxwell體本構(gòu)方程相聯(lián)系，即

(11)

式中：λ,μ為彈性拉梅常數(shù)；η為黏度(見式(5))。利用平面應(yīng)變?chǔ)舲z=0的條件，三維問題可簡化為二維問題。

2.3 數(shù)值模型

用有限元軟件COMSOL[42]數(shù)值計(jì)算上文定義的巖石圈彎曲問題。模擬的黏彈性巖石圈為2 000 km×81 km的長方形區(qū)域。除精度檢驗(yàn)外，本文結(jié)果一律基于參考模型計(jì)算。在參考模型中，采用5階三角形單元對(duì)計(jì)算域進(jìn)行剖分，共有126 520個(gè)單元，總自由度數(shù)為8 864 562。在區(qū)域的右半部分最大單元尺寸小于2 km，而x

圖4 部分模型域的計(jì)算網(wǎng)格剖分Fig.4 Computational mesh in part of the modeling domain

2.4 計(jì)算精度檢驗(yàn)

為了驗(yàn)證數(shù)值模擬的計(jì)算精度，我們比較了不同網(wǎng)格尺寸、單元階次和時(shí)間步長的模擬結(jié)果。

在網(wǎng)格單元尺寸的檢驗(yàn)中，在參考模型網(wǎng)格的基礎(chǔ)上，將最大網(wǎng)格尺寸加倍，單元數(shù)從參考模型的126 520減少到31 906，自由度從8 864 562減少到2 237 512。在單元離散階次的檢驗(yàn)中，將單元階次從參考模型的5階減小到4階，自由度從8 864 562減少到5 826 450。在計(jì)算時(shí)間步長的檢驗(yàn)中，將參考模型每一階段采用的時(shí)間步長都增加到2.5倍，4階段的最大步長依次為250、2 500、25 000和250 000 a。

以上3種檢驗(yàn)?zāi)Ｐ偷慕Y(jié)果與參考模型比較由圖5示出。圖中對(duì)比的是上表面t=100 Ma時(shí)刻x=0處的垂向位移。網(wǎng)格變疏、離散階次減小、時(shí)間步長增大的結(jié)果在圖4中依次由加號(hào)、叉號(hào)、星號(hào)標(biāo)記對(duì)應(yīng)的線表示。結(jié)果顯示4條曲線幾乎完全重合，網(wǎng)格變疏、離散階次減小、時(shí)間步長增大的3種檢驗(yàn)?zāi)Ｐ拖鄬?duì)于參考模型相對(duì)誤差最大值的量級(jí)分別為10-6、10-8和10-9，證明本文得到的結(jié)果對(duì)計(jì)算網(wǎng)格、單元階次及時(shí)間步長的選取不敏感，具有足夠的計(jì)算精度。

圖5 數(shù)值模擬計(jì)算精度檢驗(yàn)Fig.5 Accuracy tests for numerical modeling results

3 模擬結(jié)果與討論

3.1 模擬結(jié)果

本節(jié)給出的模擬結(jié)果均基于參考模型。巖石圈上表面的垂向位移隨著載荷加載時(shí)間的變化如圖6所示。圖6(a)與6(b)分別繪出t=0+～9 Ma與t=9～100 Ma兩個(gè)時(shí)間段的垂向位移隨x坐標(biāo)的變化。圖6(a)中的t=0+解是程序給出的瞬時(shí)彈性位移解(圖8中t=0+結(jié)果表示瞬時(shí)彈性應(yīng)力)。結(jié)果表明，前一時(shí)間段外緣隆起向內(nèi)(加載一側(cè))運(yùn)動(dòng)，而后一時(shí)間段外緣隆起向外移動(dòng)，即隨著載荷作用時(shí)間的增加，外緣隆起與載荷中心處的距離先減小后增大。

圖6 巖石圈上表面的垂直方向位移隨載荷作用時(shí)間的變化Fig.6 Vertical displacement of the upper surface of the lithosphere as a function of loading time

上表面載荷中心處(x=0)的垂向位移，以及上表面x=250 km處的垂向位移與水平位移隨時(shí)間的連續(xù)變化如圖7所示。圖7除清楚表明載荷中心與外緣隆起相反的垂向運(yùn)動(dòng)外，也清楚地反映了外緣隆起區(qū)域存在先內(nèi)后外的水平運(yùn)動(dòng)特征。

圖7 計(jì)算位移的時(shí)間演化Fig.7 Evolution of the computed displacement

數(shù)值模擬得到的黏彈性應(yīng)力松弛過程由圖8給出。圖中繪出的是x=0處水平和垂向正應(yīng)力σxx,σyy以及x=250 km處的水平方向正應(yīng)力σxx和剪應(yīng)力τxy在不同時(shí)刻隨深度的變化。圖8表明，巖石圈中應(yīng)力的黏彈性松弛是一個(gè)持續(xù)的動(dòng)力學(xué)過程。在加載之后，巖石圈整個(gè)厚度均產(chǎn)生彈性響應(yīng)，水平正應(yīng)力隨深度近乎線性分布，與彈性薄板假設(shè)很接近。隨著載荷作用的持續(xù)，巖石圈深部的應(yīng)力逐漸松弛，而淺部的正應(yīng)力持續(xù)增加，巖石圈的彈性厚度因而持續(xù)減小。外緣隆起下方區(qū)域的剪應(yīng)力呈現(xiàn)先增加后減小的趨勢。這樣的剪應(yīng)力演化過程與外緣隆起區(qū)域先內(nèi)后外的水平運(yùn)動(dòng)特征是一致的。

圖8 巖石圈黏彈性應(yīng)力松弛過程Fig.8 Processes of viscoelastic stress relaxation in the lithosphere

3.2 討論

Turcotte和Schubert[1]給出利用巖石圈彈性彎曲產(chǎn)生的外緣隆起位置反推出彈性巖石圈厚度的經(jīng)典公式

(12)

由外緣隆起的位置xb可計(jì)算彈性巖石圈厚度Te；Te隨xb的單調(diào)變化圖像如圖9(a)所示。圖9(b)給出本研究計(jì)算得到的外緣隆起位置隨時(shí)間的演化規(guī)律。圖9(b)表明，外緣隆起與載荷中心的距離先是持續(xù)減小，直到t=9 Ma達(dá)到最小值242 km。而后，外緣隆起與載荷中心的距離又隨時(shí)間持續(xù)單調(diào)遞增直到1 000 Ma。

結(jié)合圖9(a)和9(b)，即將本文模擬得到的外緣隆起位置變化由公式(12)轉(zhuǎn)換成彈性巖石圈厚度，得到圖10所示的Te隨加載時(shí)間的演化曲線。由于外緣隆起區(qū)域先內(nèi)后外的水平運(yùn)動(dòng)導(dǎo)致圖10中Te先隨加載時(shí)間減小而后又隨時(shí)間增加。顯然，彈性厚度隨時(shí)間增加階段的結(jié)果是違反物理規(guī)律的。因此，可以得到結(jié)論：前人廣泛采用的利用外緣隆起位置反演彈性巖石圈厚度的做法具有嚴(yán)重缺陷，在載荷作用時(shí)間超過9 Ma以后會(huì)嚴(yán)重高估彈性巖石圈厚度。

事實(shí)上，即使在載荷作用時(shí)間小于9 Ma的時(shí)間內(nèi)，利用外緣隆起位置反演得到的彈性巖石圈厚度與根據(jù)松弛時(shí)間定義的彈性厚度也有明顯的區(qū)別。根據(jù)Turcotte和Schubert[1]的公式

(13)

我們計(jì)算了本文所用黏度剖面對(duì)應(yīng)的深度yEL隨松弛時(shí)間的變化關(guān)系，并繪制在圖10中。如果用松弛時(shí)間大于載荷作用時(shí)間作為彈性巖石圈的定義，顯然松弛時(shí)間曲線對(duì)應(yīng)于彈性厚度的另一種定義。

圖10清楚地表明，兩種定義給出的彈性厚度僅在0.2～5.0 Ma的時(shí)間范圍內(nèi)較為符合(差異不超過0.64 km， 且在t=0.30 Ma和t=4.5 Ma時(shí)兩線重合)。在載荷作用時(shí)間小于0.2 Ma時(shí)，兩種定義有顯著差別但總體趨勢一致。在t=0.001 Ma時(shí)兩條曲線差距最大(4.4 km)，接著差異慢慢減小到零(t=0.001 8 Ma)，隨后逐漸增大到3.64 km(t=0.008 Ma)后又逐漸減小。在t>5 Ma時(shí)，兩種定義給出的彈性巖石圈厚度差異隨時(shí)間快速增大。由于松弛時(shí)間定義的彈性厚度yEL是直接由模型給定的光滑單調(diào)變化的黏度剖面計(jì)算得到的，圖10中兩條曲線間差距的復(fù)雜交替變化完全由巖石圈黏彈性松弛導(dǎo)致的外緣隆起移動(dòng)特征造成。由圖9(b)可以看出，外緣隆起位置隨加載時(shí)間對(duì)數(shù)的變化先后經(jīng)歷了斜率增加到接近常斜率、而后斜率又逐漸減小為零、最后斜率反向并逐漸增加的復(fù)雜演變。正是巖石圈黏彈性松弛的復(fù)雜性導(dǎo)致彈性巖石圈兩種定義給出的厚度曲線出現(xiàn)多次交叉的現(xiàn)象。

圖10 彈性巖石圈兩種傳統(tǒng)定義給出的厚度對(duì)比Fig.10 Comparison between the thickness values calculated based on two classical definitions of elastic lithosphere

4 結(jié)論

本文利用變黏度可壓縮黏彈性有限元模型定量研究巖石圈彎曲的黏彈性松弛問題。模擬結(jié)果顯示，外緣隆起位置具有先向彎曲中心運(yùn)動(dòng)后向彎曲外側(cè)移動(dòng)的特點(diǎn)，反向運(yùn)動(dòng)發(fā)生在載荷作用約9 Ma后。由于經(jīng)典彈性板彎曲理論給出外緣隆起與彈性巖石圈厚度的單調(diào)正相關(guān)關(guān)系，據(jù)此反演的彈性巖石圈厚度不能應(yīng)用于載荷作用超過9 Ma情形。否則，會(huì)出現(xiàn)彈性巖石圈厚度隨加載時(shí)間增加的荒謬結(jié)果，違反巖石圈黏彈性松弛的物理規(guī)律。進(jìn)一步與基于松弛時(shí)間的彈性巖石圈厚度定義進(jìn)行對(duì)比，本文模擬結(jié)果表明，由基于彈性板彎曲解利用外緣隆起位置約束的彈性巖石圈厚度僅在載荷持續(xù)0.2～5.0 Ma的時(shí)間段內(nèi)能夠給出較為合理的結(jié)果，誤差不超過0.64 km。在載荷作用時(shí)間小于0.2 Ma時(shí)，基于兩種傳統(tǒng)定義給出的彈性巖石圈厚度總體趨勢一致，但差異最大可達(dá)4.4 km。而在載荷作用時(shí)間大于5 Ma時(shí)，兩種定義給出的彈性厚度差異隨時(shí)間快速增大。兩種傳統(tǒng)定義給出的彈性巖石圈厚度曲線出現(xiàn)多次交叉的現(xiàn)象，明確揭示巖石圈彎曲的黏彈性松弛是一個(gè)復(fù)雜的動(dòng)力學(xué)過程，因此在確定彈性巖石圈厚度的地球動(dòng)力學(xué)相關(guān)研究中應(yīng)當(dāng)充分考慮黏彈性松弛的影響。