内陆河流域二维地下水运动模拟模型研究

内陆河流域二维地下水运动模拟模型研究

在干旱和半干旱地区，水资源是最宝贵的自然资源。这不仅是沙漠绿洲形成、发展和稳定的基础，也是环境的必要和有机组成部分，也是生态环境的制约因素。而地下水作为水资源的重要组成部分和天然植被赖以生存的重要水源,有时甚至是唯一水源,在干旱区水资源开发利用中占据重要的地位。随着干旱内陆河流域水土资源矛盾的日益凸显,再加之地表水资源在时间和空间上分布极不均匀,地下水资源已经成了干旱内陆河流域社会经济发展和生态环境建设必不可少的部分,是干旱区绿洲得以维持和存在的重要依靠,也是研究内陆河流域水资源统一调度管理和优化配置的关键要素。由此可见,地下水资源的可持续利用,紧密关系到国民经济的可持续发展和干旱区生态环境的重建。然而如何合理利用好地下水资源,是一个很值得深入探讨的问题。由于内陆河流域从河流源头到下游地区可划分为高山冰雪冻土带、山区植被带、平原绿洲带和下游荒漠带,而降水主要集中在高山冰雪冻土带和山区植被带,山前绿洲带和下游荒漠带降水稀少,基本不产流,因此,出山口径流基本代表了流域的水资源量,且形成与出山口以上的高山冰雪带和山区植被带。山区径流出山口后,在天然情况下部分洪水通过漫溢渗漏补给地下水,其余的通过河道渗漏补给河道两侧的地下水,最终消失于荒漠或尾闾湖中。在水资源大规模开发利用的情况下,人类活动增加了中下游的河道两侧径流的侧枝循环,通过地表径流灌溉和地下水开采,多次垂向交换,最终河水趋于消失。然而随着计算机技术的高速发展及其在各个学科领域的广泛应用,地下水数值模拟技术也得到了长足进展,国内外已经建立了许多城市和地区的地下水数值模拟模型和管理模型,例如MODFLOW、MT3D99、VizualMODFLOWF、VizualGroundwater和FEFLOW等,然而对干旱区内陆河流域地下水系统进行数值模拟相对较少。利用数值模型对地下水流和溶质运移问题进行模拟的方法以其有效性、灵活性和相对廉价性逐渐成为干旱区地下水研究领域的一种不可或缺的重要方法,在水资源评价问题、地下水污染问题、地下水管理与合理开发以及地下水-地表水联合评价等方面受到越来越大的重视和广泛的应用。因此,本文将FEFLOW模型引用到干旱内陆河流域,选取新疆台兰河流域作为研究区,开展干旱内陆流域地下水系统数值模拟研究。从而为内陆河流域地下水的模拟提供一种可借鉴的方法。为了更准确地评价研究区地下水资源,预测地下水系统状态及其变化趋势,本文在详细分析研究区地质和水文地质条件的基础上,建立研究区水文地质概念模型和地下水流模拟模型,采用基于有限元方法的FEFLOW(Finiteelementsubsurfaceflowsystem)对所建立的二维地下水流数值方程进行求解,完成模型的识别和验证,从而有助于认识台兰河流域在人类活动影响下地下水资源的发展趋势,并预测给定条件下台兰河流域地下水流场变化情况,为干旱区地下水资源的合理开发利用提供科学依据,旨在促进干旱内陆河流域地下水系统的研究。1地形地貌和地貌台兰河流域位于新疆阿克苏地区温宿县境内,地理坐标为80°21′44″～81°10′14″E,40°41′41″～42°15′13″N,地处北温带干旱区,气候干燥,日照充足,多风沙,降水稀少,蒸发较大,昼夜温差大,多年平均降水量为177mm,多年平均蒸发量为1800mm,属典型干旱大陆性气候。从地形地貌上看,流域分为丘陵区和冲洪积平原区,总体地势北高南低,可明显分为3个带,即上游冲洪积倾斜平原带、中游冲洪积扇细土平原带和下游冲洪积扇细土平原荒漠地带。根据岩性结构特征及相对关系和所处的地貌部位划分,主要为上更新统冲洪积砂砾石层全新统冲洪积物。2水文干预模型的建立2.1地下水含水层分布根据台兰河流域水文地质资料和不同时期地下水等水位线图,确定以北部的冲洪积倾斜平原带到南部的冲洪积扇的细粒土平原荒漠地带为研究区。分析现有的探孔和生产井资料,台兰河流域属潜水-承压水系统。自前山带至溢出带为潜水含水层。扇后缘靠山前带为粗颗粒的砂卵石,至扇缘前为粗中颗粒相砂砾石,向平原过渡为细颗粒相沉积物。台兰河冲洪积扇第四系含水层均大于200m。但由于其上游补给源不同,各冲洪积扇的潜水含水层的富水程度差异也很大。潜水的埋藏深度,在前山带大于100m,向扇缘方向逐渐变为50～30m,30～0m。自扇缘带以南到南部沙漠地带为潜水-承压含水层,在200m深度内,可见2～4个含水层。含水层岩性由砂砾石逐渐相变为含砾粗砂、粉细砂、亚砂、亚粘及粘土互层。南部沙漠边缘一带,含水层岩性为粉细砂,冲洪积成因的含水层厚度变薄。研究区内广泛分布第四系松散层孔隙潜水,含水层为砂碎石及中细砂层。由于受地质构造、地形地貌和沉积条件的制约,自北而南沉积物颗粒渐细,地下水位埋深逐渐变浅,富水性渐弱,含水层由北部单一的潜水含水层,向南逐渐过渡为具有统一水力联系的双层或多层结构的潜水-承压含水层。由于沉积环境的变化交替,粘性土的分布在水平方向上都有较大的变化,隔水层呈断续状或透镜体分布,很难在区域内找到一层分布比较稳定的隔水层,故承压水和潜水在区域上有不可分割的水力联系,构成了统一的含水层系统,可视为统一的多层介质的潜水含水系统。本次模拟将研究区概化成非均质各向同性二维非稳定流潜水含水层。2.2流量边界及边界根据所确定的研究范围及研究区地下水体系特征,模型的边界条件(图1)概化为:AB、MN、NA为第二类边界(补给边界),通量根据水文资料确定;HI为第二类边界(排泄边界),由于边界处为平原荒漠地带,可以认为不受抽水井的影响,但由于缺少观测井长期观测资料,因此,很难确定荒漠中地下水位的变化,所以把该边界定为定流量边界(定流量值为0.0038m/d)。其余假设为第二类零流量边界,概化后的计算面积为38.71×108m2。其北边界以地下水埋深50m为界,南至荒漠边缘为界,东西分别以邻县界为准。研究区共自动剖分三角形单元1675个,结点897个。2.3基于fewelle-qued的调整根据收集地矿、新疆生产建设兵团和水利部门钻孔抽水试验资料,确定研究区渗透系数K、给水度μ的初始值,结合研究区水文地质图和含水层剖面图,直接由FEFLOW调用调试,调试时K和μ值的改变仍然采用ArcView方便改动,并作出K和μ的分区图(图2)。由于研究资料缺少,无法准确确定隔水层底板高程,但考虑到研究区域大,水平落差小(1/1000～1/1500),进行概化认定潜水层是100m的含水层,所以把模拟区的潜水含水层隔水底板高程全部给定为一个固定值,即区域地面最低处高程减去100m的值:110m。3数学模型的建立3.1h3水层隔水板内h∂∂x(Κh∂Η∂x)+∂∂y(Κh∂Η∂y)+W(x,y,t)=μ∂Η∂tΗ(x,y,t)|t=0=Η0(x,y)(x,y)∈DΗ(x,y,t)|Γ1=Φ1(x,y,t)(x,y)∈Γ1∂Η∂n|Γ2=q2(x,y,t)(x,y)∈Γ2(1)式中,H为水位(或称水头)[L];H0:水头初始条件[L];h:潜水含水层隔水底板至潜水自由面的距离[L];K:渗透系数[L/T];D:研究区域;μ:给水度;Г1:第一类边界;Г2:第二类边界;Φ1:第一类边界条件,为给定边界的水头[L];q2:第二类边界条件,为给定边界上的单宽流量[L/T];n:为第二类边界Г2的外法线方向;W(x,y,t):地下水源汇项,即垂直方向上地下水的抽出和入渗补给[L/T]。3.21补充排水项目研究区补给项主要包括大气降水凝结水入渗、渠系水入渗、田间灌溉水入渗及地下水侧向径流;排泄项主要包括地下水开采、侧向流出、蒸发蒸腾。3.2.1凝结水入渗量台兰河流域多年平均降水量177mm,并且降水次数少、有效的降水量低。因此,降水入渗对地下水的补给必须考虑地下水埋深因素,并采取分配在某些时段的方法。根据本区已有的研究,降水、凝结水入渗补给地下水限于地下水埋深小于5.0m的地区。降水入渗以全年每次降水大于10mm的降水量之和为基数(有效降水量),台兰河流域平原区多年平均为62mm,降水入渗系数取0.03～0.20(根据地下水埋深不同有所变化),分配在降水集中的7至9三个月。凝结水在地下水位埋深1～3m区全年为27mm,埋深3～5m区为35mm。凝结水主要发生在秋季温差大、湿度大的降温过程中,因此凝结水入渗的分配是:地下水埋深1～3m区从10-05至11-05的31d中每天入渗1mm,埋深3～5m区从10-05至11-10的36d中每天入渗1mm。3.2.2灌溉水入渗补给从理论上讲,田间灌水入渗量的大小与灌溉定额、灌水量和灌溉地段的包气带岩性、厚度、湿度及气象等多种因素有关。本文用入渗系数法计算入渗量。由于灌溉水入渗补给主要发生在耕地上,将2005年台兰河流域土地利用格局图输入ArcViewGIS中,将耕地图和乡级行政图迭加,按照种植面积来确定灌溉制度,依各作物面积和灌溉制度算出每次灌水时流域某种作物的灌井水量和水库水量,然后根据渠系利用系数算出应抽井水量。同时根据不同地下水埋深耕地的面积,对迭加的图进行识别、剖分,并结合不同地下水埋深的入渗系数,从而确定接受补给耕地单元面积,由此计算田间灌溉水的入渗量(表1)。3.2.3乡单元的抽水量渠系渗漏补给量指渠系渗漏补给地下水的水量。将流域内渠系分布图输入到ArcViewGIS中,划分网格,确定渠系所在单元格,以备FEFLOW调用。在流域内按每眼机井抽水量200m3/h,每天抽12h,由分布在作物面积上的机井数量算出乡单元机井的日抽水强度,从而得到在乡单元上每天的抽水量,进一步根据式(2)计算渠系渗漏水入渗补给量:Q=qw×(1-N′)×0.9/SU(2)式中‚Q为渠系渗漏补给量‚qw为乡单元某种作物灌水时机井每天抽水量(m3/d),Ν′为渠系有效利用系数,0.9是渠系渗漏水的90%入渗补给地下水‚SU为乡单元耕地面积(m2)。3.2.4地下水侧向径流补给和排泄量计算根据文献和《阿克苏地区台兰河灌区规划报告》中计算地下水侧向流入(出)量‚计算公式为∶Q=K×B×H×I×cosα(3)式中,K为地下水渗透系数(m/d),B为过水断面宽度(m),H为地下水含水层厚度(m),I为地下水水力坡度(‰),α为地下水流向与断面线垂线夹角。3.2.5机井抽水量的测定台兰河流域机井数目多,是模拟的难点所在。按每眼机井抽水量200m3/h,每天抽12h,由分布在作物面积上的机井数量算出乡单元机井的日抽水强度,从而得到在乡单元上每天的抽水量。根据台兰河流域地下水开采强度的分区,同时将流域地下水开采强度分区图输入到ArcViewGIS中,根据强度分区将地下水的开采量分布到相应的网格单元上,以备FEFLOW调用,再结合乡单元上每天的地下水抽取量,从而计算研究区的实际地下水开采量。3.2.6蒸发量的季节变化土面蒸发指的是没有植被覆盖或植被盖度很小地区的蒸发。在稳定蒸发条件下即为潜水蒸发,可以按照不同地下水埋深计算潜水蒸发量。台兰河流域多年平均蒸发量最大在7月,其次在6、5月,然后在8、4月。据此把地下水埋深在1～3m区的陆面蒸发量按月分配,地下水埋深3～5m区的陆面蒸发量分配在7月。根据台兰河流域的气候,平原区多年平均降雨量为62mm,植物生长旺盛季节在3—10月的8个月中,因此,此间的蒸散也最强烈,但不同季节的变化很明显。因此林木草地蒸发引起的损失,在各时间段按比例分配到植物生长的几个月中(表2)。4实测结果及分析以台兰河2005年度的土地利用类型图及2004—2007年地下水位的观测数据调试模型参数并校正模型,通过多次运行模型、调参、优选,使得模型的计算结果与实测地下水数据变化趋势尽可能一致,最终确定模型选用的参数,利用10个长期观测井的地下水位实测值与模拟值进行对比。用2008年度地下水位的观测数据来验证模型。K和μ的调试结果见图2。图3是5个分布在流域上、中、下游的长期观测井(井号1、2、5、8、10)2008年度的地下水位实测值与模拟值的对比。可以看出,观测井地下水位的模拟值与实测值的变化趋势基本吻合,说明模型拟合的地下水位变化趋势与研究区的水文地质情况一致,模型的效果是较好的,可以用来模拟区域地下水的运动状态。从模拟结果来看,台兰河流域实测地下水位与模拟地下水位拟合过程比较理想,平均绝对误差0.05m,最大误差±0.5m。其它5个井的拟合情况见表3。根据各源汇项的计算结果,上述对地下水运动的数值模拟和模型校验是建立在水均衡原理上的。对整个流域来说,流域进水-流域消耗=流域排水。因此,流域消耗项即垂向源汇项的处理是地下水运动的模型校验和数值模拟的关键。其中流域消耗项主要包括渠系入渗、田间入渗、降雨入渗、潜水蒸发、植物蒸腾及地下水开采。利用2004—2007年的资料计算的平均值与利用2008年的地下水位观测资料验证模型时的结果进行比较,两者的均衡量比较接近,均衡差率为-7.45%,绝对值小于10%,表明其结果符合规定的差率范围(表4)。5建立水文模型利用ArcViewGIS与FEFLOW软件,通过对内陆河流域的潜水-微承压水系统选用二维地下潜水流模型,基于水量均衡来处理地下水的源汇项,模拟了内陆河流域地下水的演化特征。从地下水流数值模拟结果可以看出,所建立的水文地质概念模型和数学模型是正确的,选取的水文地质参数和计算的源汇项基本合理,符合台兰河流域地下水的实际情况,可用于该流域流场研究和水源地的开采规划,在此基础上可以为流域地