依兰陨石坑形成过程数值模拟研究*

任健康，张庆明，刘文近，龙仁荣，龚自正，张品亮，宋光明，武 强，任思远

（1.北京理工大学 爆炸科学与技术国家重点实验室，北京 100081；2.北京卫星环境工程研究所，北京 100094）

小行星以超过10 km/s 的速度撞击地表会引起热辐射、地震和海啸等环境效应[1]。陨石坑是小行星撞击地球的直接证据和记录。根据陨石坑的形貌特征，一般将陨石坑分为简单坑、复杂坑和撞击盆地3 种[2]。目前，关于地球和其他类地行星陨石坑的形貌和形成条件的研究已有很多。例如：Saito 等[3]研究了小行星和彗星的尺寸、撞击速度以及角度对撞击陨石坑尺寸和形貌的影响；Ivanov[4-5]研究了简单坑和复杂坑的形成过程及机理，模拟了希克苏鲁伯陨石坑的形成机制，数值模拟结果与地质勘探数据吻合较好；Halim 等[6]复现了月球南极沙克尔顿陨石坑的形成条件及演化过程；Yue 等[7]基于数值模拟，研究了岫岩陨石坑熔化层的分布。缺乏对陨石坑底部熔化层质量分布规律以及撞击因素对成坑半径影响规律的研究。陈鸣等[8]在依兰陨石坑底部发现了冲击形成的花岗岩质撞击角砾岩，证实了依兰陨石坑是由小行星撞击形成的；Chen 等[9]研究了依兰陨石坑的地形地貌、地质特征和冲击变质特征。尚缺少对依兰陨石坑的成坑过程、形貌特征、底部熔化层分布规律的研究。

本文中，以依兰陨石坑为研究对象，针对上述研究存在的问题，利用iSALE-2D 仿真代码模拟依兰陨石坑成坑撞击条件、形成过程，研究其形成机理和坑底熔化层的分布规律。

1 依兰陨石坑

依兰陨石坑是中国第一大陨石坑，位于我国黑龙江省哈尔滨市的依兰县，地处小丘陵地区，该陨石坑坐落在白垩纪花岗岩基岩上。依兰陨石坑最终直径Df（坑缘直径，以坑缘隆起为基准测量的陨石坑的直径）为1 850 m，坑缘坑深dr（坑缘到坑底的深度）平均为260 m。根据坑中心东南方向处400 m 钻探的岩心（图1(b)中黄色矩形区域）材料分析结果，在坑底堆积的花岗岩质撞击角砾岩（图1(b)中红色矩形区域内存在熔化材料）中发现石英多组面状变形页理（planar deformation foliations, PDFs），提供了确凿的冲击变质证据，证实了该坑为小行星撞击地表形成的简单陨石坑[9]。依兰陨石坑全景照片和坑形示意图如图1 所示，图中Da为陨石坑的表观直径，dt为真实坑深，da为表观坑深。

图1 依兰陨石坑[8]Fig.1 Yilan crater[8]

2 形成过程

目前研究陨石坑形成机理主要有3 种方法：实验室尺度下小尺寸超高速撞击体撞击靶体实验研究、理论研究和数值模拟研究。由于实验成本高且未能反映真实撞击尺度、理论研究未能再现成坑动态过程等局限性，因此本文中采用数值模拟方法，探究依兰陨石坑的形成条件、过程以及熔化层分布规律。本文中，基于iSALE-2D 中的Euler 算法对依兰陨石坑展开数值模拟。iSALE-2D 是一个基于SALE[10]计算程序，支持Lagrange、Euler 和ALE 算法，包含多种适用于地质材料的强度模型和状态方程，广泛应用于撞击陨石坑模拟研究[6,11-12]。该版本程序只能计算二维轴对称模型，因此在本文讨论的范围内，撞击体均垂直撞击靶体，在将来授权的iSALE-3D 版本中可进一步讨论撞击角度对依兰陨石坑形成影响。

2.1 数值模型

球形小行星垂直撞击地面二维轴对称计算模型如图2 所示。靶体宽和高分别为3.1 和1.8 km。为了容纳撞击成坑过程产生的反溅抛射物，靶体上方建立宽3.1 km、高1.4 km 的真空区。模型采用Euler 网格计算。为了节约计算时间和成本，整个模型网格分为密分网格区和粗分网格区。Davison 等[13]已研究表明，当球形撞击体半径内超过20 个单元，撞击体半径内的网格数增加只会增加计算成本，对陨石坑形态影响很小。因此，细分网格区单元大小不变，固定为5 m，粗分网格区的单元尺寸以1.05 倍速度增长。靶体上部真空区在垂直和水平方向分别有60 和400 个细分单元网格。靶体细分网格区是由400 个水平单元和140 个垂直单元组成。在水平、垂直方向上，真空区、靶体的粗分网格区网格数均为50 个。同时，为了记录整个陨石坑内部冲击熔体的形成与演化过程，在模型的细分网格区布设Lagrange示踪点。

图2 模型工况网格参数Fig.2 Grid parameters of the model

模型上边界条件为out flow（允许物质流出边界）；下边界为no slip（垂直、平行于边界方向上材料的速度为零）；左边界为对称边界；右边界为free slip（垂直于边界方向上材料的速度为零）。

根据陨石坑成坑模型律计算[14]估计依兰陨石坑形成的撞击条件，可能是由直径约为100 m的小行星以15 km/s 的速度垂直撞击地表形成的。本文中共模拟8 组工况，具体参数如表1 所示，其中dp为撞击体直径，u为撞击速度，t为计算时间。

表1 模拟工况Table 1 Simulation conditions

2.2 重力初始化

为了更好地模拟地球陨石坑的形成过程，在整个模型上施加重力场，重力加速度为9.81 m/s2，方向与撞击速度同向。重力初始化在撞击之前完成，目的是计算靶体在不同深度的层压和密度值，并将计算结果作为初始状态赋给靶体材料。iSALE-2D 代码是从靶体表面开始计算岩石层压、靶体材料密度等参数，基于上层计算的结果，后续每个单元网格的密度均由对应岩石层压下的密度给出，得出靶体材料参数随着深度变化的分布，完成靶体重力场的初始化。在模拟撞击的初始时刻，模型加入重力场后，靶体材料参数重力初始化结果如图3 所示，图3(a)中靶体材料的岩石静压值随着靶体深度的增加而线性增大，靶体材料的密度值也相应初始化到对应靶体深度的状态，随靶体深度的增大而增大（见图3(b)）。至此，iSALE-2D 在初始时刻完成了1.8 km 高的靶体重力初始化。

图3 重力初始化靶体岩石静压和靶体密度Fig.3 Gravity-initialized lithostatic pressure and density of target

2.3 材料模型及参数

由于依兰陨石坑位于白垩纪花岗岩上，因此本文的模拟工况选择花岗岩作为靶体材料。撞击依兰陨石坑的小行星的材质是未知的，为计算简便，撞击体与靶体材料均选择花岗岩。小行星撞击地球速度在10 km/s 以上，陨石坑形成过程中的材料会发生熔化、气化[14]。ANEOS 状态方程是一种由计算机代码计算得到的列表形式的数据集合，它能较好地表征复杂材料在液相/气相转变、气相分子团和高压下的相变[15-16]，因此对撞击体和靶体材料均采用花岗岩ANEOS 状态方程[16]。

花岗岩材料的强度模型选用iSALE 中的ROCK 强度模型[17]，强度随材料累积应变而降低，ROCK 强度模型将材料屈服强度Y定义为：

数值模拟过程中材料损伤的表征对撞击过程中靶体强度描述起关键作用，由于材料能够恢复弹性应变，因此损伤的度量主要由塑性应变来确定，当物质完全损伤时将不能承受剪切力，此时定义物质损伤度D为1[18]。本文中采用Ivanov 损伤模型[19]，损伤（D取值范围为[0,1]）是塑性应变的函数：

式中：εp为总塑性应变的累积量，εf为损伤时的塑性应变，εfb为低压状态下的最小失效应变，pc为压缩失效应力，B为大于0 的常数。花岗岩的Ivanov 损伤模型具体参数：εfb=10-4,B=10-11Pa-1,pc=300.0 MPa。

小行星撞击地表时部分花岗岩材料的温度会接近其熔化温度，当材料温度达到熔化温度时剪切强度下降到零，因此需要在模型计算过程中动态调整高温区材料的剪切强度和内摩擦力，更好地反映成坑过程中坑壁的坍塌和坑底的物质运动。本文中使用Ohnaka 热软化模型[20]，利用脆性到脆-塑性转变状态下的剪切破坏强度关系式来近似这种高温区材料的剪切强度和内摩擦力动态调整行为：

式中：Y为环境温度下的材料屈服强度（包括损伤的影响），Yt为热软化强度，T为环境温度，Tm为材料的熔化温度，ξ 为材料常数，Tm0为材料在零压下的熔化温度，a0和c为材料常数。花岗岩的Ohnaka 热软化模型具体参数：Tm0=1 673.0 K,a0=6.0 GPa,c=3.0, ξ=1.2。

2.4 成坑过程分析

8 组工况成坑过程相似，因此，以直径120 m、速度12 km/s 的花岗岩质小行星垂直撞击地表为例，介绍简单坑的形成过程。简单坑形成过程主要分为3 个阶段：接触与压缩阶段、开坑阶段和后期调整阶段[2]。为了便于观察撞击过程中材料的运动和坑形演化过程，后处理中通过Python 脚本在模型对称轴右侧100 m 处单独显示间隔50 m 的6 个绿色Lagrange 示踪点，并更改模型细分网格区设置的示踪点显示方式，如图4 所示材料上的示踪网格（并非本文中使用的Euler 算法网格）。

(1)接触和压缩阶段：如图4(a)所示，撞击体接触靶体表面开始向下侵彻，在撞击区域形成高温高压区，温度远超过1 500 K（见图4(b)），陨石坑半径和坑深逐渐增大，直到撞击体本身完全卸载，形成1～2 倍撞击体直径的陨石坑，接触压缩阶段结束。这一阶段持续时间约为40 ms。

(2)开坑阶段：靶体中的材料在剩余速度的作用下促使撞击坑不断扩大，形成反溅的抛射物，抛射物的抛射角度从90°（抛射物垂直靶面）逐渐减小到60°（见图4(c)～(d)）。当陨石坑底部停止增长，形成一个碗状的最大瞬时坑，开坑阶段结束。如图5 所示，陨石坑在t=4 s 时坑深达到最大值为440 m，此时陨石坑表观直径达到1 170 m，坑内温度达到1 000 K。

(3)后期调整阶段：在开坑阶段结束后，陨石坑内的温度依然很高，陨石坑壁面的材料由于热软化导致材料剪切强度很低，受到重力发生滑落，汇聚在坑底，导致坑底不断升高形成隆起（见图4(e)），升高到一定程度后，底部隆起又在重力作用下崩塌，如此循环（见图5 中3 个虚线框区）。循环过程中陨石坑内的温度降低，材料的剪切强度不断增大，陨石坑坑内隆起高度逐渐降低。如图5 陨石坑坑深的时程曲线表明，在t=4～150 s 时，陨石坑坑深呈阶梯状减小，在t=20, 80, 140 s 时，陨石坑底部出现隆起。当t＞150 s 时，陨石坑坑深基本保持不变，陨石坑成坑过程基本完成（见图4(f)）。图4 中6 个Lagrange 示踪点运动状态表明坑底堆积物在后期调整过程中发生堆叠。

图4 直径120 m、速度12 km/s 的花岗岩质小行星撞击成坑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 的花岗岩质小行星撞击坑的半径和坑深随时间的演化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 陨石坑熔化层分布

陨石坑成坑过程中，撞击点处为高温高压区，在此区域内的部分撞击体与靶体材料会发生相变，即熔化和气化。熔化材料的形成和分布规律可以作为陨石坑形成机理和勘探的参考依据，而考虑材料气化会增加计算成本，因此在程序计算过程中材料密度低于5 kg/m3时会自动删除，本文中不考虑材料的气化。

3.1 熔化判据

花岗岩材料受到冲击后，从初始状态跳跃到Hugoniot 曲线上对应压力下的状态，后续由于稀疏波沿着等熵线卸载，当材料的峰值压力超过临界压力值时，才会发生熔化。材料从固态沿着等熵线卸载，等熵线与初始熔化相线相交时，材料开始发生初始熔化；材料压力继续卸载，若等熵线与完全熔化相线相交时，材料则完全熔化。因此，依据靶体的峰值压力判断材料是否发熔化或完全熔化，将花岗岩的初始和完全熔化的压力阈值（46、56 GPa）作为材料熔化判据[7,16]写入iSALE-2D 后处理程序中，材料在整个撞击时间历程中峰值压力超过给定的完全熔化压力阈值（56 GPa）时，就认为该材料发生完全熔化，并在数据处理过程中以红色示踪点表示其分布位置，未发生完全熔化的不作标记。

3.2 熔化层分布规律

图6 是工况7 条件下陨石坑熔化层的形成过程与分布图。

图6 直径120 m、速度12 km/s 的花岗岩质小行星撞击陨石坑形成过程花岗岩熔化层分布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 阶段形成的。熔体连续地分布在坑内表面（见图6(c)），小部分熔体随抛射物离开陨石坑（见图6(d)）。由于重力作用，大部分熔体随着坑壁的坍塌沉积陨石坑底部（见图6(e)），后期调整阶段陨石坑底部不断调整，底部出现熔体堆叠（见图6(f)）。

8 组工况完全熔化的花岗岩质量统计和坑直径方向400 m 处岩心熔化层深度（dm）分布如表2 所示，工况1～3、5 未出现完全熔化的花岗岩，其余工况均出现，深度均小于等于50 m。撞击体与靶体材料完全熔化的质量（mt）由撞击体质量（mp）无量纲化，形成无量纲参数mt/mp。在12、15 km/s 撞击速度条件下，mt/mp随着撞击体直径的增大而增大，在15 km/s 撞击速度条件下其值约为12 km/s 的1.7 倍（见图7）。

图7 无量纲化的花岗岩完全熔化质量与撞击体直径的关系Fig.7 Dimensionless completely-melted granite mass varied with projectile diameter

表2 陨石坑熔化层分布Table 2 Distribution of melted layers in craters

依兰坑勘验数据为坑直径方向400 m 处湖相沉积物底下108～127 m 处出现花岗岩的冲击熔化产物，模拟结果均比勘探数据小，可能是由于后期坑壁岩石的坍塌在钻孔处的堆积，造成了此处模拟的熔化层深度与勘探深度的差异。根据熔化层分布区间19 m 可知，工况4、6～8 吻合较好。

4 成坑模型

4.1 最终陨石坑形态

图8 是8 组工况撞击体撞击地表形成陨石坑轮廓剖面对比图。从图8 可以看出，不同撞击条件下成坑形貌基本一致。成坑直径和深度随着撞击体直径的增大和撞击速度的提高而增大。依兰陨石坑的勘探数据[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 节中讨论的熔化层分布可知，最接近依兰陨石坑形成条件的为工况7。

表3 8 组模拟工况在150 s 坑形数据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和δ 分别是小行星半径、速度和密度，µ 是与材料相关的常数。µ 的理论值在1/3 和2/3 之间。对密实材料冲击实验结果表明µ 的值在0.55～0.6 之间，非常高孔隙率材料的 µ ≈0.4[22]。ν 的值与材料特性无关，不同材料ν 的取值均为0.4[23]。靶体的密度为 ρ、抗拉强度Y、黏度为 η ，重力加速度为g，撞击体与靶体剩下的材料参数与物质常数均由材料抗拉强度、黏度、质量密度3 个基本物质常数无量纲化，记作 Πm。假设小行星垂直撞击地表，成坑半径R可由下式表示：

撞击体与靶体材料保持不变，可忽略常数项 Πm，许多与成坑有关的地质材料都被认定是与率无关的，因此材料没有黏度单元物质常数 η。根据点源假设，成坑半径R对撞击条件和材料特性的依赖关系：

根据陨石坑成坑半径的增长是在靶体强度作用下停止还是在靶体重力作用下停止，可以将陨石坑成坑机理分为强度机制和重力机制。对于大型陨石坑或者强度非常低的材料，成坑半径的增长是在重力的作用下停止，靶体强度的影响可以忽略。在重力机制控制下的成坑半径相似关系：

式中：m为小行星质量，m=4πδa3/3 ；H1为常数，Y为靶体强度。

在强度机制下，成坑增长是在靶体强度的作用下停止，重力的影响可以忽略。在强度机制控制下的成坑半径相似关系：

在本文研究的依兰陨石坑撞击尺度条件下，成坑半径相似关系式中强度项大于重力项，即Y/(ρu2)＞ga/u2，根据式(9)～(11)的推导可知，依兰陨石坑的形成属于强度机制：由上述点源成坑相似律分析可知，依兰陨石坑半径与撞击速度的 µ 次方、靶体材料强度的 - µ/2 次方成正比。因本文的撞击体与靶体密度相同，故无量纲化的密度项为1。根据本文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 数据的对比Fig.9 Comparison of the crater radius fitting curve and the data of condition 7

5 结 论

采用iSALE-2D 仿真代码对依兰陨石坑的形成条件和过程进行模拟。对不同撞击条件下的成坑尺寸、撞击过程中的熔化材料的分布进行了统计分析。模拟结果与依兰陨石坑勘探结果比对表明，依兰陨石坑很有可能是由直径为120 m 的球状花岗岩小行星以12 km/s 的垂直撞击地表形成的，形成一个最终直径1 840 m、坑缘坑深为263 m 的简单陨石坑。

依据数值模拟结果复现了依兰陨石坑动态成坑过程，得出撞击过程中完全熔化的材料的分布规律，其主要分布在陨石坑底，呈堆叠状分布，少量的离散分布在靶体表面。统计得出模拟工况7 产生的熔化层质量约为5.64×1010kg，约为撞击体质量的24 倍。

在本文的研究范围内，由点源成坑相似律模型分析得出强度机制下依兰陨石坑的成坑半径的无量纲关系式，ΠR=0.943 77Π1-0.22257，模拟工况7 计算坑直径结果与坑直径拟合关系式结果相对误差10.3%。

对Collins G、Wünnemann K、Elbeshausen D 等iSALE-2D 的开发者，以及软件后处理工具pySALEPlot 的开发者Davison T 等谨表感谢。