基于VMOD模型的若尔盖泥炭沼泽地下水数值模拟

鲁瀚友，李志威，2①，胡旭跃，3（.长沙理工大学水利工程学院，湖南长沙 404；2.青海大学省部共建三江源生态和高原农牧业国家重点实验室，青海西宁 8006；3.水沙科学与水灾害防治湖南省重点实验室，湖南长沙404）

泥炭湿地萎缩现象近年来在国内外普遍发生，造成湿地生态、水土资源面临威胁［1-2］。若尔盖泥炭沼泽曾达4 600 km2［3］，是全球最大的内陆高寒泥炭分布区。如今该地区自然沟道和人工沟渠密布，致使地下水位下降、地表水流失、泥炭层退化、有机碳氧化分解，宏观上表现为泥炭沼泽面积大幅度萎缩，现存面积仅约2 200 km2［4］。影响泥炭湿地萎缩趋势的重要指标是地形地貌［1］。2011—2017年若尔盖沼泽湿地的野外调查表明，若尔盖高原的泥炭地有2种典型的自然沟道：沟道切穿泥炭层和沟道未切穿泥炭层。因此，研究不同沟道对泥炭地地下水运动的影响对于揭示泥炭沼泽萎缩机制具有重要科学意义。

国内学者对若尔盖高原的地下水变化也有一定关注，如王焱等［5］在朗州公路附近进行地下水的抽水实验，通过有限差分法计算若尔盖地下水的平均渗透系数与入渗率；李丽等［6］证实了地下水位降低是若尔盖有机碳与全氮流失的重要原因；LI等［7］建立In⁃teractive Groundwater模型研究沟渠排水对湿地退化的影响。对于若尔盖泥炭沼泽萎缩及地下水疏干，自然沟道的排水作用同样较显著［8］。泥炭地含水率高，在自重作用下沿沟道发生垮塌和崩岸，使泥炭地的沟道横向侵蚀速率高，在一定程度上造成泥炭沼泽面积减少［9］。另一方面，沟道纵向的溯源下切，向泥炭沼泽内部延伸和发展，从而加速封闭泥炭沼泽的地表水及地下水流失［8］。自然沟道对于泥炭地的排水作用既是泥炭湿地萎缩的驱动力，也是当前若尔盖流域侵蚀过程与物质输运的集中体现。然而，受降雨频率与大小的影响，泥炭地的地下水变化较快，波动幅度相对较大，加上缺少地下水的原位监测，使得泥炭地的地下水数据缺失，其地下水运动的定量分析较困难。为此，采用野外泥炭地的地下水原位观测和地下水数值模拟相结合的方法，是研究泥炭地沟道对泥炭沼泽地下水运动的影响是一个可行的途径。Visual MODFLOW（VMOD）软件是一款用于模拟地下水运动的专业软件，在国际上广泛应用于地下水及溶质的运移研究［10-11］，如MELIKADZE等［12］用此软件计算地下水在含水层中的滞留时间；CHETHAN等［13］通过建模模拟降雨对缺水地区地下水的改变。VMOD软件在国内也有一些应用，如刁维杰等［14］运用此软件对地下水资源进行了优化配置。该软件也用于泥炭地及其沟道的数值模拟，如BUJAKOWSKI等［15］运用VMOD软件研究了泥炭覆盖对河谷水文地质条件的影响；CRUZ⁃FUENTES等［16］采用此模型研究了西班牙某干旱地区沟渠灌溉的方案。因此，采用VMOD软件研究自然沟道对于泥炭地的地下水运动具有可行性。

笔者以若尔盖高原黑河上游麦曲支流局部泥炭地的地下水野外观测为基础，通过野外实测数据，采用VMOD Flex 2015软件建立切穿泥炭层和未切穿泥炭层2种情况下的泥炭地地下水模型。根据实际情况设计了不同地形坡度，共计20组工况，模拟沟道对泥炭地地下水运动的影响及水力梯度，有助于认识2种典型自然沟道对泥炭地的地下水及湿地萎缩的影响。

1 研究区域与研究方法

1.1 研究区域与野外观测

研究区位于黄河源支流黑河流域的哈曲上游，属于红原县色地镇附近的泥炭沼泽湿地。该地区多年平均降雨量为765 mm，地势低平，降雨汇流于此，排泄不畅形成湿地。1950年代以来，由于人工开挖沟渠和自然沟道的溯源侵蚀，大片湿地被疏干变成草原，沼泽面积与历史时期相比萎缩近50%。黄河源干流的溯源下切起于180 Myr BP，至今仍保持1～2 mm·a-1的下切速率［17-18］，黄河河道以“U”型弯道切穿若尔盖盆地，促进若尔盖黑河、白河流域中自然沟道的横向侵蚀和溯源下切，从而导致地表水及地下水横向排水，形成沟道两侧的疏干带［7］，这一点为实地考察所证实。研究区域选取自然沟道边缘，观测点1（32°59′11.130″N，102°59′17.922″E）和观测点2（32°59′16.212″N，102°59′14.250″E）分别代表自然沟道切穿泥炭层和未切穿泥炭层（图1）。研究区覆盖约1.5 m厚的草本质泥炭层，是密集植被根系的浅变质炭化的多孔介质，干密度一般为0.2～0.4 g·cm-3，泥炭的质量含水率约为200%，泥炭下层为粉砂层，50%粒径（d50）约为0.039 mm。不同的沟道深度导致2种沟道的排水能力及其对地下水运动的影响也不同。

2个观测点的原位地下水观测设计相似，均由3排共9个点组成。这些测量孔3排平行于沟道，3排垂直于沟道，且受局部地形影响，每个观测区9个点间距离各不相同。测压管由一个1.5 m长的铁管组成，铁管外径2 cm，内径1.8 cm，在管底部10 cm范围内有4排穿孔以便地下水通过。每个测量管分别设置于地表以下100、70和40 cm位置，在每个深度插入钢卷尺进测量管直至管内水面，以此测量地下水的水头（water head），水头测量误差控制在3 mm以内。每测1次，测压管垂直移动后要等8～17 h，以待水头达到稳定。通过野外多次观测，8 h后该地区水头可达到平衡状态。测量水头之后，通过测量钻孔内地表到地下水的距离，记录地下水位值。测量方法与FISHER等［19］相似，水位值记录3次，1次降雨前和2次降雨后。泥炭地表面地形变化较大，使用全站仪测量当地9个位置点的高程。每个观测区共有27个地下水头值、9个高程值，以及3个不同时刻共27个水位值。

图1 研究区域原位测量位置示意Fig.1 Study area and location of two measuring points

除此之外，采用蠕动泵测量了35、72和100 cm的渗透系数K，渗透系数采样位置在图1标注的观测点1与观测点2，测量渗透系数的方法为CHASON等［20］于1986年基于静水时差方法提出。先记录起始水头，使用蠕动泵抽水，抽水量约为87 cm3，抽水深度约10～15 cm，抽水时间取决于水头高度，2～4 min不等。抽完水后，连续测量水头变化直至回升到起始高度，通过水头变化时间计算渗透系数［20］。

式（1）中，K为渗透系数，m·s-1；r为测压管内半径，m；T0为水头恢复总时间，s。有时候水头并不能完全恢复，所以测量的恢复时间T0可能会小于实际所需时间。考虑到这一点，初始阶段水位回升遵循线性趋势，对后面的点进行拟合线性回归，然后用线性模型预测总恢复时间，带入式（1）计算渗透系数［21］。各深度渗透系数均测量3遍并取平均值以减少误差（表1）。

表1 渗透系数测定Table 1 Determination of hydraulic conductivity m·s-1

根据2016—2017年野外实测结果，泥炭层的渗透系数约为 3.89×10-6～3.91×10-5m·s-1，而根据ZHENG等［22］的研究，粉砂渗透系数经验值范围为1×10-9～2×10-5m·s-1，故粉砂层渗透系数取值为10-8m·s-1。

1.2 模型建立与工况设计

1.2.1 模型建立

在VMOD模型中导入实测地形点，采用Natural Neighbors差值法生成地形面，通过定义模型结构体功能把地形面生成土体，最后定义土体的各项参数。对于该文讨论的物理情景需要建立2种模型进行比较。以观测点1的实测数据为基础，建立沟道切穿泥炭层模型，以观测点2的实测数据建立沟道未切穿泥炭层模型（图2）。

图2（a）模型宽9 m，长10 m，沟道宽2 m，深2m，沿x轴正方向坡度为0.5。表层为1.5 m的泥炭沼泽，其下为粉砂的弱透水层，弱透水层是含水层重要组成部分，具有低渗透性和高储水性，其渗透系数一般小于 10-8m·s-1［22-23］，降雨入渗后水流主要通过泥炭层运动。图2（b）沟道深度0.5 m且未切穿泥炭层，9个观测点的位置受地形影响不同。

泥炭层的地下水水面一般在地表附近，该模型研究以潜水流动为主。二维潜水流方程是根据Du⁃puit假设和质量守恒原理推导的潜水非恒定流动Boussinesq方程。由此基础上不采用Dupuit假设可得到三维流潜水流动的一般方程为

式（2）中，Kx、Ky和Kz为渗透系数分量，m·s－1；H为水头，m；Ss为贮水率，m-1；w为源汇项，d-1；t为时间，s。

通过VMOD模型计算20种工况条件，这些工况都是根据实际观测改变系数实现。实际工况有2种，即切穿泥炭层和未切穿泥炭层情况。这2种工况由于观测地点不同，井的观测位置改变（图3），同样常水头与沟道水头的值也随实测值发生变化。式（2）给出了地下水潜水流运移方程，但并不能确定具体运动状态，需要加上定解条件。

图2 泥炭地的三维数值模型示意Fig.2 Sketch of 3D simulation model in peatland

图3 边界条件设定Fig.3 Definition of boundary conditions

1.2.2 边界条件设定

边界条件是定解条件的一种。侧向常水头边界是模型的上边界，沟道内有水流，属于模型的下边界，表达式为

式（3）中，φ1（x，y，z，t）为边界Γ1上的已知函数，该模型中的边界均属于第1类边界条件。

模型区域的边界条件包含沟道与常水头，沟道宽2 m，方向向右。常水头边界在坡顶与远离沟道的侧边，代表地下水来水方向，取值为。坡度方向与沟道同向，2种情况下9个观测井位置受地形影响各不相同（图3）。坡脚位置为地下水去水方式之一，沟道切穿泥炭层情况下，地下水位取值在其地表以下，沟道切穿泥炭层情况下，地下水位取值在其地表以下。其余未设置边界条件的边缘，默认为0流量边界［24］。下文研究以模型内9个观察点数据展开，此小区域内有流量的出入，且距坡脚相对较远，故能使用以上边界条件分析瞬时局部水流趋势。

非稳定流问题需要初始条件，对于地下水模型来水，初始条件指0时刻水头值，模型中默认的初始水头为0 m。

地形的影响主要指坡度i改变对地下水的影响，实测点的坡度i=0.5，这里选取以0.01为等距间隔，即从0～0.1的坡度共10组模拟地下水变化情况。

综合有无沟道、坡度的不同组合，此次共模拟计算了20种工况。每种工况除自变量参数，其他参数与原始值相同，不考虑河岸形状变化对河岸的影响。此次模拟计算的时间为稳态。有限差分网格划分为50×45个，并对观测井附近网格加密1倍。

2 结果与讨论

2.1 实际工况下模型模拟结果

图4为沟道切穿泥炭层观测点1与未切穿泥炭层观测点2的VMOD模拟运行结果。图4（a）中地下水在岸上递减相对较缓慢，观测区域水头每米约下降0.1 m。靠近边坡水头递减迅速，沟道切穿泥炭层后地下水只能通过边壁向下汇流或通过弱透水层入渗。这使得沟道边壁附近地下水的水力梯度达到0.89，等水头线密集程度是岸上9倍。图4（b）水头变化较均匀，观测区水头约每米下降0.047 m，沟道未切穿泥炭层时沟道影响相对较小。以图4（a）中5号点为例，该点约66%的地下水偏向于垂直沟道方向，而图4（b）中5号点这个数值约为59%，证明了沟道的存在促使泥炭层产生更大的横向水力梯度，使沟道具有排水的作用，而且沟道切穿泥炭层后排水效果更加明显。

图4 地下水水头等值线图Fig.4 Contour plot of water head of groundwater in peat layer

图5是沟道附近某点的来水路径。泥炭地地下水流运动受沟道与地形影响。地形影响下图5（a）的迹线有向坡脚运动的趋势，沿地形坡度方向水力梯度为0.043。在沟道影响下地下水偏向沟道方向运动，沿垂直沟道方向水力梯度为0.084。图5（b）沿地形坡度方向水力梯度为0.031，沿垂直沟道方向水力梯度为0.044。图5（a）沟道内水位较低，且并未与泥炭层相连，岸上水头较高，在边坡上有一个陡降的过程，地下水位在边壁下降近1 m。图5（b）沟道内水位相对较高且与泥炭层连接，迹线较均匀。

图5 地下水流动轨迹追踪Fig.5 Path line of groundwater movement

建立模型采用的地形和边界条件都需简化且形式不定。数值模型边缘会选择特征边界，如河流、山脊、湖泊等以减少不确定因素，但边界条件往往是根据需要和对水流运动的判断建立的，具有一定的不确定性［23］。因此，需要把建立的泥炭地地下水模型的模拟结果与实测值进行对比，以检验其可靠性。图6是观测值与数值模拟的计算值的对比。图6（a）的1、2、3、4、6号点观测值与运算结果相对吻合，5、7、8、9号点相对离散，观测值大于计算值。分析可能的原因是：（1）受地质、降雨、植物根系等影响，实际情况下沟道附近地下水运动更为复杂；（2）模型是在实际自然条件简化后建立的，并不能完全模拟地下水的实际流态；（3）泥炭地的地下水水位及水头是时变的，缺少在时间上高分辨率数据作为输入条件。对于沟道未切穿泥炭层的情况，图6（b）的结果类似，但误差在0.1 m以内，小于图6（a）。图6（a）的平均残差为-0.15 m，均方根误差为0.21 m；图6（b）的平均残差为0.022 m，均方根误差为0.045 m，说明模拟的计算值较为正确，可反映泥炭层地下水观测的真实状况。

图6 观测值与模拟值的比较Fig.6 Comparison between observed water head and calculated water head

2.2 地形对地下水的影响

在原始工况的条件下，按表1从0.01～0.10取10个坡度。对研究区9个位置地下水位进行计算，2、5、8号观测点为垂直于沟道方向断面（图7），4、5、6号观测点为地形坡度方向断面（图8）。

地形坡度的变化能引起地下水流动的改变，2、5、8号观测点是沟道方向的断面上3个点（图7）。图7（a）的3个点坡度每增加0.01，水头平均下降0.4 cm，2号点与5号点的平均水头差为7.49 cm，5号点与8号点平均水头差11.39 cm。这表明沿坡度方向，越往坡脚流向，沟道的地下水水头下降比例越大，且不受坡度改变的影响。对于沟道未切穿泥炭层，图7（b）的线间距同步扩大，从0.01坡度时的2 cm，到0.1坡度时扩大至7 cm。说明坡头到坡脚流向沟道的地下水比例变化不大，但随坡度的增加地下水向沟道出流加快。

对比图8（a）垂直沟道方向的3个点，2点间水头差均为0.07 cm，正如前文描述的沟道切穿泥炭层时岸上递减缓慢且均匀，在边壁附近陡降。沟道未切穿泥炭层时，地下水沿垂直沟道方向先快速后缓慢下降，0.05坡度时5～6号点地下水下降速度比4～5号点慢35%左右。

图7 垂直于沟道方向的坡度-水头关系曲线Fig.7 Slope‐head relation curve in perpendicular to gully direction

2.3 水力梯度对比

表2是由插值得到的3个点在不同方向上的水力梯度。2号点远离沟道，8号点靠近沟道，两者地形坡度方向上的水力梯度之差在沟道切穿泥炭层情况下是0.023 2，沟道未切穿泥炭层情况下0.003 8，前者情况下沟道排水能力更强。垂直沟道上3个点水力梯度值变化不大，切穿时距沟道3 m的点水力梯度均在0.85左右，未切穿时距沟道2.5 m的点则在0.45附近，沟道对等距点的影响相同。5号观测点位于观测区域中心，对比5号点横纵水力梯度得到图9。切穿沟道泥炭层的5号点在垂直方向的水力梯度偏小，表明其位置地下水相比4号更偏向于地形坡度方向。未切穿沟道泥炭层的5号点在地形坡度方向的水力梯度偏小，也说明其位置地下水相比2号更偏向于垂直沟道方向。

图8 地形坡度方向坡度水头关系曲线Fig.8 Slope‐head relation curve along terrain slope direction

表2 典型方向上水力梯度值Table 2 Hydraulic gradient in typical directions

2种典型沟道情况下垂直沟道方向的水力梯度均大于地形坡度方向的水力梯度（图9）。沟道切穿泥炭层时2个方向水力梯度差值平均为0.04，沟道未切穿泥炭层时这个差值为0.014，表明沟道可加速两岸沼泽的疏干，而且沟道切穿泥炭层时其作用更强。坡度对水力梯度的影响在2种情况下是不同的，沟道切穿泥炭层时垂直沟道方向的水力梯度在0.01～0.1坡度范围内只减少0.003，而沟道未切穿泥炭层在同等情况下增加0.052。沟道切穿泥炭层且未与两岸泥炭层直接接触，地下水在边壁附近陡降，通过边壁流入沟道，地形坡度对水力梯度影响小。沟道未切穿泥炭层时与泥炭层连接，沟道底部以下1 m依旧是泥炭层，地下水的变动受坡度影响较大。

图9 垂直于沟道方向和沿地形坡度方向的水力梯度对比Fig.9 Comparison of hydraulic gradient in perpendicular to gully and along the terrain slope

结合实测资料，利用VMOD模拟结果的分析表明，沟道附近泥炭地地下水水力梯度偏向沟道方向，证明沟道确实具有疏干地下水的能力。这与LI等［7］的结论一致，认为排水沟渠对湿地退化有重要影响。研究定量分析了2种典型沟道附近垂直沟道方向水力梯度的大小，沟道切穿泥炭层时5号点位置水力梯度是沟道未切穿泥炭层的1.79倍，证明沟道从未切穿泥炭层情况逐渐溯源下切直至切穿泥炭层的过程会加速地下水的流失［8］。通过0.01～0.10地形坡度的对比证明沟道未切穿泥炭层情况受地形影响大，而沟道切穿泥炭层后地下水受地形影响小，这与沟道内水流与泥炭层是否相连有关。

该研究还存在一些不足之处，例如0流量边界使计算存在误差，常水头位置仅代表瞬时水位，并非实际的动水头边界。但前人证明VMOD模拟泥炭地地下水有一定可行性［15，25］。这些问题有待增加实测野外数据和深入理论分析以进一步解决。

3 结论

基于若尔盖局部泥炭沼泽地的地水下实测数据和VMOD模型，建立当地2种典型自然沟道影响下泥炭地的地下水数值模型，对比了若尔盖2种典型沟道的排水能力。结果表明，湿地的保护可以针对排水作用更强的沟道切穿泥炭层区域进行。主要结论如下：

（1）观测区超过50%的地下水流入沟道，沟道具有局部疏干排水的能力。若尔盖泥炭地的自然沟道不仅排走部分地表水，还在非降雨期疏干沟道两侧的地下水，使其两侧的地下水潜水位降低，从而形成泥炭沼泽疏干带。

（2）自然沟道切穿泥炭层比未切穿泥炭层的距沟道2 m位置垂直于沟道的水力梯度增大79%，即排水能力显著增强。在沟道溯源下切作用下，未切穿沟道逐渐向切穿沟道发展，疏干带在沟道两侧放射式扩张。沟道的存在及发育促使若尔盖泥炭沼泽的地下水不断流失。

（3）泥炭地的沟道排水能力始终大于小幅度地形坡度影响，即自然沟道的排水功能基本不受局部地形变化而改变。这使得在沟道密布的若尔盖高原自然沟道排水现象普遍存在，其两侧的疏干带连结成网络，从而大范围降低若尔盖泥炭湿地的地下水位。