數值模擬華北克拉通巖石圈減薄的一種可能機制——下地殼榴輝巖重力失穩引起的拆沉

喬彥超，郭子祺，石耀霖

1 中國科學院遙感應用研究所，北京 100101

2 中國科學院研究生院，北京 100039

1 引 言

大陸巖石圈一直被認為是相對穩定的構造單元，然而并不是所有的前寒武紀克拉通從它們形成以后都是處于穩定狀態.近期的研究發現，如北美的Wyoming克拉通自中元古代以來發生了減薄和破壞［1－2］，南非的 Kapvaal克拉通在3.1Ga 時由于巖石圈地幔的拆沉作用而被改造［3］，Seber等研究發現地中海Alboran海盆及其周圍摩洛哥北部Rif造山帶和西班牙南部Betic造山帶（兩者均屬阿爾卑斯造山期產物）巖石圈地幔現今正在發生拆沉作用［4］，很多結果也表明［5－18］華北克拉通在古生代再活化而發生減薄，克拉通巖石圈的改造和破壞已經成為大陸動力學研究的一個熱點.也正是如此，華北克拉通巖石圈破壞已成為繼青藏高原和大別—蘇魯超高壓變質帶之后，我國又一在國際地球動力學界引起廣泛關注的重大科學問題，已成為國際地球科學前沿研究領域.

關于華北克拉通巖石圈減薄機制還是有比較多的爭議［6－8，12，17－20］，高山等［17－18］通過研究在華北克拉通發現高鎂安山巖、英安巖、和埃達克巖，并認為這些巖石是部分熔融的榴輝巖和地幔巖石相互反應生成，是華北克拉通拆沉的主要證據.拆沉作用泛指由于重力的不穩定性導致巖石圈地幔、大陸下地殼或大洋地殼沉入下伏軟流圈或地幔的過程.俯沖大洋巖石圈中榴輝巖的拆沉再循環是板塊構造的直接產物［21－23］.

拆沉理論認為巖石圈加厚是大陸巖石圈和下地殼拆沉的前奏.巖石圈加厚促使基性下地殼轉變為榴輝巖［17－18］；榴輝巖的密度比巖石圈地幔的密度高200～400kg／m3［24－26］.這種差異將造成重力上的不穩定性，使榴輝巖拆沉再循環進入地幔中［27－32］.

本文研究通過數值計算的方法，結合華北克拉通已有的地質成果，模擬華北克拉通巖石圈拆沉機制，通過研究討論能夠發生拆沉需要榴輝巖的厚度，不同規模榴輝巖造成拆沉的不同模式，以及對比熱侵蝕和拆沉所造成的地表不同反應，來揭示華北克拉通巖石圈不同減薄機制的動力所對應的不同地質構造特征.

2 模 型

2.1 數值計算的基本原理

我們曾對計算原理、使用公式以及所使用程序的驗證進行過論述［23］，本文不再詳述.研究基于粘滯系數、密度與溫度相關的二維上地幔小尺度地幔對流［33－40］，其關系式為

其中η0是等效粘滯系數，b是經驗系數，Ttop和Tbottom分別是計算區域上下邊界的溫度，ρ0是等效密度，α是膨脹系數.

我們使用了有限差分方法，并結合MIC技術［34－36］求解連續性方程，動量守恒方程和能量守恒方程.計算中使用無限大Prandtl數假定，這樣忽略動量方程中的慣性項，假定流體是不可壓縮的，并且引入 Boussinesq假設［41，42］.

計算的二維模型是長700km，寬700km的方盒子，網格為71×71，共有400×400個示蹤點.邊界條件：所有的邊都為自由邊界條件；左右兩個邊界是絕熱邊界條件 （?T／?x ＝0）；頂邊界溫度為Ttop＝273K；底邊界溫度為Tbottom＝1773K.

克拉通拆沉過程中的第一步是巖石圈加厚玄武巖相變為密度比較大的榴輝巖.計算中我們使用由Wylllie［43］提出，Anderson［44］修改的玄武巖－榴輝巖相變圖，如圖1.因此大于50km的地殼被認為是榴輝巖質的，小于25km是玄武質，深度在25～50km之間且溫度達到Teclo（z）（公式3，圖1）的玄武巖轉變為榴輝巖.

圖1 巖石圈加厚玄武巖相變為榴輝巖相變圖Fig.1 Phase transition of basalt to eclogite when lithospheres thicken

當溫度達到固相線時巖石開始融化，計算溫度，選擇 McKenzie［45］公式（圖1）：

本文使用的是溫度場定義的巖石圈厚度.巖石圈加厚前的地溫曲線為圖1中黑線，上地殼地溫梯度為25℃／km，則20km處為500℃；巖石圈100km厚，溫度為1350℃［46］；最底邊200km 處溫度為1400℃，如圖1中Told.地殼厚度為35km，即圖中原來Moho的深度.假定巖石圈快速擠壓增厚一倍到200km深，那么新的地溫曲線為圖1中綠色曲線，上地殼變為40km，溫度為500℃；巖石圈變為200km，溫度為1350℃，如圖1中Tnew.Tnew與Telo的交點深度約為34km，加厚的Moho深度為70km，那么從34km到70km都為相變區域，我們首先考慮極端情形，假定這一區域都變為榴輝巖.

2.2 模型初始溫度場和物質組成

我們計算區域為700km×700km，溫度場主要參考Tnew的值，在底邊界700km處溫度為1500℃，如圖2.關于物質組成，我們考慮極端的情形，數值模擬中我們往往對現實的模型或者資料進行很多簡化或者合理的假設，從而專注于我們研究的問題，并且數值模擬的一大優勢是，在保證正確的前提下，我們可以進行，不同參數，不同條件，不同模型的討論.使我們得到研究問題的全面資料.本研究中我們假定初始地殼厚度為70km，榴輝巖厚度為34～70km，首先計算橫向榴輝巖規模200～500km，如圖3.我們將討論橫向榴輝巖不同規模，如300～400km，200～500km，100～600km，即水平尺度分別為100km，200km，300km時所造成不同的拆沉模式.

圖2 計算的初始溫度場Fig.2 Initial thermal state

圖3 初始物質組成Fig.3 Initial density

2.3 參數

計算中使用的參數如表1.計算中關鍵的參數是不同組分密度的選取，榴輝巖的密度取3500kg／m3［47］，地殼的密度取2800kg／m3，地幔的密度取3340kg／m3［48］.

表1 計算中用到的參數Table 1 Parameters used in computation

3 計算結果

榴輝巖的厚度為34～70km，我們橫向上采用三種不同規模的榴輝巖分別為200～500km，300～400km，100～600km，討論了不同規模的榴輝巖造成的不同拆沉結果.其溫度場等值線不同時刻結果如圖4.

圖4中第一橫排為中等規模榴輝巖（200～500km）情況下，不同時刻溫度場拆沉的等溫線結果，依次第二橫排為較小規模榴輝巖（300～400km），第三橫排為較大規模榴輝巖（100～600km）.通過對比這三種情況，我們可以發現，榴輝巖的規模越大拆沉的越劇烈，最終能發生拆沉的區域也越大.為了更好的觀察拆沉的過程，同時我們給出了不同模型密度的隨時間的變化過程如圖5.

圖5的橫排次序如圖4.從圖5中可見，當榴輝巖規模較大時，拆沉開始階段的過程更傾向于雙管道向下拆沉，模型較短時為單管道拆沉的模式.為了得到更精確的定量數據，我們給出了三種模型的榴輝巖最終時刻的溫度彩圖疊加速度矢量圖.其結果如圖6—8.

圖6中三條黑色線從上到下分別為500，1000，1350℃時的溫度等值線圖，綠色線為拆沉前的巖石圈底邊界，即200km處.因此從綠色線到最下面的黑線之間的區域為華北克拉通巖石圈能夠發生拆沉的區域.通過計算我們能夠知道減薄后的巖石圈最大厚度為64km，發生拆沉的最大厚度為136km.從流體場可以看出最終時刻發育兩個比較完整的對流環.整個流體場比較穩定.

圖7中的線代表的意義如圖6.通過計算我們能夠知道減薄后的巖石圈最大厚度為95km，發生拆沉的最大厚度為105km.可以看到最終時刻1350℃等溫線以上兩端局部有小的對流環，中心區域則是比較穩定，這樣兩端可能會熱流更高，對應著比較強烈，時間比較長的巖漿活動.同華北克拉通巖漿分布地質資料［51］相比也比較符合.

圖8中的線代表的意義如圖6.通過計算我們能夠知道減薄后的巖石圈最大厚度為108km，發生拆沉的最大厚度為92km.可以看到最終時刻流體場發育多管道向下拆沉.

圖4 三種模型不同時刻溫度場等值線圖Fig.4 Temperature contours of three models at different times

圖5 三種模型不同時刻密度場變化圖Fig.5 Density variations of three models at different times

圖6 較長模型時最終時刻拆沉的結果Fig.6 Delamination result of the big model at last

圖7 模型為中等時最終時刻拆沉的結果Fig.7 Delamination result of the middle model at last

圖8 模型為最小時最終時刻拆沉的結果Fig.8 Delamination result of the small model at last

4 拆沉減薄與熱侵蝕減薄的對比討論

通過對華北克拉通巖石圈熱對流減薄機制的計算［32］及本文對華北克拉通巖石圈拆沉減薄機制的研究對比，我們可以知道在力學機制上這兩種機制都是可能發生的，即華北克拉通巖石圈由于熱擾動引起對流減薄和由于榴輝巖密度不穩定引起對流減薄都是可能發生的.但是通過細致的對比如圖9，我們可以發現不同的減薄機制的結果是不一致的，對應的地表地質特征也完全不同.

圖9 熱對流減薄和拆沉減薄結果對比圖Fig.9 Comparison of heat convection thinning and delamination

從圖9中左圖為熱對流減薄的結果，右圖是拆沉減薄的結果，從圖中藍色箭頭標志我們可以發現，熱對流減薄，中心區域變熱變薄，中心區域是一個伸展區域，發展一個大的變質核雜巖區，巖漿的發展時序是從中心向兩側對稱分布，并且向兩側減弱；熱拆沉減薄結果中，中心區域受到擠壓加厚，反而變冷，但是我們知道，對于地質歷史上發生過的加熱事件通過地質手段，如探針等方法是可以得到信息的，對于變冷的過程我們是沒有辦法得到信息的.在中心區域的兩側發生兩個巖漿活動中心，巖漿的發展時序是從兩側向中心減弱發展.

通過上面討論，何種減薄機制在這一過程中起主導作用，要更多的分析地表的地質資料，來找到足夠的證據.

5 結 論

本文研究表明，中生代華北拉通巖石圈受擠壓，玄武巖相變為密度較大的榴輝巖，榴輝巖重力失穩能夠引起下地殼的拆沉，造成巖石圈的大規模減薄.對比不同規模的榴輝巖減薄方式發現，當榴輝巖較大規模時更可能發生雙管道拆沉，而規模小時發生單管道拆沉.榴輝巖模型最大時候（100～600km）減薄后的巖石圈最大厚度為64km，發生拆沉的最大厚度為136km；規模中等時候（200～500km）巖石圈最大厚度為95km，發生拆沉的最大厚度為105km；規模比較小時候（300～400km）巖石圈最大厚度為108km，發生拆沉的最大厚度為92km.

通過對比中生代華北克拉通巖石圈熱減薄和拆沉減薄機制可知，從力學機制上都能造成巖石圈減薄，但他們的結果對應著不同的地質特征.對比我們的計算，熱侵蝕最多能減薄100km［32］，拆沉最多能夠減薄136km，但是很多資料都表明華北東部最多減薄了120km以上，所以單從減薄量上來看拆沉減薄機制更符合實際地質資料.根據華北克拉通巖漿分布規律［7，49－50］我們可以發現有東北部和南部兩個比較大的中心，這與我們文中討論的中等規模的榴輝巖拆沉減薄機制能夠在中心兩側發展兩個巖漿活動中心比較一致.同時如果我們單純考慮熱侵蝕減薄機制，它所對應的底部熱擾動太過于偶然性，為什么不在其他地方，而只在華北克拉通，相比較于熱侵蝕，由于拆沉理論得到了很多地質學的證據，如高山等的地球化學證據證明華北發生拆沉進入地幔的榴輝巖存在，所以我們更傾向于榴輝巖重力失穩引發拆沉是主導華北克拉通巖石圈減薄的機制.

（References）

［1］Sleep N H.Evolution of the continental lithosphere.Annual Review of Earth and Planetary Sciences，2005，33：369－393.

［2］Carlson R W，Pearson D G，Physical，chemical and chronological characteristics of continental mantle.Reviews of Geophysics，2005，43：RG1001，doi：10.1029／2004RG000156.

［3］Moser D E，Flowers R M，Hart R J.2001.Birth of the Kaapvaal tectosphere 3.08billion years ago.Science，291：465－468.

［4］Seber D，Barazangi M，lbebrahim A.Geophysical evidence for lithospheric delamination beneath the Alboran Sea and Rif－Betic mountains.Nature，1996，379：785－790.

［5］Chen L，Zheng T，Xu W.A thinned lithospheric image of the Tanlu Fault Zone，eastern China：Constructed from wave equation based receiver function migration.Journal of Geophysical Research，2006，111：B09312.

［6］高山，章軍鋒，許文良，等.拆沉作用與華北克拉通破壞.科學通報，2009，54：1962－1973.Gao S，Zhang J F，Xu W L，et al.Delamination and destruction of the North China Craton.Chinese Sci Bull，2009，54，doi：10.1007／s11434－009－0395－9.

［7］吳福元，徐義剛，高山，等.華北巖石圈減薄與克拉通破壞研究的主要學術爭論.巖石學報，2008，24：1145－1174.Wu F Y，Xu，Y G，Gao S，et al.Controversy over studies of the lithospheric thinning and craton destruction of North China.Acta Petrol Sin （in Chinese），2008，24：1145－1174.

［8］郝曉光，方劍，劉根友等.巖石圈拆沉還是軟流圈上涌？——正常密度理論與華北克拉通破壞重力學機制研究.地球物理學進展，2011，26（4）：1137－1140，doi：10.3969／j.issn.1004－2903.2011.04.001.Hao X G ，Fang J，Liu G Y ，et al.North China craton destruction，delamination or thermal erosion？Chinese J.Geophys.（in Chinese），2011，26（4）：1137－1140，doi：10.3969／j.issn.1004－2903.2011.04.001 .

［9］許衛衛，鄭天愉，趙亮.華北地區410km間斷面和660km間斷面結構—克拉通活化的地幔動力學狀態探測.2011，中國科學：地球科學，41：678–685.Xu W W，Zheng T Y，Zhao L.Mantle dynamics of the reactivating North China Craton：Constraints from the topographies of the 410km and 660km discontinuities.Sci China Earth Sci，2011，54：881－887，doi：10.1007／s11430－010－4163－0.

［10］張旗.從"巖石圈減薄"到"克拉通破壞".地球物理學進展，2011，26（6）：2262－2269，doi：10.3969／j.issn.1004－2903.2011.06.047.Zhang Q.The evolution from"lithospheric thinning"to"cratonic destruction" .Chinese J.Geophys.（in Chinese），2011，26（6）：2262－2269，doi：10.3969／j.issn.1004－2903.2011.06.047.

［11］Fan W M，Zhang H F，Baker J.On and off the north China craton：Where is the Archaean keel？Journal of Petrology，2000，41：933－950.

［12］Xu Y G. Thermo－tectonic destruction of the Archean lithospheric keel beneath the Sino－Korean Craton in China：Evidence，timing and mechanism.Physics and Chemisty of the Earth，2001，26：747－757.

［13］Liu Y S，Gao S.Melt－peridotite interactions：Links between garnet pyroxenite and high－Mg signature of continental crust.Earth and Planetary Science Letters，2005，234：39－57.

［14］Zheng J P，Griffin W L，O′Reilly S Y.Late Mesozoic－Eocene mantle replacement beneath the eastern North China Craton：Evidence from the Paleozoic and Cenozoic peridotite xenoliths.International Geology Review，2005，47：457－472.

［15］Zheng J P，Sun M，Zhou M F.Trace elemental and PGE geochemical constraints of Mesozoic and Cenozoic peridotitic xenoliths on lithospheric evolution of the North China Craton.Geochimica et Cosmochimica Acta，2005，69：3401－3418.

［16］Zheng J P，Griffin W L，O′Reilly S Y.Mineral chemistry of garnet peridotites from Paleozoic，Mesozoic and Cenozoic lithosphere：Constraints on mantle evolution beneath eastern China.Journal of Petrology，2006，47：2233－2256.

［17］Gao S，Rudnick R L，Carlson R W.Re－Os evidence for replacement of ancient mantle lithosphere beneath the North China craton.Earth and Planetary Science Letters，2002，198：307－322.

［18］Gao S，Rudnick R L，Yuan H L.Recycling lower continental crust in the North China craton.Nature，2004，432：892－897.

［19］Zhang H F，Sun M，Zhou X H.Mesozoic lithosphere destruction beneath the North China Craton：Evidence from major，trace element，and Sr－Nd－Pb isotope studies of Fangcheng basalts.Contributions to Mineralogy and Petrology，2002，144：241－253.

［20］Zhang H F，Sun M，Zhou X H.Secular evolution of the lithosphere beneath the eastern North China Craton：Evidence from Mesozoic basalts and high－Mg andesites.Geochimica et Cosmochimica Acta，2003，67：4373－4387.

［21］Sobolev A V，Hofmann A W，Sobolev S V，et al.An olivine－free mantle source of Hawaiian shield basalts.Nature，2005，434：590－597.

［22］Sobolev A V，Hofmann A W，Kuzmin D V，et al.The amount of recycled crust in sources of mantle－derived melts.Science，2007，316：412－417.

［23］郭慧麗，徐佩芬.地震層析成像在華北克拉通地區的研究進展.地球物理學進展，2011，26（5）：1557－1565，doi：10.3969／j.issn.1004－2903.2011.05.007.Guo H L，Xu P F.Progress of Seismic tomography applied in the North China Craton.Chinese J.Geophys.（in Chinese），2011，26（5）：1557－1565，doi：10.3969／j.issn.1004－2903.2011.05.007.

［24］Jull M，Kelemen.On the conditions for lower crustal convective instability.Journal of Geophysical Research，2001，106，6423－6446.

［25］Anderson D A.Large igneous provinces，delamination，and fertile mantle.Elements，2006，1：271－275.

［26］Anderson D A.Speculations on the nature and cause of mantle heterogeneity.Tectonophysics，2006，146：7－22.

［27］Lustrino M.How the delamination and detachment of lower crust can influence basaltic magmatism. Earth－Science Review，2005，72：21－38.

［28］Bédard J H. A catalytic delamination－driven model for coupled genesis of Archaean crust and sub－continental lithospheric mantle.Geochimica et Cosmochimica Acta，2006，70：1188－1214.

［29］Arndt N T，Goldstein S L.An open boundary between lower continental crust and mantle：Its role in crust formation and crustal recycling.Tectonophysics，1989，161：201－212.

［30］Kay R W，Kay S M.Creation and destruction of lower continental crust.Geologische Rundschau，1991，80：259－278.

［31］Rudnick R L.Making continental crust.Nature，1995.378：571－577.

［32］喬彥超，郭子祺，石耀霖.數值模擬華北克拉通巖石圈熱對流侵蝕減薄機制.中國科學：地球科學，2012，（待發表）.Qiao Y C，Guo Z Q，Shi Y L.Thermal Convection Thinning of the North China Craton：Numerical Simulation.Sci China Earth Sci.（in press）.

［33］Huang J S，Zhong S，van Hunen J.Controls on sublithospheric small－scale convection.Journal of Geophysical Research，2003，108（B8），2405.Doi：10.1029／2003JB002456.

［34］Gerya T V，Yuen D A.Characteristics－based marker－in－cell method with conservative finite－differences schemes for modeling geological flows with strongly variable transport properties.Physics of the Earth and Planetary Interiors，2003，140：295－320.

［35］Gerya T V， Maresch W V， Willner A P.Inherent gravitational instability of thickened continental crust with regionally developed low－timedium－pressure granulite facies metamorphism.Earth and Planetary Science Letters，2001，190：221－235.

［36］Gerya T V，Perchuk L L，Maresch W V.Thermal regime and gravitational instability of multi－layered continental crust：implications for the buoyant exhumation of high－grade metamorphic rocks.European Journal of Mineralogy，2002，14：687－699.

［37］葉正仁，王建.上地幔變黏度小尺度對流的數值研究.地球物理學報，2003，46（3）：335－339.Ye Z R，Wang J.A numerical research on the small scale convection with variable viscosity in the upper mantle.Chinese J.Geophys.（in Chinese），2003，46（3）：335－339.

［38］傅容珊，常筱華，黃建華，劉文忠.區域重力均衡異常和上地幔小尺度對流模型.地球物理學報，1994，37（S2）：249－258.Fu R S，Chang X H，Huang J H，et al.Regional gravity isostatic anomaly and small scale convection model in upper mantle.Chinese J.Geophys.（in Chinese），1994，37（S2）：249－258.

［39］熊熊，傅容珊，許厚澤等.增厚大陸巖石層熱邊界層對流剝離的數值模擬.地球物理學報，1998，41（S1）：33－40.Xiong X，Fu R S，Xu H Z，et al.Numerical Simulation of the convective removal of the thickened thermal boundary layer of continental lothosphere.Chinese J.Geophys.（in Chinese），1998，41（S1）：33－40.

［40］傅容珊，黃建華，董樹謙.利用地震層析成象數據計算地幔對流新模型的探討.地球物理學報，2003，46（6）：772－778.Fu R S.A new mantle convection model constrained by seismic tomography data.Chinese J.Geophys.（in Chinese），2003，46（6）：772－778.

［41］Boussinesq J.Theorie analytique de la chaleur mise en harmonie avec la thermodynamique et avec la theorie mecanique de la lumiere.Gauthier－Villars Paris，1903，2：157－176.

［42］Tackley P. Effect of strongly temperature－dependant viscosity on time－dependent 3－dimensional model of mantle convection.Journal of Geophysical Research，1993，20（20）：2187－2190.

［43］Wyllie P.The Dynamic Earth.New York：Wiley，1971：416.

［44］Anderson D L. Tectonics and composition of Venus.Geophys.Res.Lett.，1980，1：101－102.

［45］McKenzie D.The generation and compaction of partially molten rock.Journal of Petrology，1984，25（3）：713－765.

［46］McKenzie D，Roberts J M，Weiss N O.Convection in the earth＇s mantle：Towards a numerical simulation.Journal of Fluid Mechanics，1974，6：465－538.

［47］Anderson D L.The upper mantle transition region：eclogite？Geophys.Res.Lett.，1979，6：433－435.

［48］Dupeyrat L.Thermal and chemical convection in planetary mantles.Journal of Geophysical Research，1995，100：497－520.

［49］吳福元，李獻華，楊進輝等.花崗巖成因研究的若干問題.巖石學報，2007，23（4）：1217－1238.Wu F Y，Li X H，Yang J H，et al.Discussions on the petrogenesis ofgranites.Acta Petrol.Sin.（in Chinese），2007，23（4）：1217－1238.

［50］嵇少丞，王茜，許志琴.華北克拉通破壞與巖石圈減薄.地質學報，2008，82（2）：174－193.Ji S C，Wang Q，Xu Z Q.Break－up of the North China Craton through lithospheric thinning.Acta Geologica Sinica（in Chinese），2008，82（2）：174－193.致 謝 感謝鄭亮先生提供并維護并行機器.