地下水模拟读文献总结

地下水模拟中拟合的内涵_(许天福1992)

在进行地下水资源评价和地下水资源管理过程中，必须首先对建立的数学模型进行校核，使数学模型真正的刻画出水文地质实体的内部结构、功能以及受外部环境的影响。模型的校核通常采用拟合的方法，拟合就是数学模型计算的地下水位值与实际观测的地下水位值进行对比。既然拟合是模拟的地下水位与实测的地下水位之间的对比，研究拟合问题离不开地下水位。地下水位是一个复杂的自然过程在某一特定的水文地质实体中内部结构和外部环境多因子相互作用、相互影响综合的，直观的，活跃的外部表现形式，它反映了水文地质实体内物质、能量和信息的不断输人和输出及其演变关系，地下水位的某一微小的变化都是其中某一因子变动的必然结果。利用数值模拟方法研究地下水位能够从四维时空域中、从水文地质实体中的地下水与环境因子多方面联系及转换上揭示地下水位的变化规律。

从客观的复杂的水文地质实体到水文地质概念模型，再到数学模型，这一过程经过概化处理、抽象加工，最后建立的数学模型是否真正的刻画了水文地质实体的原形，通常的方法是靠拟合的效果来检验。拟合是地下水模拟计算的一个重要的不可愈越的阶段，不管采用集中参数的随机模型，还是采用分布参数的确定性模型，都要对地下水位进行拟合。拟合的效果好，说明数学模型正确，真实的揭示了水文地实体内部结构和外部环境因子相互作用，相互影响的规律，可以用建立的数学模型进行预测，以达到地下水资谏评价和合理开发利用、管理的目的。拟合效果的好坏有着丰富的深刻的内涵，告诉人们许多有用的信息。拟合的结果是对水文地质结构参数、边界条件、初始条件、源汇项等多种因素认识的集中、综合反映，拟合能对水文地质实体和外界环境的物质、能量和信息转换认识的正确与否进行高层次综合的整体的检验，是定性分析走向定量评价的一个特殊结合部位。拟合是对水文地质实体的结构、边界条件、初始条件等内外因素进行再认识的过程。地下水模拟模型是水文地质实体原形的数学抽象，而拟合效果是检验地下水模拟模型是否真实反映了客观水文地质实体的标淮，是对水文地质条件认识的正确与否的信号。地下水位是地下水系统内外诸因素相互作用、影响的综合表征，根据地下水位拟合的结果，可以对原来对水文地质条件认识的正确与否进行分析判断，对水文地质条件进行重新认识，甚至可以指导水文地质勘探，以最小的投人，获得最多的有用信息，使勘探工作目标明确。这样从实践到认识多次反复，最后会得到令人满意的拟合效果，建立的数学模型真正是客观水文地质实体的复制品，把研究水文地质实体的问题转化为研究数学模型的问题。

如采用数值模拟方法时，在开采条件下，拟合时的计算水位普遍高于实侧水位时，说明原给定的给水度..........这时就要对给水度...............水度.....................反复后就可以对给水度............识，’确定恰当的给水度分布值。当然，引起计算水位偏高还有其它原因，如初始水位给定的偏大。在拟合时，第二类边界点附近的计算水位总高于实测水位时，在其它因素正确的情况下，可能是边界侧向补给量大的缘故，从对拟合效果的分析，结合实际水文地质条件，可以得到正确的侧向补给量值，必要时可以做些补充勘探工作。拟合效果仅是事物的现象，要透过现象看事物的本质，对拟合结果进行分析判断，深人到水文地质实体实质中去以及与其相关的外部环境。从拟合结果这一线索出发，，可捕捉大量有用信息，加深对水文地啄实体深层机制的认识。拟合本身是一种数学方法，‘已的潜在意义远远超出数学本身，是客观水文地质实体的行为表现，有着具体形象的物理意义，不能把拟合作为孤立的数学过程，不要在拟合表面上大作文章，玩弄数学游戏，甚至采用“凑”的方法，单纯追求好的拟合效果，这是毫无意义的。拟合阶段是进行地下水资源评价和资源管理的基础，没有好的拟合效果，资源评价和资源管理只能是空中楼阁，拟合效果又是对水文地质实体历史行为认识的检验手段，而资源评价和管理是对水文地质实体未来行为的预侧和优化调度。好的拟合效果必须建立在对水文地质条件和外部因素正确认识的基础上。在拟合过程中，单纯的调整某一因素，对单一因素进行重新认识是不正确的，用这种方法也会得到好的拟合效果，这种好的拟合效果有时是一种假象，因为数学模型有时存在多解性。以往人们在拟合过程中，用调整参数的办法来达到好的拟合效果，即所谓的“调参”，我认为这种做法不完全正确，这种办法确定的水文地质参数有时不一定符合水文地质实体的结构。拟合过程中，对参数进行调整是必要的，不要忘记拟合效果受水文地质实体内外多因素的影响，还必须对其它因素反复的调整认识，这样的拟合才具有真实意义，这样建立的粉学模型真正是水文地质实体的复制品。拟合不仅可反映对水文地质条件认识的正确与否，而且可对数学模型中采用的地下水运动方程的墓本类型作出合乎实际的判断。以上的讨论都抛开了数值模拟方法的局限性，实际上数值模拟方法对拟合误差有一定影响，不过这种误差是有规律的系统误差，这儿不做详细讨论。客观存在的水文地质实体是错综复杂千变万化的，人们对水文地质客体的认识受技术、方法、装备及认识客观世界水平等的局限，要想一次性的，真正的、完全的、准确的把握客观水文地质实体是不可能的，数学模型是水文地质实体抽象概化的产物，并不是水文地质实体的绝对仿真，加之数值计算技术的局限性，不可能得到百分之百的拟合效果，计算水位与实侧水位绝对吻合是不可能的。只能在一定的技术经济条件下、一定精度要求的范围内，尽可能使计算水位与实际水位吻合。

地下水流系统概念模型研究_（戴振学）

地下水流系统是经过漫长地质时代的演化而形成的,其内部结构的不均匀性及外部环境因素作用程式的复杂性，使得系统的输入、输出在时空域上的变化相当复杂。在现有的经济、技术条件下，要对地下水流系统进行数学模拟与管理，必须对实体系统进行适当的概化，提炼出其本质性的东西，建立能反映实体系统外部特征和内部结构的概念模型(ConceptualModiel)，然后用数学语言加以插述，得到地下水流系统模拟的数学模型。因此，概念模型是水文地质条件的高度概括与总绪，是数学模型物理意义的表述。本文将从分析概念模型研究的历史与现状入手，论述概念模型的建模原则、方法及步骤。
 

 

一、概念模型研究的历史与现状

随着社会、经济的发展,地下水开采规模的扩大，以及对地下水流系统研究的不断深入，概念模型经历了一个由简到繁的发展过程。
19世纪中叶至本世纪初,地下水流模拟研究刚刚起步，地下水开发利用的规模比较小,开采量小于天然补给量，人们采用稳定流模型(J。Dupuit,1863)来模拟地下水流。本世纪30年代,地下水开采规模日益扩大，不少地区的地下水位持续下降，呈现出不稳定性。这时,由泰斯(C. v. Theis,1935)提出了非稳定流模型，泰斯的工作考虑到了地下水位的动态变化规律，将地下水流模拟向前推进了一大步。到了本世纪50年代,工农业生产的高速发展，对地下水需求量的增加迫使人们开始开采深部承压水，其可采资源评价又涉及到多层含水层相互越流问题。为了解决这-一问题，雅可布(C. E.Jacob)、汉土什(M.s.antush,1956)等人分别推导出了越流模型的解析
解。至此,对地下水流系统概化从单层稳态模型发展成多层越流动态模型,但对含水层内部结构与不规则
边界条件的研究仍无明显的进展。
本世纪t0年代以来，数值计算方法在水文地质学中的应用和电子计算机技术的推广使用,给区域地下水运动模拟带来了新的发展契机。计算机模拟节省了大量的人力和计算时间，使一些复杂的地下水流系统的模拟成为可能，人们开始考虑含水介质的非均质性和各向异性。复杂的越流系统和具不规则形状的各类边界条件，同时对多相流、双重介质理论也作了深入研究。这时对地下水流实体系统的概化相对要少一些，在概念模型中更多地保留了实体系统的自然特性。
但是，概念模型的研究在某种程度上被人们所忽视。在某些具体的地下水流系统模拟中,条件分析与水量计算由不同人员承担，搞条件分析的人嬴意在地下水流系统内郁绪构和边界条件上精雕细刻,从事水量计算的人鹦在数学模型及其求解方法上反复推敲,常常忽视了概念模型的研究，使得数学模型与水文地质条件相脱离,大大地降低了地下水模拟的精度，以至于有人开始怀疑数值计算的准确性。其实地下水泷数值计算已是一门比较成熟的理论，算法本身的误差是很小的，也就是说地下水流模拟失误的原因往往不在计算方法本身，而是对水文地质条件认识不清和概化不合理。因此,加强概念模型的研究非常重要。
 

二、概念模型的建模原则

通过分析概念模型发展的历史，结合我们实际工作的体会，总绪出越立概念模型应满足实用性、完整性和经济性三条原则。
1。实用性;地下水流模拟是一实用性很强的技术,解决现实问题是它的根本目的。因此,建立概念模型须与一定时期的科学技术水平相适应，能用于解决社会、经济发展中所面临的地下水模拟与管理问题。
2。完整性:概念模型必须尽可能真实全面地反映实体系统的内部结构与动态特征，保证模型在理论上的完整性,以提高地下水流系统模拟的精度。
3。经济性;概念模型须尽可能地简炼;回避枝节问题,以便用经济合理的勘探工作量来控制,用简洁的数学模型描述，且便于选择一种最优的计算方法求解数学模型,使整个地下水模拟过程在经济上最省。
上述三条相互联系.相互制约.实用性和完整性是经济性的前提条件，完整性又是实用性的技术保证s完整性与经济性相互对立，但最终统一在实用性上。
 

三、概念模型的建模方法及步异

水文地质条件是概念模型的基础。在建立概念模型之前,必须认真收集、整理和分析已有的水文地质资料,确定模拟的目的层,进而勾出职r水公体#统的内部结构与边界条件,然后才犬炉对次A法湖行概化。这里以焦作岩溶水系统为例来论述概念模型·
的建模方法及步骤。
一)地下水流系统内部结构概化
地下水流系统内部结构包括含水介质的形态、地下水运动状态、参数时空分布规律。
'
1。含水介质
 

达西定律及地下水流基本渗流方程都是在多孔介质中推导出的,要用现有的渗流模型描述各种介质中的地下水流,首先必须将含水介质概化为多孔介质。
 

将孔隙介质概化为非均质、各向同性或各向异性多孔介质已为人们所接受,但裂隙、岩溶含水介质就要视具体情况而定了。就焦作岩溶水系统来说,含水层由奥陶系碳酸盐岩构成，含水介质以溶隙、溶孔为主。

其次为溶洞。无可否认,在局部溶洞发育处,岩溶水运动为非达西流(即非线性流和紊流),但实际模拟经验告诉我们在大区域上，北方岩溶水运动近似地满足达西定律。因此，我们将本例的含水介质概化为非均质、各向异性的连续介质。
 

2地下水运动状态
自然界一切现象都是在三维空间里发生的,
地下
水也不例外。但考虑到模型的实用性和经济性,
在给
定假设（如裘布衣假设)的基础上,常常将地下水流概化为平面二维流。在本例中，虽然奥陶系岩溶发育具有多层性,但区域性的断裂及岩溶裂隙在垂向上的导水作用,使各岩溶发育层之间的连通性较好,
且具
有统一的地下水位，同时考虑到勘探工作量的限制,
将本岩溶含水层概化为单层结构,即岩溶水为平面二维流。
 

对于多层含水层越浇系统,在有勘探工程控制的情况下，常来用准三维流来模拟,即将地下水泷概化为:在含水层中为平面二维流,在越流层中为垂向一维流。当矿床开采前进行大降深疏干含水层时，地下水为明显的三维流，并且在疏干前大都作过大降深、大流量的抽（放）水试验,丰富的水文地质资料便于采用三雄流模型进行疏干预报。
 

3。水文地质参数的时空分布规律
水文地质参数是时间和空间坐标的的数，受目前的勘探技术条件和参数识别方法的限制,我们还不能反演出这种复杂的函数关系。但是，水文地质参数是慢时变的。在一定时期和外部条件下可以近似地看作恒定不变，因此，建立概念模型时,我们几乎都是将参数概化为时不变的。对于参数的空间分布规律,
常
采用离散化的参数概化方法（即参数分区或参数化)
 

来确定。一个好的参数分区能防止用于参数识别的非线性规划收敛于局部最优解,减少以至消除参敷识别问题的不适定性。
参数分区的依据可以归纳为以下几点:
(1)计算区单孔抽水试骏资料的计算结果,包括渗透系数、储水系数、给水度及单位涌水量。
(2〉含水层分布规律，即埋深、厚度和岩性组合特征。
( 3)地下水天然流场、人工干扰流场、水化学场和温度场。
()构道条件及岩溶发育规律(限于岩溶含水层)。
通过对以上四方面的资料分析，我们划出了焦作岩溶水系统参数分区的界线，并给出了各分区参数值的上、下限、初始值及主渗透方向，为利用最优化算法辨识水文地质参数打下了基础。

(二）边界条件的确定与概化

根据边界条件的性质可将其分为三类:
第一类(Dirichlet型):又称水位边界。选取水位边界应注意以下几点:(1)水位边界的位置应尽可能地远离计算区内的源（汇)项,绝对不能置抽(注)水井于水位边界上,〈2〉水位边界处要有观测点控制，以确定边界水位值( 3)在榄型城中至少应有一个水位边界结点,这对保证数值棋型解的唯一性是必要的。
第二类(Noumann型)。又称流量边界。选取二类边界应以隔水边界和弱透水边界为主，尽量不用人工划成的大流量边界。边界流量大都是根据达西定律估算出的，边界处的水力坡度和导水系数的估值误差将会导致边界流量的较大误差，并且在数值模型中处理大流量边界容易遣成边界附近结点的水位异常。
第三类:又称混合边界条件。由于边界条件中的两个常系数较难准确估值,在实际工作中应用的较少。
 

在确定焦作岩溶水系统边界时，我们主要选取了确定性的自然边界作为计算边界。东部青羊口断裂和‘南部朱村断裂均为区域性的深大断裂,构成隔水边界,北部以地表及地下分水岭作为零流量边界:西部在三姑泉一带划出一段水位边界，考虑到边界水位随季节变化，选用定水位边界是不合适的,故根据观测点长期观测资料，采用富氏级数模拟并预报边界水位随时间的变化规律,与数值模型联立求解。
三)地下水系统的输入与输出
地下水系统中所有源（汇)项,即外部环境对地下水系统的补给量与排泄量统称为输入,而相应的水位动态变化称为输出。对输入银件的概化，用来确定地下水流模型中源（汇)项的大小、时空变化规律p对输出状态的研究以帮助分析地下水系统的动态特征，同时选定数学模型的初始流场。
焦作岩溶水系统的输入主要由大气降水人渗量和旷坑排水量构成,分别用注水井和抽水井来描述。由于大气降水址是一个随机变量，在数值预报模型中，我们根据1959年至1989年的气象资料，采用随机模型预报了本区未来的降水量，从而提高了数值模型的预报精度。

(四）概念模型的校正

概念模型很难一次建立好，往往需要在模型识别过程中加以校正。概念模型的校正实质上是对水文地质条件的再认识，根据模型识别时的水位拟合误差分布规律,对参数分区、边界条件及源〈汇)项进行调整。概念模型的校正要以实际水文地质资料为依据,如果仅根据观测孔水位拟合绪果随意调整概念模型,就会陷入“数值游戏”的泥潭，即使这样能将观测孔水位拟合好,但所得到的模型也不具准确预报的能力。因此，概念模型的校正必须谶慎从事。
本例在模型识别中，采用修正的高斯-牛顿法进行参数优选,由于对水文地质条件进行了认真的概化,槟型识别工作进行得比较顺利,，仅对参数分区作过小的调整，就识别出了模型的水文地质参数，得到了焦作岩溶水系统的模拟模型。
地下水流系统概念模型是连接地下水实体系统与数学模型的桥梁。加强概念模型的研究,有助于加强基础水文地质理论的研究，改进传统的水文地质勘探方法。使之能满足地下水流数值模拟的技术要求，提高数学模型的可靠性和预报精度,加强概念模型的研究，将推动具有实用意义的岩溶、裂隙水运动模型的研究，不断完善和提高地下水流系统的模拟技术。

数值模拟是反映客观规律和定量评价的重要手段

70年代初以来，随着电子计算机和数值计算方法的发展,数值模拟逐渐取代传统的模拟技术成为定量评价的主要手段:水量计算、评价地下水资源遘都离不开它,还可以用它来研究地下水的合理开发和地下水资源的科学管理,但其发展趋势已远远超出作为一种计算手段的原有范畴，逐渐成为槟拟,再现自然界发生、发展的一些水文地质工程地质过程的重要手段。如一些大型自流水盆地内地下水的径流过程,向含水层内注人热水（冷水)后热量的传输过程，抽取地下水与地面沉降的关系,海岸带海水入侵含水层过程中海水与淡水的混合及过渡带（混舍带）的演化。

发展过程，污染物在地下水中的运移过程与溶质在含水层内行为（如氮素进入地下水后，在氧化条件下如何由氨转变为亚硝酸盐，最终成为硝酸盐的过程）等等诸如此类，复杂的水动力学过程、水化学过程现在都有可

能通过数值模拟再现出来。·

例如通过一个三维海水入侵模型，我们再现了1989~1990年山东省龙口市黄河营地区海水入侵的全过程与咸、淡水界面的运移规律。随着时间的推移，能够模拟、再现的内容越来越多，精度越来越高。

定性分析与数值揍拟结合的具体困难和解决办法（吕康林）

定性分析就是对水文地质侧绘、钻探物探、动态观侧，抽水试验等资料进行综合分析，研究地下水的补给、径流和排泄规律及含水层的结构特征娜。它是进行各种水文地质工作的签础，当然也是进行数值模拟的基础。因此，无论怎么强调它的重要性也是无可非议的。然而，在数值裤拟的实跳中我们体会到，即使有这种意识，也难以使定性分析与数值撰拟充分地结合。现就这一问题谈一下自己粗浅的肴法。首先，数值模拟的复杂性限制了它们的绪合。由于数值模拟技术比较复杂，涉及的知识面较广，目前真正攀握这一技术的专家主要集中在高等院校和科研单位。而对水文地质条件最熟悉的、理解最深刻的往往是长期工作于现场的技术人员或传统的水文地质专家。这样来，就使“梅拟技术”和“水文地质条件”分别掌握在不同的两批人手中。这种状况消弱了定性分析对数值模拟的指导作用。当然，一些模拟专家也是优秀的定性分析专家。但是由于时间和所提供资料的限制，他们对水文地质条件的认识也随之受到了限制。例如，生产单位提供的资料常常是一本报告和大批的观侧数据.............容易形成所谓的“数字游戏”申..........份报告或现场人员所介绍的情况建立模型的。为了使定性分析与数值模拟更好地结合，就必须推广数值法的应用。使那些具有丰言经雏的现场工作者和传统的水文地质专家能够直接或里多地参与拟合。这里的..................论和算法语官，而是通过开发软件的方践来实现。现在的软件技术发展很快，功胡越来越强。开发一个“像瓜”数值棋拟软件是完全可能的。使用这种软件，没有必要弄懂内部的计算方法和穆序结构，只要根据屏幕上的“莱单，提示和简单的使用指导书就可以进行数值椒拟.........................软件那样...................条件最熟的人进行数值拟合，有利于加强定性分析的指导作用，有利于提高模拟精度。影响定性分析与数值模拟维合的第二个因晰是精确数学方法的局限性。这种局限性表现在两个方面.定性分析结果难以准确地定堆化.........性和计算方铸的高度精确性的矛盾。大家知道，定性分析得出的给论常常是“导水性强、补给源丰富、有水力联系或关系密切”等诸如此类的模栩的相对的概念。由于各人对啊强弱.....受不一样，尚使用精确数攀方法，这些相对概念在数值上的变化可能很大。也批悬说这些给论难以..因而也鱿难以对数健润拟存寒斌性的匆导京义。如常见的河流与地下水的关系，，妞设定性分析得出的结论是“河流补给地下水，，我们如何用这钾结论指导数值模拟........................也是补给...................子可以看出，一个定性分析结果能引导出多个数值给果。这种现象多了，筑增大了数值解的非唯一性。遗成这种现象最根本的原因是精确数学方法。因为它要求一个变最只能俄一个“精确”悠而定性分析却不能准确地给出这个值。事实上，定性分析给果反映的“量”往往是一个值域...............可能性较大的数值，此数值只有对水文地质条件理解较深的人才能感受到，单凭文字摘述是反映不出来的。

地下水数值模拟模型拟合效果的评价（潘国营）

白沙河平原地下水系统数值模拟（邱汉）

郑州市地下水渗流场的数值模拟和优化管理1_（宫 辉 力）

0 引 言

郑州市长期过量开采地下水资源,致使区域地下水位大幅度下降,由此引发了一系列水环境问题[1]:(1)区域地下水降落漏斗不断向深陡扩展;(2)城市供水难度和费用大幅度增加;(3)导致地下建筑工程充水;(4)诱发地面塌陷和沉降。本项研究在建立地下水渗流场的数值模拟模型和优化管理模型及其耦合模型的基础上,通过模型运行优化地下水开采管理方案,进而调控地下水渗流场来解决上述不合理开采地下水资源所引发环境问题。

1 水文地质概念模型

郑州市地下水系统为多层含水系统(约900 km2)。在埋深小于250m范围内,可以概化为[1]:(1)主要由第四纪全新统和上更新统地层组成的浅层地下水系统(底板埋深60 m);(2)由中更新统和下更新统地层组成的中深层地下水系统(深度80~250m)。两者通过弱透水层越流系统发生水力联系,构成统一、非均质各向同性的准三维(浅层和中深层地下水系统为平面二维流,弱透水层越流系统为垂向一维流)地下水渗流场。数值模拟计算域分别由已知水头边界(Dirichlet型)和已知流量边界(Neumann型)[2]构成。系统的宏观输入、输出水量实测获得。

2 数值模拟模型

根据水文地质概念模型,构建描述计算域水位分布的数学模型[2,3]。当不考虑弱透水层贮水能力时,有如下模型。浅层地下水系统数学模型(Ⅰ):

中深层地下水系统数学模型(Ⅱ):

式中K为浅层地下水系统的平均渗透系数(m/d);T为中深层地下水系统导水系数(m2/d);h,H为浅层地下水系统和中深层地下水系统水位标高(m);L,S为给水度和贮水系数(无量纲);k′,M′为间隔弱透水层的厚度和垂向渗透系数(m/d);q1为浅层地下水系统垂直补给强度代数和(1/d),取正值;q2为中深层地下水系统垂直补给强度代数和(1/d),取负值;b1,b2为弱透水层的顶、底部边界标高(m);D为浅层地下水系统和中深层地下水系统同时确定的平面渗流域(m2);#(1)1,#(2)1—渗流计算域D的已知水位边界;#(1)2,#(2)2为渗流计算域D的已知流量边界;n为二类边界内法线向量。

水资源系统模拟模型及其应用_（张建中）

地下水流数值模拟的研究现状和发展趋势（魏林宏）

20多年来，随着电子计算机和数值方法的发展，数值模拟逐渐取代传统的模拟技术，成为研究地下水运动规律和定量评价地下水资源的主要手段，而且其发展趋势已远远超出作为一种计算手段的原有范畴，成为模拟一些水文地质过程发生、发展的主要手段。1 地下水流数值模拟的研究现状无论从理论上还是实践上，对于多孔介质饱和带中地下水流运动机理的认识是明确的，模拟饱和带地下水流的模型是成熟的，应用也是最广泛的（周仰效，1995）。已有许多使用灵活的计算机软件可供使用。近年来对于饱和带地下水流模拟的研究主要集中在：①三维流模型开发；②流速场与流线的计算方法；③非均质参数的区域概化；④繁杂数据的优化处理。目前进行区域三维地下水流的分析的主要软件有HST3D（Kipp1987），SWIFT（Granwell和Reeves，1981）以及目前世界上最流行的软件VisualMOD-FLOW（McDonald和Harbangh，1994）。VisualMODFLOW用于地下水流模拟，具有以下几个主要的功能：其一是水质点的向前、向后示踪流线模拟研究。根据地下水稳定流数据模拟结果，Modpath可方便地计算出三维流线分布和任意时间水质点的移动位置；其二是任意水均衡域的水均衡项的研究，ZoneBudge是用于计算任意水均衡域均衡要素的专门模块。用户根据自己的研究目的和需要，可在模拟区域任意选定水均衡计算的均衡区，通过执行ZoneBudge模块，可方便地得到所选定均衡域和整个模型区的所有水均衡信息。该功能可帮助地下水资源研究者直接确定研究区的均衡项及其数值的大小。此外，它的另一个重要用途就是通过断裂构造所在均衡域的水均衡计算，预测断裂的水通量；其三是允许用户直接接受地理信息系统（GIS）的输出数据文件和各种图形文件，这个功能对于充分发挥具有强大空间信息处理与分析功能的GIS技术在数值模拟中的作用意义重大，也为开发研制地下水三维数值模拟与GIS耦合评价奠定了坚实的基础。对于流速场与流线的计算方法，Pollock（1998）开发了计算流线的方法，相应的软件MODPATH在实践中已得到了广泛的应用。对于有限元法，由于获得的流速场在单元边界上的不连续性，导致计算流线的误差。目前已提出了克服这一问题的几种方法，一是基于流速分量的有限元法（Zijl，1984），另一种方法是流函数模型（Schmid，1981；Willms，1987；Frind等，1985），然而前者工作量大大增加，后者无法处理面壮的源、汇项。1992年Cordes与Kinzlbach提出了改善流速场不连续的方法。水文地质参数在地下水数值模拟中的地位是举足轻重的，模型中参数的不确定性将导致计算的水头、流速的不确定性，并影响到地下水资源评价结果的可靠性（束龙仓等，2000）。饱和带地下水流模拟的主要问题仍然是非均质参数的区域概化。由于勘探孔的数量有限，得出的含水层参数是局部的，不能以此来概化整个含水层的参数分布。理论研究注重于确定区域渗透系数与局部渗透系数的关系（Dagan，1986；Gelhar，1986）；渗透系数和孔隙度的空间分布与地质结构的关系（Byers和Stephens1983，Hoeksema和Kitanidis1985，Smith1981，Shclicdy1986）；导水系数与单位涌水量的关系（PeterM．Meier，JesusCarrear和XavierSanchez-vila，1998）。在模型实践方面，Sun和Yeh（1985）最先在地下水流反问题中考虑参数的模型结构识别，他们开发了一种优化方法不仅可用来识别模型的结构，而且可用来识别含水层的参数。但这种方法最终并未运用到实践中。当观测数据有限时，增加非均质含水层的分区数会降低识别参数的可靠性（Shah等1978，Yeh和Yoon1981），因此采用优化分区模型方法（YunjungHyun和Kang－kunlee1997）进行参数识别。该方法应用4个识别标准（Carrera和Neu-man1986a）和最大概率理论（Crrera和Neuman1986b）求得参数模型结构的最优解，提高了含水层参数的精度。另外近年来发展起来的人工神经网络技术在地下水流模拟中得到了应用，Abdul－Aziz和Wong于1992年开发了人工神经网络技术在估算含水层导水系数、贮水系数、越流系数方面的应用。用这种方法确定的含水层参数比传统的克里格（Kring）方法更准确（AmitabhaMukhopadhyay1998），关于人工神经网络技术Lawrance已在有关文献中详细说明（1994）。数值法庞大数据的处理一直是水文地质工作者最头疼的事。Roaza和Mcgrath于1992年对一地下水模型进行改善，使其能够直接读取地理信息系统所用的文件；Richard等于1993年在分析某一地区水流和污染质运移问题中，把地理信息系统和有限元法进行了成功的结合；Kuniansky和Lowther（1993）应用地理信息系统研制了一种自动生成有限元网格的数值模拟软件。这些方法和软件的出现，免去了文件格式间的转换，同时也避免了数据转换带来的误差，提高了模拟结果的精度。综上所述，人们对三维水流模拟软件的开发、对流速场和流线的计算方法的研究、对水文地质参数的概化，以及开发软件来解决繁杂数据的优化处理问题，无疑大大地提高了地下水流数值模拟的精度。2 水文地质数值模型基础的薄弱性与相邻学科，如水文学、水力学相比，水文地质学的发展相对而言是不成熟的，在构造模型，涉及多学科的全球变化研究中必须注意到水文地质数值模型的理论和数学基础的薄弱性。这表现在某些基本的地下水流动问题未能得到很好的解决。例如，作为非稳定流理论基础的泰斯模型，其原型在野外很少存在。即使汉士什，雅各布以及其它模型，与实际的物理模型相比，仍存在很大差距，可以说，未将地层的力学性质引入承压水流动模型，不考虑地下水运动的不可逆过程是水文地质数值模型的一个重要缺陷（汪民，1992）；连续性原理和达西定律是地下水动力学中两个最重要、最基本的定律，但它们是建立在连续介质和渗流两个对含水层进行了抽象的概念之上的。严格地说含水层是一个非连续介质，达西定律的应用也有一定的局限性，在非层流地区，误差是不可避免的。因此，以这两个定律为基础所推导的运动方程来描述地下水运动过程本身就存在固有的误差。水文地质数学模型数学基础薄弱性表现在模型反演求参时，解的不唯一问题。水文地质数值模型中，通常含有两个以上参数（如潜水含水层的渗透系数、给水度，承压含水层的导水系数、弹性储水系数等），从数学角度看，反演求参时，除非增加一些额外的约束条件，否则其解不唯一，有时甚至不稳定。从而造成一些人甚至对数值模拟方法本身的科学性产生了怀疑。另外，地下水数值模型建立时的诸多假设也会对地下水流数值模拟结果带来影响。3 地下水流数值模拟的拟合地下水流数值模拟的核心是拟合，正确理解和进行拟合对于提高数值模拟的质量是至关重要的。地下水流数值模型是描述地下水流的形成、机制、规律或其中某些侧面的摹本。它是人们对地下水流客观认识的一种表达形式。具体的地下水流是模型的原型。拟合则是要求数学模型与水流原型达到逼真的相似，即达到与原型等效，进行拟合的过程就是不断改正错误认识，不断修正数学模型。对于地下水流问题，要使模型和原型拟合很好，必须对原型有足够的认识，即对含水层的结构、边界条件，地下水的形成机制以及运动规律等有全面了解，这样才容易拟合。数值模拟要通过人来完成。现在有些模拟者，既不尊重事实，也不尊重科学，甚至弄虚作假，在实测数据和模拟结果的选用上大作文章，任意剔除不“理想”的数据和结果。基于上述情况必须要求模拟者有良好的素质，首先，必须尊重事实，尊重科学。其次，必须在业务上具备坚实的地质、水文地质知识和野外地质工作能力，同时具备扎实的数学基础和计算机编程和实际操作技能。4 地下水流数值模拟的发展趋势有问题才有发展，相信在21世纪地下水流数值模拟将在以下几个方面取得重大突破：数值模型基础理论研究将有重大突破；引进新的思维方法，新的数学工具，合理地描述含水层系统中大量的不确定性和模糊因素；三维流模型所要求的勘探工作量在目前的一般工程中还难以承担，所以现在一般都建立二维流或准三维流模型，随着勘测手段的提高，三维流模型将得到普及；国际上流行的VisualMODFLOW、ASMWIN等数值计算软件将得到普及，使人们从繁杂的手工数据整理、输入等工作中解放出来，会大大提高数值模型的使用效率；数值模型在裂隙发育区、非饱和带、海岸带的应用将逐渐成熟，在这些区域依靠数值模型预测水流的主要问题是含水层特征参数的确定，随着观测手段的改进和新的数学方法运用，含水层参数的精度将有所提高，因此这些地区的数学模型将会得到不断的完善。另外，随着计算机科学的飞速发展，以及遥感、地理信息系统和全球定位系统在地下水数值模拟中的应用，人们不仅可以直观地模拟地下水流，而且可以实时地监测地下水的动态，地下水水流的数值模拟将进入一个崭新的时代。

由MODFLOW浅谈地下水流数值模拟软件的发展趋势（吴剑锋）

MODFLOW（ModularThree-dimensionalFi-niteDifferenceGroundwaterFlowModel）是由美国地质调查局（U．S．GeologicalSurvey）的McDonald和Harbaugh于80年代开发出来的一套专门用于孔隙介质中三维有限差分地下水流数值模拟的软件［1］。自从它问世以来，MODFLOW已经在全美甚至在全世界范围内，在科研、生产、环境保护、城乡发展规划、水资源利用等许多行业和部门得到了广泛的应用，成为最为普及的地下水运动数值模拟的计算软件。笔者认为，MODFLOW之所以得到如此广泛的推广应用，是因为它代表了未来地下水流数值模型大发展趋势，有很强的实用性，具体包括程序设计结构的模块化，离散方法的简单化和求解方法的多样化。1 程序结构的模块化MODFLOW的一个显著特点就是它采用了模块化结构（ModularStructure），包括一个主程序（MainProgram）和若干个相对独立的子程序包（Package）。每个子程序中有数个模块（Module），每个模块用以完成数值模拟的一部分。例如“河流子程序包”用来模拟河流对含水层的影响；“井流子程序包”用来模拟抽水井和注水井对含水层的影响。用户可以根据实际工作需要选用其中某些相关的子程序包对地下水进行数值模拟。从辨证的角度来说，任何事物的发展都不能一次至臻完美，需要经过不断地补充、修改和完善过程，而MODFLOW的这种模块化结构使得其程序易于理解、便于操作、修改以至添加新的子程序包。事实上，MODFLOW问世以来，已经有许多新的子程序包被开发出来，用来解决一些原来MODFLOW不能解决的问题，例如用来模拟抽水引起地面沉降的子程序包（Leake和Prudic，1998），用来模拟水平流动障碍（Horizontalflow-bar-rier）的子程序包（Hsieh和Freckleton，1993）等。这些新子程序的加入，很大提高了MODFLOW的应用范围。据统计，自1992年底以来，美国地质调查局运行的22种有关地下水流水量水质数值模拟计算的程序中（表1），MODFLOW约占总应用次数的41.56％，这还不包括附属于它MODFLOWP和MODPATH，而其它的绝大多数模型的使用率都不超过5.0％，可以说MODFLOW具有绝对的权威性。如今，MODFLOW不但被广大的水文地质工作者所接受，而且已被其政府部门的认同。运用MOD-FLOW的计算结果，在大量有关地下水资源超量开采、地下水污染等诉讼案件中，甚至成为一个有力的法律依据。

回顾我国自80年代推广应用数值模拟方法解决地下水运动问题以来，虽然建立了很多模型，取得了不小的成就，但迄今为止，尤其是90年代中期以来，随着计算机技术迅猛发展，我们还没有一个通用的权威计算软件，这不能不说是一大遗憾。造成这一尴尬局面既有客观上的原因，也有主观上的原因。客观上我国的经济不发达，而要设计一个大型的计算软件，又需要大量的资金，再说，投入到这一领域的资金又不能产生“立杆见影”的经济效益，自然地国家、企业、个人都不愿意把资金投入到这一科研领域。主观上：（1）各个水文地质工作者存在门户之见，各人只认自己“门户”的计算模型和结果，互相不信赖，这样，即使有些好的模型也难于推广应用；（2）缺少多学科的协作精神，大型的地下水流模拟软件，不仅牵涉到水文地质学科，而且涉及到应用数学、计算机科学等领域，因此需要有多个学科的专业人员协作研究，MODFLOW的设计人员中，包括很多具有丰富实际经验的水文地质学家、数学家和计算机科学家。事实上，设计一个通用的结构化程序，虽然一次性投入的资金较大，但从长期来看，它的产出应该远大于投入，多数人还没有认识到这一点。这一问题类似于我国地下水资源评价与管理的脱节［2］，归根到一点还是眼光短浅，缺乏高瞻远瞩的大家风范。2 离散方法的简单化MODFLOW采用有限差分法而不是有限单元法对地下水流进行数值模拟。我们知道，有限差分方法曾是解决地下水流运动占统治地位的一种数值方法，近二、三十年来，有限单元法受到越来越多人的注意，被认为是本世纪最有影响的数值方法，并大有取代有限差分法之势。不幸的是，这种趋势的发展并不是基于对这两种数值方法理论上的对比分析，而更多的是出于对有限单元法中广泛使用的比较灵活的网格剖分的爱好，以及对复杂而抽象的数学推导所产生的数值方法的期望。实际上，有限单元法和有限差分法没有太大的差别，两者可以统一起来，但是，有限单元法解决非稳定的地下水流运动问题，在时间步长Δt较小的情况下，某些单元可能出现质量不守恒，因此会引起个别点的水头反常［3，4，5］。这些水头异常点往往出现在变化较为剧烈的源汇项附近，而这些地方又经常是我们重点关心的地带。相对于此，有限差分法就不存在这种问题。由此说来，MOD-FLOW采用有限差分方法更为合理可靠。在空间上，MODFLOW对含水层采用等距或不等距正交的长方体剖分网格。这种网格的优点在于用户易于准备数据文件，便于输入文件的规范化。其最大缺点在于增加了许多额外的计算单元（如图1），（a）为含水层空间剖分平面图的一部分，如果我们需要将图中的单元网格A加密，结果不但使A分化为更小的单元，而且使得横向和纵向上与单元A平行的B、C、D、E、F等单元都被加密（图1（b））。而这些单元我们主观上并不想加密，这样就增加了许多剖分单元，延长了程序的运行时间。采用不规则网格剖分能够有效地解决这一问题，并且对于单一的具体问题更为灵活方便，现阶段我国的水文地质工作者多采用这一方法。然而，不规则网格却不利于程序的规范化和普及化。随着计算机科学的迅猛发展，计算机受网格数量的限制越来越小，MODFLOW解决地下水流运动问题已经将含水层剖分到多达360×360×18个网格单元。在时间上，MODFLOW引进了应力期（StressPeriod）概念，它将整个模拟时间分为若干个应力期，每个应力期又可再分为若干个时间段（TimeStep）。在同一应力期，各时间段既可以按等步长，也可以按一个规定的几何序列逐渐增长。而在每个应力期内，MODFLOW规定所有的外部源汇项的强度应保持不变。这样做不但简化规范了数据文件的输入，而且使得物理概念更为明确。

3 求解方法的多样化求解的方法可以分为直接求解方法和迭代求解方法。MODFLOW原来含有两种迭代求解子程序包：SIP（StronglyImplicitProcedure）方法或称为强隐式法，SOR（SuccessiveOverrelaxation）方法或称为逐次超松弛迭代法。由于MODFLOW的模块化结构，MaryHill于1990年设计增加了一种新的迭代子程序包：PCG2子程序包，该子程序包采用PCG（Pre-conditionadConjugateGradient）方法或称为预调共轭梯度法迭代求解。Harbaugh也于1996年增加了一种直接求解子程序包，但由于直接求解需要大量的由图2可知，SIP法和PCG法精度高，而SOR方法求得的结果较理论值明显偏小，具体原因还不十分清楚，可能与MODFLOW采用SOR方法是一种片状超松弛迭代法（SSOR）有关。该方法将模型的网格分割为一系列“片”（Slices），并将计算单元按所在内存，现阶段还不具有实用价值。MODFLOW的多个求解子程序包，一方面，用户可以根据问题的实际情况选用比较合适的求解方法；另一方面，对于某一特定的实际问题，由于水文地质条件的复杂，用户选择不同的求解子程序包，可能都会收敛，也可能只收敛于一种（或几种）求解方法而不收敛于另一种（或几种）求解方法。作者通过分析选用不同求解子程序包的多次计算结果，认为SIP法和PCG法实用可靠，而运用SOR子程序包求解的结果精度低，不宜采用。究其原因，我们可以通过下面的简单实例得到论证。［实例］ 设某一含水层均质等厚，各向同性，水平侧向无限延伸，天然状态下水力坡度为零，其导水系数T＝2.3×10－3m2／s，贮水系数μ*＝7.5×10－4。该含水层中有一口完整井进行抽水试验，抽水量Q＝4.0×10－3m3／s，持续时间t＝86400s，距抽水井20m处有一观测孔。试求抽水井停抽时刻的距离-降深曲线以及抽水过程中观测孔的时间-降深曲线。根据问题的实际情况，作者利用MODFLOW，以抽水井为中心取一14955×14955m2的区域作为研究对象，剖分为121×121个网格。抽水井处网格边长5m，往外围逐渐加大至10，20，30，…，300，500m。为了便于比较，收敛指标均定为0.001m，选用不同的子程序包，可以分别求出对应的结果。数值计算结果与Theis公式［6］求得的解析解比较如图2所示。

由图2可知，SIP法和PCG法精度高，而SOR方法求得的结果较理论值明显偏小，具体原因还不十分清楚，可能与MODFLOW采用SOR方法是一种片状超松弛迭代法（SSOR）有关。该方法将模型的网格分割为一系列“片”（Slices），并将计算单元按所在的片进行组合，求解过程按片进行，对每一片所包含的计算单元利用Gauss消元法直接求解，每次求得的结果为各点的水头变化，所有的片都计算完就称为一次迭代。一次迭代完以后，MODFLOW将最大的水头变化值与收敛指标相比较，如果本次迭代得到的水头变化值仍大于收敛指标，则进行下一次迭代；反之则结束迭代过程。这种方法虽然使用了直接求解的Gauss消元法，但总的求解仍具有迭代的性质，而每次使用Gauss消元法都是在近似解的基础上进行的，因此有可能产生较大的误差。随着应用数学学科的发展，再加上MODFLOW程序结构的易添加性，可以预见，MODFLOW的求解子程序包必将更加多样化，应用范围也更为广泛。大量实际工作表明，只要恰当使用，MODFLOW不但用于模拟孔隙介质地下水的运动，而且可以用来解决裂隙介质中的地下水流动问题。不仅如此，经过合理的概化，MODFLOW还可以用来解决空气在土壤中的流动（Guo，1995）；将MODFLOW与溶质运移模拟的软件结合起来，还可以用来模拟诸如海水入侵等地下水密度为变量的问题（Guo和Benett，1997）［7］。

4 结论通过以上对MODFLOW的剖析，我们知道，MODFLOW程序结构合理，易于理解，便于操作，是一种较为权威的地下水流数值模拟软件，具有广泛的使用价值。它的权威性说到底是由其实用性所决定的。由此，笔者认为，我国的水文地质工作者应该消除“门户之见”，联合其他相关学科的专业技术人员，以实用性作为指导原则开发出自己的有关地下水流数值模拟的权威软件。而以上对MODFLOW有关缺陷的讨论，也正是我们在操作MODFLOW和设计软件时所应该注意的问题。可以说，MODFLOW为我们今后设计地下水流数值模拟软件，缩小同发达国家间的差距指明了正确的发展方向。美国Missimer国际软件公司的郭卫星博士为本文的撰写提供了数据资料，在此表示感谢！

图1 含水层在空间平面上的矩形离散示意

参考文献

［1］ McDonaldG．MichaelandArlenW．Harbaugh．Amodularthree-dimensionalfinite-differenceground-waterflowmodel，UnitedStatesGovernmentPrintingOffice，Washington，1988．

［2］ 吴剑锋．对地下水资源评价与管理现状有关问题的探讨．水文地质工程地质，1997，（5）．

［3］ 张宏仁．偏微分方程数值法在地学应用中的对比分析．地质学报1993，（3）．

［4］ 张宏仁．解渗流问题数值方法对比．水文地质工程地质，1984，（4）．

［5］ 张宏仁．有限单元法的水头反常问题．水文地质工程地质，1992，（5）．

［6］ 薛禹群．地下水动力学原理．北京：地质出版社，1986．

［7］ 郭卫星，卢国平编译．MODFLOW三维有限差分地下水流模型，南京大学地球科学系，1998．

西北黑河额济纳盆地水资源管理研究———三维地下水流数值模拟（武选民1）

　额济纳盆地地处我国西北黑河流域的下游．近年来由于黑河中游盆地用水量的不断增加流入额济纳盆地的水量减少使得额济纳盆地水文地质环境发生了极大变化泉水枯竭湖泊干枯沙尘暴频发绿洲萎缩植被生态环境严重恶化［1］ ．因此如何科学调配流域有限的水资源长期以来一直是人们普遍关注的热点问题之一．黑河中游向额济纳盆地的输入水量决定着地下水位的变化趋势潜水水位的埋藏深度及其变化直接控制着绿洲植被生态环境的演化．盆地内若能不断得到上游充足的水量输入潜水位的埋藏深度就会向浅埋的方向发展对于盆地内植被的良性生长演化和抑制生态环境的恶化无疑是有利的．但是大量的河水输入下游盆地一是会加剧中游地区用水的矛盾二是地下水位埋藏较浅时会增大地下水蒸发的无效消耗量造成宝贵水资源的浪费．因此确定合理的输水量使其既能满足绿洲正常生长演化用水量的需要又不致造成大量的无效蒸发损失对于缓解河流中下游用水矛盾抑制下游盆地绿洲植被生长演化状态恶化促进下游盆地生态环境向良性发展是十分必要的．本文通过额济纳盆地地下水系统的数值模拟对在不同输水量情况下额济纳盆地区域地下水位的变化趋势做出预测并提出合理的输水方案试图为额济纳绿洲生态环境的保护和建设提供科学依据．1　水文地质概况额济纳盆地南、西、北三面为低山所环抱东接巴丹吉林沙漠孔隙含水系统总面积约3．4万km2．据额济纳旗气象站1957－1999年的气象资料多年平均降水量38．2mm一次降雨在10mm 以上者甚少．黑河是流入额济纳盆地唯一的河流．黑河起源于祁连山流经甘肃省河西走廊的张掖、临泽盆地于高台县境内的正义峡穿过北山（ 合黎山） 流出走廊约180km 进入额济纳盆地．黑河在额济纳盆地内的总流长约240km．黑河从鼎新进入额济纳盆地流经约80km 于盆地内的狼心山分为东、西2个支流（ 河） 分别流向盆地北部的东、西居延海即黑河的尾阈．额济纳盆地为一相对完整的地下水流盆地含水层系统主要由第四系含水层（ 组） 、侏罗系碎屑岩类裂隙－孔隙含水层和基岩裂隙含水层所组成．侏罗系碎屑岩类裂隙孔隙水和基岩裂隙水呈条带状分

布于盆地周边，总面积约1500km2，水量贫乏．第四系含水层（ 组） 是盆地内含水层系统的主体，它是本次数值模拟的主要对象．第四系含水层（ 组） 包括潜水含水层、第一、二承压含水层以及第一和第二弱含水层［2，］ ．根据第四系含水层间的迭置关系，在平面上由南向北可将其划分为单一结构含水层区、双层结构含水层区和多层结构含水层区（ 图1） ．潜水含水层在盆地南部厚度可达到200m以上，一般地区厚度为30～50m，含水层底面埋深40～60m；第一承压含水层厚度一般为50～70m，顶、底面埋深分别为50～70m，110～150m；第二承压含水层厚度一般为40～60m，顶、底面埋深分别为160～170m，00～210m．　　额济纳盆地地下水系统的补给主要是黑河水季节性的垂向渗漏，其次是大气降水的入渗、相邻鼎新盆地的侧向径流、盆地东南部巴丹吉林沙漠潜水的侧向径流，以及外围基岩山区裂隙孔隙水及基岩裂隙水的侧向径流补给．地下水的排泄主要有：潜水的蒸发排泄、植被的蒸腾排泄及工农业生产和居民生活对地下水的开采．根据对盆地不同含水层地下水流场的对比分析可知［2，］ ，各含水层间的水力关系是：在盆地南部及周边地带，地下水从单一结构区的潜水含水层流向双层结构区的第一承压含水层；在盆地中部，从双层结构区的第一承压含水层流向多层结构区的第二承压含水层．地下水在不同的层中主要以水平径流的方式由补给区向盆地中心流动．在盆地中心水头的垂向分布为从第二承压含水层－第一承压含水层－潜水含水层，地下水头总体上依次由高变低，说明地下水存在着由下向上的越流，即第二承压含水层流向第一承压含水层，第一承压含水层流向其上的潜水含水层，最终在潜水含水层中以垂向蒸发蒸腾形式排泄于大气之中［3］ ．

1. 水文地质概念模型的建立

2．1　含水层系统的概化

为了实现对全盆地地下水流系统的整体模拟，充分利用盆地的自然边界，避免人为边界对模型计算造成的困难，正确反映各含水层（ 体） 相互间的联系，因此将盆地内不同成因类型的含水层视为统一含水系统，在模拟计算时用水文地质参数加以区别．

2．2　地下水系统边界条件的概化

（1） 底面边界：额济纳盆地第四系含水层的底面为侏罗系、第三系的泥岩及砂质泥岩，局部地段为震旦界大理岩，据钻孔资料，其透水性微弱，将其视为额济纳平原含水层系统的统一隔水底板．（2） 顶面边界：额济纳盆地地下水系统顶部边界为气－土界面，该界面与潜水面之间的非饱和带是联系大气降水、地表水与地下水的纽带．通过顶面边界进入和排出系统的包括：①黑河水的渗漏补给：根据黑河河道分布、河道的实际输水时间、输水量及不同河道段渗漏量的实测资料，将其概化为随时空变化的线性定量补给项．②灌溉水的回归入渗补给：灌溉水的入渗一是开采地下水灌溉农田时的入渗，二是在黑河河道输水期间引水灌溉草场、绿洲时的回归入渗补给．前者根据该区农作物的种植面积、灌溉制度和农灌开采的总水量等因素，依据同类地区灌溉水回归入渗系数，将其概化为季节性的局部面状补给项；对草场、绿洲灌溉水的回归入渗，则根据不同引灌地点、引水量、引灌时间等分区概化为季节性的面状补给项．③大气降水入渗补给：在十分干旱的额济纳地区，尽管一次降水对含水层饱和带的补给无意义，但一次降水的发生对抑制潜水的蒸发具有不可忽视的作用，因此将大气降水作为气－土界面上随时空变化的面状补给项．④潜水的蒸发排泄：蒸发排泄是指区内无植被覆盖的裸露地带的潜水蒸发排泄．据对该区所做的地下水均衡计算结果，潜水的蒸发量占地下水排泄总量的80％左右，植被蒸腾量占20％左右．可见潜水的蒸发是研究区内地下水排泄的主要项．根据地下水的埋深、潜水面以上包气带，的岩性、包气带水分运移规律等因素，将其概化为面，状排泄项．⑤植被的蒸腾消耗：额济纳盆地内绿洲植，被生长所需的水分主要依靠地下水的供给，因此植，被的蒸腾也是该区地下水的主要消耗项之一．根据，研究区绿洲的分布面积、植被根系发育带土壤结构、，植被的种类、植被在年内不同生长季节耗水量等要，素，将其概化为与潜水蒸发排泄不重合的、随季节变，化的面状排泄项．，（3） 侧向边界的概化：额济纳盆地第四系含水层，系统的四周均可视为第二类边界．其中东部、东南部，巴丹吉林沙漠区的地下水径流补给相对显著，作为，非零流量的第二类边界；西部、北部及东北部，第四，系含水层与山体间为山足面和断层接触，山区的侏，罗系碎屑岩类裂隙水和基岩裂隙水对平原区的补给，量甚微，将其作为零流量的第二类边界．，额济纳盆地与鼎新盆地接壤处的黑河河床，位，于模拟区的最上游，长度对整个盆地来说是一个很，小的尺度，为了保证模拟的仿真程度将其视为已知，水头边界．，2．3　源汇项的概化，源汇项主要是地下水的人工开采．开采地点主，要集中于平原区南部单一潜水含水层区和北部的绿，洲区．根据开采井位的分布和地下水的开采量作为，点状排泄项．，3　地下水流系统的数学模型及其解法根据以上对额济纳盆地地下水流系统所概化出的水文地质概念模型,将采用饱和－非饱和水流三维数学模型来描述其地下水流状态［4～8］ ,其数学模型的控制方程和定解条件可写为：控制方程式中：h 是压力水头（m） ；t 是时间（d） ；Ks 是饱和渗透系数张量；Kr 是相对渗透系数；z 是位置水头（m） ；q 是源汇项；F 是在饱和条件下为单位储水系数（1／m） ,在非饱和条件下为单位容水系数；θ是含水量；R 是计算区域；hi 是初始水头（m） ；h D 是狄里克雷边界BD 上的水头（m） ；qn 是底面边界Bn 上的单位面积流量,本模型中恒为零；qc 是侧向边界Bc上的单位面积流量；qp 是补给期间顶部边界Bv 上的综合补给强度；qe 是排泄期间顶部边界Bv 上的综合蒸发蒸腾强度．上述控制方程（1） 和定解条件（2） ～式构成了额济纳盆地地下水流的数学模型．4　数学模型的识别4．1　计算域的剖分对计算域采用了以六面体为基本单元进行剖分．垂直方向上剖分为5个六面体层,各六面体层由上而下依次代表饱和－非饱和潜水层、第一弱透水层、第一承压含水层、第二弱透水层、第二承压含水层．平面剖分是在额济纳盆地1∶50万综合水文地质图上进行的,单层剖分出四边形单元5514个,结点5738个,代表实际面积33987．5km2．因此全区共剖分出六面体总数27570个,结点总数34428个． 4．2　计算时限及步长根据实际获得的河流流量测定、地下水位观测孔水位变化、植被蒸腾强度等动态资料,模型识别的时限确定为1996年1月－1999年12月．考虑到计算工作量计算时段以月（30d） 为基本步长,在时间上的间隔共分为48个时段,总时间为1440d ． 4．3　初始压力水头的确定以不同含水层各个观测孔1996年1月的平均水位值为基础,应用趋势面加残差的方法,模拟出不同含水层各个剖分结点处的平均总水头值,进而得到相应结点的压力水头．地表属非饱和结点层,因此各结点的压力水头为负值． 4．4　顶面及侧向边界补给、排泄量通过计算区顶面的输出、输入项包括：河水渗漏补给、大气降水入渗、潜水蒸发和植被蒸腾．1996年1月1日至1999年12月31日黑河河水的渗漏补给总量为21．21亿m3,大气降水入渗年平均补给总量为1．25亿m3,潜水通过裸地区的年蒸发总量为5．72亿m3,绿洲区植被对潜水的年蒸腾消耗总量为2．9亿m3．全盆地地下水侧向径流年补给总量为1．4亿m3． 4．5　识别依据的选择根据实测获得的河流渗漏补给量、大气降水量、地下水位等相关的动态信息,模型识别的时间确定为1996年1月至1999年12月,以该时间段内实测的13个观测孔的地下水位动态序列作为拟合的目标． 4．6　非饱和土壤水分运动参数模拟区顶部为非饱和土壤带,为了获得非饱和土壤水分运动的相关参数,在盆地内不同岩性分布区采取样品进行了土壤水分特征曲线的测定［9］ ,其结果列于表1． 4．7　模型识别结果反演所确定的模型参数列于表2．各观测孔实测的与计算的历时水位对比示于图2．由各观测孔地下水位动态拟合曲线图可以看出,尽管计算曲线与实测曲线存在一定的偏差,但数值模拟计算的结果与实测值在总体上基本拟合,说明所建的数值模型能够较好地刻画额济纳盆地地下水流运移的基本规律．5　模型预报本次预报从1999年12月开始,以一个月为一个时段步长,预报总时长5年共60个时段． 5．1　黑河输水量的方案设计模型预测设计以下4种不同输水方案：方案一年输水总量10亿m3,类似于1996年度的实际输水总量；方案二年输水总量7．5亿m3,类似于1998年度的实际输水总量；方案三年输水总量6亿m3,类似于1996－1999年的实际输水总量的平均值；方案四年输水总量4．5亿m3,类似于1997年度的实际输水总量． 5．2　各输水方案的预测结果（1） 输水量为10亿m3／a 时地下水位变化趋势．在上游输水量为10亿m3／a 的条件下,盆地内地下水位普遍呈逐年上升的趋势．其中古日乃湖盆以西的冲洪积戈壁区地下水位以约1cm／a 的速率上升,古日乃湖盆中心水位上升速率约9cm／a ,湖盆东侧巴丹吉林沙漠边缘水位上升速率4．5cm／a ,进素土海子一带水位上升速率约3cm／a ．盆地西北隅狼心山以下东河上游一带地下水位上升速率约10cm／a ,东河中游约为6cm／a ,额济纳县城一带为8cm／a ,东居延海以南地区为14cm／a ．狼心山以下西河中游和赛汗陶莱为7cm／a ,赛汗陶莱以北地区为9．5cm／a ．（2） 输水量为7．5亿m3／a 时地下水位变化趋势．输水量为7．5亿m3／a 的情况下,古日乃湖盆以西的冲洪积戈壁区地下水位在预测的起始和末时刻基本相等；古日乃湖盆西侧计算末期相比初期水位下降了0．25m,但在预测的后三年当中水位没有出现持续下降的趋势；古日乃湖盆中心及东侧水位则有逐年上升的趋势,其速率分别约为6cm／a和3．5cm／a ；进素土海子一带水位出现逐年下降的趋势,下降速率为4．25cm／a ．盆地西北隅地下水位总体上处于稳中有升的状态,狼心山以下东河上中游及西河中游一带地下水位逐年上升,其上升速率为3～4．5cm／a ；额济纳及赛汗陶莱地区的水位则处于基本稳定的状态．（3） 输水量为6亿m3／a 时地下水位变化趋势．在输水量为6．0亿m3／a 的情况下,全区地下水位总体上处于缓慢持续下降的趋势．古日乃湖盆以西的冲洪积戈壁区水位下降速率为5．5cm／a ,古日乃湖盆西侧为4．8cm／a ,古日乃湖盆中心约为4cm／a ,湖盆东侧水位基本处于稳定状态,进素土海子一带水位下降速率为5．2cm／a ．盆地西北隅东河上中游地区水位接近平衡,额济纳县城一带水位的下降速率为3cm／a ,东居延海以南地区水位下降速率约为5cm／a ；西河中游一带水位下降速率为4cm／a ,向下游到赛汗陶莱水位下降速率增大为5～7cm／a ．（4） 输水量为4．5亿m3／a 时地下水位变化趋势．当上游向额济纳盆地的输水量减少到4．5亿m3／a 时,盆地内的地下水位将以较大的速率持续下降．在古日乃湖盆以西的冲洪积戈壁区地下水位下降速率约为8cm／a ,古日乃湖盆西侧为7．5cm／a ,湖盆东侧为3．5cm／a ,古日乃湖盆中心8．25cm／a ,进素土海子一带为8cm／a 以上．狼心山以下东河上游一带水位下降速率为4cm／a ,东河中游一带为6．5cm／a ,西河中游一带为7cm／a ,额济纳县城及赛汗陶莱地区水位的下降速率达到15cm／a 以上．由以上4种输水方案的预测结果可以看出,抑制额济纳盆地地下水位持续下降的合理输水总量应为7．5亿m3／a ,其中通过狼心山水文站的输水量不应小于5．5亿m3／a ．6　结论（1） 本次数值模拟是第一次在大范围区域进行的地下水流系统饱和－非饱和三维数值法计算,较为成功地尝试了对大型盆地进行数值模拟．（2） 从数值模拟结果可以看出,本次所建的额济纳盆地地下水流系统的概念模型是合理的,对各类边界条件的概化是恰当的,较好地刻画了额济纳盆地水文地质条件．（3） 为抑制额济纳盆地绿洲生态环境的进一步恶化,保证绿洲正常生长演化的生态用水量是必需的前提条件．模型预报的结果显示,黑河上中游向额济纳盆地合理的输水总量应为7．5亿m3／a ,其中通过狼心山水文站的水量应不小于5．5亿m3／a ．（4） 本次数值模拟是在大范围进行的,由于受工作精度和相关资料的限制,可能存在不足之处,今后有待进一步完善．

参考文献：
［1］ 陈梦熊．西北干旱区水资源开发与防止生态环境恶化［J］ ．水文地质工程地质199724（1） ：4－6．
　　CHEN M X．Water resource utilization and prevention of
eco-environment from deterioration in the arid region of
northwestern China ［ J］ ．Hydrogeology and Engineering
Geology 199724（1） ：4－6．

［2］ 武选民．西北黑河额济纳盆地地下水利用与生态环境保护研究［ D］ ．武汉：中国地质大学2000．　　WU X M．Groundwater exploitation and ecological environmentprotectionin Ejina basin of Hei River watershed northwestern China ［ D］ ．Wuhan ：China University of
Geosciences 2000．

［3］ 武选民史生胜黎志恒等．西北黑河下游额济纳盆地地下水系统研究（ 下）［J］ ．水文地质工程地质200229（2） ：30－33．　　WU X MSHI S S LI Z Het al ．Groundwater system in Ejina basin lower reach of Hei River northwestern China （the second part） ［ J］ ．Hydrogeology and Engineering Geology 200229（2） ：30－33．

［4］ Wang H F Anderson M P ．Introduction to groundwater modeling ［ M］ ．San Francisco ：W H Freeman and Company1982．

［5］ 李俊亭．地下水数值模拟［ M］ ．北京：地质出版社1988．
　　LI J T ．Numerical si mulation of groundwater ［ M］ ．Beijing
：Geological Publishing House 1988．
［6］ 陈崇希唐仲华．地下水流动问题数值方法［ M］ ．北京：
中国地质大学出版社1990．
　　CHEN C XTANG Z H．Problems in numerical si mulation
of groundwater flow ［ M］ ．Beijing ：China University
of Geosciences Press 1990．
［7］ Sukhija B S ．Groundwater recharge in semi-arid regions
of India ：An overview of result obtained using traces ［J］ ．
Hydrogeology Journal 19964（3） ：50－71．
［8］ Xu C YSingh V P ．A review of monthly water balance
models for water resources investigation ［ J］ ．Water Resource
Management 199812（1） ：31－50．
［9］ 许兆义．包气带水文地质专论［ M］ ．北京：地震出版社
1993．
　　XU Z Y．Unsaturated hydrogeology ［ M］ ．Beijing ：Seismological
Publishing House 1993．

VisualMODFLOW在地下水模拟中的应用———以河北省栾城县为例（贾金生）

近半个世纪以来水文地质学大致经历了两个发展阶段：从20世纪50年代到70年代中期是奠基阶段主要接受苏联思想的影响；从20世纪70年代中期到90年代是发展阶段随系统科学、环境科学、现代应用数学与计算机技术等新思想、新理论与新技术的输入使水文地质学的基本概念与研究范畴发生了巨大变革从定性研究步入了定量研究阶段纳入了系统工程的轨道［1］。随着工农业生产的发展和人民生活水平的提高地下水资源的供需矛盾日益突出这就对地下水资源的评价与管理提出了更高的要求即要从定量角度对地下水资源进行预测和评价建立合理开发利用方案。但水文地质条件本身固有的复杂性限制了用地下水动力学中建立的解析法解决问题的广泛性［2］。20世纪70年代中期以后数值模拟技术的应用使解决这一问题成为可能。地下水系统数学模型的建立是定量研究地下水运动规律的关键［3］。经过近30年的发展地下水模拟在我国经历了从无到有、从简单的水流模型到比较复杂的物质和热量运移模型从仿制到独立研制最后走向世界的艰难发展历程。现在我国已经建立了几乎囊括国际地下水模拟中心（IG-WMC）P．venderHeijde分类中所有的模型。研究范围涉及饱和带、非饱和带及饱和－非饱和带［4］。但是到目前为止我国还没有一个通用的、具有知识产权的权威计算软件这不能不说是一大遗憾［5］。近年来对于饱和带地下水流模拟的研究主要集中在三维流模型开发、流速场与流线的计算方法、非均质参数的区域概化和繁杂数据的优化处理。目前进行区域三维地下水流分析的主要软件有HST3D（Kipp1987）、SWIFT（Granwell和Reeves1981）以及世界上最流行的VisualMODFLOW（McDonaldandHarbangh1994）［6］。本文首先对VisualMODFLOW软件和MODFLOW基础理论作了简单介绍在此基础上应用VisualMODFLOW软件模拟了河北省栾城县的地下水流情况。通过对模型的校正与识别得出了比较满意的结果所建数学模型能够较好地反映当地的实际水文地质条件可以用来进行预测和对地下水资源的管理。1 VisualMODFLOW软件介绍MODFLOW（Themodularfinite－differencegroundwaterflowmodel）是由美国地质调查局（USGS）开发的用来模拟地下水流动和污染物迁移等特性的计算机程序（McDonaldandHarbaugh1988；HarbaughandMcDonald1996）。目前MODFLOW是全世界范围内模拟地下水流的应用程序。MODFLOW如此受欢迎归功于它具有以下特点：①MODFLOW所用的有限差分方法容易理解并且适用于许多现实条件。②MODFLOW可以用于一维、二维、准三维和三维模型。③数据输入格式、基本理论和每一个模块都经过了广泛验证。④模块化结构便于用户根据实际需要添加程序、完善功能和其他应用软件如Surfer、Excel等的结合。⑤MODFLOW模拟的结果可以用许多软件如Surfer、AutoCAD等显示和处理而且其自己的三维可视化结果处理也很便于用户理解和应用［78］。MODFLOW用来模拟的含水层系统应具有以下特点：饱和流状态、适合Darcy定律、地下水密度保持恒定以及水平水力传导率和导水系数的主流方向在整个含水层系统中保持不变。许多需要研究地下水流和污染物迁移的地下水含水层系统中都满足这些条件。MODFLOW可以模拟的水力过程及影响地下水系统的水力要素如图1所示。MODFLOW可以模拟潜水、承压水和隔水层中的稳定流与瞬变流的情况。许多影响因素和水文过程如河流、溪流、排水沟、泉水、水库、作物蒸散量、降雨和灌溉入渗补给等都可以用MODFLOW来模拟。MODFLOW提供了求解地下水流有限差分公式的很多种方法如强隐式迭代法SIP、逐次超松弛迭代法SOR、预调共轭梯度迭代法PCG2、SSOR等。用户可以根据自己研究的实际情况选择适合的有效求解方法。由于实际地质及水文地质条件的差异选择不同的求解程序包所得的结果是不一样的。几种求解方法可能都收敛也可能只收敛于一种或几种求解方法。MODFLOW的求解过程中引入了应力期（StressPeri-od）概念它将整个模拟时间分为若干个应力期每个应力期又可再分为若干个时段（TimeStep）。在同一个应力期各时间段既可以按等步长也可以按一个规定的几何序列逐渐增长。而在每个应力期内MOD-FLOW规定所有的外部源汇项的强度应保持不变。这样做不但简化规范了数据文件的输入而且使得物理概念更为明确。MODFLOW应用有限差分方法模拟含水层中的地下水流情况。求解过程中MODFLOW从空间上对含水层采用等距或不等距正交长方体剖分网格（图2）。这种网格的优点在于用户易于准备数据文件便于输入文件的规范化。其最大缺点在于增加了许多额外计算单元。但是随着计算机科学的迅猛发展计算机受网格数量的限制越来越小MODFLOW解决地下水流运动问题时可以将含水层剖分多达360×360×18个网格单元。VisualMODFLOW是目前国际上最盛行、且被各国同行一致认可的三维地下水流和溶质运移模拟评价的标准可视化专业软件系统。该系统是由加拿大WaterlooHydrogeologicInc．在MODFLOW软件基础上应用现代可视化技术开发研制的1994年8月首次在国际上公开发行。高度集成的软件包包括了用于地下水流模拟的MODFLOW、粒子运动轨迹和传播时间模拟的MODPATH、污染物在地下水中输移过程模拟的MT3D以及用于水文地质参数估计与优化的PEST并且具有直观的、强有力的图形交互界面。新颖的菜单结构便于用户对研究区离散及选择有效计算单元、确定边界条件与参数赋值、运行及校正模型以及用等值线或颜色阴影实现结果的可视化真正实现了人机对话。在模型的开发及结果显示过程中模型网格、输入参数和模拟结果都可以用剖面图或平面图显示。这个软件系统的最大特点是将数值模拟过程中的各个步骤天衣无缝地连接起来从开始建模、输入和修改各类水文地质参数与几何参数、运行模型、反演校正参数一直到显示输出结果使整个过程从头到尾系统化、规范化。这些特点是目前我国乃至世界上同类软件所不具备的。

2 模型的理论基础

MODFLOW是一个三维有限差分地下水流动模型它基于以下基本方程［9］：

式（1）中：Kxx、Kyy、Kzz为沿x、y、z坐标轴方向的水力传导率（LT－1）h是水头（L）W是在非平衡状态下通过均质、各向同性土壤介质单位体积的流量表示地下水的源和汇（T－1）Ss表示多孔介质的比贮水系数（L－1）t是时间（T）。

对三维稳定流动MODPATH的质量平衡方程可用有效孔隙率和渗流流速表示如式（2）：

式（2）中：Vx、Vy、Vz表示线性流动流速矢量在各坐标轴方向的分量（LT－1）n是含水层有效孔隙率（％）W表示由含水层内部单位体积源和汇产生的水量（T－1）。污染物运移模型MT3D的基本方程为：

式（3）中：C表示溶于水中的地下水污染物浓度（CL－1）t是时间（T）xi表示沿坐标轴各方向的距离（L）Dij表示水力扩散系数Vi表示地下水渗流速度（LT－1）qi表示源和汇的单位流量（T－1）Cs表示源和汇的浓度（CL－1）P表示含水层孔隙率（％）∑Rk表示化学反应项。

3 栾城县地下水数值模拟模型3．1 自然及社会经济概况栾城县位于河北省西南部北纬37°37′～38°东经114°30′～114°50′之间是太行山东麓与华北平原相交接的山前倾斜平原地区。地面标高在45～66m之间地形从西北缓缓向东南倾斜地面坡度介于1／500～1／1000之间。现状（1990年）栾城县境内只有氵交河和总退水渠———石家庄市废污水及雨水排放河道穿过县境。氵交河发源于鹿泉市南部的五峰山栾城县境内20．5km。氵交河上游多支流为季节性河流现在主要承接石家庄市的废污水及雨水。栾城县属暖温带半湿润地区大陆性气候明显温差较大多年平均气温12．2℃最冷月1月份平均气温－3．7℃最热月7月份平均气温26．0℃。≥10℃积温4251℃无霜期191d年日照时数2544h相对湿度65％多年平均年蒸发量1644．5mm年平均风速2．6m／s多年平均降雨量483．5mm（1949—2000年）且分布不均80％的降雨集中在7～9月份。栾城县属石家庄市管辖全县土地面积397km2。据1990年统计资料全县耕地面积31333．3hm2总人口328609人。1990年全县社会总产值69963万元（1980年不变价）工农业总产值60884万元（1980年不变价）其中工业总产值45622万元、农业总产值15262万元。改革开放后栾城县工业发展迅速主要行业分为建材、化工、食品、纺织、机械制造、电子设备、印刷业等。乡镇企业发展迅速主要集中在建材、造纸、化工、机械制造、纺织、食品及畜牧业等行业1990年乡镇企业完成产值33422万元（1980年不变价）。栾城县农业发达是我国著名的粮食产区。全县耕地土质良好且全部为水浇地主要靠机井灌溉沿氵交河两岸还利用一部分废污水实行井渠双灌。农业生产主要以种植业为主实行冬小麦—夏玉米一年两熟的轮作制度耕地复种指数高达1．68［10］。3．2 水文地质概况3．2．1 含水层主要特征 栾城县属太行山山前倾斜平原水文地质单元的一部分。由滹沱河、槐沙河冲洪积扇组成［10］。地下水主要赋存于第四系松散岩层孔隙中。第四系按岩性和沉积年代、含水层与隔水层的分布情况以及水动力条件等特征可划分为4个含水岩组。第1含水岩组地层上相当于全新统（Q4）第2含水岩组地层上相当于上更新统（Q3）第3含水岩组地层上相当于中更新统（Q2）第4含水岩组地层上相当于下更新统（Q1）。在水平方向上含水层由西部向东部单层厚度逐渐变厚粒度由细变粗层次由少增多富水性由弱变强。在垂直方向上上部及下部砂层粒度较细、厚度较小中部砂层粒度较粗含水层厚度较大。第1、2含水组发育较好两组之间没有连续的隔水层和弱透水层水力联系十分密切具有统一的地下水位因此可以视为一个含水层是目前的主要开采层段。该含水层以上更新统和全更新统的下段为主岩性为粗砂、中粗砂夹卵砾石。含水层底板埋深60～120m单井涌水量北部为50～70m3／（h·m）中部为30～50m3／（h·m）南部及西部边缘小于30m3／（h·m）栾城县钻孔统计资料如表1。3．2．2 地下水形成特征 栾城县地下水补给源主要有大气降水入渗、河流渗漏、灌溉回归以及侧向补给。大气降水入渗补给是栾城县地下水最主要的补给方式由于包气带岩性以砂和亚粘土为主入渗系数较大多在0．25左右地形坡度小为降水入渗提供了良好的条件。河渠渗漏补给主要来源于氵交河及总退水渠河道没有防渗设施漏水严重每年大约有0．2×108～0．5×108m3的水量入渗补给地下水。栾城县农田灌溉以地下水为主灌溉方式多为大水漫灌入渗量约占农田灌溉用水量的15％。侧向补给主要来源于西部元氏县及鹿泉市的山前侧向入渗多年平均值为0．3×108m3左右。栾城县地下水排泄以人工开采为主。农业灌溉用水占总用水量的88％左右且主要集中在冬小麦大量耗水的3～6月份。由于连年超采地下水埋深已达20m左右地下水蒸散量几乎为零。

3．2．3 地下水动态特征 栾城县地下水动态属降水、灌溉入渗补给———农业开采型地下水位变化主要受农业开采量及降雨量大小的控制。年初1～2月份地下水停采地下水位开始逐渐回升达最高水位。春灌开始后地下水位逐渐下降直至雨季开始农业停止开采地下水位开始逐渐回升如果秋冬季灌溉用水量较大地下水位也会暂时下降但总的趋势是回升直至第2年3月份达到最高值。目前栾城县地下水总体流向为自西北流向东南地下水水力坡度与地面坡度基本一致。3．2．4 水资源及开发利用现状 栾城县地表水资源缺乏目前只有石家庄市总退水渠和氵交河河段常年有明流通过。但是由于石家庄市废污水处理力度不够生物耗氧量、化学需氧量、悬浮物、硫化物、氨氮、氰化物、挥发性酚等污染物严重超标不具备利用价值只是在沿河两岸部分耕地利用废污水灌溉。栾城县地处太行山中段山前倾斜平原地带地下水资源十分丰富水质良好均适于工农业及生活利用。多年平均地下水资源量为0．81×108m3左右总补给量为1．0×108m3左右。近年来由于持续干旱和农业水利化程度的提高农业用水量逐年增加工业、乡镇企业的迅速发展以及人民生活水平的提高都对地下水资源造成了沉重负担。由于地下水连年超采引起地下水位的持续下降全县平均地下水埋深已由1949年的3～4m下降到1990年的15m左右。到1990年全县共有机井6754眼其中农村用井6703眼县级以上工业用井45眼城镇生活用井6眼。1990年全县平均降水量为583．8mm属于偏丰水年。全县开采地下水1．55×108m3左右其中农村供水系统1．44×108m3城镇自来水供水0．015×108m3县以上工业用水0．095×108m3。1990年全县地下水总补给量为1．05×108m3全县缺水量为0．5×108m3超采后使地下水位下降0．55m左右［10］。

3．3 水文地质概念模型3．3．1 地下水系统的基本分析 从以上的水文地质条件和当地地下水资源的利用情况分析现阶段栾城县地下水系统属于降水、灌溉入渗补给———农业开采型研究区基本上构成了地下水补、径、排的完整系统。其主要补给源有二：一是侧向补给即西北部山前侧向补给；二是垂向补给包括降雨入渗、灌溉回归及河渠入渗补给。地下水的排泄以农业开采为主其东南方向的侧向流出量很小。由于地下水埋深较大蒸发量接近于零。地下水流向为由西北流向东南水力坡度与地形坡度大体一致。而在地下水大规模开采前的天然状态时补给项主要为河流入渗、降雨入渗、侧渗及灌溉回归（地表水灌溉）排泄项主要为侧向流出及潜水蒸发。3．3．2 水文地质条件概化 计算区的南北边界与地下水位等值线近似直交虽有水量流进流出但流量很小视其为零流量边界；西北部边界侧向流入量较大且研究程度很高可以给出较为可靠的侧向流量数据因此作为第二类边界（流量边界）处理。东南边界近似与地下水位等值线平行水力坡度变化不大视其为变流量（依靠水位确定流量）边界［11］。计算区的底边界为一相对隔水的粘土层将其视为隔水边界计算区的包气带岩性为亚砂土及中细砂局部有亚粘土含水层接受降雨、灌溉及河渠的入渗补给。计算区地下水赋存于第四系上更新—全新统松散地层中含水层为非均质各向异性的潜水含水层主要开采层段第1、2含水层具有统一的地下水位属于混合开采按二维流考虑。3．4 数值模拟模型根据上述水文地质条件可以得到研究区的地下水数学模型如下［711］：

式中：h为潜水面高程（m）B为含水层底部高程（m）q（xyt）为边界单宽流量（m2／T）W为源汇项的代数和（m3）h0（xy）为初始水位（m）KxKy为分别为沿X轴、Y轴的含水层渗透系数（m／d）Sy为给水度（无量纲）Γ1为流量计算边界n为内法线t为时间。3．4．1 模型求解 求解过程采用WaterlooHydrogeologicInc．公司推出的VisualMODFLOW2．8．1（1999年）软件它以美国地质调查局（U．S．GeologicalSurvey）开发的MODFLOW（“ModularThree－dimensionalFinite－differenceGround－waterFlowModel”）为内核将MODFLOW、MODPATH、MT3D和PEST几个计算模块高度集成是当前地下水模拟软件中具有很高权威的专业软件。本次计算将第Ⅰ、Ⅱ含水组合并为一层计算区域划分为30行、30列共900个单元（图3）。其中：无效计算单元390个有效计算单元510个流量计算单元26个。本次计算调用的外应力子程序包为水井子程序所有补给项和排泄项都用注水井和抽水井的方式来表示。在分析降雨资料的前提下通过实地调查利用已有的观测资料确定河流来水量、农田灌溉的时间及灌水量最终确定注水井及抽水井的工作制度。模型计算过程则用WaterlooHydrogeologicInc．公司最新推荐的WHI子程序包。整个模拟过程分为12个应力期每个应力期为1个月时间步长为10d。模型运行时选择了MOD-FLOW和ZoneBudget2个选项对研究区域的地下水流和区域水平衡进行了模拟。图3 栾城县研究区网格剖分图Fig．3 DiscretizationmeshofresearchzoneinLuanchengcounty3．4．2 模型校正 利用建立的数学模型计算观测孔所在单元的水头并和实测的水头进行对比从而反求有关的水文地质参数。本次计算选用了5个观测孔都分布在第Ⅰ、Ⅱ含水层中。校正时间从1990年1月1日到1990年12月31日时间步长为30d共12个时段通过反演求参得出了含水层渗透系数和给水度的大小及分区结果（表2、图4）。5个观测孔的拟合情况与统计参数见图5、表3。

由图5和表3可以看出5个观测孔的平均水位误差最大值为0．14m正负误差比较均衡系统稳定性较好表明所建的数学模型、对水文地质条件的概化、边界条件的确定都与研究区域实际情况吻合较好。3．5 模型识别利用1990年的实测资料进行了模型校正求出了水力传导率与给水度的参数值与分区情况。但是经校正后的数学模型是否适于研究区地下水位与水量的近期及远景预测还必须通过进一步的实测资料来验证其可靠性。模型验证取1991年1月至1996年12月水位观测资料共分66个时段进行计算地下水开采量采用实际统计资料各项补给量与排泄项的计算方法与模型校正时一样。模型验证的可靠程度通过水位过程曲线拟合图（图6）来体现它反映了1991年1月至1996年12月梯度场的拟合情况。在选定的4个地下水位观测孔中从1991年1月至1996年12月模拟期中总共有265个模拟值在这些模拟值中计算误差大于0．5m的值有48个占总数的18．1％其中大于1m的值有3个占1．1％。这些误差主要产生在降雨偏少的1994年和降雨偏多的1996年表明在地下水位变化幅度较大时模型敏感度不够。在4个地下水位观测孔中误差最小的是0604＃观测井它位于栾城县中部受河流入渗及边界条件影响较少模拟结果稳定。误差较大的是0606－2＃与0607－1＃观测井两井都距离退水渠与氵交河较近受河流入渗影响剧烈模拟结果较差。但总的看4个观测孔的总体拟合程度较高。在1990年至1996年期间降雨类型包括了丰水年（1996年）降雨量为923．8mm平水年（1991年）降雨量为485．5mm枯水年（1994年）降雨量为260．5mm。模型的模拟结果表明模型是可靠的可以用于研究区地下水位与水量的近期及远景预测。

4 结论与讨论通过本研究结果可以获得了以下几点认识。目前我国在地下水数值模拟领域与其他先进国家相比还是有差距的。我国科技工作者应该加紧努力争取早日推出具有知识产权的地下水数值模拟权威软件。对水文地质条件的充分认识是地下水模拟工作的基础和关键所在。模型建立者应该对研究区有全面的了解包括水文地质状况、水资源开发利用历史与现状、工农业生产状况、农业种植结构、水力工程的分布与利用等等。在山前平原地区第1、2含水组具有密切水力联系、属于混合开采的区域可以将二层合并为一层考虑。建立模型过程中首先要力求地面高程与含水层底板高程输入的精度与准确度这是计算机计算的基础。在准备输入模型的各项补给量与排泄量时一定要紧扣当地实际情况这是反演求参的前提。只有在这些工作取得满意精度的情况下所求得的参数才能切实反映研究区含水层的实际情况校正与识别后的模型才可以用于研究区地下水资源的管理。参考文献：［1］ 陈梦熊．中国水文地质环境地质问题研究［M］．北京：地震出版社1998．［2］ 王文科李俊亭．地下水流数值模拟的发展与展望［J］．西北地质199516（4）：52－56．［3］ 李颖．水资源系统中地下水系统数学模型研究［J］．世界地质199817（2）：41－46．［4］ 薛禹群吴吉春．面临21世纪的中国地下水模拟问题［J］．水文地质工程地质1999（5）：1－3．［5］ 吴剑峰朱学愚．由MODFLOW浅谈地下水流数值模拟软件的发展趋势［J］．工程勘察2000（2）：12－15．［6］ 魏林宏束龙仓郝振纯．地下水流数值模拟的研究现状与发展趋势［J］．重庆大学学报（自然科学版）2000（10）：50－52．［7］ 武强董东林．水资源评价的可视化专业软件（VisualMODFLOW）与应用潜力［J］．水文地质工程地质1999（5）：21－23．［8］ 周念清朱蓉朱学愚．MODFLOW在宿迁市地下水资源评价中的应用［J］．水文地质工程地质2000（6）：9－13．［9］ 何杉．ProcessingMODFLOW软件在地下水污染防治中的应用［J］．199957（3）：16－18．［10］ 栾城县水政水资源综合管理办公室．石家庄市栾城县水资源开发利用现状调查报告［R］∙1993．［11］ 唐益群叶为民．地下水资源概论［M］．上海：同济大学出版社1998．