地质构造变形数值模拟研究的原理、方法及相关进展.pdf

第20卷第5期 2005 年5月 地球科学进展 ADVANCESIN EARTH SCIENCE Vol .20 No.5 May.，2005 文章编号 1001- 8166 （2005 ）05- 0549- 07 地质构造变形数值模拟研究的原理、 方法及相关进展 * 林 舸 1，赵重斌1，2 ，张晏华 1，2，王岳军1，单业华1，3 （1.中国科学院广州地球化学研究所， 中国科学院边缘海地质重点实验室， 广东 广州 510640 ； 2. CSIRO Exploration ＆ Mining ， Australia Perth W A 6151 ； 3.中南大学地学与环境工程学院计算地球科学研究中心， 湖南 长沙 430012 ） 摘 要 地质构造变形数值模拟实验研究是遵循自然科学规律， 依据数学、 物理学、 化学等学科中已 为人们认识的规律为原理； 依靠计算机先进的综合处理系统； 以固体地球科学资料为基础所构建的 地质模型为实验研究本体。因而， 地质构造变形数值模拟是采用数学物理方法对相关地质学问题 的科学描述， 并采用理论分析和数值模拟实验方法， 对所描述的相关地质学过程进行定量化求解。 由于数值模拟可以综合利用其它学科的研究成果， 建立和模拟不受时空限制的各种地质模型， 为定 量化地解决一些复杂的地质学问题， 提供有效的分析工具， 已经引起了国内外地学界的重视和深入 研究。数值模拟实验研究方法的发展， 将极大地促进我国地学界从采用了多年的传统的描述性研 究方法向科学的定量化研究方法的转变。 关 键 词 地质； 构造变形； 数值模拟； 原理方法 中图分类号 P54 文献标识码A 0 引 言 现今获取和观察到的大量地质信息和实例， 都 是地质构造演化过程中某一幕或某一阶段的地质表 现， 是地质构造演化过程最后一幕的综合结果或某 一次地质事件的具体产物。然而， 早期的地质信息 或地质作用所发生的过程可能被后期地质作用所改 造、 叠加或抹失。因此， 完整地观察或再现地质构造 变形演化过程， 热动力驱动和构造变形之间的耦 合和叠加关系， 俯冲和碰撞带的形成与演化发展， 造 山带伸展垮塌与陆内裂解作用等事件的演变过程及 其内在动力机制显得尤为困难。 传统的地质地球物理地球化学理论和分析 方法， 无疑是再现地质构造变形演化过程的重要手 段之一， 并已取得长足的发展； 由于不同的学者从不 同的侧重点来研究地质构造演化的动力学机制， 以 及对地质过程理解方面的局限性， 这些手段在研究 地球动力学机制方面难免存在一定程度的局限性。 物理模拟实验可以帮助人们理解构造变形过程和动 力学的作用过程， 但也存在严重的时空尺度的局限 性。数值模拟实验研究可以综合利用其它方法 （如 地质、 地球物理、 地球化学 ）的研究成果， 建立和模 拟不受时空限制的各种地质模型， 而成为现代地球 科学研究的重要方法之一。通过地质演化过程的正 演模拟， 重建地球或一个区域的地质演化过程， 从而 得到深部动力学与浅部构造响应的重要信息， 以帮 助人们认识与理解地质构造演化过程。通过数值模 拟实验研究， 可以理解和认识地质历史时期热动力 作用与构造变形之间的耦合与叠置关系， 地质历史 时期内岩石圈构造、 热结构、 应力应变场的演变， 岩 石圈内部不同层次或尺度物质能量的传输、 转换 及其调整过程， 岩石圈演变与重要构造形迹的形成 * 收稿日期2004- 03-08 ； 修回日期2004- 09-28. *基金项目 国家重点基础研究发展规划项目 “主碰撞带地壳挤压增厚和缩短阶段的构造变形和成矿响应” （编号2002 CB412601 ）； 中 国科学院广州地球化学研究所所长基金项目 “壳幔过渡带的形成及揭示的动力学模式 ” （编号030513 ）联合资助. 作者简介 林舸 （1948-），男，湖南洞口人，研究员，主要从事大地构造与成矿学研究. E- m ail gelin ＠ gig. ac. cn 过程及其形成机理等。 构造变形数值模拟实验研究的发展主要在于 ①它是一种大学科的交叉研究和综合分析研究方 法， 所涉及到的主要学科包括地质学、 地球物理、 地 球化学、 数学、 计算力学和计算机技术等。科学上的 重大突破， 新的生长点乃至新学科的产生，往往在 相邻学科彼此交叉和相互渗透的过程中形成。特别 是地学中 “以问题为导向”的研究更体现了这一 点 ［1］。②将常规的广为大地构造学或地质演化过 程研究所采用的归纳、 演绎的逻辑推测方法， 向通过 由计算机算法和技术表示的数学逻辑方法推进， 从 而为综合分析地质资料提供有效的研究方法； 同时， 也为实现地球动力学研究的创新从理论上和方法上 提供可靠的支持。 1 地质构造变形数值模拟的原理与方法 地质构造变形数值模拟实验的基本原理， 是遵 循自然科学规律， 如物质不灭定律、 能量守恒定律 等； 并依据数学、 物理学、 化学等学科中已为人们认 识的基本规律为原理； 依靠计算机先进的综合处理 系统； 并以固体地球科学资料为基础所构建的地质 模型为实验研究本体。地质构造变形数值模拟实验 研究， 不是单纯的数据处理和图像分析， 它是采用数 学物理方法对相关地质学问题的科学描述， 并采用 理论分析和数值模拟实验方法， 对所描述的相关地 质过程进行定量化的求解。因而， 对能充分反映岩 石圈及地壳变形的基本物理化学过程 （如流体流动 过程、 热传递过程、 介质变形过程、 热液流体的混合 化学反应过程和水岩化学反应过程等 ）及其相关的 基本科学问题的认识和知识积累， 是探讨地质构造 变形等复杂地质学问题的基本前提条件。一般地 说， 地质构造变形数值模拟实验研究先从简单模型 做起， 简单的几何形态是数值模拟常采用的方法， 能 有效地认识一些基本的科学问题， 如温度场的改变， 导致流体的对流及含矿热液系统化学反应的发生， 再现成矿物质的溶解与沉淀过程。由于地质学问题 的复杂性， 既使对这些具有简单几何构形， 但能充分 反映影响岩石圈及地壳变形的基本问题的求解， 采 用纯数学分析的方法常常也是不太可能的。因此， 研究和发展高级数值计算方法， 及数值模拟技能的 积累无疑也是为解答上述的这些基本科学问题进行 求解和探讨地质构造变形而不可缺少的研究内容。 为了探讨地质构造变形， 必须深入研究与岩石 圈及地壳变形机制相关的基本科学问题。如孔隙流 体在岩石圈及地壳中的对流， 这一对流作用可导致 类似地幔柱式的地温重分布； 由于组成岩石圈及地 壳的材料强度随温度变化而变化， 而这一变化又可 导致组成岩石圈及地壳材料的屈服破坏。因而， 这 种地温的重分布无疑又对地壳的构造变形有着极为 重要的影响。又如在岩石圈及地壳中发生的不同流 体的混合化学反应及水岩化学反应， 这些化学反应 可导致岩石圈及地壳中骨架材料的溶解或沉积， 这 种由于化学反应引起的溶解或沉积作用， 又将改变 岩石圈及地壳中骨架材料的结构及孔隙分布， 从而 明显地改变岩石圈及地壳材料的强度及传热特性， 无疑又将导致组成岩石圈及地壳材料的屈服破坏， 结果也是严重地影响岩石圈及地壳的构造变形。由 此可见， 为了深入地研究和探讨构造变形问题， 我们 必须对诸如上述的一些基本科学问题进行深入的探 讨。为了实现这一目标， 我们需要对一些能反映这 些基本科学问题的概念和理论模型进行详细的 研究。 构造变形数值模拟实验研究的研究方法， 主要 包括三大部分 一是建立科学的高级数值计算算法， 这是数值模拟实验研究的核心。只有依据现已为人 们所发现和认识的自然科学规律建立的数学物理方 程组， 才有可能定量化模拟相关的地质学过程， 并逐 步将地质科学纳入定量的科学理论范畴。由于问题 的复杂性， 有限元和有限差分是解决这类耦合问题 最有效的方法； 二是将建立的高级数值计算算法转 换成计算机语言， 这一部分工作将主要由计算机方 面的专家来做； 三是依据已有的地质资料， 以观察到 的地质现象为依据， 选择参数进行数值模拟实验 研究。 依据观测或推测的边界条件、 实验参数， 来建立 实验研究的初始条件、 边界条件， 并不断修正、 更改 实验参数、 初始条件、 边界条件进行模拟实验， 逐步 使实验结果趋近于或等于观察到或科学推测的地质 构造模型， 从数学理论和方法上证实或证伪已有的 地质模型， 建立符合区域地质实况的地质构造模式。 因此， 综合分析已有资料， 构建合理的地质模型， 从 不同方面选择实验研究的初始条件、 边界条件及实 验参数， 并选择不同的数学模型进行实验研究， 并从 地质、 地球物理、 地球化学等方面提供约束条件， 限 制实验结果的多解性， 并使实验结果趋近于地质现 状， 是数值模拟的主要研究内容与方法。 2 国外地质构造数值模拟研究的现状 基于数值模拟研究的优越性， 近20年， 那些已 055 地球科学进展 第20卷 在传统构造地质学、 成矿学、 大地构造学研究领域卓 有成就的专家、 学者， 如澳大利亚的 Bruce Hobbs ， 英国的 England ， Housem an ， 德国的 Cee Pashan ，加 拿大的 Paul W illiam s ， 瑞士的 Neil Mancktelow 等， 都 已致力于数值模拟方面的研究； 目前许多国家 （如 美国、 加拿大、 英国、 法国、 德国、 日本、 瑞士、 澳大利 亚 ） 也都相继建立起计算机模拟研究实验室或研究 中心。国外研究者相继利用数值模拟手段研究了碰 撞与碰撞后过程深部动力学机制、 浅部响应和远程 地质构造间的热应力传输与耗损过程、 板块俯冲对 俯冲带热动力的调整及其对岩浆作用源区性质、 俯 冲带形变特征及大陆增生等的影响、 大陆俯冲陆 陆碰撞过程下插板片断离作用的成因机制及其对同 碰撞和碰撞后岩浆作用与大陆构造变形的影响与制 约作用、 陆内岩浆作用与伸展构造、 多陆块相互作用 与侧向构造逃逸间耦合作用机制、 岩石圈变形温 度流体相互作用机制等； 如 Bott ［2］模拟了板缘力 驱动下板内应力应变的演化； Tapponnier ［3］对板块俯 冲过程及俯冲带中变形的研究； Mckenzie ［4］对大陆 引张和盆地形成过程的模拟； England 等［5］在1982 年首次将黏性薄层流变模型运用到以印度板块和欧 亚板块的碰撞为代表的大陆形变， 在一矩形区域内 对青藏高原的挤压隆升演化进行了数值模拟； Housem an 等 ［6～8 ］成功模拟了印度板块和欧亚板块 碰撞后板块边缘变形与地幔对流之间的耦合关系， 他们的实验研究结果表明 地壳厚度变化是调整碰 撞动力学环境变化的主要地质构造效应， 走滑剪切 断裂只是部分地缓解了这种应力场的变化， 地壳厚 度变化产生的向上浮力， 导致陆块间的侧向运行和 地壳厚度变化。 Zhang 等 ［9，10］用模拟手段成功地解 释了澳大利亚东部的岩石圈构造和异常应力场的成 因 澳大利亚东部被动板缘或陆/洋壳交接带， 有一 高密度区指向大陆的上地幔契， 正是这一高密度区 的存在， 导致了澳大利亚东部的岩石圈构造和异常 应力场的发生和地壳构造运动。 Zhao 等 ［11］用变 形流体热力学化学反应全耦合的模拟方法实 现对构造控矿理论的研究和找矿勘探的应用。 国际上近几年的研究表明， 构造变形动力学数 值模拟研究的方法是成熟的。目前的发展主要在数 学公式的推导与高级计算机算法的设计， 并主要利 用有限元和有限差分方法。模拟研究也充分考虑到 岩石圈的特性， 加强了在流变学理论方面的研究； 更 大程度上考虑到实际地质现况。近几年， 国外研究 者对能充分反映岩石圈及地壳变形的基本物理化学 过程 （如流体流动过程， 热传递过程， 介质变形过 程， 热液流体的混合化学反应过程和水岩化学反应 过程等 ）及其相关的基本地质科学问题的认识和发 展高级数值计算方法， 进行了初步研究。澳大利亚 的 Zhao 等 ［12］在近几年来， 已研究和发展了可用于 求解岩石圈及地壳中介质变形过程， 流体流动过程， 热传递过程， 热液混合化学反应和水岩化学反应过 程的一些高级计算机算法， 利用简单的几何模型进 行了一些有效的研究工作， 为定量化地解决许多复 杂的地质学问题， 提供了新的思路和研究方法。 3 地质构造变形数值模拟研究在我国 的进展 近几年， 我国地质学者采用不同的方法， 在构造 变形动力学数值模拟研究方面做了一些工作， 并有 一批成果报道 ［13～25 ］。其中， 宁杰远等［13］对俯冲带 波速结构的数值模拟；毛兴华等 ［14］根据弹性本构 关系， 用平面应变有限元方法计算了相变引起的体 积变化、 板块内部温度差、 密度异常及边界力产生的 应力场分布情况， 对俯冲带深部应力场开展了二维 弹性有限元数值模拟； 张贵宾等 ［15］利用中英合作完 成的有限元粘弹塑构造数值模拟技术 FEVPLIB ， 模 拟研究了青藏高原西部横过西昆仑和塔里木结合带 剖面的这一动力学演化过程；何丽娟 ［16］采用粘弹 性动力学模型模拟研究了岩石圈流变性对拉张盆地 构造热演化的影响。李祖宁等 ［17］将印度板块持续 地向北推进、 下伏地幔小尺度对流对增厚大陆岩石 层的搬离作用以及剥蚀作用视为形成现今东亚大陆 形变和应力场格局的主驱动力， 在一梯形区域内， 利用数值模拟的方法，研究了东亚大陆在不同的边 界条件、 不同的剥蚀率系数及不同的岩石力学参数 条件下的形变及应力场格局。青藏高原的隆升及其 动力学机制， 一直是地球科学家关注的热门课题。 李祖宁等 ［18］修改了 England 和 Mckenzie 的黏性薄 层流变模型中控制大陆形变的连续性方程，将剥蚀 作用对高原隆升演化的影响直接引入该方程，并考 虑下伏地幔小尺度对流对增厚岩石层的搬离作用对 高原隆升演化后期的影响， 用有限差分法直接模拟 青藏高原隆升过程。数值模拟结果所显示的高原隆 升演化过程与实际观测资料吻合较好， 揭示了高原 隆升演化过程的非平稳和多阶段的特性。毕思 文 ［19］根据山地系统岩石圈运动的主要因素构造 带和断裂带， 应用变形体模拟构造带，采用摩擦模 型模拟断裂带， 并根据周边位移情况给出边界条件。 155 第5期 林 舸等 地质构造变形数值模拟研究的原理、 方法及相关进展 应用 ABAQUS 有限元分析软件和位移加载， 开展了 青藏高原山地系统岩石圈应力位移场的数值模拟 分析， 给出了应力和位移的分布规律。 中国科学院广州地球化学研究所地质构造变形 数值模拟实验研究小组， 近年来与国外学者进行合 作研究， 初步地研究了岩石圈与地壳中流体对流的 可能性， 揭示了流体对流在岩石圈及地壳中发生的 基本前提条件。尽管地学界通过岩石探针，从物质 成分分析对比， 普遍认为在岩石圈板块中存在着大 规模的对流。但是并不完全清楚在不同的边界条件 下， 对流在岩石圈及地壳中发生的基本前提条件。 通过对岩石圈及地壳不同尺度的流体对流过程的理 论分析和数值模拟研究， 已认识到 岩石圈与地壳中 流体对流对岩石圈及地壳中的应力分布及其演化过 程有着极为重要的影响。岩石圈及地壳在不同的边 界条件下， 本身具有一个 “系统临界地温温度梯 度 ” ， 当岩石圈及地壳中的实际地温温度梯度大于 或等于这一相关的 “系统临界地温温度梯度 ”时， 对 流就可能在岩石圈及地壳中发生 ［26～37 ］。这一科学 认识从根本上揭示了流体对流在岩石圈及地壳中发 生的基本前提条件， 对探讨华北中生代岩石圈减薄 的动力学机制有着非常重要的理论指导意义。 地质构造变形数值模拟实验的主要内容之一是 研究构造应力与构造变形， 即应力与应变的关系问 题。因此， 开展构造地质基础的数值模拟实验， 是地 质构造变形数值模拟实验的主要内容之一。同时， 数值模拟不但可以在华北、 华南陆块这类板块尺度 大区域的实验研究； 而且对矿区这类小尺度的实验 研究同样是有效的研究方法。在构造地质学的基础 理论方面， 数值模拟还可重点围绕着应力反演 （ stress inversion ）和构造复原 （ restoration ）， 以及褶 皱/断裂构造导致流体的运移、 集中、 混合和沉淀富 集开展有关的创新性研究。 应力反演法是利用断层数据来反演当时原地古 应力状态的一种方法。目前已提出的多种方法都带 有不同程度的局限性， 表现为它们一般不能够有效 地区分开不同构造期的断层数据， 从而导致错误的 应力估计。这种局限性的产生原因主要在于应力反 演问题的非线性。然而， 经过简单的非线性变换， 在 Fry ［38］提出的特征空间 （ sigm aspace ）里， 应力反演 问题可以转变为线性问题， 即同构造期断层数据的 分布呈现出统一的超平面， 而不同构造期断层数据 的分布呈现出不同的超平面。在此基础上， 引入模 糊线性聚类方法， 不仅考虑到数据测量误差给应力 估计带来的不确定性， 而且还能够根据观察数据的 内在特征合理地确定出最佳的构造期数目 ［39～43 ］。 构造复原或恢复是构造地质学中研究构造变形场的 一种重要手段。通常的做法属于不连续恢复， 即先 将离散的有限单元去应变， 然后再拼贴所有的去应 变单元， 最终得到复原的结果。这种人为的分步计 算的不足在于复原结果往往取决于所采用的拼贴方 案， 要求输入较多的控制参数， 且算法不易编程实 现。基于有限元分析， 提出的连续恢复算法， 能够较 好地克服以上的不足之处。目前， 已将它运用到裂 谷盆地三维断裂构造面的复原， 探讨花岗岩侵入体 的侵位空间和侵位机制 ［44～46 ］。小构造 （褶皱、 断 裂、 节理等 ）的活动， 与成矿过程具有密切的联系。 对褶皱/断裂活动导致的变形场的模拟显示， 褶皱/ 断裂活动造成的位移场均可导致形应变和体应变 （岩石体积增加 ） 在褶皱轴部或断裂尖端附近集中， 并进一步引起岩体渗透率的增加和流场流压的降 低， 导致成矿流体在这些位置运移、 集中、 混合， 沉淀 富集； 构造运动及其对流体成岩/成矿耦合作用的 数值模拟结果显示， 在不同的构造环境下， 成矿流体 的流速、 流线和温度分布， 成矿元素的溶解和沉淀类 型及分布区域存在着明显的差异， 成矿流体的混合 反应是成矿作用的重要机制之一 ［47～49 ］。 研究组分别运用 FLAC 软件和薄席片粘性模 型， 研究了可能导致华北中生代岩石圈减薄的动力 学过程。研究结果显示， 在块体挤压应力引起的岩 石圈/地壳加厚背景下， 重力伸展垮塌作用是导致 华北陆块盆山体系形成的重要因素。热应力模拟 表明， 热侵蚀作用是导致厚且冷的克拉通岩石圈减 薄的主导因素， 而机械伸展引起的减薄幅度有限； 华 北陆块中生代岩石圈减薄的过程是地热热流与应力 共同作用的结果 ［50，51］。该项成果的创新主要在3 个方面 ①地球热力学模拟本构方程的建立。本项 研究建立的模型全面考虑了岩石圈热传导特性， 岩 石放射性生热的能力以及源自地幔深处的热流的影 响。这在国际上现时进行的地球热力学研究中， 是 同行公认的较全面的数值建模方案。引导出的模型 能真实地反映岩石圈的热演化。②本项研究对热侵 蚀减薄和机械构造伸展减薄进行了综合考虑和耦合 模拟。国际上类似的模型多将两种机制单独考虑。 因此， 本模型较真实地再造大地构造事件的实际过 程。③对华北岩石圈的减薄过程进行了初步的定量 描述。这种定量的结论代表着对推理性的、 模糊的 华北热构造事件常规描述的进一步的认识。 255 地球科学进展 第20卷 4 问题与讨论 我国在大陆地质构造变形动力学数值模拟研究 方面的工作虽已经开始， 但只能说还在探索阶段。 在碰撞与碰撞后过程深部动力学机制、 俯冲带形变 特征及大陆增生等的影响、 大陆俯冲陆陆碰撞过 程下插板片断离作用的成因机制及其对同碰撞和碰 撞后岩浆作用与大陆构造变形的影响与制约作用、 陆内岩浆作用与伸展构造、 多陆块相互作用与侧向 构造逃逸间耦合作用机制等数值模拟研究方面， 尚 未开展深入研究。一方面是我们起步晚， 另一方面 是尚未能真正认识到数值模拟研究的重要性及其可 能的前景。 由于地质学问题的复杂性， 要对哪些具有简单 几何构形， 但能充分反映影响岩石圈及地壳变形的 基本科学问题的求解， 必须研究和发展高级数值计 算方法。这就需要大学科的交叉和科学家们的合 作， 特别需要数学家， 物理学家和化学家的参与。如 果地质学家想用一组数理方程式来描述岩石圈及地 壳变形的过程， 没有数学家， 物理学家和化学家的参 与几乎是不可能的。有了他们的参与， 加上地质学 家对地质科学资料的积累和认识及其正确的逻辑思 维方法， 利用的现代高科技手段， 研究和发展一些高 级计算机算法， 为定量化地解决许多复杂的地质学 问题， 甚至用一组数理方程式来描述岩石圈及地壳 变形的过程， 也不是不可能的。 参考文献 （ References ） ［1］ Sun Chengquan ， Qu J iansheng. Status and trends ofthe interna - tional Earthsciencesstudies ［J］. Advancesin Earth Science ， 2002 ，17（3）344- 347.［孙成权，曲建升.国际地球科学发展 态势 ［J］.地球科学进展，2002 ，17（3） 344- 347.］ ［2］ Bott M H ， Dean D S. Stresssystem s atyoung continental m argins ［J］. Nature ，1972 ，235 23- 25. ［3］ Tapponnier P. Slip- line fieldtheory and lar ge- scale continental tectonics ［J］. Nature ，1976 ，264 319- 324. ［4］ Mc Kenzie D P. Som erem arks onthe m ovem entofs m allm eltfrac- tionsinthe m antle ［J］. Earth Planet Sciences Letter ，1989 ，95 53- 72. ［5］ England P C ， Mckenzie D A . Thin viscoussheet m odelfor cont i- nental de ation ［J］. J ournal of Geophysics Research ，1982 ， 70295- 321. ［6］ Philip England ， Gregory Housem an. Role of Lithosphericstren gth heterogeneitiesinthe tectonics of Tibet and neighbouring regions ［J］. Nature ，1985 ，315 （6 017 ）297-301. ［7］ England P ， Housem an G A . Extension during continental co nver- gence ， with applicationtothe Tibetan Plateau ［J］. J ournalof Ge- ophysics Rresearch ，1989 ，9417 561- 17 579. ［8］ Gregory Housem an ， Philip England. A Crustal - Thickening Model forthe Ind- Asian Collision ［M］. Los Angeles ， USA ， the Rubey Colloquium - Tectonic Evolution of Asia ，1994. ［9］ Zhang Y ， Scheibner E ， Ord A ， etal . Num ericalm odelling ofcrus- talstressesinthe eastern Australian passiv e m argin ［J］. Austral- ian J ournal of Earth Sciences ，1996 ，43161- 175. ［10］ Zhang Y ， Scheibner E ， Hobbs B E ， etal . Lithosphericstructure in SE Australia A m odelbased on gravity ， geoid and m echanical m odeling ［A］. In Braun J， Dooley J， Goleby B ， et al ， eds. Structure and Dynamics ofthe Australianlith osphere ［C］. AGU Monograph ， Geodynamics Series ，1998 ，2689- 108. ［11］ Zhao Chongbin ， Hobbs B E ， Mhlhaus H B . Finite elem entm od- elling oftem perature gradient Drivenrock al teration and minerali- zationin porousrock m asses ［J］. Com puter s in Applied Mechanicsand Engineering ，1998 ，165 175- 187. ［12］ Zhao Chongbin ， Hobbs B E ， Mhlhaus HB ， et al . Num erical m odelling ofdouble Diffusion drivenreactiv eflowtransportin de- ablefluid- saturated porous m edia with p articular considera- tion oftem perature- dependentchemicalreac tionrates ［J］. Inter- national J ournalfor Com puter - Aided Enginee ring and Software Engineering Com putations ，2000 ，17367-385. ［13］ Ning J ieyuan ， Zang shaoxian. Num erical Sim ulation of P- w av e velocitystructure ofsubductionzones ［J］. Chinese J ournalof Ge- ophysics ，2001 ，44（2）190- 198.［宁杰远，臧绍先.俯冲带波 速结构的数值模 拟 ［J］.地球 物理学 报，2001 ，44（2）190- 198.］ ［14］ Mao xinhua ， Liu Yaqing ， Ye Guoyang ， et al. 2- D elastic FEM sim ulation onstressstateinthe deep part of a subductedslab ［J］. Acta Seism ologica Sinica ，2002 ，24（3）277-284.［毛兴 华， 刘亚静， 叶国扬， 等.俯冲带深部应力场的二维弹性有限 元数值模拟 ［J］.地震研究，2002 ，24（3）277-284.］ ［15］ Zhang Guibin ， Martin H PBott . Modelling the evlution ofa hal f graben basin using elasto- viscoplasticfini te elem entanalysis ［J］. Earth Science Frontiers ，2000 ，7（4）442- 448.［张贵宾， Martin HPBott .半地堑盆地演化机制的粘弹塑数 值模拟 ［J］.地学 前缘，2000 ，7（4）442-448.］ ［16］ He Li juan. Effects oflithosphericrheology ontherm al - m echanical m odelling ofextension basins ［J］. Chinese J ournalof Geophysics ， 2002 ，45（1）49- 55.［何丽娟.岩石圈流变性对拉张盆地构造 热演化历 史的 影响［J］.地 球物 理学 报，2002 ，45（1）49- 55.］ ［17］ Li Zuning ， Fu Rongshan ， Zheng Yong ， et al. Evolution ofthe de ation and stress patterns ofthe East As ia continent under m ulti- drivingforce ［J］. Acta Seism ologica Sinica ，2002 ，24（1） 17- 26.［李祖宁， 傅容珊， 郑勇， 等.多种驱动力作用下东亚大 陆形变及应力场演化 ［J］.地震学报，2002 ，24（1）17- 26.］ ［18］ Li Zuning ， Fu Rongshan ， Huang J ianhua. Num ericalsim ulation ofthe Qinghai- Xizang Plateau uplift underth e effectof denuda- tion and m antle convection ［J］. Chinese J ournal of Geophysics ， 2002 ，45（4）516-524.［李 祖宁， 傅容珊， 黄建华.剥蚀及地 幔作用下青藏高原隆升过程的数值模拟 ［J］.地球物理学报， 355 第5期 林 舸等 地质构造变形数值模拟研究的原理、 方法及相关进展 2002 ，45（4）516-524.］ ［19］ Bi Siwen. Study on dynamic num ericalsim ulat ion of m ountain systemin Tibet Plateau ［J］.J ournalof Mountain Science ，2001 ， 19（4）289- 298.［毕思文.青藏高原山地系统动力学数值模 拟研究 ［J］.山地学报，2001 ，19（4）289-298.］ ［20］ Fu Rongshan ， Xu Yaomin ， Huang J ianhua ， et al. Num erical sim ulation ofthe com pression uplift ofthe Qi nghai - Xizang Plat- eau ［J］. Chinese J ournalof Geophysics ，2000 ，43（3）346- 355. ［傅容珊， 徐耀民， 黄建华， 等.青藏高原挤 压隆升过程的数 值模拟 ［J］.地球物理学报， 2000 ，43（3）346-355.］ ［21］ Fu Rongshan ， Huang J ianhua ， Xu Yaoming ， et al . Num erical sim ulation ofthe collision betw een Indian an d Eurasian plates and the de ations ofthe present Chinese conti nent ［J］. Acta Seis- m ologica Sinica ，2000 ，22（1）1-7.［傅容 珊，黄建华，徐耀 民，等.印度与欧亚板块碰撞的数值模拟和现代中国大陆形 变 ［J］.地震学报，2000 ，22（1）1- 7.］ ［22］ Zhu J ieshou. Structure oflower m antle and co re- m antle boundary region and its geodynamics ［J］. Advancesin Earth Science ， 2000 ，15（2）139- 142.［朱介寿.下地幔及 核幔边结构及地 球动力学 ［J］.地球科学进展，2000 ，15（2）139- 142.］ ［23］ He Li juan. Advanceintectono- therm al m odel ling ofsedim entary basins ［J］. Advancesin Earth Science ，2000 ，15（6）661- 665. ［何丽娟.沉积盆地构造热演化模拟的研究进展 ［J］.地球 科学进展，2000 ，15（6）661-665.］ ［24］ Zhu J ieshou ， Cao J iamin ， Cai Xuelin ， et al. Study forthree- di- m ensionalstructure of Earthinterior and geo dynamicsin China and ad jancentand searegions ［J］. Advancesin Earth Science ， 2003 ，18（4）497- 503.［朱介寿，曹家敏，蔡学林，等.中国 及邻区陆域海域地球内部三维结构及动力学 研究［J］.地球 科学进展，2003 ，18（4）497- 503.］ ［25］ Zheng Guizhou ， Shen Yongli.3- D analysys of geological chara c- teristics and status research of 3- Dgeology m odelling ［J］. Ad- vancesin Earth Science ，2004 ，19（2）218-223.［郑贵州，申永 利.地质特征三维 分析及三维模拟现 状研究 ［J］.地球科学 进展，2004 ，19（2）218-223.］ ［26］ Chongbin Zhao ， Hobbs B E ， Mhlhaus H B ， etal. Finite elem ent m odelling ofthree- dim ensionalconvection p roblem sinfluid- satu- rated porous m edia heatedfrombelow ［J］. Com m unications in Numerical sin Engineering ，2001 ，17101- 114. ［27］ Chongbin Zhao ， Ge Lin ， Hobbs B E ， etal . Finite elem entm od- elling ofheat Transf