克拉通巖石圈對流減薄的數(shù)值模擬

程華冬， 黃金

1 中國科學技術(shù)大學地球和空間科學學院地震與地球內(nèi)部物理實驗室， 合肥　230026 2 蒙城地球物理國家野外科學觀測研究站， 合肥　230026

克拉通巖石圈對流減薄的數(shù)值模擬

1 中國科學技術(shù)大學地球和空間科學學院地震與地球內(nèi)部物理實驗室， 合肥230026 2 蒙城地球物理國家野外科學觀測研究站， 合肥230026

摘要采用二維有限元數(shù)值模擬的方法研究了巖石圈的對流減薄過程，特別是克拉通巖石圈的對流減薄過程.模型的主要參數(shù)包括增厚巖石圈的寬度x、增厚倍數(shù)γ、以及與巖石圈組分變化導致的黏性和密度變化密切相關(guān)的黏性比(ηc)和浮力數(shù)(B)或等效密度變化(Δρtc).數(shù)值計算結(jié)果顯示，地幔對流將逐漸減薄增厚的巖石圈部分，(1)當B=0和ηc=1時，即對一般地幔巖石圈，增厚巖石圈對流減薄的時間可表示為0.0073γ0.70x0.26.將數(shù)值結(jié)果應(yīng)用于地球，意味著增厚到300 km的巖石圈，如寬度為300 km，對流移除增厚部分回到初始平衡厚度120 km大約需要225 Ma；如寬度為1500 km，移除增厚部分大約需要342 Ma.(2)當B和ηc較小，克拉通巖石圈對流減薄過程與一般加厚巖石圈的對流減薄過程類似，但減薄時間受克拉通組分浮力和黏性比的影響而顯著增長，克拉通巖石圈對流減薄的時間可表示為γ0.78ηc-0.36x0.04.因而，對300 km厚的克拉通巖石圈，如克拉通巖石圈的密度比周圍地幔的密度低0.4%(即B=0.1)，寬度1500 km，若克拉通巖石圈黏性因組分影響比普通地幔巖石圈大10倍，其被對流減薄到120 km大約需要1.18 Ga.(3)當B和ηc增大到一定量時(如B≥0.2且ηc>10)，克拉通巖石圈被移除的過程將發(fā)生變化，由于組分浮力的影響，對流主要不是將克拉通巖石圈帶到軟流圈地幔中，而主要是將較厚的巖石圈物質(zhì)向兩邊推送.在此情況下，克拉通巖石圈能長時間(>3 Ga)保持穩(wěn)定.

關(guān)鍵詞克拉通巖石圈； 巖石圈減薄； 巖石圈減薄移除時間

1引言

克拉通是大陸上相對穩(wěn)定的構(gòu)造單元，包括地殼和巖石圈地幔(Rudnick and Fountain, 1995;Schubert et al., 2001).地質(zhì)、地球物理和地球化學研究顯示：克拉通有200多公里的巖石圈根，有的厚達300到400 km；克拉通區(qū)域表面熱流值小；巖石圈根相對較冷、含水量低，并且顯示出是部分熔融分異后的殘留物；巖石圈根的密度相對較小、黏性較高；克拉通形成后可長期保持穩(wěn)定，有的可長達40億年(Boyd et al., 1985;Carlson et al., 2005;King, 2005;Pearson, 1999;Pearson et al., 1995).但是有些克拉通(如美國的Wyoming克拉通、我國的華北克拉通)在其形成后遭遇過改造(Carlson et al., 2004;Gao et al., 2002).其中華北克拉通是目前研究最為廣泛和深入的地區(qū).在華北克拉通發(fā)現(xiàn)了38億年的陸殼殘留物(Liu et al., 1992).在經(jīng)歷了早元古代(約25—18億年間)的增生聚合碰撞事件后隨著巖石圈厚度增大,華北克拉通作為一個完整的大陸塊體長時期處于相對穩(wěn)定的狀態(tài)(Zhao et al., 2001).金伯利巖研究顯示約在早古生代(～4.6億年前)華北克拉通地幔巖石圈根的厚度約為200 km(路鳳香等, 1991).大約在中侏羅世(～1.6億年)中國東部開始大規(guī)?；鹕絿姲l(fā).目前的初步研究認為華北克拉通巖石圈減薄發(fā)生在晚中生代，在早白堊世(1.2—1.3億年)達到高潮，巖石圈的厚度從200 km減薄到80～120 km左右.伴隨著這一過程，大量的火山活動使地幔巖石浸入地殼，地殼運動非?；钴S，導致古老的華北克拉通的解體或破壞(Griffin et al., 1998;Menzies et al., 1993;Wu et al., 2005;Xu, 2001;鄧晉福等, 2003;范蔚茗和Menzies, 1992;李江海等, 2004;邵濟安等, 2002;吳福元等, 2003;許文良等, 2004;趙越等, 2004).

盡管目前對造成克拉通巖石圈減薄的機制和過程仍有許多爭論，但普遍認為對流減薄在克拉通巖石圈減薄過程中具有重要影響.克拉通巖石圈減薄的概念性模型包括拆沉、對流侵蝕、熱侵蝕、橄欖體和熔體的相互作用、巖漿提取、巖石圈地幔水化、以及機械拉張等(吳福元等, 2008).不同機制導致的巖石圈減薄后的后果和地表反應(yīng)存在區(qū)別(吳福元等, 2008;朱日祥等, 2012).各種概念模型都有相應(yīng)的地質(zhì)、地球物理和地球化學觀測等數(shù)據(jù)支撐(Carlson et al., 2005;Duggen et al., 2003;King, 2005;Schurr et al., 2006;Zhao and Zheng, 2007;吳福元等, 2008;許文良等, 2004).但對華北克拉通，目前的研究結(jié)果顯示，由于巖石圈俯沖的影響，華北克拉通區(qū)域及其下覆地幔受到影響，軟流圈地幔和巖石圈地幔的相互作用是導致華北克拉通被破壞的主因(朱日祥等, 2011, 2012).因而研究軟流圈地幔對流對巖石圈地幔的影響具有重要意義(朱日祥等, 2012).

研究地幔對流對巖石圈的影響一直是地球動力學研究所關(guān)心的話題(Buck and Toks?z, 1983;Davies, 1994;Houseman et al., 1981;Huang et al., 2003;Marotta et al., 1999;Morency et al., 2002;Olson et al., 1988;Schott and Schmeling, 1998;Shapiro et al., 1999;Yuen and Fleitout, 1985;喬彥超等, 2012, 2013).數(shù)值模擬結(jié)果顯示，加厚巖石圈減薄的過程和時間與模型參數(shù)密切相關(guān).不同模型由于關(guān)注的對象不同，模型參數(shù)各異，導致得到的巖石圈減薄時間尺度變化很大，從～10 Ma(Marotta et al.,1999)到幾百個Ma(Buck and Toks?z, 1983;Houseman et al., 1981)、甚至能夠長期保持穩(wěn)定(Shapiro et al., 1999).因而，定量研究加厚巖石圈的對流減薄和模型參數(shù)的關(guān)系，顯得十分重要.

Morency等(2002)通過數(shù)值計算和理論分析的方法，定量分析了加厚巖石圈根在對流地幔中被移除的時間.其結(jié)果顯示，加厚巖石圈根的移除時間與其大小(即其二維模型中的加厚巖石圈根的厚度和寬度)和黏性參數(shù)密切相關(guān)，并給出了估算公式.這對討論增厚巖石圈的減薄過程及其對地表的影響等具有重要意義.然而在Morency等(2002)的模型中并未考慮克拉通地幔巖石圈和普通地幔巖石圈之間可能存在的由化學組分不同引起的密度和黏性差異，而這正是克拉通巖石圈和一般地幔巖石圈間的差異所在.現(xiàn)有研究都顯示，克拉通巖石圈的密度和黏性差異對其穩(wěn)定性具有重大影響(Lenardic and Moresi, 1999;Lenardic et al., 2003;O′Neill et al., 2008;Wang et al., 2014;Yoshida, 2012).因而本文將在Morency等(2002)的定量化方法的研究基礎(chǔ)上，通過在數(shù)值模型中增加考慮由化學組分引起的巖石圈的密度和黏性變化的影響，進一步研究克拉通巖石圈在對流地幔中的減薄過程和減薄時間，為討論克拉通巖石圈的對流減薄提供數(shù)值參考模型.

2數(shù)值模型

(1)

其中，u，P，η，T和H分別為流體速度矢量、壓力、黏性系數(shù)、溫度和內(nèi)生熱，C代表組分場，初始時刻，巖石圈內(nèi)C=1，其余部分C=0.ez為z方向上的單位矢量，Ra為Rayleigh數(shù)，B是組分浮力數(shù)，其定義分別為：

(2)

其中，ρ0是地幔的參考密度，α是熱膨脹系數(shù)，g是重力加速度，ΔT是模型上下邊界溫度差，D是模型深度，κ是熱擴散系數(shù)，η0是底部邊界的參考黏性，Δρc=ρ0-ρlith是組分密度差，ρlith指克拉通巖石圈的密度.方程組(1)無量綱化時的參數(shù)為：長度[L]=D；時間[t]=D2/κ；黏性[η]=η0；溫度[T]=ΔT；內(nèi)生熱率[H]=(CpκΔT)/D2.

模型假定地幔黏性與溫度相關(guān)，其無量綱形式為：

(3)

式中，E是無量綱的活化能，與活化能E*的關(guān)系：E=E*/(RΔT)，T0=Ts/ΔT是無量綱的表面溫度，ηz表示黏性隨深度的變化，這里主要用于表示軟流圈.當無量綱化的深度Z<0.66時ηz=1/30，其他區(qū)域為ηz=1.ηc是組分黏性比，表示巖石圈組分變化導致的黏性變化.C=0的組分ηc=1，C=1的區(qū)域為巖石圈，通過改變ηc來模擬克拉通巖石圈的黏性變化.

數(shù)值模型是3×1的二維模型，方程(1)采用二維有限元程序Citcom解算.模型劃分為289×97個有限單元網(wǎng)格，在靠近上下邊界處，適度進行了網(wǎng)格加密.上表面給定溫度，下表面熱流為0，上表面固定，下表面和側(cè)面邊界自由滑移.本文模型深度取值1000 km，一方面能反映出軟流圈與下地幔的黏性跳變，另一方面則是增大分辨率，模型選用的其他參數(shù)值見表1.根據(jù)(2)式，模型瑞利數(shù)Ra=1.3×107.由于我們研究的是擬穩(wěn)狀態(tài)下的對流模型，考慮到地球的冷卻效應(yīng)，為保證地幔溫度保持大約1250 K，模型選取的內(nèi)生熱略大于地球內(nèi)生熱(Turcotte and Schubert, 2002).

表1　模型參數(shù)值

為了定量研究克拉通巖石圈根的對流減薄時間，這里參考Morency等(2002)的做法，即：(1)首先計算一個參考模型，使其達到擬穩(wěn)狀態(tài)，此時模型的表面熱流、內(nèi)部溫度以及均方根速度等基本保持不變.由于E很大，模型表現(xiàn)為近似剛性蓋層下的小尺度對流形式(圖1(a—d)).該邊界層厚度，Zq，用以刻畫巖石圈厚度，其值按橫向平均溫度T≤0.9Ti計算，Ti定義為無量綱模型深度在0.3到0.8之間的溫度平均值.(2)克拉通巖石圈根通過增加參考模型局部區(qū)域巖石圈厚度來模擬(圖1(e—h)).具體方法是在模型中部選取一定寬度(x)的區(qū)域?qū)⒃搮^(qū)域厚度從Zq增加到Z0，Z0=(γ+1)Zq.由于邊界層基本處于熱傳導狀態(tài)，內(nèi)部溫度近似成線性關(guān)系，增厚邊界層內(nèi)的溫度隨深度的變化關(guān)系和增厚前一樣，只是系數(shù)減小1/(γ+1).此時，增厚巖石圈相比于初始時，厚度增加了γZq，這里我們將γ稱之為增厚倍數(shù).(3)本文的一個重要任務(wù)是確定巖石圈

圖1　參考模型在擬穩(wěn)狀態(tài)下的等溫線和流線(a)以及橫向平均溫度(b)、黏性(c)和均方根速度Vrms(d).模型caseB3在初始狀態(tài)下的等溫線和流線(e)以及橫向平均溫度(f)、黏性(g)和Vrms(h)(f—h)中虛線是整個模型的橫向平均值，實線是巖石圈加厚部分的橫向平均值.(i—l)顯示加厚巖石圈演化過程的溫度場和流場.(a，e，i—l)中，實線代表等溫線，相鄰等溫線間距為無量綱溫度0.1，黑色點線表示流函數(shù)值小于零的流線，黑色虛線表示流函數(shù)值大于等于零的流線.Fig.1　Isotherm, streamline (a) and radial dependences of horizontal average temperature (b), viscosity (c) and Vrms (d) at the quasi stable state for the reference model. Isotherm, streamline (e) and radial dependences of horizontal average temperature (f), viscosity (g) and Vrms (h) at the initial state of the Model caseB3 The dashed line are displayed for the horizontal average value of the whole model and the solid line are displayed for the horizontal average of the thickened lithosphere in (f—h). (i—l) The temperature and flow field evolution of the thickened lithosphere. Solid line are displayed for isotherm, the spacing of adjacent isotherm is dimensionless temperature 0.1, black dotted line for positive streamline and black dashed line for negative streamline in (a，e，i—l).

對流減薄的時間.圖1(e，i—l)分別給出了模型caseB3在巖石圈加厚后的溫度場和流場的時間演化過程.從圖中可以看到，經(jīng)過一定時間的演化，增厚巖石圈厚度從(γ+1)Zq基本重新減薄到Zq.增厚巖石圈根移除時間的計算方法仍采用Morency等(2002)的方法，即計算每一時間步加厚巖石圈區(qū)域在深度Z=0.9Zq處的橫向平均溫度進而給出時間變化曲線(圖2)，然后利用溫度演化公式(Morency et al., 2002)：

(4)

對模型數(shù)據(jù)進行擬合，得到Ta、ta、Tst和A.其中ta是延遲時間，表示溫度變化只有在此時間之后才符合(4)式；Ta相當于ta時間時，加厚巖石圈區(qū)域在Z=0.9Zq處的平均溫度；Tst則相當于達到平衡時加厚巖石圈區(qū)域在Z=0.9Zq處的平均溫度.巖石圈根移除時間定義為：

(5)

必須注意的是，當γ較小時，延遲時間ta近似為零(圖2中的實線)，但當γ較大時，延遲時間ta就增大(圖2中的虛點線).依此方法確定的模型巖石圈根移除時間計算結(jié)果見表2和表3.

表2　模型參數(shù)和結(jié)果

表3　模型參數(shù)和結(jié)果

續(xù)表3

注：*增厚的巖石圈根不能被完全移除，沒有計算地幔內(nèi)部溫度Ti及移除時間t2D；#增厚的巖石圈根被移除掉，由于浮力數(shù)較大，0.9Zq處的橫向平均溫度變化復雜，不能使用上面的方法得到移除的具體時間，但能夠大致得到移除時間在0.07以上.

圖2　深度Z=0.9Zq處的橫向平均溫度的時間演化曲線不同線型(點線除外)表示在不同的增厚倍數(shù)γ下的數(shù)值計算值，而通過對應(yīng)線型的點線代表對應(yīng)數(shù)值計算值的擬合曲線.Fig.2　Horizontal average temperature evolution with time at Z=0.9ZqThe line (except dotted line) in different linear are displayed for the numerical values of the different thickening factors γ and the dotted line corresponds to the line in different linear are displayed the fitting curves corresponds to the numerical values.

3數(shù)值結(jié)果

增厚巖石圈的對流減薄與巖石圈的物性參數(shù)密切相關(guān)(Huang et al., 2003;Morency et al., 2002).瑞利數(shù)和活化能等對巖石圈減薄過程具有重要影響(Huang et al., 2003;Morency et al., 2002)，但由于Morency等(2002)在其工作中詳細分析和討論了Ra(瑞利數(shù))和E(活化能)對增厚巖石圈移除時間的影響，為能充分討論克拉通巖石圈組分差異導致的密度和黏性的影響，本文將不再考慮Ra和E變化的影響，而是將其作為模型給定參數(shù)(參見表1).根據(jù)表1的模型參數(shù)，本文參考模型在處于擬穩(wěn)狀態(tài)時的熱邊界層厚度Zq=0.12.下面首先將給出B=0和ηc=1的一組模型的結(jié)果，該組模型忽略克拉通巖石圈固有的密度和黏性差異.然后通過改變B和ηc值研究克拉通巖石圈固有的密度和黏性差異的影響.即在討論并總結(jié)一般巖石圈對流減薄規(guī)律的基礎(chǔ)上，進一步討論分析克拉通巖石圈對流減薄的規(guī)律.

3.1加厚巖石圈根的對流減薄

對一般加厚巖石圈，在本文模型中即假定B=0和ηc=1.該組模型與Morency等(2002)的模型基本一致.但應(yīng)該注意，與Morency等(2002)模型不同的是，為和實際地球更加接近，本文考慮了660 km處的黏性躍變(Hager and Richards, 1989).表2給出了γ變化從0.5到2和x變化從0.4到3的共20個算例的結(jié)果.

圖1(e，i—l)顯示了模型caseB3的演化過程，該模型增厚巖石圈的寬度為1.5，增厚倍數(shù)為1.0.圖1e是模型caseB3在初始時刻的等溫線和流線.從圖中可以看到，增厚巖石圈的減薄過程受到側(cè)面和底部對流的影響(圖1i和1j).在兩側(cè)以及底部對流的同時作用下，增厚的巖石圈逐步被移除、抹平(圖1(i—k)).在0.0073時刻，增厚部分的巖石圈基本被移除(圖1k).這與表2給出的、依據(jù)公式(4)和(5)計算得到的模型caseB3的移除時間0.0072基本一致.下面將詳細分析增厚巖石圈移除時間與γ和x的關(guān)系.

圖3a給出了對流移除時間與增厚倍數(shù)γ的對數(shù)關(guān)系圖.圖中相同形狀表示x相同的模型，連接相同形狀的直線是通過線性擬合得到的相同x值下增厚巖石圈移除時間和增厚倍數(shù)的關(guān)系.從圖中可以看到：對相同的x，對流移除時間與增厚倍數(shù)基本成冪函數(shù)關(guān)系，對流移除時間隨增厚倍數(shù)增大而增大；對不同的x，擬合直線的斜率略有差異，但基本一致.計算結(jié)果顯示：

(6)

圖3b給出了對流移除時間與增厚寬度x的對數(shù)關(guān)系圖.圖中相同形狀表示γ相同的模型.連接相同形狀的直線是線性擬合得到的相同增厚倍數(shù)時，移除時間和增厚巖石圈寬度的關(guān)系.從圖中可以看出，對流移除時間和增厚寬度也基本成冪函數(shù)關(guān)系，具體計算結(jié)果顯示：

(7)

綜合(6)式和(7)式的分析，我們得到對流移除時間與增厚倍數(shù)γ和增厚寬度x之間的量化關(guān)系：

(8)

圖3c給出了根據(jù)(8)式計算的數(shù)值和表2中給出的模型結(jié)果的比較.從圖中可以看到，式(8)可以很好地表示模型的計算結(jié)果，只是在移除時間很大和很小時略微存在差異.

根據(jù)(8)式可以方便地估算一個增厚巖石圈被對流減薄的時間.對于初始平衡厚度為120 km的巖石圈，將其厚度增加到300 km(即γ=1.5)，當增厚巖石圈寬度為300 km(x=0.3)時，增厚巖石圈根被移除的時間為225 Ma；而當加厚巖石圈寬度為1500 km(x=1.5)時，則其增厚巖石圈根的移除時間是342 Ma；而對于更大的加厚巖石圈寬度，如10000 km(x=10.0)，其根的移除時間是559 Ma.

盡管本文模型的黏性結(jié)構(gòu)與Morency等(2002)的模型略有差別，但兩者的結(jié)果仍具有一致性.Morency等(2002)通過數(shù)值計算和理論分析認為，加厚巖石圈的減薄包括側(cè)面和底部的侵蝕影響，當加厚巖石圈寬度很小時，以側(cè)面侵蝕為主；而當加厚巖石圈寬度很大時，以底部侵蝕為主.忽略次要影響的前提下，他們給出的無量綱關(guān)系式為：

當x

(10)

式中，x0是反映底部和側(cè)面侵蝕作用的特征寬度，

(11)

式(9)—(11)中的變量與本文一致，其中V是無量綱的活化體積(無量綱參數(shù)為[V]=(RΔT)/(ρ0gD))，ΔTe是Morency等(2002)引進的流變溫度(Davaille and Jaupart, 1994;Grasset and Parmentier, 1998;Solomatov and Moresi, 2000)，即

(12)

式中

(13)

如前面提到的，本文采用在深度大于660 km的下地幔黏性增加30倍來表示黏性隨深度的變化，

圖3　增厚巖石圈對流移除時間與增厚倍數(shù)γ(a)和寬度x(b)的關(guān)系；(c) 模型計算結(jié)果和擬合結(jié)果的比較Fig.3　Convective removal duration of the thickened lithosphere scale with the thickening factor γ (a) and width x (b)；(c) Comparison of the model results and the fitting results

而Morency等(2002)采用活化體積V*來表示黏性隨深度的變化.要達到巖石圈下覆軟流圈的黏性與下地幔的平均黏性相差約30倍，V*的取值大約應(yīng)為0.75 cm3·mol-1，相應(yīng)的無量綱的活化體積V=2.1621.根據(jù)表1及表2的數(shù)據(jù)，通過式(12)和(13)計算得到ΔTe=0.183.將此V及ΔTe值和本文模型的其他參數(shù)值代入(9)—(11)式，并取γ=1.5可計算得到無量綱的時間：

(14)

根據(jù)前面給出的時間和長度的量綱，式(14)量綱化，即為：

(15)

式中x*以km計.這意味著，對于初始平衡厚度為120 km的巖石圈，將其厚度增加到300 km(即γ=1.5)，當增厚巖石圈寬度為300 km時，其增厚巖石圈根被移除的時間為121 Ma；而當加厚巖石圈寬度為1500 km，則其根的移除時間是428 Ma.而對于更大的加厚巖石圈寬度，如10000 km，根據(jù)(15)式，其根的移除時間仍為428 Ma.

從上面的結(jié)果中可以看到，Morency等(2002)的結(jié)果與本文模型結(jié)果存在一定差異.這主要是由于Morency等(2002)在分析數(shù)據(jù)時的理論模型過于簡單所導致的.如圖1所示，增厚巖石圈根的移除是側(cè)面和底部對流侵蝕同時作用的結(jié)果，但在Morency等(2002)的模型中，以特征寬度為界限，寬度小時忽略了底部的作用，而寬度大時又忽略了側(cè)面的作用，所以他們的結(jié)果顯示在寬度小時，移除時間與寬度成正比，而寬度較大時，移除時間與寬度無關(guān).

3.2克拉通巖石圈根的對流減薄

本文關(guān)注的焦點是克拉通巖石圈根的對流減薄，這就涉及到克拉通巖石圈的密度和黏性參數(shù)的變化.因此，研究克拉通巖石圈根的對流減薄，即使不考慮Ra和E的影響，模型至少也包含有與克拉通巖石圈根的大小(γ,x)和化學組分導致的密度與黏性變化(B,ηc)有關(guān)的四個參量.前面已經(jīng)研究了不考慮化學組分導致的密度與黏性變化情況下的增厚巖石圈根在熱對流下的對流減薄問題，并得到了對流減薄時間和增厚的寬度以及厚度的量化關(guān)系.這里將在此基礎(chǔ)上，通過改變B和ηc來模擬化學組分變化導致的克拉通巖石圈的密度和黏性變化，以討論它們對克拉通巖石圈根的對流減薄的影響.表3給出了這組模型的計算結(jié)果.

比較表2和表3中移除時間的大小，我們發(fā)現(xiàn)在模型中增加克拉通巖石圈的密度和黏性變化后，增厚巖石圈根更難以被移除.但是，即使考慮克拉通巖石圈的密度和黏性的影響，當B和ηc較小時，巖石圈根的移除過程和3.1節(jié)模型所示的普通巖石圈根移除過程基本類似(如圖4).巖石圈根的移除仍然是側(cè)面和底部對流減薄共同作用的后果.圖4顯示了模型caseB3_7的溫度、流場和克拉通巖石圈組分隨時間演化的過程.

從圖4可以看到模型溫度場的演化和普通增厚巖石圈模型的(圖1(e，i—l))演化相似，底部和側(cè)面的對流侵蝕著巖石圈根，使其慢慢變薄，直至重新達到穩(wěn)定.從組分場(圖4(g—l))也可以看到，克拉通巖石圈組分被側(cè)面和底部的對流帶到了下部地幔并混入其中，增厚的克拉通巖石圈根逐漸被移除.表3顯示caseB3_7的移除時間為0.0266，從圖4中也可以看到，圖4f和圖4l對應(yīng)的時間為0.0259，此時加厚的克拉通巖石圈根幾乎被完全移除.

但當B或ηc較大時，移除時間將顯著增長(表3).圖5顯示的是caseB3_14、caseB3_15、caseB3_16、caseB3_20和caseB3_23的組分時間演化圖.該組模型的克拉通巖石圈根大小(γ=1.0,x=1.5)相同，克拉通組分變化導致的黏性變化ηc以及浮力數(shù)B的數(shù)值如圖所示.從圖5中可以明顯看到，當B=0.2時，隨著ηc的增大，增厚的巖石圈根的移除時間不斷加長(圖5中的第1,2,3行)，這也清楚地顯示在巖石圈內(nèi)溫度場的演化上(圖6).圖6給出模型caseB3_14-17、caseB3_20和caseB3_23在深度Z=0.9Zq處的橫向平均溫度演化曲線.一般來說，該處的溫度演化基本滿足(4)式的指數(shù)規(guī)律(圖2).圖6中顯示，在B=0.2且黏性比較小時(caseB3_14和caseB3_15)，平均溫度曲線基本滿足(4)式.擬合得到的增厚的巖石圈根移除時間分別為0.0282(894 Ma)和0.0503(1.60 Ga)(表3).但當黏性比較大時，在0.1的時間內(nèi)，Z=0.9Zq處的橫向平均溫度一直處于上升階段(圖6中caseB3_16，ηc=10).這表明，在0.1的時間內(nèi)克拉通巖石圈并沒有被充分移除，這也反映在圖5的組分場上(圖5中的第3行).

圖5還顯示，當B或ηc足夠大時(如ηc=10；B=0.3和0.4)，克拉通巖石圈的對流減薄過程與前述過程表現(xiàn)的非常不一樣(圖5中的第4,5行).圖5顯示，此時很少克拉通巖石圈組分被帶入其下的地幔中，克拉通巖石圈組分在漫長的演化過程中，向兩邊擴散，這也會使得克拉通巖石圈得到部分減薄(圖5中的第4,5行)，但這一過程非常漫長.從模型caseB3_20和caseB3_23在深度Z=0.9Zq處的橫向平均溫度演化曲線(圖6中的較粗的實線和虛線)上看，平均溫度在0.1時間內(nèi)仍處于上升階段.很顯然，如果平均溫度一直處于上升，就無法利用(4)式進行擬合，從而無法確定對流減薄時間.無量綱0.1的時間大約相當于3.17 Ga，這個時間非常長，如果模型在這么長時間內(nèi)無法計算對流移除時間，我們就不將這類模型納入后面的統(tǒng)計分析中.

圖4　caseB3_7在不同時刻的(a—f)等溫線及流線和(g—l)組分場演化(a—f)中實線代表等溫線，相鄰等溫線間距為無量綱溫度0.1.黑色點線表示流線為負值，黑色虛線表示流線為正值.Fig.4　The isotherm, streamline (a—f) and composition (g—l) evolution of the model caseB3_7, displayed at different time Solid line for isotherm, the spacing of adjacent isotherm is dimensionless temperature 0.1.Black dotted line for positive streamline and black dashed line for negative streamline in (a—f).

圖5　不同組分黏性比以及浮力數(shù)下，增厚的克拉通巖石圈根的組分場演化.增厚克拉通巖石圈根的增厚倍數(shù)γ=1.0，增厚寬度x=1.5Fig.5　The composition of the thickening cratonic lithosphere, after stretching factor γ=1.0, width x=1.5, displayed at different viscosity ratio and buoyance number

圖6　模型caseB3_14-17、caseB3_20和caseB3_23在深度Z=0.9Zq處的橫向平均溫度演化曲線實線、虛線、點劃線、點線以及較粗的實線、虛線分別表示caseB3_14—17、caseB3_20和caseB3_23.Fig.6　Horizontal average temperature evolution with time at Z=0.9ZqSolid，dashed，dot-dashed，dotted line and the thicker solid，dashed line are displayed for caseB3_14—17，caseB3_20和caseB3_23.

(16)

(17)

其中，ΔTe是通過式(12)和(13)計算得到的，表3中的Ti略有變化，因此熱化學對流模型的ΔTe也會變化.下面將仿照3.1節(jié)模型數(shù)據(jù)的分析方法對表3中數(shù)據(jù)進行分析.在我們的熱化學對流模型中，除了增厚倍數(shù)γ和增厚區(qū)域?qū)挾葂外，還包括組分浮力數(shù)B和組分變化導致的黏性變化ηc.將組分浮力數(shù)B按照(17)式轉(zhuǎn)化為Δρtc，則模型的四個參數(shù)為γ、x、ηc和Δρtc.參考3.1節(jié)的分析結(jié)果，假定：

(18)

下面將根據(jù)表3的結(jié)果確定(18)式中的冪指數(shù).

圖7a給出的是增厚區(qū)域?qū)挾纫欢l件(x=1.0)下對流移除時間隨著巖石圈增厚倍數(shù)γ的對數(shù)關(guān)系圖，其中相同形狀表示相同的組分黏性比ηc(圓圈表示ηc=1，方形表示ηc=3，菱形表示ηc=10)，而相同的線型則代表著相同的浮力數(shù)B(實線表示B=0，點劃線表示B=0.05，虛線表示B=0.1)，通過相同形狀的線條是其線性擬合.從圖中能夠看到：總體的趨勢是對流移除時間隨著增厚倍數(shù)的增大而增大，這和熱對流的結(jié)果一致.同時，還可以看到，圖7a中實線相連的菱形(caseB4_2、caseC4_2和caseD4_2)，點劃線相連的菱形(caseB4_4、caseC4_5和caseD4_5)以及虛線相連的菱形(caseB4_7、caseC4_8和caseD4_8)對應(yīng)的三條擬合直線幾乎平行，其斜率相差不大，平均值約為0.33.實線相連的方形(caseB4_1,caseC4_1和caseD4_1)，點劃線相連的方形(caseB4_3,caseC4_4和caseD4_4)以及虛線相連的方形(caseB4_6,caseC4_7和caseD4_7)對應(yīng)的三條擬合直線也幾乎平行，其斜率相差不大，平均值約為0.56.實線相連的圓圈(caseA4、caseB4、caseC4和caseD4)以及虛線相連的圓圈(caseB4_5、caseC4_6和caseD4_6)對應(yīng)的兩條擬合直線幾乎平行，其斜率相差不大，平均值約為0.76.這說明對相同的ηc，對流移除時間隨著增厚倍數(shù)的增大而增大，對流移除時間與增厚倍數(shù)基本成冪函數(shù)關(guān)系，B的影響不大(圖7a中的相同形狀的實線、虛線和點劃線的斜率大致相同)；但對流移除時間與增厚倍數(shù)的關(guān)系的冪指數(shù)隨著ηc的增大而減小(圖7a中圓圈、方形和菱形對應(yīng)的斜率，即θ，具體值見表4)，這說明對流移除時間與增厚倍數(shù)的關(guān)系與ηc相關(guān).根據(jù)表4的結(jié)果，可以得到對流移除時間與增厚倍數(shù)的關(guān)系的冪指數(shù)θ與黏性變化ηc的關(guān)系(圖7b).

圖7　(a)對流移除時間和巖石圈根增厚倍數(shù)γ的量化關(guān)系；(b)量化關(guān)系參數(shù)θ隨著組分黏性比ηc的變化關(guān)系；(c)去除了增厚倍數(shù)γ的影響后，對流移除時間隨著組分黏性比ηc的量化關(guān)系；(d)去除了增厚倍數(shù)γ和組分黏性比ηc的影響后，對流移除時間隨著Δρtc的量化關(guān)系；(e)去除增厚倍數(shù)γ，組分黏性比ηc以及Δρtc的影響后，對流移除時間隨著增厚寬度x的量化關(guān)系；(f) 模型計算結(jié)果和擬合結(jié)果的比較Fig.7　(a) Convective removal duration scale with the lithospheric root thickening factor γ; (b) Scaling law parameters θ varies with the composition viscosity ratio ηc; (c) Convective removal duration scale with the composition viscosity ratio ηc without the effect of the lithospheric root thickening factor γ; (d) Convective removal duration scale with the Δρtc without the effect of the lithospheric root thickening factor γ and the composition viscosity ratio ηc; (e) Convective removal duration scale with the width x without the effect of the lithospheric root thickening factor γ, the composition viscosity ratio ηc and the Δρtc; (f) Comparison of the model results and the fitting results

表4　γ的量化關(guān)系參數(shù)

(19)

即

(20)

由圖7a我們知道增厚倍數(shù)γ對于對流移除時間的影響很大，且當巖石圈的組分黏性比ηc不同時，影響的幅度也不同，所以在下面的分析數(shù)據(jù)中，我們都將增厚倍數(shù)對移除時間的影響利用(20)式去除后再來研究.去除增厚倍數(shù)的影響后，我們發(fā)現(xiàn)巖石圈的組分黏性比對移除時間的影響占主導地位(圖7c—7e).圖7c給出的是去除了增厚倍數(shù)γ的影響后，對流移除時間隨著組分黏性比ηc的對數(shù)關(guān)系圖，圖中的圓圈包含表3中所有的B≤0.1的模型和表2中的caseB1、caseB2、caseB3、caseB4、caseC4和caseD4.圖7c中包含18組數(shù)據(jù).每組增厚倍數(shù)、增厚區(qū)域的寬度以及等效密度變化都相同.從圖7c可以看到增厚倍數(shù)、增厚區(qū)域的寬度以及等效密度變化的影響幾乎難以區(qū)分.因此在獲取對流移除時間和ηc的冪指數(shù)時，我們采取將所有相同的ηc的模型的對流移除時間求平均，然后進行擬合，據(jù)此得到α=0.52，即

(21)

和3.1節(jié)一樣，在討論對流移除時間與等效密度變化Δρtc的關(guān)系時，也先將對流移除時間中包含的γ和ηc影響先行去除.圖7d給出了去除了增厚倍數(shù)γ以及組分黏性比ηc影響后，對流移除時間隨著等效密度變化Δρtc的變化規(guī)律，圖中的圓圈同樣包含了表3中所有的B≤0.1的模型和表2中的caseB1、caseB2、caseB3、caseB4、caseC4和caseD4.如圖7d所示，隨著熱邊界層總的等效密度變化的增大(浮力數(shù)的減小)，整體上的變化趨勢是相同的，都表現(xiàn)出了對流移除時間在減小，這里采用對所有數(shù)據(jù)點進行線性擬合的方法確定(18)式中的β.擬合結(jié)果給出β=-0.21，故

(22)

最后，在得到(20)—(22)式的基礎(chǔ)上，通過將對流移除時間進行歸算，即扣除增厚倍數(shù)γ、組分黏性比ηc以及等效密度變化Δρtc的影響后來研究它和增厚寬度x的關(guān)系，如圖7e所示.圖7e給出的是去除了增厚倍數(shù)γ、組分黏性比ηc以及等效密度變化Δρtc的影響后，對流移除時間隨著增厚區(qū)域?qū)挾葂的變化規(guī)律，圖中的圓圈同樣包含了表3中所有的γ=1.0且B≤0.1的模型和表2中的caseB1、caseB2、caseB3和caseB4.從圖7e中可以看到：利用(20)—(22)式歸算后，對固定的x，對流移除時間基本一致.因此，在獲取(18)式的μ值時，我們再一次采用將相同x時的數(shù)據(jù)取平均的方法.根據(jù)平均值數(shù)據(jù)擬合得到μ=0.04，即

(23)

μ值很小，說明了增厚寬度x對減薄時間的影響很小，這與一般巖石圈對流減薄關(guān)系(8)式不同.這可能主要是由于在此種情況下，側(cè)面對流侵蝕的影響顯著減小.綜合(20)—(23)式的分析，最后得到：

(24)

圖7f給出了根據(jù)(24)式計算的移除時間和數(shù)值模型計算得到的移除時間之間的比較，從圖7f中可以看到，(24)式可以很好地估算模型的計算結(jié)果.

式(24)可用來估算克拉通巖石圈被對流減薄的時間.采用3.1節(jié)相同的假定，即假定300 km(即γ=1.5)厚的克拉通巖石圈，在對流的作用下減薄到初始平衡厚度120 km.若克拉通巖石圈的密度比周圍地幔的密度低0.4%(即B=0.1或Δρtc=0.082)且組分黏性比ηc=10，在克拉通巖石圈寬度為300 km(x=0.3)時，其對流減薄時間大約為1.11 Ga；在克拉通巖石圈寬度為1500 km(x=1.5)時，其對流減薄時間增加到大約1.18 Ga.對比3.1節(jié)結(jié)果可以看到，0.4%的組分密度變化和10倍的組分黏性變化，使得巖石圈的移除時間從225～342 Ma增加到1.11～1.18 Ga.這顯示，克拉通巖石圈組分變化對移除時間具有重要影響.

保持其他參數(shù)不變，僅改變克拉通巖石圈組分浮力或黏性比，可以估算組分浮力或黏性比的影響.如仍假定300 km厚(即γ=1.5)、1500 km寬(x=1.5)、組分黏性比ηc=10，但若克拉通巖石圈的密度只比周圍地幔的密度低0.2%(即B=0.05或Δρtc=0.132)，對流減薄時間略微減小為1.07 Ga.這說明組分密度的變化影響不大.但如仍假定300 km厚(即γ=1.5)、1500 km寬(x=1.5)、克拉通巖石圈的密度比周圍地幔的密度低0.4%(即B=0.1或Δρtc=0.082)，而組分黏性比ηc=1，則依據(jù)(24)式，對流減薄時間約為426 Ma.這說明組分黏性的變化具有很大影響.

同樣也可估算巖石圈厚度變化的影響.如假定240 km(即γ=1.0)厚的克拉通巖石圈，在對流的作用下減薄到初始平衡厚度120 km，根據(jù)(24)式，在克拉通巖石圈寬度為1500 km(x=1.5)和克拉通巖石圈的密度比周圍地幔的密度低0.4%(即B=0.1或Δρtc=0.082)時，如組分黏性比ηc=10，克拉通巖石圈對流減薄時間大約是1.03 Ga.這顯示，克拉通巖石圈的厚度也具有一定影響.

早期的研究顯示，克拉通巖石圈根僅僅只有化學浮力并不能保持長期穩(wěn)定，其根的強度(即克拉通巖石圈根的黏性，本文指組分黏性比ηc)是其保持長期穩(wěn)定的關(guān)鍵因素(Lenardic and Moresi, 1999).最近，Wang等(2014)在研究克拉通的穩(wěn)定性時，也指出在組分黏性比ηc=3時，具有浮力的克拉通巖石圈根就能在對流的地幔中保持長達億年不受侵蝕；而當組分黏性比ηc增大到10，即使克拉通巖石圈具有的浮力更小抑或沒有本身固有的浮力，都能在地幔中長期的穩(wěn)定下來.這些都與式(24)的結(jié)果相一致.

根據(jù)目前的研究結(jié)果，華北克拉通巖石圈減薄發(fā)生在晚中生代，在早白堊世(120—130 Ma)達到高潮(吳福元等, 2008)，并且?guī)r石圈的厚度從200多公里減薄到80～120 km左右，說明華北克拉通的減薄持續(xù)時間約為140 Ma.為能應(yīng)用(24)式，這里假定華北克拉通從約200 km(即γ=0.67)厚經(jīng)過對流減薄后減薄到120 km，空間尺度為1500 km(即x=1.5).如果根據(jù)一般克拉通都具有高黏性、低密度的特點，假定華北克拉通巖石圈的組分黏性比ηc=10(華北克拉通巖石圈的黏性是周圍地幔黏性的10倍)，密度比周圍地幔的密度低0.4%(即B=0.1或Δρtc=0.082)，那么根據(jù)(24)式，對流減薄時間約897 Ma.這個時間明顯比140 Ma長.華北克拉通要在140 Ma得到減薄，必須發(fā)生交代置換或水化等作用使得密度或黏性產(chǎn)生變化.若僅發(fā)生密度變化，如交代作用使得克拉通巖石圈密度和周圍地幔的密度相同(即B=0或Δρtc=0.182)，則對流減薄時間約為760 Ma；即使外推到克拉通巖石圈的密度比周圍地幔高0.4%(即B=-0.1或Δρtc=0.282)，對流減薄時間也需約692 Ma.這說明如果化學交代等作用僅使克拉通巖石圈密度增大，就不足于在約140 Ma的時間內(nèi)使得華北克拉通巖石圈的厚度減薄到120 km.但若考慮交代、熔融或水化等作用可以減小克拉通巖石圈的組分黏性比時，華北克拉通就可以在200 Ma內(nèi)得到減薄.根據(jù)(24)式，即使假定華北克拉通巖石圈的密度仍比周圍地幔的密度低0.4%，當組分黏性比減小到ηc=3時，對流減薄時間減小到約446 Ma；當組分黏性比減小到ηc=1，對流減薄時間減小到約227 Ma；如果組分黏性比降到ηc=0.5，那么克拉通巖石圈對流減薄時間就只需145 Ma.這與目前研究得到的華北克拉通可能由于太平洋的俯沖，經(jīng)過脫水作用，華北克拉通巖石圈經(jīng)歷了強烈的水化，而大大降低了巖石圈的強度的結(jié)果一致(Xia et al., 2013;朱日祥等, 2012).

4結(jié)論

本文采用二維有限元數(shù)值模擬的方法，通過增厚(γ)一定寬度(x)內(nèi)的擬穩(wěn)狀態(tài)下的模型邊界層厚度來模擬加厚的巖石圈，然后考慮由化學組分引起的巖石圈的密度(B或Δρtc)和黏性變化(ηc)，定量研究了增厚的克拉通巖石圈在對流地幔中的移除時間，結(jié)果顯示：

(1) 當B=0和ηc=1，即忽略克拉通巖石圈固有的密度和黏性差異時，加厚巖石圈的移除時間基本上是隨著增厚巖石圈的寬度(x)增大而增大，同時也隨著增厚巖石圈的厚度(γ)增大而增大，并且通過計算分析，得到移除時間和增厚巖石圈的寬度以及厚度的量化關(guān)系t2D=0.0073γ0.70x0.26.這顯示，將一個初始厚度為120 km的巖石圈增厚到300 km，根據(jù)增厚的寬度(300 km或1500 km)的不同，其巖石圈根的對流移除時間在225 Ma或342 Ma不等.

(3) 對華北克拉通，如果其200 km厚的巖石圈是在140 Ma減薄到120 km厚，那么需要交代、熔融或水化作用減小其黏性.

致謝作者感謝李曙光老師以及Shijie Zhong教授對本研究提出的寶貴建議.

References

Boyd F R, Gurney J J, Richardson S H. 1985. Evidence for a 150～200 km thick Archaean lithosphere from diamond inclusion thermobarometry.Nature, 315(6018): 387-389.

Buck W R, Toks?z M N. 1983. Thermal effects of continental collisions: thickening a variable viscosity lithosphere.Tectonophysics, 100(1-3): 53-69.Carlson R W, Irving A J, Schulze D J, et al. 2004. Timing of Precambrian melt depletion and Phanerozoic refertilization events in the lithospheric mantle of the Wyoming craton and adjacent Central Plains Orogen.Lithos, 77(1-4): 453-472.

Carlson R W, Pearson D G, James D E. 2005. Physical, chemical, and chronological characteristics of continental mantle.Rev.Geophys., 43(1), doi: 10.1029/2004RG000156.

Davaille A, Jaupart C. 1994. Onset of thermal convection in fluids with temperature-dependent viscosity: application to the oceanic mantle.JournalofGeophysicalResearch, 99(B10): 19853-19866.

Davies G F. 1994. Thermomechanical erosion of the lithosphere by mantle plumes.JournalofGeophysicalResearch, 99(B8): 15709-15722.

Deng J F, Su S G, Zhao H L, et al. 2003. Deep processes of Mesozoic Yanshanian lithosphere thinning in North China.EarthScienceFrontiers(in Chinese), 10(3): 41-50.

Duggen S, Hoernle K, Van Den Bogaard P, et al. 2003. Deep roots of the Messinian salinity crisis.Nature, 422(6932): 602-606.

Fan W M, Menzies M A. 1992. Destruction of aged lower lithosphere and accretion of asthenosphere mantle beneath eastern China.GeotectonicaetMetallogenia(in Chinese), 16: 171-180.

Gao S, Rudnick R L, Carlson R W, et al. 2002. Re-Os evidence for replacement of ancient mantle lithosphere beneath the North China craton.EarthandPlanetaryScienceLetters, 198(3-4): 307-322.

Grasset O, Parmentier E M. 1998. Thermal convection in a volumetrically heated, infinite Prandtl number fluid with strongly temperature-dependent viscosity: Implications for planetary thermal evolution.JournalofGeophysicalResearch, 103(B8): 18171-18181.Griffin W L, Andi Z, O′reilly S Y, et al. 1998. Phanerozoic evolution of the lithosphere beneath the Sino-Korean Craton.∥ Mantle Dynamics and Plate Interactions in East Asia.AmericanGeophysicalUnion, 27: 107-126. Hager B H, Richards M A. 1989. Long-wavelength variations in earth′s geoid-physical models and dynamical implications.Philos.Trans.Roy.Soc.A, 328(1599): 309-327.

Houseman G A, Mckenzie D P, Molnar P. 1981. Convective instability of a thickened boundary layer and its relevance for the thermal evolution of continental convergent belts.JournalofGeophysicalResearch, 86(B7): 6115-6132.

Huang J S, Zhong S J, Van Hunen J. 2003. Controls on sublithospheric small-scale convection.JournalofGeophysicalResearch, 108(B8), doi: 10.1029/2003JB002456.

King S D. 2005. Archean cratons and mantle dynamics.EarthandPlanetaryScienceLetters, 234(1-2): 1-14.

Lenardic A, Moresi L N. 1999. Some thoughts on the stability of cratonic lithosphere: Effects of buoyancy and viscosity.JournalofGeophysicalResearch, 104(B6): 12747-12758.

Lenardic A, Moresi L N, Muhlhaus H. 2003. Longevity and stability of cratonic lithosphere: Insights from numerical simulations of coupled mantle convection and continental tectonics.JournalofGeophysicalResearch, 108(B6), doi: 10.1029/2002JB001859.

Li J H, Niu X L, Kusky T, et al. 2004. Neoarchean plate tectonic evolution of North China and its correlation with global cratonic blocks.EarthScienceFrontiers(in Chinese), 11(3): 273-283.

Liu D Y, Nutman A P, Compston W, et al. 1992. Remnants of ≥3800 Ma crust in the Chinese Part of the Sino-Korean craton.Geology, 20(4): 339-342.

Lu F X, Han Z G, Zheng J P, et al. 1991. Characteristics of Paleozoic mantle-lithosphere in Fuxian, Liaoning Province.GeologicalScienceandTechnologyInformation(in Chinese), 10(S1): 2-20.

Marotta A M, Fernàndez M, Sabadini R. 1999. The onset of extension during lithospheric shortening: a two-dimensional thermomechanical model for lithospheric unrooting.Geophys.J.Int., 139(1): 98-114.

Menzies M A, Fan W M, Zhang M. 1993. Palaeozoic and Cenozoic lithoprobes and the loss of >120 km of Archaean lithosphere, Sino-Korean craton, China.GeologicalSociety,London,SpecialPublications, 76(1): 71-81.Morency C, Doin M P, Dumoulin C. 2002. Convective destabilization of a thickened continental lithosphere.EarthandPlanetaryScienceLetters, 202(2): 303-320. Moresi L N, Solomatov V S. 1995. Numerical investigation of 2D convection with extremely large viscosity variations.Phys.Fluids, 7(9): 2154-2162.

Olson P, Schubert G, Anderson C, et al. 1988. Plume formation and lithosphere erosion: A comparison of laboratory and numerical experiments.JournalofGeophysicalResearch, 93(B12): 15065-15084.

O′Neill C J, Lenardic A, Griffin W L, et al. 2008. Dynamics of cratons in an evolving mantle.Lithos, 102(1-2): 12-24.

Pearson D G, Snyder G A, Shirey S B, et al. 1995. Archean Re-Os age for Siberian eclogites and constraints on Archaean tectonics.Nature, 374(6524): 711-713.

Pearson D G. 1999. The age of continental roots.Lithos, 48(1-4): 171-194.

Qiao Y C, Guo Z Q, Shi Y L. 2012. Numerical simulation of a possible thinning mechanism of the North China Craton-Gravitational instability delamination induced by lower crust eclogite.ChineseJ.Geophys. (in Chinese), 55(12): 4249-4256, doi: 10.6038/j.issn.0001-5733.2012.12.036.

Qiao Y C, Guo Z Q, Shi Y L. 2013. Thermal convection thinning of the North China Craton: Numerical simulation.ScienceChina:EarthSciences, 56(5): 773-782.

Rudnick R L, Fountain D M. 1995. Nature and composition of the continental crust: a Lower crustal perspective.Rev.Geophys., 33(3): 267-309. Schott B, Schmeling H. 1998. Delamination and detachment of a lithospheric root.Tectonophysics, 296(3-4): 225-247.

Schubert G, Turcotte D L, Olson P. 2001. Mantle Convection in the Earth and Planets. Cambridge, UK: Cambridge University Press, 940.

Schurr B, Rietbrock A, Asch G, et al. 2006. Evidence for lithospheric detachment in the central Andes from local earthquake tomography.Tectonophysics, 415(1-4): 203-223.

Shao J A, Zhang L Q, Li D M. 2002. Three Proterozoic extensional events in North China Craton.ActaPetrologicaSinica(in Chinese), 18(2): 152-160.

Shapiro S S, Hager B H, Jordan T H. 1999. Stability and dynamics of the continental tectosphere.Lithos, 48(1-4): 115-133.

Solomatov V S, Moresi L N. 2000. Scaling of time-dependent stagnant lid convection: Application to small-scale convection on Earth and other terrestrial planets.JournalofGeophysicalResearch, 105(B9): 21795-21817.

Turcotte D L, Schubert G. 2002. Geodynamics. Cambridge: Cambridge University Press.

Wang H L, Van Hunen J, Pearson D G, et al. 2014. Craton stability and longevity: The roles of composition-dependent rheology and buoyancy.EarthandPlanetaryScienceLetters, 391: 224-233.

Wu F Y, Ge W C, Sun D Y, et al. 2003. Discussions on the lithospheric thinning in eastern China.EarthScienceFrontiers(in Chinese), 10(3): 51-60.

Wu F Y, Lin J Q, Wilde S A, et al. 2005. Nature and significance of the Early Cretaceous giant igneous event in eastern China.EarthandPlanetaryScienceLetters, 233(1-2): 103-119.

Wu F Y, Xu Y G, Gao S, et al. 2008. Lithospheric thinning and destruction of the North China Craton.ActaPetrologicaSinica(in Chinese), 24(6): 1145-1174.

Xia Q K, Liu J, Liu S C, et al. 2013. High water content in Mesozoic primitive basalts of the North China Craton and implications on the destruction of cratonic mantle lithosphere.EarthandPlanetaryScienceLetters, 361: 85-97.

Xu W L, Wang Q H, Wang D Y, et al. 2004. Processes and mechanism of Mesozoic lithospheric thinning in eastern North China Craton: evidence from Mesozoic igneous rocks and deep-seated xenoliths.EarthScienceFrontiers(in Chinese), 11(3): 309-313.

Xu Y G. 2001. Thermo-tectonic destruction of the Archaean lithospheric keel beneath the Sino-Korean Craton in China: Evidence, timing and mechanism.Phys.Chem.EarthA, 26(9-10): 747-757.

Yoshida M. 2012. Dynamic role of the rheological contrast between cratonic and oceanic lithospheres in the longevity of cratonic lithosphere: A three-dimensional numerical study.Tectonophysics, 532-535: 156-166. Yuen D A, Fleitout L. 1985. Thinning of the lithosphere by small-scale convective destabilization.Nature, 313(5998): 125-128.

Zhao G C, Wilde S A, Cawood P A, et al. 2001. Archean blocks and their boundaries in the North China Craton: lithological, geochemical, structural and P-T path constraints and tectonic evolution.PrecambrianRes., 107(1-2): 45-73.

Zhao L, Zheng T Y. 2007. Complex upper-mantle deformation beneath the North China Craton: implications for lithospheric thinning.Geophys.J.Int., 170(3): 1095-1099.

Zhao Y, Xu G, Zhang S Z, et al. 2004. Yanshanian movement and conversion of tectonic regimes in East Asia.EarthScienceFrontiers(in Chinese), 11(3): 319-328.

Zhong S J, Zuber M T, Moresi L, et al. 2000. Role of temperature-dependent viscosity and surface plates in spherical shell models of mantle convection.JournalofGeophysicalResearch, 105(B5): 11063-11082.

Zhu R X, Chen L, Wu F Y, et al. 2011. Timing, scale and mechanism of the destruction of the North China Craton.ScienceChina:EarthSciences, 54(6): 789-797.

Zhu R X, Xu Y G, Zhu G, et al. 2012. Destruction of the North China Craton.ScienceChina:EarthSciences, 55(10): 1565-1587.

附中文參考文獻

鄧晉福, 蘇尚國, 趙海玲等. 2003. 華北地區(qū)燕山期巖石圈減薄的深部過程. 地學前緣, 10(3): 41-50.

范蔚茗, Menzies M A. 1992. 中國東部古老巖石圈下部的破壞和軟流圈地幔的增生. 大地構(gòu)造與成礦學, 16: 171-180.

李江海, 牛向龍, Kusky T等. 2004. 從全球?qū)Ρ忍接懭A北克拉通早期地質(zhì)演化與板塊構(gòu)造過程. 地學前緣, 11(3): 273-283.

路鳳香, 韓柱國, 鄭建平等. 1991. 遼寧復縣地區(qū)古生代巖石圈地幔特征. 地質(zhì)科技情報, 10(S1): 2-20.

喬彥超, 郭子祺, 石耀霖. 2012. 數(shù)值模擬華北克拉通巖石圈減薄的一種可能機制——下地殼榴輝巖重力失穩(wěn)引起的拆沉. 地球物理學報, 55(12): 4249-4256, doi: 10.6038/j.issn.0001-5733.2012.12.036.

喬彥超, 郭子祺, 石耀霖. 2013. 數(shù)值模擬華北克拉通巖石圈熱對流侵蝕減薄機制. 中國科學: 地球科學, 43(4): 642-652.

邵濟安, 張履橋, 李大明. 2002. 華北克拉通元古代的三次伸展事件. 巖石學報, 18(2): 152-160.

吳福元, 葛文春, 孫德有等. 2003. 中國東部巖石圈減薄研究中的幾個問題. 地學前緣, 10(3): 51-60.

吳福元, 徐義剛, 高山等. 2008. 華北巖石圈減薄與克拉通破壞研究的主要學術(shù)爭論. 巖石學報, 24(6): 1145-1174.

許文良, 王清海, 王冬艷等. 2004. 華北克拉通東部中生代巖石圈減薄的過程與機制: 中生代火成巖和深源捕虜體證據(jù). 地學前緣, 11(3): 309-313.

趙越, 徐剛, 張拴宏等. 2004. 燕山運動與東亞構(gòu)造體制的轉(zhuǎn)變. 地學前緣, 11(3): 319-328.

朱日祥, 陳凌, 吳福元等. 2011. 華北克拉通破壞的時間、范圍與機制. 中國科學: 地球科學, 41(5): 583-592.

朱日祥, 徐義剛, 朱光等. 2012. 華北克拉通破壞. 中國科學: 地球科學, 42(8): 1135-1159.

(本文編輯何燕)

Numerical simulation on convective thinning of the cratonic lithosphere

1LaboratoryofSeismologyandPhysicsofEarthInterior,SchoolofEarthandSpaceSciences,UniversityofScienceandTechnologyofChina,Hefei230026,China2MengchengNationalGeophysicalObservatory,Hefei230026,China

AbstractConvective removal on lithospheric root, especially on cratonic lithosphere root is studied via numerical simulations in 2D models. The control parameters in the models include the width, x, stretching factor, γ, viscosity ratio, ηc, and chemical buoyancy number, B, or effective density contrast, Δρtc. Numerical results show that mantle convection thins the thickened lithosphere, and (1) when B=0 and ηc=1, i.e., for general lithosphere, the root removal duration is scaled as 0.0073γ0.70x0.26, which means, for a 300 km thickness lithosphere with an initially equilibrium thickness of 120 km, and a root width of 300 km or 1500 km, it takes 225 Ma or 342 Ma to thin by sublithospheric convection to its equilibrium; (2) for small B and ηc, the process to thin cratonic lithosphere is similar to that to thin general lithosphere, but the root removal duration is increased significantly, and the root removal duration is scaled as With this scaling law, for a 300 km thickness lithosphere with an initially equilibrium thickness of 120 km, and a root width of 1500 km, the root removal duration is 1.18 Ga if ηc=10; and (3) for large B and ηc(B≥0.2 and ηc>10), the process to thin cratonic lithosphere is very different. Because of the influences of chemical buoyancy and viscosity, the cratonic lithosphere is spread to adjacent lithosphere instead of mixing with underneath mantle. In these cases, the root removal duration is very long (>3 Ga).

KeywordsCratonic lithosphere; Lithospheric thinning; Lithospheric root removal duration

基金項目國家自然科學基金項目(41474082,91014005)資助.

作者簡介程華冬，男，1986年生，安徽省肥西縣人，主要研究方向：地球動力學.E-mail:chengyd@mail.ustc.edu.cn *通訊作者黃金水，男，教授，主要研究方向：地球動力學.E-mail:jshhuang@ustc.edu.cn

doi:10.6038/cjg20160412 中圖分類號P315

收稿日期2015-01-27，2016-03-04收修定稿