地下水流数值模型设计与应用

1、 地下水流数值模型设计与应用陈崇希 中国地质大学(武汉)地下水流数值模拟大纲前言第1章 地下水流定解问题概述第2章 地下水流有限差分法第3章 地下水流有限元法（限于迦辽金有限元法）第4章 边界条件 源汇 等特殊问题的模拟与处理 4.1 混合抽水井的模拟 4.2 自流井的模拟 4.3 水平井的模拟 4.4 混合观测孔水位的模拟 4.5 抽水井水位的校正方法 4.6 井孔表皮效应（井周扰动效应）的处理 4.7 面井内抽水井水位降深的计算 4.8 非完整抽水井附近观测孔井水位的校正 4.9 泉的模拟4.10 非完整河流与地下水间补给、排泄的模拟4.11 大气降雨、河渠、湖库入渗补给滞后性的刻画4.1

2、2 潜水蒸发排泄的处理4.13 初始水头条件的模拟4.14 含水系统海边界问题4.15 数值模型中流速和流量的计算问题4.16 岩溶管道-裂隙-孔隙三重介质地下水流模拟4.17 地下水开采-地面沉降模拟第5章 地下水流模型的拓宽 5.1 地下水饱和-非饱和流耦合模型 5.2 地下水-地表水流耦合模型 5.3 分布式水文模型简介第6章 反求水文地质参数的数值方法第7章 数值模型设计应用及水文地质勘查 6.1 数值模拟设计应用主要步骤 6.2 水文地质概念模型的设计 6.3 数值模型的建立 6.4 预测模型若干问题 6.5 关于基岩含水系统建模的特殊问题 6.6 数值模型水文地质勘查等若干问题第8

3、章 数值模拟实例 附件:三维流多边形（棱柱体）网格 有限差分方程的建立 我们已经学习过基于矩形网格的差分方法。不难看出，无论是等格距还是变格距矩形网格，这种剖分是有局限性，通常不能满足实际问题的需要。常见的实际问题大多是：含水层渗流区边界形状不规则，非均质参数分区的界线不规则，抽水井群和观测孔群的布置一般也是任意的。这些问题若用矩形网格剖分则很不方便,特别是井孔及内部与地下水有水力联系的河流等源汇点/线。若将所有井孔及其它源汇点都与格点一致,精度较高,但许多网格没有必要地被加密；若网格密度合理,则可能井孔及其它源汇点不能与格点一致而丢失精度。为此，我们介绍一种以辅助三角形剖分为基础的不规则网格

4、差分法，通常被称之为(任意)多边形网格法。 目前看来,这种方法是能使边界、内部界线及井孔等源汇点/线与网络格点、格线最佳逼近,而又尽量少增加格点数目的网络格式,是一种比较实用的方法。因此,本教材特别推荐此法(PGMS)。 Istok(1989)认为:“有限差分法难以处理不规则边界、非均质和各向异性介质或倾斜含水层”。把有限差分法限于矩形网格是不妥的。 多边形网格有限差分法其差分方程的建立,我们采用直接根据达西定律和水均衡原理建立差分方程。形成多边形棱柱均衡网格的方法ijkpqacbaboyxijkijkijkABC图2-6-1多边形网格均衡系统示意图2001多边形棱柱均衡网格的差分方程的建立

5、对于第m层格点 i 为中心的网格 ，它在平面上的投影区如图2-6-1B，根据达西定律和水均衡原理建立m层格点 i 的差分方程。 首先求出单位时间内通过 的周边流入 内的水量。由图2-6-1B可见，可以先求出通过三角形 ijk 内的 pb 边和 qb 边流入内的水量，然后对格点周围各三角形做类似的计算，并求和，则可得到通过的周边流入内的水量。由达西定律有： （2-6-1）式中： 根据 （2-6-1）式的计算方法，对m层格点周围的所有三角形做类似的计算，并求和得另外，容易求得通过均衡体顶底面流入的水量为（2-6-3） （2-6-2） 若记Qwi为其源汇项，包括开采井、泉流量、河流的补给、排泄、降雨

6、入渗等，则方程可以写为，（2-6-4） 这就是多边形网格地下水三维流方程。 引 言 地下水资源评价、管理的核心是地下水动态的预测，当前地下水动态预测最有效的手段，甚至可以说唯一的手段，是数值模拟。国际上流行的地下水模拟专业软件很多，如Visual MODFLOW、FEFLOW等，它们在对诸水文地质模拟要素的刻画上大同小异。当前，Visual MODFLOW在国内外被普遍采用，它在可视化方面做得比较精緻*，然而该软件在地下水流的模拟上尚不尽人意，特别是对于普遍存在的混合抽水井流量分配、混合观测孔的水位、自流井的流量动态及降雨渠道等入渗补给地下水的滞后性等模拟要素，或者缺少模拟功能或者存在明显的缺

7、陷， 在泉、地下水-河流补给/排泄等方面也存在待改进部分； 当然，MODFLOW采用的矩形网格剖分方法也要损失模拟精度和/或增大运算量。 基于上述认识，于上世纪90年代申请国家自然科学基金多项以研究上述主要问题的流动机理及模拟方法*。做好基础准备后,于20032005年立项研制地下水软件。中国地调局和水文地质环境地质研究所要求“软件的核心部分对水文地质条件和井孔的刻画要求优于MODFLOW软件（地下水三维流矩形有限差分模拟软件）” 我们深刻体会到 ：数值模拟的核心是“防止模拟失真，提高仿真性”。因此，努力分析流动机理并用于数值模拟。这些成果已经集成为一个基于多边形网格的三维地下水流有限差分模拟

8、系统（简称PGMS，即Polygon-grid finite-difference groundwater modeling system）。 PGMS软件选择多边形网格有限差分数值模拟，就 剖分 而言，它比MODFLOW 矩形网格有限差分更能适应复杂的边界条件，容易地将水文地质诸模拟要素置于格点处和格线逼近水文地质界线；就精度而言，与常用的有限元相当，且物理意义明确。 PGMS软件基本模型已与一、二、三维地下水流解析解做过模拟对比，模拟结果良好。 PGMS软件已在近几年地质调查项目（新疆渭干河流域, 陕西渭北岩溶水, 甘肃河西走廊疏勒河流域 及 黑河中游地区 等项目)地下水模型中得到全面的应

9、用与检验，个别水文地质模拟要素，如混合抽水井和混合观测孔和降雨渠道等入渗补给滞后性等，则被应用更多次。 尽管现今的PGMS软件属1.0版,但已具有很好的应用前景。我们将进一步发展、完善该软件。致谢 PGMS模型的研制获得了中国地质调查局和水文地质环境地质研究所地质大调查项目及多项国家自然科学基金项目的支持; 广西地矿局科技处和北海地矿公司专门设置混合抽水试验场为本项目提供现场试验及滨海区地下水潮汐效应的系统观測，中国地质大学（武汉）渗流实验室为模型的验证试验提供了实验室条件，甘肃地调院、新疆一水、甘肃二水,陕西二水和山东烟台地质队以及水利系统等单位为软件的应用提供了野外现场的系统资料，国内外专

10、家学者对PGMS的开发给予了热情的鼓励和支持，在此一并表示感谢！4.1 混合抽水井的模拟 4.1.1 概述 什么是混合井孔?依原来的定义,贯穿两个或两个以上含水层的井孔称为混合井孔(图4-1-1)。 这种传统上的定义实际上隐含着一个假定，即含水层内地下水是水平二维流动。 4.1 混合抽水井的模拟 4.1.1 概述 图4-1-1 混合井孔的传统定义4-1-1 混合抽水井、混合观测孔的原定义(陈崇希等,1992) 我们认为:三维流场中，常规的（非点状滤管的）抽水井或观测孔，不管其置于多层结构还是均质单层结构的含水系统，都属于混合井孔(图4-1-2缺) 。这是由于三维流场中,井孔滤管中不同深度的水头

11、是不相等的，因此井管中的水要发生垂直流动，即使在不抽水的观测孔中也一样。这是混合井孔的本质所在(陈崇希等,1998b)。 由于自然界地下水普遍存在三维流(图4-1-3)，而井孔，特别是抽水井几乎都是非点状（滤管）井，因此混合抽水井普遍存在。 图4-1-3 自然界地下水普遍的三维流模式(陈崇希等,2005b) 开采地下水的井孔多采用混合井（除非水质不符合要求），以简化成井工艺，增加出水量和降低成本。因此，几乎没有无混合井的地区，特别是民间可以凿井的我国（包括台湾大量的是混合井）。既使现在开始控制打井，控制开采，但以往的混合抽水井转变为混合“观测孔”。 混合井在地下水流模型中如何刻画，国外尚未很好

12、解决。目前主要采用由美国地调局推出，并在国际上广泛流传和应用的三维有限差分地下水流模型MODFLOW软件处理。 MODFLOW建议：“多层井的流量必须以某种方式人为地分配给每一单层，把井流量按每层的导水系数大小分配，即 Qi /Qw=Ti /T ” (1988、1996、2000)。 为了便于讨论，又不失其一般，我们以贯穿两个含水层的混合抽水井为例进行讨论。如此，上式可表述为 Q1 /Q2=T1 /T2 MODFLOW对混合井的这种处理方法没有给以理论上任何的分析、说明，缺乏理论依据，应用中与实际也不符(陈崇希等, 1992)。这是因为 (陈崇希等,2003a,2007,2008;黎明等,20

13、08)： （1）含水层导水系数对井孔流量的影响，不会如此简单。例如,混合井附近岩性(渗透系数)发生变化，甚至混合井打在岩性透镜体上，怎样影响井管的流量分配？含水层厚度发生变化又如何改变流量的分配？ K3K1K2 混合井管贯穿岩性透镜体图4-1-3 混合井管贯穿岩性透镜体上 （2）含水层的参数影响混合井流量的分配，导水系数只是其中一个因素，含水层的弹性给水度（储水系数）就不起作用了吗? （3）井管的流量分配不仅与含水层的水文地质参数分布有关，还与外边界条件、井径、水泵吸水口的位置及其它抽水井的干扰等有关， （4）按MODFLOW的“流量预分法” ，一个混合抽水井的Q1/Q2 始终是个常量(因为T

14、1/T2为常量)。然而在模拟过程中，特别是预测时，该混合抽水井附近可加入或关闭 抽水井，在这种井群干扰下，原混合井的流量比还会保持常量吗？ （5）混合观测孔是混合抽水井的特殊情况（Qw =0）(图4-1-4)，对于两层混合的观测孔，其孔中水位(混合水位)必界于两含水层水位之间，即混合观测孔对于其中一含水层(例如2含水层)起抽水作用（Q20）, 对于另一含水层（1含水层）起注水作用（Q10）。如此，Q1/Q20。如此，两个比值怎能相等？！ 图4-1-4 混合观测孔中的各分层的流量及漏斗曲线 笔者将包括上述分析、批评MODFLOW的稿件投到美国WRR杂志,其某审稿者云：“MODFLOW并没有说上述

15、比例关系可用于混合观测孔”而否定此稿件。实际上，当抽水井的流量Qw足够小时，不会改变Q20和Q10的情况,仍存在比值Q1/Q2为负值的情况，即出现负值(Q1/Q2)等于正值(T1 /T2)的荒谬。 从理论上说，一种正确的方法，当蜕化为其特殊条件时，也必定是正确的。科学研究中，人们正是通过许多特殊来认识一般的；反过来，通过特殊来鉴别某理论的正确与否。这是常用的科学的方法,怎么WRR杂志审稿者,?! （6）按各层导水系数T的比例来预先人为划分各层的流量是缺乏理论依据的，因为这种做法要求各分层的有效井径相等和井壁处的水力坡度上下处处相等，这两个条件不可能人为控制，预先也不得而知。 基于上述分析，用导

16、水系数的比值分配各层的流量是不妥的。实质上，这种方法不是模拟，而是“处理”，一种与机理不符的“处理”。 的确，Javandkl and Witherspoon(1969) 在多层井地下水流数值模拟研究中认为：当流动稳定时,在所研究的条件下,沿井壁的水力梯度是相等的。由此可得出,混合井各层的进水流量与其导水系数成正比。 但哪篇论文的井流基本模型是：和(实质上)的圆岛模型,它的结论不能随意地拓宽到一般条件。 笔者最近获得的几个混合井流解析模型表明(陈崇希等,2011; 陈崇希,2012)：对于傍河稳定井流；只有当各层分别为均质、等厚含水层,外边界条件相同，即两层抽水井至河流的距离相等（ ）,及内边

17、界条件相同，即两层抽水井的有效井径（非公称井径）相等（ ），并允许忽略井中水流的水头损失时，各分层流量的分配才与导水系数成正比。 对于不稳定井流；则还要求两含水层的压力传导系数相等 ，且在混合井凿成（尚未抽水）足够长的时间后 (陈崇希等,2011; 陈崇希,2012)。 如此多的假定条件的要求，对于实际的水文地质条件是不可能满足的。因此，数值模型中采用“流量预分法”处理混合抽水井是不对的。北海抽水试验场的研究(陈崇希等,1992; 陈崇希等,1998)也表明其错误。 另一处理混合抽水井的方法是将混合井两含水层之间的井筒给予一个很大的垂向渗透系数Kz（模拟井筒效应），通过均衡计算确定各层流量(W

18、en-Hsing Chiang,2001），垂向渗透系数可表示为(Neville等,2004) 式中 和 分别是水的密度和动力粘滞系数; 是重力加速度; 井管内径。大垂向渗透系数法比“流量预分法”有所改进,但用圆管层流公式刻画抽水井井管的垂向渗透系数仍有很大的局限性(将在4.3中的物理模型一节中分析)。 “防止模拟失真，提高仿真性”是水文地质模拟工作者的核心任务( 陈崇希,2003)，而地下水流系统中普遍存在的混合井，是当前国内外模拟失真的主要问题之一，应当引起我们的足够重视。 混合抽水井模拟在水质模型中的重要意义*: 上述问题促使我们去思考地下水动力学的井流分析中引入热传导源汇理论的适用性。

19、传统地下水流模型都是纯渗流模型，因此以井孔的滤管井壁为其边界( 因为水进入滤管井壁后已不是渗流)，为此必需给出滤管井壁的边界条件流量分布或水头分布。这是出现上述问题的根源。 4.1.2 混合抽水井模拟方法“渗流-管流耦合模型”的提出 我们的思路是(1992)：为了避免难以给定的滤管井壁的边界条件，将井管内的水流与井周围的地下水流一起作为模拟对象，从而把模型的边界从滤管井壁移至井口(严格地讲是水泵进水口), 采用真实抽水流量或井中水位；滤管井壁的流量分布和水头分布是模拟的结果，而无需人为预先给定。所建的新模型称为“渗流-管流耦合模型” (1992)，其中渗流刻画地下水流，管流刻画井管中的水流。

20、对于“渗流-管流耦合模型”的求解，提出两种方法： “第一种是分别建立渗流模型和管流模型，然后采用迭代方法求解” (陈崇希等,1992)； 第二种是把管流部分视为渗透系数很大的圆柱形透镜体，并对管流引入形式上服从达西定律的“等效渗透系数”，从而将渗流-管流模型视为含有圆柱形“透镜体”的新的“渗流”系统(陈崇希等, 1992) 。对于这个新“渗流”系统，关键问题是如何确定井管的“等效渗透系数”。 下面介绍第二种方法的具体推导. 从流体力学圆管中水流的水头损失H(伯努里)方程出发 （4-1-2） 式中：H 是水头损失(m)；f是摩擦系数；l是管长(m)；d是管内直径(m)；u是管内平均流速(m/d)

21、；g是重力加速度(m/d2)。当管流为层流时 （4-1-3） 式中：Re是雷诺数。而 （4-1-4） （4-1-5）式中：是流体运动粘度；是流体动力粘度；是流体密度。将（4-1-3）至（4-1-5）式代入（4-1-2）式， （4-1-6） 式中：J是水力坡度；是流体重度。将上式写成渗透流速的形式，对于管流，其空隙率n=1，则渗透流速V=nu=u，故有 （4-1-7） 将此式与达西定律 对比，可得出圆管“层流状态等效渗透系数KL”的表达式： （4-1-8） 其中，KL定义为层流状态下管流的等效渗透系数，方程（4-1-8） 分别由陈崇希（1966）和J.Bear（1972）分别获得。当管流呈紊流状

22、态时，方程（4-1-2）可改写为(1992) （4-1-9）如果定义紊流状态下管流的等效渗透系数KN为 （4-1-10) 则 （4-1-11)该式也具有达西定律的形式。 本来，紊流的运动规律有别于达西定律，而且不同流态区具有不同的形式。引入等效渗透系数后，将5个流态分区（1个层流区和4个紊流区）的运动规律统一为达西定律形式，且与地下水渗流定律一致。即 V=KeJ （14） 其中（4-1-12） 为了求解水头分布，必须已知管流的水头损失H，依（4-1-2）式水头损失H又依赖于f 和u，f值又取决于Re(尼古拉池试验曲线 见图4-1-5 )，Re却又与u有关，即（4-1-4）式诸要素间互相依赖,因

23、此用迭代方法求解。为了保证必要的模拟精度，我们建议图4-1-6所示的模拟方法。图4-1-5 Nikuradse试验曲线 图4-1-6 混合井模拟流程图 本来,涉及混合抽水井的地下水流问题属于线性-非线性流问题,其定解问题的表述是此较复杂的,引入等效渗透系数Ke 后,只要将传统的纯地下水渗流问题的数学描述,将其中的渗透系数K用等效渗透系数Ke取代,就可以了,表述简单明辽。 4.1.3 北海混合抽水试验场 混合抽水试验确定分层水文地质参数实例(1989-1992) 上述方法首次用于广西北海孔隙含水系统混合抽水试验场(陈崇希等,1992; 陈崇希等,1998) . 试验场为三个含水层夹两个弱透水层的

24、含水系统.该试验场由间距为1.9m的双混合抽水井(和承压含水层)及8个观测孔(和承压含水层各3个,和承压含水层混合孔1个,潜水层1个)组成。观测孔至混合抽水井距离为6.9437.14m, 混合观测孔至混合抽水井距离为10.96m。试验时间为1988年12月4日至22日: (1) 410日为自然水位观测阶段,以掌握自然水位下降速率和潮汐效应变化,观测时间间隔为1h;(2)1116日为定流量抽水时间,两井流量之和为3394m3/d,延续时间4.5d; (3)1722日为停止抽水水位恢复阶段。 试验场为三个含水层夹两个弱透水层的含水系统，所建的为准三维流模型，并采用三角单元剖分的有限元数值方法。所有

25、抽水井和观测孔都设置结点，而观测孔的实测水位经过天然动态的校正后,再用以拟合。(详见4.13) 经过调试、优选，在不改变现场水文地质工程师提供的参数分区及参数值上下限的前提下进行拟合求参。最终的拟合动态曲线示于图4-1-7,拟合误差统计示于表4-1-1。表4-1-2 给出混合抽水井孔分层流量Q1、Q2的变化及流量以及导水系数之比的一览表,就此实例而言，MODFLOW流量分配法的误差达70%85%。这么大的误差是不可接受的。 尽管“渗流-管流耦合模型”已用于多个实例,但由于它的重要性,我们仍对其做了室内物理模拟及其对应的数值模拟研究(陈崇希等,2004b;Chen等,2003*) ,以进一步检验

26、该模型的仿真性。研究成果是肯定的。 图4-1-7 北海混合抽水试验 水头降深拟合曲线(据陈崇希等,1992;1998) 1实测水头降深曲线；2模拟水头降深曲线 表4-1-1 模拟水头与观测孔水位降深误差统计表 (据陈崇希 等,1992;1998)表4-1-2 混合抽水井孔分层流量Q1、Q2的变化及流量以及导水系数比的一览表(据陈崇希等,1998)黎明的博士论文第7章 数值模型设计、应用及水文地质勘查 大多数值模拟教材及一些商业软件，基本上限于模拟方法，而不涉及模型设计。实际上，数值模拟应用中，首先碰到的问题是模型设计。失误在模型设计上失误在数值模型应用起跑线上的数值模拟的案例并非个別,包括自认

27、为站在我国“数值模拟制高点” 的学者。下文将分析其案例*。 为此,本书增加模型设计一章，除了简述应用步骤外,主要针对普遍存在的问题提出笔者的见解,其中有些是建议性的，有些是讨论性的，有些还难以形成“要求”。本书采用较多的正反面实例都是宝贵的学术财富，来说明的问题，提供读者讨论。7.1 数值模型设计主要步骤 数值模拟方法求解地下水流问题通常分为两个大步骤:(1)反演问题，即模型识别（或模型校正）问题: 包括模型类型的检验、含水系统参数的确定和边界条件的校正等等； (2)正演问题，即预测问题: 是种种条件改变下（地下水开采的分布、强度的变化；水库的修建，渠系运行等地表水调度，灌溉制度的变化等水利工

28、程运转；矿山地下水的疏干排水）的地下水的水头、流量及水均衡诸要素等动态的预测。地下水动态预测是地下水资源-环境优化管理的基础 *。7.1.1 明确目标任务数值模拟的目标任务不同，其模型设计等也会有差异。当前使用较多的，例如：区域（盆地级，流域、大支流级）地下水资源评价；水源地级的地下水动态预测；地下水开采-地面沉降；地下水开采-海水入侵；地下水开采-生态环境保护；矿坑排水疏干及水利工程相关问题 等。7.1.2 水文地质资料搜集 广泛搜集水文地质资料，全面分析含水系统的结构、岩性，地下水的种种补给、排泄的方式与地下水的动态，并分析由此控制的区域地下水径流特征及动态规律。“要尽量搜集已有的水文地质

29、资料”，这是后续工作的基础。这一点似乎不会有异议，但不然，例如某区域性地面沉降模型，在众多的潜水观测孔中才搜集14个孔的水位，还以此作全区17000km2面积的时空插值，这怎么控制得了！且不说潜水面的分布远比承压水头复杂的多。数值模型设计者不应限于室内搜集资料，还应现场重点踏勘，充分倾听野外水文地质工作者的见解并与之交换意见。这一点很重要。7.1.3 水文地质模型的设计（在7.2节重点分析）地下水赋存、运动于复杂的地质体中，并与大气及地表水存在多种形式交换，因此自然界地下水流不可能、也没有必要完全按原型设计，即水文地质条件要概化为模型，称为水文地质(概念)模型；但模型又不能失真,即模型的运转要

30、能够基本再现地下水的动态。模型应用中，在水文地质模型设计之后，要梳理涉及的诸地下水的模拟要素，以便寻找合适的软件，或者明确自已编制的相关程序应包含哪些模拟要素。水文地质模型的设计是模型应用中最重要的环节。 7.1.4 数学模型的建立 当前涉及的地下水模拟任务，基本上都是井孔-含水系统问题，而这类问题的数学模型目前主要有两类：第一类是传统的纯渗流模型（例如MODFLOW、FEFLOW等软件），这种模型将抽/注水井视为线汇/线源，其流量/水头采用人工分配的方法，因而混合观测孔的水位数据不能利用，混合抽/注水井、自流井和水平井的处理（非模拟）不符合水流机理，若免强采用，可能导致明显的差误。 第二类是

31、 “渗流-管流耦合模型”属于Darcy-非Darcy流模型,该模型由陈崇希等（1992, 1998b)提出后,做了大量模型和引用,如实例和理想模型（陈崇希等，1995b,2003b,2003c,2003d,2004a,2004b,2005b；成建梅等,1998；万军伟，2002；刘文波,2003；胡立堂，2004；黎明，2005；韩巍, 2005；Chen C X等,1999,2003b；Cheng J M等,1998,2005；Hu Li-tang等,2011；Mohamed 等, 2006；wang wei等, 2007；地大 200? ；Simon 等2008;Seung-Hyun Ki

32、m 等,2008）及 PGMS软件(陈崇希等，2005a；2007)。 这类模型的优点：刻画井孔-含水系统克服了传统纯渗流模型以井壁为边界条件,需要人工给定流量分配/水头分布的缺陷；另外，这类模型容易发展为“岩溶管道-裂隙-孔隙三重介质模型”(陈崇希，1995a)。缺点：模拟量比较大，但当前的PC机完全可以承担模拟一般的实际地下水流问题。基于上述，本教材推荐后者。7.1.5 数值模型的建立(在7.3详细分析)由上述建立的水文地质概念模型和数学模型,选用第2、3章的基本数值方法及第4章诸相关地下水要素的模拟或处理方法，将时空离散化,选用合适的软件或自编程序,依据搜集到的水文地质数据,通过模型识别

33、,建立数值模型。此部分是实施模型设计,确立数值模型过程最重要的一步,也是最化时间的一个流程。具体內容、步骤及注意问题在7.3中详细分析。7.1.6 数值模拟预测 (在7.4详述)建立数值模型的主要目的之一是:确定可持续开采量及地下水资源-环境（优化）管理。地下水资源-环境优化管理的基础 是地下水动态预测。所谓“最优化模型”，实际上是地下水的预测模型与数学上（运筹学）最优化程序的耦合。因此，水文地质模拟者的任务是提供一个优质的预测模型 *。7.2 水文地质(概念)模型设计地质体及其中的地下水流十分复杂，因此所建立的模型不可能，也没有必要完全按原型设计，即水文地质条件要概化(概念)模型；但模型不能

34、失真,即模型的运转要能够基本再现地下水的动态。这是设计水文地质模型的两原则。目前的一种倾向是随意地简化模型，或由于不会模拟某些要素而被“简化”掉。分析研究区地下水流的基本特征，建立地下水流概念模型是模型设计的核心，往往也是困难的一步。要求设计者具较高的综合能力和高度的责任心。不要轻视对水文地质条件的分析，许多概念模型的错误就是由此引起。重视对水文地质条件的分析，这不是口头主张与否的问题,几乎所有的文献都如此主张,甚至十分强调对水文地质条件的分析，然而重要的是看其如何具体分析水文地质条件、流动特征而建立相应的仿真模型。根据搜集的水文地质资料，进行地下水流模型的初步设计，主要内容包括：模型类型的选

35、定，模型空间范围的划定，边界条件等的确定关于边界条件先说明一点：地下水流概念模型设计、建立过程中，不应只说明研究区存在第一、二类边界等，这仅仅停留在教科书上的分类，而应涉及研究区的具体条件的分析与设置。提出源汇、边界条件等诸水文地质模拟要素的模拟方法如果某商业软件或他人的软件能够满足研究区全部水文地质要素的模拟，且模拟方法无所改进时，也可直接使用，但应指明所采用的软件名称、版本和所有者等。当然,应该将数值模型中有关边界条件、源汇等诸要素采用哪一种模拟方法加以说明, 以表明使用的正确性。在上述基础上提出资料清单，如果必要和可能的话，还要做补充水文地质勘查设计。补充水文地质勘查工作完成之后，根据新

36、信息再修订模型设计。7.2.1 模型基本类型的论证地下水流模型的基本类型通常是指模型的维数，各向同性/各向异性，稳定流/不稳定流,当然,Darcy流/非Darcy流/Darcy-非Darcy流也属于模型的基本类型。“当前数值模拟的核心问题是防止模拟失真，努力提高仿真性。”“模拟失真是多方面的，其中最主要的是概念模型失真。一旦水文地质概念模型不能刻画研究区地下水的基本运动规律，其后的工作再精细，也不可能弥补建立概念模型的失误*。因此，要十分重视水文地质概念模型的建立”。 (陈崇希,1999c,2002a,2003a) 上世纪70年代中后期之前,我国地下水流计算主要是解析法。哪时笔者有一句“经典的

37、话” “后期对数据的精心分析也难以弥补设计上的错误”。例如对下图所示水文地质条件(夹角为120度的两隔水边界形成的扇形含水层)做抽水试验设计，以求取水文地质参数。“后期对数据的精心分析也难以弥补设计上的错误”。观测孔的设置也有讲究。“然而，这个问题并没有引起水文地质数值模拟工作者的普遍重视。我们要避免对水文地质条件、地下水流动基本特征未加精细分析就摆出模型。这种模型的预测结果是难以可信的，”尽管有的模型其短时间的拟合度还不是很差。然而“拟合好或差不是检验建模合理性的唯一标准，更重要的也是最根本的是视水文地质实体的概化是否合理，模拟方法是否符合机理，正因为如此，我们十分强调对水文地质条件的正确分

38、析和建模论证。”“建模是需要论证的。”“ 在建模论证中最主要手段是流网分析。未恰当建立概念模型的实例，不是个别现象” (陈崇希等,1999c,2002a,2003a) 7.2.1.1 建模维数问题及其不同意见实例群1 一个地区建立数值模型，首先的一个问题是建立几维的模型。这是个有争议的问题。上世纪70年代来，直到2003-2005年我国新一轮地下水资源及其环境问题调查建立地下水模型之前，我国所做的地下水数值模型估计有数百个，大部分是二维模型，部分是准三维流模型，很少量是三维流模型。1. 建立盆地-流域级模型的维数问题存在种种不同见解：“地下水资源评价问题用二维流区域模型就足够了”；“建三维模型

39、资料要求高, 难以满足，建准三维模型足矣”“区域地下水流模型，平面尺度为几十至几百公里计，而厚度仅几十至几百米，差两三个数量级，没有必要做三维模型”；“资料充分时做三维模型；资料不足时,可做二维模型或准三维模型”；等等。2. 笔者认为*： 实际区域地下水流基本上都属于三维流：由于补给(降雨、河渠、灌溉等入渗)、排泄(蒸发蒸腾、泉及非完整河等)都发生在潜水面上或其下浅部，而含水系统的厚度大多为百米级，这是区域地下水形成三维流的基本原因（见图4-1-3）。特别是“重点目标区”往往在冲洪积扇的前缘和冲积平原区，哪里的地下水流速的垂直分量更是不可忽略。（如黑河中游和下游地区，疏勒河流域，渭干河流域，渭

40、北岩溶水地区，晋陕蒙天桥泉域，等）。另外，地下水开采、矿坑排水及水库、渠系等水利工程的运作，大多起到加强三维流的作用。在这些条件下，二维流模型难以满足研究目标的需要。 图4-1-3 自然界地下水普遍的三维流模式(陈崇希等,2005b) 三维流问题，只要含水层和弱透水层的渗透系数相差两个数量级以上，即可用准三维模型刻画吗？Neuman 等（1969）研究认为，当含水层和弱透水层的渗透系数相差两个数量级或更大时，此时忽略弱透水层内水平流速分量而导致目的层内水头误差小于5%，即三维流问题可以近似采用准三维模型刻画。对于实际问题，含水层与弱透水层的渗透系数之比往往达三个数量级以上，另外，准三维流模型对

41、水文地质资料要求、数据整理和模型运行工作量等诸方面都比三维流模型简单、省时得多，基于这两个原因，许多国内外地下水模拟专家主张采用准三维流模型。然而，我们注意到:自然界弱透水层中较普遍地呈现“千层饼状”的薄互层岩性，存在“米级”，“厘米级”乃至更簿的互层状沉积物，或者含水系统概化时，由于模型分层数远小于剖面上砂性土/粘性土的层数而概化的“弱透水层”中已夹有砂性土，它们在模型的“层”中表现出各向异性的特征。那么，这种情况下准三维流模型的相对误差还是小于5%吗？笔者等通过数值模拟研究后得出结论（陈崇希等2001）： 含水层和弱透水层的渗透系数相差两个数量级或更大的条件，三维流问题可以近似采用准三维模

42、型的误差小于5%的结论，仅适用于弱透水层为各向同性介质； 对于弱透水层的各向异性比为10时(实际条件通常明显地大于此值)，误差已超过27.7%，而且随着弱透水层单位储水系数 的减小和模拟时间的延长，误差还要增大，已存在超过30.5%的情况。”（陈崇希等2001）基于此，普遍的情况不宜概化为准三维流模型。 二维模型能正确地确定可持续开采量吗？ 区域地下水资源评价的核心是确定可持续开采量，可持续开采量评价准则是“补给增量与排泄减量之和”。（陈崇希 1966,1978,1982；Bredhoft 2002）。然而,补给增量与排泄减量是与开采井的层位密切相关，在含水系统的浅部与深部开采地下水的效果（补

43、排的增减量）显然是不同的。因此，仅就“可持续开采量评价”而言，也应按地下水流原本的面目来建模，简化为二维模型是不妥的。 区域水平尺度/垂向尺度影响建模维数? 区域地下水流模型平面尺度大，垂向尺度小，并不能说明流速的垂直分量可以忽略不计。黑河中游地区张液市东25-30km老寺庙农场27#孔(孔深136.25m),水头高出地面7.02m，而同处的潜水位却埋藏深度大约为1上下，两者相差约8m。若采用二维流模型识别，拟合哪个水位水头？如果说承压水头与潜水位相差8m并不普遍,那么两者相差23m，能被忽略吗? 资料少先建二维模型? 关于资料占有程度与建模维数问题，笔者认为：只要大体掌握含水系统的结构、岩性

44、规律(依地貌类型)，补给/排泄位置、方式等，就可建三维模型，我国特有的1/20万的区域水文地质报告为此提供基础，更何况目前绝大部分地区具有更详细的水文地质资料。拟合孔不足是实际情况，这只能说模型识别的精度不高，却与勘查阶段相适应。如果资料较少的勘查初级阶段先建二维模型，那么： 这意味着:将降雨、河渠、灌溉等的水入渗到潜水面的补给（这是地下水的主要补给来源）平均地分配到整个含水系统的厚度上;将潜水面的蒸发量平均地分配到整个含水系统的厚度上;将出露于地表的泉流量平均地分配到整个含水系统的厚度上;将非完整抽水井的流量平均地分配到整个含水系统的厚度上;将非完整河视为完整河，拟合水位用的观测孔，不管它的

45、垂直位置是深部承压水层还是潜水层；我们能接受吗？! 一旦未来地质孔多了，抽水试验多了，拟合用的观测孔多了，认为可以建三维模型时，那么原先二维模型识别的参数则要“推倒重来”。相反，若一开始就建三维模型，当增加几个地质孔时，只要局部修改含水系统的分层标高；当增加几个观测孔时，只要进行模型再识别，改变部分参数，提高精度， 无需“推倒重来”! 如果开始模型识别资料不多时就建三维模型，则可通过所进建模型的灵敏度分析或模型识别过程中的信息来指导未来勘查工作的重点所在，以用较少的勘查量获取更多有效的信息，为下次模型识别提高质量打下基础。上述,就理论上讨论了关于建模维数问题的不同见解。下面以实例说明。例1-1

46、 天桥泉域岩溶地下水流模型 (陈崇希等,1995c)研究区位于晋陕蒙接壤区黄河谷地寒武系中统-奥陶系中统岩溶水天桥泉域，南北长约200km，东西宽约100km，面积约1.4万km2。地势东高、西低；南、北高，中间低。图7-2-1晋陕蒙天桥泉域位置图 (陈崇希等,1995c)区域地质: 研究区处于山西台背斜与鄂尔多斯台向斜的交接处。地层出露较全，属于我国华北典型的地层岩性。地质构造:研究区位于背斜与向斜的交接部位。岩层走向NNE，倾向NWW，产状平缓，黄河以西岩层倾角变陡。本区岩溶含水系统指寒武系中统张夏组至奥陶系中统上马家沟组的碳酸盐岩含水岩系。以奥陶系中统峰峰组岩层（白云质泥灰岩或泥质白云岩

47、）及石炭系（页岩、泥岩及砂岩为主的含煤地层）为其隔水顶板。以寒武系中统徐庄组（泥质灰岩、砂岩及页岩） 及寒武系下统石英砂岩为其隔水底板。在研究区的东北部，含水岩系裸露；东南部，含水岩系常有第四系或与第三系覆盖；中部黄河东岸多有石炭系覆盖；黄河西岸岩溶含水岩层倾角逐渐变大，超过45，埋深迅速变大，上部多为石炭系及二迭系覆盖（图7-2-）。图7-2-晋陕蒙天桥泉域地质略图(陈崇希等,1995c)图7-2-晋陕蒙天桥泉域地质详图(陈崇希等,1995c)黄河于天桥背斜轴部（在山西铁匠铺北数公里处）流过，河床切入奥陶中统上马家沟地层，地下水集中由此“天窗”排入黄河，形成著名的“天桥泉” *（严格地讲并非

48、泉，而是地下水较集中地从黄河河床底面排泄入河），流量达12.46m3/s。边界条件:岩溶含水系统的北界为含水层尖灭或与花岗岩接触；南界（除有1km开口外）及东南界为与变质岩接触；东北及东界为地下分水岭；上述北、东、南处理为隔水边界或零流量边界。西部寒武系中统-奥陶系中统碳酸盐岩层已深埋，岩溶不会发育，取一个足够大的埋深（800m）近似作为隔水的西界。上世纪90年代前后几年内，已做过几次数值模型研究,都是平面二维流模型，且对泉流量采取人工给定*（当时对泉的模拟尚缺乏经验）。 上世纪90年代初，笔者承担国家“八五”科技攻关项目的专题，重新建立天桥泉域的地下水流模型。建模最为关键的是：本区寒武系中统

49、张夏组奥淘系中统马家沟组岩溶水能否用平面二维来刻画？图7-2-3a由占孔水位所做的等水位线,黄河东西侧地下水都流向天桥泉。通过天桥泉东西向剖面（图7-2-a）所示，“泉域西界以埋深800m为隔水边界”。既然是隔水边界，它便不能作为流线的起点(平面二维流将其作为水流来源的边界)，此处无源头，补给源在天桥泉排泄区的东侧，西侧只能起到岩溶水“过路”的作用。因此，流线应当如图7-2-所示：在排泄区东侧近黄河地段，上部地下水宏观上大体水平地直接流向黄河，在黄河河床底附近转为垂直向上流。而在远离排泄区东侧的补给区，如地下水分水岭附近，地下水以垂直向下流为主，而后沿隔水底板流动，越过黄河河床底流向排泄区的西

50、侧，受西侧隔水边界阻挡，再沿隔水顶板向东流排泄入黄河底，显然存在“绕流”现象笔者称之为“绕流模型”。这种流动形式是不可能由平面二维流模型来刻画的。（陈崇希等,1995c) 图7-2-3 天桥泉域地下水绕流剖面示意图(陈崇希等,1995c) 另外，平面二维流模型必须满足裘布依假定，然而水源地水文地质勘查表明，当占孔揭穿顶板峰峰组之后，地下水水头一般高出地面，且随占孔的加深，水头不均匀地增高。终孔时，水头一般高出地面1623m。这一现象从另一侧面证明该处地下水的垂直流速不能忽略不计，裘布依假定不能成立。因此本区地下水属三维流系统，这是水文地质概念模型的根本所在，否则，将导致模拟失真。”(陈崇希等,

51、1995c) 上面分析天桥泉域地下水流存在“绕流”现象,论证了建立三维流模型的必要性。 该区曾做过极珍贵的大型抽水试验(最大流量达3.09m3/s, 抽水试验持续时间46d),但可惜的是由于含水层埋深大，观测孔虽有24个，但主要沿黄河岸边浅埋部位展布，使得做初始等水头线提供初始条件遇到困难。这个问题不解决，数值模型不可能建立。我们根据地下水动力学中Theis公式讨论之6(陈崇希,1983,59页；或陈崇希等,1999,82页；2011，115页)将以水头H为因变量的数学模型改为以降深s为因变量的数学模型,从而将原来的初始条 改为 ,即全区抽水试验初始时刻的降深处处为零而得以解决!当时所具有的资

52、料，能够建立三维流模型吗?基岩含水系统建模的另一个基本问题是各向异性问题, 本区的放水试验,其等降深线具椭圆形,证明了这一特征。由于基岩含水系统各向异性的主方向受构造控制,一般不可能刚好为正东、正北和正上(倾角90)方向，因此水流的控制方程不能写为：而应写为：7-2-17-2-2(7-2-1) (7-2-2) 一些文献忽视了这一点。例如：某数值模拟论文述：“泉域区内寒武-奥陶系碳酸岩含水岩组中广泛分布的裂隙主要有两组,倾向分别为318和312”吴吉春等,2001,然而却将控制方程写成(7-2-1)式的平面二维形式。笔者并不认为该文作者理论上不理解这一点,而是疏忽。但此疏忽却造成不良影响。 天桥

53、泉域放水试验形成的椭圆形等降深线,长轴为NNE向约70, 根据椭圆长、短半轴之比,依地下水动力学理论(陈崇希,1983，(6-2-32)式；陈崇希等,2011,(7-2-32)式), 各向异性渗透系数之比约为32。这一点为模型识别求取水文地质参数提供重要信息。最终,我们建立的天桥泉域张夏组马家沟组岩溶地下水流的模型,是各向异性、三维“绕流”不稳定流模型。后话 今天返回去看“天桥模拟报告”,似乎存在几点不足:应设法建立预测天桥黄河正常水位时的泉流量(上 述12.46m3/s是黄河断流时的流量)；预测两水源地运行条件下天桥泉流量的衰减动态。例1-2 柳林泉域岩溶水流的模型山西省柳林岩溶水泉域处于吕

54、梁复背斜-鄂尔多斯台向斜的交接处，为一单斜构造。地层出露较全，属于我国华北典型的地层岩性，寒武系中统-奥陶系中统碳酸盐岩层构成岩溶含水系统(下页示意图)。东侧为潜水，西侧为承压水。柳林泉出露在奥陶系与石炭系的接触带地段，这是典型的“单斜含水层-泉流系统”。要正确建立地下水流基本模型，必需对地下水流（流网）的基本特征作正确的分析。“单斜含水层-泉流系统（图7-2-a）是自然界最常见的地下水系统之一。我国南方小型构造盆地（如湖南省的斗笠山盆地，恩口盆地，桥头和盆地等）和北方大型盆地（如鄂尔多斯盆地等）周边的泉（群），大多可概化为单斜含水层-泉流系统。遗憾的是，有的模型设计者采用平面二维流模型刻画这

55、类含水系统”(陈崇希等,2002a)。图7-2-单斜含水层-泉流系统(陈崇希等,2002a)1. 潜水面;2. 隔水边界;3. 岩溶不发育地下水滞留区;4. 流线;5.泉;6. 观测孔;7. 降雨补给单斜含水层-泉流系统分承压与无压两区。垂向上的上边界：承压区的顶面隔水边界可依地层岩性很容易确定，无压区以有降雨入渗补给的潜水面、泉(或加非完整河)为其顶面边界，其底面隔水边界除根据地层岩性外，还要依其岩溶发育随深度的衰减情况分析确定。平面上的边界随具体条件而定。上述的是含水层系统框架，它的划定一般并不困难。“重要的问题是，如何分析这种条件下地下水流动的基本特征，或者说，怎么勾画流线”。一些地下水

56、流模拟者不会绘制复杂条件下的流线，导致基本模型类型设计的错误不是个别例子。(陈崇希等,2002a) “20世纪60年代初，笔者带学生在湖南省斗笠山、恩口实习，提出应画成如图7-2-b所示的流线。如此勾画的流线反映了单斜含水层-泉流系统（以前曾称为“自流斜地含水系统”）的基本特征。另外，在泉口附近打孔，随着孔深的加大，水头上升，这一现象对此也是一个佐证。可见，单斜含水层-泉流系统的水文地质概念模型是三维流。”(陈崇希等,2002a)图7-2-单斜含水层-泉流系统(陈崇希等,2002a)1. 潜水面;2. 隔水边界;3. 岩溶不发育地下水滞留区;4. 流线;5.泉;6. 观测孔;7. 降雨补给基于

57、上述水文地质条件的流线（流网）分析，若将其概化为水平二维流模型，就意味着：地下水向“泉”的流动概化为向“完整井”的流动；然而所模拟的柳林岩溶水泉域的“寒武-奥陶系含水岩组”(吴吉春等,2001)厚度依地层表估计超过600m,即意味着模型中的抽水井都超过600 m。读者能接受吗? 概化为水平二维流模型你能接受吗?经常可以从文献中看到：拟合较好，说明所建模型合理、可靠。要明确，拟合求参一般是多解的。水文地质模拟者应当下大力气去合理概化水文地质模型,并采用正确的仿真技术。只有在此基础上,才有真正意义上的“拟合”。 (陈崇希,2003a)例1-3 含混合抽水井孔的三维流问题，因故降为二维流模型。(1)

58、认为研究区存在大量的混合井孔“形成人工天窗，使层与层之间互相沟通，故所测水头已不再是某层的水头，而是潜水层和承压水层的综合反映，概化为平面二维流混合层”模型。(桑学锋等，2007）这种见解是比较普遍的。如此分析是否合理？(2)某专著（2007，259页）认为：在三维流水文地质条件下,如果“水头是在装有长过滤器的井中测量的，但三维模型要求的是点上的水头值，所以这种水头值只适用于二维区模型。” （薛禹群,2007，259页）读者如何认识? 如何评价?这是认识到建模区存在混合观测孔，因不会模拟混合观测孔，无奈而建为二维流模型。但此条件下建为二维流模型就正确了吗？难道存在混合观测孔的三维流，流场会转变

59、成二维流？!(3)有的虽然也建为三维流模型（如某海水入侵模型），但论文、专著中不提该区存在混合井孔，却悄悄地加以处理。这就不属于技术问题了，而是学术诚信问题！例1-4 郑州黄河滩地三个含水层混合抽水试验确定水文地质参数这是一个笔者未参与设计和现场试验工作，只是后来见到混合抽水试验求取水文地质参数的技术报告因存在混合抽水而采用平面二维流模型，我们用以重新建模混合抽水试验三维流模型确定分层水文地质参数(陈崇希等,1985b,1998b)*。试验区位于郑州市北郊，北依黄河，南连黄泛平原，限于黄河与大堤之间的滩地（图7-2-4），面积约30km2。混合抽水试验场涉及7080m以浅的三个孔隙含水层及其间

60、的两个弱透水层。黄河水是本区地下水最主要的稳定的补给来源，此段黄河属于“悬河”（图7-2-5），以“注水式”补给地下水（见4.10），而近水平地向南流出本区。抽水井组北距黄河1300m，如图7-2-4所示。3个贯穿、含水层的混合抽水井位于边长15m的等边三角形顶点。三个抽水井中间设一个3层混合观测孔（中心观测孔），其它观测孔分布在外围。共有18个观测孔，其中5个为三层混合观测孔，其他13个为、含水层混合孔。三个混合抽水井的流量达 1. 38万m3/d，抽水历时7d。图7-2-4试验场位置图图7-2-5 地下水补给和排泄示意图依以上述条件建立水文地质概念模型为：非完整河旁（黄河仅切至含水层），三