基于Modflow-Hydrus耦合模型的改良盐碱土水平井排水效果分析

代锋刚，王晓燕,谷明旭，余婉露，连鹏达

(1. 河北地质大学河北省水资源可持续利用与开发重点实验室, 河北 石家庄 050031； 2. 河北地质大学河北省水资源可持续利用与产业结构优化协同创新中心， 河北 石家庄 050031; 3. 河北省高校生态环境地质应用技术研发中心， 河北 石家庄 050031; 4. 中国地质调查局水文地质环境地质调查中心，河北 保定071051； 5. 河北省地质环境监测院河北省地质资源环境监测与保护重点实验室, 河北 石家庄050021)

利用盐碱地的治理改良来增加补充耕地资源，是国家从战略高度坚守耕地红线和保障粮食安全的重大举措之一[1].盐碱地改良技术方法有物理、化学、生物、工程技术等，其中化学法见效快，但环境污染的风险大，成本高.物理法、生物法用客土、覆盖等技术存在返盐风险.从盐碱地的综合整治方法、排盐效果评价、投资成本、风险分析、推广使用等方面比较，工程排水洗盐是目前重要的水利措施之一.国内外不少学者对盐碱地治理改良进行了研究，日趋成熟的暗管排水排盐技术在盐碱地治理改良中得到较为广泛的应用[2].国外在通过控制排水降低非点源污染[3]、控制排水智能设备、提高水资源利用效率等方面成果显著[4].国内学者盐碱土治理改良方面，排水系统结合区域(新疆[5]、滨海地区[6]等)水文气象、土壤类型、作物种类进行综合调控，为排水技术的推广打好了一定的基础.其中水平井排水技术与排水沟、排水暗管、井灌井排等方法相比具有较大优势[7-8],其机理是灌溉与降水淋滤洗盐作用，利用水平井能有效控制土壤毛细水上升高度和地下水位的临界深度，避免改良的盐碱土壤再次返盐.

文中通过现场、室内试验及数值模拟技术，建立饱和-非饱和土壤水Modflow-Hydrus耦合模型，对不同水文气象环境下水平井的排水效果进行分析，为滨海盐碱地改良、排水工程设计提供理论参考.

1 材料与方法

1.1 试验区概况

试验区设在冀中平原河北省沧州市沧县境内典型盐碱地区，为冲积平原向滨海平原过渡带，地势总体趋势是由西南向东北微倾斜，海拔5～7 m，自然坡降约0.000 08.区内多年平均降水量约600 mm，6—9月高水位期的降水量约占全年降水总量的78%，降水量年内分布极为不均，年平均蒸发量971 mm.根据试验区钻孔岩心编录表，现场采集不同层土样进行颗粒分析、容重测试等室内试验，获得不同层土壤的容重ρ、孔隙度Φ及颗粒级配数据，试验区土壤物理参数见表1，表中d为土层深度，θ为粒组含量.

表1 试验区土壤物理参数Tab.1 Soil physical parameters in test area

1.2 试验设计

试验设计如图1所示.

图1 试验设计示意图(单位：m)Fig.1 Schematic diagram of test design(unit: m)

试验区为300 m×200 m的矩形场地，东西设有2条排水沟，距离280 m,其中间设有3条水平井，间距110 m,水平井布设与排水沟平行，水平井井管和滤水管长度分别为180 m和100 m，内径210 mm，地下埋深9m，所在位置为承压含水层，区内有东西2条排水沟，3眼监测孔和3条水平井分别位于中间.

1.3 土壤水分参数测定

1) 颗粒分析

颗粒分析是指测定田间土壤中不同粒径颗粒质量与该土总质量的百分比.野外土样采集并编号封存，进行土壤颗粒分析，粒径小于0.075 mm时采用密度计法，粒径大于0.075 mm时采用筛分法，两者相结合，土样定名主要依据土工试验规范和岩土工程勘察规范.

2) 土壤容重

环刀法测定土壤容重，即根据环刀内湿土质量与干土质量的比值计算所测土壤含水率，由于实测值具有一定精度，因此作为典型参考值用于率定岩土水分特征参数.

3) 水文地质参数

利用单孔抽水试验来获取试验区含水层的导水系数和释水系数等水文地质参数.本次抽水试验目标含水层为承压含水层，岩性以粉砂为主，地下水埋深8～30 m.抽水孔井型结构为完整井，埋深8～30 m为滤水管,其余为不透水管.本次选用Diver水位监测仪，监测频次高，数据密度大，用于观察记录地下水位的动态变化，记录频次为1次/min.抽到4 h基本稳定，继续稳定16 h后停抽并开始水位恢复，得到可靠的观测资料，利用水位恢复试验数据推求参数.

2 Modflow-Hydrus耦合模型

2.1 Modflow-Hydrus耦合模型简介

Hydrus模型通常仅局限于小范围包气带土壤水计算问题，不建议用于大区域范围地下水计算问题.虽然模拟软件Modflow在区域尺度地下水计算方面具有优越性，但是由于潜水面不断变化，Modflow处理其边界条件不够准确，如果直接利用Modflow计算分析地下水浅埋区的水文地质问题，那么会产生较大的误差.故采用Modflow-Hydrus耦合模型，各取所长，弥补两者不足.首先，对Hydrus模型计算结果进行统计分析，确定综合补给强度-埋深分段函数.其次，估算试验区不同时段不同区域的综合补给强度.最后，用BASIC语言将其编成Modflow软件面状补给包的数据文件格式Rch.dat，导入Modflow软件进行计算地下水位埋深等值线，重复利用综合补给强度-埋深函数计算地下水补给强度，直到模型收敛.计算得到rch1.dat和rch0.dat分别乘以一个权重系数，得到新源数据文件导入Modflow进行模拟计算，能够加快收敛速度.通过重复迭代计算，Modflow-Hydrus耦合模型体现Hydrus在描述降水入渗-潜水蒸发自适应源汇项时的高精度及Modflow具备大区域地下水流计算问题的优点，同时降低了地下水变源汇项确定的难度.

Modflow(2005)UZF1是 Visual Modflow三维地下水模拟软件系统中新增模块，其中地下水补给过程和蒸散发模块可以用于模拟非饱和带土壤水运动问题.UZF1模块替换了传统Modflow中蒸发蒸腾模块和面状补给模块.UZF1模块采用非饱和带土壤的入渗速率，并不是饱和渗透系数，饱和渗透系数与入渗速率紧密相关.UZF1模块与蒸发蒸腾模块不同，UZF1模块中蒸散发量计算时，起初考虑计算蒸散发极限埋深上部的包气带，忽略地下水极限埋深，如果水量计算不能满足蒸散发量，直接从地下水中扣除.另外，UZF1模块与蒸散发模块，如果水位高于地面，地下水直接流向地表.UZF模块与传统Visual Modflow相比，在面状补给的处理方面，能更精细地刻画非饱和带水文运动过程，但是UZF模块仅考虑了重力势而忽略了毛管势，因此对非饱和带土壤水分运移过程的呈现与实际还存在偏差.

与两者相比，Modflow-Hydrus耦合模型能实时计算包气带和饱水带的水分交换量，降水入渗或灌溉入渗补给量、潜水蒸发量具有自适应性，因此耦合模型提高了非饱和带和饱和带上边界源汇项的计算精度，增加了数值模拟的仿真度[9].

2.2 数学模型

根据试验设计情况，将整个入渗过程概化为以水头为变量一维变饱和带水分运移问题，其数学模型可以用Richards方程[10]描述，即

(1)

式中：h为压力水头，cm，非饱和带内为负值，饱和带内为正值；K为非饱和带水力传导度，cm/d，随含水率变化；C为容水度，cm-1，相当于饱和含水层的储水率；s为源汇项，cm/d.

初始条件和边界条件为

(2)

式中：h0为模型剖面初始水位，cm；上边界条件为面状补给或排泄，包括大气降水、潜水蒸发、灌溉排水等，qs为地表土壤水流量,cm/d；下边界条件为地下水潜水面，d为地下水位埋深，m.

Hydrus模型剖面包括整个非饱和带范围，模拟计算不同埋深非饱和带水分运移.地面以下30 m为剖面模型范围，模型以5种不同岩性概化为5层，模型垂向按2 cm间距进行剖分.

模拟应力期设为10a(3 600d)，采用变化的时间步长，软件根据数值解收敛迭代次数自动调整时间步长.设初始时间步长为0.001 d，步长区间设为0.000 1 d～0.200 0 d.

据试验区水文地质条件，地表以下30m内含水层概化为2层.上层厚度约4 m，岩性以粉土为主，夹薄层粉黏、粉砂层，概化为潜水含水层；中间夹层为埋深5～8 m黏土，视为弱透水层；下层厚度22 m，岩性为粉砂，概化为承压含水层；底部为黏土层，构成相对隔水底板.模型的范围以试验区为中心向四周外扩大1.2 km，通过模型试算，在5 m抽水降深时，3个水平井水位降幅小于1 cm，忽略不计，视为无穷远边界，其余边界均定为隔水边界.源汇项根据试验区水文气象环境设定.根据上述水文地质条件，计算区地下水数学模型可描述为潜水含水层地下水微分方程和下部承压含水层地下水流微分方程，即

(3)

(4)

上述式中：H1为潜水水位，m；H2为承压水水位，m；T1为潜水含水层导水系数,m2/d；T2为承压水含水层导水系数,m2/d；μ为潜水给水度；μ*为承压水弹性释水系数；Г1，Г2为潜水与承压水含水层边界；Q为承压水开采量，m3/d；σ′为潜水与承压水之间的垂向越流系数，d-1；w为源汇项，降水-气象环境及灌溉环境下的综合补给强度,m/d.采用Visual Modflow系统对水平井工程排水效果进行分析模拟，有效计算范围为矩形，长为2 650 m，宽为2 500 m，均匀剖分网格，分别用50 m×50 m方格空间离散，试验区采用加密剖分，间距10 m，其邻近部分用25 m格距过渡.

2.3 土壤水分特征参数

根据现场原位测试和室内试验的数据，分析确定土壤水分参数，包括饱和渗透系数、颗粒分析、干容重、孔隙度.利用SSCBD模型可以对水分特征参数进行预测，同时根据各层土壤的粒组含量和容重确定土壤的水分特征参数(田间持水率θr，饱和持水率θs)，根据室内试验及野外试验获取的数据进行修正.结合土壤水分参数经验值，构建了土壤水分参数数据库(见表2).

表2 土壤水分特征参数Tab.2 Soil moisture characteristic parameters

对耦合模型进行识别和检验，计算地下水位等值线与实际观测地下水等值线的年内变化趋势进行分析，用地下水位的拟合均方差小于允许误差作为数值解收敛的判别标准.地下水流场拟合精度以模拟区地下水最高水位和最低水位差作为参照，拟合误差δ≤10%，水位模拟动态和实际动态基本相似，水位变幅误差Δ≤10%，因此模型基本可靠.

3 结果与分析

3.1 模拟情景设置

试验区设置2种情景进行模拟，天然气象情景和典型灌溉情景，不同情景分持续排水和断续排水，根据不同年排水量(Ma)、不同工程管理方式等条件分析水平井有效控制范围(单水平井A1为SA1，多井并排A2为SA2)和水平井排水效果.

天然气象情景下，水平排水井的排水强度设为0.40、持续排水0.30、断续排水0.75，单位为m3/(d·m).不同排水强度的单井控制范围各异，水平井两侧地下水位降幅大于2 m(见表3)，说明天然气象环境下，不同的年排水量单排水井的控制范围不同，单井控制范围和年排水量呈正相关，对比持续排水，断续排水的控制范围显著减小，本次模拟的排水制度设定为每年3月、5月、10月进行排水.

典型灌溉情景下，水平排水井排水强度设为0.54、持续排水0.45、断续排水0.90，单位为m3/(d·m).不同排水强度情况下单井控制范围不同，水平井两侧地下水位降幅达2 m以上(见表3)，对比天然气象环境，单排水平井控制范围呈减小趋势.因此，不同的灌溉排水制度、田间工程排水方式、天然气象环境等对单水平排水井的控制范围均有影响，本次分别模拟以上不同情景水平井排水效果和地下水位等值线.

表3 各种环境下不同年排水量单井控制范围Tab.3 Single well control range of displacement in various environments

3.2 持续排水条件下模拟分析

持续排水条件下，对天然气象条件、典型灌溉的排水过程分别进行模拟分析，天然气象情景和典型灌溉条件下，预测分析试验区水位疏干至稳定后地下水位变化过程，计算平水期地下水位埋深等值线，如图2所示.

图2 试验区持续排水潜水埋深等值线Fig.2 Continuous drainage diving buried depth equivalent line in test area

图2a表明在排水强度为0.3 m3/(d·m)持续排水条件下，3个月后在排水井周围300 m以内，潜水含水层地下水位埋深有所增加，从现状1.25 m增至2.00 m以上，区内土壤水盐状况改善明显.在第3—4年，抽水基本达稳定时，从图2b和观测孔数据看出，排水强度为0.4 m3/(d·m)时，排水井周围400 m处平均水位降至0.75 m，即地下水位平均埋深达2.00 m以上，可以明显改善土壤水盐状况，单井有效控制范围大于800 m；排水沟两侧500 m处地下水水位从初始值降至0.45 m.

由流场叠加原理可知，若平行布置排水的水平井，水位降深叠加后大于0.8 m，水平排水井两侧有效控制范围相对较大，单井有效控制的范围可达到1 000 m.试验区的控制范围可以延伸外围300 m，由此可推算，1眼排水水平井就可以基本满足盐渍化改良的要求.由模拟计算结果可知，当地下水位埋深较小时，降水季节变化严重影响降水入渗补给量；当地下水位埋深较大时，随包气带厚度增加，降水入渗补给强度与降水季节密切相关，补给量峰值稍有延迟和滞后；当地下水位埋深较大时，在包括丰平枯的水文年内大气降水入渗补给强度大小和时间弱相关，年降水量平均分配在各个不同时间段.地下水位埋深较小时，包气带厚度较小，季节变化对地下水综合补给强度影响明显，若地下水位埋深达6 m以上，各季节时段地下水综合补给强度基本达到极值.即可认为当地下水位埋深较大时，通过厚层非饱和带降水入渗补给基本稳定，因此研究地下水流问题时常用入渗系数概化降水入渗过程.天然气象情景和典型灌溉情景2种环境下，在持续排水时试验区潜水位大幅度下降，形成以试验区为中心的地下水位降落叠加漏斗，明显疏干了潜水含水层，可以明显改善盐碱地的水文地质条件.持续排水条件下，水平排水井排水效果比较明显，有效控制范围比较大，但是田间工程管理不利于操作，建议盐碱土水利工程改良实践采用断续排水方式为宜.

3.3 断续排水条件下模拟分析

考虑田间灌排工程管理操控方便，针对天然气象条件、典型灌溉情景、断续排水过程分别进行模拟计算.在断续排水条件下，天然气象条件下和典型灌溉情景下，排水期间水平井有效控制范围内，潜水含水层地下水位下降明显，以试验区为中心形成地下水位降落漏斗.天然气象情景下，断续排水条件下，3个水平排水井不同排水时段计算的潜水位埋深等值线及承压水位埋深等值线，如图3所示，蓝色细虚线表示初始平均地下水位埋深线和排水疏干后潜水含水层地下水位埋深线，蓝色细实线表示承压含水层地下水位埋深线.

图3 试验区断续排水条件下不同排水期地下水水位埋深曲线Fig.3 Buried depth of groundwater level in different drainage periods under intermittent drainage conditions in the test area

图3表明第1年3月地下水位下降最大值为5.47 m，5月增至5.87 m，第2年5月增至5.92 m，天然条件下的排水时段，潜水含水层疏干效果明显，水文地质条件有利于盐碱土改良.统筹考虑水文气象情景及典型灌溉条件，探索田间水平井的科学合理排水模式，尤其是在蒸发强烈季节和容易返盐时段，有效控制潜水含水层地下水位，试验区因地下水位埋深小而导致返盐的现象可以减少.根据区域内主要农作物灌溉季节和返盐时段，设定每年排水时段为3月、5月、10月，其余月份不排水，通过模拟计算得到，疏干稳定后天然气象情景条件下，水平井排水强度为1.0 m3/(d·m)地下水位埋深等值线如图2c所示，典型灌溉条件下，排水强度为1.2 m3/(d·m)地下水位埋深等值线如图2d所示.天然情景下断续排水情况下，不同排水时段潜水含水层地下水位埋深等值线和承压含水层地下水位埋深等值线如图2c，2d所示，单水平井有效控制周围400 m范围，断续排水时段潜水含水层地下水位下降幅度较大，疏干效果显著.

天然气象环境下，在水平井排水前，人为干扰较少，试验区潜水蒸发和大气降水入渗补给影响不明显，数值基本相等，因此选择潜水蒸发量和降水入渗补给量作为指标量，在试验区排水水平井断续排水后，划分不同区域对潜水蒸发量与降水入渗补给量比值分别进行统计，然后进行试验区排水影响分区，如图4所示.

图4 试验区断续排水影响分区图Fig.4 Zone of intermittent drainage effect in test area

试验区排水后地下水位下降，影响毛细水上升高度，减小潜水蒸发量约占大气降水入渗补给量10%，改变土壤包气带水分运动形式，有效改善土壤的水盐状况；强烈影响区，地下水位下降幅度大，潜水蒸发占大气降水入渗补给量30%左右，其影响范围大小为4倍试验区面积；明显影响区，包气带土壤水毛细水上升高度受影响，因为其与潜水面有密切水力联系，随潜水面改变而变化，潜水蒸发量约占大气降水入渗补给量50%，其影响范围大小为3倍试验区面积.由此可推测试验区进行断续排水，以试验区为中心，面积为8倍试验区范围利于盐碱土改良的水文地质条件均可以得到有效改善，使非饱和带水分总体向下运移，浅部潜水含水层水位呈现下降趋势，甚至小于临界水位，有效控制毛细水上升高度，避免土壤返盐.

4 结 论

1) 通过野外和室内试验确定土壤水分特征参数，利用Hydrus软件建立试验区饱和-非饱和带土壤水分运移剖面数值模型.对试验区局部范围土壤水分运移过程进行分析，考虑试验区地层岩性、气象水文、灌溉排水条件，分析确定了不同时段地下水综合补给强度-埋深分段函数，可以借鉴到大区域盐渍土水利工程改良效果评价和大气降水入渗补给系数确定的研究.

2) 利用Modflow-Hydrus耦合模型分析不同情景下试验区水平井的地下水位控制范围及排水效果，在天然气象环境持续排水条件下，单个水平井能够控制800 m范围，断续排水条件下，单个水平井能够控制200 m范围；在典型灌溉持续排水条件下，单个水平井能够控制500 m范围，断续排水条件下，单个水平井能够控制100 m范围.平行布置水平井时，每条水平排水井的控制范围会叠加增大.单个水平井的控制范围影响因素较多，控制范围与水平井的年抽水强度大小呈正相关.

3) 水平井适宜渗透性小的土层结构，且排水效果明显.水平井具有一定埋深，大强度排水时地下水调节空间较大，水平井滤水管和土壤含水层充分接触，增加汇水面积，同时排水量可以用作盐碱土改良淋滤洗盐的补充水源.因此浅层水平排水井用于盐碱土水利工程改良具有较好的应用前景.