大沽河流域土壤水-地下水流耦合模拟及补给量估算*

2019-02-26张旭洋黄修东徐绍辉

土壤学报 2019年1期

张旭洋林青黄修东徐绍辉†

（1 青岛大学环境科学与工程学院，山东青岛 266071）

（2 青岛市水文局，山东青岛 266071）

大沽河水源地是青岛市主要的供水水源地之一，由于大规模开采地下水，造成地下水位大幅下降、海水入侵等一系列环境问题，给当地社会经济发展造成很大影响［1-2］。准确计算大沽河流域的地下水补给量，对保障该地区经济可持续性发展和地下水资源的合理开发利用具有重要的理论和实际意义。土壤水是大气降水、地表水、土壤水和地下水“四水转化”中的一个非常重要的环节［3-5］，与地下水之间存在着密切联系，且在一定条件下可以相互转化，综合考虑水分在包气带和地下含水层的运动过程，可以更加准确地估算地下水补给量。

估算地下水补给的方法大致可以分为物理法、化学法和数值模拟方法。物理方法主要包括水均衡法［6-7］、零通量面法、地中渗透仪法［8］、地下水位波动法［9］等；化学方法主要有人工示踪剂法［10］和氯离子质量平衡法［11-12］；数值模拟方法可以分为：基于土壤水动力学过程的渗流带模型［13］（Richards 方程）和地下水流模型［14］。由于研究对象和原理的限制，传统方法大多只单独从地下水或者包气带的角度考虑问题，忽略了土壤水和地下水之间的相互联系，导致计算结果存在一定误差。

将土壤水与地下水作为一个整体来进行研究，可以更加全面地认识土壤水、地下水的运动规律。韩双平等［15］通过调节试验区的潜水埋深，研究了在不同水位埋深下，冬小麦和夏玉米生长期内土壤水和地下水的转化关系；Adomako等［16］通过测量和分析降水、土壤水和地下水中δD和δ18O的含量差异，表明降水、土壤水和地下水之间关系密切；牛赟等［17］选取黑河荒漠过渡带作为实验地点，分析了该地区降水、土壤水和地下水的相关性，并建立三者之间的回归模型；窦超银等［18］通过观测土壤水分、盐分以及地下水位在一年内的变化情况，对地下水浅埋条件下土壤水分的运动进行了研究；卢小慧等［19］以中国科学院栾城农业生态试验站的地下水位观测资料及气象资料为基础，综合运用降水、蒸发、土壤水、地下水动态观测资料，利用EARTH 模型计算了河北平原地下水垂向入渗补给量。Chen和Hu［20］在考虑水分在包气带和地下水运动的基础上，开发出一款土壤水文模型，并利用Nebraska Sand Hills的实测数据进行验证，研究表明受到地下水位影响的土壤含水量可以提高21%，蒸发量也相应提高7%～21%。

已有研究未从水动力学角度将土壤水和地下水真正作为一个整体看待，也未建立流域尺度上的土壤水-地下水耦合模型，并以此计算大气降水和灌溉水等对地下水的补给。本研究参照Seo等［21］的HYDRUS package for MODFLOW软件的原理，结合GIS技术，建立流域尺度土壤水-地下水流耦合模型，利用研究区夏玉米生长期内的降水、蒸发、地下水位和土壤水分含量等数据对模型进行校正，并用2012—2013年的资料进行验证。

1 材料与方法

1.1 大沽河流域土壤水-地下水耦合模型的原理

土壤水-地下水耦合模型包含两大部分：HYDRUS子模块和MODFLOW主模块。HYDRUS子模块基于Šimůnek等［22］编译的HYDRUS-1D源程序，考虑了土壤水分运动的主要过程和影响因素，如降水、蒸发、渗透、毛细上升、植物根系吸水、地表积水和径流、土壤水储量等。MODFLOW模块基于Harbaugh等［23］编译的MODFLOW-2000源程序，主要用来模拟各种条件下水流在地下含水层中的运动。

HYDRUS子模块和MODFLOW模块耦合的原理主要基于Seo等［21］的HYDRUS package for MODFLOW软件，将对土壤水分运动的模拟作为一个子程序耦合到地下水运动模型MODFLOW模块中。地下水流区域的空间离散参照Harbaugh等［23］的离散方法，将含水层离散成若干网格，将地下水运动简化为平面二维流动；土壤包气带根据植被类型、地下水埋深、土壤类型以及降水分布等因素的不同，将其划分为N个不同的区域（N≤地下水网格单元），每个区域包含一个或多个地下水网格单元，该分区内所有地下水网格单元的地下水位平均值将作为土壤底部边界的水头值，土壤包气带中水分的运动则简化为垂向一维流动，这样整个研究区从地表至隔水底板就被划分为N个柱体，其中的水流运动被简化为准三维流动。每个土壤柱以地下水位为界，地下水位以下，用MODFLOW模块求解；地下水位以上，用HYDRUS-1D模块求解，在MODFLOW的每一个时间步长内，HYDRUS子模块通过多次迭代求解Richards方程得到土壤剖面底部流量，MODFLOW模块获得该流量经计算得到一个新的水位，作为下一计算步长内土壤剖面底部边界的水头值。HYDRUS模块和MODFLOW模块耦合计算的过程如图1所示。

HYDRUS-MODFLOW耦合模型中，水在土壤和地下含水层中的流动分别采用不同的方程，土壤水分运动采用修正的Richards方程，忽略土壤水侧向流动，仅考虑一维垂向运移和根系吸水项，方程求解采用伽辽金线性有限元法。方程具体形式如下：

式中，C(h)为比水容量（cm-1），C(h)=dθ/dh，θ为土壤体积含水量（cm3·cm-3）；h为压力水头（mm）；K(h)为非饱和土壤导水率（cm·d-1）；S(z, t)为t时刻z深度处耗水速率，取该处作物根系吸水率（cm3·d-1）；t为时间（d）；z为土壤深度（cm），坐标向下为正。

图1 土壤水-地下水耦合模型流程图Fig. 1 Schematic diagram of soil water/groundwater coupled model

地下水流动数值模拟的理论基础是孔隙介质中地下水二维运动方程、定解条件及数值计算方法。运动方程根据质量守恒定律导出，定解条件由研究区域地下水的初始条件与计算区域的边界条件给定，地下潜水非稳定流的二维运动偏微分方程为：

式中，K为渗透系数(m·d-1)；h为潜水层的厚度（m）；H为含水层任一点的水头标高（m）；W为单位时间单位面积上的垂直水量交换（m·d-1）；μ为给水度；t为时间（d）；q(x, y, z)为侧向补给量；φ为Γ2-1上的已知函数；Γ2-1为第一类边界，定水头边界；Γ2-2为第二类边界，隔水边界；Γ2-3为二类边界，侧向补给边界。

1.2 研究区概况

大沽河是山东半岛主要河流之一，位于东经120°03′—120°25′，北纬36°10′—37°12′之间。青岛市境内流域总面积为4 781 km2（见图2），属华北暖温带沿海湿润季风气候，多年平均气温12.3℃，多年平均降水量685.4 mm，降水主要集中在6—9月份，约占全年降水量的69.2%。多年平均蒸发量为983.9 mm，是平均降雨量的1.45倍。流域内农作物主要为冬小麦—夏玉米轮作。

1.3 大沽河流域水文地质概念模型

图2 大沽河流域的地理位置Fig. 2 Geographical location of the Dagu River Basin

土壤水-地下水流耦合模型的上边界条件为大气边界，以大气降水量、农田灌溉水量、地表蒸发量和植物蒸腾量作为模型的主要源汇项；耦合模型的下边界条件设定为隔水边界，这主要是由于大沽河流域地下含水层下伏地层为渗透能力极差的基岩；将研究区的东西两侧视为隔水边界，主要原因是研究区内的含水砂层沿大沽河周围分布、向东西两侧尖灭，且东西两边也是流域地表分水岭边界；研究区北部的含水层在大、小沽河入境处与外界的含水层相互连通，故将研究区的北部边界视为透水边界；南部边界根据课题组多年收集的地下水位实际观测资料，取多年平均水位作为定水头边界。

大沽河流域地下水的补给来源主要包括大气降水和农田灌溉；含水层的排泄项主要是潜水蒸发及人工开采。

1.4 模型离散化

研究区总面积为4 781 km2，将渗流区剖分为113行×77列，每个单元格为1 000 m×1 000 m，共划分为8 701个正方形单元，其中有效单元4 750个，非活动单元3 951 个。根据以往对大沽河水源地，即地下水库的勘察结果［24-25］，将整个流域的含水层参数分为4个区，其中地下水库区域由于含水层给水度大且差异明显，细分为3个区；地下水库以外的区域，因透水性差，故概化为一个区。给水度和渗透系数分区一致，如图3a所示。含水层垂向上简化为一层，为潜水含水层。

流域内土壤类型根据本课题组多年采样的测定结果，并通过kriging插值获得黏粒和砂粒在整个流域内不同深度的分布情况，得到相应的土壤质地；流域内的土地利用类型依据欧盟联合研究中心(Joint Research Centre，JRC)空间应用研究所(Institute of Space Applications，SAI)2000年全球土地覆盖数据（2000 Global Land Cover Data，GLC2000），将流域土地利用类型简化为耕地（3 461 km2）和城镇居民用地（1 289 km2）；由于大沽河流域地下水埋深大都小于6 m，将地下水埋深划分为4个等级：0～2 m、2～4 m、4～6 m和＞6 m；鉴于流域内降雨量空间分布的不一致性，考虑到研究区内有蓝村、岚西头、郭庄等8个气象站，采用泰森多边形方法进行分区，每一分区内的降雨条件视为均一，分区内的降雨由所属气象站数据代表。

针对流域内土壤类型、土地利用类型、地下水埋深以及降水分布等具体情况，将整个流域划分为47个区域如图3b所示，每个区域包含一个或多个网格单元，土壤包气带中的水分运动简化为一维流动。

图3 模型空间离散示意图Fig. 3 Schematic diagrams of spatial discreteness

根据模拟期间地下水埋深的变化情况，为使不同分区土壤剖面的底部边界在模拟期间始终处于地下水位以下，以含水层的隔水底板作为土壤剖面的底部边界。

1.5 模型相关参数及初始条件

以地下水的埋深为界，将土壤剖面简化为2层：下层土壤质地为砂土，代表地下水位波动带的土壤质地；上层根据实测数据细分为2～3种不同的质地，土壤水力学参数以课题组多年实测数据为准。

表1 不同类型土壤水力学参数Table 1 Soil hydraulic parameters relative to type of the soil

研究区地下水开采量参照青岛市水利局发布的2012、2013年《青岛市水资源公报》的统计结果。青岛地区的地下水开采主要集中在地下水库范围内，地下水库以外由于没有好的含水层，其渗透系数和给水度较小，开采量较低。因此，根据不同行政区浅层淡水的总开采量求得不同行政区在地下水库内、外每平方千米的开采量，将其作为耦合模型中每个单元格的地下水开采量，以抽水井的形式加入到模型中，流域内每个单元井的开采量如表2所示。

表2 大沽河流域年均供水量统计结果Table 2 Statistics of annual water supply of the Dagu River Basin

根系吸水模型采用van Genuchten宏观吸水模型，根系吸水参数选取HYDRUS-1D数据库中的玉米（wesseling，1991）和小麦(wesseling，1991)经验值［22］，由于夏玉米生长期内，大约91.6 %的根系分布在0～1 m的土壤层内，故模型中仅在0～1 m的土层设置根系吸水项。降水、灌溉水和田间蒸散发构成作物生长期内土壤水分运动模型的上边界条件和源汇项，农业灌溉量根据当地农民的实际灌溉时间和强度以降雨的方式加入到模型中，在夏玉米整个生长期内，灌溉主要集中在播种时；冬小麦生育期内灌溉一次，主要集中在拔节期前后（生育期第140天左右），每次的灌水量80 mm［26-27］。降水和蒸发数据主要来自于青岛市水利局2012年和2013年的气象统计结果。

大沽河流域土壤水-地下水耦合模型的初始流场，以模拟时间开始时的地下水位的监测资料作为依据，在流域的上、中、下游一共选取35口水位监测井的实测数据，采用内插法形成耦合模型的初始地下水流场。

由于流域内地势平坦，地表导水率较大，故模拟中不考虑地表径流。模型模拟期内，大沽河流域降水量偏小，河流大部分时间处于干涸状态，故不考虑河流入渗对地下水的补给。

2 结果与讨论

2.1 校正期土壤剖面含水量和地下水位的拟合结果

利用土壤水-地下水流耦合模型和上述参数对2013年6月16日—2013年9月29日夏玉米生长期内（共105 d）土壤水和地下水的运动过程进行模拟，将土壤剖面含水量、地下水位的模拟结果分别与实测数据进行拟合分析，采用均方根误差（RMSE）和决定系数（R2）来评价耦合模型的模拟效果。

根据已有的监测资料，随机选取大沽河流域上中下游6个土壤水分监测点（5#、3#、2#、8#、1#、10#）在7月12日的实测数据与耦合模拟所得的土壤剖面含水量进行拟合分析，除个别点外，土壤剖面含水量的R2在0.65～0.91之间，RMSE在0.005～0.01之间，因篇幅所限只给出5#、2#和10#三个监测点的拟合结果，如图4所示。由于模拟时间处于夏玉米生长旺期，玉米蒸散发量较大，而且7月9日至7月12日，大沽河流域经历了较强的降水过程，导致整个土壤剖面含水量的变化比较剧烈，模拟结果较差。5#和10#观测点0.4～0.6 m 处的实测土壤含水量明显偏高，这是因为该层土壤粘粒含量较高，持水性好，造成水分在该层集聚，从而对水分的下渗补给过程产生重要影响。土壤剖面含水量的拟合结果表明：耦合模型模拟出的土壤剖面含水量与研究区实测点的实测数据的变化趋势基本相同，模拟得到的土壤含水量与实测值总体上差别不大，说明耦合模型能较准确地模拟玉米生长期内水分在包气带中的运动过程。

图4 校正期土壤剖面含水量的拟合结果Fig. 4 Fitting of soil profile water content in the calibration phase

将耦合模拟所得的地下水位和实测水位进行拟合分析，随机选取大沽河上中下游有地下水位观测井的三个网格，绘制地下水位随时间的变化趋势线，如图5所示，地下水位的R2在0.52～0.68之间，RMSE在0.12～0.38之间。从图5中可以看出，实测地下水位在7月16日前后有明显的升高，根据青岛市水利局统计的2013年降水数据显示，在7月9日至7月12日，青岛市有明显的降雨过程，降雨量达到100 mm，由于大沽河流域平均地下水埋深在4～6 m左右，而且可能存在大孔隙优先流，较强程度的降水可以迅速入渗补给到地下水，在短时间内导致地下水位明显上升。而模型中没有考虑优先流，对于地下水位升高的现象，模拟结果滞后于实际观测结果，并且模拟曲线比实测曲线更加平缓。地下水位的拟合结果表明：土壤水-地下水流耦合模型模拟得出的地下水位与研究区实测地下水位的变化趋势基本相同，耦合模拟得到的地下水位与实测地下水位总体上差别不大，说明耦合模型能较好地模拟玉米生长期内水分在地下含水层中的运动过程。

图5 校正期地下水位的拟合结果Fig.5 Fitting of groundwater table in the calibration phase

图6为2013年8月10日研究区范围内地下水位实测值和模拟值的等值线图。从图中可以看出，在整个研究区范围内，模拟地下水位和实测地下水位大体相同，证明耦合模型可以准确地模拟流域范围内地下水位的变化趋势。

图6 地下水位模拟值和实测值等值线图Fig. 6 Contour map of simulated and measured groundwater tables

基于对土壤剖面含水量和地下水位拟合结果的分析，可以证明土壤水-地下水流耦合模型的概化及选取的土壤水力学参数和水文地质参数是合理的，校正后的土壤水-地下水流耦合模型能够比较准确地模拟大沽河流域夏玉米生长期内水分在包气带和地下含水层中的运动情况，因此耦合计算所得的地下水补给量具有较高的准确性。识别后的含水层给水度µ和渗透系数K如表3所示。

2.2 验证期土壤剖面含水量和地下水位的拟合结果

利用校正后的耦合模型对大沽河流域土壤水和地下水的运动过程进行模拟。模拟时间从2012年6月16日至2013年6月16日，模拟时间根据夏玉米和冬小麦的生长过程分为两个应力期，夏玉米期从2012年6月16日至2012年9月29日；冬小麦期从2012年9月30日至2013年6月16日。因篇幅限制，随机选取大沽河流域上中下游（K094029A、K094007A和K0920430）三个水位观测点的实测数据说明地下水位的拟合情况；随机选取上中下游（5#、2#和1#）三个土壤水分监测点在2012年7月14日和2013年3月28日的实测数据说明土壤剖面含水量的拟合结果。

表3 不同分区给水度µ和渗透系数K值表Table 3 Specific yield ( μ ) and permeation coefficient (K ) relative to zone

由图7和图8可知：验证期耦合模拟的地下水位、土壤剖面含水量与实测数据的拟合度较高：地下水位的模拟值与实测值的R2在0.56～0.81，RMSE为0.17～0.19，拟合结果较好。2012年7月14日土壤剖面含水量的R2为0.53～0.86，RMSE为0.006～0.011；2013年3月28日土壤剖面含水量的R2在0.73～0.92，RMSE为0.008～0.011。冬小麦生长期内土壤剖面含水量的拟合结果明显比夏玉米生长期的拟合结果好，这是因为在夏玉米生长期，降水和植物蒸散发作用强烈，土壤剖面含水量变化剧烈，导致拟合结果较差。验证期内地下水位和土壤剖面含水量的拟合结果表明，土壤水-地下水耦合模型能够准确的模拟研究区土壤水和地下水的运动过程，模拟所得的地下水补给量准确度和可信度较高，可以用于计算不同降雨水平年条件下地下水的补给量。

根据耦合模型的输出文件，得到含水层的水量平衡关系，进而求取地下水的补给量，地下含水层的水量平衡关系如表4所示。校正期大沽河流域大气降水量为328 mm，灌溉量为80 mm，耦合模拟计算出的地下水补给量为3.15×109m3；验证期大沽河流域大气降水量为653 mm，灌溉量为160 mm，土壤水-地下水耦合模型计算出的地下水补给量为4.77×109m3，可以为制定流域水资源优化配置方案提供科学依据。

3 结论

基于HYDRUS package for MODFLOW软件的原理，结合GIS技术，建立大沽河流域土壤水-地下水流耦合模型。根据2013年夏玉米生长期内的降水、蒸发、地下水位以及土壤剖面含水量等资料，构建2013年夏玉米生长期内土壤水-地下水耦合运动模型，经过模型校正，地下水位、土壤剖面含水量的模拟值和实测值吻合度较好，在夏玉米生长期内，整个流域范围内地下水补给量为3.15×109m3。校正期的拟合结果表明模型的概化及选取的土壤水力学参数和水文地质参数是比较合理的，校正后的耦合模型基本可以反映大沽河流域土壤水和地下水的运动状况，耦合计算所得的地下水补给量可信度较高。利用校正后的模型及参数，对2012年6月16日—2013年6月16日期间的土壤水和地下水运动情况进行模拟，结果表明：土壤水和地下水流耦合模型能够较好的反映研究区土壤水和地下水的时空变化，耦合计算所得的地下水补给量为4.77×109m3。土壤水-地下水流耦合模型能够较好地模拟大沽河流域土壤水和地下水的时空变化，并且能够准确地计算地下水垂向入渗补给量，可以为流域尺度制定更加合理的水资源优化配置方案提供科学依据。