依蘭隕石坑形成過程數(shù)值模擬研究*

任健康，張慶明，劉文近，龍仁榮，龔自正，張品亮，宋光明，武 強，任思遠

（1.北京理工大學 爆炸科學與技術國家重點實驗室，北京 100081；2.北京衛(wèi)星環(huán)境工程研究所，北京 100094）

小行星以超過10 km/s 的速度撞擊地表會引起熱輻射、地震和海嘯等環(huán)境效應[1]。隕石坑是小行星撞擊地球的直接證據(jù)和記錄。根據(jù)隕石坑的形貌特征，一般將隕石坑分為簡單坑、復雜坑和撞擊盆地3 種[2]。目前，關于地球和其他類地行星隕石坑的形貌和形成條件的研究已有很多。例如：Saito 等[3]研究了小行星和彗星的尺寸、撞擊速度以及角度對撞擊隕石坑尺寸和形貌的影響；Ivanov[4-5]研究了簡單坑和復雜坑的形成過程及機理，模擬了?？颂K魯伯隕石坑的形成機制，數(shù)值模擬結果與地質(zhì)勘探數(shù)據(jù)吻合較好；Halim 等[6]復現(xiàn)了月球南極沙克爾頓隕石坑的形成條件及演化過程；Yue 等[7]基于數(shù)值模擬，研究了岫巖隕石坑熔化層的分布。缺乏對隕石坑底部熔化層質(zhì)量分布規(guī)律以及撞擊因素對成坑半徑影響規(guī)律的研究。陳鳴等[8]在依蘭隕石坑底部發(fā)現(xiàn)了沖擊形成的花崗巖質(zhì)撞擊角礫巖，證實了依蘭隕石坑是由小行星撞擊形成的；Chen 等[9]研究了依蘭隕石坑的地形地貌、地質(zhì)特征和沖擊變質(zhì)特征。尚缺少對依蘭隕石坑的成坑過程、形貌特征、底部熔化層分布規(guī)律的研究。

本文中，以依蘭隕石坑為研究對象，針對上述研究存在的問題，利用iSALE-2D 仿真代碼模擬依蘭隕石坑成坑撞擊條件、形成過程，研究其形成機理和坑底熔化層的分布規(guī)律。

1 依蘭隕石坑

依蘭隕石坑是中國第一大隕石坑，位于我國黑龍江省哈爾濱市的依蘭縣，地處小丘陵地區(qū)，該隕石坑坐落在白堊紀花崗巖基巖上。依蘭隕石坑最終直徑Df（坑緣直徑，以坑緣隆起為基準測量的隕石坑的直徑）為1 850 m，坑緣坑深dr（坑緣到坑底的深度）平均為260 m。根據(jù)坑中心東南方向處400 m 鉆探的巖心（圖1(b)中黃色矩形區(qū)域）材料分析結果，在坑底堆積的花崗巖質(zhì)撞擊角礫巖（圖1(b)中紅色矩形區(qū)域內(nèi)存在熔化材料）中發(fā)現(xiàn)石英多組面狀變形頁理（planar deformation foliations, PDFs），提供了確鑿的沖擊變質(zhì)證據(jù)，證實了該坑為小行星撞擊地表形成的簡單隕石坑[9]。依蘭隕石坑全景照片和坑形示意圖如圖1 所示，圖中Da為隕石坑的表觀直徑，dt為真實坑深，da為表觀坑深。

圖1 依蘭隕石坑[8]Fig.1 Yilan crater[8]

2 形成過程

目前研究隕石坑形成機理主要有3 種方法：實驗室尺度下小尺寸超高速撞擊體撞擊靶體實驗研究、理論研究和數(shù)值模擬研究。由于實驗成本高且未能反映真實撞擊尺度、理論研究未能再現(xiàn)成坑動態(tài)過程等局限性，因此本文中采用數(shù)值模擬方法，探究依蘭隕石坑的形成條件、過程以及熔化層分布規(guī)律。本文中，基于iSALE-2D 中的Euler 算法對依蘭隕石坑展開數(shù)值模擬。iSALE-2D 是一個基于SALE[10]計算程序，支持Lagrange、Euler 和ALE 算法，包含多種適用于地質(zhì)材料的強度模型和狀態(tài)方程，廣泛應用于撞擊隕石坑模擬研究[6,11-12]。該版本程序只能計算二維軸對稱模型，因此在本文討論的范圍內(nèi)，撞擊體均垂直撞擊靶體，在將來授權的iSALE-3D 版本中可進一步討論撞擊角度對依蘭隕石坑形成影響。

2.1 數(shù)值模型

球形小行星垂直撞擊地面二維軸對稱計算模型如圖2 所示。靶體寬和高分別為3.1 和1.8 km。為了容納撞擊成坑過程產(chǎn)生的反濺拋射物，靶體上方建立寬3.1 km、高1.4 km 的真空區(qū)。模型采用Euler 網(wǎng)格計算。為了節(jié)約計算時間和成本，整個模型網(wǎng)格分為密分網(wǎng)格區(qū)和粗分網(wǎng)格區(qū)。Davison 等[13]已研究表明，當球形撞擊體半徑內(nèi)超過20 個單元，撞擊體半徑內(nèi)的網(wǎng)格數(shù)增加只會增加計算成本，對隕石坑形態(tài)影響很小。因此，細分網(wǎng)格區(qū)單元大小不變，固定為5 m，粗分網(wǎng)格區(qū)的單元尺寸以1.05 倍速度增長。靶體上部真空區(qū)在垂直和水平方向分別有60 和400 個細分單元網(wǎng)格。靶體細分網(wǎng)格區(qū)是由400 個水平單元和140 個垂直單元組成。在水平、垂直方向上，真空區(qū)、靶體的粗分網(wǎng)格區(qū)網(wǎng)格數(shù)均為50 個。同時，為了記錄整個隕石坑內(nèi)部沖擊熔體的形成與演化過程，在模型的細分網(wǎng)格區(qū)布設Lagrange示蹤點。

圖2 模型工況網(wǎng)格參數(shù)Fig.2 Grid parameters of the model

模型上邊界條件為out flow（允許物質(zhì)流出邊界）；下邊界為no slip（垂直、平行于邊界方向上材料的速度為零）；左邊界為對稱邊界；右邊界為free slip（垂直于邊界方向上材料的速度為零）。

根據(jù)隕石坑成坑模型律計算[14]估計依蘭隕石坑形成的撞擊條件，可能是由直徑約為100 m的小行星以15 km/s 的速度垂直撞擊地表形成的。本文中共模擬8 組工況，具體參數(shù)如表1 所示，其中dp為撞擊體直徑，u為撞擊速度，t為計算時間。

表1 模擬工況Table 1 Simulation conditions

2.2 重力初始化

為了更好地模擬地球隕石坑的形成過程，在整個模型上施加重力場，重力加速度為9.81 m/s2，方向與撞擊速度同向。重力初始化在撞擊之前完成，目的是計算靶體在不同深度的層壓和密度值，并將計算結果作為初始狀態(tài)賦給靶體材料。iSALE-2D 代碼是從靶體表面開始計算巖石層壓、靶體材料密度等參數(shù)，基于上層計算的結果，后續(xù)每個單元網(wǎng)格的密度均由對應巖石層壓下的密度給出，得出靶體材料參數(shù)隨著深度變化的分布，完成靶體重力場的初始化。在模擬撞擊的初始時刻，模型加入重力場后，靶體材料參數(shù)重力初始化結果如圖3 所示，圖3(a)中靶體材料的巖石靜壓值隨著靶體深度的增加而線性增大，靶體材料的密度值也相應初始化到對應靶體深度的狀態(tài)，隨靶體深度的增大而增大（見圖3(b)）。至此，iSALE-2D 在初始時刻完成了1.8 km 高的靶體重力初始化。

圖3 重力初始化靶體巖石靜壓和靶體密度Fig.3 Gravity-initialized lithostatic pressure and density of target

2.3 材料模型及參數(shù)

由于依蘭隕石坑位于白堊紀花崗巖上，因此本文的模擬工況選擇花崗巖作為靶體材料。撞擊依蘭隕石坑的小行星的材質(zhì)是未知的，為計算簡便，撞擊體與靶體材料均選擇花崗巖。小行星撞擊地球速度在10 km/s 以上，隕石坑形成過程中的材料會發(fā)生熔化、氣化[14]。ANEOS 狀態(tài)方程是一種由計算機代碼計算得到的列表形式的數(shù)據(jù)集合，它能較好地表征復雜材料在液相/氣相轉(zhuǎn)變、氣相分子團和高壓下的相變[15-16]，因此對撞擊體和靶體材料均采用花崗巖ANEOS 狀態(tài)方程[16]。

花崗巖材料的強度模型選用iSALE 中的ROCK 強度模型[17]，強度隨材料累積應變而降低，ROCK 強度模型將材料屈服強度Y定義為：

數(shù)值模擬過程中材料損傷的表征對撞擊過程中靶體強度描述起關鍵作用，由于材料能夠恢復彈性應變，因此損傷的度量主要由塑性應變來確定，當物質(zhì)完全損傷時將不能承受剪切力，此時定義物質(zhì)損傷度D為1[18]。本文中采用Ivanov 損傷模型[19]，損傷（D取值范圍為[0,1]）是塑性應變的函數(shù)：

式中：εp為總塑性應變的累積量，εf為損傷時的塑性應變，εfb為低壓狀態(tài)下的最小失效應變，pc為壓縮失效應力，B為大于0 的常數(shù)。花崗巖的Ivanov 損傷模型具體參數(shù)：εfb=10-4,B=10-11Pa-1,pc=300.0 MPa。

小行星撞擊地表時部分花崗巖材料的溫度會接近其熔化溫度，當材料溫度達到熔化溫度時剪切強度下降到零，因此需要在模型計算過程中動態(tài)調(diào)整高溫區(qū)材料的剪切強度和內(nèi)摩擦力，更好地反映成坑過程中坑壁的坍塌和坑底的物質(zhì)運動。本文中使用Ohnaka 熱軟化模型[20]，利用脆性到脆-塑性轉(zhuǎn)變狀態(tài)下的剪切破壞強度關系式來近似這種高溫區(qū)材料的剪切強度和內(nèi)摩擦力動態(tài)調(diào)整行為：

式中：Y為環(huán)境溫度下的材料屈服強度（包括損傷的影響），Yt為熱軟化強度，T為環(huán)境溫度，Tm為材料的熔化溫度，ξ 為材料常數(shù)，Tm0為材料在零壓下的熔化溫度，a0和c為材料常數(shù)?；◢弾r的Ohnaka 熱軟化模型具體參數(shù)：Tm0=1 673.0 K,a0=6.0 GPa,c=3.0, ξ=1.2。

2.4 成坑過程分析

8 組工況成坑過程相似，因此，以直徑120 m、速度12 km/s 的花崗巖質(zhì)小行星垂直撞擊地表為例，介紹簡單坑的形成過程。簡單坑形成過程主要分為3 個階段：接觸與壓縮階段、開坑階段和后期調(diào)整階段[2]。為了便于觀察撞擊過程中材料的運動和坑形演化過程，后處理中通過Python 腳本在模型對稱軸右側(cè)100 m 處單獨顯示間隔50 m 的6 個綠色Lagrange 示蹤點，并更改模型細分網(wǎng)格區(qū)設置的示蹤點顯示方式，如圖4 所示材料上的示蹤網(wǎng)格（并非本文中使用的Euler 算法網(wǎng)格）。

(1)接觸和壓縮階段：如圖4(a)所示，撞擊體接觸靶體表面開始向下侵徹，在撞擊區(qū)域形成高溫高壓區(qū)，溫度遠超過1 500 K（見圖4(b)），隕石坑半徑和坑深逐漸增大，直到撞擊體本身完全卸載，形成1～2 倍撞擊體直徑的隕石坑，接觸壓縮階段結束。這一階段持續(xù)時間約為40 ms。

(2)開坑階段：靶體中的材料在剩余速度的作用下促使撞擊坑不斷擴大，形成反濺的拋射物，拋射物的拋射角度從90°（拋射物垂直靶面）逐漸減小到60°（見圖4(c)～(d)）。當隕石坑底部停止增長，形成一個碗狀的最大瞬時坑，開坑階段結束。如圖5 所示，隕石坑在t=4 s 時坑深達到最大值為440 m，此時隕石坑表觀直徑達到1 170 m，坑內(nèi)溫度達到1 000 K。

(3)后期調(diào)整階段：在開坑階段結束后，隕石坑內(nèi)的溫度依然很高，隕石坑壁面的材料由于熱軟化導致材料剪切強度很低，受到重力發(fā)生滑落，匯聚在坑底，導致坑底不斷升高形成隆起（見圖4(e)），升高到一定程度后，底部隆起又在重力作用下崩塌，如此循環(huán)（見圖5 中3 個虛線框區(qū)）。循環(huán)過程中隕石坑內(nèi)的溫度降低，材料的剪切強度不斷增大，隕石坑坑內(nèi)隆起高度逐漸降低。如圖5 隕石坑坑深的時程曲線表明，在t=4～150 s 時，隕石坑坑深呈階梯狀減小，在t=20, 80, 140 s 時，隕石坑底部出現(xiàn)隆起。當t＞150 s 時，隕石坑坑深基本保持不變，隕石坑成坑過程基本完成（見圖4(f)）。圖4 中6 個Lagrange 示蹤點運動狀態(tài)表明坑底堆積物在后期調(diào)整過程中發(fā)生堆疊。

圖4 直徑120 m、速度12 km/s 的花崗巖質(zhì)小行星撞擊成坑3 個階段材料和溫度分布Fig.4 Material and temperature distributions in three stages of crater formation induced by the impact of a 120-m-diameter granitic asteroid with the velocity of 12 km/s

圖5 直徑120 m、速度12 km/s 的花崗巖質(zhì)小行星撞擊坑的半徑和坑深隨時間的演化Fig.5 Time histories of the radius and depth of the crater induced by the impact of a 120-m-diameter granite asteroid with the velocity of 12 km/s

3 隕石坑熔化層分布

隕石坑成坑過程中，撞擊點處為高溫高壓區(qū)，在此區(qū)域內(nèi)的部分撞擊體與靶體材料會發(fā)生相變，即熔化和氣化。熔化材料的形成和分布規(guī)律可以作為隕石坑形成機理和勘探的參考依據(jù)，而考慮材料氣化會增加計算成本，因此在程序計算過程中材料密度低于5 kg/m3時會自動刪除，本文中不考慮材料的氣化。

3.1 熔化判據(jù)

花崗巖材料受到?jīng)_擊后，從初始狀態(tài)跳躍到Hugoniot 曲線上對應壓力下的狀態(tài)，后續(xù)由于稀疏波沿著等熵線卸載，當材料的峰值壓力超過臨界壓力值時，才會發(fā)生熔化。材料從固態(tài)沿著等熵線卸載，等熵線與初始熔化相線相交時，材料開始發(fā)生初始熔化；材料壓力繼續(xù)卸載，若等熵線與完全熔化相線相交時，材料則完全熔化。因此，依據(jù)靶體的峰值壓力判斷材料是否發(fā)熔化或完全熔化，將花崗巖的初始和完全熔化的壓力閾值（46、56 GPa）作為材料熔化判據(jù)[7,16]寫入iSALE-2D 后處理程序中，材料在整個撞擊時間歷程中峰值壓力超過給定的完全熔化壓力閾值（56 GPa）時，就認為該材料發(fā)生完全熔化，并在數(shù)據(jù)處理過程中以紅色示蹤點表示其分布位置，未發(fā)生完全熔化的不作標記。

3.2 熔化層分布規(guī)律

圖6 是工況7 條件下隕石坑熔化層的形成過程與分布圖。

圖6 直徑120 m、速度12 km/s 的花崗巖質(zhì)小行星撞擊隕石坑形成過程花崗巖熔化層分布Fig.6 Distribution of granite melt layers during the crater formation induced by the impact of a 120-m-diameter granite asteroid with the velocity of 12 km/s

如圖6(a)～(b)所示，花崗巖熔體主要是在成坑的第1 階段形成的。熔體連續(xù)地分布在坑內(nèi)表面（見圖6(c)），小部分熔體隨拋射物離開隕石坑（見圖6(d)）。由于重力作用，大部分熔體隨著坑壁的坍塌沉積隕石坑底部（見圖6(e)），后期調(diào)整階段隕石坑底部不斷調(diào)整，底部出現(xiàn)熔體堆疊（見圖6(f)）。

8 組工況完全熔化的花崗巖質(zhì)量統(tǒng)計和坑直徑方向400 m 處巖心熔化層深度（dm）分布如表2 所示，工況1～3、5 未出現(xiàn)完全熔化的花崗巖，其余工況均出現(xiàn)，深度均小于等于50 m。撞擊體與靶體材料完全熔化的質(zhì)量（mt）由撞擊體質(zhì)量（mp）無量綱化，形成無量綱參數(shù)mt/mp。在12、15 km/s 撞擊速度條件下，mt/mp隨著撞擊體直徑的增大而增大，在15 km/s 撞擊速度條件下其值約為12 km/s 的1.7 倍（見圖7）。

圖7 無量綱化的花崗巖完全熔化質(zhì)量與撞擊體直徑的關系Fig.7 Dimensionless completely-melted granite mass varied with projectile diameter

表2 隕石坑熔化層分布Table 2 Distribution of melted layers in craters

依蘭坑勘驗數(shù)據(jù)為坑直徑方向400 m 處湖相沉積物底下108～127 m 處出現(xiàn)花崗巖的沖擊熔化產(chǎn)物，模擬結果均比勘探數(shù)據(jù)小，可能是由于后期坑壁巖石的坍塌在鉆孔處的堆積，造成了此處模擬的熔化層深度與勘探深度的差異。根據(jù)熔化層分布區(qū)間19 m 可知，工況4、6～8 吻合較好。

4 成坑模型

4.1 最終隕石坑形態(tài)

圖8 是8 組工況撞擊體撞擊地表形成隕石坑輪廓剖面對比圖。從圖8 可以看出，不同撞擊條件下成坑形貌基本一致。成坑直徑和深度隨著撞擊體直徑的增大和撞擊速度的提高而增大。依蘭隕石坑的勘探數(shù)據(jù)[9]：最終撞擊坑直徑為1.85 km，坑緣坑深約為260 m。表3 中得出撞擊體直徑90、100 m 的4 個工況最終撞擊坑直徑均小于1.85 km，撞擊體直徑110、120 m 的4 個工況最終撞擊坑直徑分別為1 770、1 930、1 840 和2 040 m，坑緣坑深分別為251、317、263 和359 m。120 m、12 km/s 的工況成坑的最終直徑為1 840 m，坑緣坑深為263 m，由此可得該工況撞擊條件下形成的隕石坑和依蘭隕石坑坑形最接近，結合第3 節(jié)中討論的熔化層分布可知，最接近依蘭隕石坑形成條件的為工況7。

表3 8 組模擬工況在150 s 坑形數(shù)據(jù)Table 3 Crater shape data for eight simulation conditions at 150 s

圖8 不同模擬工況在計算結束時間150 s 所對應的隕石坑剖面輪廓Fig.8 The crater profiles corresponding to the different simulation conditions at the end time of 150 s

4.2 點源成坑相似律模型

點源相似律模型廣泛用于將實驗室尺度的沖擊結果外延至行星級別的撞擊結果[21-23]。點源假設認為，當碰撞效應范圍遠大于小行星本身時，可以將小行星等效一個點源[21]。小行星對撞擊結果的影響僅由一個變量決定C=auμδν，而非獨立的3 個a、u和δ 自變量[21]。a、u和δ 分別是小行星半徑、速度和密度，μ 是與材料相關的常數(shù)。μ 的理論值在1/3 和2/3 之間。對密實材料沖擊實驗結果表明μ 的值在0.55～0.6 之間，非常高孔隙率材料的 μ ≈0.4[22]。ν 的值與材料特性無關，不同材料ν 的取值均為0.4[23]。靶體的密度為 ρ、抗拉強度Y、黏度為 η ，重力加速度為g，撞擊體與靶體剩下的材料參數(shù)與物質(zhì)常數(shù)均由材料抗拉強度、黏度、質(zhì)量密度3 個基本物質(zhì)常數(shù)無量綱化，記作 Πm。假設小行星垂直撞擊地表，成坑半徑R可由下式表示：

撞擊體與靶體材料保持不變，可忽略常數(shù)項 Πm，許多與成坑有關的地質(zhì)材料都被認定是與率無關的，因此材料沒有黏度單元物質(zhì)常數(shù) η。根據(jù)點源假設，成坑半徑R對撞擊條件和材料特性的依賴關系：

根據(jù)隕石坑成坑半徑的增長是在靶體強度作用下停止還是在靶體重力作用下停止，可以將隕石坑成坑機理分為強度機制和重力機制。對于大型隕石坑或者強度非常低的材料，成坑半徑的增長是在重力的作用下停止，靶體強度的影響可以忽略。在重力機制控制下的成坑半徑相似關系：

式中：m為小行星質(zhì)量，m=4πδa3/3 ；H1為常數(shù)，Y為靶體強度。

在強度機制下，成坑增長是在靶體強度的作用下停止，重力的影響可以忽略。在強度機制控制下的成坑半徑相似關系：

在本文研究的依蘭隕石坑撞擊尺度條件下，成坑半徑相似關系式中強度項大于重力項，即Y/(ρu2)＞ga/u2，根據(jù)式(9)～(11)的推導可知，依蘭隕石坑的形成屬于強度機制：由上述點源成坑相似律分析可知，依蘭隕石坑半徑與撞擊速度的 μ 次方、靶體材料強度的 - μ/2 次方成正比。因本文的撞擊體與靶體密度相同，故無量綱化的密度項為1。根據(jù)本文8 組工況的計算結果，以及額外追加計算的4 組工況（直徑110 m，撞擊速度10、11、13 和14 km/s）的模擬結果來擬合式(11)，可得H2=0.943 77，- μ/2 =-0.222 57，將工況7 代入擬合曲線中，可得隕石坑最終坑直徑為2030 m，與該工況的模擬結果相對誤差為10.3%，工況7 和擬合曲線對比如圖9 所示。

圖9 成坑半徑擬合曲線與工況7 數(shù)據(jù)的對比Fig.9 Comparison of the crater radius fitting curve and the data of condition 7

5 結 論

采用iSALE-2D 仿真代碼對依蘭隕石坑的形成條件和過程進行模擬。對不同撞擊條件下的成坑尺寸、撞擊過程中的熔化材料的分布進行了統(tǒng)計分析。模擬結果與依蘭隕石坑勘探結果比對表明，依蘭隕石坑很有可能是由直徑為120 m 的球狀花崗巖小行星以12 km/s 的垂直撞擊地表形成的，形成一個最終直徑1 840 m、坑緣坑深為263 m 的簡單隕石坑。

依據(jù)數(shù)值模擬結果復現(xiàn)了依蘭隕石坑動態(tài)成坑過程，得出撞擊過程中完全熔化的材料的分布規(guī)律，其主要分布在隕石坑底，呈堆疊狀分布，少量的離散分布在靶體表面。統(tǒng)計得出模擬工況7 產(chǎn)生的熔化層質(zhì)量約為5.64×1010kg，約為撞擊體質(zhì)量的24 倍。

在本文的研究范圍內(nèi)，由點源成坑相似律模型分析得出強度機制下依蘭隕石坑的成坑半徑的無量綱關系式，ΠR=0.943 77Π1-0.22257，模擬工況7 計算坑直徑結果與坑直徑擬合關系式結果相對誤差10.3%。

對Collins G、Wünnemann K、Elbeshausen D 等iSALE-2D 的開發(fā)者，以及軟件后處理工具pySALEPlot 的開發(fā)者Davison T 等謹表感謝。