外源流水对流域地貌的重塑：河流下切与流域袭夺*

李伊菲，石耀霖，张怀

(中国科学院大学地球与行星科学学院 中国科学院计算地球动力学重点实验室， 北京 100049) (2021年2月8日收稿； 2021年3月23日收修改稿)

高海拔地区，由于高原内部通常位于地形屏障的背风侧，因此具有干旱到半干旱气候的特征，流域降水较少，盆地演化缓慢，有利于内流盆地的形成[1]。除此之外，在低海拔地区，也分布着内流盆地。在新近纪晚期和第四纪早期，许多内流盆地已经转换成了外流盆地[2]，如河套盆地、共和盆地、Mianeh盆地、Ebro盆地[1,3-5]。内流盆地转换成外流盆地后，给下游带来了大量的外源流水。内流湖泊外泄在地质上属于瞬时事件，在相对比较稳定的构造演化背景下，会对流域或者盆地的构造演化产生深远影响。这种地质过程现今仍在进行，例如，青藏高原可可西里卓乃湖在2011年9月溃决后,引发了学术界关于其流域下游一系列盐湖发生顺序湖水外泄，逐步贯通，并发展成为长江支流的可能性的讨论[6-7]。

湖水的外流机制分为两种：可能是内流盆地湖水溢出导致的，也可能是相邻的河流溯源侵蚀连通内流盆地的结果[8]。伊朗高原的Mianeh盆地在～4—7 Ma是内流盆地，～4—5 Ma存在一个1 km高0.5 km深的Mianeh古湖，高原边缘的高水位和平缓地形的共同作用导致湖泊溢出，进而导致Qezel Owzan河在～4 Ma的时间开始加速下切，至今切割出一个约1 km深的基岩峡谷[1]。河套盆地在上新世至晚更新世早期是内流盆地，由于下游的原晋陕黄河的溯源侵蚀，切穿了晋陕峡谷北部的分水岭，沟通了河套盆地。河套盆地与晋陕黄河贯通，快速下切，形成现在的晋陕峡谷[9]。

无论是哪种机制，湖水的外泄都会给下游的河流和流域带来大量的外源流水，引起流域地貌的重塑。外源流水的注入如何影响流域地貌的形态和特征？这些特征又受哪些构造、地表过程的影响？为了回答这些问题，本文建立二维地表演化数值模型，对外源流水注入过程对河流地貌的影响进行系统的数值模拟研究。

1 模型与数值方法

1.1 控制方程

外源流水对河流地貌的重塑是外源流水注入过程和河流地貌演化过程共同作用的结果，因此需要模拟外源流水注入过程和地表演化过程两部分。其中地表演化过程由河流下切过程和山坡扩散过程组成，在构造活动地区高程变化还受到地壳垂直运动的影响。我们在不规则的Voronoi多边形网格上建立了一系列二维有限体积计算格式，对控制方程进行相应离散化并得到数值解[10]。

模型中河流下切过程对地表的侵蚀作用，假设河流下切速率是河道剪应力的幂函数[11]：

(1)

式中：z[L]为基岩高程，t[T]为时间，ke为剪应力系数，τb[M·L-1·T-2]为河道剪应力，τc[M·L-1·T-2]为河流下切剪应力阈值，a是无量纲的指数系数。

在长时间尺度的计算中，假设水流是准稳态的，质量和动量守恒，结合曼宁定律(Manning equation)，可以得到河道剪应力

(2)

其中:ρ是水密度，g是重力加速度，Q[L3·T-1]是水流量，w[L]是河道宽度，S是河道梯度，N是曼宁系数，α和β是无量纲的指数系数。

由于水流量Q的计算代价太大，因此引入水力几何关系，将计算水流量替换为计算汇流面积[12]

Q=kqAc,

(3)

(4)

其中:kw和kq是带量纲的系数，b和c是无量纲的指数系数，A[L2]是汇流面积。令τc=0，联立上述方程可以得到

(5)

山坡上蠕变、雨水飞溅和生物输运都会引起扩散过程，对凸面山顶的发育起重要的控制作用[13]。这里，考虑山坡扩散过程产生物质搬运在长时间尺度下的平均作用，采用线性扩散过程来进行模拟[14-15]。

(6)

式中Kd为山坡扩散系数。所有的算例中都取Kd=1×10-2m2·a-1。

将隆升过程简化为在时空上是不变的，因此给定隆升速率即可对全区域进行计算。

(7)

式中，U[L·T-1]为流域隆升速率。

1.2 外源流水注入和径流计算

模型中的径流包括外源流水注入产生的径流和研究区域降水产生的径流两部分。我们假设在研究区域的二维空间上，降水在空间和时间上是均一的。由于自然界中地下水的流动具有随机性[16]，不考虑地下水和蒸发对径流的影响，降水全部转化成径流。径流选择单向流模型，即水流和水流携带的物质流向高程小于本单元且梯度最大的相邻单元。研究区域内降水产生的节点的汇流面积等于上游所有节点的面积之和。采用CASCADE算法，对径流流向和节点的汇流面积进行计算[17]。

Pelletier[18]在讨论Grand Canyon晚新生代演化的过程中，假设在湖水外泄前，区域内没有下切过程。湖水外泄后，开始下切，来模拟湖水外泄对流域的影响。在我们的模型中，不讨论某个具体的流域，也不讨论外源流水的来源。而是讨论不同的外源流水注入流量和外源流水注入位置对流域的重塑作用。因此，通过给指定单元赋外源流水注入流量的方式来模拟外源流水注入的过程。

假设外源流水的注入口只有一个，研究外源流水对下游流域及相邻流域的重塑作用。考虑外源流水注入后对河流地貌在长时间尺度下的影响，假设外源流水注入的过程中的水流仍是准稳态的过程。不考虑外源流水注入后，注入产生的径流对下游流域的瞬时作用。根据汇流面积A和水流量Q的幂率关系，即式(3)，给定外源流水注入流量对应的汇流面积，在外源流水注入开始后，增加注入位置对应单元的汇流面积后，根据上述汇流面积的计算方法，计算全区域的汇流面积A。

1.3 模型与初边值条件

图1展示了模型的初始和边界条件，其中黑色和洋红色线条是初始河网(汇流面积A>1×108m2)，黑色线条是主河道(汇流面积A>1×1010m2)。边界条件用黑色直线表示，其中虚线表示开放边界，实线表示封闭边界。红色圆点表示供选择的外源流水注入的位置。初始地形是保持边界条件不变，给定带有低缓坡度和随机扰动的初始地形(0≤z≤5)，对该区域施加均匀的隆升速率、河流下切系数和山坡扩散系数的情况下，演化1 Ma后得到的2条准垂直于底边界的山脊山谷。需要注意的是，在后文中讨论参数的影响时，每个算例在外源流水注入前后的河流下切系数和隆升速率是不变的。因此，每个算例的初始地形高度可能不同，但是流域的河网分布是相似的。

图1 模型的初始与边界条件Fig.1 Initial and boundary conditions

2 模拟结果

我们的模拟，旨在通过对比不同参数演化得到的外源流水注入后河流地貌的演化过程和特征，归纳这种情况下流域袭夺、河流深切的机制和控制因素。算例所用参数参考前人研究成果[13,19-21]，见表1。

表1 数值模型采用的参数Table 1 Parameters used in the numerical experiments

2.1 流域袭夺与河流深切

图2展示了河流下切系数Kc=1×10-4a-1，扩散系数Kd=1×10-2m2·a-1，外源流水注入流量对应的汇流面积A=3×1011m2，隆升速率U=1×10-3m·a-1，外源流水注入位置为图中红点所示的情况下，外源流水注入后下游流域(左侧流域)和相邻流域(右侧流域)河流地貌的地形、河流下切速率和分水岭处河网的演化过程。

图2 外源流水注入后河流地貌的演化过程Fig.2 Fluvial landscape evolution after the exogenous inflow

结果显示，外源流水注入后初期，左侧连接外源流水注入点和主河道的支流、左侧主河道的河流下切速率快速增加(图2(d))，与外源流水注入点相连的支流迅速转换为主河道，左侧主河道处高程迅速下降，河道被深切(图2(a))。图2(g)展示了部分初始分水岭位置(图2(d)中黑色方框)的河流下切速率的变化量和河网。其中，背景颜色表示河流下切速率和初始河流下切速率的差值：蓝色表示河流下切速率加快，红色表示河流下切速率减慢。蓝绿色表示左侧流域的河网，洋红色表示右侧流域的河网，不同颜色河网之间是分水岭所在位置。演化100 a时，外源流水注入的影响还没有传播到分水岭处(图2(g))。

随着左侧主河道不断快速下切，主河道高程下降(图2(b))，即支流的侵蚀基准点下降。支流与主河道之间的梯度增加，导致山坡处支流的河流下切作用加强，河流下切增强的位置逐渐从主河道位置传播到山坡的位置(图2(e))。初始分水岭左侧的山坡加速下切，分水岭右侧流域还没有被影响(图2(h))。左侧主河道经过调整，高程和梯度降低，河流下切的速率逐渐恢复至外源流水注入前的大小(图2(e))，主河道高程逐渐趋于平稳。

分水岭左侧的山坡受加速的支流下切的影响，高程不断下降(图2(c))，导致分水岭右侧的单元的径流流向从向右流入右侧主河道转为向左，归为左侧流域。左侧流域逐渐对右侧流域进行袭夺，分水岭向右移动(图2(f)和2(i))。同时，右侧流域的部分单元逐渐被袭夺到左侧流域，促使左侧流域山坡处单元的支流流量增加，支流加速下切；右侧流域山坡处单元的支流流量减小，河流下切减缓(图2(i))，单元高程增高，成为新的分水岭。新的分水岭形成后，随着左侧流域的不断扩张，右侧流域不断被袭夺，周而复始，分水岭不断向右移动。左侧流域逐渐增大，右侧流域逐渐减小(图2(h)和2(i))。

2.2 外源流水注入影响下的河流地貌特征

2.2.1 外源流水注入流量对河流地貌特征的影响

通过改变外源流水注入流量对应的汇流面积，讨论外源流水注入流量对河流地貌特征的影响。图3(a)～3(c)展示了其他参数与算例1相同，外源流水注入流量对应汇流面积增加到A=3×1013m2后，河流地貌的演化结果。外源流水注入流量增加后，左侧主河道的河流下切速率增加得更明显(图3(a))，主河道被强烈地深切。左侧流域袭夺右侧流域的速度更快，分水岭处的河流下切速率的变化量更大(图3(b))。演化1 Ma后左侧流域袭夺了更多的右侧流域(图3(c))。

图3(d)～3(f)展示了其他参数与算例1相同，当外源流水注入流量对应汇流面积减小到A=3×109m2后，河流地貌的演化结果。外源流水注入给主河道增加的径流量较小，引起的河流下切作用的增幅也较小(图3(d))。因此，主河道和两侧山坡间的梯度的增加量较小。与算例1相比，河流下切增强的位置从主河道位置传播到山坡的速度较慢，105a时还未传播到分水岭的位置(图3(e))。演化1 Ma后分水岭基本保持在原来的位置，左侧流域基本还未开始袭夺右侧流域(图3(f))。

2.2.2 外源流水注入位置对河流地貌特征的影响

通过改变汇流面积的输入单元，讨论外源流水注入位置对河流地貌特征的影响。首先讨论外泄位置垂直于主河道流动方向移动的情况。图4(a)～4(c)展示了其他参数与算例1相同，外源流水注入位置从左侧主河道正上方(400,150)转换到流域左上角(400,200)后河流地貌的演化特征。结果显示，外源流水注入位置决定了外源流水注入位置附近的河道分布情况，原先的主河道慢慢演化成支流，外源流水注入带来的大量径流使原先的支流被快速下切成为主河道(图4(c))。外泄位置的改变，对外泄位置下游的流域袭夺和河流深切过程影响不大(图4(a)～4(c))。

接着，讨论外泄位置平行于主河道流动方向的情况。图4(d)～4(i)展示了其他参数与算例4(本节中上一算例)相同，外源流水注入位置由左侧流域左上角(400,200)转换到流域中部左侧(200,200)后，河流地貌的演化特征。外泄位置下游部分(x<200 km)，左侧流域袭夺右侧流域的过程和算例4基本相同(图4(d),4(f))。对外源流水注入位置的上游(x>200 km)，主河道没有直接接收外源流水注入带来的径流，因此没有立刻开始加剧下切。下游主河道被深切后，上游主河道因为和下游主河道间梯度的增加，逐渐开始加剧下切，并向两侧山坡传播(图4(d))。由于上游主河道的高程是逐渐降低的，因此和算例1 相比，左侧上游流域袭夺右侧上游流域的开始时间更靠后，袭夺的速度也更慢(图4(e))。演化1 Ma后，左侧流域袭夺右侧流域的面积较算例4较小，上游河流的深切量较小(图4(f))。

图3 不同外源流水注入流量下的河流地貌特征Fig.3 Fluvial landscape features with different exogenous inflow discharge

图4 不同外源流水注入位置下的河流地貌特征Fig.4 Fluvial landscape features with different exogenous inflow position

2.3 地表过程对河流地貌特征的影响

2.3.1 河流下切系数对河流地貌特征的影响

采用不同的河流下切系数，讨论河流下切系数对河流地貌特征的影响。结果显示，河流下切系数减小到Kc=1×10-5a-1后，外源流水注入后，河谷被深切，河流下切加强区域从主河道向两侧山坡传播的速度减慢(图5(a))。算例1中，3×104a时，河流下切增强的区域已经传播到分水岭处(图2(h))。但在本例中，105a时，河流下切的增强区域仍未传播到分水岭处(图5(a))，分水岭位置保持不变(图5(b))。1 Ma时，左侧流域对右侧流域的袭夺不明显，分水岭位置基本保持不变(图5(c))。

河流下切系数增大到Kc=3×10-4a-1后，河流下切加强区域从主河道向两侧山坡传播的速度增快(图5(d))。算例1中，105a时，河流下切的增强区域传播刚刚越过初始分水岭(图2(i))。但在本例中，105a时，河流下切的增强区域早已越过初始分水岭(图5(d)和5(e))。左侧流域袭夺右侧流域的面积更大，分水岭位置早已右移(图5(e))。1 Ma时，左侧流域对右侧流域的袭夺很明显(图5(f))。

2.3.2 隆升速率的影响

此例采用不同的隆升速率，图6展示了隆升速率U=3×10-3m·a-1和U=3×10-4m·a-1的情况下，河流地貌的演化过程和特征。结果显示，当隆升速率增加到U=3×10-3m·a-1后，河流下切的加强区域从主河道向两侧山坡传播的速度与算例1基本相同，初始分水岭附近的河流下切的增强和减弱区域基本一致，河流下切的改变量明显增加(图6(a)和6(b))。演化1 Ma后，左侧流域对右侧流域的袭夺基本一致，区域高程明显增加(图6(c))。

当隆升速率减小到U=3×10-4m·a-1后，河流下切的加强区域从主河道向两侧山坡传播的速度与算例1基本相同，初始分水岭附近的河流下切的增强和减弱区域基本一致，河流下切的改变量明显减小(图6(d)和6(e))。演化1 Ma后，左侧流域对右侧流域的袭夺基本一致，区域高程明显减小(图6(f))。

图5 不同河道下切系数下的河流地貌特征Fig.5 Fluvial landscape features with different incision coefficient

图6 不同隆升速率下的河流地貌特征Fig.6 Fluvial landscape features with the different uplift rate

隆升速率与演化得到的地形高度、局部梯度呈正相关。河流下切过程、山坡扩散过程和局部梯度呈正相关。因此，隆升速率变化时，地表演化过程、流域袭夺过程一起变化。因此，隆升速率影响流域的平均高程和局部梯度，不影响流域袭夺的进程。

3 讨论与启示

3.1 流域袭夺

Whipple等[22]在探讨分水岭迁移时提出了正反馈和负反馈两套体系。侵略方流域的高侵蚀率驱动了分水岭的移动，分水岭移动后，流域面积加大，进而促进了侵略方流域侵蚀率的增加。这种正反馈会促使流域分水岭不断移动。负反馈则是指，高侵蚀率驱动分水岭移动后，河道纵剖面的不均衡增加，产生了局部的隆升和侵蚀间的不平衡，河道剖面调整，使得这种纵剖面的不均衡减小，分水岭的移动被抑制。我们在这种框架下，通过分析流域袭夺面积随时间的变化，探讨不同因素在外源流水注入引起的流域袭夺过程中，所起的作用和动力学原因。

图7(a)展示了不同外源流水注入流量下流域袭夺面积随时间的变化，外源流水注入流量决定了左侧流域是否对右侧流域袭夺以及袭夺的速度。外源流水注入是直接作用在左侧流域主河道上的，直接影响左侧主河道的河流下切速率并呈正相关。因此，外源流水注入流量增大，主河道高程降低速率增大。山坡上支流的侵蚀基准点高程下降得更快，引起的支流纵剖面的不均衡更大，支流更快地进行调整，因此流域的袭夺也就更快。同理，当外源流水注入流量较小时，引起的支流纵剖面的不均衡，不足以驱动对相邻流域的袭夺，则可能导致外源流水注入后对相邻流域不产生影响。由于外源流水注入流量不直接作用在山坡和支流上，因此河道剖面的不均衡传递到分水岭的时间差异不大。这导致不同外源流水注入流量引起的左侧流域开始袭夺的时间差异不大(图7(a))。

图7(b)展示了不同外源流水注入位置下流域袭夺面积随时间的变化，图中的外泄位置对应的是图1(b)所示节点号。当外泄位置在准垂直于主河道的方向上改变时，外源流水注入位置不影响其下游流域是否袭夺相邻流域、袭夺的速度和开始袭夺的时间。但是，当外泄位置沿主河道方向改变时，外泄位置上游流域的袭夺将会减缓，也造成整体流域的平均袭夺速度减缓。

图7 不同参数下的袭夺流域面积Fig.7 Captured drainage area with different parameters

图7(c)展示了不同河流下切系数下流域袭夺面积随时间的变化，河流下切系数与流域袭夺的速度呈正相关，与流域开始袭夺的时间呈负相关。河流下切系数控制整体流域的河流下切过程。河流下切系数增大，主河道高程降低速率增大，山坡上支流的侵蚀基准点高程下降得更快，引起的支流纵剖面的不均衡增大。由于支流也被更大的河流下切系数控制，因此更快地进行调整以适应纵剖面的不均衡，将其传递到初始分水岭，导致流域袭夺的时间更早。高河流下切系数的算例中，流域袭夺在3×104a左右基本停止，而非无穷无尽地进行下去(图7(c))，这是因为隆升与侵蚀之间的负反馈机制抑制了流域袭夺的进一步发展。

图7(d)展示了不同隆升速率下流域袭夺面积随时间的变化。隆升速率增加，同时增加地形梯度和高程。河流下切过程、山坡扩散过程都与地形梯度成正比。因此，隆升速率的增加主要影响流域地形高程和梯度。但与流域袭夺的开始时间和流域袭夺的速度无关(图7(d))。

3.2 低起伏面的形成

Yang等[23]指出，构造变形引起的河流重组可能会导致流域缺失。缺失流域面积的流域河流下切减弱，无法平衡构造隆升，进而会导致低起伏面的形成。在我们的模拟中，外源流水注入后可能会引发下游流域袭夺相邻流域，造成相邻流域的流域面积缺失。因此，我们通过提取模型中的纵剖面，讨论外源流水注入会如何影响流域的地形和起伏，以及是否会引起相邻流域向低起伏面发展。

图8(a)和8(b)展示了算例10的初始地形河网和演化1 Ma后的地形河网，外源流水注入位置为图8(a)和8(b)中点B。在本例中，选择较高的河流下切系数和外源流水注入流量，因此左侧流域对右侧流域的袭夺面积较大，右侧主河道明显变短。分别选择左侧流域支流最高点A、外源流水注入点B、右侧流域被袭夺过的支流最高点C、右侧流域未被袭夺过的支流最高点D(图8(a)～8(b))。提取节点A～D到主河道的支流和对应主河道的纵剖面(图8(c)～8(f))。横轴表示节点沿河道到出水口的距离，纵轴表示高程，黑点表示支流和主河道的交汇点，即支流的侵蚀基准点。

图8(d)代表左侧流域接收外源流水注入的支流和主河道。主河道快速地被深切，达到新的平衡。图8(c)代表左侧流域中，除了外源流水注入点连接主河道支流外，其他支流及相连左侧主河道的演化结果。左侧主河道的快速下切导致支流的侵蚀基准点快速下降，促使支流快速下切。支流的纵剖面上有明显的因为侵蚀基准点下降导致的裂点，裂点(图8(c)中绿色点)逐渐沿支流向上游移动，最终达到新的平衡。

图8 支流和主河道纵剖面Fig.8 Longitudinal profiles of tributaries and main channels

右侧流域中未被袭夺的山坡单元(图8(f))，主河道在右侧流域被袭夺前，没有变化(演化100 a的纵剖面与初始纵剖面重合)。随着流域不断被袭夺，右侧流域面积缺失，主河道高程逐渐上升，纵剖面的局部梯度逐渐减小(图8(f))。由此可以推断，外源流水注入导致的流域袭夺，也可能会导致被袭夺的流域形成低起伏的地形。图8(e)展示了右侧流域中，被左侧流域袭夺的支流和主河道的演化过程。右侧流域被袭夺前，河道纵剖面的演化过程与图8(f)相似。被袭夺后，新的支流连通左侧高程更低的主河道，支流和主河道间梯度增大，逐渐加速支流的下切(图8(e))。

3.3 地质启示

现今的汾河是黄河的一条支流，发源于宁武县管涔山脉，流经太原和临汾盆地，于河津汇入黄河，干流全长694 km[24]，与汾河毗邻的晋陕黄河北起托克托县河口镇，南至河津，全长700多km[25]。张珂[9]指出，新近纪晚期曾有从滹沱盆地南流至晋中盆地的古河，晋陕黄河和汾河曾经长度相近，是吕梁山隆升产生的长度相近的两条相似的河流，汾河的规模可能更大，最终汇入侯马、运城及渭河等盆地。托克托台地位于晋陕峡谷上峡口，蒋复初等[3]认为，河套古湖的湖水于100 kaB. P.开始外泄，河套—晋陕峡谷贯通。湖水迅速下泄到晋陕峡谷中，晋陕河谷迅速下切[9]。河套古湖的湖水下泄，可能促使了晋陕黄河流域对汾河流域的部分袭夺，使原先发源于滹沱盆地的汾河上游被部分袭夺，汾河长度减小，其演化过程可能与图8相似。汾河的纵剖面存在明显的裂点[26]，除构造因素的影响外，可能也与晋陕黄河接收河套古湖湖水后加速下切，汾河的侵蚀基准点下降有关(与图8(c)中裂点迁移相似)。因此，河套古湖的湖水外泄可能在晋陕黄河和汾河流域的演化过程中，起到了不可忽视的作用。新近纪晚期和现今2条河流的差异，可能与河套古湖湖水下泄有关。

本文主要对山地区域外源流水对流域地貌产生的影响进行了讨论。对平原区域的河流，由于河道高差降低，河道梯度减小，径流的增加对河道的下切加速作用没有山地区域强。另外，由于径流速度减慢，因此上游带来的泥沙无法被径流快速带走，而是可能沿程淤积，甚至引发溃堤。例如历史上黄河曾经多次夺淮入海：黄河中下游多沙少水的特征导致其夺淮后，不仅给淮河带去了大量的径流，还给淮河带去了大量的泥沙。而淮河原本的水力条件无法输运如此大量的径流和泥沙，导致泥沙沉积在河道中，甚至堵塞淮河的入海口，进一步导致泥沙的溯源淤积[27-29]。因此对于平原区域的外源流水问题，需要进一步的研究和讨论。

4 结论

利用二维有限体积数值模型，系统地研究了外源流水注入对河流地貌的影响，其结论包括以下3个方面：

1)外源流水注入后，接收其注入径流的河道必然会经历河流快速下切的过程。连接外源流水注入位置到主河道的支流逐渐从支流变成主河道，原先的部分主河道可能会退化为支流。河道被深切后，其支流受侵蚀基准点高程下降的影响，开始加速下切，进而整体流域高程下降。流域的平均高程依赖于地表隆升速率、河流下切系数和外源流水注入流量。

2)当外源流水注入流量较大时，下游流域会袭夺相邻流域。河流下切系数控制流域开始袭夺的时间。河流下切系数和外源流水注入流量主要控制流域袭夺速度。相邻流域被袭夺后，该流域的河流下切过程减缓，可能促使低起伏地形的形成。外源流水注入引起的流域袭夺过程和地表隆升速率无关。

3)河套古湖湖水下泄，导致了晋陕黄河的深切。晋陕黄河流域可能已经袭夺了古汾河的部分流域。