1. 从核心需求出发为什么要在Python里跑Modflow6地球物理与水文地质圈这两年有个明显趋势大家不再满足于把MODFLOW当“黑盒”点两下界面出张等值线图而是想把建模流程纳入更完整的科学计算链路里。MODFLOW 6这个版本的最大变化是整个代码框架用面向对象思路重写底层的求解器、时间步进、边界条件都变成了一套可以“勾搭”进Python生态的模块。配合flopy工具库你几乎可以在Jupyter Notebook里完成从网格生成、参数赋值、运行模拟到读取结果的完整闭环。我在实际项目里最强烈的体感是过去用传统界面建模调一次参数要反复点菜单、导出、再导入现在只需要写好一份可复现的Python脚本改参数等于改变量批量跑情景模拟变得极其顺手。尤其是溶质运移这块MODFLOW 6把地下水流模块GWF和溶质运移模块GWT分开建模两者通过 exchange 机制耦合逻辑比旧版本清晰太多。这篇内容就是围绕一个典型的含承压含水层的溶质运移案例从建模思路、Python源码结构、参数计算到调试排错完整走一遍。不管你是刚接触数值模拟的研究生还是工作中需要用模型辅助决策的工程师按照这个流程都能把模型跑起来并且能看懂每一行代码在干什么。内容中涉及的所有脚本结构和参数设计均来自我在项目中的实际使用方案亲测可跑、可复现。2. MODFLOW 6的Python建模体系全景2.1 从传统建模到Python API的范式转变传统的MODFLOW建模流程大致是划分网格→分配水文地质参数→设置边界条件→运行计算→后处理。这套流程在GMS或ModelMuse里做当然没有问题但它最大的隐患是“过程不可追踪”。模型文件是谁生成的、参数是哪个版本改的、边界条件具体怎么设的时间一长很容易说不清。换成Python脚本之后整个过程变成了一条可版本管理的“代码生产线”。flopy作为USGS官方维护的Python工具包天然支持MODFLOW 6全部关键字你可以直接在脚本里new一个Simulation对象往里面塞GWF模型、GWT模型、离散化信息、应力期和输出控制。整份脚本本身就是模型说明书以后想回溯哪个版本改了什么都清楚这对我这种经常要跟多个项目并行的人来讲价值太大了。2.2 flopy与MODFLOW 6的模块化架构要写对源码先搞清楚MODFLOW 6的模块划分。顶层是Simulation负责管理整个模拟过程Simulation下面挂GWF模型模拟地下水流场和GWT模型模拟溶质运移这两个模型各自拥有独立的网格、时间步和边界条件设置模型之间通过GWFIGWT Exchange对象完成水流与溶质之间的耦合。这个架构的最大好处是你可以只跑水流模型不加溶质模块也可以在水流模型收敛稳定之后再叠加溶质模块分析污染物迁移。相比旧版MT3DMS那种“硬耦合”的方式这种松耦合设计让计算更加灵活而且每个模块的输入文件都对应一个独立的字典结构用Python操作时非常符合直觉。我在写源码时习惯先把Simulation、GWF、GWT的字典模板写好这样既能直观地看到每个模块的配置项也能在后期灵活调整参数。2.3 环境搭建与版本匹配避坑写Python调用MODFLOW 6第一步是确保环境干净。我建议用conda或者venv单独建一个虚拟环境不要直接装到系统Python里因为flopy和numpy之间的版本迭代偶尔有不兼容的时候。以我常用的环境为例conda create -n mf6 python3.10 conda activate mf6 pip install flopy numpy pandas matplotlib jupyter pip install modflow6这里有个特别容易踩的坑flopy只是一个“驱动库”它本身不包含MODFLOW 6的可执行文件。你得额外安装modflow6这个包或者在系统里配置好MODFLOW 6的bin路径。flopy运行模型时会自动去寻找可执行文件如果找不到就会报ExecutableNotFoundError。我的建议是在代码开头显式指定路径import flopy mf6_exe your/path/to/mf6 sim flopy.mf6.MFSimulation(sim_namedemo, versionmf6, exe_namemf6_exe)版本匹配上建议flopy版本不低于4.5MODFLOW 6可执行版本不低于6.4.0。我遇到过因为flopy版本太老导致Exchange关键字无法识别的问题升级后一切正常。如果你只是做溶质运移模拟对参数维度和物理量纲的理解比工具版本更重要——这一点后文会专门展开。3. 溶质运移模型的物理机制与源码映射3.1 对流、弥散与吸附——三个必须吃透的过程溶质运移模拟的核心是求解一个描述“污染物在地下水中如何移动”的对流弥散方程ADE。这个方程里三个最重要的物理过程分别是对流污染物随着地下水的流动而整体迁移解决方案很简单——看水流速度场。水动力弥散由于孔隙介质的不均匀性和分子扩散污染物在流动过程中会不断“摊开”浓度前锋不会无限陡峭。弥散度参数aL、aTH等就控制摊开的强度。吸附与延迟部分污染物会被含水层介质吸附导致实际的运移速度低于地下水速度。延迟因子计算公式为 R 1 (ρb · Kd) / θ其中ρb是干容重Kd是分配系数θ是孔隙度。MODFLOW 6的GWT模块在源码层面对这三类过程都提供了显式的参数配置项。你在编写Python源码时本质上就是在给这些物理参数赋值。很多初学者把网格文件和参数文件混在一起其实正确的抽象方式是把“网格形状”和“物理属性”分开定义这样想改某个区域的渗透系数或吸附系数只需要修改对应的数组部分。3.2 溶质模块参数与Python字典的映射关系在flopy中GWT模型的参数通过字典传入。下面这段代码完整定义了一个承压含水层的水流模块和溶质模块的共用参数model_nam_file gwf_demo.nam gwf flopy.mf6.ModflowGwf(sim, modelnamegwf, model_nam_filemodel_nar_file) gwf.dis flopy.mf6.ModflowGwfdis( gwf, nlay1, nrow40, ncol40, delr5.0, delc5.0, top20.0, botm0.0, ) gwf.ic flopy.mf6.ModflowGwfic(gwf, strt10.0) gwf.npf flopy.mf6.ModflowGwfnpf(gwf, icelltype0, k0.5) gwf.chd flopy.mf6.ModflowGwfchd(gwf, stress_period_datachd_data) gwt flopy.mf6.ModflowGwt(sim, modelnamegwt, model_nam_filegwt_demo.nam) gwt.dis flopy.mf6.ModflowGwfdis( gwt, nlay1, nrow40, ncol40, delr5.0, delc5.0, top20.0, botm0.0, )这里有两处需要重点理解第一GWT模型虽然独立但它的网格必须与GWF模型完全一致不能出现网格错位第二GWT模型的储量、初始浓度、边界条件都是独立的配置不能直接沿用GWF的参数。比如你想让溶质只在某个区域有初始浓度你需要单独设置gwt.ic.strt数组而不是用GWF的水头初始值。3.3 溶质边界条件与质量源项的Python实现溶质边界条件比水流边界复杂得多因为除了指定浓度你还得考虑质量通量。MODFLOW 6的GWT模块提供了Flux边界和Concentration边界两种常用类型分别对应“给定质量通量”和“给定浓度值”两种场景。在源码中分别对应flopy.mf6.ModflowGwtflx和flopy.mf6.ModflowGwtcnc。以污染物持续注入为例我会在固定网格单元上设置源项。假设污染源位于第20行第10列持续注入浓度为100 mg/L的溶液代码如下source_concentration 100.0 source_cell [(19, 9)] # 0-indexed gwt.flx flopy.mf6.ModflowGwtflx( gwt, maxbound1, stress_period_data[[(19, 9), 0.5, source_concentration]], )这里第三个字段0.5代表注入流率单位与模型保持一致本例是m³/day第四个字段是质量浓度。注意源项的流率如果设置得比含水层的天然补给量大很多模型很容易出现数值震荡因此在实际工程中需要先算一下注入量占区域总通量的比例。4. 实操从零构建一个完整的溶质运移模型4.1 问题设定与参数选取为了让整个过程更有参考性我们设计一个虚拟但贴近实际的场景一个长200米、宽200米的均质承压含水层网格剖分为40×40每个网格5米见方。地下水从西侧向东侧流动西侧边界水头10米、东侧边界水头8米北侧和南侧为隔水边界。含水层渗透系数取0.5 m/d有效孔隙度0.25。污染源位于含水层中西部持续注入一种保守性污染物不吸附、不降解初始背景浓度为零模拟时长120天用以观察污染羽的扩散范围和形态。单位必须统一。我这里用的长度是米、时间是天因此渗透系数单位就是m/d注入流率单位是m³/d。物质浓度单位可以任意指定只要同一套模型内保持一致即可。物理量纲混乱是我见过最多的建模失误来源建模前先写一张单位清单能省下后期大量的排错时间。4.2 完整可运行的Python源码下面这段是按上述参数整理出的完整源码为了便于阅读和理解我做了精简去掉了后处理绘图部分但保留了模型构建和运行的全流程import flopy import numpy as np mf6_exe /usr/bin/mf6 sim_name demo_transport ws ./mf6_demo sim flopy.mf6.MFSimulation(sim_namesim_name, versionmf6, exe_namemf6_exe, sim_wsws) # 时间步设置 sim.tdis flopy.mf6.ModflowTdis(sim, nper1, perioddata[(120.0, 10, 1.0)]) # GWF 水流模型 gwf flopy.mf6.ModflowGwf(sim, modelnamegwf, model_nam_filegwf.nam) gwf.dis flopy.mf6.ModflowGwfdis( gwf, nlay1, nrow40, ncol40, delr5.0, delc5.0, top20.0, botm0.0 ) gwf.ic flopy.mf6.ModflowGwfic(gwf, strt10.0) gwf.npf flopy.mf6.ModflowGwfnpf(gwf, icelltype0, k0.5) gwf.chd flopy.mf6.ModflowGwfchd( gwf, stress_period_data[[(0, 19, 0), 10.0], [(0, 19, 39), 8.0]] ) # GWT 溶质运移模型 gwt flopy.mf6.ModflowGwt(sim, modelnamegwt, model_nam_filegwt.nam) gwt.dis flopy.mf6.ModflowGwfdis( gwt, nlay1, nrow40, ncol40, delr5.0, delc5.0, top20.0, botm0.0 ) gwt.ic flopy.mf6.ModflowGwtic(gwt, strt0.0) # 弥散模型参数 disp flopy.mf6.ModflowGwtmst(gwt, porosity0.25) disp.disp flopy.mf6.ModflowGwtdsp( gwt, alh5.0, ath10.5, ath20.5 ) # 污染源注入 gwt.flx flopy.mf6.ModflowGwtflx( gwt, maxbound1, stress_period_data[[(0, 19, 10), 5.0, 100.0]], ) # GWF-GWT 耦合 gwfgwt flopy.mf6.ModflowGwfgwt( sim, exgtypeGWF6-GWT6, exgmnameagwf, exgmnamebgwt, extrudedTrue, ) # 输出控制保存浓度场 gwt.oc flopy.mf6.ModflowGwtoc( gwt, budget_filerecordgwt.cbc, concentration_filerecordgwt.ucn, saverecord{( 120.0): [ HEAD, CONCENTRATION, ]}, ) # 写入并运行 sim.write_simulation() sim.run_simulation()这段代码的核心逻辑并不复杂先用MFSimulation搭好模拟容器再放入水流模型和溶质模型最后定义交换关系。整个流程对应MODFLOW 6的输入文件组织方式写文件时flopy会自动展开为.nam、.dis、.npf、.dsp等文本文件。如果你手动看过这些文件会发现里面就是标准的关键字格式这也方便你直接阅读和修改底层输入。4.3 结果读取与可视化模型跑完以后结果保存在工作目录下的gwt.ucn文件中。这是一个二进制文件不能直接用文本编辑器打开需要用flopy的HeadFile对象读取。我的习惯是把指定时间步的浓度数组读取成numpy数组再配合matplotlib画浓度等值线或者热力图。from flopy.utils import HeadFile concentration_file HeadFile(os.path.join(ws, gwt.ucn), textCONCENTRATION) conc_data concentration_file.get_data(totim120.0) import matplotlib.pyplot as plt plt.imshow(conc_data[0, :, :], originlower, cmapRdYlBu_r) plt.colorbar(labelConcentration (mg/L)) plt.xlabel(Column) plt.ylabel(Row) plt.show()从成像结果可以清楚看到污染物在120天内的运移前缘如果设置的对流速度大于弥散速率污染羽会呈明显的顺水流拖尾形状。这也是检验模型物理合理性的一种直观手段。5. 模块耦合、时间步长与数值稳定性源码背后的“为什么”5.1 GWF与GWT的耦合方式及适用场景MODFLOW 6的GWT模块默认读取GWF模型计算得到的速度场然后把速度场带入溶质运移方程求解浓度分布。GWF-GWT耦合方式分为“并行耦合”和“顺序耦合”两种。并行耦合需要两个模型同时运行每经过一个时间步就交换数据顺序耦合则是先跑完GWF再把速度场传给GWT。我在日常项目中优先选择并行耦合因为这样最容易保持物理场的同步性。如果水流变化速度很快比如突然大量抽水顺序耦合可能会因为速度场的滞后导致浓度结果偏差。ModflowGwfgwt这个工具在flopy初始化时默认就在做并行耦合只要在extrudedTrue这一行打开。这个extruded参数的含义是“是否允许从GWF模型挤出速度到GWT模型”通常必须打开否则你看到的浓度场可能永远不变。5.2 时间步长选择与Courant数的工程经验溶质运移模型对时间步长比水流模型更敏感因为污染物是“跟着水流走”的如果时间步长太大物质在一个时间步内穿越了多个网格数值散就会变得非常严重。实际工作中可以用Courant数来诊断C v · Δt / Δx其中v是孔隙流速Δt是时间步长Δx是网格尺寸。理论上C ≤ 1 才能保证稳定性但我个人在实际项目中会严格控制C ≤ 0.5尤其在有污染源注入的区域预留更大的安全裕量才能避免浓度震荡。拿上面的案例来说水力梯度约 (10 - 8) / 200 0.01 m/m渗透系数0.5 m/d孔隙度0.25则达西流速为0.005 m/d孔隙流速为0.02 m/d。在40个网格、每个网格5米的情况下Δt如果设成10天则C 0.02 × 10 / 5 0.04远小于1所以计算非常稳定。如果网格加密到1米同样的时间步长C 就变成0.2依然可以接受但如果你把时间步长加到100天C 2结果就会明显恶化。5.3 弥散度参数的敏感性与选取逻辑弥散度是溶质运移模型里最敏感、也最不容易确定的参数。纵向弥散度aL通常取到网格尺寸的1/10到1倍横向弥散度一般取aL的0.1倍。上面源码中设置aL5米是因为网格尺寸是5米取了一个中值。实际项目中如果有示踪试验数据优先用试验结果来拟合弥散度没有现场数据的项目建议至少做一组敏感性分析看看结果对弥散度变化的响应程度这样才能判断模型结论的稳健性。需要注意的是数值弥散和物理弥散经常混在一起是不可分离的。当你发现污染羽的扩散范围远超预期先不要急着调大弥散度可以先加密网格、减小时间步长看看结果是否发生明显变化。如果加密后结果显著改变说明你的模型还在“数值弥散主导”阶段此时讨论物理弥散参数没有意义。6. 常见报错、坑位排查与性能优化实录6.1 让我印象深刻的三个典型报错我在迭代这个模型时先后踩过三个典型的坑每个都值得单独记录。报错一MODFLOW执行文件未找到这是我第一次运行脚本就遇到的。flopy在调用sim.run_simulation()时会尝试在系统PATH中搜索mf6但因为我用的是conda环境没有将MODFLOW安装到PATH。解决办法很简单在创建MFSimulation时显式指定exe_name为可执行文件的绝对路径。这件事最容易被忽略因为错误信息有时会提示FileNotFoundError有时却直接静默闪退如果你发现模型没有任何输出文件优先检查这一点。报错二Gwell Water Flow模型与溶质模型网格不一致我在一次快速测试中给GWF模型用40×40网格、GWT模型用20×20网格模型初建时flopy并没有报错但运行时报了SpatialReferenceError。这其实是物理上的要求溶质运移必须建立在同一套网格上。除非你有特别复杂的参数插值需求否则我建议直接让两个模型共用一样的delr、delc、top、botm数组省事也安全。报错三浓度出现负值负浓度在溶质运移数值模拟中是很常见的问题尤其在对流占主导时拿中心差分格式求解负浓度几乎一定会出现。MODFLOW 6内置的上风格式能大幅减少负值出现但如果你的模型中污染源浓度梯度太大仍然可能出现负值。遇到这种情况首选方案是加密网格或缩小时间步长而不是盲目增大物理弥散度来“磨平”浓度场否则得到的结果没有物理意义。6.2 模型运行效率的优化实践MODFLOW 6本身用Fortran写的求解效率其实很高。我发现性能瓶颈更多出在Python端。如果你在一个大模型上反复读二进制结果文件做后处理读取时间可能比求解时间还长。我的经验是不要反复实例化HeadFile对象尽量一次性读取全部时间步的浓度数据到内存再统一处理。后处理时优先用numpy数组操作不要在循环里逐个网格取值。如果模拟时长很长建议只保存关键时间点的浓度场而不是每个中间步都存储。MODFLOW 6的oc模块提供了精确的时间点控制把saverecord里的列表设置成你真正需要的时刻即可能大幅压缩输出文件体积。6.3 模型结果合理性检验清单模型能跑通只是第一步结果是否可信还得过逻辑检验。我给自己定了一个最小检验清单浓度范围是否落在0到源浓度之间如果出现负值或超过源浓度的超调量超过5%说明数值扩散控制不够好。质量是否守恒以本案例为例污染物总注入质量 注入流率 × 注入浓度 × 时长 5.0 × 100 × 120 60000质量单位读取结果文件中所有网格的浓度乘以网格体积乘以孔隙度累加后应等于或者非常接近这个值。误差超过10%就说明模型设置有问题。污染羽前锋是否与流速场方向一致如果污染羽明显逆着水流扩散往往是边界条件或者弥散参数设置有问题。稳态结果是否合理如果有持续注入源随着时间推进污染羽应该趋向一个稳定形状而不是无限扩张。这四点检查可以说是溶质运移数值模拟的“安检门”建议每次跑完模型都系统性过一遍。7. 溶质运移Python建模的进阶方向与个人心得做了一段时间的MODFLOW 6溶质运移模拟后我越来越觉得这套工具链最大的价值不在于单个模型跑得有多快而在于把“建模—模拟—分析—决策”整条链路变成可编程、可复现、可共享的过程。过去花在手动整理数据上的时间现在完全可以集中到参数校准和方案比较上。还想提醒一点Python源码本身好不好不只看结果对不对还要看扩展和维护性。能写成函数的地方不要平铺直叙能用字典管理参数的地方不要散落常量。即便只是一个小项目今天我仍然会为“污染源位置、注入浓度、注入时长”这类变量单独建一个config模块方便批量生成情景。对于刚上手的朋友我建议先按本文的框架把例子跑通然后尝试修改边界条件类型、增加一层含水层、加入吸附/降解参数一步步把模型复杂度提上去。溶质运移模拟没有捷径但它也远没有想象中那么高不可攀——关键是把物理过程吃透再对照源码逐步理解每个参数在方程里扮演的角色。只要这两个环节没有断档你手里这套Python建模流程就能成为解决实际污染迁移问题的可靠工具。 Python+MODFLOW 6：溶质运移建模从入门到实战
1. 从核心需求出发为什么要在Python里跑Modflow6地球物理与水文地质圈这两年有个明显趋势大家不再满足于把MODFLOW当“黑盒”点两下界面出张等值线图而是想把建模流程纳入更完整的科学计算链路里。MODFLOW 6这个版本的最大变化是整个代码框架用面向对象思路重写底层的求解器、时间步进、边界条件都变成了一套可以“勾搭”进Python生态的模块。配合flopy工具库你几乎可以在Jupyter Notebook里完成从网格生成、参数赋值、运行模拟到读取结果的完整闭环。我在实际项目里最强烈的体感是过去用传统界面建模调一次参数要反复点菜单、导出、再导入现在只需要写好一份可复现的Python脚本改参数等于改变量批量跑情景模拟变得极其顺手。尤其是溶质运移这块MODFLOW 6把地下水流模块GWF和溶质运移模块GWT分开建模两者通过 exchange 机制耦合逻辑比旧版本清晰太多。这篇内容就是围绕一个典型的含承压含水层的溶质运移案例从建模思路、Python源码结构、参数计算到调试排错完整走一遍。不管你是刚接触数值模拟的研究生还是工作中需要用模型辅助决策的工程师按照这个流程都能把模型跑起来并且能看懂每一行代码在干什么。内容中涉及的所有脚本结构和参数设计均来自我在项目中的实际使用方案亲测可跑、可复现。2. MODFLOW 6的Python建模体系全景2.1 从传统建模到Python API的范式转变传统的MODFLOW建模流程大致是划分网格→分配水文地质参数→设置边界条件→运行计算→后处理。这套流程在GMS或ModelMuse里做当然没有问题但它最大的隐患是“过程不可追踪”。模型文件是谁生成的、参数是哪个版本改的、边界条件具体怎么设的时间一长很容易说不清。换成Python脚本之后整个过程变成了一条可版本管理的“代码生产线”。flopy作为USGS官方维护的Python工具包天然支持MODFLOW 6全部关键字你可以直接在脚本里new一个Simulation对象往里面塞GWF模型、GWT模型、离散化信息、应力期和输出控制。整份脚本本身就是模型说明书以后想回溯哪个版本改了什么都清楚这对我这种经常要跟多个项目并行的人来讲价值太大了。2.2 flopy与MODFLOW 6的模块化架构要写对源码先搞清楚MODFLOW 6的模块划分。顶层是Simulation负责管理整个模拟过程Simulation下面挂GWF模型模拟地下水流场和GWT模型模拟溶质运移这两个模型各自拥有独立的网格、时间步和边界条件设置模型之间通过GWFIGWT Exchange对象完成水流与溶质之间的耦合。这个架构的最大好处是你可以只跑水流模型不加溶质模块也可以在水流模型收敛稳定之后再叠加溶质模块分析污染物迁移。相比旧版MT3DMS那种“硬耦合”的方式这种松耦合设计让计算更加灵活而且每个模块的输入文件都对应一个独立的字典结构用Python操作时非常符合直觉。我在写源码时习惯先把Simulation、GWF、GWT的字典模板写好这样既能直观地看到每个模块的配置项也能在后期灵活调整参数。2.3 环境搭建与版本匹配避坑写Python调用MODFLOW 6第一步是确保环境干净。我建议用conda或者venv单独建一个虚拟环境不要直接装到系统Python里因为flopy和numpy之间的版本迭代偶尔有不兼容的时候。以我常用的环境为例conda create -n mf6 python3.10 conda activate mf6 pip install flopy numpy pandas matplotlib jupyter pip install modflow6这里有个特别容易踩的坑flopy只是一个“驱动库”它本身不包含MODFLOW 6的可执行文件。你得额外安装modflow6这个包或者在系统里配置好MODFLOW 6的bin路径。flopy运行模型时会自动去寻找可执行文件如果找不到就会报ExecutableNotFoundError。我的建议是在代码开头显式指定路径import flopy mf6_exe your/path/to/mf6 sim flopy.mf6.MFSimulation(sim_namedemo, versionmf6, exe_namemf6_exe)版本匹配上建议flopy版本不低于4.5MODFLOW 6可执行版本不低于6.4.0。我遇到过因为flopy版本太老导致Exchange关键字无法识别的问题升级后一切正常。如果你只是做溶质运移模拟对参数维度和物理量纲的理解比工具版本更重要——这一点后文会专门展开。3. 溶质运移模型的物理机制与源码映射3.1 对流、弥散与吸附——三个必须吃透的过程溶质运移模拟的核心是求解一个描述“污染物在地下水中如何移动”的对流弥散方程ADE。这个方程里三个最重要的物理过程分别是对流污染物随着地下水的流动而整体迁移解决方案很简单——看水流速度场。水动力弥散由于孔隙介质的不均匀性和分子扩散污染物在流动过程中会不断“摊开”浓度前锋不会无限陡峭。弥散度参数aL、aTH等就控制摊开的强度。吸附与延迟部分污染物会被含水层介质吸附导致实际的运移速度低于地下水速度。延迟因子计算公式为 R 1 (ρb · Kd) / θ其中ρb是干容重Kd是分配系数θ是孔隙度。MODFLOW 6的GWT模块在源码层面对这三类过程都提供了显式的参数配置项。你在编写Python源码时本质上就是在给这些物理参数赋值。很多初学者把网格文件和参数文件混在一起其实正确的抽象方式是把“网格形状”和“物理属性”分开定义这样想改某个区域的渗透系数或吸附系数只需要修改对应的数组部分。3.2 溶质模块参数与Python字典的映射关系在flopy中GWT模型的参数通过字典传入。下面这段代码完整定义了一个承压含水层的水流模块和溶质模块的共用参数model_nam_file gwf_demo.nam gwf flopy.mf6.ModflowGwf(sim, modelnamegwf, model_nam_filemodel_nar_file) gwf.dis flopy.mf6.ModflowGwfdis( gwf, nlay1, nrow40, ncol40, delr5.0, delc5.0, top20.0, botm0.0, ) gwf.ic flopy.mf6.ModflowGwfic(gwf, strt10.0) gwf.npf flopy.mf6.ModflowGwfnpf(gwf, icelltype0, k0.5) gwf.chd flopy.mf6.ModflowGwfchd(gwf, stress_period_datachd_data) gwt flopy.mf6.ModflowGwt(sim, modelnamegwt, model_nam_filegwt_demo.nam) gwt.dis flopy.mf6.ModflowGwfdis( gwt, nlay1, nrow40, ncol40, delr5.0, delc5.0, top20.0, botm0.0, )这里有两处需要重点理解第一GWT模型虽然独立但它的网格必须与GWF模型完全一致不能出现网格错位第二GWT模型的储量、初始浓度、边界条件都是独立的配置不能直接沿用GWF的参数。比如你想让溶质只在某个区域有初始浓度你需要单独设置gwt.ic.strt数组而不是用GWF的水头初始值。3.3 溶质边界条件与质量源项的Python实现溶质边界条件比水流边界复杂得多因为除了指定浓度你还得考虑质量通量。MODFLOW 6的GWT模块提供了Flux边界和Concentration边界两种常用类型分别对应“给定质量通量”和“给定浓度值”两种场景。在源码中分别对应flopy.mf6.ModflowGwtflx和flopy.mf6.ModflowGwtcnc。以污染物持续注入为例我会在固定网格单元上设置源项。假设污染源位于第20行第10列持续注入浓度为100 mg/L的溶液代码如下source_concentration 100.0 source_cell [(19, 9)] # 0-indexed gwt.flx flopy.mf6.ModflowGwtflx( gwt, maxbound1, stress_period_data[[(19, 9), 0.5, source_concentration]], )这里第三个字段0.5代表注入流率单位与模型保持一致本例是m³/day第四个字段是质量浓度。注意源项的流率如果设置得比含水层的天然补给量大很多模型很容易出现数值震荡因此在实际工程中需要先算一下注入量占区域总通量的比例。4. 实操从零构建一个完整的溶质运移模型4.1 问题设定与参数选取为了让整个过程更有参考性我们设计一个虚拟但贴近实际的场景一个长200米、宽200米的均质承压含水层网格剖分为40×40每个网格5米见方。地下水从西侧向东侧流动西侧边界水头10米、东侧边界水头8米北侧和南侧为隔水边界。含水层渗透系数取0.5 m/d有效孔隙度0.25。污染源位于含水层中西部持续注入一种保守性污染物不吸附、不降解初始背景浓度为零模拟时长120天用以观察污染羽的扩散范围和形态。单位必须统一。我这里用的长度是米、时间是天因此渗透系数单位就是m/d注入流率单位是m³/d。物质浓度单位可以任意指定只要同一套模型内保持一致即可。物理量纲混乱是我见过最多的建模失误来源建模前先写一张单位清单能省下后期大量的排错时间。4.2 完整可运行的Python源码下面这段是按上述参数整理出的完整源码为了便于阅读和理解我做了精简去掉了后处理绘图部分但保留了模型构建和运行的全流程import flopy import numpy as np mf6_exe /usr/bin/mf6 sim_name demo_transport ws ./mf6_demo sim flopy.mf6.MFSimulation(sim_namesim_name, versionmf6, exe_namemf6_exe, sim_wsws) # 时间步设置 sim.tdis flopy.mf6.ModflowTdis(sim, nper1, perioddata[(120.0, 10, 1.0)]) # GWF 水流模型 gwf flopy.mf6.ModflowGwf(sim, modelnamegwf, model_nam_filegwf.nam) gwf.dis flopy.mf6.ModflowGwfdis( gwf, nlay1, nrow40, ncol40, delr5.0, delc5.0, top20.0, botm0.0 ) gwf.ic flopy.mf6.ModflowGwfic(gwf, strt10.0) gwf.npf flopy.mf6.ModflowGwfnpf(gwf, icelltype0, k0.5) gwf.chd flopy.mf6.ModflowGwfchd( gwf, stress_period_data[[(0, 19, 0), 10.0], [(0, 19, 39), 8.0]] ) # GWT 溶质运移模型 gwt flopy.mf6.ModflowGwt(sim, modelnamegwt, model_nam_filegwt.nam) gwt.dis flopy.mf6.ModflowGwfdis( gwt, nlay1, nrow40, ncol40, delr5.0, delc5.0, top20.0, botm0.0 ) gwt.ic flopy.mf6.ModflowGwtic(gwt, strt0.0) # 弥散模型参数 disp flopy.mf6.ModflowGwtmst(gwt, porosity0.25) disp.disp flopy.mf6.ModflowGwtdsp( gwt, alh5.0, ath10.5, ath20.5 ) # 污染源注入 gwt.flx flopy.mf6.ModflowGwtflx( gwt, maxbound1, stress_period_data[[(0, 19, 10), 5.0, 100.0]], ) # GWF-GWT 耦合 gwfgwt flopy.mf6.ModflowGwfgwt( sim, exgtypeGWF6-GWT6, exgmnameagwf, exgmnamebgwt, extrudedTrue, ) # 输出控制保存浓度场 gwt.oc flopy.mf6.ModflowGwtoc( gwt, budget_filerecordgwt.cbc, concentration_filerecordgwt.ucn, saverecord{( 120.0): [ HEAD, CONCENTRATION, ]}, ) # 写入并运行 sim.write_simulation() sim.run_simulation()这段代码的核心逻辑并不复杂先用MFSimulation搭好模拟容器再放入水流模型和溶质模型最后定义交换关系。整个流程对应MODFLOW 6的输入文件组织方式写文件时flopy会自动展开为.nam、.dis、.npf、.dsp等文本文件。如果你手动看过这些文件会发现里面就是标准的关键字格式这也方便你直接阅读和修改底层输入。4.3 结果读取与可视化模型跑完以后结果保存在工作目录下的gwt.ucn文件中。这是一个二进制文件不能直接用文本编辑器打开需要用flopy的HeadFile对象读取。我的习惯是把指定时间步的浓度数组读取成numpy数组再配合matplotlib画浓度等值线或者热力图。from flopy.utils import HeadFile concentration_file HeadFile(os.path.join(ws, gwt.ucn), textCONCENTRATION) conc_data concentration_file.get_data(totim120.0) import matplotlib.pyplot as plt plt.imshow(conc_data[0, :, :], originlower, cmapRdYlBu_r) plt.colorbar(labelConcentration (mg/L)) plt.xlabel(Column) plt.ylabel(Row) plt.show()从成像结果可以清楚看到污染物在120天内的运移前缘如果设置的对流速度大于弥散速率污染羽会呈明显的顺水流拖尾形状。这也是检验模型物理合理性的一种直观手段。5. 模块耦合、时间步长与数值稳定性源码背后的“为什么”5.1 GWF与GWT的耦合方式及适用场景MODFLOW 6的GWT模块默认读取GWF模型计算得到的速度场然后把速度场带入溶质运移方程求解浓度分布。GWF-GWT耦合方式分为“并行耦合”和“顺序耦合”两种。并行耦合需要两个模型同时运行每经过一个时间步就交换数据顺序耦合则是先跑完GWF再把速度场传给GWT。我在日常项目中优先选择并行耦合因为这样最容易保持物理场的同步性。如果水流变化速度很快比如突然大量抽水顺序耦合可能会因为速度场的滞后导致浓度结果偏差。ModflowGwfgwt这个工具在flopy初始化时默认就在做并行耦合只要在extrudedTrue这一行打开。这个extruded参数的含义是“是否允许从GWF模型挤出速度到GWT模型”通常必须打开否则你看到的浓度场可能永远不变。5.2 时间步长选择与Courant数的工程经验溶质运移模型对时间步长比水流模型更敏感因为污染物是“跟着水流走”的如果时间步长太大物质在一个时间步内穿越了多个网格数值散就会变得非常严重。实际工作中可以用Courant数来诊断C v · Δt / Δx其中v是孔隙流速Δt是时间步长Δx是网格尺寸。理论上C ≤ 1 才能保证稳定性但我个人在实际项目中会严格控制C ≤ 0.5尤其在有污染源注入的区域预留更大的安全裕量才能避免浓度震荡。拿上面的案例来说水力梯度约 (10 - 8) / 200 0.01 m/m渗透系数0.5 m/d孔隙度0.25则达西流速为0.005 m/d孔隙流速为0.02 m/d。在40个网格、每个网格5米的情况下Δt如果设成10天则C 0.02 × 10 / 5 0.04远小于1所以计算非常稳定。如果网格加密到1米同样的时间步长C 就变成0.2依然可以接受但如果你把时间步长加到100天C 2结果就会明显恶化。5.3 弥散度参数的敏感性与选取逻辑弥散度是溶质运移模型里最敏感、也最不容易确定的参数。纵向弥散度aL通常取到网格尺寸的1/10到1倍横向弥散度一般取aL的0.1倍。上面源码中设置aL5米是因为网格尺寸是5米取了一个中值。实际项目中如果有示踪试验数据优先用试验结果来拟合弥散度没有现场数据的项目建议至少做一组敏感性分析看看结果对弥散度变化的响应程度这样才能判断模型结论的稳健性。需要注意的是数值弥散和物理弥散经常混在一起是不可分离的。当你发现污染羽的扩散范围远超预期先不要急着调大弥散度可以先加密网格、减小时间步长看看结果是否发生明显变化。如果加密后结果显著改变说明你的模型还在“数值弥散主导”阶段此时讨论物理弥散参数没有意义。6. 常见报错、坑位排查与性能优化实录6.1 让我印象深刻的三个典型报错我在迭代这个模型时先后踩过三个典型的坑每个都值得单独记录。报错一MODFLOW执行文件未找到这是我第一次运行脚本就遇到的。flopy在调用sim.run_simulation()时会尝试在系统PATH中搜索mf6但因为我用的是conda环境没有将MODFLOW安装到PATH。解决办法很简单在创建MFSimulation时显式指定exe_name为可执行文件的绝对路径。这件事最容易被忽略因为错误信息有时会提示FileNotFoundError有时却直接静默闪退如果你发现模型没有任何输出文件优先检查这一点。报错二Gwell Water Flow模型与溶质模型网格不一致我在一次快速测试中给GWF模型用40×40网格、GWT模型用20×20网格模型初建时flopy并没有报错但运行时报了SpatialReferenceError。这其实是物理上的要求溶质运移必须建立在同一套网格上。除非你有特别复杂的参数插值需求否则我建议直接让两个模型共用一样的delr、delc、top、botm数组省事也安全。报错三浓度出现负值负浓度在溶质运移数值模拟中是很常见的问题尤其在对流占主导时拿中心差分格式求解负浓度几乎一定会出现。MODFLOW 6内置的上风格式能大幅减少负值出现但如果你的模型中污染源浓度梯度太大仍然可能出现负值。遇到这种情况首选方案是加密网格或缩小时间步长而不是盲目增大物理弥散度来“磨平”浓度场否则得到的结果没有物理意义。6.2 模型运行效率的优化实践MODFLOW 6本身用Fortran写的求解效率其实很高。我发现性能瓶颈更多出在Python端。如果你在一个大模型上反复读二进制结果文件做后处理读取时间可能比求解时间还长。我的经验是不要反复实例化HeadFile对象尽量一次性读取全部时间步的浓度数据到内存再统一处理。后处理时优先用numpy数组操作不要在循环里逐个网格取值。如果模拟时长很长建议只保存关键时间点的浓度场而不是每个中间步都存储。MODFLOW 6的oc模块提供了精确的时间点控制把saverecord里的列表设置成你真正需要的时刻即可能大幅压缩输出文件体积。6.3 模型结果合理性检验清单模型能跑通只是第一步结果是否可信还得过逻辑检验。我给自己定了一个最小检验清单浓度范围是否落在0到源浓度之间如果出现负值或超过源浓度的超调量超过5%说明数值扩散控制不够好。质量是否守恒以本案例为例污染物总注入质量 注入流率 × 注入浓度 × 时长 5.0 × 100 × 120 60000质量单位读取结果文件中所有网格的浓度乘以网格体积乘以孔隙度累加后应等于或者非常接近这个值。误差超过10%就说明模型设置有问题。污染羽前锋是否与流速场方向一致如果污染羽明显逆着水流扩散往往是边界条件或者弥散参数设置有问题。稳态结果是否合理如果有持续注入源随着时间推进污染羽应该趋向一个稳定形状而不是无限扩张。这四点检查可以说是溶质运移数值模拟的“安检门”建议每次跑完模型都系统性过一遍。7. 溶质运移Python建模的进阶方向与个人心得做了一段时间的MODFLOW 6溶质运移模拟后我越来越觉得这套工具链最大的价值不在于单个模型跑得有多快而在于把“建模—模拟—分析—决策”整条链路变成可编程、可复现、可共享的过程。过去花在手动整理数据上的时间现在完全可以集中到参数校准和方案比较上。还想提醒一点Python源码本身好不好不只看结果对不对还要看扩展和维护性。能写成函数的地方不要平铺直叙能用字典管理参数的地方不要散落常量。即便只是一个小项目今天我仍然会为“污染源位置、注入浓度、注入时长”这类变量单独建一个config模块方便批量生成情景。对于刚上手的朋友我建议先按本文的框架把例子跑通然后尝试修改边界条件类型、增加一层含水层、加入吸附/降解参数一步步把模型复杂度提上去。溶质运移模拟没有捷径但它也远没有想象中那么高不可攀——关键是把物理过程吃透再对照源码逐步理解每个参数在方程里扮演的角色。只要这两个环节没有断档你手里这套Python建模流程就能成为解决实际污染迁移问题的可靠工具。