滦河流域内蒙段地下水资源模拟评价分析

'灌溉排水学报
2008年12月
第27卷第6期
JournalofIrrigationandDrainage文章编号:1672-3317(2008)06-0123-04*滦河流域内蒙段地下水资源模拟评价分析李凤玲,朝伦巴根,高瑞忠(内蒙古农业大学水利与土木建筑工程学院,呼和浩特010018)摘
要:以滦河流域上游闪电河子流域地下水资源评价为例,阐述了地下水资源评价中水文地质实体辨识,含水层概念模型选择,边界条件的确定,初始水头、含水层底板高程、水文地质参数空间变异性分析及用可视模块化的三维有限差模型(VISUALMODFLOW)评价地下水资源的全过程。基于地质统计学插值、估值理论,对于MODF-LOW减小误差和简化参数的调节试算过程提出了新的计算思路。关
键
词:地下水资源;模拟评价;滦河流域;MODFLOW中图分类号:P641.8


文献标志码:A1
前
言美国水文地质专家McDonaldMG和HarbaughAW(1988)的
地下水流的模块化三维有限差(MOD-[1]FLOW)
软件问世以后,20世纪90年代开始在我国得以应用。1994年加拿大WaterlooHydrogeologics公司在MODFLOW程序的基础上推出VISUALMODFLOW以来,在我国有许多成功应用的先例。杨青[2][3][4]春、杨树青、金咪用VISUALMODFLOW评价了区域地下水资源,但他们在输入初始水头、渗透系数、底板高程前,没对这些变量进行大尺度的无偏、估计方差最小的估值,因此输入的这些参数很有限。然而,VISUALMODFLOW在单元接单元的(Cel-lby-Cell)水均衡计算与水头计算时,用与待模拟水头的单元块格点最近的5个点的数据,按距离平方倒数分别对初始水头、含水层底板高程和渗透(导水)系数进行估值。显然,根据有限的初始值,对大区域内的节点值进行这样的估值误差大,影响评价精度。在输入初始水头、渗透系数、底板高程前,以无偏、估计方差最小原则的地质统计学方法,估计出足够的单元节点的这些参数值后,再输入到模型,因而保证了模拟计算的精度。2
评价区水文地质实体辨识及含水层的概念模型评价区位于滦河流域内蒙段上游的闪电河流域。开采含水层分布在闪电河河谷平原的漫滩和一级阶地上。经地质钻探、水文物探查明,研究区含水层上部为孔隙潜水含水层、下部为孔隙承压含水层,孔隙承压含水层底板岩性为上新统巴彦花群砂砾岩。评价区目前开采孔隙潜水。潜水含水层岩性主要为粉细砂,其间夹有很薄的粘土层。鉴于水资源评价区孔隙潜水含水层岩性基本为粉细沙,其概念模型为均质各向同性非稳定流模型。2.1
初始水头分布及其空间变异性分析VISUALMODFLOW3.1已把GSILIB(GEOSTATISTICALSOFWARELINBRARY)中简单、普通克里格模型嵌入到模型中,变差函数参数给出标准值。然而,水文地质参数的变差函数参数随水文地质条件的不同而变化,标准值不能涵盖具体水文地质条件对应的变差函数参数,因此,在分析初始水头分布及其空间变异性时,在VISUALMODFLOW之外采用了泛克里格法。为此,按地质统计学布设采样网的要求,在评价区选不规则分布的87个与潜水含水层平均厚度一致的开采井为初始水位采样井。以2002年9月平均*收稿日期:2007-12-10基金项目:国家自然科学基金重点项目(50139040)作者简介:李凤玲(1964-),女,呼和浩特人,副教授,博士研究生,从事水资源评价和优化配置的教学与科研工作。123 水位为初始水位。考虑采样井水位埋深受开采过程的影响,适度加密了地下水埋深观测时间。在2002年9月的1,5,10,15,20,25,30日观测了87个采样孔地下水埋深。87个地下水埋深观测井地面高程用高精度GPS测量,并用1
25000地形图内插高程值校核确认。泛克里格法用66
94(间距h=1000m)网格对初[3]始地下水位进行了无偏估计方差最小估计。边界内各节点水位就是输入到VISUALMODFLOW模型的初始水位。2.2
含水层边界与VISUALMODFLOW模型单元块边界VISUALMODFLOW模型以单元块为边界,而单元块是长方块。含水层的边界不可能与VISUALMODFLOW模型单元块边界相重合,因此用无效单元来消除落到含水层边界外的单元块的作用。将含水层以240m
300m单元块剖分。总行数(NROW)=200;总列数(NCOL)=200;总层数(NLAY)=1。2.3
含水层渗透系数及弹性贮水系数的确定及其空间变异性分析

为了确定VISUALMODFLOW模型计算需要的渗透系数K和弹性储水率SS的初值,在闪电河中、下游做了3个定流量非稳定流抽水试验。图1为其中的一个抽水试验[5][6]场图。用BP-NUEMAN、BP-BOULTON、RAGA-[7]NEUMAN法估计出K、SS。评价区还有水文地质队、水利勘测设计院利用稳定流抽水试验估计的K、SS结果。把稳定流抽水试验估计参数的结果统一到非稳定流抽水试验结果后,共有34个抽水试验求参结果。对K、SS按66
94(间距h=1000m)网格节点用泛克立格法进行无偏、估计方差最小的估值。边界内各节点的K、SS就是输入图1
抽水试验抽水孔及观测孔布置示意图VISUALMODFLOW模型的含水层系统的K、SS。2.4
含水层底板高程变异性分析评价区还有水文、煤田地质队、水利勘测设计院和课题组打的勘探孔38个。用普通克里格法按70
65(间距h=1000m)网格节点,对含水层底板高程进行无偏估计方差最小的估值。含水层边界内各节点的含水层底板高程就是输入VISUALMODFLOW模型的含水层系统含水层底板高程。2.5
定流量边界的处理确定由定流量边界流入评价区的流量时,把无流边界与VISUALMODFLOW的
井软件包
相结合,把无流边界内结点视为注水井,并把定流量分配给它们。3
VISUALMODFLOW模型模拟计算3.1
变水头单元的水头计算对于格点(i,j,k),除了有6个断面的流入流出量外,也可能有外部应力项(汇源项)。评价区的应力项有降雨补给、潜水蒸发、地下水开采、水库渗漏、河道渗漏补给。这些项有的与水头有关,有的与水头无关。必须对每一个格点都列出方程,得到一个线性方程组,用强隐式法模拟不同应力期不同时段的水头分布和水均衡计算。3.1.1
应力项的确定1)降雨应力项-降雨补给量(ai,j,k,1)的确定ai,j,k,1=Ii,j
DELRj
DELCi(1)3式中:ai,j,k,1为格点(i,j,k)的单元降雨补给量(L/T);Ii,j为面积(DELRj
DELCi)上的补给通量(L/T);DELRj,DELCi为列和行距。用安设在评价区2m高处的AZW-001自动气象站观测的降水量、TDR测量的土壤含水量、SEBA地下水位自动测仪观测的地下水位系列数据,进行多元回归得到的降雨入渗补给量与降水量、土壤含水率间的关系式为:y=-1.55+0.002304P+0.0088
40-0.0107
100+0.064
180

(R=0.76)(2)式中:P为时段降雨量(mm);
40,100,180分别为40,100,180cm处的土壤体积含水率。通过入渗系数
=Y/P124 计算出Ii,j,并输入给VISUALMODFLOW。2)蒸发应力项(ai,j,k,2)的确定用RE,i,j,k表示从单元面积DELRj
DELRi内的每单位地下水面积、单位时间蒸发损失量(L/T),把不同情况下的RE,i,j,k输入给VISUALMODFLOW,用下试计算蒸发应力项。ai,j,k,2=RE,i,j,k
DELRj
DELRi(3)

3)河床补给应力项(ai,j,k,3)的确定蓝旗电厂铁路桥到双山水库的闪电河河床有河床补给量,计算河床补给量式为:ai,j,k,3=CRIVi,j,k,3(Ri,j,k-Hi,j,k)(4)2式中:CRIVi,j,k,3为河床水力传导系数(L/T),Hi,j,k为河床平均水位。将CRIVi,j,k,3和Ri,j,k输给VISUALMODFLOW程序。4)地下水开采应力项(ai,j,k,4)的确定单元中井的抽水量或注水量与单元面积、水头相独立。ai,j,k,4等于单元中所有抽水井抽水量之和。ai,j,k,4以负值输入程序。3.1.2
模型的识别与检验评价区3个长期观测孔有水位计自动观测数据,其余4个长期观测孔有39个月每5日1次的人工观测数据。模型主要通过试算对渗透系数、储水率、降雨入渗系数、潜水蒸发系数、河床平均面积及河床推积物的平均厚度进行识别。模型识别阶段选在2002年9月1日至2004年12月31。对所有格点组成的方程组通过应力期、时间步长和迭代循环过程求解个点水头分布。识别准则为:
水位观测值与模拟值的相关系数21/2RE大、点距分布在95%的置信区间内。
平均残差平方的根RMS=最小,其中RE=Hcal-Hob;REn为残差;Hcal为水位模拟计算值;Hob为水位观测值;n为观测井数。
残差方差最小。当识别的判别准则不满足时,微调参数代入方程组重新计算。在识别阶段,7个观测孔实测水位与模拟水位的RMS301=0.393;RMS901=0.337;并且7个观测孔的水位观测值与模拟值点距都分布在95%的置信区间内;残差分别为0.199m和0.163m。上面的判定数据表明,模型识别结果可行。识别后的参数如表1所示。表1
MODFLOW模型识别参数渗透系数/弹性释水降雨入渗潜水蒸发河谷平均河床堆积平均位置给水度(m
d-1)系

数系

数系

数宽度/m厚度/m河谷180.00010.120.150.250.3漫滩130.000080.10.20.18

阶地80.000080.080.120.18



模型的检验时段选择2005年1月1日~12月31日,以模型识别校正后的参数进行模拟计算。2005年2月18日和9月1日长期观测孔的计算水位与实测水位对比,RMS分别为0.38和0.451,残差分别为0.155m和0.204m,可以得出,二者具有较好的拟合关系,从而说明了模型的有效性。3.2
水均衡计算VISUALMODFLOW模型单元水均衡计算和单元变化水头的计算彼此独立,前者可作后者的独立检验。水均衡计算是一个单元接一个单元的计算,计算结果存到磁盘以便应用。单元接单元计算的流入流出单元的水量有4种:
应力流量,它是外应力作用到单元后产生的流量、流入单元的应力水量为正,流出的为负;
储存量,它是水头变化引起的储存量释放或储存量得到的净水量。储存量释放的净水量为正,反之为负;
内部流量,它们分别是流过(i,j,k)和(i+1,j,k)、(i,j,k)和(i,j+1,k)、(i,j,k)和(i,j,k+1)单元界面的流量,流向与行、列和厚度增加方向一致的为正,相反的为负;
流量,从定水头单元流入变水头单元的流量为正,反之为负。水均衡计算的误差用下试计算。100(IN-OUT)D=(5)(IN+OUT)/2式中:D,IN,OUT分别为相对误差、总流入、总流出量。MODFLOW对相对误差D没有给出衡量计算好坏的参照值,根据海河流域水资源评价细则,取D
10%。2002~2005年水均衡计算结果列入表2。125
















3表2
2002~2005年水均衡计算结果万mOUTIN贮存量的降
雨定水头河
道定流量边贮存量的项
目蒸发量抽水量变

化补给量补给量补给量界补给量变

化6169.451144.364715.656810.5624.373.80426.355462.12合
计12029.4613397.16(IN或OUT)D/%2.69

从表2看出,D远小于10%,满足均衡计算精度。因而,从另一个角度验证了前面计算的水头分布是合理的。4
结
论通过利用MODFLOW对滦河流域闪电河子流域地下水资源的模拟评价可知,现状开采条件下研究区33多年平均净补给量为569.88万m,仍有93.06万m的开采潜力。当用水量超过此开采潜力时可通过减小排泄量和激化补给量的方法增加开采量。限于水位观测时间有限,水均衡模拟计算时间只有2.4年。随着水均衡模拟计算时间增加,相对误差会变得更小。参考文献:[1]
McDonaldMG,HarbaughAW,ChapterAI.Amodularthree-dimensionalfinite-differenceground-waterflowmodel,Book6[M].Washington:UnitedStatesGovernmentPrintingoffice,1988.[2]
杨青春,卢文喜,马洪云.VISUALMODFLOW在吉林省西部地下水数值模拟中的应用[J].水文地质工程地质,2005,32(3):67-69.[3]
杨树青,史海滨,杨金忠,等.干旱区微咸水灌溉对地下水环境影响研究[J].水利学报,2007,38(5):565-574.[4]
金咪,耿庆文.VISUALMODFLOW在淄博高新区地下水流场及水质模拟中的应用[J].华北水利水电学院学报,2003,24(3):56-59.[5]
高瑞忠,朝伦巴根,朱仲元,等.BP神经网络在Boulton法确定潜水含水层参数中的应用[J].水文地质工程地质,2005,32(3),42-45.[6]
朝伦巴根,高瑞忠,朱仲元,等.基于Neuman理论的人工神经网络优化确定含水层参数[J].工程勘察,2006,(2):19-26.[7]
高瑞忠,朝伦巴根,朱仲元,等.确定潜水含水层系统参数的Boulton-RAGA方法[J].农业工程学报,2004,(6):23-26.SimulationandAssessmentofGroundwaterResourcesintheInnerMogolianSectionofLuanheRiverBasinLIFeng-ling,Chaolunbagen,GAORu-izhong(CollegeofHydraulicandCivilEngineering,InnerMongoliaAgriculturalUniversity,Huhhot010018,China)Abstract:Thecomprehensiveassessmentofgroundwaterresourceshasbeenfinishedintermsofecologyandindustryoverallarrangementbytheresearchjobofscientificallocationofwaterresourcesintheeco-logicalconstructiveareabetweenBeijingandInnerMongolia.ThispapertakestheworkofsimulationandassessmentofgroundwaterresourcesintheInnerMongoliansectionofLuanheRiverbasinasaexampletosetforthoverviewsabouthowtoidentifyhydrogeolocialregime,tochoosegroundwatermodel,todeter-mineboundarycondition,tostudyspacevariationofinitialwatertableelevation,lowerconfiningbedele-vation,andhydrogeologicalparametersandtouseVISUALMODFLOWwhenassessmentofgroundwaterresources.Keywords:groundwaterresources;simulationandassessment;LuanheRiverbasin;MODFLOW126'