• ISSN 1673-5722
  • CN 11-5429/P

郯庐断裂带宿迁段合欢路土壤氡分布特征与迁移特征的数值模拟

戴波 赵启光 张敏 张扬 周晓成

娄世平, 杨玉永, 刘瑞峰, 徐秀杰, 贾荣光, 董翔. 基于Android智能终端的地震现场应急指挥技术系统运维信息管理平台[J]. 震灾防御技术, 2018, 13(3): 727-735. doi: 10.11899/zzfy20180324
引用本文: 戴波,赵启光,张敏,张扬,周晓成,2021. 郯庐断裂带宿迁段合欢路土壤氡分布特征与迁移特征的数值模拟. 震灾防御技术,16(1):220−228. doi:10.11899/zzfy20210123. doi: 10.11899/zzfy20210123
Lou Shiping, Yang Yuyong, Liu Ruifeng, Xu Xiujie, Jia Rongguang, Dong Xiang. The Information Operation and Maintenance Management Platform of Earthquake Field Emergency Command Technology System Based on Android Intelligent Terminal[J]. Technology for Earthquake Disaster Prevention, 2018, 13(3): 727-735. doi: 10.11899/zzfy20180324
Citation: Dai Bo, Zhao Qiguang, Zhang Min, Zhang Yang, Zhou Xiaocheng. Numerical Simulation of Radon Concentration in Soil of Hehuan Road in Tan Lu Fault Suqian Segment[J]. Technology for Earthquake Disaster Prevention, 2021, 16(1): 220-228. doi: 10.11899/zzfy20210123

郯庐断裂带宿迁段合欢路土壤氡分布特征与迁移特征的数值模拟

doi: 10.11899/zzfy20210123
基金项目: 宿迁市活动断层探测与地震危险性项目(SQZBTB20110162);江苏省产业前瞻与关键核心技术-竞争项目(BE2020116);2021年度第一批江苏省地方标准项目(124)
详细信息
    作者简介:

    戴波,男,生于1987年。工程师。主要从事断裂带地球化学与地震活动性方面的研究。E-mail:okman310@foxmail.com

Numerical Simulation of Radon Concentration in Soil of Hehuan Road in Tan Lu Fault Suqian Segment

  • 摘要: 沿郯庐断裂带宿迁段合欢路布设土壤氡测线并对测得的数据进行分析研究,结果表明研究区域内土壤氡探测结果对断裂带位置、断层类型和特征有较好的指示性。在已知场地地质构造参数、地层岩性特征的基础上,建立土壤氡迁移模型。建立模型时识别了断裂类型和覆盖层分层特征,并讨论了断层内部结构与氡浓度取值,使模型与实际情况相等。数值模拟结果印证了土壤氡异常峰区间在圈定断裂位置和断层破碎带中的重要作用,揭示了氡在不同覆盖层中的迁移速度,并定量解释了合欢路场地断层土壤氡异常的成因。
  • 起伏界面波场复杂,在不规则界面上既产生波的反射、散射和绕射,又产生面波,更有受界面起伏特征影响的散射波等(孙建国,2007)。一般情况下,为处理复杂边界,需增加很多甚至是数倍于规则边界的计算量。为精确模拟起伏界面条件下的地震波传播特征,需解决两个关键问题——计算方法选择与模型网格精确剖分。近年来,提出的数值模拟计算方法包括有限差分法(FEM)(Levander,1988Olsen等,19951997Graves,1996Olsen,2000Moczo等,2001)、伪谱法(PSM)(Patera,1984Furumura等,1998)、有限元法(FDM)(Smith,1975Koketsu等,2004;)、谱元法(Patera,1984Seriani,1998Komatitsch等,1998199920002002)。网格剖分方式可分为网格法和非网格法两大类(孙建国,2007)。网格法分为规则网格(Moczo,1989Jastram等,1994)、非结构网格(Tessmer等,1992Zahradník,1993)和贴体网格(垂向变换)(Tessmer等,1992Hestholm等,1994)等。近年来,国内外对于网格法开展了系列研究工作。规则网格差分法最早由Alterman等(1968)引入地震波数值模拟计算中,并利用规则网格差分法完成了层状介质中二阶弹性波波动方程的模拟,随后Alford等(1974)分析了网格选取与模拟精度之间的关系。Madariaga(1976)提出了交错网格差分格法,并很好地应用于正演模拟中,但该方法使用低阶差分,精度较低,随后Dablain(1986)提出了任意阶差分精度算法,Crase等(1990)实现了高阶差分算法,因此有限差分法数值模拟精度被大大提高。但采用规则网格剖分时,如果界面起伏较大,会导致模型界面阶梯状离散,进而产生数值散射(Xue,2014)。六面体(四边形)单元的谱元法已成功应用于模拟波的传播(Priolo等,1994Komatitsch等,20022005),在处理不规则界面时,显示均匀介质内差异过大的网格会激发强的网格反射波(Zhou等,2010),由于其均基于规则网格实现,无法从根本上消除阶梯状网格引起的虚假绕射波,为此学者们发展了不规则网格有限差分法。

    不规则网格有限差分法在有限元法中的非结构网格适用性较好,易生成网格,但对精度有一定影响,三角形单元增加了存储量,降低了效率和准确性(Komatitsch等,2000)。贴体网格(垂向变换)思路相近,依照实际地表起伏剖分曲线网格,使用较方便,适用于同位网格和交错网格(Hestholm等,2002Zhang等,2006),但对网格变形范围具有一定限制(Tessmer等,1992)。Pitarka(1999)提出了各向同性介质中矩形不规则网格有限差分法。Hestholm等(2002)和Tarrass等(2011)提出了以贴体网格的方式对起伏地表进行剖分。Zhang等(2006)和Lan等(2011)利用同位网格的差分方式对曲坐标系下的波动方程进行差分求解,实现了正演模拟,且取得了较高精度。Zhang等(2019)以不同网格分辨率的模拟,强调了复杂情况下大规模、高分辨率地震波模拟的重要性。

    本文通过谱元法分析六面体单元剖分起伏界面的不同方式对地震波场的影响。对于地下介质中传播的地震波,可由波动方程、初始条件和边界条件进行描述,即在$ \mathrm{\Omega } $中地震波位移满足以下方程:

    $$ \begin{array}{c}\rho \ddot{u}=\nabla \cdot{\boldsymbol{\sigma}} +f\\{\boldsymbol{ \sigma}} =\mathit{C}:\varepsilon \\ \varepsilon =\dfrac{1}{2}\left[\nabla{\boldsymbol{ u}}+(\nabla {\boldsymbol{u}}{)}^{{\rm{T}}}\right]\end{array} $$ (1)

    式中,$ {\boldsymbol{u}} $表示位移向量,$ {\boldsymbol{\sigma}} $表示应力张量,$ \boldsymbol{C} $表示刚度张量,$ \rho $表示密度,$ f $表示外力,“:”表示张量相乘。

    在谱元法中,方程(1)点乘任意试验权函数$ w $,并进行分步积分,可得到:

    $$ \begin{array}{c}{\displaystyle\int }_{\mathrm{\Omega }}\rho w\cdot \ddot{{\boldsymbol{u}}}{\rm{d}}\Omega +{\displaystyle\int }_{\mathrm{\Omega }}\nabla w:{\boldsymbol{C}}:\nabla {\boldsymbol{u}}{\rm{d}}\Omega ={\displaystyle\int }_{\mathrm{\Omega }}w\cdot f{\rm{d}}\Omega +{\displaystyle\int }_{{\Gamma }_{\text{abs}}}w\cdot t{\rm{d}}\Gamma \end{array} $$ (2)

    式中,$ \mathrm{\Omega } $表示物理计算域,$ \Gamma $为其边界,包括自由地表边界$ {\Gamma }_{{\rm{f}}} $和人工吸收边界$ {\Gamma }_{{\rm{abs}}} $

    由于地表应力为零,$ {\Gamma }_{\mathrm{f}} $在部分过程中消失,因此式(2)自然满足地表条件。

    在参考单元Λ中,引入N阶多项式的基函数,其在单元$ {\mathrm{\Omega }}_{\mathrm{e}} $上对应1组节点,$ {u}_{N}^{{\rm{e}}} $$ {w}_{N}^{{\rm{e}}} $分别为$ u $$ w $在这些节点上的Legrange插值。这些节点在[−1,1]区间内取为Gauss-Lobatto-Legendre点,是下式的根:

    $$ \begin{array}{c}\left(1-{\xi }^{2}\right){P}_{N}^{\mathrm{{'}}}\left(\xi \right)=0\end{array} $$ (3)

    式中,${P}_{N}^{\mathrm{{'}}}\left(\xi \right)$为Legendre N阶多项式的偏导数。

    在参考单元Λ中,$ {u}_{N}^{e} $表示为:

    $$ \begin{array}{c}{u}_{N}^{e}\left(\xi ,\eta ,\gamma \right)=\displaystyle\sum\limits_{p=0}^{N}{u}_{N}^{e}\left({\xi }_{p},{\eta }_{p},{\gamma }_{p}\right)h\left({\xi }_{p}\right)h\left({\eta }_{p}\right)h\left({\gamma }_{p}\right)\end{array} $$ (4)

    式中,$ h\left({\xi }_{p}\right) $表示$ p $项一维Legrange插值。

    引入多项式近似,方程(2)的积分可利用Gauss-Lobatto-Legendre积分规则得到:

    $$ \begin{array}{c}\displaystyle\int {u}_{N}{w}_{N}d\Omega =\displaystyle\sum\limits_{e=1}^{{n}_{el}}{\displaystyle\int }_{{\Omega }_{e}}{u}_{N}^{e}{w}_{N}^{e}d\Omega \approx \displaystyle\sum\limits_{e=1}^{{n}_{el}}\displaystyle\sum\limits_{i=0}^{N}{w}_{i}\displaystyle\sum\limits_{j=0}^{N}{w}_{j}\displaystyle\sum\limits_{k=0}^{N}{w}_{k}{\mathit{J}}_{\mathit{e}}\left({\xi }_{i},{\eta }_{j},{\gamma }_{k}\right)\cdot {u}_{N}^{e}\left({\xi }_{i},{\eta }_{j},{\gamma }_{k}\right)w\left({\xi }_{i},{\eta }_{j},{\gamma }_{k}\right)\end{array} $$ (5)

    式中,$ {w}_{i} $为积分系数,$ {\boldsymbol{J}}_{i} $为Jacobian矩阵,由影射$ {\mathrm{\Omega }}_{\mathrm{e}} $到Λ的关系$ {\Gamma }_{e} $决定。空间离散化后,进行时间离散化,用表示介质位移未知量,得到:

    $$ \begin{array}{c}M\ddot{U}+KU=F\end{array} $$ (6)

    式中,$ \mathit{M} $为质量矩阵,$ \mathit{K} $为刚度矩阵,为源项。

    Legendre谱元的不重叠单元的空间离散,以Gauss-Lobatto-Legendre节点插值多项式,结合Gauss-Lobatto-Legendre积分规则,使质量矩阵$ \mathit{M} $为对角线矩阵,$ \mathit{K} $为带宽矩阵。$ \mathit{M} $的对角性使谱元法计算效率较有限元法有了较大提高。

    常规的网格剖分方式中多采用规则的阶梯状网格离散地质模型(Dablain,1986Levander,1988董良国等,20002004),如图1所示。在进行数值模拟计算时,当遇到有起伏界面的地质模型时,直接使用矩形网格离散界面会导致阶梯状的空间离散。如果网格剖分过于稀疏,如图2(a)所示,单元尺寸差异很大,在计算中会导致发散;如果网格剖分过于密集,如图2(b)所示,在加密网格后采用台阶状网格近似描绘起伏地表(Robertsson,1996)时,会大大增加计算量,占用计算机内存,且无法彻底消除阶梯的散射问题。为解决上述问题,须寻找使网格尽可能保持均匀、剖分尽可能贴合起伏界面的剖分方式。

    图 1  矩形网格剖分起伏界面
    Figure 1.  Rectangular mesh generation of undulating interface
    图 2  矩形网格剖分示意图
    Figure 2.  Schematic diagram of rectangular mesh generation

    (1)问题一:深度方向(z轴方向)不同单元个数之间的过渡

    取起伏界面处的1个台阶状网格进行研究(图3),为使网格贴合起伏界面,并保证均为正交或略倾斜的六面体单元,需考虑深度方向(z轴方向)不同单元个数之间的过渡。由图3可知,二维矩形网格剖分产生了阶梯状界面。

    图 3  台阶状网格截面示意图
    Figure 3.  Schematic diagram of the stepped mesh section

    三维模型在地形等高线图上所处位置如图4(a)所示,其中x轴为南北向,y轴为东西向,z轴为深度方向,模型东西长6 000 m,南北长5 000 m,深度为1 100 m。起伏界面位于深度600 m处,最深处为700 m,起伏界面将模型分为上下两部分,图4(a)中红色外框对应上半部分模型深度为600 m的范围,红色内框对应上半部分模型深度为700 m的范围。图4(b)中z轴方向上红色外框范围的单元较红色内框范围单元少1个,产生了阶梯状网格。考虑到三维模型在2个方向(xz平面与yz平面)上均有阶梯状网格(图5),因此不同网格个数之间的过渡在2个方向上均需进行。当2个方向上单元过渡后,将拐角处(图5中点A、B、C、D)的单元进行了2次叠加二合一处理,导致单元发生扭曲畸变。

    图 4  三维模型在地形中所处位置示意图
    Figure 4.  Location of the 3D calculation model in terrain
    图 5  起伏界面上与下的阶梯状网格完成单元二合一处理后变为斜面网格示意图
    Figure 5.  The stepped mesh above and below the undulating interface become inclined grids using two-in-one unit method

    (2)问题二:扭曲单元的衔接

    当2个方向(xz平面与yz平面)均完成问题一中深度方向上的单元个数过渡后,需考虑2个方向拐角处因进行了2次过渡叠加而发生扭曲的单元衔接问题。

    本文引入TrueGrid网格前处理软件,针对问题一,对含起伏界面的研究区域进行六面体单元网格剖分,在固有台阶状网格处使用六面体单元二合一的方式,如图6(a)、6(b)所示,将1个台阶状单元网格通过扯动单元的方式,通过1个变形的六面体单元进行过渡衔接,将深度方向(z轴方向)上的不同单元个数进行过渡衔接。在上下台阶网格处使用该方法,将图3所示台阶状网格变为图6(c)所示的连续网格。

    图 6  台阶状网格处使用单元二合一的方式处理过程示意图
    Figure 6.  process of using two-in-one unit method on stepped mesh

    完成单元二合一处理后的网格模型如图7所示,由图7可知,模型消除了台阶状网格,将台阶状网格变为了斜面网格,可将斜面投影至起伏界面,使网格完全贴合界面。但在4个拐角处,由于叠加了2次二合一处理,使拐角处单元发生扭曲畸变,因此,需单独处理4个拐角处的单元,即解决问题二。

    图 7  起伏界面上与下的阶梯状网格完成单元二合一处理后变为斜面网格的模型示意图
    Figure 7.  The stepped mesh above and below the undulating interface become inclined grids using two-in-one unit method

    本文采用3种剖分方式,以对比不同方式的影响。由于4个拐角处理方式一致,本文仅以1个拐角为例。

    (1)剖分方式一:在拐角处删除四纵列单元

    在网格剖分时单元个数变化策略中,将16个单元合并为12个是网格变化策略中常用的方式(蒙茂洲,2010),如图8所示,该方法可使网格处理较规整。拐角处单元发生的扭曲畸变如图9(a)所示,删除四纵列扭曲单元如图9(b)所示,采用将16个单元合并为12个单元的方法调整其余单元位置,补充至删除的四纵列单元的位置上。为满足计算要求,略调整了12个单元的中心位置,如图9(c)所示,至此完成了剖分。

    图 8  16个单元合并为12个单元的过程示意图
    Figure 8.  The process of merging 16 units into 12 units
    图 9  剖分方式一的过程示意图
    Figure 9.  The process of mesh generation of the first method

    (2)剖分方式二:在拐角处删除一纵列单元

    将16个单元合并为12个单元的网格数量调整方法删除网格数量较多,且调整后倾斜的单元数量较多,采用单元个数变化策略中的将4个单元合并为3个单元的方式,如图10所示,局部调整方式可减少倾斜的单元个数(罗宏伟等,2013)。拐角处的单元发生扭曲畸变,如图11(a)所示,删除一纵列扭曲单元后如图11(b)所示。采用将4个单元合并为3个单元的方法,调整其余单元的位置,补充至删除的一纵列单元的位置上,如图11(c)所示。

    图 10  4个单元合并为3个单元的过程示意图
    Figure 10.  The process of merging 4 units into 3 units.
    图 11  剖分方式二的过程示意图
    Figure 11.  The process of mesh generation of the second method

    (3)剖分方式三:将图7中的斜面部分单独剖分为1个过渡环

    过渡环内、外块体相差1个单元数,即1个台阶。过渡环内、外侧网格充分贴合完成剖分,如图12所示,可知过渡环的设置使拐角处单元变形较小。

    图 12  剖分方式三的过程示意图
    Figure 12.  The process of mesh generation of the third method

    由TrueGrid软件可进行3种剖分方式生成的整体网格质量评价,主要指标包括单元扭曲度及正交性等,评价结果如图13所示。由图13可知,网格质量良好,基本为正交的六面体单元。扭曲度是检查相对的2个角之间的差值,该值越小,单元扭曲度越小,由图13(a)可知,大部分单元相对2个角之间的差值为0。正交性是检查每个单元节点处的角度与90°的差值,该值越小,单元正交性越好,由图13(b)可知,大部分单元节点处的角度与90°的差值为0。

    图 13  整体网格质量评价结果
    Figure 13.  Evaluation results of overall mesh quality

    本文选择的研究区域(图14红色框)位于京津冀地区的第四系地层局部起伏界面,附近地震活动频繁,是强震灾害多发区。根据物探及钻孔资料(金安蜀等,1980国家地震局,1986孙若昧等,1996赵金仁等,1999王夫运等,2005),在既有北京盆地模型(高孟潭等,2002)的基础上叠加了地形条件,建立了北京盆地含起伏地形的地下介质三维速度结构模型(付长华等,2012)。x轴为南、北方向,y轴为东、西方向,z轴为深度方向,整体深度为1100 m,盆地最深处为700 m,深度分布如图15所示,剖分网格的单元长度为100 m,表1所示为介质和震源模型参数,震源时间函数通过主频为4 Hz的Ricker子波表示,模拟的整个地震波传播时长为4 s。在三维地震波场模拟中计算量仍较大,且对储存能力的要求较高,为此,使用天河二号超级计算机进行并行计算,计算区域被剖分为9个进程。整个区域除自由地表外,其余表面均设置为吸收边界,以消除地震波在区域边界的反射。

    图 14  研究区域第四系沉积地层地形
    Figure 14.  Topographic map of quaternary sedimentary strata in the study area
    图 15  研究区域深度分布
    Figure 15.  Depth distribution map of the study area
    表 1  模型基本参数
    Table 1.  Basic parameters of the model
    模型尺寸/m 总步长/步 格点间距/m 时间步长/s 震源时间函数主频/Hz 震源时间函数
    6 000×5 000×1 100(长×宽×深) 4 000 100 0.001 4 主频为4 Hz的Ricker子波
    下载: 导出CSV 
    | 显示表格

    观测点布置位置如下:①测点位置一如图16中红色点所示,3种剖分方式最大的区别在于拐角处理,在3种方式网格模型中六边形单元畸变最大的位置,即在拐角处布置4组测点如图17所示,在3种模型中取同一位置,以反映在拐角连接处不同剖分方式在拐角处的单元模拟结果的差异;②测点位置二如图16中蓝色点所示,在建立的三维网格模型中水平地表布置测点,在3种模型中取同一位置,以反映在水平地表处不同剖分方式对距拐角不同距离处的单元模拟结果的差异。

    图 16  模型中测点位置的俯视图
    Figure 16.  Top view of test points position in the model
    图 17  模型中测点位置一的截面图
    Figure 17.  Section view of test point position one in the model

    3种剖分方式在测点位置一处的模拟结果对比如图1819所示。根据Ricker子波理论波形可知,起伏界面位于模型600 m深度附近,当波传至起伏界面时,剖分方式二、三波形较规则,剖分方式一在0.2、2 s左右波形出现了不规则散射,剖分方式二、三时程曲线较吻合。当模型取均匀介质,即介质波速比为1∶1时,波的持时为2.7 s左右;当模型取非均匀介质,即介质波速比为0.8∶1时,波的持时为3.3 s左右,这符合软弱介质使持时增长的原理。

    图 18  均匀介质中3种剖分方式在测点位置一处四纵列单元的模拟时程曲线结果对比
    Figure 18.  Comparison of the simulation time history curves of three mesh generation methods in a uniform medium with four rows of elements at the test point one
    图 19  非均匀介质中3种剖分方式在测点位置一处四纵列单元的模拟时程曲线结果对比
    Figure 19.  Comparison of the simulation time history curves of three mesh generation methods in a four-row unit at a test point in heterogeneous medium

    图18可知,在深度700 m台阶网格上下附近,波形出现不规则的散射。在网格中取相同位置的单元直接使用阶梯状网格模拟的时程曲线如图20所示,对比图1819可知,剖分方式二、三改善了台阶处网格在数值模拟时产生的散射现象。

    图 20  非均匀介质中直接使用台阶状网格在测点位置一处四纵列单元的模拟时程曲线结果对比
    Figure 20.  Comparison of simulation time history curve results of four rows of units at a test point directly using stepmesh in non-uniform medium

    测点位置二如图16中蓝色点所示,分布于模型水平面的8个点位,每个点位间隔100 m,且从1点至8点离拐角单元距离越来越近。3种剖分方式在测点位置二处的模拟结果对比如图21所示,由图21可知,对于水平地表单元而言,当距拐角900 m左右时,3种剖分方式时程曲线基本相同,且距拐角距离越远,时程曲线越吻合;剖分方式二、三时程曲线较一致,而当距拐角距离<300 m时,剖分方式一存在较大差异。

    图 21  3种剖分方式在测点位置二的时程曲线对比
    Figure 21.  Comparison of the time history curves of the three mesh generation methods at the second test point

    采用谱元法进行数值模拟的优势在于计算精度高、可较好地适应起伏界面、可克服有限元计算量大的缺点,因此被广泛使用。在以往的研究中,通常缺少灵活的六面体单元贴合起伏界面的剖分方式。利用谱元法规则六面体单元进行网格剖分时,界面起伏较大处会导致阶梯状离散,进而产生数值散射。为消除阶梯网格对起伏界面地震动模拟计算结果的影响,本文采用的剖分方式使网格尽可能保持均匀,且剖分尽可能贴合起伏界面。本文取得的研究进展主要有:

    (1)本文基于TureGrid软件提出了模拟计算使用的起伏界面处地形六面体单元网格剖分方式,该剖分方式在固有台阶状网格处使用六面体单元二合一的方式将深度方向上的不同单元个数进行过渡衔接。采用Fortran语言编写了应用程序,通过该应用程序可快速建立起伏界面处的地形六面体网格模型。考虑到三维模型在2个方向(xz平面与yz平面)上均有台阶状网格,本文采取了删除四纵列拐角单元、删除一纵列拐角单元及构造单元过渡环3种剖分方式,解决2个方向上单元二合一后拐角处产生扭曲单元的问题,做到了六面体单元尽可能正交,网格尽可能均匀,剖分尽可能贴合起伏界面。

    (2)本文通过3种剖分方式建立整体六面体网格模型,通过TrueGrid软件进行整体网格质量评价,单元扭曲度及正交性评价结果显示网格质量良好,基本为正交的规则六面体单元,符合网格剖分要求,满足谱元法计算需求。

    (3)3种模型通过谱元法进行数值模拟计算,验证了,剖分方式的正确性。对于拐角处的单元而言,剖分方式一会引起波形散射,而剖分方式二、三均为规则波形。对于远离拐角处的单元而言,当距拐角距离>900 m时,3种剖分方式模拟出的结果较一致;当距拐角距离为300~900 m时,剖分方式一与剖分方式二、三存在较小差异;当距拐角距离<300 m时,剖分方式一与剖分方式二、三存在较大差异,剖分方式二、三差异始终较小。这表明删除拐角处一纵列单元方式与设置过渡环的方式均可使用,而删除四纵列拐角单元方式不推荐使用。

    在起伏界面处,本文涉及的网格剖分后的剖分处理方式优势包括网格优势和程序优势。网格优势体现在:①剖分使用六面体网格,六面体网格在计算精度及存储空间等方面具有显著优势(王瑞等,2020);②改善了规则网格剖分时利用台阶状网格近似描绘起伏界面的问题,消除了数值模拟时阶梯状边界引起的绕射和散射等问题;③剖分的网格尺寸尽可能保持均匀,避免因网格尺寸差异导致计算过程中出现发散;④剖分贴合起伏界面,提高了数值模拟的精度。程序优势体现在:基于TrueGrid软件编写了应用程序,提出了模拟计算用的起伏界面处地形六面体单元网格剖分方式,通过该程序可快速建立起伏界面处的六面体网格模型。

    在今后的工作中,可将剖分方式二或三改为循环程序,通过输入起伏界面不同深度的边界坐标可得到该界面内三维精细化六面体网格剖分结果,可将程序运用到大尺度、大范围地形效应的模拟中。通过精准的网格剖分投影贴合起伏界面,以达到较高的计算精度,从而解决用六面体网格离散起伏界面模型中产生的阶梯状界面,提高谱元法处理起伏界面问题的灵活性。

    致谢 国家超级计算天津中心提供了超算计算平台,在此表示感谢。

  • 图  1  活断层上氡迁移的物理模型

    Figure  1.  Physical model of radon migration in active fault

    图  2  郯庐断裂带宿迁段地质构造图

    Figure  2.  Geological structure map of the Tanlu fault Suqian segment

    图  3  合欢路场地位置图

    Figure  3.  Location of Hehuan road

    图  4  合欢路测线土壤氡浓度、CO2浓度

    Figure  4.  Soil radon and carbon dioxide survey line of Hehuan road

    图  5  合欢路土壤氡迁移模型

    Figure  5.  Radon migration mode of Hehuan Road

    图  6  F5-2断层土壤氡迁移数值模拟

    Figure  6.  Numerical simulation of soil radon migration in fault F5-2

    图  7  F5-1断层土壤氡迁移数值模拟

    Figure  7.  Numerical simulation of soil radon migration in fault F5-1

    表  1  覆盖层中氡迁移相关参数

    Table  1.   Parameters related to radon migration in overburden

    土壤性质孔隙度e/%扩散系数D/10−2 cm2·s−1有效扩散系数D*/cm2·s−1对流速度V/10−4 cm·s−1
    40.0 4.50~7.00 0.110~0.175 8.0
    疏松沉积物 20.0 2.00~2.50 0.100~0.125 6.8
    白黏土 59.3 1.53 0.023 4.0
    砂质黏土 10.8 1.09 0.100 4.2
    下载: 导出CSV

    表  2  土壤氡测线探测结果与浅层人工地震探测结果

    Table  2.   Detection results of soil radon line and shallow seismic

    土壤氡测线
    分段名称
    测线长度
    /m
    最大值
    /kBq·m−3
    最小值
    /kBq·m−3
    平均值
    /kBq·m−3
    异常下限
    /kBq·m−3
    异常
    形态
    异常峰区
    间间距/m
    异常区间
    /m
    浅层人工
    地震测线
    分支断裂
    编号
    断点位置
    /m
    断错
    性质
    倾向
    倾角/°
    上断点
    埋深/m
    SG13601.091.099.4023.46单峰1490~1545SF5-21508正断E706
    SG23600.990.994.117.77双峰231970~2022SF5-11988正断W6518
    下载: 导出CSV

    表  3  断点上覆盖层参数

    Table  3.   Overburden parameters on fault point

    断点断点上覆盖层参数
    覆盖层深度/m厚度/m描述对流速度/cm·s−1有效扩散系数/cm2·s−1
    断裂F5~2断点 1 0~4 4 灰黑、灰黄色黏土,质软 4.5 0.06
    2 4~10 6 黄棕、青灰色条纹砂质黏土,含有锰结核和姜结石 6.0 0.11
    3 10~24 14 黄棕、青灰色黏土 4.0 0.05
    4 24~32 8 黄棕、灰黄色砂质黏土 5.0 0.08
    断裂F5-1断点 1 0~2.0 2 灰黑色,黑色黏土,质软 4.5 0.06
    2 2.0~3.8 1.8 棕红色,局部有灰色粗砂质黏土,含有钙质结核 6.0 0.11
    3 3.8~4.5 0.7 青灰色粉砂质黏土,质软 5.0 0.08
    4 4.5~5.8 1.3 棕色、褐棕色黏土,含有铁锰结核 4.0 0.04
    下载: 导出CSV
  • [1] 晁洪太, 李家灵, 崔昭文等, 1994. 郯庐断裂带中段全新世活断层的特征滑动行为与特征地震. 内陆地震, 8(4): 297—-304.

    Chao H. T., Li J. L., Cui Z. W., et al., 1994. Characteristic slip behavior of the holocene fault in the central section of the Tanlu fault zone and the characteristic earthquakes. Inland Earthquake, 8(4): 297—304. (in Chinese)
    [2] 车用太, 刘耀炜, 何钄, 2015. 断层带土壤气中H2观测——探索地震短临预报的新途径. 地震, 35(4): 1—10.

    Che Y. T., Liu Y. W., He L., 2015. Hydrogen monitoring in fault zone soil gas—a new approach to short/immediate earthquake prediction. Earthquake, 35(4): 1—10. (in Chinese)
    [3] 戴波, 赵启光, 张敏等, 2020. 土壤氡对郯庐断裂宿迁段F5断裂探测和活动性的研究. 地震工程学报, 42(6): 1479—1486.

    Dai B., Zhao Q. G., Zhang M., et al., 2020. Detection and activity of the fault F5 in Suqian segment of the Tanlu fault by using soil radon. China Earthquake Engineering Journal, 42(6): 1479—1486. (in Chinese)
    [4] 杜建国, 宇文欣, 李圣强等, 1998. 八宝山断裂带逸出氡的地球化学特征及其映震效能. 地震, 18(2): 155—162.

    Du J. G., Yu W. X., Li S. Q., et al., 1998. The geochemical characteristics of escaped radon from the Babaoshan fault zone and its earthquake reflecting effect. Earthquake, 18(2): 155—162. (in Chinese)
    [5] 付晓飞, 方德庆, 吕延防等, 2005. 从断裂带内部结构出发评价断层垂向封闭性的方法. 地球科学——中国地质大学学报, 30(3): 328—336.

    Fu X. F., Fang D. Q., Lü Y. F., et al., 2005. Method of evaluating vertical sealing of faults in terms of the internal structure of fault zones. Earth Science—Journal of China University of Geosciences, 30(3): 328—336. (in Chinese)
    [6] 葛良全, 邹功江, 谷懿等, 2012. 非稳态条件下壤中氡浓度数理模型探讨. 成都理工大学学报(自然科学版), 39(3): 323—327.

    Ge L. Q., Zou G. J., Gu Y., et al., 2012. Research on mathematic-physical model of radon migration in soil under unsteady conditions. Journal of Chengdu University of Technology (Science & Technology Edition), 39(3): 323—327. (in Chinese)
    [7] 柯云龙, 刘耀炜, 张磊等, 2018. 川滇地震预报实验场高精度氢观测台阵建设分析. 地震, 38(3): 35—48.

    Ke Y. L., Liu Y. W., Zhang L., et al., 2018. Establishment and analysis of the high-precision hydrogen observation array in China earthquake science experiment site. Earthquake, 38(3): 35—48. (in Chinese)
    [8] 刘保金, 酆少英, 姬计法等, 2015. 郯庐断裂带中南段的岩石圈精细结构. 地球物理学报, 58(5): 1610—1621. doi: 10.6038/cjg20150513

    Liu B. J., Feng S. Y., Ji J. F., et al., 2015. Fine lithosphere structure beneath the middle-southern segment of the Tan-Lu fault zone. Chinese Journal of Geophysics, 58(5): 1610—1621. (in Chinese) doi: 10.6038/cjg20150513
    [9] 刘洪涛, 2018. 土壤氡迁移数值模拟及土壤氡对流速度的研究. 北京: 中国地质大学(北京).

    Liu H. T., 2018. Study on the numerical simulation of soil radon migration and soil radon convective velocity. Beijing: China University of Geosciences (Beijing). (in Chinese)
    [10] 刘菁华, 2006. 活断层上覆盖层中氡迁移的数值模拟及反演拟合. 长春: 吉林大学.

    Liu J. H., 2006. Numerical simulation, inversion fitting of radon migration in the overburden above active fault. Changchun: Jilin University. (in Chinese)
    [11] 刘耀炜, 任宏微, 2009. 汶川8.0级地震氡观测值震后效应特征初步分析. 地震, 29(1): 121—131. doi: 10.3969/j.issn.1000-3274.2009.01.016

    Liu Y. W., Ren H. W., 2009. Preliminary analysis of the characteristics of post-seismic effect of radon after the Wenchuan 8.0 earthquake. Earthquake, 29(1): 121—131. (in Chinese) doi: 10.3969/j.issn.1000-3274.2009.01.016
    [12] 邵永新, 杨绪连, 李一兵, 2007. 海河隐伏活断层探测中土壤气氡和气汞测量及其结果. 地震地质, 29(3): 627—636.

    Shao Y. X., Yang X. L., Li Y. B., 2007. The result and measurement of soil gas radon and soil gas mercury in the exploration of Haihe hidden fault. Seismology and Geology, 29(3): 627—636. (in Chinese)
    [13] 汪成民, 李宣瑚, 1991. 我国断层气测量在地震科学研究中的应用现状. 中国地震, 7(2): 19—30.

    Wang C. M., Li X. H., 1991. Applications of fracture-gas measurement to the earthquake studies in China. Earthquake Research in China, 7(2): 19—30. (in Chinese)
    [14] 王江, 李营, 陈志, 2017. 口泉断裂断层气地球化学变化特征及断层活动性. 地震, 37(1): 39—51.

    Wang J., Li Y., Chen Z., 2017. Gas geochemistry and activity of the Kouquan fault in Shanxi province. Earthquake, 37(1): 39—51. (in Chinese)
    [15] 吴慧山, 1995. 氡测量方法与应用. 北京: 原子能出版社.

    Wu H. S., 1995. Methods and applications of radon. Beijing: Atomic Energy Press. (in Chinese)
    [16] 伍剑波, 张慧, 苏鹤军, 2014. 断层气氡在不同类型覆盖层中迁移规律的数值模拟. 地震学报, 36(1): 118—128.

    Wu J. B., Zhang H., Su H. J., 2014. Numerical simulation for migration rule of fault gas radon in different overburden. Acta Seismologica Sinica, 36(1): 118—128. (in Chinese)
    [17] 许汉刚, 范小平, 冉勇康等, 2016. 郯庐断裂带宿迁段F5断裂浅层地震勘探新证据. 地震地质, 38(1): 31—43.

    Xu H. G., Fan X. P., Ran Y. K., et al., 2016. New evidences of the Holocene fault in Suqian segment of the Tanlu fault zone discovered by shallow seismic exploration method. Seismology and Geology, 38(1): 31—43. (in Chinese)
    [18] 张慧, 苏鹤军, 李晨桦, 2013. 合作市隐伏断层控制性地球化学探测场地试验. 地震工程学报, 35(3): 618—624.

    Zhang H., Su H. J., Li C. H., 2013. Field test on the geochemical detection of concealed fault in Hezuo city. China Earthquake Engineering Journal, 35(3): 618—624. (in Chinese)
    [19] 周慧玲, 苏鹤军, 张慧等, 2018. 基于孕震物理模式的断层气流动监测网络布设技术. 地震工程学报, 40(5): 1052—1060.

    Zhou H. L., Su H. J., Zhang H., et al., 2018. Mobile monitoring network layout technique for fault gas based on seismogenic mode. China Earthquake Engineering Journal, 40(5): 1052—1060. (in Chinese)
    [20] 周晓成, 王传远, 柴炽章等, 2011. 海原断裂带东南段土壤气体地球化学特征. 地震地质, 33(1): 123—132.

    Zhou X. C., Wang C. Y., Chai Z. Z., et al., 2011. The geochemical characteristics of soil gas in the southeastern part of Haiyuan fault. Seismology and Geology, 33(1): 123—132. (in Chinese)
    [21] Abdoh A., Pilkington M., 1989. Radon emanation studies of the Ile Bizard Fault, Montreal. Geoexploration, 25(4): 341—354. doi: 10.1016/0016-7142(89)90005-7
    [22] Baubron J. C., Rigo A., Toutain J. P., 2002. Soil gas profiles as a tool to characterise active tectonic areas: the Jaut Pass example (Pyrenees, France). Earth and Planetary Science Letters, 196(1—2): 69—81.
    [23] Giammanco S., Gurrieri S., Valenza M., 1998. Anomalous soil CO2 degassing in relation to faults and eruptive fissures on Mount Etna (Sicily, Italy). Bulletin of Volcanology, 60(4): 252—259. doi: 10.1007/s004450050231
    [24] Ho C. K., 2008. Analytical risk-based model of gaseous and liquid-phase radon transport in landfills with radium sources. Environmental Modelling & Software, 23(9): 1163—1170.
    [25] Iakovleva V. S., Ryzhakova N. K., 2003. Spatial and temporal variations of radon concentration in soil air. Radiation Measurements, 36(1—6): 385—388.
    [26] Janik M., Bossew P., 2016. Analysis of simultaneous time series of indoor, outdoor and soil air radon concentrations, meteorological and seismic data. Nukleonika, 61(3): 295—302. doi: 10.1515/nuka-2016-0049
    [27] King C. Y., King B. S., Evans W. C., et al., 1996. Spatial radon anomalies on active faults in California. Applied Geochemistry, 11(4): 497—510. doi: 10.1016/0883-2927(96)00003-0
    [28] Lombardi S., Voltattorni N., 2010. Rn, He and CO2 soil gas geochemistry for the study of active and inactive faults. Applied Geochemistry, 25(8): 1206—1220. doi: 10.1016/j.apgeochem.2010.05.006
    [29] Toutain J. P., Baubron J. C., 1999. Gas geochemistry and seismotectonics: a review. Tectonophysics, 304(1—2): 1—27.
    [30] Wakita H., Nakamura Y., Notsu K., et al., 1980. Radon anomaly: a possible precursor of the 1978 Izu-Oshima-Kinkai earthquake. Science, 207(4433): 882—883. doi: 10.1126/science.207.4433.882
    [31] Walia V., Yang T. F., Lin S. J., et al., 2013. Temporal variation of soil gas compositions for earthquake surveillance in Taiwan. Radiation Measurements, 50: 154—159. doi: 10.1016/j.radmeas.2012.11.007
    [32] Yakovleva V. S., Parovik R. I., 2010. Solution of diffusion–advection equation of radon transport in many-layered geological media. Nukleonika, 55(4): 601—606.
    [33] Zhang W., Zhang D. S., Wang X. F., et al., 2014. Analysis of mathematical model for migration law of radon in underground multilayer strata. Mathematical Problems in Engineering, 2014: 250852.
  • 加载中
图(7) / 表(3)
计量
  • 文章访问数:  147
  • HTML全文浏览量:  17
  • PDF下载量:  15
  • 被引次数: 0
出版历程
  • 收稿日期:  2020-07-07
  • 网络出版日期:  2021-07-12
  • 刊出日期:  2021-03-01

目录

/

返回文章
返回