常布辉 李根东 苏小飞 王会永 王艳明 温乐

关键词：ＭＩＫＥＳＨＥ；生态环境约束；水权交易；沈乌灌域；河套灌区

０引言

水权交易的本质是对水资源进行时空的二次分配［１－２］，因此可能会引起一系列生态负外部效应，尤其是在降水量少、以地表水为主要水源的西北内陆灌区［３］。在水资源短缺加剧的背景下必然存在着复杂的社会利益冲突［４－５］，因此国内外学者对水权交易的外部生态环境效应进行了大量研究。

张忠［６］结合东阳义乌地区水权交易的典型案例，围绕水权交易过程中产生的外部效應进行了阐释，研究结果表明当初始水权已经被不同的水资源利用主体完全分配，针对生态安全等新增的生态用水问题就必须通过水权的转让或交易解决。基于生态经济学的Ｅｍｅｒｇｙ理论和外部理论，Ｌｙｕ等［７］在分析农业水权交易导致出让水权者产生生态经济损失后，提出了一种包括社会、经济和生态环境补偿的突发事件量化方法，同时发现中国水权交易市场的生态经济补偿呈持续上升趋势，２０１３年至２０１７年由９４００万元上升到了１．１４亿元。

２０世纪８０年代国外出现了水权交易与生态环境的关系研究［８］。例如：美国加利福尼亚州的ＳａｎＪｏａｑｕｉｎ河谷将农业水权转向工业，导致灌溉补水不足，使得砷等有毒物质得不到有效的降解，最终导致Ｋｅｓｔｅｒｓｏｎ国家野生动物园庇护所发生污染并关闭［９］。洛杉矶通过水权交易获得蒙诺湖支流取水权后，导致该湖的水资源量大幅度下降，生态平衡被打破，出现了大量野禽死亡的现象［１０］。Ｅｔｃｈｅｌｌｓ等［１１］发现水资源在维持河流生态环境方面的价值相当于或大于其消费性利用的价值，水权交易不仅降低了水资源的质量，而且可能引起地下水位下降、土壤盐碱化等问题，从而降低水资源的再生能力。Ｌａｄｓｏｎ等［１２］从环境角度解释了改变水权交易取水的时间和地点所带来的外部环境效应：在空间上，由于不同取水点和运输渠道的漏损率不同，因此在进行水权交易时对第三方产生的影响程度也不同；在时间上，干旱期取水比非干旱期取水所带来的外部影响大。Ｈｅａｎｅｙ等［１３］考虑时间因素的影响，把水权交易外部效应分为供给可靠性、取水时效性、储蓄水及取水费用、水质４类。Ｈａｎａｋ［１４］讨论美国加利福尼亚州地下水区域间交易对第三方回流的影响，认为出口禁令可以有效限制地下水输出，减轻地下水位的下降程度。

基于以上研究可知，虽然目前国内外学者在水权确权方法、水权交易模式和水权交易外部影响等研究领域已经获得了较为丰富的研究成果，但在水权交易外部影响方面，国内对水权交易的经济社会效应关注较多，对生态环境效应的研究相对较少，水权交易对生态环境的作用机理、影响程度及应对措施等仍需要深入研究［１５］。

１研究区概况与数据来源

１．１研究区概况

河套灌区位于内蒙古自治区巴彦淖尔市，是全国３个特大型灌区之一，也是我国最大的一首制自流引黄灌区。沈乌灌域位于三盛公水利枢纽西北部，乌兰布和沙漠穿沙公路以北，磴口县与杭锦后旗行政边界以南，河套总干渠及乌拉河干渠以西，狼山冲洪积坡地边界以东，总面积约１８．６８万ｈｍ２。沈乌灌域属于温带大陆性干旱气候区，降水稀少，蒸发强烈，干燥多风，昼夜温差大，日照时间长。２０１９年沈乌灌域通过实施节水改造向外转让水权１．２亿ｍ３，现状引水由原来的４．２６亿ｍ３下降到３．０６亿ｍ３。

１．２数据来源

１）地下水埋深。课题组在沈乌灌域布设了４７眼地下水埋深监测井，监测获取了２０１６—２０１９年的五日地下水埋深。

２）归一化植被指数（ＮＤＶＩ）。选用２０１６年８月２１日Ｌａｎｄｓａｔ＿５和２０１９年９月４日Ｌａｎｄｓａｔ＿７卫星数据，分辨率为３０ｍ，云量少，使用ＥＮＶＩ５．２经过辐射定标、大气校正、裁剪等预处理后计算得到研究区自然植被分布区的ＮＤＶＩ值。

２技术路线

技术路线见图１。

３地下水埋深与ＤＮＶＩ之间的关系

３．１数据处理

ＮＤＶＩ数据是３０ｍ×３０ｍ分辨率的栅格数据，而地下水埋深数据来源于分散在研究区的地下水监测井，为分析二者之间的关系，需要对数据进行预处理，以保证两个变量数据的匹配性。数据预处理的方法有３种。

１）空间网格作为基本单元。本研究分析空间网格设置为１ｋｍ×１ｋｍ，需要利用ＡｒｃＧＩＳ工具将３０ｍ×３０ｍ的ＮＤＶＩ数据重采样至１ｋｍ×１ｋｍ。同时，利用克里金方法插值到对应空间网格。通过匹配共获取５１７３６组数据样本，基于空间网格的地下水埋深与ＮＤＶＩ散点图（图略）可以看出，两者之间的关系没有明显规律。

２）泰森多边形法划分基本单元。根据离散分布的地下水埋深监测井，得到地下水埋深监测井的泰森多边形，如图２所示。利用ＡｒｃＧＩＳ的空间统计工具计算每个泰森多边形内ＮＤＶＩ平均值，２０１６—２０１９年３９个（部分监测井数据不连续）多边形区域一共确定１５６组样本。

利用泰森多边形获取的１５６组样本绘制地下水埋深与ＮＤＶＩ关系图（见图３），由图３可以看出，地下水埋深与ＮＤＶＩ成较明显的线性关系，随着地下水埋深的增大，ＮＤＶＩ明显下降。

３）土质单元划分基本单元。土壤类型是影响植被生长的重要因子，根据世界土壤数据库土质类型，沈乌灌域可划分为６６个土质单元。利用ＡｒｃＧＩＳ的空间统计工具计算每个土质单元ＮＤＶＩ平均值和地下水埋深，２０１６年各土质单元ＮＤＶＩ平均值和地下水埋深如图４所示。

２０１６—２０１９年６６个土质单元（其中１个土质单元无ＮＤＶＩ数据）一共确定２６０组样本。图５为２６０组样本的地下水埋深与ＮＤＶＩ关系散点图，由图５可知，地下水埋深小于３．４ｍ时，地下水埋深与ＮＤＶＩ没有关系；地下水埋深大于３．４ｍ时，地下水埋深与ＮＤＶＩ成明显的线性关系。

３．２结果分析

通过数据预处理可知，基于泰森多边形基本单元和土质基本单元均得到地下水埋深与ＮＤＶＩ的线性关系，利用皮尔逊相关系数法和多元线性回归法分别对其进行相关关系分析。

１）泰森多边形基本单元的分析结果。２０１６—２０１９年３９个泰森多边形基本单元的地下水埋深与ＮＤＶＩ的皮尔逊相关系数为－０．３７５，且通过显著性水平为１％的显著性检验。地下水埋深与ＮＤＶＩ的相关性属于负相关，即地下水埋深越小，泰森多边形基本单元内的ＮＤＶＩ平均值越大。

以泰森多边形基本单元的ＮＤＶＩ平均值作为被解释变量，将地下水埋深、土壤含盐量、土壤沙含量、土壤淤泥含量、土壤黏土含量、土壤有机碳含量、土壤电导率、土壤碎石体积百分比等要素因子作为解释变量，构建多元线性回归分析模型，得到以下结论：多元线性回归模型的拟合优度为０．１４，表明选取的解释变量对ＮＤＶＩ值的解释程度为１４％；地下水埋深的检验统计量为－３．８２，Ｐ值小于０．０５（Ｐ在ｔ检验中用于对各变量系数显著性进行检验，Ｐ＜０．０５表示显著，Ｐ＞０．０５表示不显著），表明泰森多边形单元内的地下水埋深对ＮＤＶＩ平均值具有显著的影响；地下水埋深的偏回归系数为－０．００８６，说明在其他解释变量不变的条件下，地下水埋深每下降１ｍ，泰森多边形单元内ＮＤＶＩ平均值会降低０．００８６。

２）土质基本单元的分析结果。当地下水埋深小于３．４ｍ时，土质基本单元地下水埋深与ＮＤＶＩ没有关系；当地下水埋深大于３．４ｍ时，土质基本单元地下水埋深与ＮＤＶＩ成显著线性关系。因此，将地下水埋深３．４ｍ作为临界值，将２０１６—２０１９年６６个土质单元的２６０组样本分别进行分析。在２６０组样本中，有２２０组样本地下水埋深小于３．４ｍ，记为Ａ组；４０组样本地下水埋深大于３．４ｍ，记为Ｂ组。

Ａ组地下水埋深与ＮＤＶＩ的皮尔逊相关系数为０．０２，且未通过显著性检验，表明当地下水埋深小于３．４ｍ时，地下水埋深与ＮＤＶＩ不存在相关性。Ａ组多元线性回归的结果显示，地下水埋深的检验统计量为０．０５（Ｐ＞０．０５），没有通过显著性检验；而土壤盐分和土壤沙含量、土壤黏土含量等变量通过了显著性检验（Ｐ＜０．０１）。结果表明：当地下水埋深小于３．４ｍ时，地下水埋深对ＮＤＶＩ没有显著的影响，而土壤盐分和土壤质地对ＮＤＶＩ具有显著的影响。

Ｂ组地下水埋深与ＮＤＶＩ的皮尔逊相关系数为－０．５５９，且通过了显著性水平为１％的显著性检验，表明当地下水埋深大于３．４ｍ时，地下水埋深与ＮＤＶＩ存在显著负相关性。Ｂ组多元线性回归模型的拟合优度为０．５８，表明选取的解释变量对ＮＤＶＩ值的解释程度为５８％；地下水埋深的检验统计量为－２．３５（Ｐ＜０．０５），表明当地下水埋深大于３．４ｍ时，地下水埋深对土质基本单元ＮＤＶＩ平均值具有显著的影响；地下水埋深的偏回归系数为－０．０８０７，说明在其他解释变量不变的条件下，当地下水埋深大于３．４ｍ时，地下水埋深每下降１ｍ，土质基本单元ＮＤＶＩ平均值会降低０．０８０７。

通过对比泰森多边形基本单元和土质基本单元Ｂ组的多元线性回归结果可知，土质基本单元Ｂ组的多元线性回归模型的拟合优度远高于泰森多边形基本单元。由土质基本单元Ａ组的多元线性回归结果可知，土壤质地在一定条件下对ＮＤＶＩ具有显著影响，因此泰森多边形基本单元的分析结果准确性会受区域内土壤质地差异性的影响。同时考虑泰森多边形基本单元多元线性回归模型的拟合优度相对较低，本研究认为土质基本单元的分析结果具有更高的可信度。

综上所述，利用皮尔逊相关系数法和多元线性回归法分析地下水埋深与ＮＤＶＩ的关系，得到以下结论：当地下水埋深小于３．４ｍ时，地下水埋深与ＮＤＶＩ没有相关性；当地下水埋深大于３．４ｍ时，地下水埋深与ＮＤＶＩ具有显著负相关性，在其他因素不变的条件下，地下水埋深每下降１ｍ，土质基本单元的ＮＤＶＩ平均值会降低０．０８０７。

４模型构建及验证

結合沈乌灌域水循环特点，构建基于ＭＩＫＥＳＨＥ模型的沈乌灌域分布式水文模型。鉴于灌域内水分水平运动比较弱，而垂向运动较强的特点，模型主要由气象模块、地表模块（简化处理）、不饱和带、饱和带以及一维河流模块组成。

４．１模型结构及原理

ＭＩＫＥＳＨＥ模型由数个独立的且相互联系的基于物理过程的模块构成，每个子模块用于一个主要的水文过程的描述，根据不同的模拟要求，这些模块可以互相分离也可以综合起来应用，分离开来可以分别描述水文循环的各个过程，综合起来可以描述整个流域的水文循环过程。

１）蒸散发计算模块。在ＭＩＫＥＳＨＥ蒸散发模块的计算中，需要运用气象和植被覆盖数据对降雨蒸发量进行预测。其中包括植物冠层对降雨的截留量、由植物冠层降落到土壤层的水量以及植物冠层、土壤表面和根系的蒸发量等决定因素。

２）坡面流模块。当流域地表的净雨强度超过土壤的下渗能力时，就会产生地表径流，地表径流沿坡面汇入河网，最后由河网内的流域出口断面流出。而径流在河网内的流经路径及径流量主要由流域地形、地表的摩阻力以及通过蒸散发和下渗失去的水量决定。在ＭＩＫＥＳＨＥ的坡面流模块中，主要采用扩散波对圣维南方程组进行近似模拟，运用有限差分的方法对方程组进行计算。

３）不饱和带模块。不饱和带是流域水文系统中一个至关重要的部分，同时也是ＭＩＫＥＳＨＥ的一个核心计算模块，在模型模拟应用中起着重要作用。受重力作用影响，不饱和带出流在下渗过程中主要作垂直方向上的运动。ＭＩＫＥＳＨＥ计算模块包含了反复迭代的计算过程，用于校正土壤含水量以及计算土壤浅层地下水位的动态变化。

ＭＩＫＥＳＨＥ在计算不饱和带模块时主要采用３种方法：Ｒｉｃｈａｒｄｓ方程、模拟简单的重力水出流过程、两层水量平衡模型。

４）饱和带计算模块。在ＭＩＫＥＳＨＥ的水流运动模块中，饱和带模块用于计算流域中的饱和地下水流运动。ＭＩＫＥＳＨＥ可以描述三维水流在异质含水层不定边界上的运动。模型对于时空变换的独立参数，采用三维达西方程进行模拟，并用隐式有限差分进行数值求解。

ＭＩＫＥＳＨＥ对上述地下水模拟提供兩种解决办法，一种是逐次超松弛法，另一种是前承条件共轭梯度法。

５）河流计算模块。模型中的不同水量最后通过河网进行汇流，最后到达流域出口。坡面流、壤中流和基流通过格网向流域较低的方向传播，最后到达流域相邻格网的最低点，也就是河流的位置，然后汇入河流，进行下一步的河流汇流计算。由于ＭＩＫＥＳＨＥ本身不包含河道汇流的计算，因此需要耦合ＤＨＩ公司研发的另外一个水动力演算模型ＭＩＫＥ１１进行这一部分的计算。耦合了ＭＩＫＥ１１水动力模型之后，能够对地表和地下的水流动力过程进行完整的、动态的耦合模拟。

ＭＩＫＥ１１河流模型采用的是水动力学模型，即明渠不稳定流隐格式有限差分解，其差分格式采用六点中心隐式差分格式，其数值计算采用传统的“追赶法”。ＭＩＫＥ１１计算参数包括两类：数值参数，主要是方程组迭代求解时的有关参数，如迭代次数及迭代计算精度；物理参数，主要是河网的阻力参数。

４．２模型构建

１）地表高程。本文使用的地表高程数据来源于地理空间数据云提供的ＳＲＴＭＤＥＭＵＴＭ３０ｍ分辨率数字高程数据产品。

２）气象数据。搜集沈乌灌域附近气象站数据，用Ｐｅｎｍａｎ－Ｍｏｎｔｅｉｔｈ公式计算研究区域的参考蒸散发量。

３）土地利用。通过遥感影像解译，获取区域及各分干渠范围内耕地、水域、草地、居民区、道路、裸地等土地利用类型。土地利用中不同下垫面对应的植被叶面积指数、根系深度以及蒸散发参数都是不同的。叶面积指数的获取在没有实测数据的前提下，借鉴已有研究中ＮＤＶＩ与叶面积指数之间的经验关系。

４）农田灌溉。灌溉模块是对灌域从引水到田间灌水过程的模拟，引水经过渠道逐级损失进入田间，一部分形成下渗补给土壤水和地下水，一部分消耗于蒸散发。由于每年种植结构以及总引水量受人为影响大，总引水在各分干渠及主要支渠间的配水过程非常复杂，因此没有客观规律可循。为了尽可能真实反映灌溉过程中的渠系渗漏、田间灌溉和渗漏，在渠道引水测流数据的基础上，结合渠道渗漏试验和田间试验确定的渠系水利用系数、田间水利用系数，分渠段计算得出相关渠段控制范围内的渠系渗漏水量和田间灌溉水量。

５）农田排水沟。模型中河流模块用于模拟排水沟的排水过程。排水沟位置确定后，通过实测的沟底高程以及末端测流断面特征和水位流量关系，结合排水沟比降，每２ｋｍ设置一个计算断面，以满足计算精度的需要。最终通过河流模块与ＭＩＫＥＳＨＥ地表、土壤、地下水模块之间的耦合完成设置。

６）土壤数据及参数。土壤数据主要来源于ＦＡＯ网站上的ＨＷＳＤ数据库。数据库中土壤分为上层（０～３０ｃｍ）和下层（３０～１００ｃｍ）。运用ＳＰＡＷ软件的ＳＷＣ模块，结合研究区土壤栅格属性表，按对应的属性寻找出变量。根据属性表中查出的变量，运用ＳＰＡＷ软件计算得出ＭＩＫＥＳＨＥ模型需要的土壤参数。

４．３模型参数率定及验证

２０１５年为率定期，２０１６年为验证期。通过调整地下水模型结构与含水层参数，利用不同时期地下水实测数据及二维流场，对模拟水位和流场进行校验。通过修正模型结构与参数，控制并再现地下水径流条件，使地下水达到了较好拟合。由拟合结果（图略）可以看出，模拟结果与实测数据的拟合情况较好。对已有的４４眼地下水观测井拟合结果进行统计，见表１。由表１可以看出，平均误差在０．５ｍ以内的井数为３４眼，占比７７．２７％。综上所述，模型精度满足要求，可以用于地下水埋深变化预测。

５地下水埋深预测

通过模型模拟，对沈乌灌域在３．０６亿ｍ３总引水量下的地下水埋深进行了预测分析。２０２５年沈乌灌域地下水埋深预测结果见图６。由图６可以看出，２０２５年地下水埋深在３．４ｍ以上的区域明显扩大，由２０２０年的４．８６万ｈｍ２增加到了９．３４万ｈｍ２。由图７可以看出，到２０３０年，沈乌灌域地下水埋深持续增加，大于３．４ｍ的区域持续扩大，仅在东南部、北部和西部局部范围内存在小于３．４ｍ的区域。大于３．４ｍ的区域由２０２５年的９．３４万ｈｍ２持续增加到了２０３０年的１４．５６万ｈｍ２，占整个灌域面积（１８．６８万ｈｍ２）的７８％。

６结果分析讨论

通过上述分析可知，在保持３．０６亿ｍ３引水量的情况下，沈乌灌域的地下水位会持续下降。预测结果显示，到２０２５年将有５０％的区域地下水埋深超过３．４ｍ，到２０３０年将有７８％的区域地下水埋深超过３．４ｍ。在地下水埋深超过３．４ｍ的区域，天然植被状况将会受到明显影响，主要体现在植被覆盖度下降。

７结论

基于ＭＩＫＥＳＨＥ分布式水文模型，采用实测水文数据和基于遥感的ＮＤＶＩ数据，运用数值模拟和统计分析等手段，对沈乌灌域出让水权后的地下水埋深和ＮＤＶＩ进行了预测分析，得到的结论如下。

当地下水埋深小于３．４ｍ时，ＮＤＶＩ与土壤质地存在一定的相关性；当地下水埋深大于３．４ｍ时，地下水埋深与ＮＤＶＩ具有显著负相关性，在其他因素不变的条件下，地下水埋深每增加１ｍ，土质基本单元的ＮＤＶＩ平均值会降低０．０８０７。

根据模型预测结果，相较于２０２０年，２０２５年地下水埋深在３．４ｍ以上的区域将增加到９．３４万ｈｍ２；到２０３０年，沈乌灌域地下水埋深将持续增加，大于３．４ｍ的区域将增加到１４．５６万ｈｍ２，仅在东南部（引黄总渠首）、北部（地下水排泄区）和西部（山洪补给区）局部范围内存在小于３．４ｍ的区域。针对天然植被可能存在的状况，灌域管理单位要加强灌域历年凌汛期以及灌溉间隙的生态补水调度，针对地下水埋深增加区域的补水水体进行补水级别划分，以减缓重点区域地下水埋深增加趋势。

农业作为黄河流域的用水大户，随着黄河流域水资源短缺与经济社会和生态环境发展矛盾的加剧，水权交易作为推动区域水资源优化配置的有力措施，需要在注重经济社会效益的同时，更加重视生态环境的响应，全面、系统地做好水权交易引起的水资源重新配置后的外部效应将是今后一个重要研究方向。

本文的不足之处在于地下水观测数据序列不长，对灌域内天然植被种类的变化考虑不足。下一步将在灌域地下水埋深变化频繁的区域内设置典型区域生态样方，针对天然植被覆盖度变化以及植被种类更替规律等进行系统监测研究。