基于HydroGeoSphere的河谷地区地表水地下水水流水质联合模拟

张将伟，卢文喜，陈 末，林 琳，李久辉，崔尚进

(吉林大学环境与资源学院，吉林 长春 130012)

从水循环的角度来看，地表水与地下水均是水文循环的重要一环，二者在水量与水质上相互作用，紧密联系。为了更加合理的反映流域水文循环实际过程，应当将地下水与地表水系统作为一个整体进行模拟分析[1]。自20世纪70年代Pinder等对耦合模型进行了研究探索[2]后，地下水地表水耦合模型迅速发展，Govindaraju等建立地表渠系水流和地下水流耦合模型；Maxwell等将地下水地表水耦合模型ParFlow与气象模型ARPSZ结合[3]；王中根等耦合地表水模型SWAT与地下水模型MODFLOW，应用于海河流域[4]。但上述模型仅仅考虑水流耦合模拟，而忽视了地表水地下水在水质方面的联系，现代社会水量与水质均是水资源开发利用的核心制约[5]，进行地表水地下水水流水质联合模拟是发掘区域水资源系统潜力的必要手段。

目前已有多种模型可以进行地表水地下水联合模拟，其中HydroGeoSphere模型从水文循环物理过程出发，通过耦合方程将地表水地下水数学模拟模型平列起来，同时求解，可以更客观地描述水分转化关系与溶质运移过程[6]。本文综合考虑了坡面漫流、河道流与包气带水、饱和带水的相互联系，利用HydroGeoSphere软件建立地表水地下水水量水质联合模拟模型，藉此得到降雨量变化情况下研究区内水流和总氮的分布特征与变化规律，为今后进行地表水地下水联合调度、污染预报等方面研究提供科技支撑，能够更好地发挥区域水资源的整体性功能。

1 假想例简介

建立了一个南北长500 m，东西宽500 m的河谷假想例作为研究区，如图1与图1.2所示。研究区中部存在一河流，流向为北至南流动，经入流断面流入研究区(如图2中Γ4所示)，由出流断面流出(如图2中Γ2所示)。研究区东西两侧山脊线为分水岭(如图2中Γ1，Γ3所示)。研究区多年平均降雨量为800 mm/a，水面蒸发量为1 000 mm/a。研究区内土地利用类型有草地和滩地，如图2所示，其地表水水力参数取经验值如表1。地表水系统源汇项主要有河口流入量、大气降水补给、地表水体蒸发、地表水下渗补量以及河口流出量。

图1 研究区三维图Fig.1 3-D diagram of research area

图2 研究区平面图Fig.2 Planimetric map of research area

分区曼宁粗糙系数洼地储量/m消减储量/m耦合长度/m草地0.300.250.0020.01滩地0.030.250.0020.01

研究区含水层以第四系沉积物为主，水文地质参数取值如表2所示。含水层南侧和北侧(如图2中Γ2,Γ4所示)设为定水头边界，水头值分别为7 m和9 m，东侧和西侧边界(如图2中Γ1,Γ3所示)为隔水边界。研究区地下水系统源汇项有北部边界流入量、大气降水、地表水下渗量、抽水井抽水量及南部边界流出量。

表2 研究区水文地质参数表Tab.2 Value of hydrogeology parameter of study area

研究区南部存在露天垃圾堆放场，污染物以氮元素污染为主。设置一东西向监测线横穿垃圾场，监测污染物下渗情况(如图2中A-A'所示)。场地内纵向弥散系数10 m，横向弥散系数1 m。研究区东侧、北侧和西侧(如图2中Γ1,Γ3,Γ4所示)为零溶质通量边界；南侧水流交换通量边界(如图2中Γ2所示)概化为第三类边界条件。研究区南部拟建一抽水井，抽水方案如表3所示。

表3 抽水方案表Tab.3 Scheme of pumping

2 数学模拟模型

2.1 水流数学模型

2.1.1 地表水水流数学模拟模型

地表水的运动采用二维平均深度圣维南(Saint Venant)方程来描述，如模型(1)：

(1)

式中：do为地表径流深度，[L]；Kox和Koy分别为X轴和Y轴方向上的等效水力传导系数，[LT-1]；h0(x,y)为(x,y)点初始水位，[L]；Q为地表水的源汇项，为单位面积体积流量，[LT-1]；Γgs为地表水与地下水之间交换通量，[T-1]；φ0为等效地表孔隙度，无量纲；Γ1,Γ3为分水岭，Γ2,Γ4入流断面和出流断面。

2.1.2 地下水水流数学模拟模型

研究区内地下水为三维变饱和地下水流，使用Richard方程描述其水流运动，如模型(2)：

(2)

式中：Kxx,Kyy,Kzz为渗透系数，[LT-1]；Krw为相对渗透率，无量纲；H为地下水水头，[L]；W为地下水的源汇项，为单位体积流量，[T-1]；Γgs为地表水与地下水之间交换通量，[T-1]；φ为孔隙度，无量纲；Sw为饱和度，无量纲；Ss为储水系数，[L-1]；qG为边界上的水流通量，[LT-1]；H0为地下水初始水头，[L]；Γ1,Γ3为隔水边界， 已知流量边界；DG为研究区地下水区域。

2.1.3 地表水-地下水水流模型耦合

本次研究采用双重节点法对地下水地表水模型进行耦合。这种方法假设地表水与地下水之间存在一个水量交换层，地表水与地下水交换过程符合达西定律。并在达西定律的基础上引入相对渗透系数，反映包气带对地表水地下水水量交换的影响，如(3)式所示：

(3)

式中：qgs为地下水和地表水之间单位面积体积流量[LT-1]；do为地表径流深度[L]；Γgs为地表水与地下水之间交换通量[T-1]；krw为相对渗透率，无量纲；Kzz为地下水表层介质渗透系数[LT-1]；H表示地下水水头[L]；h表示地表水水头[L]；lexch为地表水与地下水之间的耦合长度[L]。

2.2 溶质运移数学模型

2.2.1 地表水溶质运移数学模拟模型

地表水中溶质二维运移过程数学模拟方程如模型(4)：

(4)

式中：qo为地表水流速，[LT-1]；c为地表水溶质浓度，ML-3；Do为地表水水动力弥散系数，L2T-1；Ωs为地表水溶质的源汇项，ML-3T-1；Ωgs为地下水地表水之间溶质交换量，ML-3T-1；Ro地表水溶质阻滞系数，无量纲；Γ1,Γ2,Γ3,Γ4为纽曼边界。

2.2.2 地下水溶质运移数学模拟模型

使用三维变饱和溶质运移模型刻画地下水溶质变化，控制方程如模型(5)：

(5)

式中：q为地下水流速，量纲[LT-1]；C为地下水溶质浓度，[ML-3]；D为水动力弥散系数，[L2T-1]；Ωg为地表水地下水间溶质的交换量，[ML-3T-1]；Ωg为地下水溶质的源汇项，[ML-3T-1]；n为孔隙度，无量纲；R为阻滞系数，无量纲；Γ1,Γ3为零通量纽曼边界，Γ2,Γ4第三类边界条件。

2.2.3 地表水-地下水溶质运移模型耦合

地表水与地下水溶质运移模型耦合方法还处于探索发展阶段[7]，本文假设地表水地下水间的溶质运移过程符合Fick定律,使用一阶多项式描述地表水地下水间的溶质运移过程如式(6)所示：

(6)

式中：Ωgs为地表水地下水的溶质交换量，[ML-3T-1]；Γ0为上游节点流量，[T-1]；C为地下水溶质浓度，[ML-3]；c为地表水溶质浓度，[ML-3]。

3 模型的离散

本文在耦合地表水模型与地下水模型后，对耦合模型进行统一的时间空间离散。

3.1 时间间离散

地表水地下水联合模拟的关键在于时间尺度[9]，采取自适应时间步长方法对联合模拟模型进行时间离散可以在满足模型精度同时有效提高计算效率。

本次研究模拟期为120 d。自适应时间步长法各项参数如下所示。初始时间0 s，初始时间步长0.5 s，单个时间步长内允许最大水头变化0.5 m，允许最大饱和度变化0.1，最小步长递增倍数0.5，最大步长递增倍数2，最大时间步长由经验公式[8](7)计算得到：

(7)

式中：S为贮水系数；a为流域边长；T为导水系数。据此计算得最大时间步长为62 500 s。

3.2 空间离散

研究区空间离散使用三角网格剖分，平面上共剖分出989个节点，1 888个单元，网格平均分布。垂向上将研究区分为15层，为研究地表水与地下水的交互作用，将地表水地下水交换界面附近剖分得更精细。如图3所示。

图3 研究区空间离散图Fig.3 Spatial discrete graph of study area

4 结果分析

抽水井以方案1方案2方案3抽水时，地表水地下水联合模拟模型运行结果如下。

方案1:模拟期末刻，抽水井内水位为7.04 m，污染物随降水下渗进入地下水，下渗至地面以下4 m，垂向污染范围约410.79 m2，如图4(a)所示。地表水出流断面径流量789.54 m3/d，污染物浓度为0.0043 g/L。研究区地表水地下水补排关系为地表水补给地下水。

方案2：模拟期末刻，抽水井水位降至6.73 m。污染物下渗4.75 m，垂向污染范围约430.29 m2，如图4(b)所示。地表水出流断面径流量为747.37 m3/d，污染物浓度为0.004 g/L。研究区地表水地下水补排关系为地表水补给地下水。

方案3：模拟期末刻，抽水井水位为6.23 m，污染物下渗5.3 m，垂向污染范围约472.52 m2，如图4(c)所示。地表水出流断面径流量为705.39 m3/d，污染物浓度为0.003 8 g/L。研究区地表水地下水补排关系为地表水补给地下水。

图4 各方案污染羽垂向分布图Fig.4 Vertical distribution of pollution plume in each scheme

5 结论与建议

5.1 结 论

(1)模拟时段内研究区地表水地下水补排关系为地表水补给地下水。随抽水井抽水量增加，地下水水位持续下降，大量地表水下渗补给地下水。当抽水井抽水量为43 m3/d时，河口地表径流量减少5.34%，当抽水井抽水量为86 m3/d时，河口地表径流量减少10.66%。

(2)地下水的开采会导致地下水污染加重。随抽水量增加，地表水地下水水头差增大，地表水下渗量增大，地表水污染物通过对流不断向地下水下渗。抽水井抽水量为43 m3/d时，地下水垂向污染范围增大4.14%。抽水井抽水量为86 m3/d时，地下水垂向污染范围增大15.03%。

5.2 建 议

(1)地表水-地下水联合模拟还处于发展阶段，但其应用前景极为广阔。可结合3S技术研究水循环过程，也可针对生态基流，水库最小下泄流量，水资源评价中地表水地下水资源重复量等专门问题进行研究。

(2)地表水-地下水联合模拟模型对水文地质资料与地表水水力资料要求高。可利用灵敏度分析性分析等方法确定敏感度较大参数，简化非敏感参数，从而达到降低资料要求同时保证模拟精度的要求。

(3)本文在运用地下水-地表水联合模拟的方法，研究地表水与地下水的水质水量的变化。但由于技术上的不成熟，仅考虑了溶质的对流弥散作用，未考虑可能发生的化学反应及生物化学反应，致使模型与实际系统仍存在一定误差。

□

[1] Maxwell R M, Miller N L. Development of a coupled land surface and groundwater model [J]. Journal of Hydrometeorology, 2005,(6):233-247.

[2] Pinder G F, Sauer S P. Numerical simulation of flood wave modification due to bank storage effects [J]. Water Resources Research, 1971, 7(1):63-70.

[3] Maxwell R M, Miller N L. Development of a coupled land surface and groundwater model [J]. Journal of Hydrometeorology, 2005,(6):233-247.

[4] 王中根,朱新军,李 尉,等.海河流域地表水与地下水耦合模拟[J].地理科学进展,2011,(11):1345-1353.

[5] 王建华,肖伟华,王浩等.变化环境下河流水量水质联合模拟与评价[J].科学通报,2013,58(12):1 101-1 108.

[6] 凌敏华,陈 喜,程勤波,等.地表水与地下水耦合模型研究进展[J].水利水电科技进展,2010,(4):79-84.

[7] 曾献奎.基于HydroGeoSphere的凌海市大、小凌河扇地地下水-地表水耦合数值模拟研究[D]. 长春：吉林大学,2009.

[8] 薛禹群,谢春红.地下水数值模拟[M].北京:科学出版社,2007:248-253.

[9] 邓 洁,魏加华,邵景力. 河渠与地下水相互转化耦合模型研究进展[J]. 南水北调与水利科技,2008,(2):75-79.