1、Modflow 四小时速成四小时速成1 水流问题有限差分法2 Modflow3 Processing Modflow4 PM Path5 例题练习主要参考书目：1 地下水流的数学模型和数值法，孙讷正，1981.5，地质出版社。2 地下水动力学（第二版），薛禹群，1997.9，地质出版社。3 地下水流动问题数值方法，陈崇希等，1990.6，中国地质大学出版社。4 地下水流数值模拟，1989.9,李俊亭，地质出版社 教学中涉及软件：1 Processing Modflow2 Surfer(V.6或更高)3 Auto CAD 2000Modflow 4小时速成小时速成数值解与数值方法数值解数值解 用
2、有限个离散点的函数（水头）值，近似代替函数（水头）的连续分布，称数值解。数值解属于近似解，在地下水渗流计算中，是最“精确”的（概化误差最小）地下水计算误差地下水计算误差 概化误差(要命的误差)、测量误差、截断误差等数值方法数值方法 求解微分方程数值解的方法，常用方法有:有限差分法(FDM)、有限单元法(FEM)二维渗流问题矩形网格有限差剖分示意图二维渗流问题矩形网格有限差剖分示意图有限差剖分结点法有限差剖分结点法二维渗流问题矩形网格有限差剖分示意图二维渗流问题矩形网格有限差剖分示意图有限差剖分格点法有限差剖分格点法二维问题矩形网格有限差解之几何意义,1,111,1,111,1,1,221,1,
3、1,11,1,1,1,22 ()()()()i j m ki j m kS i j mkij m ki j m kij m ki j m kij mij mi jm ki j m ki jm ki j m ki jmi jmHHSx y ztx zkHHkHHyy zkHHkHHx 1,1,1,11,1,1,1,221,2 ()()x y zi j mki j m ki j mki j m ki j mi j mi j m kx ykHHkHHzW 三维问题的“简隐式简隐式”有限差分方程(Modflow用该格式)1 不仅计算含水层水位，还要计算弱透水层中的水位不仅计算含水层水位，还要计算弱透水
4、层中的水位 2 根据需要，含水层、相对隔水层还要进一步分成若干层根据需要，含水层、相对隔水层还要进一步分成若干层简隐式特点简隐式特点：1 差分格式无条件稳定与收敛；差分格式无条件稳定与收敛；2 时间步长随意；时间步长随意；3 只要时、空剖分步长只要时、空剖分步长“充分小充分小”，就足够，就足够“精确精确”4个下标含义1,111111/2,1,1/2,1,1111,1/2,1,1/2,1,1/2()()()()(kki j mi j mi j mkkkkkij mij mi j mij mij mi j mkkkki jmi jmi j mi jmi jmi j mi j mHHSxytxTHH
5、THHyyTHHTHHxk i i,j j,mmH H11,1/2,1,1/2,1,1,2,1,)()2 kki j mi j mi j mi j mi j mi j ki j mij mij mi j mij mHHkxyMMWxyTTTTT 调 和 平 均 导 水 系 数:相 邻 格 点 只 要 有 一 个 导 水 系 数 为 零，其“平 均 值”就 为 零这 一 特 点，可 满 足 自 动 隔 水 边 界 的 要 求！k k+1 1k k+1 1-1 1i i,j j,mm+1 1H H准三维问题的准三维问题的“简隐式简隐式”有限差分方程有限差分方程 1 仅求解各含水层（组）水位，弱透水
6、层不计算仅求解各含水层（组）水位，弱透水层不计算 2 一般情况下，单一含水层（组）不再进一步分层一般情况下，单一含水层（组）不再进一步分层(?)1 那些结（格）点需要求解那些结（格）点需要求解 内格点=需要（列出差分方程式）流量边界（第二类）格点=需要（列出差分方程式）水位边界（第一类）格点=不需要（因为水位已知啦）外格点=不需要2 通用差分方程的特点通用差分方程的特点 a 所有内格点的方程式都是相同的（给编写程序带来方便）b 当导水系数用调和平均值时（假设:计算区域之外的导水系数为零），隔水边界（零流量边界）的差分式与内结点相同。c 流量边界“等价于”隔水边界+“注(抽)水井”，注水井流量的
7、大小 等于边界流入量。有限差分方程中那些格点需要计算有限差分方程中那些格点需要计算1 概化水文地质条件数学模型2 离散化数值模型3 制备数据文件集合(BAS.dat,BCF.dat,Well.dat)4 运行模拟程序Modflow.exe5 输出模拟结果文件集合 Heads.dat,Drawdown.dat6 解读模拟结果(均衡分析、绘图、表等)前处理软件前处理软件协助自动或半自动制作数据文件集合，具有校对、图形显示、人机交互特征 后处理软件后处理软件具有自动整理结果数据文件为简洁易懂的图表 等功能数值模拟步骤前处理后处理地下水模拟程序 美国地质调查局(UGS):MODFLOW 地下水流模拟程
8、序 MT3D、MT3DMS溶质运移模拟程序集成前后处理程序 Visual ModFlow Processing ModFlow(PM)Groundwater Modeling System(GMS)Visual Modflow 数值模拟：常用商业软件Modflow剖分特点与格元编号Modflow格元均衡示意图Modflow格元侧面流量计算表达式垂向传导（越流）系数的计算三维流垂向传导系数Cv11211112212kkkkkkkkkkvvKVvvKVKVr cCVKVv “准三维流”垂向越流系数越流系数VCont无夹层VCont有夹层VCont1122kkCVVCr c Modflow通用差分方
9、程式通用差分方程式Modflow88模块化程序结构递推求解模式时间剖分方式时间剖分问题时间剖分问题1 总模拟期分成若干“应力期”2 每个“应力期”分解为若干时间段3 每个应力期内各模块包参数为 “常数”，不同应力期间参数可以 变化，以实现随时间动态变化。Processing Modflow集成前后处理系统介绍File文件菜单Processing Modflow集成前后处理系统介绍File/Preferences模拟程序设置Processing Modflow集成前后处理系统介绍Grid剖分菜单Processing Modflow集成前后处理系统介绍Parameters参数菜单Processin
10、g Modflow集成前后处理系统介绍Models模型菜单Processing Modflow集成前后处理系统介绍Tools前后处理工具菜单基本模型基本模型无任何外应力包无任何外应力包基本基本ModflowModflow模型模型，仅含BAS、BCF、SIP与SSOR能够模拟以下简单渗流问题能够模拟以下简单渗流问题1 形状任意（马赛克排列成任意近似形状）2 地层参数（渗透系数、导水系数、储水系数、给水度等）总体非均质，各格元内部均质，不同格元允许参数不同。3 稳定流、非稳定流4 仅允许含有隔水边界与定水头边界（永远不随应力期变化）井模块包井模块包Well源汇项强度与水位无关源汇项强度与水位无关在
11、基本Modflow模型的基础上，增加模拟井功能1 模拟抽、注水井。注水为正，抽水为负2 模拟流量边界 流量边界隔水边界抽（注）水井3 每个格元只能有一个抽（注）水井，其值应为位于该格元 各井流量的代数和4 包含信息 含井总格元数 各各应力期的：井位置、井流量（不同应力期可以不同）“面状补给面状补给”模块包Recharge(RCH)源汇项强度与水位无关源汇项强度与水位无关增加模拟“面状补给面状补给”功能1 模拟大气降水入渗、凝结水入渗等功能2 每格元只能有一个补给强度值，其值为位于该格元代数和3 包含信息 含面补总格元数？补给加载到哪一层？各各应力期：补给强度（注意，与抽水井流量单位不同）排水沟
12、模块排水沟模块包Drain(DRN)源汇项强度与水位源汇项强度与水位“有关有关”！1 可模拟排水沟、泉水溢出等功能（只出不进，类似二极管）（只出不进，类似二极管）2 只能向“排水沟”排，不能由排水沟进入含水层3 排水沟水力传导系数(Drain Hydraulic Conductance)描述阻力大小，近似模拟向 排沟溢出的通畅程度4 排水沟排水高程 Elevation of the Drain：泉口高程；排水沟内水位;4 包含信息 各应力期的：排水沟位置、哪层、排水高程、排水沟水力传导系数,()0 i j ki j ki j ki j ki j ki j ki j ki j ki j kQDC
13、DHDHDQDHD当当Drain Hydraulic Conductance在该格元所控制的“排水沟”中，单位水头差所溢出的水量 CdElevation of the Drain该格元所控制的“排水沟”，平均排水高程Di,j,k。所谓水头差，是指格元地下水位与该排水高程之差。,()0 i j ki j ki j ki j ki j ki j ki j ki j ki j kQDCDHDHDQDHD当当用逐格元输入法输入排水沟参数：Equivalent Hydraulic Conductivity排水沟在该节点所控制单位长度、单位水头差所溢出的水量 KElevation of the Drain
14、排水沟在该节点所控制排水高程 Di,j,k用折线输入法输入排水沟参数：,()()0 i j ki j ki j ki j ki j ki j ki j ki j ki j ki j ki j ki j ki j kQDKLHDCDHDHDQDHD当当,i j ki j ki j kCDKL自动化程度高！在某地下水模型中，选用了在某地下水模型中，选用了“排水沟排水沟”模块包模块包问题：模拟程序运行后，发现排水沟出水太多问题：模拟程序运行后，发现排水沟出水太多了，怎么办？了，怎么办？解决方法解决方法1 沟排水太通畅了，把系数沟排水太通畅了，把系数Cd或或K减小一点，把地减小一点，把地下水憋一下。或
15、者下水憋一下。或者2 沟的排水高程沟的排水高程 D 太低了，高点试试太低了，高点试试3 两种方法一起上，再试试！两种方法一起上，再试试！河流模块河流模块包River(RIV)_1,()i j ki j ki j ki j ki j ki j kQRivCRivHRivHHHRiv 当当地下水向河流排泄时,()()i j ki j ki j ki j ki j ki j ki j ki j ki j ki j ki j ki j kHBrivQRivCRivHRivHHBrivQRivCRivHRivBriv当当当河流向地下水渗漏时河流模块包River(RIV)_21 可模拟河流渗漏与溢出2 河
16、底水力传导系数为等效阻力，近似模拟水流的通畅程度3 河水位指河水水面高程4 包含信息 各应力期的：河流位置、河水高程、河底（等效）阻力层底面高程 河底阻力层等效水力传导系数用逐格元输入法输入River模块参数：,()i j ki j ki j ki j ki j ki j kQRivCRivHRivHHHRiv 当当地下水向河流排泄时,()()i j ki j ki j ki j ki j ki j ki j ki j ki j ki j ki j ki j kHBrivQRivCRivHRivHHBrivQRivCRivHRivBriv当当当河流向地下水渗漏时Criv格元单位水头差渗漏（溢出
17、）量用折线输入法输入排水沟参数：自动化程度高！蒸腾模块包Evapotranspiration(EVT)1 可近似模拟地下水浅埋带蒸发 与 植被根系吸收地下水蒸发功能2 涉及参数：蒸发面高程、极限蒸发深度、水面蒸发强度,max,max,1,0 ,i j kMaxi j ki j ki j ki j ki j kMaxi j ki j ki j ki j ki j kEvtEvtHHsHsHEvtEvtHsDHHsDEvtHsmax,i j kDH 0max1n 柯夫达阿维里扬诺夫经验公式,max,max,1,0 ,i j kMaxi j ki j ki j ki j ki j kMaxi j k
18、i j ki j ki j ki j kEvtEvtHHsHsHEvtEvtHsDHHsDEvtHsmax,i j kDHMaxEvt,i j kHsmaxD一般用外部场文件输入通用水位边界模块包General_Head_Boundary(GHB)1 可近似模拟变水位边界、与水头有关的变流量边界等2 涉及参数：边界外侧水位Hbi,j,k、边界水力传导系数CHBi,j,k,()i j ki j ki j ki j kQGBCHBHBH应用举例1 变水头边界2 边界流量随降深变化3 阻水构造（如阻水断层），边界流量 与内外水头差成正比CHBijk单位“格元水位与外部水位”差“诱发”的流入（出）量,
19、i j kCHB,i j kHBGHBKGHBHB,i j ki j ki j kCHBKL,()i j ki j ki j ki j kQGBCHBHBH注意：1 在定水头边界定水头边界上加任何“应力包”都是徒劳的，是白增加的。相当于在大海中“抽水”或“注水”！2 在定水头边界附近抽水，可“袭夺”定水头边界“取之不尽的水源”，而边界水位永远不会下降，这种情况自然界少见。3 总之，应慎用定水头边界！4 推荐：若边界即能提出水位边界，又能出提流量边界，则应首选流量边界。否则当出现边界量失真时不易察觉。四种模型层类型矩阵数据文件编辑外部矩阵数据文件输入外部矩阵数据文件输入（覆盖）加载矢量图形文件(
20、如Autu_CAD 文件)格元图例易于看懂加载各种模块分布PMPath主界面阻滞因子及其含义阻滞因子及其含义PMPath设置“粒子”分布与数量PMPath“粒子”迹线演示PMPath“粒子”迹线演示PMPath“粒子”迹线阶段标记设置Processing Modflow 前后处理器容量限制Processing Modflow 单层数据文件格式Processing Modflow 单层数据文件格式Processing Modflow 各种包含“包”集合Pressing Modflow Pro 例例 1（三维稳定流问题）问题描述：问题描述：由两种不同岩性构成潜水含水层，各向同性，分布范围为矩形。南
21、、北边界为隔水边界；东、西边界为河流定水头边界，水头分别为8m、9m。上下两层水平渗透系数分别为8.64m/d、43.2m/d，厚度分别为4m、6m，垂向渗透系数为水平方向的10%，地面高程为10m。含水层接受面状补给强度为6.91210-4 m/d。上层西边界附近有一污染源。任务：任务：(东界附近)用一完整抽水井来截获污染源，确定抽水流量，以水力隔离污染源。建立并运行稳定渗流数值模型步骤建立并运行稳定渗流数值模型步骤 1 创建模型 2 给模型赋值 3 运行模拟程序 4 观察校核模拟结果 5 计算水均衡 6 模拟成果输出剖分伊始，用一个大“盒子”把计算区套起来Pressing Modflow
22、Pro 例例 1计算结果显示Pressing Modflow Pro 例例 1计算结果显示Pressing Modflow Pro 例 2（单层非稳定流问题）概况：概况：单层潜水含水层，粗砂，渗透系数160m/d，各向同性，给水度0.06；顶底板高程分别为25m、0m；东西界宽6000m，为隔水边界；南北向10000m，北界为定水位边界，水位15m，南界为定流量边界，单宽流量0.0672m2/d；在湿润季节(4个月),入渗补给强度为7.510-4 m/d，干旱季节(8个月)无入渗；共9眼开采井，单井开采量3888m3/d，仅在干旱季节用于灌溉开采，湿润季节不开采；任务任务:1 建立近似稳定流模
23、型，无开采，年平均入渗补给强度为2.510-4 m/d；2 以“稳定”场作为初始条件，运行10年，近似得到“拟稳定”场；3 以“拟稳定”场的湿润季节末作为初始条件，模拟干旱季节开采8个月后的流场；4 以开采8个月后的流场作为初始条件，模拟经4个月补给后的流场。5 重复3、4五次，及模拟灌溉开采五年，与第一年开采流场进行比较。构建例2 “稳定流”模型剖分与加密剖分问题稳定解流场（无抽水天然情况）例 2 枯水期末、补给期末流场怎样添加“观测孔”怎样显示“观测孔”曲线怎样制作“动画”Pressing Modflow Pro 例 3（多层稳定流问题）概况：一概况：一山间河谷，两侧为隔水山体。共有两层含
24、水层，上部潜水与下部承压水，河水与潜水水力联系密切，承压含水层底板为相对不透水岩石，潜水与承压水之间有一相对弱透水层分隔。河流参数：河流参数：上游界河水位19.4m，下游17m；河宽100m；河底“淤积”层厚度1m，渗透系数2m/d；淤积层底界高程，上游17.4m，下游15m。潜水含水层参数：潜水含水层参数：水平与垂直渗透系数分别为5m/d、0.5m/d；给水度0.06；有效孔隙度0.20弱透水层参数：弱透水层参数：水平与垂直渗透系数分别为0.5m/d、0.05m/d；厚度2m承压含水层参数：承压含水层参数：水平与垂直渗透系数分别为2m/d、1m/d；储水率510-5d-1，有效孔隙度0.25
25、。抽水井抽水井：有三眼下部承压含水层抽水井，单井抽水量500 m3/d。边界条件：边界条件：各层顶底板高程、上下游边界水位等，都存在数据文件中，建模时直接调取。任务：任务：1 构建三层结构稳定流数值模型（含河流与抽水井）；2 确定(迹线显示)抽水井的“捕获带”范围；3 分析抽水井激发夺取河水量比例。主要步骤近似定义边界形状及边界特征估算河流参数（Conductance）承压含水层稳定流场分布及其水井“捕获”带“粒子”正向追踪法，确定“污染源”污染途经有限差所计算的格点水位有限差所计算的格点水位：格元内“平均”水位1 当格元内没有抽水井时，数值解格点水位就是格元形心点形心点 的水位。2 当格元内
26、有抽水井时，格点水位相当于某虚拟“大井大井”（面积井）的井壁水位，比抽水井内实际水位高格点水位、观测井水位、开采井水位2ln()44.810exp()2ffTHQxrHQTxxxxr稳定井流蒂姆公式:差分向井单元汇流式:联立求解得虚拟大井半径:,2,2,2,()()2.25()ln42.25 ()ln4/4.81 fwfwi jrri jwRrSi jwfswi jwwi jHHsQHsQsQRTtHsQRTrSQTtHHsQTxSHH以承压水抽水井为例:,取近似公式:则式中:开采井内实际水位开采井格点(有限差)计算()-fwwRrfsQs QR rsR rs水位为井内降深曲线函数，由多落程抽水试验获得。为-区间地层水头损失与井损为-区间地层降深，可用Jacob公式近似估算，近似取持续抽水t(一般取0.51天)时的降深。开采井水位估算方法-推荐106 以上有不当之处，请大家给与批评指正，以上有不当之处，请大家给与批评指正，谢谢大家！谢谢大家！
