河西走廊系列盆地构造演化的三维数值模拟*

李蔚琳，程惠红，张 怀，石耀霖

(中国科学院大学 中国科学院计算地球动力学重点实验室, 北京 100049)

河西走廊地区地处青藏高原东北缘，位于阿拉善块体与北祁连加里东褶皱带之间，是青藏高原向东北方向推挤过程中形成的一个前陆盆地系[1]，也是青藏高原东北缘的主要汇聚带之一，见图1。一系列重大的地球物理场、地形、地貌均以此为界[2]，其活动构造过程记录了青藏高原东变形特征，是研究新生盆地构造演化的理想地区，且对认识青藏高原的扩展变形过程和力学机制有着重要意义。

蓝色曲线为北祁连山地区河流水系；①玉门盆地；②酒泉盆地；③张掖盆地；④武威盆地。图1 河西走廊及其邻区构造示意图Fig.1 Schematic map of the tectonics of the Hexi Corridor and its adjacent areas

目前，国内外对青藏高原东北缘活动构造的发育、形成和生长扩展过程已进行了大量的研究。在受印度板块与欧亚板块碰撞后向北的推挤过程中，青藏高原东北缘的构造发育基本上遵循由南向北扩展的过程，不同阶段多区域断裂经历发育、生长和发生性质转换[3-6]。且该区域的盆地演化自西向东逐渐变年轻[7-8]，通过前陆盆地向背驮式盆地以及外流水系向内流水系的转变，逐渐形成现今的青藏高原东北缘的构造格局[9-10]。此外，随着“中国地壳运动观测网络工程”项目实施，许多学者根据1999年和2001年GPS观测结果对青藏高原东北缘地壳水平状态进行研究[11-12]，得出青藏高原东北缘的西部(主要包括柴达木盆地及其北侧的祁连山地区)，构造运动特征以北东向挤压运动、地壳南北向缩短为主；而青藏高原东北缘的东部(主要包括青海东部、海原断裂带以南及六盘山断裂带以西地区)，构造变形特征以整体顺时针旋转为主。

近几十年活动构造运动特征和定量研究结果表明，青藏高原东北缘发育有近东西向-北东东向的大型左旋走滑断裂带(海原断裂断裂带、阿尔金断裂、祁连山北缘断裂带等)和北西西向的逆冲断裂带(河西走廊内部及边缘断裂、祁连山内部及边缘断裂等)。在这些主要断裂带控制下，一方面，青藏高原东北缘的祁连山呈现出山前的褶皱逆冲向北部前陆方向不断扩展生长的特征；另一方面，河西走廊地区内被一组活动性很强的北北西向的断裂切割成呈一字形排列的隆起与拗陷相间的构造格局，形成盆山相间的地貌，走廊内的系列盆地自东向西依次为武威—张掖—酒泉—玉门盆地[6, 13-14]。通过大量的野外地质考察，郑文俊[14]很好地给出该区域发展过程：30 Ma以前，祁连山东西端为大面积沉积区且与其北部的河西走廊连接行成大片宽广区域；30～20 Ma时，祁连山及河西走廊地区变形较为缓慢；大约12 Ma以前，祁连山地区山体自西向东基本连成一片，且形成一些规模较小的断裂，与此同时祁连山北部的河西走廊开始成形，且河西走廊北部山体也有了一些轻微的变形；近8 Ma时，海原断裂已经完全贯通，河西走廊前陆盆地形成，且分别位于河西走廊南北两端的两条断裂带(祁连山北缘断裂带和走廊北缘断裂带)已经形成；随后走廊内部进一步发育出一组NNW向的山间隆起，将河西走廊断陷带分隔成大小不等的4个盆地，而祁连山与河西走廊盆地内部继续发生褶皱变形。

据此，在前人已有的工作基础上，本文尝试采用三维黏弹塑性有限元数值模拟计算方法，以河西走廊及其邻区近5 Ma的构造环境为基础，综合考虑区域地质构造单元、地层结构、主要断裂带等因素，以期定量给出河西走廊及其邻区近5 Ma内的构造变形及演化过程，刻画河西走廊两侧断裂带的发展过程和系列盆地演化过程，进而定量分析探讨走廊内系列盆地形成的力学机制，对青藏高原东北缘的横向扩展过程和机理进行讨论。

1 数值计算模型

1.1 模型设置

基于青藏高原东北缘晚新生代的构造环境，综合考虑地质构造单元的相对完整性和区域断裂、盆地空间分布信息[2, 14-27]，结合GPS水平速度场[28]，搭建河西走廊及其邻区的模型(35.8～42.1°N、94.5～107°E)，见图2，包含河西走廊及其邻区内的4条主要断裂带和河西走廊内的4个次级地块(盆地)。南部为祁连山地区，西北角和北部区域分别为塔里木地块和阿拉善地块的一部分；河西走廊断陷带贯穿模型的中部，其中包括4个次一级子块体：玉门盆地、酒泉盆地、张掖盆地和武威盆地。

自古近纪早期开始，受印度板块与欧亚板块于中生代末—新生代早期的碰撞及持续至今的向北推挤作用的远程效应的影响，青藏高原持续的隆升、阿尔金断裂带大规模的左行走滑和祁连山北缘断裂带的逆冲，形成河西走廊典型的新生代走滑—压陷盆地[29]。研究区域内活动断裂系包含4条主要断裂带[21]：祁连山北缘断裂带、走廊北缘断裂带、邻近河西走廊地区的阿尔金断裂带的东段和海原断裂带。由于构成上述断裂带的断层复杂多样，在模型中把属于同一条较大断裂带的断层束、断层段统一简化为宽度为5 km的断裂带表示，参数见表1。

ATF:阿尔金断裂带；NHCF:走廊北缘断裂带；NQLF:祁连山北缘断裂带；HYF:海原断裂带； YMB:玉门盆地；JQB:酒泉盆地；ZYB:张掖盆地；WWB:武威盆地。图2 有限元模型及边界条件Fig.2 Model domain and boundary conditions

断裂带名称走向倾向倾角/(°)深度/km海原断裂带弧形展布南西8020 (上地壳)阿尔金断裂带(东段)NEE近东西7560 (Moho)走廊北缘断裂带NWW南西7550 (Moho)祁连山北缘断裂带NWW西南7560 (Moho)

注：弧形展布表示的断层走向为向北东方向凸起的弧形。

本文将研究区域的经纬度坐标系投影转换到直角坐标系进行计算(X轴东向为正，Y轴北向为正，Z轴向上为正)。由于研究区域内地壳厚度较大(最厚为62 km)和莫霍面变化梯度高[30]，模型深度取地表以下100 km，已达到岩石圈地幔。据研究区域莫霍面等深线图[30]及青藏高原东北缘地壳三维速度结构探测结果[31-32]，垂直方向上模型分成4层：上地壳范围为地表至地下20 km深度处；中地壳深度为地下20～35 km；下地壳深度为地下48～62 km；岩石圈地幔深度为莫霍面至地下100 km的深度范围。计算模型采用六面体网格，并在断裂带地区在水平方向和垂直方向进行整体加密。网格尺度在水平方向上保持一致，长宽均为10 km，深度方向上将100 km厚的岩石圈模型分成11层，模型浅部网格厚度为5 km、深部网格厚度为10 km，整体共有48 941个节点、43 935个单元。

研究区域南部受到印度板块向北的推挤作用，而其北部受到阿拉善块体和鄂尔多斯块体的阻挡作用，据此对边界条件设置：1)水平方向上：据地壳运动观测网络给出研究地区的GPS速度场[28]作为一级近似插值到4个侧边界，作为水平速度约束条件，且假定从地表到100 km深度保持一致；2)垂直方向上：顶面为自由边界，即法向应力和剪应力均为零；鉴于目前对该区域岩石圈上地幔运动状态的未知，暂且将底面垂直方向速度约束为0，而水平方向约束自由。

1.2 动力学方程

青藏高原东北缘地壳缩短和河西走廊系列盆地演化过程中，上地壳将表现为脆性，而中、下地壳和岩石圈地幔介质的流变性质，据此，模型中黏弹性介质采用Maxwell体本构关系。地壳岩石圈的静态的力平衡方程为

(1)

其中，σij为应力张量(i,j=1, 2, 3)，ρgi为体力项。有限元模型的计算每一时间步的应变增量为

{dε}={dεv}+{dεe}+{dεp}.

(2)

式中：dεv、dεe和dεp分别表示黏性、弹性和塑性的应变增量。黏性与弹性应变的关系由Maxwell体的线性黏弹性关系给出，其本构关系描述为

{dεv}=[Q]-1{σt}dt=

[Q]-1({σt-dt}+{dσ})dt，

(3)

{dεe}=[D]-1{dσ}.

(4)

式中：{σt}为t时刻的应力张量，dt为时间增量，{dσ}为应力张量增量，[D]为弹性材料矩阵，[Q]为与黏度相关的材料矩阵。

(5)

(6)

式中：η为黏滞系数，E为杨氏模量，ν为泊松比。在模拟过程中，当加载达到材料屈服极限的时候，材料开始发生塑性变形，采用Drucker-Prager塑性屈服准则模拟上地壳和嵌于地壳中的断层：

(7)

(8)

(9)

(10)

式中：I1为应力张量的第一不变量，J2为偏应力张量的第二不变量，α和β为与C(内聚力)和φ(内摩擦角)相关的材料常数，σave为平均应力。这里挤压应力为负。

综上所述，三维黏弹塑性模型材料的本构方程可以表示为

(11)

1.3 模型参数选取

据前人研究结果及地质资料显示，青藏高原东北缘地区在向北东方向扩展的过程中，上地壳表现为脆性，而中、下地壳及地幔介质的流变性质在长时间尺度的构造过程中得以体现[33]，具体表现在黏弹性介质与松弛时间上的差别。青藏高原的中、下地壳低黏度带发育普遍，平均黏滞系数比周边稳定的塔里木地块和阿拉善地块的地壳黏滞系数低[33]。孙玉军等[34]通过岩石圈地震波速度结构、热结构和岩石圈岩石组成等数据计算得到鄂尔多斯地块及其周围三维岩石圈流变结构，结果显示青藏高原下地壳黏滞系数约为(0.4～3.0)×1021Pa·s；石耀霖和曹建玲[35]利用实验室流变试验结果、以岩石圈温度和应变速率的最新研究成果为基础，计算中国大陆岩石圈等效黏滞系数，其中高原下方中地壳为(3.0～5.0)×1020Pa·s，下地壳为(3.0～6.0)×1019Pa·s；万永革等[36]对青藏高原北部进行模拟计算时采用Maxwell体黏弹性介质模型，黏弹性系数在地壳上部19 km范围内取值1.0×1021Pa·s，在中地壳19～38 km取值2.0×1019Pa·s，在下地壳取值为6.3×1018Pa·s。综上，模型中地壳力学参数取值，见表2[34-38]。同时，参考地质资料并比较模型中河西走廊地区与其邻区(塔里木、阿拉善和祁连山地区)对应地层的力学参数数值的大小关系，将模型中4个次级地块区域设置为“硬包体”，即其力学性质(杨氏模量、黏滞系数)设置为对应地层力学性质(杨氏模量、黏滞系数)参数的2倍。

根据青藏高原东北缘各区域典型地震折射剖面P波和S波速度结构及泊松比结构[39]，对研究区域内岩石圈分层结构进行分析。弹性模量E和泊松比ν与地震波传播速度有如下关系：

(12)

本研究中采用ABAQUS软件的Newton-Raphson算法进行求解。将每个分析步划分为一系列的载荷增量步，每个增量步内又进行若干次迭代。据河西走廊中地壳的黏弹性松弛时间(黏滞系数与杨氏模量之比)约为203 a，下地壳的黏弹性松弛时间约为10 a，岩石圈地幔的黏弹性松弛时间约为87 a。在数值模拟过程中，以500 a为一个时间步长，计算10 000个时间步(约5 Ma)且应力已松弛并稳定。

2 数值模拟结果与分析

2.1 区域速度场分布

将研究区域内的水平速度场数值模拟计算结果与实测的GPS观测值进行对比(图3)，可以看出模拟计算出的水平速度场与GPS观测值整体上吻合，且在数值和方向上都与GPS数据有较好

表2 有限元模型的材料参数Table 2 Material parameters in the finite-element model

蓝色箭头表示GPS速度插值所得到的边界条件，红色箭头表示GPS速度，黑色箭头为模型计算得到的河西走廊及其邻区地表速度场。图3 研究区域内数值计算地表速度场与GPS数据对比图Fig.3 Comparison of the modeled velocity field with GPS velocity field

的一致性。同时，在区域背景应力场的作用下，研究区域地表变形速度场整体上呈现出南东向旋转；水平位移速率具有非常明显的分带性，即南部速度大于北部、东部速度大于西部。

并且，在此次数值模拟过程中，我们进一步计算了河西走廊及其邻区的4条主要断裂带的平均滑动速率，分段将其与对应断裂段的地质滑动速率和GPS滑动速率进行比较，见表3。数值计算得出的4条断裂带上的滑动速率与地质、GPS滑动速率吻合较好；海原断裂带祁连段的走滑速率达到4.5 mm/a，走廊北缘断裂带嘉峪关断裂的走滑速率达到2.8 mm/a，祁连山北缘断裂带佛洞庙—红崖子断裂的走滑速率达到3.5 mm/a。另外，计算得出合黎山南缘断裂、龙首山南缘断裂、玉门断裂和张掖断裂在背景应力作用下是有一定的滑动速率(低于1～2 mm/a)，但野外地质考察上未观测到，是否是断裂闭锁还是未滑动需要进一步探索。

表3 计算得出的4条主要断裂的平均滑动速率与地质滑动速率和GPS滑动速率比较结果Table 3 Comparison of the calculated average slip rates of the four major fault zones in the tectonic process of the Hexi Corridor with the geological and GPS slip rates in segments

2.2 区域内的垂直和水平位移

图4(a)为结合地质资料[14]得出数值模拟过程中设置的研究区域内的初始地表地形图。在模拟过程中，将北祁连山地区的初始高程设置为1 500 m，走廊北缘断裂带的初始高程设置为300～600 m，其他地区的初始高程皆与水平面相同。图4(b)给出在初始地形图的基础上，经过5 Ma构造变形之后研究区域的地形图。可以看出，在整体受压扭条件下，河西走廊北部受青藏高原块体挤压而发生变形，引起阿拉善块体边缘的地壳缩短；塔里木地区和阿拉善地区上地壳相对下沉，且阿拉善南缘地区上地壳向青藏高原推挤；河西走廊北缘相对隆升，且祁连山地区隆起同时促使河西走廊内部发生大规模变形褶皱；河西走廊内部，4个次级块体(硬包体)上地壳物质相对下沉，挤压到次级块体下方中地壳，使得4个次级块体之间形成山。

图4(c)给出地表垂直抬升速率数值计算云图。研究区域内，河西走廊及其邻区的构造运动表现为山体的相对大幅度抬升与盆地的相对强烈下沉，且与前人研究结果具有较好的吻合(见表4)。在构造背景应力压扭作用下，研究区域内的构造变形以整体隆升为主，但各地块隆升速度不均匀：1)研究区域内断裂带的上盘一侧抬升速度较快; 2)与其他地区相比，祁连山地区总体上抬升较快，且北祁连山的隆升速度比南祁连山快; 3)在整体压扭作用下，祁连山北缘断裂带榆木山断裂段呈现微弱东西引张，抬升速度比其两端慢，在计算得到的地表垂直抬升速率云图中呈现类似“缺口”的形态; 4)模型中河西走廊内部由于4个硬包体的存在而出现垂直运动差异，依次出现4个成左行排列且被NW向小的隆起分割的系列盆地。

图4 研究区域内垂直变形结果Fig.4 Contour map of vertical deformation in the study area

研究区域抬升速率/(mm/a)前人结果计算结果塔里木地块0.12～0.15[33]～0.24河西走廊走廊北缘断裂带～1.4[6,14,16]～1.8走廊内部0.36±0.11[14]～0.41走廊内系列盆地0.28±0.12[17-18]～0.29祁连山地区祁连山山体0.5±0.1[14]～0.62祁连山北缘断裂带0.65±0.15[6]～0.78海原断裂带0.41±0.11[24]～0.65阿拉善地块0.27±0.13[34, 37]～0.37

3 讨论

3.1 计算参数的影响分析

在数值模拟计算中计算参数的选取往往会对计算结果有较大的影响。前面数值模型据河西走廊及其邻区构造运动和河西走廊内系列盆地分布特征，将河西走廊内的次级块体设置“硬包体”，即其力学参数(杨氏模量、黏滞系数)比周围地层的强时，在区域压扭构造环境下，河西走廊内会出现一系列大小不等的4个盆地。为进一步分析河西走廊内4个次级块体的力学参数选取是否恰当及在区域压扭作用下地壳横向不均匀性对盆地形成演化的影响，通过改变河西走廊内4个次级块体的力学性质设置不同的对比模型。1)将河西走廊内的4个次级块体设置为“软包体”，将其力学参数(杨氏模量、黏滞系数)降低为周围地层对应参数的一半时的两种情况；2) 将河西走廊内的4个次级块体设置为和周围地层对应参数一样。图5给出这两种对比模型下计算得到的地形图，同前面硬包体的计算结果(图4(b))进行比较，可以看出：当河西走廊内的次级块体为软包体(图5(a))时，数值模拟计算结果同块体为“硬包体”相反，即4个次级地块不但没有生成盆地反而有所抬升，因此，当河西走廊内的次级块体为“硬包体”时的模拟结果更符合盆山相间的实际情况；而当河西走廊内4个次级块体性质与周围地层的相同时，河西走廊地区不会出现盆地和断裂的起伏(图5(b))，垂直方向上的位移变化仅取决于该区的初始形态。

图5 次级块体(包体)力学参数不同模型计算得到的地形图Fig.5 The vertical displacement nephogram calculated using three kinds of models

上述计算结果显示：在压扭作用下，当河西走廊内的次级块体为硬包体时，河西走廊内会依次形成4个大小不等的系列盆地，那么在实际地质和地球物理的观测中，这些“硬包体”是否可以被观测得到？目前，已有不少研究指出在青藏高原地区的一些大的盆地(如塔里木盆地)的层析成像结果为地震波高速和杨氏模量相应较高的地区[40]，且其Lg波和Q较大[41]及黏滞系数较高[35]。也就是说大的“硬包体”是存在的，那么河西走廊这些小的“硬包体”是否存在呢？早期的区域层析结果[42]分辨率较低，无法识别出河西走廊内的速度差异。但是，基于现今密集台网改进了的结果[43]得到的层析成像结果，已经可以看到河西走廊内部的波速差异，这为今后的研究提供了希望，进而能够结合新的地质和地球物理的其他观测资料，构建更反映实际的模型。

3.2 黑河流动模式

在一个垂直构造运动强烈发育的地区，地貌演化与构造运动过程密切相关，特别是在河流水系的演变发育过程中表现极为敏感，并且能够记录地体隆升与挤压扩展的重要信息。据此，将北祁连山地区河网的地貌特征与计算得到的该地区的地表垂直方向的抬升速率进行对比，作为讨论计算结果在河西走廊构造演化应用的一个重要方面。

在原有祁连山至河西走廊的缓坡地貌基础上，如果北祁连山有一个窄型条带地区发生迅速隆升，那么北祁连山既有的河流在向北流的过程中，遇到山体隆升，河道会出现两种演变可能：一是河流中段隆起，造成下游河流断头、上游河流倒流改道；二是河流切割山体的速率大于山体抬升速率，河流切割隆起的山脊形成深切河谷，继续沿原河道流动。结合地质资料，通过对北祁连山地区的河流水系(图6)进行分析，可以看到：除黑河干流发源于祁连山腹地外，其余河流皆发育于祁连山山前；北祁连山地区的河网大体趋势呈平行等距的弧状排列，总体呈放射状均匀分布。根据河流走向，将该地区的河网分为3类：第1类是像北大河，从北祁连山南坡发源，深深切割山脊，向北流入河西走廊，这说明可能这些河流在原来自南向北流的过程中，遇到北祁连山隆起时，侵蚀切割山体速率高于山体抬升速率，河流切割山体继续北流。第2类是像摆浪河、马营河等，从北祁连山北坡发源，自南向北流动，却在北祁连山山前断头，这说明这些河流可能是由于北祁连山山体抬升速度高于其切割速度，则在山脊断头或倒流。第3类河流则是像黑河，汇集了北祁连山北流的河水，平行于山间谷地流动，然后在北祁连山的榆木山段找到“缺口”而切割北祁连山北流。

F1:池家刺窝断裂; F2:阿尔金断裂; F3:嘉峪关断裂; F4:金塔南山断裂; F5:慕少梁断裂; F6:阿右旗断裂; F7:文殊山断裂; F8: 合黎山南缘断裂; F9:玉门断裂; F10:榆木山断裂; F11:龙首山断裂; F12:昌马断裂; F13:佛洞庙—红崖子断裂; ①石油河；②白杨河；③北大河；④洪水坝河；⑤摆浪河；⑥马营河；⑦黑河；⑧洪水河；⑨山丹河。图6 北祁连山地区河流水系构造简图Fig.6 Schematic map of the tectonics of the river system in the North Qilian Mountains

对比北祁连山地区河流水系构造简图(图6)和计算得到的研究区域内北祁连山地区的山体抬升速率云图(图3(b))，图中祁连山地区榆木山段的“缺口”可以解释该地区的的河流流动模式：黑河在北流过程中，在北祁连山地区隆升的速度较慢的榆树山段间隙，穿过祁连山，流入河西走廊，验证了模型的正确性。从而说明本研究中的模型计算可以正确反映北祁连山山脊迅速隆起的区域变形特征。

4 结论

在现有的研究基础上，以河西走廊及其邻区晚新生代所处的构造环境为框架，综合考虑研究区域内的地层结构、主要的地质构造单元、断裂带及其几何特征和地质构造单元的相对完整性等因素，建立河西走廊地区三维黏弹塑性有限元模型。利用实测GPS速度场作为边界条件进行约束，数值分析河西走廊及其邻区在区域压扭构造作用下，近5 Ma河西走廊内系列盆地的构造演化过程，得出如下结论。

1)青藏高原扩展过程中，祁连山地区作为青藏高原东北缘向北东方向扩展生长的前缘地区，其生长是由南向北有序向外扩展的过程和模式。同时，青藏高原东北缘向北东方向的扩展是通过一系列北西西向逆冲断裂带(如祁连山北缘断裂带、海原断裂带等)而实现的。根据地质资料[44-46]，这些断裂带在具有逆冲特征的同时，仍然发生地壳缩短和大量的左旋走滑，以至于整个地块具有向东滑移的趋势，这可能也可以作为青藏高原东北缘的向北东方向扩展生长的动力来源。

2)在压扭构造作用下的地壳横向不均匀性对河西走廊构造变形有较大的影响。通过设置模型中河西走廊内次级块体的力学性质(杨氏模量、黏滞系数)，建立3组对比模型，即将河西走廊内的4个次级块体设置为“硬包体”、“软包体”和同周围岩体参数一致的“包体”。在区域压扭作用下，只有在河西走廊内存在多个“硬包体”时，河西走廊及其邻区在5 Ma的构造演化过程中才会依次出现多个成左行排列的系列盆地，且盆地间被一组NNW-NW向的隆起分割所分割。

通过对比计算得到的研究区域内5 Ma演化过程的表面抬升速率图与祁连山地区的河网分布，解释了现今北祁连山地区河流网络分布的现象，也验证了我们模拟的正确性。而新的密集台站为我们的研究在层析成像技术方面提供了新的希望，希望在未来的研究中，能够结合新的地质和地球物理的其他的观测资料，将研究继续向前推进。