高文冰，姚聪聪，夏 岩，喻晓琳，张洪波，4
（1. 中国石油长庆油田公司勘探开发研究院，陕西 西安 710021；
2. 低渗透油气田勘探开发国家工程实验室，陕西 西安 710021；
3. 长安大学水利与环境学院，陕西 西安 710054；
4. 长安大学旱区地下水文与生态效应教育部重点实验室，陕西 西安 710054）
水文观测数据是开展工程水文设计和水资源规划与管理的基础，也是国际水文界一直以来持续关注的热点。近年来，随着国家水利事业的快速发展，我国水文监测站网体系日臻完善，目前已发展到12.1 万处，仅国家基本水文站就建有3 154处，实现了对大江大河及其主要支流、有防洪任务的中小河流水文监测的全面覆盖。然而值得注意的是，由于站点建设偏向人口聚集区、重点河湖保护区、重大工程辐射区以及洪涝灾害防控区，因此在西部许多欠发达或人口相对少的地区（如黄土高原）站点密度相对稀疏，无资料或资料短缺问题仍旧普遍存在，从而对区域工程建设、水资源规划与管理以及生态环境保护等方面造成较大影响［1］，极大地限制了西部地区生态文明建设与高质量发展的进程。
黄土高原是中华文明的发祥地，也是我国黄河流域生态保护与高质量发展以及生态文明建设的主战区。然而，受独特的自然地理条件和长期经济欠发达等因素的制约，黄土高原水文站网的密度一直相对稀疏，且由于区域空间异质性强烈，使得水文数据的短缺的影响被无形放大，一定程度上阻碍了区域水利事业和社会经济的迅速崛起。在黄土高原的无资料或资料短缺区，通常采用传统水文移用方法，但受空间异质性影响，设计结果和规划值常偏差较大，进而引发部分工程或规划的效能无法有效发挥。马莲河流域是黄土高原的一个典型流域，其水系广泛发育、支流交错，但水文站点分布较少，近2 万平方米的流域面积上仅有洪德、贾桥、庆阳、雨落坪四个站点，且支流站点只有一个，现有观测水文数据无法有效覆盖的区域较大，给流域水文设计和水资源规划与管理带来了极大的不便。特别是近几十年来，气候变化和人类活动的耦合影响更趋强烈，如高强度降水、地表覆被改善、董志塬综合治理等，使得马莲河流域在双重驱动下水循环及水文要素的时变性及空间异质性更趋显著，区域水文系统的复杂性亦不断攀升［2，3］，故当前探明马莲河流域不同区域水资源赋存及动态变化，已成为进一步夯实生态文明建设成果，保障区域社会经济高质量发展的关键，且极为迫切。
在现行的技术手段中，分布式水文模型是估算流域水资源空间分布格局的有效工具，已在变化环境下流域水文过程模拟中广泛应用。近年来，作为经典的分布式水文模型，SWAT模型在黄土高原的径流模拟中应用较为普遍［4］。如李文婷等基于SWAT 模型分析了黄河源区降水量及蓝绿水量的变化趋势［5］，李小兵利用SWAT 模型模拟了黄土高原秃尾河流域的月径流量［6］，刘闻采用SWAT 模型模拟了渭河流域关中段的径流变化［7］，庞佼基于SWAT 模型分析了黄土高原安家沟流域的月径流过程［8］。通过SWAT 模型构建，有效实现了有限站点径流资料向流域尺度地表水资源量时空分布格局的扩展，减小了水文移用过程的偏差。但是上述研究中也明确指出，由于黄土高原产流机制较为复杂，地下水对径流的补给强烈［9，10］，单纯的水文模型在模拟地表水-地下水水量交换方面能力极为有限，较难获得满意的结果。因此，将其与地下水模拟较为出色的MODFLOW 模型相耦合，可进一步减小模拟误差［11］，获得地下水补给强烈地区更为合理水文模拟结果。
鉴于此，本文拟聚焦黄土高原马莲河流域，通过构建和校准流域SWAT-MODFLOW 耦合模型，开展流域尺度的水文过程模拟，细化流域无资料或资料短缺地区不同地表覆盖情景下的水文过程，并对流域地表水资源量的时空分布格局和演化特征进行评价，以期为马莲河流域水资源管理和相关涉水工程设计提供数据参考，并为黄土高原其他无资料或数据短缺地区的径流模拟与水资源评价提供技术方法支撑。
马莲河地处陇东黄土高原，由发源于麻黄山的西川和发源于白玉山西麓的东川汇合而成。河流自西北流向东南，东川与西川汇合后称环江，环江与柔远川汇合后称马莲河。主河道全长374.8 km，流域面积19 086 km2，多年平均流量13.9 m3∕s，多年平均年径流量4.39 亿m3。流域属温带大陆性气候，年降水量330～617 mm，多年平均气温8.5～9.3 ℃。汛期降水主要集中在6-9 月，约占全年总降水量的70%，同期径流量约占全年的67%，非汛期径流量相对匮乏，尤其是枯水期水量严重短缺。流域内设有上游干流洪德站、支流柔远川贾桥站、下游干流庆阳站和流域出口雨落坪站4 个水文站点，具体位置如图1 所示。自1980-2020年，由于沙漠整治、退耕还林还草政策及城市化进程的推进，流域内发生了以自耕地向草地及林地转移为主体的植被覆盖改变，同时城乡∕工矿∕居民用地面积也逐渐增加，覆被变化及人类活动对流域径流的影响日趋强烈。
图1 马莲河流域地形、河流水系及下垫面变化图Fig.1 Topography， river system map and underlying surface changes of the Malian River basin
马莲河流域系黄土高原典型区，具有地表水与地下水交互强烈的特点。已有研究表明，马莲河流域地下水对径流的补给约占多年平均径流量的60%～80%，且由于流域内黄土层厚度较大，存在一定的空间异质性，使得降水入渗包气带和地下水系统的时间并不相同，具体表现为优先流和活塞流双模式运行［12，13］。不同空间位置、不同下渗模式导致降水经地下水系统返回河道的时间呈现出明显的差异化，也引发了流域降雨-径流关系的相关性不显著。图2 显示了流域内4 个水文站面降水量-年径流量的协同演化关系，统计结果表明上游、中游及支流站点的相关系数R2不超过0.17，出口站点的相关系数R2也仅为0.372 6，从侧面验证了以上对马莲河流域三水转化的科学认知。由此也可以发现，要厘清马莲河流域地表水资源量空间分布格局和演化特征，准确模拟区域地下水对径流的补给过程至关重要。
图2 4个水文站点降雨-径流相关关系图Fig.2 The relationship between rainfall and runoff at the four hydrological stations
鉴于以上分析，本文选择SWAT-MODFLOW 耦合模型实施马莲河流域的水文模拟。为了实现模型建立和校验，采用了马莲河流域DEM数据、土地利用数据、土壤数据、气象数据以及径流数据，数据来源如表1所示。其中，庆阳和雨落坪站作为模型率定与校验站，洪德和贾桥作为流域尺度水资源时空分布拟合效果的检验站。
表1 选用数据序列类型、长度、精度及其来源Tab.1 Selection of data sequence type， length， precision and its source
2.1 SWAT模型
SWAT（Soil and Water Assessment Tool）是美国农业部（USDA）农业研究中心开发的一个分布式流域水文模型［14］，包含水文过程、土壤侵蚀过程和水质多个子模块。水文过程基于水量平衡方程进行，分为陆面部分及水面部分。陆面部分中，地表径流计算多采用SCS 曲线法计算［15］，用于模拟不同土壤或土地利用类型情景下的多尺度径流过程；
壤中流计算考虑了地形、土壤有效含水量、土壤渗透性等因素，采用动态储存方式进行管理；
地下径流计算将流域地下水划分为潜水和承压水两部分，潜水直接汇入河道，承压水假定汇入流域外，仅当潜水超过指定水位阈值时，才能给河道补给水量；
蒸发计算考虑了植物蒸腾及土壤水分蒸发，植物蒸腾通过叶面指数和潜在蒸散发的线性关系推求，潜在蒸散发多选择Penman-Monteith 法进行计算［16］。水面部分为主河道及蓄水体的汇流演算，采用曼宁公式和马斯京根法计算。
2.2 MODFLOW 模型
MODFLOW 模型是美国地质调查局开发的一款模块化三维有限差分地下水流动模型，主要用于空隙介质中地下水流运动的数值模拟［17］。通过把研究区在空间和时间上进行离散，建立每个网格点的水均衡方程式，可列出各个离散格点水头之间的大型方程组并求解，进而得出研究区的水头分布。涉及的主要方程如下：
式中：h为水头，m；
K为含水介质的水平渗透系数，m∕d；
Kz为含水介质的垂向渗透系数，m∕d；
W为含水层的源和汇，1∕d；
Ss为自由面以下含水层储水率，1∕m；
t为时间，d。
2.3 SWAT-MODFLOW 耦合模型
SWAT-MODFLOW 耦合模型主要整合了SWAT 模型和MODFLOW 模型的优势，可有效解决地下水与地表水交互强烈流域的水文模拟问题，主要原理如图3 所示。其中，SWAT 模型可模拟地表水与土壤水运动过程，MODFLOW 模型可模拟地表水和土壤水与地下水的交互过程，模型耦合主要基于地表水模型最小单元（HRUs）和地下水模型最小单元（Cells）之间的映射进行。将SWAT 模型的最小水文响应单元（HRUs）离散分解成的具有特定地理位置信息的单个多边形（DHRUs），之后与MODFLOW 有限差分网格（Cells）相映射，从而在SWAT和MODFLOW 之间传递变量。地表水模型向地下水模型输入其最小网格模拟的深层入渗量，地下水模型向地表水模型输入最小单元模拟的河流排泄量，耦合模型数据传递示意图如图4 所示。模型耦合接口可透过QGIS软件的QSWATMOD插件实现［10］。
图3 SWAT-MODFLOW 耦合模型原理示意图Fig.3 Schematic diagram of SWAT-MODFLOW coupling model
图4 SWAT-MODFLOW 耦合模型数据传递示意图Fig.4 Schematic diagram of data transfer of SWAT-MODFLOW coupling model
2.4 模型校准
SWAT 模型的校准通过SWAT-CUP 软件中的SUFI-2 算法实现，MODFLOW 模型的校准通过调整水文地质参数进行，模型的校准精度采用决定系数R²和纳什效率系数NSE评价。R²描述了模拟值与实测值之间的线性相关性，越接近1 则相关性越好［18］。NSE取值在-∞～1之间，越接近1则可信度越高；
接近0表示结果接近实测值的平均水平；
远小于0则表示结果不可信。一般认为，当R2＞0.60 且NSE＞0.50 时认为模型的模拟结果较好［19］。
决定系数R²和纳什效率系数NSE的计算公式如下：
式中：Qm，i为实测流量，m³∕s；
Qp，i为模拟流量，m³∕s；
Qm，avg为多年实测平均流量，m³∕s；
Qp，avg为多年模拟平均流量，m³∕s；
n为实测时间序列长度。
3.1 SWAT模型构建
SWAT 模型建模的基本流程如图5 所示。研究设定流域集水面积阈值为20 000 hm2，添加庆阳水文站作为控制点，洪德、贾桥水文站为观测点，设置雨落坪站为出口点，将研究区划分为56个子流域，划分结果如图6（a）所示。采用多水文响应单元划分法确定水文响应单元，设定土地利用面积阈值、土壤类型面积阈值及坡度阈值为5%，将研究区划分为519 个HRU，划分结果如图6（b）所示。加载气象数据后即可完成SWAT 模型的初步构建。
图5 SWAT水文模型建立流程Fig.5 SWAT hydrological model establishment process
图6 马莲河流域子流域和水文响应单元划分Fig.6 Sub-basin and hydrological response unit delineation in Malian River Basin
3.2 MODFLOW 模型构建
将研究区潜水含水层概化为三维非稳定流，将流域边界视为零通量边界。模型率定期与验证期与SWAT 模型相对应，设定1月为一个应力期，时间步长为1 d。按照单元中心法对研究区水文地质模型进行空间离散，划分过程采用1 000 m×1 000 m的正方形网格，最终形成有效单元格40 592 个。参考马莲河流域水文地质图，结合地质与钻孔资料，根据不同岩性渗透系数经验值、不同岩性给水度经验值以及不同岩性贮水率经验值确定研究区域参数初值，根据含水层性质对研究区域划分，最终将马莲河流域划分为4个子区域，如图7所示。
图7 马莲河流域水文地质参数分区图Fig.7 Hydrogeological parameters zoning map of the Malian River Basin
3.3 SWAT-MODFLOW 耦合模型构建
通过在QGIS 插件中构建水文模型与地下水模型最小单元之间的映射关系来构建耦合模型，主要包括：创建SWAT-MODFLOW 链接；
创建swatmf_link.txt 文件；
运行耦合模拟程序三个步骤，具体过程如图8所示。
图8 SWAT-MODFLOW 耦合模型模拟过程Fig.8 Simulation process of SWAT-MODFLOW coupling model
3.4 模型率定及验证
选取1986-1988年为SWAT模型的率定期，1991-1992年为验证期，设置预热期长度2 a，模拟尺度为月尺度，对庆阳站和雨落坪站分别进行率定，流域最敏感的参数包括CN2、SOL_K、ESCO、EPCO、CH_K2、CH_N2等。SWAT-MODFLOW 模型参数率定时段与SWAT 模型保持一致，使用耦合模型联合出流后所形成的河道断面流量对MODFLOW 模型各水文分区的水文地质参数进行反复调整，最终得到流域不同分区的水文地质参数率定结果，如表2 所示。耦合前后各站点的拟合评价结果如表3所示，实测值和模拟值的拟合情况如图9所示。
表2 马莲河流域水文地质参数率定值Tab.2 Values of hydrogeological parameters in the Malian River Basin
表3 各站点模型拟合评价结果Tab.3 Evaluation indicators of fitting results of two models at the sites
图9 水文站点实测值和模拟值拟合情况Fig.9 Fitting of measured and simulated values at hydrological sites
由SWAT 模型的模拟结果可知，单一水文模型对流域丰水期的模拟效果较好，流量峰值基本对应，但枯水期河道基流的模拟效果较差，原因在于SWAT 模型在地下水模块采用的是集总式结构，处理地下水交互的方式是一维概念性的，并未考虑到不同子流域之间的地下水运动关系，导致地下水位动态变化的模拟效果欠佳，进而影响SWAT 模型的模拟精度［20］。将其与MODFLOW 模型耦合后，枯水期月均流量的拟合效果明显变好，R2与NSE系数也进一步提高，拟合效果总体较为满意。对模型校验站洪德站和贾桥站的实测与模拟结果进行对比可知，在对应的率定期和验证期，洪德站耦合模型的NSE系数分别为0.88 及0.79，贾桥站耦合模型的NSE系数分别为0.75 及0.88。由此证明，基于SWAT-MODFLOW 耦合模型开展流域尺度的地表水资源量空间分布模拟是可行的。
4.1 出口断面流量时空分布
4.1.1 出口断面流量空间分布
采用1980 年的土地利用数据模拟马莲河流域1965-2019年近天然条件下各子流域的出口断面流量，并绘制多年平均年径流量空间分布图［图10（a）］。基于设站的水文资料，以4个水文站点实测多年平均年径流量（变异点1996年以前）为依据，根据不同区间的控制面积推求各子流域的出口断面流量，并将推算结果与耦合模型模拟结果进行对比，如图10（b）所示。
图10 各子流域出口断面多年平均径流空间分布图及不同情景径流对比图Fig.10 Spatial distribution of multi-year average runoff at each sub-basin outlet section and comparison of runoff in different scenarios
分析各子流域出口断面流量模拟结果可知［图10（a）及10（b）］，不同子流域或区域出口断面流量的空间差异性较大，尤其是支流与干流的空间异质性十分显著，整体呈现为自西北向东南逐渐递增的趋势。若在支流和干流数据间实施传统水文移用法极有可能造成较大误差，导致水利工程设施安全或区域水资源规划面临威胁。
对比面积比推算结果和耦合模型模拟结果可知［图10（b）］：马莲河上游流域的径流推算结果与耦合模型模拟结果较为一致；
中游流域出口断面流量较小的子流域径流推算结果与耦合模型模拟结果较为一致，但出口断面流量较大的子流域偏差则较明显，主要表现为面积比径流推算结果对干流水量的高估；
下游流域模拟结果偏差不明显。从与四个水文站点实测值的对比来看，耦合模型模拟结果与实测结果更为接近，而面积比推算结果明显偏大，水资源量高估较为严重。综上，认为基于SWAT-MODFLOW 耦合模型模拟的流域水资源量是相对可信的，基本能较为可靠地反映无资料或资料短缺地区的地表水资源量，而在传统水文设计中常用的面积比推算法容易出现高估，导致设计偏差。
4.1.2 出口断面流量丰枯变化
根据近天然条件下的径流模拟结果，绘制56个子流域出口断面的年径流量分布箱图，如图11所示，包括上游东西川、支流柔远川、干流环江及下游马莲河四个区域。由图11 可知，近天然条件下马莲河流域不同子流域出口断面年径流量基本呈现丰枯交替的平稳状态：上游及支流所在区域出口断面水量偏少，丰枯差异不显著；
干流及下游所在区域出口断面水量丰富，丰枯变化明显；
不同区域的出口断面水量及丰枯差异较大，具有显著的空间异质性。根据各子流域出口断面年径流量的变差系数可知，环江所在区域的变差系数高于各子流域变差系数的平均值0.46，其他区域的变差系数则基本小于该值，因此认为环江所在区域各子流域出口断面年径流量的年际变化较大，其他区域的年际变化次于该区域。
图11 56个子流域出口断面年径流量分布箱图Fig.11 Annual runoff distribution box map of 56 sub-basin outlet cross-sections
4.2 产水能力空间分布
计算近天然条件下各子流域的产水能力，并将其用径流模数M 表示，绘制多年月平均径流模数的空间分布图（图12）。在图12中，径流模数越大的子流域，单位面积上的产流量越多，颜色越蓝。由图12可知，马莲河流域各子流域或区域的产水能力具有显著的空间异质性：近天然条件下上游东川和西川干流区域的产水量最多，多年月平均径流模数大于1.5 L∕（s·km2）；
环江中上游处产水量最少，多年月平均径流模数小于0.65 L∕（s·km2）；
柔远川、环江下游及固城川、合水川、九龙河处产水量适中，多年月平均径流模数在0.65～1.5 L∕（s·km2）之间。
图12 马莲河流域近天然条件下多年月平均径流模数空间分布图Fig.12 Spatial distribution of multi-year monthly average runoff modulus under near-natural conditions in the Malian River basin
4.3 下垫面影响下的地表水资源量变化
4.3.1 土地利用演化规律
马莲河流域1980年与2020年的土地利用分布如图1所示，1980-2020 年的土地利用变化及转移情况如表4 及表5 所示。根据各土地类型的空间分布及其面积比例可知，马莲河流域耕地和草地分布均匀，林地多分布在东南部，近40 年内上中游流域林地面积明显增加（332 km²），耕地面积减少（500 km²），“退耕还林还草”政策效果显著。城镇化方面，流域城乡∕工矿∕居民用地增加了114 km²，人类活动强度加大，对流域水文循环过程的影响将显著增大。由土地利用转移矩阵可知（表5），马莲河流域内的土地流转主要表现为以耕地向草地、林地、城乡∕工矿∕居民用地的转变以及以草地向林地的转变，综合表现为流域内耕地面积减少，林地、草地、城乡∕工矿∕居民用地面积增加。
表4 马莲河流域1980-2020年土地利用类型变化表Tab.4 Changes in land use types in the Malian River Basin from 1980 to 2020
表5 马莲河流域1980-2020年土地利用类型转移矩阵Tab.5 Land use type transfer matrix in the Malian River Basin from 1980 to 2020
4.3.2 下垫面影响下的地表水资源量
将2020 年马莲河流域的土地利用数据带入SWAT-MODFLOW 耦合模型，模拟现状下垫面条件下的流域地表水资源量，并绘制各子流域出口断面多年平均年径流量空间分布图［图13（a）］，绘制各水文站不同情况下的年径流量分布箱图［图13（b）］。由图13（a）及13（b）可知，下垫面影响下流域出口断面多年平均年径流量的空间分布与近天然条件相差不大，不同子流域或区域出口断面径流量的空间差异性显著。雨落坪站近天然条件下的年径流量略大于现状下垫面条件下的年径流量，且二者均大于实测年径流量，这表明下垫面变化对马莲河径流演化的影响相对有限，而以取用水为代表的人类扰动不容忽视，应主要驱动了实测径流序列的趋势性变化。其他站点结果与雨落坪站类似。
图13 现状下垫面条件下多年平均年径流空间分布图与不同情景下的站点年径流量箱图Fig.13 Spatial distribution of multi-year average runoff under current substratum conditions and box plots of annual runoff volume distribution under different conditions
4.3.3 直接与间接人类活动影响
为便于分析，将由下垫面变化导致的径流减少归为间接人类活动的影响，由近天然条件与现状下垫面下水量的差值得出；
由直接人类取用水活动导致的径流减少归为直接人类活动的影响，由现状下垫面下的模拟水量与实测水量的差值得出。统计二者对径流减少的相对贡献见表6。由表6 可知，不同站点控制流域范围内的直接和间接人类活动影响量不同，对马莲河流域径流变化的贡献也不同。以耕地向草地和林地转移的流域植被覆盖改变，对马莲河流域地表水资源量的影响有限，未发生实质性变化，而流域直接人类活动对径流减少的影响强烈，各站直接人类活动对径流减少的贡献均达到了95%以上。
表6 马莲河流域间接和直接人类活动对年径流量减少的贡献Tab.6 Contribution of indirect and direct human activities to the reduction of annual runoff in the Malian River basin
以黄土高原地区马莲河流域雨落坪水文站以上区域为研究区，通过构建SWAT-MODFLOW 耦合模型，开展了马莲河流域地表水资源量时空分布评价及其演化分析，为该流域无资料或资料短缺地区的地表水资源量评估提供参考。具体结论如下：
（1）SWAT-MODFLOW 模型在马莲河流域具有较好的适用性，能够用于评价马莲河流域地表水资源的时空分布。SWAT模型对马莲河流域洪峰模拟较为准确，但枯季拟合偏差较大，对旱区枯季水资源利用技术支撑有限。与MODFLOW 模型耦合后，枯季径流的准确性有了较大提升，出口站点NSE系数提高至0.92。
（2）马莲河流域地表水资源的空间差异性导致数据移用存在一定偏差，而水文设计中常用的面积比推算法也可能会高估区域地表水资源量，基于SWAT-MODFLOW 耦合模型获取的地表水资源量空间分布成果更具有参考价值。
（3）流域内不同子流域出口断面径流的空间异质性显著，整体呈现自西北向东南逐渐递增的趋势；
不同子流域出口断面流量基本呈现丰枯交替的平稳状态；
环江所在区域各子流域出口断面年径流量的年际变化较大，其他区域的年际变化次于该区域。
（4）流域内不同子流域或区域产水能力的空间异质性显著，近天然条件下上游东西川干流区域产水能力最强，环江中上游区域产水能力较弱，柔远川、环江下游及固城川、合水川、九龙河处产水能力介于二者之间。
（5）近40 年来马莲河流域的覆被情形发生了耕地面积减少，林地、草地、城乡∕工矿∕居民用地面积增加的改变，以耕地向草地和林地转移的流域植被覆盖改变对马莲河流域地表水资源量的影响有限，而以取用水为代表的流域直接人类活动对径流减少的影响强烈。
猜你喜欢 马莲径流量水文 马莲，马莲散文选刊·下半月(2023年6期)2023-06-29继往开来 守正创新——河北省水文工程地质勘查院河北地质(2021年3期)2021-11-05马莲湖辽河(2021年4期)2021-08-27我家有一墩马莲辽河(2021年4期)2021-08-27水文河南水利年鉴(2020年0期)2020-06-09水文水资源管理河南水利年鉴(2020年0期)2020-06-09好学的马莲菜小学生导刊(2018年20期)2018-11-30水文比拟法在计算河川径流量时的修正水利规划与设计(2017年5期)2017-06-09水文河南水利年鉴(2017年0期)2017-05-19SCS模型在红壤土坡地降雨径流量估算中的应用水土保持通报(2014年5期)2014-06-09
