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

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

1.   依兰陨石坑

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

图  1  依兰陨石坑^[8]

Figure  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  模型工况网格参数

Figure  2.  Grid parameters of the model

下载: 全尺寸图片 幻灯片

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

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

表  1  模拟工况

Table  1.  Simulation conditions

工况	dp /m	u/(km∙s−1)	t/s
1	90	12	150
2	90	15	150
3	100	12	150
4	100	15	150
5	110	12	150
6	110	15	150
7	120	12	150
8	120	15	150

下载: 导出CSV 

| 显示表格

2.2   重力初始化

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

图  3  重力初始化靶体岩石静压和靶体密度

Figure  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为损伤系数，${Y}_{\mathrm{d}}$为损伤材料的强度，${Y}_{\mathrm{i}}$为未损伤材料的强度。损伤材料强度由Drucker-Prager关系^[17]定义：

式中：${Y}_{\mathrm{d}0}$为零压下的材料屈服强度，损伤材料的内聚力，${\mu }_{\mathrm{d}}$为损伤材料的内摩擦因数，p为压力，${Y}_{\mathrm{d}\mathrm{m}}$为高压下损伤材料强度极限。未受损伤的材料强度是由Lundborg关系^[17]定义的：

式中：${Y}_{\mathrm{i}0}$为未损伤材料的内聚力，${\mu }_{\mathrm{i}}$为未损伤材料的内摩擦因数，${Y}_{\mathrm{i}\mathrm{m}}$为未损伤材料高压强度极限。花岗岩的ROCK强度模型具体参数：Y_i0=10.0 MPa, μ_i=2.0, Y_im=2.5 GPa, Y_d0=10.0 kPa, μ_d=0.6, Y_dm=2.5 GPa。

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

式中：${\varepsilon }_{\mathrm{p}}$为总塑性应变的累积量，${\varepsilon }_{\mathrm{f}}$为损伤时的塑性应变，${\varepsilon }_{\mathrm{f}\mathrm{b}}$为低压状态下的最小失效应变，${p}_{\mathrm{c}}$为压缩失效应力，B为大于0的常数。花岗岩的Ivanov损伤模型具体参数：ε_fb=10⁻⁴, B=10⁻¹¹ Pa⁻¹, p_c=300.0 MPa。

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

式中：Y为环境温度下的材料屈服强度（包括损伤的影响），${Y}_{\mathrm{t}}$为热软化强度，T为环境温度，${T}_{\mathrm{m}}$为材料的熔化温度，$\xi$为材料常数，${T}_{\mathrm{m}0}$为材料在零压下的熔化温度，${a}_{0}$和$c$为材料常数。花岗岩的Ohnaka热软化模型具体参数：T_m0=1 673.0 K, a₀=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算法网格）。

图  4  直径120 m、速度12 km/s的花岗岩质小行星撞击成坑3个阶段材料和温度分布

Figure  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

下载: 全尺寸图片 幻灯片

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

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

图  5  直径120 m、速度12 km/s的花岗岩质小行星撞击坑的半径和坑深随时间的演化

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

3.   陨石坑熔化层分布

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

3.1   熔化判据

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

3.2   熔化层分布规律

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

图  6  直径120 m、速度12 km/s的花岗岩质小行星撞击陨石坑形成过程花岗岩熔化层分布

Figure  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处岩心熔化层深度（d_m）分布如表2所示，工况1～3、5未出现完全熔化的花岗岩，其余工况均出现，深度均小于等于50 m。撞击体与靶体材料完全熔化的质量（m_t）由撞击体质量（m_p）无量纲化，形成无量纲参数m_t/m_p。在12、15 km/s撞击速度条件下，m_t/m_p随着撞击体直径的增大而增大，在15 km/s撞击速度条件下其值约为12 km/s的1.7倍（见图7）。

表  2  陨石坑熔化层分布

Table  2.  Distribution of melted layers in craters

工况	dm /m	mt /(1010 kg)
1		2.34
2		3.86
3		3.25
4	30～50	5.62
5		4.39
6	10～50	7.65
7	30～50	5.64
8	35～50	9.65

下载: 导出CSV 

| 显示表格

图  7  无量纲化的花岗岩完全熔化质量与撞击体直径的关系

Figure  7.  Dimensionless completely-melted granite mass varied with projectile diameter

下载: 全尺寸图片 幻灯片

依兰坑勘验数据为坑直径方向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个工况最终撞击坑直径分别为1770、1930、1840和2040 m，坑缘坑深分别为251、317、263和359 m。120 m、12 km/s的工况成坑的最终直径为1 840 m，坑缘坑深为263 m，由此可得该工况撞击条件下形成的陨石坑和依兰陨石坑坑形最接近，结合第3节中讨论的熔化层分布可知，最接近依兰陨石坑形成条件的为工况7。

图  8  不同模拟工况在计算结束时间150 s所对应的陨石坑剖面轮廓

Figure  8.  The crater profiles corresponding to the different simulation conditions at the end time of 150 s

下载: 全尺寸图片 幻灯片

表  3  8组模拟工况在150 s坑形数据

Table  3.  Crater shape data for eight simulation conditions at 150 s

工况	dp/m	u/(km∙s−1)	Df/m	dr/m	dr/dp
1	90	12	1470	182	2.02
2	90	15	1660	208	2.31
3	100	12	1610	227	2.27
4	100	15	1710	269	2.69
5	110	12	1770	251	2.28
6	110	15	1930	317	2.88
7	120	12	1840	263	2.19
8	120	15	2040	359	2.99

下载: 导出CSV 

| 显示表格

4.2   点源成坑相似律模型

点源相似律模型广泛用于将实验室尺度的冲击结果外延至行星级别的撞击结果^[21-23]。点源假设认为，当碰撞效应范围远大于小行星本身时，可以将小行星等效一个点源^[21]。小行星对撞击结果的影响仅由一个变量决定$C=a{u}^{\mu }{\delta }^{\nu }$，而非独立的3个$a$、$u$和$\delta$自变量^[21]。$a$、$u$和$\delta$分别是小行星半径、速度和密度，$\mu$是与材料相关的常数。$\mu$的理论值在1/3和2/3之间。对密实材料冲击实验结果表明$\mu$的值在0.55～0.6之间，非常高孔隙率材料的$\mu \approx 0.4$^[22]。$\nu$的值与材料特性无关，不同材料$\nu$的取值均为0.4^[23]。靶体的密度为$\rho$、抗拉强度$Y$、 黏度为$\eta$，重力加速度为$g$，撞击体与靶体剩下的材料参数与物质常数均由材料抗拉强度、黏度、质量密度3个基本物质常数无量纲化，记作${\varPi }_{\mathrm{m}}$。假设小行星垂直撞击地表，成坑半径$R$可由下式表示：

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

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

式中：$m$为小行星质量，$m=4{\text{π}} \mathrm{\delta }{a}^{3}/3$；H₁为常数，Y为靶体强度。

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

式中：H₂为常数，记${\varPi }_{{R}}=R{(\rho /m)}^{1/3}$，${\varPi }_{1}=Y/\left(\rho {u}^{2}\right)$。

在本文研究的依兰陨石坑撞击尺度条件下，成坑半径相似关系式中强度项大于重力项，即$Y/\left(\rho {u}^{2}\right) {\text{＞}} ga/{u}^{2}$，根据式(9)～(11)的推导可知，依兰陨石坑的形成属于强度机制：由上述点源成坑相似律分析可知，依兰陨石坑半径与撞击速度的$\mu$次方、靶体材料强度的$-\mu /2$次方成正比。因本文的撞击体与靶体密度相同，故无量纲化的密度项为1。根据本文8组工况的计算结果，以及额外追加计算的4组工况（直径110 m，撞击速度10、11、13和14 km/s） 的模拟结果来拟合式(11)，可得H₂=0.94377，$-\mu /2$=−0.22257，将工况7代入拟合曲线中，可得陨石坑最终坑直径为2030 m，与该工况的模拟结果相对误差为10.3%，工况7和拟合曲线对比如图9所示。

图  9  成坑半径拟合曲线与工况7数据的对比

Figure  9.  Comparison of the crater radius fitting curve and the data of condition 7

下载: 全尺寸图片 幻灯片

5.   结　论

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

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

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

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