水文地质边界条件处理-陈崇希
e4常规观测孔水位的形成机理及确定方法
一般观测孔实质上是一部分滤管进水 （ 抽水） 而另 , !（ "，’，&） H’ * （ , ） #$ ) %$ %$ I 一部分滤管出水 （ 注水） ， 即观测孔并不只简单地反 式中， #$、 %$是观测孔滤管顶点、 底点的标高* 本文称 映含水层地下பைடு நூலகம்的水头， 而是兼有抽水与注水作用
上它只有一个点与含水层相连通， 反映该点的水头 值； 另一类是常规观测孔， 它有一定的长度与含水层 相连通* 通常所谓观测孔是指后者* 美国著名学者 、 前 苏 联 著 名 学 者 !"#$%$& 和 ’$&()(* 在对承压含水层非完整井流的研究中， 都认为观测 ./01234 孔中的水头降深 （ ! "， #$， %$， &） 反映该孔滤水管中各 点降深 !（ ’，&） 的平均值， 即 I "， （ ! "，#J，%J，&） (
#$ 水流机理与模拟方法
如图 # 所示， 在地下水三维流场中， 穿入一根不 抽水的滤管— — —观测孔， 原来流场中滤管位置处的 水头值不相等 （ 如图 # 所示下部水头值高于上部水 头值） ， 滤管穿入后， 由于滤管的水流阻力远远小于 原来的孔隙介质， 因此滤管中的原始水头差导致比 原来孔隙介质条件下大得多的垂向流速和井筒垂向 流量( 如此， 依水流连续性原理， 井管下部必须从含 水层中进水— — —抽水； 井管上部必须向含水层中出 水— — —注水( 滤管中水流流速 （ 流量） 的增量 （ 与设 置观测孔之前的孔隙介质比较） 和滤管与孔隙介质 间水量的交换 （ 抽水、 注水） 必导致观测孔及其周围 含水介质水头的再分布( 总的趋势是滤管内的水头 差将会缩小， 但水头值不会变成相等， 否则， 滤管内 万方数据 水流将会停滞， 违背了部分进水部分出水的水流连
第4章-2边界条件源汇等特殊问题的模拟与处理资料
rw rwe A
（4-7-4）
若将此有效井径 rw 直接替代上节井水位的校正方程 中的井径,即可将两种校正合二为一。
我们注意到，附加水位降深 sA 与抽水井流量Q成正 比，而“§4.6 抽水井水位的校正方法”中的由于径向流 的校正值也是与抽水井流量Q成正比，为此，实际应用中 可将两者合二为一,用改变抽水井所在格/结点的含水层渗 透系数K,来校正该点的水头值——抽水井水位值。PGMS 软件采用参数CKK乘以渗透系数K来实现此目的。当CKK ＜0时，抽水井水位向下校正；当CKK＞0时，抽水井水 位向上校正。
“井周扰动效应”是指由于钻进方法、 抽水井的结构或成井处理等等的不同，使 抽水井附近含水层的渗透性发生变化，而 对抽水井中水位降深产生的影响(图4-7-1)。
s=sN+ΔsA 在油藏工程中称为“井孔的表皮效应”。
图 4-7-1 井周扰动现象及其漏斗
对此，可直接利用文献(陈崇希,1983§４.７） 中已建立的有关方程。即
0.75
Q 2 T
ln
RA rw
0.75
Qt RA2 e
（4-8-2）
因此，点井中的水位降深sw与面井降s之差Δs为
s
sw
s
Q 2 T
ln
RA rw
0.75
Q 2 T
ln
0.472
RA rw
（4-8-3）
式中的“圆形含水层的半径”可按下式折算
RA
FA
（4-8-4）
对于非达西泉流的模拟, 只要将上式中的渗透系数改为等效渗透系 数和即可实现。我们曾模拟了理想模型的非达西流的泉流量动态(见图417-6)，以及广西环江县北山岩溶管道-孔隙三重介质泉流量动态(图4-92)。喀斯特岩溶裂隙泉流量模拟之所以困难,是由于泉附近局部空隙结构、 尺度变化对流量会起决定性作用，而这些多方面复杂因素在勘查中难以 获得信息，使得泉流量的拟合成为多解。正因为如此,我们更加强调:模 拟者特别需要从水文地质条件方面分析入手，把握泉水的来源及控制因 素。
关于稳定井流模型和Dupuit公_省略__对陈崇希教授_商榷_一文的答复_薛禹群
水利学报SHUILI XUEBAO 2011年10月第42卷第10期文章编号：0559-9350（2011）10-1252-05收稿日期：2011-07-17作者简介：薛禹群（1931-），男，江苏无锡人，教授，中科院院士，主要从事地下水动力学和地下水数值模拟研究。
E-mail ：yuqunx@关于稳定井流模型和Dupuit 公式的讨论——对陈崇希教授“商榷”一文的答复薛禹群（南京大学地球科学与工程学院水科学系，江苏南京210093）由于忙于其他方面的工作，没能定期阅读《水利学报》，因而不清楚《水利学报》开展了有关稳定井流模型和Dupuit 公式的讨论，陈崇希教授2010年8月发表于《水利学报》的文章（文献［1］）也迟至近日才看到，因而迟复了。
甚歉！兹就有关问题谈谈个人的看法。
1地下水的流动状态关于地下水的流动状态我们一贯认为“地下水的状态总是在不断发展、变化着。
所谓稳定只是有限时段内的一种暂时平衡现象。
当水位变化很小，即源汇项w ()x ，t →-w ()x ，渗透系数K ij →-K （式中-w ，-K 分别为w ，K 的极限值）、边界条件也不随时间变化时，便有∂H ∂t →0，可视作稳定流问题来研究”（式中H 为地下水水头，t 为时间）［2］。
在地下水动力学［3］中也贯彻这样的观点，所以在其第三章讲述稳定运动的一开头，就说明“随着抽水时间的延续，降深不断增大，漏斗不断扩展。
若没有其它补给源时，地下水向井的运动始终处于非稳定状态”（文献［3］，62页）。
但教材不同于专著，专著可以阐述作家一家之言，文献［2］就是这样做的，除了前述简单的说明外，就没有再提及稳定流了。
教材则不行，除了要考虑该学科的研究进展外，还必须考虑教学大纲的要求和广大产业部门的需要与看法。
产业部门目前还大量使用着稳定流公式，如坝下渗流、绕坝渗流、水库渗漏、干扰井群、边界附近井流、各种不完整井流、承压－潜水井流的公式等等，并不是他们不知道稳定流公式有这样那样的缺陷和假设条件，而是因为目前利用非稳定流理论来推导这类解析解公式有的还有困难，有的即使推出，公式往往过于复杂，不便应用。
地下水运动数值模拟过程中边界条件问题探讨_卢文喜
2003年3月水 利 学 报SHUILI XUE BAO 第3期收稿日期:2001-11-14作者简介:卢文喜(1956-),男,吉林德惠人,教授,博士生导师,主要从事生态水文和地下水系统数值模拟和优化管理方面研究。
文章编号:0559-9350(2003)03-0033-04地下水运动数值模拟过程中边界条件问题探讨卢文喜1(1.吉林大学环境与资源学院,吉林长春 130026)摘要:本文对地下水运动数值模拟过程中边界条件的涵义和处理方法进行了分析和讨论。
阐述了边界条件所包含的双重意义。
指出随着人类活动影响强度的日益增大,边界条件的处理要面临一些新的更为复杂的问题。
在模型预报之前必须首先对边界条件做出预报。
边界条件的预报既要考虑自然因素的作用,同时也要考虑人类活动(人工开采和人工补给)的影响及由于邻区水流条件变化而产生的耦合效应。
之后,给出了两个应用实例。
关键词:地下水;数值模拟;边界条件中图分类号:P641.2文献标识码:A在地下水运动数值模拟的过程中,模拟预报结果的正确与否与边界条件处理得是否恰当密切相关[1,2]。
尤其是在人类活动影响强度日益增大的今天,在处理边界条件时,常常会面临一些新的更为复杂的问题。
原因在于边界处的水流状况往往不仅受到自然因素的控制,而且还深受人类活动(如人工开采和人工补给)的影响[3,4],同时还可能受到邻区水流条件变化的扰动,而对于人为边界更是如此[5]。
所以必须对边界条件给予应有的重视,深入探讨其多重的内涵并研究出切实可行的处理方法。
1 边界条件涵义探讨在地下水运动数值模拟的过程中,一般都是在概念模型的基础上,建立描述地下水流的数学模型,然后再采用某种数值方法,对模型离散并求解。
对于分布参数的地下水流数学模型而言,模型主要由两部分内容组成:①描述地下水运动规律的偏微分方程;②反映地下水模拟区域具体特征的边界条件和初始条件(若为稳定运动则没有初始条件)[6,7]。
这里的边界条件具有两重意义:一是它与初始条件一起构成地下水流数学模型的定解条件,用来说明具体目标系统的边界所具有的特定状态,从而使模型的求解能够得到切合实际状况的特解。
已阅 地下水流模拟系统PGMS_1_0版_简介
2 PGMS 软件结构
110 版 PGMS 软件是运行在 Microsoft Windows 平台 上的可视化地下水三维流模拟系统 ,主要由 4 个可执 行程序组成 :平面三角形网格剖分程序 (主要进行平面 三角网络剖分) 、PGMS 前后处理过程的主控程序 (处 理三维可视化模型与核心计算程序之间的数据交换) 、 基于多边形网格的三维有限差分模型的数值计算程序 (三维饱和水流计算程序) 和含水层三维结构动态演示 程序 (显示三维实体模型) 。
(1) 无 需 人 工 分 配 混 合 抽 水 井 各 层 的 流 量 : MODFLOW[1] 建议 “: 多层井的流量必须以某种方式人 为地分配给每一单层 , ……,把井流量按每层的导水系 数大小分配”。多年的研究表明此建议不妥 。PGMS 符合机理地模拟混合抽水井 ,各层的流量是模拟的结 果 ,不可人为分配[2～3] 。
参考文献 :
[ 1 ] Mcdonald M G, Harbaugh A W. A modular Three dimensional finite2difference groundwater flow model [ R] . Techniques of Water Resources Investigations Book 6. Chapter A1. Washington D C : U S Geological Survey United States Government Printing Office , 1988.
2007 年第 6 期
水文地质工程地质
Ⅰ
地下水流模拟系统 PGMS( 110 版) 简介
-四小时学会Modflow
该格元所控制的“排水沟”， 平均排水高程Di,j,k 。 所谓水头差，是指格元地下水位 与该排水高程之差。
用折线输入法输入排水沟参数：
QDi, j,k K L (H i, j,k i, j,k i, j,k Di, j,k ) CDi, j,k (Hi, j,k Di, j,k ) 当Hi, j,k Di, j,k
QRivi, j,k CRivi, j,k (HRivi, j,k Hi, j,k ) 当Hi, j,k HRivi, j,k
当河流向地下水渗漏时
当Hi, j,k Brivi, j,k QRivi, j,k CRivi, j,k (HRivi, j,k Hi, j,k ) 当Hi, j,k Brivi, j,k QRivi, j,k CRivi, j,k (HRivi, j,k Brivi, j,k )
前处理软件－协助自动或半自动制作数据文件集合，
具有校对、图形显示、人机交互特征 后处理软件－具有自动整理结果数据文件为简洁易懂的图表
等功能
数值模拟：常用商业软件
地下水模拟程序 美国地质调查局(UGS): MODFLOW 地下水流模拟程序 MT3D、MT3DMS…溶质运移模拟程序
集成前后处理程序 Visual ModFlow Processing ModFlow(PM) Groundwater Modeling System(GMS) Visual Modflow ……
自动化程度高！
在某地下水模型中，选用了“排水沟”模块包 问题：模拟程序运行后，发现排水沟出水太多 了，怎么办？
解决方法 1 沟排水太通畅了，把系数Cd或K减小一点，把地
下水憋一下。或者 2 沟的排水高程 D 太低了，高点试试 3 两种方法一起上，再试试！
潜水完整干扰井非稳定流抽水试验计算水文地质参数
潜水完整干扰井非稳定流抽水试验直线图解法（水文地质参数计算）（发表于《内蒙古水利》（2011年第1期），荣获“2010—2011年度科学发展与构建和谐社会理论实践成果”一等奖，颁奖单位：四川西部文献编译研究中心、《环球人文地理》（理论版）编辑部、中国西部科技杂志社）内容摘要：对于潜水完整干扰井非稳定流抽水试验计算水文地质参数，可应用类似的承压--无压水定解问题的太斯解析解公式，依据“势能叠加原理”求出干扰井解析式，再进行承压水降深S 和潜水含水层天然厚度H与抽水计算时刻观测井中水柱高h的平方差（H2-h2）近似代换之后，得到相应的潜水干扰井太斯解析解公式。
然后，再进行雅可布式简化，得到便于图解计算的潜水干扰井雅可布公式，类似于承压水雅可布图解法，求得水文地质参数K的较精确值和μ的可参考值。
引言：在水文地质勘查实践中，为更清楚地了解含水层水动力学特性，能够掌握较切合实际的水文地质参数，为开采井布局和设计提供较为可靠的出水量及其与水位变化关系的依据性资料，常布置实施多主孔完整干扰井非稳定流抽水试验。
对于潜水含水层，至近年来的规范规程和教科书，都未直接给出完整干扰井非稳定流抽水试验计算水文地质参数的具体公式，在生产实践中，常需水文地质人员自行推演计算公式，缺乏表达的统一性和应用的普遍性。
本文根据地下水动力学教程提示，对潜水完整干扰井非稳定流抽水试验计算水文地质参数方法和公式进行论述。
一、水文地质模型1、含水层在平原区或其它潜水含水层分布较广、连续发育的地区，含水层厚度分布较稳定，岩性较均匀，底板倾斜度较小，适用于地下水动力学方法计算抽水试验水文地质参数。
2、边界条件边界距抽水井距离大于500m，可按无限边界处理。
图1 潜水含水层完整干扰井非稳定流抽水水文地质模型示意图二、单井数学模型及其解析解对于分布连续广泛，厚度相差不大，岩性较均匀的承压—无压含水层中完整井井流，常用如下数学模型描述：其中：t—抽水时间（T）；r—计算水位点与抽水井距离（L）；Q—井孔抽水量（L3T-1）。
水文地质边界条件处理-陈崇希
（1）含水层导水系数对井孔流量的影响， 不会如此简单。例如,混合井附近岩性(渗透系 数)发生变化，甚至混合井打在岩性透镜体上， 怎样影响井管的流量分配？含水层厚度发生变 化又如何改变流量的分配？
K1
K3
K2
混合井管贯穿岩性透镜体
图4-1-3 混合井管贯穿岩性透镜体上
（2）含水层的参数影响混合井流量的分配，导水系数只 是其中一个因素，含水层的弹性给水度（储水系数）就不起作
第5章 地下水流模型的拓展
§5.1 地下水饱和-非饱和流耦合模型 §5.2 地下水-地表水流耦合模型 §5.3 分布式水文模型简介
第6章 数值模型设计 应用及水文地质勘查
§6.1 数值模拟设计应用主要步骤 §6.2 水文地质概念模型的设计
§6.3 数值模型的建立
§6.4 预测模型若干问题 §6.5 关于基岩含水系统建模的特殊问题 §6.6 数值模型水文地质勘查等若干问题
这种井群干扰下，原混合井的流量比还会保持常量吗？
（5）混合观测孔是混合抽水井的特殊情况（Qw =0）(图4-14) ，对于两层混合的观测孔，其孔中水位(混合水位)必界于两
含水层水位之间，即混合观测孔对于其中一含水层(例如2含水 层)起抽水作用（Q2>0）, 对于另一含水层（1含水层）起注水 作用（Q1<0）。如此，Q1/Q2<0。而两含水层的导水系数的比 值肯定是正值，即T1/T2>0。如此，两个比值怎能相等？！
若记Qwi为其源汇项，包括开采井、泉流量、河流的补给、 排泄、降雨入渗等，则方程可以写为，
hn,m1 hin, m1 hkn,m1 hin, m1 hin, m11 hin, m1 Tij , m j pb Tik , m bq Ki , m 1 / 2 Ai zi , m 1 zi , m ij ik e hin, m1 hin, m11 hin, m1 hin, m Ki , m 1 / 2 Ai Qwi , m i , m Ai zi , m zi , m 1 / 2 tn 1 tn
不同地质边界条件岩溶隧道涌水量预测及展望
该方法为隧道影响范围 (或集水范围 ) 内泉水流 量的总和 , Q = ∑Q泉 。即隧道施工影响范围所有泉
i=1 n
水 ,成为隧道涌水量的补给源 。
1. 6 地下动力学法
该方法适用于均质地层 ,即在同一水文地质单元 , 利用抽水试验取得渗透系数 ( K ) , 再选择相似边界条 件的公式计算全段隧道涌水量 。 岩溶地区为非均质含水层 , 抽水试验求渗透系数 的前提不对 , 将岩溶作为均质或似均质地层对待 , 求得 隧道涌水量仅具参考价值 。
94
我国可溶岩地区占全国总面积的 1 /8, 分布十分 广泛 ,几乎遍及全国 。由于南方地区气候湿润 ,雨量充 沛 ,岩溶发育普遍 ,其中西南地区以典型裸露为特征 , 而其余地区以覆盖充填型较为普遍 。岩溶发育呈现为 非均质地层 ,岩溶地下水活动多处于非均质 、 非稳定的 紊流状态 ,边界条件极为复杂 。虽然多年研究取得很 多成果 ,但由于排水水文地质与供水观念上相差甚远 , 各种水文地质参数的来源 、 取值 、 运用及涌水量的计算 有较大的出入 ,因此预测难免出现较大的误差 。 南岭隧道施工对整个隧道范围的主要水文地质点 进行了监测 ,并对隧道涌水量进行了测定 ,收集了较为 充足的数据 ,取得了丰富的第一手资料 。通过对资料 的分析整理 ,不仅对隧道施工岩溶水文地质条件变化 及涌水曲线特征有了较为深刻的认识 , 而且对岩溶水 文地质参数及隧道涌水量预测提出新的概念及方法 。
Pred iction of Wa ter - bursti n g Volume in Karst Tunnel under D ifferen t Geolog ica l Boundary Cond ition s
地下水开采一地面沉降模型研究
水文地质工程地质
地下水开采 一地面沉降模型研究
陈崇希’ ， 装顺平， （ １ ． 中国 地质大学・ 武汉 环境地质研究所， 湖北武汉 ４ ３ ０ （ Ｙ ７ ４ ； ２ ， 中国地震局地球物理研究所， 北京 １ ０ ０ ０ ８ １ ）
误差还是小于５ ％吗？本研究带着此问题做了数值模 拟研究， 发现 ：
｝ １ Ｊ ａ ｖ ａ ｎ ｄ ｅ ｌ 和Ｗ ｉ ｔ ｈ ｅ ｓ ｐ ｏ ｏ ｎ 的“ 误差小于５ ％ ” 的结
论仅适用于弱透水层为各向同性介质； ②对于弱透水层的各向异性比为 １ ０时， 误差已超
口护
过２ ７ ． ７ ％， 而且随着弱透水层单位储水系数ｆ ， ． 的减小
和模 拟 时 间 的 延 长， 误 差 还 要 增 大， 己存 在 超 过 ３ ０ ． ５ ％的情况。 基于此研究结果， 本专题对苏州市采用三维不稳 定流模型。
３ ． ２ 土层固结与地下水流如何藕合？— 根据单位 储水系数的物理意义将两者藕合起来 １ ９ ８ ９ 年比利时地调所做的《 上海地面沉降模型） ， 其沉降模型是“ 地下水流动模型加上一维垂直流动固 结模型” 。即“ 模型分成两步” ： 先由水流模型给出每一
１ ９ ９ ５ 年分别对上海市和上海市浦东区及天津市做了
三个地下水开采一 地面沉降模型， 都是准三维流模型。 国外 ９ ０年代做的几个主要实例模型 ． 例如 Ｒ ． Ｂ ｒ ａ －
ｖ ｏ （ 美国） 等做的美国休斯敦模型， Ａ ． Ｒ ｉ ｖ ｅ ｒ ａ （ 法国） 等做 的墨西哥城模型， Ｇ ． Ｃ ａ ｍ ｂ ｏ ｌ ａ ｔ ｉ （ 意大利） 做的意大利拉 沮纳区域地下水流模型， Ｋ ． Ｄ ａ ｉ ｔ ｏ （ 日本） 做的 日本大绍 平原的模型等等， 都属于准三维流模型。 我们注意到， 长江三角洲的沉积物， 弱透水层多呈 “ 千层饼状” 的岩性， 它们在模型的“ 层” 中表现出各向
地下水动力学,稳定井流与非稳定井流
5结论
理论推导与实际算例表明文中提出的非稳定井流试验数据分析方法是可行的。与现有的方法相比较，具有：（1）不论稳定与非稳定井流的情况下均可以应用；(2）仅要求计算井损参数的情况下，可以不考虑非稳定流井函数的具体形式；（3）非线性指数不论是否已知的情况下均可以应用；（4）全部数据分析过程可以程序化，利用计算机完成全部计算过程等优点。
关键词：稳定井流；非稳定井流；泰斯公式；裘布依公式；井函数
1地下水动力学发展史
在18世纪中期开始，一些法国工程师和科学家的杰出成就，奠定了地下水渗流力学 的理论基础。这当中主要包括达西定律 、裘布依假设 、布西涅克斯潜水运动方程 。达西定律是法国工程师达西在解决城市给水问题时，根据均质砂中垂直水流实验结果，在1856年总结出线形渗流方程 ，即地下水的渗流速度v与水力梯降J成正比的线形关系；该线形渗流方程也就是著名的达西定律，它的建立是渗流力学诞生的标志。裘布依假设是法国工程师、水力学家裘布依针对缓变流动的潜水,于1863年提出用潜水位h代替侧压水头,这种处理方法使得同一剖面各点的渗透速度相等。得益于裘布依假设，达西定律在实际中被迅速推广，这也使得渗流力学得以迅速发展。布西涅克斯方程是在1904年法国数学力学家布西涅克斯在认为水是不可压缩的条件下，利用裘布依假设，给出潜水渗流运动的微分方程,这为非稳定流理论 的发展奠定了基础。他创造性地将坐标原点取在含水层底板上(以下坐标都是这种设置方法) ,这使得方程中的水位变量与潜水流厚度相等,这极大地方便了方程在实际中的应用。
4实用分析
在大部分地质报告，无论是抽水试验还是矿井涌水量预测，特别是通过抽水试验来求含水层参数，多数是采用稳定井流解析法。但是一些地质报告在采用这种方法是时却忽视对水问地质条件的分析，忽视当地是否存在稳定井流的可能性，似乎以为“稳定井流”是一种不受条件限制的、可以任意选用的计算方法，但是这是一种错误的观点。
数值模拟在地下水环境影响评价中的应用
数值模拟在地下水环境影响评价中的应用作者：王志刚等来源：《计算机光盘软件与应用》2013年第24期摘要：本文利用地下水数值模拟软件Visual Modflow，在详细研究某电厂周边水文地质条件的基础上，建立了研究区的地下水流动模型和地下水溶质运移模型，对电厂事故工况下泄露到地下水中轻质柴油的运移规律进行了模拟研究，就其对地下水的影响进行了分析评价。
结果表明：模拟期内柴油对含水层造成了污染，并且存在超标现象，但没有对保护目标造成明显影响。
结合数值模拟过程，提出了数值模拟在地下水评价应用中的一些注意事项，可为今后的地下水评价提供参考。
关键词：数值模拟；地下水；环境影响评价中图分类号：X820.3随着我国经济社会的发展，人们对环境质量的要求越来越高[1]，而同时，我国的环境污染问题，尤其是地下水环境污染的问题越来越严重。
作为重要的环境因子，地下水已经成为城市和工农业用水的主要水源[2]。
在经济社会可持续发展的前提下，防治地下水环境污染已引起人们的重视，全面、准确地评价建设项目对地下水环境的影响对于保护地下水资源意义重大[3]。
环保部2011年6月发布了《环境影响评价技术导则—地下水环境》（HJ610-2011）[4]（以下简称《导则》）正式作为我国建设项目地下水环境影响评价工作的规范和指导，这就对地下水环境影响评价提出了更高的要求。
随着科学技术的发展，数值法越来越成为解决地下水问题的重要方法，数值模拟技术以其方便、灵活和高效的特点在水文地质领域得到广泛应用[5-7]。
而《导则》要求，对于地下水环境影响一级评价，必须采用数值法进行预测和评价，对于地下水环境影响二级评价，可选择采用数值法或者解析法进行预测和评价。
可见，数值法在地下水环境影响评价中占有重要地位。
然而，由于自《导则》开始实施以来的时间较短，数值模拟在地下水环境影响评价中的应用还存在诸多方面的问题。
本文拟通过建立某电厂厂址区域的地下水数值模型，利用Visual Modflow软件模拟事故工况下污染物在地下水中的迁移规律，进而就本工程对该区域地下水环境的影响进行预测和评价。
地下水溶质运移理论及模型读书报告
式中，W为球面积；n为有效孔隙率；JD为弥散通量，且 ,Vv为均衡段空隙体积。
忽略高阶微量，化简后得：
于是该点源的定解问题可以写成：
(R≧0,t>0)
(R>0)
(t>0)
(t>0)
(t>0)（该式将点源处浓度限制在有限区域）
通过Boltzmann变换，将原来的偏微分方程定解问题转变为常微分方程定解问题，可求得空间瞬时点源的解为：
III、一、二维水动力弥散的差分解法的比较
相同点：都采用类似的差分原理进行差分，得到的差分格式的基本类型一致。
区别：一维条件下3种格式采用“追赶法”求解，的三对角线方程组。而二维条件下所给出的三种格式组成的方程组是五对角线方程组。为了避免解五对角线方程组的困难，特提出交替方向隐式法，简称ADI方法。它的优点是：不是一次对整体矩阵求逆，而是分两次对三对角线矩阵求逆，这样就把二维问题简化为多次解一维问题。
注入空间连续点源时，假定注入的是理想示踪剂。将连续点源视为无数的瞬时点源之和，直接利用空间瞬时点源的解，利用积分得出解。当时间足够长时，该问题的解为：
2、
数值解法可以应用于复杂的情况，在实际应用中起着很好的效果。现就书上出现的数值解法做一些简单介绍。
（1）有限差分法
有限差分法的基本思想是：将研究空间划分成许多小的网格，把时间分成许多小段 ，每个网格中心点处的未知变量视为该网格上的平均值，然后利用差商近似代替微商，形成研究区域上离散分布的有限个代数方程，求解方程组便可得该 时刻上各格点上的取值。然后按照一个个的 逐个往前求解。
迦辽金法属于加权剩余法，且由于其他加权剩余法，应用更普遍。
煤矿水文地质勘测中存在的问题及完善对策研究
煤矿水文地质勘测中存在的问题及完善对策研究摘要：煤矿水文地质勘测工作对于矿区安全生产发挥着重要作用。
本文在介绍煤矿水文地质勘测概念和工作内容基础上，分析了煤矿水文地质勘测过程普遍面临的问题，针对上述问题，从加强多个角度探讨了煤矿水文地质勘测工作的有效措施，期望对今后煤矿水文地质勘测实践带来有益的思考。
关键词：煤矿地质；防治水；问题；有效措施1.煤矿水文地质勘测中存在的问题1.1水文地质资料不完整在开展煤矿勘测工作之前，需要了解矿井地质条件及周边的环境，包括矿井充水源、水害分布及积水量等内容，但是一些煤矿企业没有认真研究矿井充水条件，特别是充水性图。
作为一种对井下水文地质资料实际测量并综合记录的图纸，矿井充水性能图详细记载矿区充水规律，同时也是矿井勘测及制定防水方案的重要依据。
没有经过水文地质调研，所获取的勘测信息，缺乏针对性，贸然进行煤层开采会导致作业面突水事故发生。
一些煤矿采掘企业缺少对周边环境的勘察，没有标记特殊地貌的位置、范围、积排水量、水情变化等详细数据，导致实际开采过程反复打眼钻孔确认安全方位，严重拖延工作面采掘进度，影响煤矿开采效率。
1.2重视程度不足煤炭企业最注重的就是效益，但是进行水文地质勘探会消耗大量的人力、物力、财力，且需要应用多种先进技术，会加大煤矿开采的成本。
因此，很多煤矿企业都会根据实际情况简化煤矿地质工程勘察工作，若在平坦区域进行煤矿开采，煤矿企业可能都不进行水文地质勘探。
而这些问题主要是由不重视水文地质勘探工作造成的，不仅会影响煤矿开采的安全性，也会影响煤炭企业的经济效益。
1.3缺少专业技术人才在煤矿实际开采工作中，很少有企业成立专门的防治水专业化部门，并针对矿区防治水形势提出整体解决策略。
从思想上看，煤矿开采企业不重视勘测工作，将其简单归类为辅助生产部门，缺少专业化的勘测队伍。
从管理角度看，煤矿勘测工作任务量难以量化，缺少相关管理制度考评，导致一些人员工作内容繁多而工作落实的实际效果不佳，在基层管理队伍中普遍存在畏难情绪，导致矿区勘测工作进展缓慢。
地下水环境影响评价工作常见问题及思考
地下水环境影响评价工作常见问题及思考王晓曦（辽宁省环境规划院有限公司.辽宁沈阳110165）摘要:通过调查分析公开环境影响评价资料，利用水文地质学原理，发现地下水环境影响预测工作中存在预测方法选择错误、水文地质条件分析不足、地下水流数学模型构建错误、数值模型未经充分识别等常见问题。
在今后的地下水环境影响评价工作中，建议加强水文地质条件和地下水流动机理分析，以利于提高预测评估的可靠性。
关键词:地下水环境影响评价;解析法;数值法;水文地质条件Abstract:Based on the investigation and analysis of environmental impact assessment reports and the principle of hydrogeology,founding that there are some common problems in the prediction of ground w ater environmental impact,such as the wrong selection of prediction methods,insufficient analysis of hydrogeological conditions,wrong construction of mathematical model of groundwater flow,and insufficient identification of numerical model.It is suggested to strengthen the analysis of hydrogeological conditions and groundwater flow mechanism in the future groundwater EIA work to improve the reliability of prediction and assessment.Key words：ground w ater environmental impact assessment;analytical method;numerical method;hydrogeological cond让i ons中图分类号：X820.3 文献标识码:A文章编号:1674-1021（2021）03-107-041引言环境影响评价制度是我国生态文明建设的一项重要法律制度。
《浙江省水文管理条例》广场宣传日活动成功举行
表１ 各工 况渗流 指标 统计 表
［ ３ ］张蔚榛 ． 地下水非稳定流计算 和地下水资源 评价 ［ Ｍ］ ．北京 ：
科学 出 版社 ，１ ９ ８ ３ ．
［ ４ ］陈崇希 ，林 敏 ．地下水 动力学 ［ Ｍ ］ ．武汉 ：中国地质 大学 出
版 社 ，１ ９ ９ ０ ．
［ ５ ］陶月赞 ，曹彭强 ，席道瑛 ．垂 向入渗与河渠边界影响下潜水 非 稳定流参数求 解 ［ Ｊ ］ ．水利学报 ，２ ０ ０ ６ ，３ ７（ ８ ） ：９ １ ３ —９ １ ７ ．
上 游 调 配 径 流 量 时 空分 配 。
响研究 ［ Ｊ ］ ．南水 北调与水利科技 ，２ ０ ０ ７ ，５（ ５ ） ：１ ０ ９ — １ １ ２ ． ［ ９ ］南江江 ．汾河太原城区段 治理美化二期工程南延伸段 蓄水池 ：１ ５—１ ７ ．
河道断流 时，因河道 渗漏 ，水 位逐渐 下降 ，如 同降 水头 渗
透试验 ，河水位 与周边 环境地 下水 位水头 差逐 渐减 小 ，河 道 内水量逐渐 漏光 。这 种情 况可采 用无 补给模 型 ，计算 河
［ ８ ］侯宏冰 ，刘长礼 ．黄河下游悬河温县段河水位变化对侧 渗的影
水位的变化 ，可计 算河道 水 蓄满 １次 ，水 面能 够维 持多 长 时 间；对 于枯水 期河 道水渗 漏光 ，通 过上 游水库 调配径 流 量间歇性补 给河道 蓄水 ，可采用定 流量 渗漏模 型 ，可计 算
对于径流时空 分布不 均 的河道 ，丰水期 ，镇 区段 河 道
［ ６ ］王旭东 ，殷 宗泽，宰金 珉 ．有 限区域 地下 水非 稳定 流解 析解
满蓄 ，河水位 达 到设 计 高 程 ，此 时，可采 用 定 水 头模 型 ，
主要 评价地下 水影 响边界 。对 于径 流稳定 期 ，采用 定流 量 补给模 型 ，可用 于计 算 河道 上游 补 给流 量 。对于 枯 水期 ，
潜水稳定井流的剖面二维数值模拟方法
第42卷 第4期2023年 7月 地质科技通报B u l l e t i n o f G e o l o g i c a l S c i e n c e a n d T e c h n o l o g yV o l .42 N o .4J u l . 2023王旭升,谢永桦,陈崇希.潜水稳定井流的剖面二维数值模拟方法[J ].地质科技通报,2023,42(4):27-36.W a n g X u s h e n g ,X i e Y o n g h u a ,C h e n C h o n g x i .S e c t i o n a l 2D n u m e r i c a l m o d e l l i n g m e t h o d f o r s t e a d y s t a t e w e l l -f l o w i n a n u n c o n -f i n e d a q u i f e r [J ].B u l l e t i n o f G e o l o g i c a l S c i e n c e a n d T e c h n o l o g y,2023,42(4):27-36.基金项目:国家自然科学基金项目(41972263;41772249)作者简介:王旭升(1974 ),男,教授,博士生导师,主要从事地下水动力学和水文模型研究工作㊂E -m a i l :w x s h @c u gb .e d u .c n ©Ed i t o r i a l O f f i ce of B u l l e t i n o f G e o l og i c a l S c i e n c e a n d T e ch n o l o g y .T hi s i s a n o pe n a c c e s s a r t i c l e u n d e r t h e C C B Y -N C -N D l i c e n s e .潜水稳定井流的剖面二维数值模拟方法王旭升1a,谢永桦1b,陈崇希2(1.中国地质大学(北京)a .水利部地下水保护重点实验室(筹);b .水资源与环境学院,北京100083;2.中国地质大学(武汉)环境学院,武汉430078)摘 要:经典D u p u i t 井流模型与考虑入渗的改进D u p u i t 井流模型,都受到D u pu i t 假定的影响,可能存在系统误差㊂建立反映三维流或轴对称二维流特性的潜水井流数值模型,是检验D u p u i t 模型可靠性的重要手段㊂提出一种模拟潜水稳定井流的剖面二维数值模拟方法,通过参数转换把柱坐标系的渗流方程变换为等效的直角坐标系方程,利用MO D F L OW 方体网格有限差分模型实现剖面二维流场模拟㊂数值模型以抽水井的井中水位为已知条件,渗出面的排水量按照D a r c y 定律差分公式计算,潜水面则通过MO D F L OW 对干-湿单元的处理加以识别,抽水流量经水均衡计算得到㊂通过采用精细化网格建立典型案例模型,获得模拟精度很高的结果,使反算抽水井流量的相对误差不超过0.2%㊂以此检验D u p u i t 井流模型,发现解析公式得到的水位线总体与数值模拟结果一致,仅在抽水井附近由于没有考虑水跃而偏低,且误差受到含水层渗透系数各向异性的影响㊂在有入渗的情况下,分水岭附近的渗流违反D u p u i t 假定㊂然而,改进的D u p u i t 井流公式计算的分水岭水位相对误差低于1%㊂这一数值模拟方法简单实用,但也受到MO D F L OW 本身局限性的约束㊂关键词:D u pu i t 假定;抽水井;入渗补给;渗出面;潜水面;各向异性;有限差分法2023-01-13收稿;2023-04-14修回;2023-04-19接受中图分类号:P 641 文章编号:2096-8523(2023)04-0027-10d o i :10.19509/j .c n k i .d z k q.t b 20230024 开放科学(资源服务)标识码(O S I D ):S e c t i o n a l 2D n u m e r i c a l m o d e l l i n g m e t h o d f o r s t e a d ys t a t e w e l l -f l o w i n a n u n c o n f i n e d a qu i f e r W a n g X u s h e n g 1a,X i e Y o n g h u a 1b ,C h e n C h o n gx i 2(1a .K e y L a b o r a t o r y o f G r o u n d w a t e r C o n s e r v a t i o n o f M i n i s t r y o f W a t e r R e s o u r c e s (i n p r e pa r a t i o n );1b .Sc h o o l o f W a t e r R e s o u r c e s &E n v i r o n m e n t ,C h i n a U n i v e r s i t y o f G e o s c i e n c e s (B e i j i n g),B e i j i n g 100083,C h i n a ;2.S c h o o l o f E n v i r o n m e n t S t u d i e s ,C h i n a U n i v e r s i t y of G e o s c i e n c e s (W u h a n ),W u h a n 430078,C h i n a)A b s t r a c t :[O b je c t i v e ]B o t h t h e c l a s s i c a l D u p u i t m o d e l a n d t h e m o d if i e d D u p u i t m o d e l i n c l u d i ng i n f i l t r a t i o n a r e i n f l u e n c e d b y D u p u i t s a s s u m p t i o n a n dh a v e p o t e n ti a l s y s t e m a t i c e r r o r s .B u i l d i n g a nu m e r i c a l m o d e l o f w e l l f l o w i n a n u n c o n f i n e d a q u i f e r b y c h a r a c t e r i z i n g t h e t h r e e -d i m e n s i o n a l o r a x i s ymm e t r i c t w o -d i m e n s i o n a l (2D )f l o w i s a n e s s e n t i a l a p p r o a c h t o v e r i f y t h e p e r f o r m a n c e o f D u p u i t -t y pe m o d e l s .[M e t h o d s ]I n t h i s s t u d y ,a 2D n u m e r i c a l m o d e l i s p r o p o s e df o r t h e s t e a d y s t a t e w e l l -f l o w i n a n u n c o n f i n e d a qu i f e r ,i n w h i c h t h e c o n t r o l e q u a t i o n o f s e e p a g e i n t h e c y l i n d r i c a l c o o r d i n a t e s i s t r a n s f o r m e d e q u i v a l e n t l y to t h e C a r t e s i a n c o o r d i n a t e s t h r o u g h p a r a m e t e r t r a n s f o r m a t i o n ,a n d t h e s e c t i o n a l 2D m o d e l l i n g i s i m pl e m e n t e d v i a t h e MO D F L OW f i n i t e -d i f f e r e n c e g r i d o f c u b i c b l o c k s .I n t h e n u m e r i c m o d e l ,t h e w a t e r l e v e l i n t h e p u m p i n gCopyright ©博看网. All Rights Reserved.h t t p s://d z k j q b.c u g.e d u.c n地质科技通报2023年w e l l i s a g i v e n c o n d i t i o n,t h e f l u x a c r o s s t h e s e e p a g e f a c e i s e s t i m a t e d b y d i f f e r e n c e e q u a t i o n a c c o r d i n g t o D a r c y s l a w,t h e p h r e a t i c s u r f a c e i s i d e n t i f i e d b y t h e t r e a t m e n t o f d r y a n d w e t c e l l s i n MO D F L OW,a n d t h e p u m p i n g r a t e i s d e t e r m i n e d f r o m t h e w a t e r d e b u g c a l c u l a t i o n.[R e s u l t s]F i n e g r i d s a r e c o n s t r u c t e d i n-n u m e r i c a l m o d e l s o f t y p i c a l c a s e s t o o b t a i n h i g h-p r e c i s i o n r e s u l t s,i n w h i c h t h e r e l a t i v e e r r o r o f t h e b a c k-w a r d s e s t i m a t e d p u m p i n g r a t e i s n o m o r e t h a n0.2%.T h i s n u m e r i c a l m o d e l i s u s e d t o c h e c k t h e D u p u i t-t y p e w e l l-f l o w m o d e l s.A s i n d i c a t e d,t h e g r o u n d w a t e r l e v e l e s t i m a t e d f r o m t h e a n a l y t i c a l f o r m u l a s g e n e r-a l l y a g r e e s w e l l w i t h t h e n u m e r i c a l m o d e l l i n g r e s u l t s,e x c e p t t h a t n e a r t h e w e l l,w h e r e t h e a n a l y t i c a l s o l u-t i o n u n d e r e s t i m a t e s t h e g r o u n d w a t e r l e v e l d u e t o i g n o r i n g t h e w a t e r j u m p a n d t h e e r r o r s d e p e n d o n t h e a n i-s o t r o p i c p e r m e a b i l i t y o f t h e a q u i f e r.W h e n i n f i l t r a t i o n e x i s t s,t h e f l o w i n t h e v i c i n i t y o f w a t e r s h e d d o e s n o t f o l l o w D u p u i t s a s s u m p t i o n.H o w e v e r,t h e e s t i m a t e d g r o u n d w a t e r l e v e l o n t h e w a t e r s h e d b y m o d i f i e d D u-p u i t w e l l-f l o w e q u a t i o n h a s a l o w l e v e l o f r e l a t i v e e r r o r,w h i c h i s l e s s t h a n1%.[C o n c l u s i o n]T h i s n u m e r i-c a l m e t h o d i s s i m p l e a n d p r a c t i c a l h o w e v e r i t i s a l s o i n f l u e n c e d b y l i m i t a t i o n s i n MO D F L OW.K e y w o r d s:D u p u i t s a s s u m p t i o n;p u m p i n g w e l l;i n f i l t r a t i o n r e c h a r g e;s e e p a g e f a c e;p h r e e t i c s u r f a c e;a n i-s o t r o p i c;f i n i t e d i f f e r e n c e m e t h o dR e c e i v e d:2023-01-13;R e v i s e d:2023-04-14;A c c e p t e d:2023-04-19地下水动力学中的井流模型具有轴对称特征,使用柱坐标系进行描述㊂经典D u p u i t井流模型[1]引入D u p u i t假定,把地下水流简化为单纯的径向一维流问题,便于获得控制方程的解析解㊂D u p u i t井流模型忽略了降水入渗的影响㊂为此,有研究者提出了改进模型[2-4],将入渗补给考虑在径向一维流问题的控制方程中㊂这些模型能够提供抽水井周围潜水面形态的理论公式㊂轴对称状态的潜水井流实际属于柱坐标系下的二维渗流过程,采用D u p u i t假设是将其简化为径向一维流问题的手段㊂在含水层的局部地段(如地下分水岭处),地下水的垂向流速可能显著大于水平流速,不满足D u p u i t假定㊂那么,如何检验D u p u i t假定对潜水稳定井流问题的有效性?这是地下水动力学中一个受到关注的问题㊂2020年陈崇希[4]在提出考虑入渗的改进D u p u i t模型时,也指出所建立的方程包含D u p u i t假定,可能在地下分水岭和抽水井附近带来误差,但只做了简单定性讨论,未定量分析D u p u i t假定对解析结果的影响㊂检验D u p u i t假定的影响,需要求解柱坐标系的二维潜水井流问题㊂不过,由于包含潜水面这样的自由面边界,获得解析解的难度很大㊂潜水井流问题还涉及另一个需要考察的实际因素,就是渗出面㊂地下水流向抽水井时,井壁水位高于井中水位(即水跃),介于两者之间的一段井壁就是渗出面㊂这是地表水等势体附近流线和等水头线弯曲状态(不满足D u p u i t假设)所产生的必然结果,其存在性与流场是否稳定无关[2,5-7]㊂前苏联学者已经证明[5-6]:即使存在水跃,D u p u i t模型的流量公式仍然是可靠的(无入渗补给条件),尽管在抽水井附近实际地下水位高于D u p u i t公式计算的结果㊂对于潜水井流的水跃现象,B o u l t o n[7]利用L a p l a c e方程的有限差分解法进行了系统的研究㊂张有龄[8]在阐述潜水井流问题和水跃现象时,详细引用了杨式德所做的有限差分模型,其方法与B o u l t o n所用方法类似㊂其他学者也用有限差分法[9]或有限单元法[10]对此做了数值模拟研究㊂由于非线性带来的困难,目前尚未研究得到水跃的严格解析解㊂前人所做的潜水井流数值模拟没有考虑降水入渗补给,而且需要采用柱坐标系进行网格离散,不容易在现有的通用数值模拟工具中实现㊂如果能够在直角坐标系建立与柱坐标系等价的剖面二维渗流数值模型,将为模型网格离散以及处理非均质含水层㊁入渗补给㊁渗出面等水文地质要素带来便利㊂目前国际上流行的三维地下水流有限差分模拟程序MO D F L OW[11]采用的是直角坐标系,可在一些专业软件进行快速建模㊂笔者拟利用控制方程的等效原理,在直角坐标系下采用等效方法模拟潜水稳定井流剖面二维流场,与D u p u i t模型或改进D u p u i t 模型解析解进行对比,对有关问题进行讨论㊂1直角坐标系模拟轴对称渗流方法1.1两种坐标系的渗流方程在直角坐标系中(图1-a),二维地下水稳定渗方程[12]的一般形式为:∂∂x K x∂H∂x+∂∂z K z∂H∂z=0(1)式中:H为水头;x,z分别为水平与垂直坐标,而K x和K z分别为对应坐标方向的渗透系数㊂该方程成立的条件是渗透张量的主轴与坐标轴一致,且流场内部没有源汇项㊂潜水面作为二维流场的上边82Copyright©博看网. All Rights Reserved.第4期王旭升等:潜水稳定井流的剖面二维数值模拟方法界,其稳态约束条件可表示为[12]:K x∂H ∂x ∂z w t∂x -K z ∂H ∂z+w =0(2)式中:z w t 为潜水面的高度(随x 坐标变化);w 为直角坐标系入渗补给强度㊂在柱坐标系中(图1-b ),轴对称渗流方程[13]的一般形式为:1r ∂∂r K h r ∂H ∂r +∂∂z K v ∂H∂z=0(3)式中:r 为径向距离;K h 和K v 分别为水平和垂直方向的渗透系数㊂这种柱坐标系下的潜水面稳态边界条件表示为:K h r ∂H ∂r ∂z w t∂r -K v r ∂H ∂z+εr =0(4)式中:ε为柱坐标系下的入渗补给强度㊂1.2等效变换方法已有一些研究者提出了在直角坐标系等效模拟轴对称地下水流的方法㊂S a m a n i 等[14]将直角坐标系的水平距离变换为x =l n (r )建立MO D F L OW 模型,然后按照这种对数变换,保持水平渗透系数不变,而对垂向渗透系数实施变换,产生等效结果㊂L a n ge v i n [15]提出另外一种在MO D F L OW 中转换参数的方法,他发现MO D F L OW 用对数加权平均法计算的块体单元界面等效渗透系数正好与D u -p u i t 公式的效果一致㊂L o u w y c k 等[16]建立x 方向均匀划分的MO D F L OW 模型,然后通过参数变换获得等效模拟结果,将水平渗透系数变换为:图1 直角坐标系(a )和柱坐标系(b )的潜水面示意图(z 为相对高度;x 为横向距离;r 为径向距离;后同)F i g .1 S c h e m a t i c s o f t h e p h r e e t i c s u r f a c e i n t h e c a r t e s i a n (a )a n d c y l i n d r i c a l (b )c o o r d i n a t e s ys t e m s K x =2πK hl n (r j +1/2/r j-1/2)(5)式中:j 为单元体在x 方向的离散坐标;r j +1/2,r j -1/2分别为单元外侧和内侧的径向距离㊂这种方法实际上也是取x 与l n r 成正比㊂本研究从控制方程等价的角度建立变换关系,即保持直角坐标系的水平距离x 与柱坐标系的径向距离r 相同,通过参数变换,使式(1)与式(3)等价,而式(2)与式(4)等价㊂满足这种等价关系的变换式为:K x =2πx K h Δy ,K z =2πx K v Δy ,w =2πx εΔy,x =r (6)式中:乘积因子2πx 为径向距离r 处完整过水断面的周长;肀y 为模型在x ,z 正交方向的宽度㊂在式(6)中变量与肀y 相除,可以保持参数的量纲不变㊂从物理意义上来讲,肀y 是剖面模型所代表的过水断面单位宽度㊂MO D F L OW 模型默认单元都是三维方体,因此必然会有y 方向的宽度肀y ,将其设定为一个任意单位宽度(例如1m )代入式(6)即可㊂L a n ge v i n [15]提出的方法本质也是如此,但他没有从控制方程等价的角度进行解释,也没有处理渗出面㊂1.3渗出面的处理渗出面是地下水流的特殊边界,在忽略流速水头的情况下,其水头与渗出点的位置高度相等㊂对于有入渗补给的潜水稳定井流模型,渗出面既可能出现在抽水井的井壁,也可能出现在外周含水层与大气的交界面上㊂假设井壁和含水层外周均沿铅直方向延伸,分别位于r =r w (井半径)和r =R (含水层外周半径)处,则渗出面的边界条件可表示为:H (x ,z )=z ,x =r w ,h w <z <(h w +L w )(7)H (x ,z )=z ,x =R ,h R <z <(h R +L R )(8)式中:h w 为抽水井中的水位;h R 为含水层外周水位,而L w 和L R 分别为井壁与含水层外周渗出面产生的水跃值㊂式(7,8)表明渗出面属于已知水头的边界㊂然而,在实际的模拟过程中,预先并不知道L w 和L R 的数值,因此难以在模型中直接设置某个高度段为已知水头边界㊂为了能让模型识别出渗出面,可以改用以下描述方式:V x =-K x ∂H∂x,H (r w +δx ,z >h w )>z 92Copyright ©博看网. All Rights Reserved.h t t p s://d z k j q b.c u g.e d u.c n地质科技通报2023年或H(R-δx,z>h R)>z(9) V x=0,H(r w+δx,z>h w)ɤz或H(R-δx,z>h R)ɤz(10)式中:δx>0,表示偏移一小段距离㊂采用偏移小段距离的水头大小来判断是否处于渗出面附近,符合数值模型中边界单元的特性㊂以井壁附近的单元为例,如图2-a所示,第k层和k+1层的相邻单元(宽度为肀x)计算水头分别为H(k)和H(k+1),它们均与x=r w+0.5肀x处的潜水面高度z w近似相等,但是只有第k层的水头满足H(k)>z(k)而处于渗出面附近㊂这里有δx=0.5肀x㊂图2-a中,第k层单元与渗出面接触,并且排泄的侧向流量Q x按照D a r c y定律计算,即将式(9)的差分格式与单元界面的面积相乘,计算为:Q x=-K x2(H k-z k)ΔxΔzΔy,H k>z k(11) Q x为负值表示流向与x轴相反㊂图2-a中第k+1层单元格点高于潜水面,符合式(10)所述条件,即Q x=0㊂在含水层外周,可以证明存在类似关系㊂因此,可以将渗出面的边界条件替换为第三类边界,即: Q d=C d(H k-z k),H k>z k;Q d=0,H kɤz k(12)式中:Q d为模型单元向外部排泄的流量;C d称为模型单元的排水系数㊂根据式(11),并把式(6)代入,可得:C d=2K xΔzΔyΔx=4πx K hΔzΔx,x=r w或x=R(13)a.井壁渗出面相邻单元;b.井孔饱水带相邻单元;r w为抽水井半径;Δx为单元长度;Δz为单元厚度;Δy为单元宽度;Q x为侧面流量;z w为井中水位;z(k)为单元中点高度;H(k)为单元格点水头图2含水层网格单元与井孔的关系示意图F i g.2 S c h e m a t i c s o f t h e r e l a t i o n s h i p b e t w e e n g r i d b l o c k s o f t h e a q u i f e r a n d t h e w e l l显然,排水系数的取值与模型网格剖分有关,并非自然界含水层的固定参数㊂实际上,井中水位以下的井壁(水头值为h w),也可以用类似的方法处理㊂如图2-b所示,设第k层有一个模型单元与井孔单元相邻,其宽度为肀x,2个单元交界面上的水头固定为h w,则侧向流量的差分公式为:Q x=-K x2(H k-h w)ΔxΔzΔy(14)式(14)与处理渗出面的式(11)是类似的㊂因此也可使用类似式(12)的公式计算单元的排泄流量:Q d=C d(H k-h w)(15)其中排水系数C d的计算方法与式(13)相同㊂2模型设计与MO D F L OW用法2.1含水层模拟情景本研究模拟均质含水层中的潜水井流问题,以便与D u p u i t模型和改进的D u p u i t模型[4]解析解进x,r为距离;R为圆岛半径;L R为外侧边界水跃;H(r,z)为水头分布函数;h w为井中水位;L w为井壁水跃图3潜水井流模型示意图F i g.3 D i a g r a mm a t i c m a p o f t h e m o d e l f o r t h e w e l l-f l o w i n a n u n c o n f i n e d a q u i f e r行对比㊂模型如图3所示,有一半径为R=2000m 03Copyright©博看网. All Rights Reserved.第4期王旭升等:潜水稳定井流的剖面二维数值模拟方法的圆岛状均质含水层,外周边界水位为h R=50m,圆岛中心处有一个半径为r w=0.1m的抽水井,以恒定流量Q w抽水,最终形成井中水位h w㊂含水层从隔水底板到地面的总厚度为70m,地下水接受大气降水入渗补给,强度为ε,最大可达1mm/d㊂地下水流场达到稳定,形成非线性的潜水面形态,并可能在抽水井的井壁发育长度为L w的渗出面(即水跃值为L w),在含水层外周边界发育长度为L R的渗出面㊂按照含水层渗透系数与入渗补给强度的变化,设置不同情景进行模拟㊂情景编号和参数见表1㊂表1模拟情景编号和参数T a b l e1 C o d e a n d p a r a m e t e r s o f m o d e l l i n g s c e n a r i o s模拟情景水平方向渗透系数K h/(m㊃d-1)各向异性比值K v/K h入渗补给强度ε/(mm㊃d-1)井中水位h w/mA1010.035B1011.035C100.1,10.01.035D10.11.020 2.2离散网格模拟工作是在A q u a v e o公司发布软件G r o u n d-w a t e r M o d e l i n g S y s t e m(GM S)中完成的,构建三维方体有限差分网格(y方向只取一个单元宽度),调用MO D F L OW模块处理水文地质要素㊂模型网格的初始剖分方式为:在垂向上分为70层,每层厚度1m;水平方向也进行密集划分,总体每个网格单元的宽度为2m,抽水井本身占一排宽度为0.1m的单元,其附近含水层单元加密,宽度缩小到0.01m㊂模型中井孔所在的这一列单元设置为无效单元,即并不直接参与有限差分计算,而是以边界的方式作用到模型中(具体处理见2.3节)㊂模型最右侧一排单元与含水层外周边界相邻㊂在模拟过程中,根据渗出面和潜水面分水岭的结果进行模型网格的微调,提高精度㊂每个含水层单元的水平渗透系数K x和垂向渗透系数K v,根据式(6)计算,其中x为单元形心与抽水井中心线的距离㊂2.3抽水井和外周边界的处理本次模拟并不像解析模型那样把抽水井处理为已知流量的边界,因为抽水井本身也被离散为很多个单元格(类似混合井),事先并不能准确掌握每个单元分配的流量㊂为此,本研究将抽水井的水位处理为已知条件,即给定h w的数值,然后在模拟形成稳定流之后反求出抽水井的流量㊂在井壁处,z坐标高于h w的点可能位于渗出面,其水头值与z坐标一致㊂与井壁相邻的一排模型单元,采用MO D F-L OW中的D r a i n模块[11]处理地下水向抽水井的排泄㊂高于井中水位的边界单元,排泄量用式(12)计算,以z k为D r a i n模块的排水高程㊂低于井中水位的边界单元,排泄量用式(15)计算,以h w为D r a i n 模块的排水高程㊂D r a i n模块的排水系数统一采用式(13)计算,取x=r w㊂调用MO D F L OW运行稳定流模拟之后,再根据模拟结果确定抽水井的流量Q w㊂它实际上由2个部分构成的:一部分是从渗出面排出的地下水量,一部分是从井水位下方井壁排出的地下水量㊂因此,抽水井流量的统计可以表示为:Q w=ðN k=1C d(H k-h w)+ðM k=N+1C d(H k-z k)(16)式中:N为井中水位以下的最高井壁单元层号;M 为对应渗出面的最高井壁单元层号㊂从模拟结果计算Q w的方法,是将井壁相邻的饱和带单元组合设置为一个均衡区,进行编码,调用MO D F L OW的B u d g e t模块计算均衡区流进和流出的水量㊂没有计算误差的情况下,两者的绝对值应相等,并且与Q w一致㊂如果潜水面存在分水岭,也可以将分水岭半径范围以内的全部入渗补给水量计为Q w㊂这2种计算方法的结果不一定相等,只要控制相对误差不超过1%,就可以认为数值模拟结果是可靠的㊂在模型的外周边界处,渗出面的情况与井壁渗出面是类似的,也用D r a i n模块处理㊂对高于水位h R的边界单元,排泄量用式(12)计算,以z k为D r a i n模块的排水高程,取x=R计算排水系数C d㊂低于h R的边界单元则将水头固定为h R㊂2.4潜水面的处理潜水面属于自由水面,各点的水头值与其高度相等,这可以用来判断潜水面的位置㊂与此同时,潜水面还需要满足式(2),作为渗流场的边界条件㊂MO D F L OW并不直接求解式(2),而是以入渗补给量作为源㊁建立潜水面所在网格单元的水均衡方程,与饱和带的方程组一起求解㊂从数值计算的角度来讲,这与求解式(2)是等价的㊂如果计算水头高于单元的底部,则该单元为 湿单元 ,是有效的㊂否则,就判断为 干单元 ,相当于无效单元,不再参加下一次计算,而是将潜水面移入下部的相邻单元[11,17]㊂MO D F L OW的R e c h a r g e模块可以自动将入渗补给加入到最高处的湿单元上㊂模型入渗强度w根据式(6)计算,其中x为单元形心的r值㊂模拟过程中,潜水面的形成实际上表现为网格单元的逐渐疏干过程㊂先将所有不属于边界的网格单元水头设置为模型顶部高程(70m)㊂然后启动渗流场计算和 干 ㊁ 湿 单元判断,反复迭代,逐渐将水头高于潜水面的单元疏干,而余下的网格单元13Copyright©博看网. All Rights Reserved.h t t p s ://d z k j q b .c u g.e d u .c n 地质科技通报2023年构成饱和带的渗流空间㊂迭代过程中可能会出现疏干过快或振荡现象,导致计算不稳定,因此需要对迭代算法进行控制,使疏干过程平稳进行㊂2.5求解算法和收敛精度MO D F L OW 形成有限差分方程组后,可采用多种方法进行求解㊂本次模拟采用预共轭梯度法(P C G 2模块[18]),这是一种迭代计算效果较高的算法㊂为保证数值模拟精度,P C G 2模块中水头计算的收敛标准设置为0.001m ,水均衡计算的收敛标准设置为0.01m 3/d㊂由于潜水面的形成是一个将包气带单元逐步疏干的过程,不宜太快,将P C G 2模块中迭代计算的松弛因子(r e l a x a t i o n p a r a m e t e r )从默认的1.0调整为0.5㊂3 典型模拟结果分析3.1对比D u pu i t 模型解析解表1中的情景A 是无入渗㊁均质各向同性含水层条件,剖面二维渗流模拟结果可以与D u p u i t 公式进行对比㊂图4给出了抽水井附近的模型网格与模拟结果㊂显然,潜水面的模拟水位高于解析解,在井壁处模拟得到渗出面的顶部高为40.5m ,比井中水位35.0m 高出了5.5m ,这就是水跃值㊂圆岛外周的水跃值小于1m㊂图4 情景A 模拟水位与经典D u pu i t 模型解析解的局部对比F i g .4 L o c a l c o m pa r i s o nb e t w e e n t h e m o d e l l e d g r o u n d -w a t e r l e v e l o f t h e sc e n a r i o -A a nd t he a n a l yt i c a l s o l u t i o n o f t h e c l a s s i c a l D u pu i t m o d e l 根据模拟结果,对抽水井流量进行了计算,得到Q w =4051.8m 3/d,其中从渗出面排泄的地下水量为364.2m 3/d ,占9.0%㊂模型的总体水均衡误差(流入㊁流出之差)绝对值小于0.01m 3/d,说明有很高的差分计算精度㊂解析方法计算流量的D u pu i t 公式[1]为:Q w =πK h (h 2R -h 2w )l n (R /r w )(17)把含水层参数和边界水位代入式(17),得到Q w =4044.6m 3/d㊂与之相比,数值模型得到的流量偏高0.2%㊂D u p u i t 模型的流量公式对三维流也是适用的,而且不受水跃的影响,这一点已经得到了严格的证明[5]㊂因此,上述0.2%的相对误差代表了数值模型在流量计算上的误差水平㊂3.2对比改进的D u pu i t 模型解析解表1中的情景B 是有入渗㊁均质各向同性含水层条件,剖面二维渗流模拟结果可以与改进的D u -pu i t 公式[4]进行对比(图5)㊂宏观上看,模拟的潜水面曲线与解析解很接近,在分水岭处几乎重合(图5-a)㊂在抽水井附近,解析解则明显低于模拟结果(图5-b ),渗出面产生的水跃值L w ʈ7.1m ㊂对于情景B ,抽水井流量与井中水位的解析关系为[4]:Q w =πK h (h 2R -h 2w )l n (R /r w )+επR22l n (R /r w )(18)把参数和边界条件代入式(18),得到抽水井流量的解析解为Q w =4679.0m 3/d ,高于无入渗的情景㊂数值模型反算的结果为Q w =4687.9m 3/d,其中渗出面流量510.7m 3/d ,占10.9%㊂模型的水均衡误差仍然小于0.01m 3/d,数值模拟得到的抽水流量比解析解偏大0.2%,相对误差很小㊂这种一致性也反映在分水岭上,解析解确定的分水岭在x =1220.4m 处,数值模拟的分水岭在x =1221m处,相对误差仅有0.05%,而且数值模拟的分水岭水位是50.52m ,与解析解一致㊂分水岭处地下水主要呈现垂向流动,明显违反D u p u i t 假设㊂然而,从本例对比结果来看,引入D u p u i t 假设并没有影响解析解计算分水岭水位的准确性㊂a .整体范围;b .抽水井附近图5 情景B 模拟水位与改进D u pu i t 模型解析解的对比F i g .5 C o m pa r i s o nb e t w e e n t h e m o d e l l e d g r o u n d w a t e r l e v e l o f sc e n a r i o -B a nd t he a n a l yt i c a l s o l u t i o n o f t h e m o d i f i e d D u pu i t m o d e l 3.3有入渗的各向异性含水层情景表1中的情景C 和D 都是各向异性含水层,可以考察含水层各向异性特征对潜水井流的影响㊂改23Copyright ©博看网. All Rights Reserved.第4期王旭升等:潜水稳定井流的剖面二维数值模拟方法进D u pu i t 模型[4]并没有包含对各向异性的处理㊂情景C 与情景B 的基本条件是相同的,唯一区别在于情景C 的垂向渗透系数与水平渗透系数不同㊂根据模拟结果(图6),各向异性对水位的影响主要表现在抽水井附近(图6-a ),而距离大于100m时的水位线对各向异性比(K v /K h )不敏感(图6-b)㊂分水岭位置和水位的模拟结果与情景B 几乎一致,说明各向异性对抽水井流量也几乎没有影响㊂从图6-a 可以看出,K v /K h 越小(即垂向渗透系数越低)则模拟水位线越高于改进D u pu i t 模型的解析解㊂当K v /K h =0.1时,水跃值为L w ʈ9.5m ㊂当K v /K h =10时,水跃值只有L w ʈ4.5m ㊂这说明:垂向渗透系数越大,解析公式得到的水位越逼近实际水位㊂情景D 的水平渗透系数为其他情景的1/10,而且井中水位比其他情景低了10m ㊂模拟结果见图7㊂与情景B 和C 相比,分水岭的位置(在r =643m )更加靠近抽水井,但水位更高(62.0m ),相应的抽水流量也降低到Q w =1287.0m 3/d,这主要是水平渗透性降低导致水流阻力增加的结果㊂这种情况下的水跃值也很大,达到L w ʈ33.5m ,使得在抽水井附近模拟水位远高于改进D u p u i t 公式的解析解㊂当r >117m 时,解析公式计算水位的相对误差低于1%㊂在分水岭处,解析公式计算水位为61.8m ,相对误差约0.4%㊂a .抽水井附近;b .全部范围图6 情景C 的模拟水位与改进D u p u i t 模型解析解的对比F i g .6 C o m pa r i s o nb e t w e e n t h e m o d e l l e d g r o u n d w a t e r l e v e l o f sc e n a r i o -C a nd t he a n a l y t i c a l s o l u t i o n of t h e m o d i f i e d D u pu i t m o d el 图7 情景D 的模拟结果F i g .7 N u m e r i c a l m o d e l l i n g re s u l t s of s c e n a r i o -D 4 讨 论4.1流速的换算如果要根据上述剖面二维渗流模拟结果绘制流线,需要注意2个问题㊂第1个问题是等效模型计算的流速与轴对称模型是否具有相同的方向?如果方向相同,那么流线的形态是一样的㊂第2个问题是两者计算的流速大小是否相同如果有差异,则意味着不能直接计算流线上质点的运移时间㊂下面对此进行讨论㊂轴对称坐标系与直角坐标系的D a r c y 流速分别为:V r =-K h∂H ∂r ,V v =-K v∂H∂z(19)V x =-K x∂H ∂x ,V z =-K z ∂H∂z(20)式中:V r 和V v 分别为轴对称坐标系的径向流速与垂向流速;V x 和V z 分别为直角坐标系的横向流速与垂向流速㊂把1.2节所提供的参数等效变换方法代入流速计算式,保持肀(∂H /∂r )=(∂H /∂x ),可建立2种坐标系D a r c y 流速的转换关系:V r =Δy 2πx V x ,V v =Δy 2πx V z (21)式(21)表明,2种坐标系D a r c y 流速是不同的,但V r /V v 与V x /V z 的数值是相同的,即流动方向不变㊂因此,2种坐标系的流线形态是相同的,而计算质点运移时间需要用式(21)换算流速㊂GM S 软件提供了MO D P A T H 模块与MO D F L OW 对接,可追踪粒子绘制流线,采用的是P o l l o c k 算法[19]㊂该算法33Copyright ©博看网. All Rights Reserved.h t t p s ://d z k j q b .c u g.e d u .c n 地质科技通报2023年计算粒子对流运移速度用到介质的孔隙度,即D a r -c y 流速与孔隙度的比值㊂因此,在调用MO D P -A T H 模块时,把实际孔隙度与式(21)中的比例系数(2πx /Δy )相乘作为等效孔隙度输入,即可得到符合轴对称坐标系的运移时间㊂图7展示了MO D P A T H 模块所绘制情景D 中的流线,且流线上相邻箭头之间表示一段相同的运移时间㊂潜水面上粒子的投放间距随径向距离的增大而增加,以保证任意2条相邻流线之间获得的补给流量为固定值㊂在分水岭下方,可以追踪出一根几乎垂直的流线,即为分水线(三维空间为分水面),其内侧区域的流线较为稀疏(反映抽水井流量占补给流量较小),流线上的箭头也较少(反映地下水流速较快)㊂4.2为何D u pu i t 假设未导致显著的分水岭水位误差在有入渗的情况下,只要抽水流量小于圆岛范围的入渗补给总量,就会出现分水岭㊂在分水岭处,地下水存在垂直向下的流速,而水平流速为零,显然违反D u p u i t 假定㊂如果按照直觉进行推测,D u p u i t 假定应该在分水岭处造成较大的误差㊂然而,在本研究的有入渗情景B ㊁C 和D 中,可用改进的D u p u i t 井流公式计算分水岭的水位,与数值模拟结果对比,其相对误差均低于1%,说明并没有受到D u p u i t 假定的显著影响㊂下面对于这种 意外 或 反常 现象做一个理论上的探讨㊂在剖面二维的潜水井流模型中,底部为隔水边界,因此垂向流速为零,而顶部潜水面处的边界条件由式(4)描述㊂如果忽略二阶小项,则流速的上㊁下边界条件可近似为:V v |z =h =-K v ∂H ∂z|z =h ʈ-ε,V v |z =0=0(22)式中:V v 为垂向D a r c y 流速,以向上为正;h 为潜水面相对底板的高度,即以底板为零高程时式(4)中的z w t ㊂尽管垂向流速随z 的变化是非线性的,我们可以根据边界条件用线性化公式来近似评估水力梯度,即:∂H ∂z ʈεK v zh(23)因此水头随z 变化的近似式为:H (z )ʈH z =0+ε2K v hz2或 H (z )ʈh -ε2K v h(h 2-z 2)(24)设分水岭处的水位为h d ,若采用D u p u i t 假定忽略垂向水力梯度,则底部水头的假设值为H z =0=h d ,而根据式(24)其实际值(取z =h d )近似为:H z =0ʈh d -εh d2K v(25)因此D u p u i t 假定计算含水层底部水头的相对误差为:e =|h d -H z =0|h d ʈε2K v(26)该误差水平取决于入渗强度与垂向渗透系数的比值㊂在案例情景B ,C 和D 中,ε=0.001m /d,K v =0.1~100m /d,因此理论上近似估计的相对误差不超过0.5%㊂如果要分析分水岭处水位的计算误差,则式(26)应改为:e =|h d -h D |h d(27)式中:h D 是采用改进D u p u i t 模型[4]计算的分水岭水位,即:h 2D=h 2R+εK h R 2-r 2d 2-r 2dl n R r d(28)其中r d 为分水岭的横向坐标㊂实际水位h d 目前还没有严格的解析解,所以与数值模拟结果对比是检验D u p u i t 假定的主要方法㊂为了掌握更普遍的情况,我们做理论上的近似分析㊂首先,从水均衡角度可以得到:Q (r )=2πʏhV h r d z =πε(r 2-r 2d ), r d ɤr ɤR (29)式中:Q (r )为穿过距离r 处圆环过水断面的流量(向外为正);V h 为水平方向的D a r c y 流速㊂把Da r -c y 定律形式V h =-K h ∂H /∂r 代入式(29),得到: ʏh∂H ∂r d z =dd r ʏh0H d z -h d h d r =ε2K h r 2dr-r, r d ɤr ɤR(30)在式(30)等号两边对r 积分,利用外侧边界条件但忽略渗出面,可得以下近似解:ʏhH d z -12h 2ʈ12h 2R +ε2K h R 2-r 22-r 2d l n R r,r d ɤr ɤR(31)把式(24)代入式(31),有:12-ε3K vh 2ʈ12h 2R +ε2K h R 2-r 22-r 2d l n R r,r d ɤr ɤR(32)把r =r d 代入式(32)得到分水岭处实际水位的近似解,即:43Copyright ©博看网. All Rights Reserved.。
对“再论‘单位涌水量就是导水系数’”一文的回应
对“再论‘单位涌水量就是导水系数’”一文的回应陈崇希【期刊名称】《水文地质工程地质》【年(卷),期】2017(044)005【总页数】1页(P172)【作者】陈崇希【作者单位】中国地质大学环境地质研究所,湖北武汉430074【正文语种】中文论文《再论‘单位涌水量就是导水系数’”》[1](下简称《单位涌水量》)中，质疑两本地下水动力学教材[2～3]，涉及的是基本理论问题。
《单位涌水量》指出：Thiem影响半径稳定井流模型是“假设地下水静水面为水平面，而在自然界几乎没有静水面是水平面的含水层”，这种模型、方程还能应用吗？笔者作为上述两本地下水动力学教材之一的主笔，对质疑做出回应。
(1)《单位涌水量》质疑的是影响半径稳定井流模型、公式。
笔者主笔的地下水动力学(P97～P101)[3]强调无界含水层中不可能形成稳定井流, 即明确判定Thiem 的“影响半径稳定井流模型”是错误的, 不可与Dupuit圆形定水头边界模型——“圆岛模型”混为一谈[3]。
当今的地下水动力学教材中不应再引入“影响半径稳定井流”理论系统, 更不应将它与不稳定井流理论并列为教材的两个重点[4～6]。
(2)虽然《单位涌水量》研究的是稳定井流条件, 但提出了一个很好的问题，即：在自然界几乎没有水头面水平 (处于静水状态) 的含水层，地下水动力学教材理应建立符合自然状态的井流模型。
有的地下水动力学教材 [2]认为：自然界地下水水头面坡度一般很小, 尤其在平原区, 通常为千分之几到万分之几，因此可以近似使用假设天然水力坡度为零(静止)的含水层中的理想化井流模型。
如此论证是缺乏科学性的。
运动的地下水与静止的地下水属于质的差别，运动的快与慢属于量的差别，两者怎能混淆？再者，水力坡度的大小仅仅是影响地下水径流量的因素之一，还有含水层导水系数及含水层宽度等因素也影响它。
这就是《单位涌水量》对地下水动力学教材提出质疑──缺乏科学依据地“对客观存在的地下水径流量忽略不计”。
新疆渭干河流域地下水三维流数值模拟
新疆渭干河流域地下水三维流数值模拟黎 明1a,陈崇希1b,张明江2(1.中国地质大学a.研究生院;b.环境地质研究所,武汉430074;2.新疆地质矿产局第一水文地质大队,乌鲁木齐830091)摘 要:采用任意多边形网格有限差分法求解模型,根据地下水水头长观资料识别模型的水文地质参数,建立了含混合井孔的6层模拟系统的地下水三维不稳定流模型来刻画渭干河流域地下水流动的基本规律。
运用识别后的模型进行了地下水资源评价,预测了未来10a不同开采条件下的地下水动态变化,探讨了地下水资源的合理开发模式,提出在以蒸发排泄为主的条件下,应在保证生态用水的前提下最大限度地夺取无效蒸发。
关键词:混合井孔;地下水三维流模型;数值模拟;地下水资源评价;渭干河流域;新疆中图分类号:P641.2;P641.8 文献标识码:A 文章编号:1000-7849(2005)01-0074-05 渭干河流域位于天山南麓,塔里木盆地北缘,范围涉及新疆阿克苏地区的库车、新和及沙雅三县。
流域内多年平均年降水量为65mm,多年平均年蒸发量约为2000mm,为典型的暖温带大陆性干旱气候。
长期以来,该地区只片面注重地表水的开发而忽视了地下水的开发,对水资源不合理的开发利用,不仅浪费了大量的水资源,而且已造成了河流断流、植被生态环境恶化、土地盐渍化、沙漠化等生态环境问题,因此,准确评价并合理开发利用地下水资源具有重要意义。
近年来,数值模拟成为地下水资源—地质环境保护定量评价的主要手段[1]。
笔者将从新疆渭干河流域地下水的形成、循环规律入手,建立地下水三维流数值模型,进而利用该模型评价地下水资源,研究水资源的合理开发模式。
1 水文地质条件概况渭干河流域的地势总体为北高南低,自西北向东南倾斜。
依地貌特征分为两部分:北部是以却勒塔格山为主体的构造剥蚀低山丘陵,海拔高度为1100～2060m,山体主要呈东西向展布,由第三系上新统地层构成;中部和南部分别为渭干河冲洪积平原和塔里木河泛滥平原,海拔950～1000m,村庄相对密集,为主要农业区,广泛分布第四系沉积物;其边缘地带大部分为无人类活动的荒漠地带。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
流的模拟上尚不尽人意，特别是对于普遍存在的混合抽水井流量分
配、混合观测孔的水位、自流井的流量动态及降雨渠道等入渗补给 地下水的滞后性等模拟要素，或者缺少模拟功能或者存在明显的缺
陷，
在泉、地下水-河流补给/排泄等方面也存在待改进部分；
当然，MODFLOW采用的矩形网格剖分方法也要损失模拟精度和/或 增大运算量。
我们深刻体会到 ：数值模拟的核心是“防止模拟失
真，提高仿真性”。因此，努力分析流动机理并用于数值
模拟。这些成果已经集成为一个基于多边形网格的三维地 下水流有限差分模拟系统（简称PGMS，即Polygon-grid
finite-difference groundwater modeling system）。
的三维有限差分地下水流模型MODFLOW软件处理。
MODFLOW建议：“多层井的流量必须以某种方式人为地 分配给每一单层，……把井流量按每层的导水系数大小分配， 即 Qi /Qw=Ti /ΣT ” (1988、1996、2000)。
为了便于讨论，又不失其一般，我们以贯穿两个含水 层的混合抽水井为例进行讨论。如此，上式可表述为 Q1 /Q2=T1 /T2 MODFLOW对混合井的这种处理方法没有给以理论上 任何的分析、说明，缺乏理论依据，应用中与实际也不符 (陈崇希等, 1992)。ห้องสมุดไป่ตู้是因为 (陈崇希等,2003a,2007,2008; 黎明等,2008)：
所在(陈崇希等,1998b)。
由于自然界地下水普遍存在三维流(图4-1-3)，而井孔，特别是抽水井几乎 都是非点状（滤管）井，因此混合抽水井普遍存在。
§4.1 混合抽水井的模拟 4.1.1 概述
Ⅰ-Ⅱ 混 合 观 测 孔 Ⅱ 层 观 测 孔 Ⅰ 层 观 测 孔 Ⅰ-Ⅱ 混 合 抽 水 井 Q
H2 v1 H1 Q v2 v3 Q1
§4.10 §4.11 §4.12 §4.13 §4.14 §4.15 §4.16 §4.17
非完整河流与地下水间补给、排泄的模拟 大气降雨、河渠、湖库入渗补给滞后性的刻画 潜水蒸发排泄的处理 初始水头条件的模拟 含水系统海边界问题 数值模型中流速和流量的计算问题 岩溶管道-裂隙-孔隙三重介质地下水流模拟 地下水开采-地面沉降的模拟
图4-1-4 混合观测孔中的各分层的流量及漏斗曲线
笔者将包括上述分析、批评MODFLOW的稿件投到美国 WRR杂志,其某审稿者云：“MODFLOW并没有说上述比例关 系可用于混合观测孔”而否定此稿件。实际上，当抽水井的流 量Qw足够小时，不会改变Q2>0和Q1<0的情况,仍存在比值 Q1/Q2为负值的情况，即出现负值(Q1/Q2)等于正值(T1 /T2)的 荒谬。 从理论上说，一种正确的方法，当蜕化为其特殊条件时， 也必定是正确的。科学研究中，人们正是通过许多特殊来认识 一般的；反过来，通过特殊来鉴别某理论的正确与否。这是常
i 1,2,3 N ; m 2, , N z 1
（2-6-4） 这就是多边形网格地下水三维流方程。
引 言
地下水资源评价、管理的核心是地下水动态的预测，当前地下
水动态预测最有效的手段，甚至可以说唯一的手段，是数值模拟。
国际上流行的地下水模拟专业软件很多，其中以Visual MODFLOW 应用最广泛，它在可视化方面做得比较精緻*，然而该软件在地下水
第7章
数值模拟实例
三维流多边形（棱柱体）网格 有限差分方程的建立
我们已经学习过基于矩形网格的差分方法。不难看出，无论是等
格距还是变格距矩形网格，这种剖分是有局限性，通常不能满足实
际问题的需要。常见的实际问题大多是：含水层渗流区边界形状不 规则，非均质参数分区的界线不规则，抽水井和观测孔的布置一般 也是任意的。这些问题若用矩形网格剖分则很不方便,特别是井孔及 内部与地下水有水力联系的河流等源汇点/线。若将所有井孔及其它 源汇点都与格点一致,精度较高,但许多网格没有必要地被加密；若网 格密度合理,则可能井孔及其它源汇点不能与格点一致而丢失精度。 为此，我们介绍一种以辅助三角形剖分为基础的不规则网格差分法， 通常被称之为(任意)多边形网格法。
（1）含水层导水系数对井孔流量的影响， 不会如此简单。例如,混合井附近岩性(渗透系 数)发生变化，甚至混合井打在岩性透镜体上， 怎样影响井管的流量分配？含水层厚度发生变 化又如何改变流量的分配？
K1
K3
K2
混合井管贯穿岩性透镜体
图4-1-3 混合井管贯穿岩性透镜体上
（2）含水层的参数影响混合井流量的分配，导水系数只 是其中一个因素，含水层的弹性给水度（储水系数）就不起作
若记Qwi为其源汇项，包括开采井、泉流量、河流的补给、 排泄、降雨入渗等，则方程可以写为，
hn,m1 hin, m1 hkn,m1 hin, m1 hin, m11 hin, m1 Tij , m j pb Tik , m bq Ki , m 1 / 2 Ai zi , m 1 zi , m ij ik e hin, m1 hin, m11 hin, m1 hin, m Ki , m 1 / 2 Ai Qwi , m i , m Ai zi , m zi , m 1 / 2 tn 1 tn
开采地下水的井孔多采用混合井（除非水质不符合要求），
以简化成井工艺，增加出水量和降低成本。因此，几乎没有无
混合井的地区，特别是民间可以凿井的我国（包括台湾大量的 是混合井）。既使现在开始控制打井，控制开采，但以往的混
合抽水井转变为混合“观测孔”。
混合井在地下水流模型中如何刻画，国外尚未很好解决。 目前主要采用由美国地调局推出，并在国际上广泛流传和应用
致谢
PGMS模型的研制获得了中国地质调查局和水文地质环境
地质研究所地质大调查项目及多项国家自然科学基金项目的支
持; 广西地矿局科技处和北海地矿公司专门设置混合抽水试验
场为本项目提供现场试验及滨海区地下水潮汐效应的系统观測，
中国地质大学（武汉）渗流实验室为模型的验证试验提供了实 验室条件，甘肃地调院、新疆一水、甘肃二水,陕西二水和山东
用了吗?
（3）井管的流量分配不仅与含水层的水文地质参数分布有 关，还与外边界条件、井径、水泵吸水口的位置及其它抽水井 的干扰等有关，
（4）按MODFLOW的“流量预分法” ，一个混合抽水井 特别是预测时，该混合抽水井附近可加入或关闭 抽水井，在
的Q1/Q2 始终是个常量(因为T1/T2为常量)。然而在模拟过程中，
PGMS软件选择多边形网格有限差分数值模拟，就 剖分 而言，它
比MODFLOW 矩形网格有限差分更能适应复杂的边界条件，容易地
将水文地质诸模拟要素置于格点处和格线逼近水文地质界线；就精度 而言，与常用的有限元相当，且物理意义明确。
PGMS软件基本模型已与一、二、三维地下水流解析解做过模拟
对比，模拟结果良好。 PGMS软件已在近几年地质调查项目（新疆渭干河流域, 陕西渭 北岩溶水, 甘肃河西走廊疏勒河流域 及 黑河中游地区 等项目地下水 模中得到全面的应用与检验，个别水文地质模拟要素，如混合抽水井 和混合观测孔和降雨渠道等入渗补给滞后性等，则被应用更多次。 尽管现今的PGMS软件属1.0版,但已具有很好的应用前景。我们 将进一步发展、完善该软件。
M2 +M’’+M1 M2 +M’’ M2 +M’’+M1+M’
Z H0
潜水层 M’ 1弱透水层
M1 Ⅰ承压含水层 M’’ 2弱透水层
Q2 v4 Q2
M2
M2 Ⅱ承压含水层
0
4-1-1 混合抽水井、混合观测孔的原定义(陈崇希等,1992)
(x, y)
图4-1-1
混合井孔的传统定义
图4-1-3 自然界地下水普遍的三维流模式(陈崇希等,2005b)
目前看来,这种方法是能使边界、内部界线及井孔等源汇点
/线与网络格点、格线最佳逼近,而又尽量少增加格点数目的网
络格式,是一种比较实用的方法。因此,本教材特别推荐此法。 Istok(1989)认为:“有限差分法难以处理不规则边界、非均
质和各向异性介质或倾斜含水层”。把有限差分法限于矩形网
格是不妥的。 多边形网格有限差分法其差分方程的建立,我们采用直接根 据达西定律和水均衡原理建立差分方程。
（2-6-1）
Qe Tij ,m
式中：
ij
pb Tik ,m
ik
bq
根据 （2-6-1）式的计算方法，对m层格点周围的所有三角 形做类似的计算，并求和得
Q侧 Tij ,m
e
h n,1 hin,m1 jm
ij
pb Tik ,m
hkn,1 hin,m1 m
ik
bq
地下水三维流数值模拟 边界条件 源汇等特殊问题 的模拟与处理 陈崇希
中国地质大学（武汉） 环境地质研究所
《地下水流数值模拟》大纲
前言
第1章 地下水流定解问题概述
第2章 地下水流有限差分法
第3章 地下水流有限元法（限于迦辽金有限元法）
第4章 边界条件 源汇 等特殊问题的模拟与处理
§4.1 §4.2 §4.3 §4.4 §4.5 §4.6 §4.7 §4.8 §4.9 混合抽水井的模拟 自流井的模拟 水平井的模拟 混合观测孔水位的模拟 抽水井水位的校正方法 井孔表皮效应（井周扰动效应）的处理 面井内抽水井水位降深的计算 非完整抽水井附近观测孔井水位的校正 泉的模拟
第5章 地下水流模型的拓展
§5.1 地下水饱和-非饱和流耦合模型 §5.2 地下水-地表水流耦合模型 §5.3 分布式水文模型简介
第6章 数值模型设计 应用及水文地质勘查
§6.1 数值模拟设计应用主要步骤 §6.2 水文地质概念模型的设计
§6.3 数值模型的建立