中国科学院大学学报  2019, Vol. 36 Issue (2): 196-207   PDF    
河西走廊系列盆地构造演化的三维数值模拟
李蔚琳, 程惠红, 张怀, 石耀霖     
中国科学院大学 中国科学院计算地球动力学重点实验室, 北京 100049
摘要: 河西走廊是青藏高原向东北方向扩展生长的前缘地区。在青藏高原向东北方向推挤过程中,走廊内形成盆山相间的地貌,自东向西为武威—张掖—酒泉—玉门盆地,呈一字形排列,是研究新生盆地构造演化的理想地区。基于河西走廊及其邻区晚新生代构造环境、现今盆地与断裂带的几何构造布展、GPS和历史地震资料,建立区域三维黏弹塑性有限元地质模型,以期动态刻画近5 Ma河西走廊系列盆地的构造演化过程,探讨地壳的横向不均匀性对该区域构造变形的影响。数值分析结果如下。1)在近5 Ma的区域压扭构造作用下,河西走廊多个力学性质较强的次级块体依次形成左行排列的系列盆地且被NNW-NW向的断裂带和隆起分割。2)与研究区域内的其他地区相比,河西走廊与祁连山和阿拉善地块的交界处及祁连山地区整体抬升速度较快,且北祁连山抬升速度大于南祁连山。3)祁连山北缘榆木山断裂段呈现微弱东西引张,抬升速度比两端慢,地表垂直抬升速率呈现“缺口”形态。4)与河西走廊相邻的塔里木和阿拉善地块的上地壳相对于不断隆升的青藏高原下沉,且挤压盆地中地壳,使得盆地之间出现隆起。数值分析的结果反映了河西走廊系列盆地的演化过程,解释了现今祁连山地区河流网络分布现象,也揭示了青藏高原东北缘向北东方向扩展生长过程中潜在的动力来源。
关键词: 河西走廊     青藏高原东北缘     构造演化     黏弹塑性     有限元模拟    
Three-dimensional numerical modeling of the tectonic evolution of the serial basins in the Hexi Corridor in Northwest China
LI Weilin, CHENG Huihong, ZHANG Huai, SHI Yaolin     
Key Laboratory of Computational Geodynamics of CAS, University of Chinese Academy of Sciences, Beijing 100049, China
Abstract: The Hexi Corridor, located in the northeastern margin of Qinghai-Tibet Plateau, is the leading edge of northeastern expanding of Qinghai-Tibet Plateau. In the northeastern thrusting of Qinghai-Tibet Plateau, there forms the landscape of basins and mountains in the Hexi Corridor, which is an ideal area for studying the tectonic evolution of the Cenozoic basins. The four basins, Wuwei, Zhangye, Jiuquan, and Yumen basins, lay in the Hexi Corridor from the east to the west. Based on the Late Cenozoic tectonic environment in the Hexi Corridor and its adjacent areas, the geometric structure of the present basins and fault zones, GPS, and historical seismic data, we establish a three-dimensional visco-elastic-plastic finite element model to describe the tectonic evolution process of the serial basins in the Hexi Corridor dynamically and explore the effects of lateral inhomogeneity of the crust on the tectonic deformation in this area. The results of the simulation are shown as follows. 1) The hard secondary blocks in the Hexi Corridor form four basins in a sequence of left echelon arrangement in the model and they are separated by the NNW-NW faults. 2) Compared with other regions, the Qilian Mountains and the junction areas of the Hexi Corridor with the Qilian Mountains and the Alashan block generally uplift fast and the uplifting speed of the Northern Qilian Mountains is higher than that of the Southern Qilian Mountains. 3) Under the overall compression-torsion, there appears an approximate "pull-point tectonic structure" gap on the Elm Shan fault in the northern Qilian fault zone where the uplifting speed is slower than in the adjacent part on the fault. 4) In the Tarim and Alashan blocks, the upper crust sinks correspondingly and squeezes lower blocks, and thus the tectonic framework between the basins and the mountains forms. The results of the modeling reflect the evolution of the serial basins in the Hexi Corridor, explain the distribution of the current river network in the Qilian Mountains area, and reveal the potential power source of the northeastern thrusting of the northeastern Tibetan Plateau.
Keywords: Hexi Corridor     Northeastern Tibetan Plateau     tectonic evolution     visco-elasto-plastic     finite element modeling    

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

Download:
蓝色曲线为北祁连山地区河流水系;①玉门盆地;②酒泉盆地;③张掖盆地;④武威盆地。 图 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个次一级子块体:玉门盆地、酒泉盆地、张掖盆地和武威盆地。

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

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

表 1 模型使用的断层几何参数 Table 1 Geometric parameters of faults in the model

本文将研究区域的经纬度坐标系投影转换到直角坐标系进行计算(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体本构关系。地壳岩石圈的静态的力平衡方程为

$ \frac{{\partial {\sigma _{ij}}}}{{\partial {x_j}}} + \rho {g_i} = 0 $ (1)

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

$ \left\{ {{\rm{d}}\varepsilon } \right\} = \left\{ {{\rm{d}}{\varepsilon ^{\rm{v}}}} \right\} + \left\{ {{\rm{d}}{\varepsilon ^{\rm{e}}}} \right\} + \left\{ {{\rm{d}}{\varepsilon ^{\rm{p}}}} \right\}. $ (2)

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

$ \begin{array}{*{20}{c}} {\left\{ {{\rm{d}}{\varepsilon ^{\rm{v}}}} \right\} = {{\left[ \mathit{\boldsymbol{Q}} \right]}^{ - 1}}\left\{ {{\sigma ^{\rm{t}}}} \right\}{\rm{d}}t = }\\ {{{\left[ \mathit{\boldsymbol{Q}} \right]}^{ - 1}}\left( {\left\{ {{\sigma ^{t - {\rm{d}}t}}} \right\} + \left\{ {{\rm{d}}\sigma } \right\}} \right){\rm{d}}t,} \end{array} $ (3)
$ \left\{ {{\rm{d}}{\varepsilon ^{e}}} \right\} = {\left[ \mathit{\boldsymbol{D}} \right]^{ - 1}}\left\{ {{\rm{d}}\sigma } \right\}. $ (4)

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

$ \begin{array}{*{20}{c}} {\left[ \mathit{\boldsymbol{D}} \right] = \frac{E}{{\left( {1 + \nu } \right)\left( {1 - 2\nu } \right)}}}\\ {\left( {\begin{array}{*{20}{c}} {1 - \nu }&\nu &\nu &0&0&0\\ \nu &{1 - \nu }&\nu &0&0&0\\ \nu &\nu &{1 - \nu }&0&0&0\\ 0&0&0&{0.5 - \nu }&0&0\\ 0&0&0&0&{0.5 - \nu }&0\\ 0&0&0&0&0&{0.5 - \nu } \end{array}} \right),} \end{array} $ (5)
$ {\left[ \mathit{\boldsymbol{Q}} \right]^{ - 1}} = \frac{1}{\eta }\left( {\begin{array}{*{20}{c}} {\frac{1}{3}}&{ - \frac{1}{6}}&{ - \frac{1}{6}}&0&0&0\\ { - \frac{1}{6}}&{\frac{1}{3}}&{ - \frac{1}{6}}&0&0&0\\ { - \frac{1}{6}}&{ - \frac{1}{6}}&{\frac{1}{3}}&0&0&0\\ 0&0&0&1&0&0\\ 0&0&0&0&1&0\\ 0&0&0&0&0&1 \end{array}} \right). $ (6)

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

$ F\left( {\sigma ,c} \right) = \alpha {I_1} + \sqrt {{J_2}} - \beta , $ (7)
$ \alpha = \frac{{2\sin \varphi }}{{\sqrt 3 \left( {3 + \sin \varphi } \right)}},\beta = \frac{{6{\rm{C}}\cos \varphi }}{{\sqrt 3 \left( {3 + \sin \varphi } \right)}}, $ (8)
$ {I_1} = 3{\sigma _{{\rm{ave}}}} = \frac{{{\sigma _{ii}}}}{3}, $ (9)
$ \begin{array}{l} \sqrt {{J_2}} = \\ \sqrt {\frac{1}{6}\left[ {{{\left( {{\sigma _1} - {\sigma _2}} \right)}^2} + {{\left( {{\sigma _2} - {\sigma _3}} \right)}^2} + {{\left( {{\sigma _3} - {\sigma _1}} \right)}^2}} \right]} . \end{array} $ (10)

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

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

$ \left\{ {{\rm{d}}\sigma } \right\} = \left( {\left[ {\mathit{\boldsymbol{\tilde D}}} \right] - \left[ {{\mathit{\boldsymbol{D}}_{\rm{p}}}} \right]} \right)\left\{ {{\rm{d}}\varepsilon } \right\} + \left\{ {{\rm{d}}\tilde \sigma } \right\} - \left\{ {{\rm{d}}{{\tilde \sigma }_{\rm{p}}}} \right\}. $ (11)
1.3 模型参数选取

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

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

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

$ E = \frac{{\rho V_{\rm{s}}^2\left( {3V_{\rm{p}}^2 - 4V_{\rm{s}}^2} \right)}}{{V_{\rm{p}}^2 - V_{\rm{s}}^2}},\nu = \frac{{V_{\rm{p}}^2 - 2V_{\rm{s}}^2}}{{2\left( {V_{\rm{p}}^2 - V_{\rm{s}}^2} \right)}}. $ (12)

式中:Vp是P波速度;Vs是S波速度;ρ是密度,且ρ=0.32Vp+0.77。计算中,当缺少VS的情况下,常常利用经验公式Vp/VS= $ \sqrt 3 $,由Vp求解VS。对于地幔介质,取统计关系式VS=Vp/1.86[31-32, 39]

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

2 数值模拟结果与分析 2.1 区域速度场分布

将研究区域内的水平速度场数值模拟计算结果与实测的GPS观测值进行对比(图 3),可以看出模拟计算出的水平速度场与GPS观测值整体上吻合,且在数值和方向上都与GPS数据有较好的一致性。同时,在区域背景应力场的作用下,研究区域地表变形速度场整体上呈现出南东向旋转;水平位移速率具有非常明显的分带性,即南部速度大于北部、东部速度大于西部。

Download:
蓝色箭头表示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个次级块体之间形成山。

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

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

表 4 计算的研究区域内主要地块的平均抬升速率与前人研究结果的比较 Table 4 Comparison of the calculated average slip rates of the major blocks in the tectonic process of the Hexi Corridor with the previous research results
3 讨论 3.1 计算参数的影响分析

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

Download:
图 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类河流则是像黑河,汇集了北祁连山北流的河水,平行于山间谷地流动,然后在北祁连山的榆木山段找到“缺口”而切割北祁连山北流。

Download:
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演化过程的表面抬升速率图与祁连山地区的河网分布,解释了现今北祁连山地区河流网络分布的现象,也验证了我们模拟的正确性。而新的密集台站为我们的研究在层析成像技术方面提供了新的希望,希望在未来的研究中,能够结合新的地质和地球物理的其他的观测资料,将研究继续向前推进。

参考文献
[1]
张培震, 邓起东, 张国民, 等. 中国大陆的强震活动与活动地块[J]. 中国科学:地球科学, 2003, 33(S1): 12-20.
[2]
李清河, 张元生, 涂毅敏, 等. 祁连山-河西走廊地壳速度结构及速度与电性的联合解释[J]. 地球物理学报, 1998, 41(2): 197-210. DOI:10.3321/j.issn:0001-5733.1998.02.008
[3]
Zheng W J, Zhang P Z, He W G, et al. Transformation of displacement between strike-slip and crustal shortening in the northern margin of the Tibetan Plateau:Evidence from decadal GPS measurements and late Quaternary slip rates on faults[J]. Tectonophysics, 2013, 584(1): 267-280.
[4]
Zheng D W, Clark M K, Zhang P Z, et al. Erosion, fault initiation and topographic growth of the North Qilian Shan (northern Tibetan Plateau)[J]. Geosphere, 2010, 6(6): 937-941. DOI:10.1130/GES00523.1
[5]
万景林, 郑文俊, 郑德文, 等. 祁连山北缘晚新生代构造活动的低温热年代学证据[J]. 地球化学, 2010, 39(5): 439-446.
[6]
张宁, 郑文俊, 刘兴旺, 等. 河西走廊西端黑山断裂运动学特征及其在构造转换中的意义[J]. 地球科学与环境学报, 2016, 38(2): 245-257. DOI:10.3969/j.issn.1672-6561.2016.02.012
[7]
袁道阳, 刘百篪, 吕太乙, 等. 北祁连山东段活动断裂带的分段性研究[J]. 地震工程学报, 1998(4): 27-34.
[8]
Tapponnier P, Xu Z, Roger F, et al. Oblique stepwise rise and growth of the Tibet Plateau[J]. Science, 2001, 294(5547): 1671-1677. DOI:10.1126/science.105978
[9]
Métivier F, Gaudemer Y, Tapponnier P, et al. Northeastward growth of the Tibet plateau deduced from balanced reconstruction of two depositional areas:The Qaidam and Hexi Corridor basins, China[J]. Tectonics, 1998, 17(6): 823-842. DOI:10.1029/98TC02764
[10]
Fielding E, Isacks B, Barazangi M, et al. How flat is Tibet?[J]. Geology, 1994, 22(2): 163-167.
[11]
王庆良, 王文萍, 崔笃信, 等. 青藏块体东北缘现今地壳运动[J]. 大地测量与地球动力学, 2002, 22(4): 12-16.
[12]
崔笃信, 王庆良, 王文萍, 等. 甘青块体东西部运动分界带的确定[J]. 大地测量与地球动力学, 2007, 27(1): 24-30.
[13]
陈杰, Wyrwoll K H, 卢演俦, 等. 祁连山北缘玉门砾岩的磁性地层年代与褶皱过程[J]. 第四纪研究, 2006, 26(1): 20-31. DOI:10.3321/j.issn:1001-7410.2006.01.004
[14]
郑文俊.河西走廊及其邻区活动构造图像及构造变形模式[D].北京: 中国地震局地质研究所, 2009. http://www.cnki.com.cn/Article/CJFDTotal-GJZT201003011.htm
[15]
陈柏林, 刘建生, 张永双, 等. 玉门断裂全新世活动特征及其与玉门地震的关系[J]. 地质论评, 2005, 51(2): 138-142. DOI:10.3321/j.issn:0371-5736.2005.02.004
[16]
陈柏林, 王春宇, 崔玲玲. 河西走廊西段南北向左行逆冲活动断裂的发现及其意义[J]. 地质学报, 2009, 83(7): 937-945. DOI:10.3321/j.issn:0001-5717.2009.07.004
[17]
陆洁民, 郭召杰, 赵泽辉, 等. 新生代酒西盆地沉积特征及其与祁连山隆升关系的研究[J]. 高校地质学报, 2004, 10(1): 50-61. DOI:10.3969/j.issn.1006-7493.2004.01.004
[18]
陈汉林, 杨树锋, 肖安成, 等. 酒泉盆地南缘新生代冲断带的变形特征和变形时间[J]. 石油与天然气地质, 2006, 27(4): 488-494. DOI:10.3321/j.issn:0253-9985.2006.04.008
[19]
杨树锋, 陈汉林, 程晓敢, 等. 祁连山北缘冲断带的特征与空间变化规律[J]. 地学前缘, 2007, 14(5): 213-223.
[20]
国家地震局地质研究所, 宁夏回族自治区地震局. 海原断裂带[M]. 北京: 地震出版社, 1990: 6-255.
[21]
国家地震局地质研究所, 国家地震局兰州地震研究所. 祁连山-河西走廊活动断裂系[M]. 北京: 地震出版社, 1993: 1-340.
[22]
宋春晖, 赵志军. 青藏高原北缘酒西盆地13 Ma以来沉积演化与构造隆升[J]. 中国科学:地球科学, 2001, 31(S1): 155-162.
[23]
冉勇康, 邓起东. 海原断裂的古地震及特征地震破裂的分级性讨论[J]. 第四纪研究, 1998, 18(3): 271-278. DOI:10.3321/j.issn:1001-7410.1998.03.011
[24]
邓起东, 张维岐, 张培震, 等. 海原走滑断裂带及其尾端挤压构造[J]. 地震地质, 1989, 11(1): 1-14.
[25]
崔笃信, 王庆良, 胡亚轩, 等. 青藏高原东北缘岩石圈变形及其机理[J]. 地球物理学报, 2009, 52(6): 1490-1499. DOI:10.3969/j.issn.0001-5733.2009.06.010
[26]
Lasserre C, Peltzer G, Peltzer G, et al. Interseismic strain across the Altyn Tagh Fault System (Northern Tibet), measured by SAR interferometry[C]//Agu Fall Meeting Abstracts, 2001: 1977-1985.
[27]
Lasserre C, Gaudemer Y, Tapponnier P, et al. Fast late Pleistocene slip rate on the Leng Long Ling segment of the Haiyuan fault, Qinghai, China[J]. Journal of Geophysical Research Solid Earth, 2002, 107(B11): ETG-1-ETG 4-15.
[28]
Liang S, Gan W, Shen C, et al. Three-dimensional velocity field of present-day crustal motion of the Tibetan Plateau derived from GPS measurements[J]. Journal of Geophysical Research Solid Earth, 2013, 118(10): 5722-5732. DOI:10.1002/2013JB010503
[29]
陈柏林, 刘建生. 祁连山北缘-河西走廊地区大地形变与地震的关系[J]. 地质通报, 2009, 28(10): 1439-1447. DOI:10.3969/j.issn.1671-2552.2009.10.010
[30]
刘峡, 黄立人, 杨国华, 等. 北祁连-河西走廊地区垂直形变与构造应力场关系的初步研究:三维有限元拟合结果分析[J]. 地震地质, 2003, 25(2): 307-316. DOI:10.3969/j.issn.0253-4967.2003.02.016
[31]
Liu M, Mooney W D, Li S, et al. Crustal structure of the northeastern margin of the Tibetan plateau from the Songpan-Ganzi terrane to the Ordos basin[J]. Tectonophysics, 2006, 420(1/2): 253-266.
[32]
刘明军, 李松林, 方盛明, 等. 利用地震波速研究青藏高原东北缘地壳组成及其动力学[J]. 地球物理学报, 2008, 51(2): 412-430. DOI:10.3321/j.issn:0001-5733.2008.02.014
[33]
曹建玲, 石耀霖, 张怀, 等. 青藏高原GPS位移绕喜马拉雅东构造结顺时针旋转成因的数值模拟[J]. 科学通报, 2009, 54(2): 224-234.
[34]
孙玉军, 董树文, 范桃园, 等. 中国大陆及邻区岩石圈三维流变结构[J]. 地球物理学报, 2013, 56(9): 2936-2946.
[35]
石耀霖, 曹建玲. 中国大陆岩石圈等效黏滞系数的计算和讨论[J]. 地学前缘, 2008, 15(3): 82-95. DOI:10.3321/j.issn:1005-2321.2008.03.006
[36]
万永革, 沈正康, 曾跃华, 等. 青藏高原东北部的库仑应力积累演化对大地震发生的影响[J]. 地震学报, 2007, 29(2): 115-129. DOI:10.3321/j.issn:0253-3782.2007.02.001
[37]
柳畅, 石耀霖, 郑亮, 等. 三维黏弹性数值模拟华北盆地地震空间分布与构造应力积累关系[J]. 地球物理学报, 2012, 55(12): 3942-3957. DOI:10.6038/j.issn.0001-5733.2012.12.007
[38]
朱广生, 桂志先, 熊新斌, 等. 密度与纵横波速度关系[J]. 地球物理学报, 1995, 38(S1): 260-264.
[39]
杨辉, 王勇. 青藏高原低速层对岩石圈强度影响的三维有限元模拟[J]. 大地测量与地球动力学, 2007, 27(2): 25-31.
[40]
Wei W, Zhao D, Xu J, et al. Depth variations of P-wave azimuthal anisotropy beneath Mainland China[J]. Sci Rep, 2016, 6: 29614. DOI:10.1038/srep29614
[41]
Zhao L F, Xie X B, He J K, et al. Crustal flow pattern beneath the Tibetan Plateau constrained by regional Lg-wave Q tomography[J]. Earth & Planetary Science Letters, 2013, 383(4): 113-122.
[42]
潘佳铁, 李永华, 吴庆举, 等. 青藏高原东南部地区瑞雷波相速度层析成像[J]. 地球物理学报, 2015, 58(11): 3993-4006.
[43]
潘佳铁, 李永华, 吴庆举, 等. 基于密集流动地震台阵的青藏高原东北缘及邻区Rayleigh波相速度层析成像[J]. 地球物理学报, 2017, 60(6): 2291-2303.
[44]
Burchfiel B C, Zhang P, Wang Y, et al. Geology of the Haiyuan Fault Zone, Ningxia-Hui Autonomous Region, China, and its relation to the evolution of the Northeastern Margin of the Tibetan Plateau[J]. Tectonics, 1991, 10(6): 1091-1110. DOI:10.1029/90TC02685
[45]
Meyer B, Tapponnier P, Gaudemer Y, et al. Rate of left-lateral movement along the easternmost segment of the Altyn Tagh fault, east of 96°E (China)[J]. Geophysical Journal of the Royal Astronomical Society, 1996, 124(1): 29-44. DOI:10.1111/gji.1996.124.issue-1
[46]
袁道阳, 张培震, 刘百篪, 等. 青藏高原东北缘晚第四纪活动构造的几何图像与构造转换[J]. 地质学报, 2004, 78(2): 270-278. DOI:10.3321/j.issn:0001-5717.2004.02.017