水稻灌区基于多过程耦合的分布式水转化模型

王 璞，卢 炤，霍再林

（中国农业大学 水利与土木工程学院，北京 100083）

1 研究背景

水稻种植面积占我国粮食作物种植总面积的27%，稻田灌溉用水量占我国农业总用水量65%以上，但其灌溉水利用系数仅在0.5 到0.6 之间［1-3］。由于灌溉定额大，水稻灌区灌溉水垂向及水平向转化频繁，灌区水转化模型是灌区水资源科学管理的基础［4-6］。

目前针对灌区水转化模拟的方法主要包括：（1）基于陆地水文模型的灌区水转化过程模拟。SWAT 模型已被广泛应用于灌区水转化、面源污染、农业生产力等方面的模拟［7-9］。特别是近年来，SWAT 模型与MODFLOW 模型结合，实现了灌区地表水-土壤水-地下水的耦合模拟［7］。（2）基于农田水文模型的分布式模拟。一般是对农田水文模型进行率定验证后，将其拓展到灌区尺度进行模拟［10-11］。该方法可对农田垂向水转化及作物生长过程进行详细描述，但由于缺乏对地下水水平运动的模拟以及灌溉-排水-地下水等水转化过程的表征，无法应用于不同情景下灌区水转化的模拟。（3）基于灌区分布式水平衡的模拟。一般是对灌域水平衡进行分析计算，进而将其拓展至灌区尺度。如岳卫峰等［12］针对内蒙古河套灌区义长灌域基于四水转化关系进行了灌域水均衡分析，该方法尽管可获得灌区整体水平衡，但其时空分辨率低，无法实现灌区尺度详细的水转化过程模拟［12-13］。

综上所述，现有的灌区水转化模型尚未充分体现灌区在灌溉与排水驱动下，灌溉水、土壤水及地下水在垂向与水平向频繁转化的特点。基于此，本研究旨在发展灌区尺度基于灌排-农田水文-地下水过程耦合的水转化模型，并将其应用于黑龙江省绥化市庆安县和平灌区，为灌区尺度水平衡分析提供有效方法。

2 研究区域及水文监测

2.1 研究区概况本研究以黑龙江省和平灌区为研究区，该灌区地处松嫩平原和小兴安岭余脉的交汇地带，属呼兰河流域中上游（图1）。和平灌区地理坐标东经127°21′至127°45′，北纬46°5′至46°92′；属寒温带大陆性季风气候，年均日照时数2599 h，年均气温1.69 ℃，平均年降雨量577 mm。和平灌区控制面积约170 km2，其中耕地占比90.5%，耕地内97%种植作物为水稻。灌区分为井灌区与渠灌区，其中渠灌区地表引水来源主要包括呼兰河、拉林清河与安邦河。灌区内现有干渠1 条，主要支渠20 条，主要排水沟道9 条。灌区和平渠首灌溉期多年平均灌溉引水量5.74 亿m3。灌区地下水埋深由上游的2 ～ 3 m 向下游井灌区13 ～ 18 m 逐渐增大。

图1 研究区位置

2.2 灌区水文监测与数据收集灌区水文监测数据包括2018—2019年河流水位、干支渠流量、沟道流量及地下水位，用于模型的率定与验证。河流水位监测断面包括李山屯渠首、郑文举渠首、安邦河渠首，于河流水位监测断面布设水尺并标定黄海高程；渠道及排水沟道数据除干渠流量外监测灌区上、中、下游3 条典型支渠（二支渠、六支渠、庆胜支渠）渠首流量，2 条典型排水沟道（一支排、七支排）流量，每日监测一次。渠道流量数据使用水尺监测数据结合水位流量关系曲线计算，沟道流量数据采用多普勒超声流量计测量。由于灌区内支渠众多且未实现数字化监测，本研究难以监测全部渠道每日流量。在对灌区调研后认为灌区上中下游对应的支渠流量变化趋势基本一致，研究区域内其他支渠的流量使用已知流量的3 条支渠结合渠道设计流量推求。为保证一定的准确性使用逐日观测的干渠流量进行对比验证，如支渠总流量与干渠流量不符则考虑闸门开闭、各渠道控制面积等因素加以调整。此外，研究区布设13 眼地下水观测井，采用Hobo U20-001-01 自计水位计逐日6 次测定各观测井地下水位（图2）。模型运行所需气象数据收集于毗邻灌区的铁力气象站。同时基于灌区内土壤质地调查（64 处土壤剖面），获取作物根区土壤条件及水分运动参数分布，用于灌区土壤水分运动模块（图3）。

图2 研究区DEM 数据及水尺、观测井布设位置示意

图3 研究区土壤参数空间分布

3 灌区分布式水转化模型

3.1 模型框架

3.1.1 模型结构 该模型包括灌溉、排水、农田作物生长与水转化、地下水运动4 个模块，基于Mat⁃lab R2018b 与ArcGIS 10.7 软件实现灌区供耗排耦合过程定量模拟。其中灌溉与排水模块是在对灌区灌溉渠道及排水沟道进行数字化的基础上，进行灌溉渠道入渗损失、各水文响应单元灌溉水量、井灌区抽水量以及排水单元排水量的计算；农田作物生长及水循环模块旨在耦合模拟农田一维水转化及作物生长过程；地下水运动模块主要模拟计算地下水的垂向及水平向运动。模型将水文响应单元划分为均一网格，灌溉模块将渠道来水分配至水文响应单元，农田作物生长与水转化模块对各水文响应单元作物生长及水分转化过程进行迭代计算，输出作物生长指标及水分分量运算结果，为排水模块及地下水模块提供地表产流量及地下水补给量。各模块进行时空耦合，最终实现灌区供耗排过程的耦合模拟。模型具体结构及关系见图4。

图4 模型整体结构图

3.1.2 边界条件 依据研究区水文地质条件及地表河流分布，西北以呼兰河为定水头边界；东北以安邦河为定水头边界；南部以柯木克河为定水头边界；西部为依据呼兰河与柯木克河插值确定的定水头边界；东南部为定流量边界（图5），边界条件中水头及流量均采用实测值。

图5 研究区地下水隔水顶板高程及边界条件分布示意

3.1.3 水文响应单元划分 基于ArcGIS 10.7 软件对研究区进行水文响应单元划分，水文响应单元大小为250 m×250 m，灌区内有效水文响应单元共2715 个。

3.2 各过程定量表征方法

3.2.1 农田水转化与作物生长过程 模型中作物生长过程采用EPIC（Erosion Productivity Impact Calcu⁃lator）模型进行描述［14-15］。为便于准确量化稻田灌溉-土壤水-地下水转化过程，模型将稻田在垂向分为5 层，由上而下分别为储水区、泥浆层、传输层、饱和层与不透水层（图6）。储水区为稻田水层，当储水区高度超过田埂高度时会产生田间排水；泥浆层为犁底层之上的土壤；传输层为犁底层之下，地下水位之上的土壤，一般透水性强；饱和层为地下水所在，饱和层与传输层相邻界面为地下水位；不透水层为隔水底板。

图6 土壤垂向分层示意

水稻生育期内，稻田土壤水与田面水层逐日均有补给与消耗。在稻田存在水层与无水层时，稻田水量平衡总方程分别可表示为：

稻田存在水层：

稻田不存在水层：

式中：S 为田面水层深度，mm；P 为降雨量，mm；I 为灌溉量，mm；ET 为蒸散发量，mm；D 为排水量，mm；J 为深层渗漏量，mm；m 为土壤含水量，mm；t 为日序数。

依据作物系数与水分胁迫系数计算农田作物实际耗水量，其中参考作物耗水量ET0的计算采用FAO-56 推荐的Penman-Monteith 公式，作物潜在耗水量ETp依据参考作物耗水量与作物系数计算得出，同时依据植物叶冠生长程度划分土壤潜在蒸发Ep与作物潜在蒸腾Tp，分别计算水分胁迫下的实际蒸发、蒸腾量Ea与Ta，最终获得农田实际蒸散发量［16-19］。在此过程中所需的作物生长指标如叶面积指数等由EPIC 模型计算提供。稻田深层渗漏量依据田面水层深度及土壤质地等值计算而得。

当田面存在水层时，认为稻田深层渗漏由田面水层深度与土壤质地决定，其经验公式［20］可表示为：

式中a、b 为取决于稻田土质的参数。

当田面不存在水层且稻田土壤含水率大于田间持水率时，稻田深层渗漏依据式（4）计算［21］：

式中： ms为饱和含水率； md 为残余含水率； ks为饱和导水率；C 为常数； RDmx为作物最大根深，mm； RDt为作物当日根深，mm； mgt为泥浆层当日含水率。

模拟时将泥浆层与储水区视为整体进行水量分配。同时依据气象、土壤质地参数等数据计算逐日传输层向饱和层的渗漏量以及储水区与泥浆层蒸散发量。基于逐日的水分分量与灌溉制度迭代计算地表径流量与灌溉量，以更新下一日模型水分分量。计算流程见图7。

图7 农田水转化模拟流程图

3.2.2 地下水运动过程 各水文响应单元逐日地下水位变化量包括该水文响应单元垂直方向与地上部分的水量交换和与相邻水文响应单元的水平径流交换。地下水位的垂直补给与消耗一般取决于河流及渠道补给、深层渗漏等，对于井灌区地下水开采量为地下水重要的汇项。同时研究区地下水水力梯度较大，水平流动明显。根据裘布依假设与达西定律，地下水水平径流量可表示为［22］：

式中：JL为目标水文响应单元相邻四个方向任意一个单元向该目标单元输入的逐日单位面积地下水径流量，m/d；K 为目标单元的渗透系数，m/d；h 为目标单元潜水含水层厚度，m；B 为目标单元边长，m；Lx为目标单元周围四个方向上任意单元的地下水位，m；L 为中部目标单元的地下水位，m；X 为中部目标单元中心点与周围相邻单元中心点的距离，m。

3.2.3 灌溉输配水过程 在ArcGIS 10.7 软件中将灌区主要渠道进行数字化处理，划分渠道控制区域（图8）。根据水稻灌溉的实际情况，按照灌溉制度确定田间灌溉水量，并基于渠道来水量对各渠道控制的水文响应单元进行灌溉水量分配。同时利用经验公式估算渠道输水渗漏损失，渠道渗漏量公式如下［23］：

图8 研究区进水渠道、排水沟道数字化示意

式中：σ为每千米渠道输水损失（以渠道净流量百分数计）；A 和m 分别为渠床土壤透水系数、渠床土壤透水指数，二者依据渠道所在区域的土壤质地确定；Qn为渠道净流量，m3/s；Q1为渠道渗漏流量，m3/s；L 为渠道长度，km。

依据公式计算渠道渗漏的同时，根据渠道灌溉水量与渗漏损失量，实时更新渠道净流量。

3.2.4 排水过程 在ArcGIS 10.7 软件中将灌区主要排水沟道进行数字化处理，划分沟道控制范围（图8）。研究区为水稻灌区，生育期水量充沛，渠道来水很大一部分在补给农田灌溉后仍有剩余，故灌区排水水量由两部分组成：一部分为田间水层高度超过田埂后产生的田间排水，另一部分为渠道剩余的水量产生的退水。田间排水由各个水文响应单元稻田水平衡模型计算而得，退水产生的排水由水文响应单元所属支渠灌溉之后的剩余水量组成。

4 模型率定与验证

4.1 作物生长过程分别采用2018年与2019年灌区定位监测水稻叶面积指数（LAI）与株高数据进行作物生长模型部分参数（表1）率定与验证（图9）。选取R2与RMSE 作为模拟结果评价指标，结果表明，EPIC 模型可以较好地模拟研究区水稻生长过程（表2）。 EPIC 模型基于积温模拟作物生长，因此区域上的积温不确定性会导致定位监测数据与模拟结果存在一定偏差。如2019年叶面积指数模拟最大值出现时间先于实测值，且株高模拟值大多高于实测值应是由积温误差导致。

图9 2018、2019年作物生长指标模拟结果与实测值对比

表1 农田水转化模型参数初始值及率定值

表2 农田水转化模型率定及验证效果评价指标

4.2 灌区蒸散发选取SEBAL 模型的灌区蒸散发遥感反演值进行区域蒸散发部分参数的率定验证，将蒸散发遥感反演值与模型模拟结果二者的逐日区域平均值进行对比，率定蒸散发模块参数（田间持水率、最大作物系数、蒸腾胁迫曲线影响因子等）。2018、2019年均为丰水年，降雨较多，由于云层覆盖等原因造成遥感数据缺失严重，故选取2018 与2019年两年生育期内影像较为完整时段加以对比（图10），分别为2018年83 d，2019年65 d。使用2018年数据进行模型率定R2=0.78，RMSE=0.98 mm/d；2019年数据进行模型验证R2=0.75，RMSE=0.86 mm/d。认为模拟结果符合精度要求。

图10 2018、2019年区域蒸散发模型模拟值与遥感反演值对比

4.3 地下水埋深采用2018年和2019年地下水埋深数据用于地下水运动模块参数（导水率、给水度）率定验证（表3），取2018年地下水埋深数据用于模型率定，2019年地下水埋深数据用于模型验证。进行地下水运动模块模拟时，参考黑龙江水利勘测设计院编写的庆安县1∶5 万水文地质勘察报告（1987年）与黑龙江水文地质图（1∶350 万）对研究区域断裂情况的描述，将研究区划分为7 个水文地质分区进行模拟（图11）。模拟结果如图12，可见各点均匀分布在1∶1 相关线两侧附近，可认为模型无系统性错误。由表4 可见各种模型评价指标均在允许范围内，地下水埋深模拟值与实测值吻合较好，模型运行稳定、可靠，可用于灌区供耗排耦合过程定量模拟。2018、2019年模拟结果评价指标在灌区范围内呈现出下游优于上游的情况，其原因在于研究区由上游向下游地下水埋深递增，埋深较浅的地下水位观测井数据更易受到偶然因素干扰，但埋深较深的地下水位观测井实际情况与模型假设较为一致，因此模拟结果较为理想。2018、2019 两年由于灌区生育期内持续的田间深层渗漏以及渠道入渗补给，灌区整体地下水位呈上升趋势，但井灌区由于存在抽取地下水用于灌溉的原因，其地下水位升高幅度低于渠灌区。

表4 地下水运动模块率定及验证效果评价指标

图12 率定期、验证期地下水埋深计算值与观测对比

表3 潜水给水度Sy 、渗透系数ks2 初始值与率定值

图11 地质分区示意

4.4 排水流量2018年对一支排及七支排沟排水流量数据进行监测，对两条排水沟排水流量模拟结果与实测结果进行了对比（图13），其中一支排沟R2=0.67，RMSE=0.38 m3/s，七支排沟R2=0.61，RMSE=0.15 m3/s。2019年由于仪器故障，缺测该年排水沟道流量数据。结果表明，模型能较好地模拟灌区排水过程，但峰值时误差较大，分析认为误差来源是将排水沟道控制范围内的田间排水与渠道退水等统一划归到该排水沟，未考虑沟渠联通等情况。

图13 2018年一支排、七支排流量模拟结果与实测值对比

5 讨论

通过作物生长指标、地下水位、蒸散发及排水流量等验证表明，本文所建立的水稻灌区基于多过程耦合的分布式水转化模型具有较高精度，可以客观反映水稻灌区复杂的水转化过程。相对于目前常用的灌区水转化模型来说，该模型针对水稻灌区灌排复杂及农田灌溉水与地下水转化频繁的特点，实现了灌溉过程、排水过程、农田水分运动及作物生长过程、地下水运动过程的动态耦合。尽管众多学者尝试将流域水文模型（如SWAT 模型）与地下水运动模型（如MODFLOW 模型）进行耦合用以模拟灌区水转化过程，但由于流域水文模型对于土壤垂向水分运动的概化，难以真实反应农田灌溉水入渗与蒸发过程［7］。特别对于稻田具有明显水分分层特性条件下，本研究中所提出的稻田水转化过程具有明显的优势。进一步来说，稻田分层水转化模型较传统基于动力学过程对土壤水分运动模型参数依赖性减小，对灌区分布式水转化模型具有较强的适用性。

此外，该模型将灌排渠沟进行数字化，突出了灌区灌排过程的时空差异性。尤其是对渠系输配水过程中渗漏损失及对地下水的补给这一灌区水转化过程特有的环节得以量化表征，一定程度提高了模型精度。相比较而言，目前基于农田水文模型的分布式模拟及灌区水平衡分布式模拟模型均对复杂渠沟空间分布及其输配水及排水过程缺乏必要的考虑［6，12-13］。从模型结构来看，该模型由灌溉模块、排水模块、农田作物生长与水转化模块以及地下水模块耦合而成，模型中农田水转化过程、地下水运动过程参数较现有基于动力学过程的模拟方法大大减少，特别是模型所采取的土壤水与地下水过程统一的网格化水文响应单元提高了灌区尺度分布式模型的模拟效率。

需要指出的是，本文所建立的水稻灌区分布式水转化模型依赖于灌区渠沟分布、高分辨率土地利用及作物种植分布等空间数据。同时，相对于灌区水转化客观物理过程，模型对于排水过程模拟尚未考虑通过田埂侧向渗流与沟渠联通情况，一定程度引起了排水模拟的误差；模型中灌溉模块所涉及的灌溉制度需进一步基于渠道来水及作物需水进行优化设计。

6 结论

本文针对水稻灌区水转化特点，以黑龙江省和平灌区为例，建立了灌区基于多过程耦合的分布式水转化模型。基于灌区实测水文过程对模型加以检验，结果表明模型可以较好模拟灌区主要水转化过程。与传统灌区水转化模型相比，该模型实现了灌溉过程、排水过程、农田水分运动及作物生长过程、地下水运动过程的动态耦合，可以有效表征灌区多尺度复杂水转化过程。特别在农田尺度，通过将稻田水转化模型与作物生长模型耦合，克服了传统灌区水转化模型对作物生长过程模拟不足的劣势。模型中将灌排渠沟进行数字化，突出了灌区灌排过程的时空差异性，实现了灌区水转化的分布式模拟。该模型对于灌区用水效率时空动态评估、灌溉水资源管理等具有重要意义。后续研究中将对模型灌溉模块和排水模块优化提升，进一步提升模型的合理性及普适性。