引言

断层气是地球内部气体及挥发组分不断沿着活动板块和活动块体边界以及其它活动性断裂等地壳薄弱地带向地表迁移和释放的气体. 断裂带上的微量气体组分(如Rn，Hg，H₂，CO₂等)的变化，主要受断层和断层活动性的影响． 活动断层是地表与地下深部联系的通道，断层气从通道中溢出地表并携带大量的构造活动信息，是很好的构造指示剂(汪成民等，1991; Walia et al，2009； 吴华平等，2009; Voltattorni，Lombardi，2010)． 断层气在地震科学研究领域中应用并取得显著成果，大致归纳为以下几个方面： 探索地震，评价断层活动性，寻找隐伏断裂的位置，地震后效与震区烈度评定，以及灾害预测研究(陈万春，1996; Ioannides et al，2003)． 氡气作为断层气的一种成分，因其具有惰性、 迁移速度快、 穿透能力强、 对地震前兆反应明显及野外取样和测试简单易行，是目前探测隐伏断裂位置与评价断层活动性的一种有效手段．

氡的迁移机制与规律及浓度分布的研究始于20世纪早期，Flügge和Zimens(1939)最早提出氡运移是由扩散引起的. 它是建立在Fick定理基础上，用于解释氡气运移的一种理论(Султанходжаев и др，1979； Semkow，Parekh，1990; 贾文懿等，2000)． 随着后期不断的实践，研究人员又提出了对流作用，认为扩散与对流是氡迁移最重要的两个运移机制(Fleischer et al，1980; Malmquist et al，1989; Ioannides et al，2003; Swakon et al，2004; Iskandar et al，2005)． 因此对氡浓度的测量是一个复杂的过程，影响其含量的因素主要有： ① 岩土中铀、 镭等放射性元素的含量； ② 岩土的射气系数； ③ 岩土层的成分特征和覆盖层的厚度； ④ 岩土层的密度、 孔隙度等物理性质； ⑤ 气体的扩散系数(国家地震局科技监测司，1985)．

对于氡迁移的数值模拟研究领域，国外主要有： Abdoh和Pilkington(1989)基于对流和扩散机制理论，推导出氡在均匀覆盖层中迁移的三维浓度方程，并且可把三维的问题转化为二维来解决； Morin等(1993)推导出了三维扩散的解析解； Kohl等(1994)基于扩散-平流-衰减、 扩散-平流以及扩散-衰减等3种不同机制，模拟了氡从地表下到建筑物的迁移过程． 国内有刘菁华等(2007)根据Abdoh和Pilkington(1989)推导出的二维氡迁移浓度方程，运用有限差分法，计算对流-扩散机制下的氡浓度，模拟氡在断层上方均匀覆盖层中浓度的分布． Abdoh和Pilkington(1989)以及刘菁华等(2007)针对对流-扩散机制的各向同性和均匀覆盖层中氡的迁移问题，研究和探讨得都比较深入． 但从数值模拟角度，探讨氡在非均质覆盖层中的迁移规律、 浓度分布，并解释实测氡中遇到的典型现象与问题等方面，尚未深入研究．

关于断层上方覆盖层中氡浓度的分布规律，Swakon等(2004)、 张慧等(2005，2010)和张新基等(2005)研究发现的特殊现象主要有： ① 地表氡浓度异常与断层带位置不同步； ② 氡曲线出现单峰、 多峰及锯齿状等形态； ③ 断层上、 下盘的上方覆盖层氡浓度异常分布不对称； ④ 不同土壤属性的覆盖层出现不同的曲线形态． 针对上述现象本文分别提出了覆盖层内部含裂隙、 覆盖层下边界含裂隙及非均质覆盖层等3种物理模型，并从数值模拟角度切入，进行了讨论．

1.   氡迁移的基本物理模型和模拟原理

Abdoh和Pilkington(1989)提出的断层带上方覆盖层中氡的二维迁移模型，有3个假定条件： ① 氡迁移机制只有扩散和对流两种； ② 迁移的覆盖层为均质各向同性； ③ 土壤覆盖层中无放射性并处于稳定平衡状态． 本文所研究的非均质，其实质是在覆盖层中选取特殊区域，每个区域内部是均质和各向同性的，不同区域的同种参数设定不同的取值(如有效扩散系数、 氡气从土壤向空气的散失因子)，所以这种情况下的非均质也满足上述假设条件，Abdoh和Pilkington(1989)提出的迁移模型理论也适用于本文，基本物理模型如图 1所示．

图  1  含氡源断层带上方覆盖层的物理模型(引自Abdoh， Pilkington， 1989)

Figure  1.  The physical model of fault zone that contains a radon source and is covered by overburden (after Abdoh, Pilkington, 1989)

下载: 全尺寸图片 幻灯片

根据上述条件，推导出覆盖层中氡的浓度方程为(Abdoh，Pilkington，1989)

式中，D为扩散系数，e为覆盖层介质孔隙度，v为对流速率，λ为氡衰变系数(取值为2.07×10^-6 s^-1)，N=N(x，y)为氡浓度. 将有效扩散系数D^*=D/e代入式(1)中则为

相应的边界条件为(Abdoh，Pilkington，1989)

式中，a为土壤覆盖层宽度，b为土壤覆盖层厚度，h为氡气从土壤向空气的散失因子．

本文模拟氡在断层上方土壤覆盖层中的迁移规律，是由Matlab软件实现的． 利用偏微分求解工具pdetool绘制模型、 设定边界条件，并网格化输出所需的变量. 根据不同的物理模型改变相应的参数，最后采用具有Armijo-Goldstein线性搜索策略的阻尼高斯-牛顿迭代法(调用函数pdenonlin)求解氡迁移的二维浓度方程．

2.   数值模拟与讨论

2.1   模型1: 覆盖层内部含裂隙

2.1.1   模型建立与数值模拟

野外实际测量工作中，我们常会遇到氡气异常点与断层带位置点不相符的现象，如张慧等(2005)在金城关断裂布设的青白石测线与安宁碱水沟测线，测得氡异常与断裂带位置不同步． 为了解释这种现象，我们以裂隙倾角变化、 断裂带与裂隙偏移距离变化为变量，基于理想情况下的边值条件，建立了两个物理模型裂隙倾角变化模型(图 2a)和裂隙偏移断裂带模型(图 2b)． 假设在单个断层带上的覆盖层Ⅰ中含有裂隙Ⅱ，Ⅱ为断裂活动过程中形成的裂隙，裂隙Ⅱ在后期被土壤充填． Ⅰ与Ⅱ区域之间组构性质不同，但区域内部为各向同性均质的土壤，氡气向水平和垂直方向迁移，具体边界条件和物理模型如图 2所示．

图  2  断层带上方含裂隙覆盖层的物理模型 (a) 裂隙倾角变化模型; (b) 裂隙偏移断裂带模型

Figure  2.  (a) The model with variable dip-angle of the fracture; (b) The model with variable fracture deviation from fault zone

下载: 全尺寸图片 幻灯片

图 2a，b模型中W_f为断层带宽度，W₁为区域Ⅱ宽，t为倾角，h₂为散失因子，D^*₂为有效扩散系数. 另外，图 2b模型中W₂为区域Ⅱ高，W₃为地面上的区域Ⅱ中心到覆盖层中心的距离． 两模型的边界条件按式(3)来设定，N(x，0)=f(x)可写为

基于图 2a以裂隙的倾角t为变量，即t=30°，45°，65°和90°，区域Ⅱ连通覆盖层上下表面，且下端固定于断层带口，其它给定的边界条件与参数取值为a=20 m，b=4 m，W_f=2 m，N₁=200 Bq/m³，N₀=1500 Bq/m³. 区域Ⅰ中参数: D^*₁=0.02 cm²/s，v₁=0.0001 cm/s，h₁=0.02 cm^-1/s； 区域Ⅱ中参数: D^*₂=0.04 cm²/s，v₂=0.004 cm/s，h₂=0.2 cm^-1/s，W₁=2 m． 根据以上数据模拟结果如图 3所示．

图  3  裂隙倾角(t)变化对氡迁移的影响

Figure  3.  Effect of the changes in dip-angle(t) of the fracture on radon migration in overburden

下载: 全尺寸图片 幻灯片

基于图 2b以断裂带与裂隙的偏移距离W₃为变量，即W₃=3.5，5.5，7.5和9.5 m，其它给定的边界条件与参数取值为a=20 m，b=4 m，W_f=2 m，N₁=200 Bq/m³，N₀=1500 Bq/m³，t=45°. 区域Ⅰ中参数: D^*₁=0.01 cm²/s，v₁=0.0001 cm/s，h₁=0.02 cm^-1/s； 区域Ⅱ中参数: D^*₂=0.02 cm²/s，v₂=0.002 cm/s，h₂=0.2 cm^-1/s，W₁=1 m，W₂=3 m． 根据以上数据模拟结果如图 4所示．

图  4  裂隙偏移断裂带的距离(W₃)变化对氡迁移的影响

Figure  4.  Effect of the distance that the fracture deviates from the fault zone on radon migration in overburden

下载: 全尺寸图片 幻灯片

2.1   .2 模拟结果分析

1)裂隙的倾角变化对氡迁移的影响． 从图 2a的模拟结果(图 3)可以看出，在裂隙区域Ⅱ内氡的富集程度明显高于区域Ⅰ，随着裂隙倾角大小的变化，氡浓度富集的方向也在改变． 倾角依次从30°，45°，65°和90°变化，上方裂隙口位置氡浓度的异常峰值大小依次降低． 由以上分析得出，覆盖层中裂隙起到了加快氡迁移和促进氡富集的作用； 裂隙倾角的变化，引起氡迁移与富集的路径变化．

2)裂隙与断裂带的偏移距离对氡迁移的影响． 图 2b的模拟结果(图 4)表现出了4种偏移距离的氡迁移特征，它们之间氡浓度分布特征的区别很明显． 裂隙与断裂带偏移量越小，表现出的氡异常越明显； 随着偏移距离的增加，在覆盖层表面表现的异常峰值就越远离断层带地面投影位置，但是异常强度没有明显的改变． 由此可以得出，异常峰值位置变化与裂隙偏移断裂带的距离大小有关，异常强度大小则与其无关．

3)裂隙对氡异常峰值与断裂带位置同步性的影响． 由上述讨论及图 3、 图 4可以得出，氡异常峰值位置随着裂隙倾角、 断裂带偏移裂隙距离的改变而变化． 在该情况下就无法根据异常峰值确定隐伏断裂带的位置，二者在位置上表现出不同步现象． 模拟出现这种结果，是由于我们给定区域Ⅱ的有效扩散系数(D^*₂)、 对流速率(v₂)和散失因子(h₂)均高于区域Ⅰ所致． 由此可推断，土壤裂隙为氡气的运移提供了通道，氡气在此处富集运移，影响了地表氡气的测量，从而造成浓度曲线异常点与断层点不相符． 本文模拟结果与张慧等(2005)和张新基等(2005)野外测得的不同步现象一致，从数值模拟角度给出了地表氡异常位置与隐伏断裂位置不同步现象的解释．

2.2   模型2： 覆盖层下边界含裂隙系

2.2.1   模型建立与数值模拟

断裂是由于岩石受到的构造应力和能量超过了临界点促使岩石发生破裂并相对运动而形成的． 逆冲断层上盘活动性强，其节理和裂隙较下盘发育，导致上盘的氡异常强度高于下盘，且上下两盘氡浓度分布不对称(汪成民等，1991； 张新基等，2005)． 在实测中还会遇到氡异常峰值出现单峰、 多峰等不同情况，而且各峰值强度大小不同． 为了分析覆盖层对氡迁移规律的影响以及更好地解释实测结果，本文建立一个物理模型(图 5)，在均质各向同性的覆盖层下方有一个断层带，断层上盘有3条构造裂隙C₁，C₂和C₃． 图 5模型对应的边界条件为

图  5  裂隙系-断层带上方覆盖层的物理模型

Figure  5.  The model of fracture system-fault zone covered by overburden

下载: 全尺寸图片 幻灯片

f(x)的表达式如式(5)给定，是与断层带的位置相关的函数，在裂隙位置为二次函数，断层带处为定值N₀，其它范围为一次函数． 由于构造应力和断层活动性强度在断层带附近达到最大，所以假设距离断层带越近的裂隙，其宽度就越大，富集氡的能力就越强(富集过程不受其它自然因素影响)． 3条裂隙C₁，C₂和C₃的宽度分别为0.5，1.0，1.5 m，三者之间的等间距为1.5 m，C₃与断层带间距为1.5 m． 给定的参数值为a=30 m，b=6 m，W_f=4 m，N₀=3000 Bq/m³，v=0.001 cm/s，D^*=0.1 cm²/s，h=0.05 cm^-1/s． 假设在理想情况下，断层上方土壤覆盖层为均质各向同性，且对氡含量测定没有其它自然因素的影响. 该模型的数值模拟结果如图 6与图 7所示．

图  6  裂隙系-断层带上方覆盖层中氡迁移与断层位置关系图

(a) 氡迁移规律模拟三维立体图; (b) 氡迁移浓度分布剖面图

Figure  6.  Plot of radon migration with fault position in the overburden on the fracture system-fault zone

(a) 3D numerical simulation of radon migration； (b) Distribution of radon concentration

下载: 全尺寸图片 幻灯片

图  7  断层-裂隙系上方的覆盖层厚度(b)与地表氡浓度峰值关系图(氡浓度单位： Bq/m³)

Figure  7.  Relationship of thickness (b) of the overburden and radon concentration peak on surface above the fracture system-fault zone (radon concentration units: Bq/m³)

下载: 全尺寸图片 幻灯片

2.2   .2 模拟结果分析

图 6a，b显示，在覆盖层下边界，氡浓度含量与断层带位置有关． 裂隙C₁，C₂和C₃氡浓度在图中表现为： 图 6a中3个峰值高度依次增高，图 6b中三者的氡浓度等值线晕圈依次增大． 其它非裂隙范围的上方覆盖层中氡浓度分布以断层带边界的中线为对称轴向两边递减，与断层带的距离越远，浓度值越低．

1)覆盖层厚度对氡浓度曲线形态和迁移规律的影响． 为了研究覆盖层厚度对氡迁移的影响，在图 6基础上，对模拟覆盖层厚度分别为3.0，1.0，0.65 m的氡迁移作比较分析，结果如图 7所示． 在同样的裂隙系-断层带边值条件下，根据图 7中异常晕圈的晕色深浅和个数可以看出，厚度为3.0，1.0，0.65 m的覆盖层中氡浓度在地表的异常强度依次增大，异常峰值个数分别是1，2，3，即地表表现为单峰、 双峰以及锯齿状的氡浓度曲线形态．

结合图 6b和图 7分析得出，覆盖层厚度越小，地表表现的氡异常强度就越大； 厚度增大时，氡浓度曲线异常低缓，等值曲线宽度加大． 由图 6b可看出，覆盖层内部氡异常明显，而在其顶部氡曲线没有异常峰值． 而实测中也有类似情况，如张慧等(2005)在金城关断裂布设的安宁关山沟测线与西固蛤蟆滩测线，断层带上方覆盖层表面没有氡异常，不能根据其测量数据确定断层带位置． 张慧等(2010)通过野外场地试验得出的结论表明，正常情况下，覆盖层厚度不大时峰值突出，异常明显； 当厚度增大时异常形态低缓、 宽度加大． 因此模拟结果符合实测结果．

2)覆盖层下边界的裂隙对氡异常峰值对称性的影响． 据图 6和图 7可知，由于断层上盘含有3条裂隙导致断层带上方覆盖层上、 下盘氡异常峰值分布不对称，上盘的异常强度以及异常带宽度均大于下盘． 模拟结果符合前人实测结论： 土壤气正异常区位于断裂上盘的地表投影，而断层上断点的地表投影一般为土壤气浓度值由低到高的转折陡变部分(汪成民等，1991； 张新基等，2005； 张慧等，2010)．

2.3   模型3： 非均质覆盖层

2.3.1   模型建立与数值模拟

上述两个模型分别讨论了内部含裂隙、 下边界含裂隙系的覆盖层中氡的迁移规律. 为了使模型和计算更切合实测，本文在Abdoh和Pilkington(1989)及刘菁华等(2007)的基础上提出非均质覆盖层研究模型． 假设一个非均质模型(图 8)，覆盖层由致密性和疏松性两种土壤组成，致密性土壤夹有一层疏松性土壤，前者的土壤类型是含少量砂质的黏土，后者的是砂砾质疏松沉积物，致密层和疏松层内部结构为均质各向同性． 模型的边界条件如下:

图  8  断层上方非均质覆盖层物理模型

Figure  8.  The model of the fault zone covered by nonhomogeneous overburden

下载: 全尺寸图片 幻灯片

土壤的疏松层与致密层对氡迁移机制的影响程度不同，主要来自于有效扩散系数和对流速率两个因素． 对于不同属性的土壤，其有效扩散系数、 对流速率的取值至今没有一个定性的理论标准． 表 1是根据吴慧山等(1995)给出的资料，得出部分土壤的有效扩散系数(D^*)、 孔隙度(e)及扩散系数(D)之间的关系. 孔隙度越大扩散系数就越大，但由D^*=D/e可知，有效扩散系数的值取决于二者的大小． 在模拟中，假设疏松层有效扩散系数和对流速率都比致密层大，给定的参数值有： a=20 m，b= 6 m，W_f=4 m，N₀=1500 Bq/m³，N₁=200 Bq/m³，h=0.05 cm^-1/s. 致密层参数: v₁=0.0002 cm/s，D^*₁=0.02 cm²/s; 疏松层参数: v₂=0.0007 cm/s，D^*₂=0.2 cm²/s． 而各向同性均质的覆盖层的参数有： D^*₁=D^*₂=0.02 cm²/s，v₁=v₂=0.0002 cm/s，运用氡气二维迁移方程模拟结果如图 9与图 10所示．

表  1  氡气在不同土壤性质下的有效扩散系数(D^*)、 扩散系数(D)和孔隙度(e)三者之间的关系(据吴慧山等， 1995)

Table  1.  The relationship among effective diffusion coefficient (D^*), diffusion coefficient (D) and porosity (e) of radon in different types of soils (Wu et al, 1995)

土壤性质	孔隙度e	扩散系数D/(10-2 cm2·s-1)	有效扩散系数D*/(cm2·s-1)
砂子	40%	4.5—7.0	0.11—0.175
疏松沉积物	20%	2.0—2.5	0.1—0.125
白黏土	59.3%	1.53	0.023
砂质黏土	10.8%	1.09	0.1

下载: 导出CSV 

| 显示表格

图  9  非均质覆盖层(a)与均质覆盖层(b)氡迁移规律模拟三维立体图

Figure  9.  3D plot of numerical simulation on radon migration in nonhomogeneous (a) and homogeneous overburdens (b)

下载: 全尺寸图片 幻灯片

图  10  非均质覆盖层与均质覆盖层氡迁移浓度分布等值线对比图

Figure  10.  The contour plots of radon concentration distribution in nonhomogeneous and homogeneous overburdens

下载: 全尺寸图片 幻灯片

2.3   .2 模拟结果分析

1)非均质覆盖层中氡的浓度分布规律． 图 9和图 10表明，疏松层与致密层之间氡的迁移规律明显不同，非均质与均质之间不同． 受到对流速率和有效扩散系数的影响，在疏松层与致密层交界处氡气迁移发生明显的转变. 位于断层上方覆盖层2—4 m的疏松层，其氡气浓度值在水平和垂直两个方向变化梯度较致密层低一些，迁移强度大于致密层．

图 10显示，非均质与均质覆盖层的迁移规律有很大的区别． 在覆盖层厚度3 m处左右，氡浓度等值的蓝线与绿线相交，交点所在深度以上部分表现为非均质覆盖层中的氡迁移强度在水平和垂直方向上都大于均质覆盖层，相同浓度等值线非均质的宽度和高度都大于均质，并且每个深度的横剖面曲线形态均类似于地表下2 m处的氡浓度曲线(图 11)． 而交点以下部分则恰好相反． Swakon等(2004)通过氡测量，验证了不同属性土壤中氡的浓度分布也不同，疏松层覆盖层的氡扩散和对流作用更强． 本文的模拟结果与该理论相符合．

图  11  断层带上方地表下2 m 处的氡浓度曲线

Figure  11.  Radon concentration curve 2 m beneath the ground across the fault zone

下载: 全尺寸图片 幻灯片

2)非均质覆盖层对氡浓度曲线形态的影响． 假设在非均质和均质覆盖层同一深度分别布设一条测线(图 10)，其测量曲线如图 11所示． 从图 11可以看出，本模型中非均质与均质覆盖层在该深度处的氡浓度曲线均为单峰，异常点与断层带位置相符合，非均质覆盖层的氡浓度曲线明显加宽，即非均质覆盖层中的迁移强度大于均质． 上述提及的图 10中蓝、 绿色等值线相交点以上部分，非均质覆盖层中的氡浓度曲线与均质相比都有所加宽，加宽量大小随深度变化而不等． 根据假设条件，不难得出加宽的缘由是非均质覆盖层中疏松层的扩散-对流强度较大，对氡迁移起到了促进作用. 这一结果验证了Swakon等(2004)的理论．

3.   结论

根据氡迁移二维浓度方程，基于Matlab平台模拟了不同条件下覆盖层中氡的迁移规律，从模拟结果分析得出如下结论:

1)通过对内部含裂隙的覆盖层中氡迁移的数值模拟，得出土壤裂隙的存在为氡气的运移提供了通道，造成地表氡异常位置与断层带位置不同步. 该模拟结果与张慧等(2005)和张新基等(2005)的实测现象相符合．

2)通过对下边界含裂隙系的覆盖层中氡迁移的数值模拟，得出断层带上盘裂隙发育，氡异常强度高于下盘； 断层上下盘氡浓度分布不对称； 地表测得氡浓度曲线随着覆盖层厚度的变化会出现单峰、 双峰以及锯齿等形状．

3)通过对非均质覆盖层中氡迁移的数值模拟，得出覆盖层土壤结构属性对氡迁移有明显影响． 由于对流-扩散机制的影响，疏松层与致密层相比，氡浓度梯度变化小； 水平与垂直两个方向迁移强度高． 在理想模式下，通过与均质覆盖层对比，分析了非均质覆盖层对氡气曲线形态的影响．

4)通过对3种模拟结果总结得出，土壤裂隙为氡气的运移提供了通道，促使氡气富集； 覆盖层对氡异常有削弱作用，尽管在覆盖层下部分氡异常强度比较大，但当厚度超过一个上限时，在断层带地表投影处的氡含量不会出现异常值．

5)根据本文提出的3种模型，对断层气氡在覆盖层中迁移规律以及实测中的典型现象与问题进行了讨论． 数值模拟结果与实测有很好的一致性，表明本文的数值模拟方法，在利用氡气的运移判断隐伏断裂位置等方面的运用，具有一定的积极意义．