MT3D-USGS
摘要：MT3D-USGS: 使用 SFR/LAK/UZF 包进行传输 (SFT/LKT/UZT)，以及化学反应 (RCT) 本教程更全面地演示了如何设置一个使用 MT3D-USGS 首次发布中包含的所有新包的 MT3D-USGS 模型，并包括化学反应 (RCT)。 问题描述: 模型大小: 300 行 x 300 列 x 3 层 x 2 应力期模型。 水流模型: 使用 SFR（地表水流路径）、LAK（湖泊,MT3D-USGS
MT3D-USGS: 使用 SFR/LAK/UZF 包进行传输 (SFT/LKT/UZT)，以及化学反应 (RCT)
本教程更全面地演示了 怎样设置一个使用 MT3D-USGS 首次发布中包含的所有新包的 MT3D-USGS 模型，并包括化学反应 (RCT)。
难题描述:
模型 大致: 300 行 x 300 列 x 3 层 x 2 应力期模型。 水流模型: 使用 SFR（地表水流路径）、LAK（湖泊）和 UZF（非饱和区流）包， 并且三者之间存在连接。 溶质传输模型:
模拟河流流量传输 (SFT)，并连接到一个单一的湖泊 (LKT)。 模拟地表径流和泉水排放 (UZT) 到地表水网络。
化学反应: 包含化学反应 (RCT)。
首先导入一些库:
Python
import os import sys from pathlib import Path from tempfile import TemporaryDirectory import git import tplotlib as mpl import numpy as np import pandas as pd import pooch import flopy print(sys.version) print(f"numpy version: {np.__version__}") print(f" tplotlib version: {mpl.__version__}") print(f"flopy version: {flopy.__version__}")输出:
3.12.10 | packaged by conda-forge | ( in, Apr 10 2025, 22:21:13) [GCC 13.3.0] numpy version: 2.2.5 tplotlib version: 3.10.3 flopy version: 3.9.3创建 MODFLOW 模型并将其存储在变量 'mf' 中。模型名称将作为所有 MODFLOW 文件的名称。exe_name 应该是 MODFLOW 可执行文件的名称。在这种情况下，我们希望使用版本：'mfnwt'，即 MODFLOW-NWT。
Python
# 临时目录 temp_dir = TemporaryDirectory() model_ws = temp_dir.name # 模型 职业空间路径 # 检查我们是否在仓库中并定义数据路径。 try: root = Path(git.Repo(".", search_parent_directories=True).working_dir) except: root = None data_path = root / "examples" / "data" if root else Path.cwd() modelpth = os.path.join(model_ws, "no3") # 具体模型文件存放的子目录 modelname = "no3" # 模型名称 mfexe = "mfnwt" # MODFLOW 可执行文件名称 (MODFLOW-NWT) mtexe = "mt3dusgs" # MT3D-USGS 可执行文件名称 # 确保 modelpth 目录存在 if not os.path.isdir(modelpth): os. kedirs(modelpth, exist_ok=True) # 如果目录不存在则创建它 # 在 flopy 中实例化 MODFLOW 对象 mf = flopy.modflow.Modflow( modelname=modelname, exe_name=mfexe, model_ws=modelpth, version="mfnwt" )设置模型离散化
Python
Lx = 90000.0 # X 路线长度 (米) Ly = 90000.0 # Y 路线长度 (米) nrow = 300 # 行数 ncol = 300 # 列数 nlay = 3 # 层数 delr = Lx / ncol # 列宽 (米) delc = Ly / nrow # 行宽 (米) x x = ncol * delr # X 路线最大坐标 y x = nrow * delc # Y 路线最大坐标 X, Y = np.meshgrid( # 创建网格坐标数组 np.linspace(delr / 2, x x - delr / 2, ncol), np.linspace(y x - delc / 2, 0 + delc / 2, nrow), )实例化 MODFLOW-NWT 的输出控制 (OC) 包
Python
oc = flopy.modflow.ModflowOc(mf) # 创建输出控制包，使用默认设置实例化 MODFLOW-NWT 的求解器包
Python
# Newton-Raphson 求解器：创建 flopy nwt 包对象 headtol = 1.0e-4 # 水头收敛容差 fluxtol = 5 # 流量收敛容差 xiterout = 5000 # 最大外部迭代次数 thickfact = 1e-06 # 厚度乘数，用于调整相邻单元格之间的流量 linmeth = 2 # 线性方程组求解 技巧。2 表示 GMRES iprnwt = 1 # 打印 NWT 模块诊断信息。1 表示在列表文件中打印 ibotav = 1 # 计算层底部高程的平均 技巧。1 表示使用下层底部高程的平均值 nwt = flopy.modflow.ModflowNwt( mf, headtol=headtol, fluxtol=fluxtol, xiterout= xiterout, thickfact=thickfact, linmeth=linmeth, iprnwt=iprnwt, ibotav=ibotav, options="SIMPLE", # NWT 的选项，"SIMPLE" 一个预设选项 )实例化 MODFLOW-NWT 的离散化 (DIS) 包
Python
# 第 1 层顶面高程使用 GW Vistas 确定并本地存储 fname = "grnd_elv.txt" folder_name = "mt3d_example_sft_lkt_uzt" pooch.retrieve( # 使用 pooch 下载文件 url=f"https://github.com/modflowpy/flopy/raw/develop/examples/data/{folder_name}/dis_arrays/{fname}", fname=fname, path=data_path / folder_name / "dis_arrays", known_hash=None, ) elv_pth = data_path / folder_name / "dis_arrays" / fname grndElv = np.loadtxt(elv_pth) # 加载顶面高程数据 # 第 1 层底面高程也通过 GUI 确定并本地存储 fname = "bot1.txt" folder_name = "mt3d_example_sft_lkt_uzt" pooch.retrieve( # 使用 pooch 下载文件 url=f"https://github.com/modflowpy/flopy/raw/develop/examples/data/{folder_name}/dis_arrays/{fname}", fname=fname, path=data_path / folder_name / "dis_arrays", known_hash=None, ) bt1_pth = data_path / folder_name / "dis_arrays" / fname bot1Elv = np.loadtxt(bt1_pth) # 加载第 1 层底面高程数据 bot2Elv = np.ones(bot1Elv.shape) * 100 # 第 2 层底面高程设置为 100 bot3Elv = np.zeros(bot2Elv.shape) # 第 3 层底面高程设置为 0 botm = [bot1Elv, bot2Elv, bot3Elv] # 所有层底面高程列表 botm = np.array(botm) # 转换为 NumPy 数组 Steady = [False, False] # 每个应力期是否为稳定流。两个应力期都为瞬态流 nstp = [1, 1] # 每个应力期的 时刻步数 (这里每个应力期只有一个 时刻步，通常瞬态流不止一个) t ult = [1.0, 1.0] # 时刻步乘数 (这里每个 时刻步长度不变) # 应力期 (Stress periods) perlen = [9131.25, 9131.25] # 每个应力期的长度 (天) # 创建离散化对象 # itmuni = 4 (天); lenuni = 1 (英尺) dis = flopy.modflow.ModflowDis( mf, nlay, # 层数 nrow, # 行数 ncol, # 列数 nper=2, # 应力期数 delr=delr, # 列宽 delc=delc, # 行宽 top=grndElv, # 顶面高程 botm=botm, # 各层底面高程 laycbd=0, # 无隔水层 itmuni=4, # 时刻单位：4 代表天 (days) lenuni=1, # 长度单位：1 代表英尺 (feet) steady=Steady, # 稳定流/瞬态流标志 nstp=nstp, # 每个应力期的 时刻步数 t ult=t ult, # 时刻步乘数 perlen=perlen, # 每个应力期的长度 )实例化 MODFLOW-NWT 的上游加权 (UPW) 流包
Python
# UPW 必须在 DIS 实例化之后实例化。否则，在 mf.write_input() 经过中，flopy 将崩溃。 # UPW 输入文件的第一行是：IUPWCB HDRY NPUPW IPHDRY hdry = -1.00e30 # 干燥单元的水头值 iphdry = 0 # 确定 怎样处理干燥单元的标志 # 接下来的变量是：LAYTYP, LAY G, CHANI, LAYVKA, LAYWET laytyp = [1, 3, 3] # 层类型：>0 表示可变层 (非承压或承压-非承压混合) # laytyp[0]=1: 第一层是可变层（非承压） # laytyp[1]=3: 第二层是第三种类型（可变层） # laytyp[2]=3: 第三层是第三种类型（可变层） layavg = 0 # 0: 调 安宁均 (用于计算层间垂向导水率) chani = 1.0 # >0: CHANI 是整个层的水平各向异性 (水平 Kx/Ky 比率)。1.0 表示各向同性 layvka = 0 # =0: 表示 VKA 是垂向导水率 (VK)。=1 表示 VKA 是垂向导水率的乘数 laywet = 0 # 在 UPW 包中始终设置为零 hk = 20 # 水平导水率 (m/d 或 ft/d，取决于 lenuni) # hani = 1 # 如果 CHANI > 1，则不需要 vka = 0.5 # 垂直导水率 (m/d 或 ft/d) ss = 0.00001 # 比储存量 (1/ft 或 1/m) sy = 0.20 # 比出水度 (无量纲) upw = flopy.modflow.ModflowUpw( mf, laytyp=laytyp, layavg=layavg, chani=chani, layvka=layvka, laywet=laywet, ipakcb=53, # 将单元格流量保存到单元号 53 的文件 hdry=hdry, iphdry=iphdry, hk=hk, vka=vka, ss=ss, sy=sy, )实例化 MODFLOW-NWT 的基本 (BAS 或 BA6) 包
Python
fname = "ibnd_lay1.txt" ibnd1_pth = data_path / folder_name / "bas_arrays" / fname pooch.retrieve( url=f"https://github.com/modflowpy/flopy/raw/develop/examples/data/{folder_name}/bas_arrays/{fname}", fname=fname, path=data_path / folder_name / "bas_arrays", known_hash=None, ) ibnd1 = np.loadtxt(ibnd1_pth) # 加载第一层的 ibound 数据 ibnd2 = np.ones(ibnd1.shape) # 第二层 ibound 均为 1 (活动单元) ibnd3 = np.ones(ibnd2.shape) # 第三层 ibound 均为 1 (活动单元) ibnd = [ibnd1, ibnd2, ibnd3] # 所有层的 ibound 列表 ibnd = np.array(ibnd) # 转换为 NumPy 数组 fname = "strthd1.txt" StHd1_pth = data_path / folder_name / "bas_arrays" / fname pooch.retrieve( url=f"https://github.com/modflowpy/flopy/raw/develop/examples/data/{folder_name}/bas_arrays/{fname}", fname=fname, path=data_path / folder_name / "bas_arrays", known_hash=None, ) StHd1 = np.loadtxt(StHd1_pth) # 加载第一层的初始水头数据 fname = "strthd2.txt" StHd2_pth = data_path / folder_name / "bas_arrays" / fname pooch.retrieve( url=f"https://github.com/modflowpy/flopy/raw/develop/examples/data/{folder_name}/bas_arrays/{fname}", fname=fname, path=data_path / folder_name / "bas_arrays", known_hash=None, ) StHd2 = np.loadtxt(StHd2_pth) # 加载第二层的初始水头数据 fname = "strthd3.txt" StHd3_pth = data_path / folder_name / "bas_arrays" / fname pooch.retrieve( url=f"https://github.com/modflowpy/flopy/raw/develop/examples/data/{folder_name}/bas_arrays/{fname}", fname=fname, path=data_path / folder_name / "bas_arrays", known_hash=None, ) StHd3 = np.loadtxt(StHd3_pth) # 加载第三层的初始水头数据 strtElev = [StHd1, StHd2, StHd3] # 所有层的初始水头列表 strtElev = np.array(strtElev) # 转换为 NumPy 数组 hdry = 999.0 # 无水头单元的水头值 bas = flopy.modflow.ModflowBas(mf, ibound=ibnd, hnoflo=hdry, strt=strtElev)实例化 MODFLOW-NWT 的一般水头边界 (GHB) 包
Python
# GHB 边界位于区域的顶部 (北) 和底部 (南) 边缘，所有 3 层。 elev_stpt_row1 = 308.82281 # 第 1 行 (北边界) 的起始高程 elev_stpt_row300 = 239.13811 # 第 300 行 (南边界) 的起始高程 elev_slp = (308.82281 - 298.83 9) / (ncol - 1) # 高程坡度 (这里可能数据有误，通常应该是两个点的差值除以距离) # 注释：原始代码中的 elev_slp 计算可能有一个小错误，它使用的是两个高程值和 ncol-1，但这两个高程值看起来不是直接对应于边界两端的。 # 实际含义可能是：高程变化量 / 列数 - 1 sp = [] # 用于存储 GHB 数据的列表 # 遍历所有层 for k in [0, 1, 2]: # 层索引 (0-based) # 遍历边界行 for i in [0, 299]: # 行索引 (0-based)，代表最北和最南的行 # 遍历所有列 for j in np.arange(0, 300, 1): # 列索引 (0-based) # 跳过不满足 下面内容条件的单元格 # 这里的条件似乎是为了排除某些特定的单元格，而不是所有的边界单元格都设置GHB if (i == 1 and (j < 27 or j > 31)) or (i == 299 and (j < 26 or j > 31)): if i % 2 == 0: # 如果是偶数行 (这里 i=0 是偶数行) # 添加 GHB 数据: [层, 行, 列, 水头, 电导] sp.append([k, i, j, elev_stpt_row1 - (elev_slp * (j - 1)), 11.3636]) else: # 如果是奇数行 (这里 i=299 是奇数行) sp.append( [k, i, j, elev_stpt_row300 - (elev_slp * (j - 1)), 11.3636] ) # 为南边界的特定列范围添加 GHB (补充之前循环中可能遗漏的边界) for k in [0, 1, 2]: # 所有层 for j in np.arange(26, 32, 1): # 列 26 到 31 sp.append([k, 299, j, 238.20, 3409.0801]) # 添加 GHB 数据到最南边的特定列 ghb = flopy.modflow.ModflowGhb(mf, stress_period_data=sp) # 创建 GHB 包对象注释：这里的 GHB 定义比较复杂，包含了一些特定的行、列和计算逻辑。关键在于每个 GHB 单元格需要提供其层、行、列索引，以及对应的水头和电导值。
实例化 MODFLOW-NWT 的地表水流路径 (SFR2) 包
Python
# 使用 numpy.genfromtxt() 将预先准备好的河段数据读取到 numpy recarrays 中 # 请记住，存储在预先准备好的 NO3_ReachInput.csv 文件中的单元格索引是基于 0 的。 # Flopy 在写入文件时会将其转换为基于 1 的索引。 fname = "no3_reachinput.csv" rpth = data_path / folder_name / "sfr_data" / fname pooch.retrieve( url=f"https://github.com/modflowpy/flopy/raw/develop/examples/data/{folder_name}/sfr_data/{fname}", fname=fname, path=data_path / folder_name / "sfr_data", known_hash=None, ) reach_data = np.genfromtxt(rpth, delimiter=",", names=True) # 加载河段数据 reach_data # 显示加载的河段数据 (recarray 形式) # 使用 numpy.genfromtxt() 将预先准备好的河段数据读取到 numpy recarrays 中 fname = "no3_segmentdata.csv" spth = data_path / folder_name / "sfr_data" / fname pooch.retrieve( url=f"https://github.com/modflowpy/flopy/raw/develop/examples/data/{folder_name}/sfr_data/{fname}", fname=fname, path=data_path / folder_name / "sfr_data", known_hash=None, ) ss_segment_data = np.genfromtxt(spth, delimiter=",", names=True) # 加载分段数据 # segment_data 一个字典，键是应力期索引，值是该应力期的分段数据 segment_data = {0: ss_segment_data, 1: ss_segment_data} segment_data[0][0:1]["width1"] # 访问第一个应力期第一个分段的 width1 (宽度) nstrm = len(reach_data) # 河段数 nss = len(segment_data[0]) # 分段数 (取第一个应力期的分段数) nsfrpar = 0 # SFR 参数数量 (这里没有使用参数化， 因此为 0) const = 128390.4 # 曼宁方程常数，单位为 cfs (立方英尺/秒) dleak = 0.0001 # 排水泄漏系数 (单位通常是 L/T，如 1/d) ipakcb = 53 # 将 SFR 输出写入单元格流量文件 (单元号 53) istcb2 = 37 # 将 SFR 输出写入文 这篇文章小编将件 (单元号 37) isfropt = 1 # SFR 选项。1 表示使用 reachinput 格式 # dataset 5 (请参阅在线指南) (ITMP, IRD , IPT ) dataset_5 = { 0: [nss, 0, 0], # 第一个应力期：nss 分段，不读/不打印文件 1: [-1, 0, 0], # 第二个应力期：-1 表示沿用上一个应力期的设置 } # 输入参数通常遵循在线 MODFLOW 指南中定义的变量名称 sfr = flopy.modflow.ModflowSfr2( mf, nstrm=nstrm, # 河段数 nss=nss, # 分段数 const=const, # 曼宁方程常数 dleak=dleak, # 排水泄漏系数 ipakcb=ipakcb, # 流量输出到单元格流量文件 istcb2=istcb2, # 流量输出到文 这篇文章小编将件 isfropt=isfropt, # SFR 选项 reachinput=True, # 使用 reachinput 格式 reach_data=reach_data, # 河段数据 segment_data=segment_data, # 分段数据 dataset_5=dataset_5, # dataset 5 unit_number=15, # SFR 包的单元号 )
好的，我将继续准确翻译并添加注释。
实例化 MODFLOW-NWT 的湖泊 (LAK) 包
Python
# 读取预先准备好的湖泊数组 fname = "lakarr1.txt" LakArr_pth = data_path / folder_name / "lak_arrays" / fname pooch.retrieve( # 使用pooch下载文件 url=f"https://github.com/modflowpy/flopy/raw/develop/examples/data/{folder_name}/lak_arrays/{fname}", fname=fname, path=data_path / folder_name / "lak_arrays", known_hash=None, ) LakArr_lyr1 = np.loadtxt(LakArr_pth) # 加载第一层的湖泊定义数组 (包含湖泊单元格编号) LakArr_lyr2 = np.zeros(LakArr_lyr1.shape) # 第二层没有湖泊单元，全部为0 LakArr_lyr3 = np.zeros(LakArr_lyr2.shape) # 第三层没有湖泊单元，全部为0 LakArr = [LakArr_lyr1, LakArr_lyr2, LakArr_lyr3] # 湖泊定义数组，每层一个 LakArr = np.array(LakArr) # 转换为NumPy数组 nlakes = int(np. x(LakArr)) # 湖泊数量，根据LakArr中的最大值确定 ipakcb = ipakcb # 从上方继承的输出控制标志，表示将湖泊流量保存到单元格预算文件 theta = -1.0 # 隐式 技巧 (-1.0 表示显式，-1.0 to 1.0 之间) nssitr = 10 # 牛顿法求解平衡湖泊水位的最大迭代次数 sscncr = 1.000e-03 # 平衡湖泊水位解的收敛准则 surfdep = 2.000e00 # 湖底拓扑变化的微小高度 (用于计算与地下水的交换面积) stages = 268.00 # 运行开始时每个湖泊的初始水位 # ITMP > 0, 读取湖泊定义数据 (由lakarr提供) # ITMP1 ≥ 0, 读取每个湖泊新的补给、蒸发、径流和取水数据 (由flux_data提供) # LWRT > 0, 抑制湖泊包的输出打印 bdlknc_lyr1 = LakArr_lyr1.copy() # 湖泊-地下水连接标志， 第一层湖泊定义 bdlknc_lyr2 = LakArr_lyr1.copy() # 假设第二层和第一层有相同的连接模式 bdlknc_lyr3 = np.zeros(LakArr_lyr1.shape) # 第三层没有直接湖泊连接 # 需要将 bdlknc_lyr1 中非零值向各个 路线 (左/右和上/下) 扩展 1 个单元格 # 这是为了确保湖泊边界处的单元格能够正确地与湖泊连接 for i in np.arange(0, LakArr_lyr1.shape[0]): for j in np.arange(0, LakArr_lyr1.shape[1]): im1 = i - 1 ip1 = i + 1 jm1 = j - 1 jp1 = j + 1 # 检查上方单元格 if im1 >= 0: if LakArr_lyr1[i, j] == 1 and LakArr_lyr1[im1, j] == 0: bdlknc_lyr1[im1, j] = 1 # 如果当前是湖泊且上方是非湖泊，则将上方标记为连接单元 # 检查下方单元格 if ip1 < LakArr_lyr1.shape[0]: if LakArr_lyr1[i, j] == 1 and LakArr_lyr1[ip1, j] == 0: bdlknc_lyr1[ip1, j] = 1 # 如果当前是湖泊且下方是非湖泊，则将下方标记为连接单元 # 检查左侧单元格 if jm1 >= 0: if LakArr_lyr1[i, j] == 1 and LakArr_lyr1[i, jm1] == 0: bdlknc_lyr1[i, jm1] = 1 # 如果当前是湖泊且左侧是非湖泊，则将左侧标记为连接单元 # 检查右侧单元格 if jp1 < LakArr_lyr1.shape[1]: if LakArr_lyr1[i, j] == 1 and LakArr_lyr1[i, jp1] == 0: bdlknc_lyr1[i, jp1] = 1 # 如果当前是湖泊且右侧是非湖泊，则将右侧标记为连接单元 bdlknc = [bdlknc_lyr1, bdlknc_lyr2, bdlknc_lyr3] # 所有层的湖泊-地下水连接标志列表 bdlknc = np.array(bdlknc) # 转换为NumPy数组 # 流量数据：{应力期索引: [[湖泊索引, 流量类型索引, 流量值]]} # [湖泊索引, 流量类型索引, 流量值, ...] # 流量类型索引: 1=降雨, 2=蒸发, 3=径流, 4=取水 # 这里只有湖泊 0 的降雨和蒸发 flux_data = { 0: [[0.0073, 0.0073, 0.0, 0.0]], # 应力期 0 的流量数据 (lake 0, 降雨, 蒸发, 径流, 取水) 1: [[0.0073, 0.0073, 0.0, 0.0]] # 应力期 1 的流量数据 } lak = flopy.modflow.ModflowLak( mf, nlakes=nlakes, # 湖泊数量 ipakcb=ipakcb, # 输出控制标志 theta=theta, # 时刻权重 nssitr=nssitr, # 牛顿法迭代次数 sscncr=sscncr, # 收敛准则 surfdep=surfdep, # 湖底拓扑变化高度 stages=stages, # 初始水位 lakarr=LakArr, # 湖泊定义数组 bdlknc=bdlknc, # 湖泊-地下水连接标志 flux_data=flux_data, # 流量数据 unit_number=16, # LAK 包的单元号 )实例化用于 MODFLOW-NWT 包的测站 (GAGE) 包
Python
gages = [ [1, 225, 90, 3], # [测站编号, 行, 列, 层] [2, 68, 91, 3], [3, 33, 92, 3], [4, 165, 93, 3], [5, 123, 94, 3], [6, 77, 95, 3], [7, 173, 96, 3], [8, 328, 97, 3], [9, 115, 98, 3], [-1, -101, 1], # 这 一个 独特的测站定义，可能是针对整个流域或其他汇总输出 ] # gages = [[1,38,61,1],[2,67,62,1], [3,176,63,1], [4,152, ,1], [5,186,65,1], [6,31,66,1]] # 注释掉的旧测站定义 files = [ "no3.gag", # 主测站输出文件 "seg1_gag.out", # 各分段测站输出文件 "seg2_gag.out", "seg3_gag.out", "seg4_gag.out", "seg5_gag.out", "seg6_gag.out", "seg7_gag.out", "seg8_gag.out", "seg9_gag.out", "lak1_gag.out", # 湖泊测站输出文件 ] numgage = len(gages) # 测站数量 gage = flopy.modflow.ModflowGage(mf, numgage=numgage, gage_data=gages, filenames=files)实例化 MODFLOW-NWT 的非饱和带流 (UZF) 包
Python
nuztop = 2 # 非饱和带顶层标识。2表示顶层是 的，即降雨和蒸发直接 影响于该层。 iuzfopt = 2 # 非饱和带计算选项。2表示启用地表径流和泉水排放。 irun = 1 # 径流计算标志。1表示计算地表径流。 iet = 0 # 蒸发计算标志。0表示不计算蒸发 (通常在MODFLOW中由EVT包处理)。 iuzfcb = 52 # 将 UZF 输出保存到单元格预算文件 (单元号 52)。 iuzfcb2 = 0 # 不将 UZF 输出保存到辅助文 这篇文章小编将件。 ntrail2 = 20 # 跟踪粒子数量的尾迹单元格数。 nsets2 = 20 # 每次调用 NWT 求解器时计算流量的 UZF 单元格组数。 nuzgag = 2 # UZF 测站数量。 surfdep = 2.0 # 地表凹陷深度 (用于计算地表径流的起始水深)。 eps = 3.0 # 地表径流计算中的一个参数 (通常与粗糙度有关)。 thts = 0.30 # 饱和含水量。 thti = 0.13079 # 初始含水量。 fname = "iuzbnd.txt" fname_uzbnd = data_path / folder_name / "uzf_arrays" / fname pooch.retrieve( url=f"https://github.com/modflowpy/flopy/raw/develop/examples/data/{folder_name}/uzf_arrays/{fname}", fname=fname, path=data_path / folder_name / "uzf_arrays", known_hash=None, ) fname = "irunbnd.txt" fname_runbnd = data_path / folder_name / "uzf_arrays" / fname pooch.retrieve( url=f"https://github.com/modflowpy/flopy/raw/develop/examples/data/{folder_name}/uzf_arrays/{fname}", fname=fname, path=data_path / folder_name / "uzf_arrays", known_hash=None, ) iuzfbnd = np.loadtxt(fname_uzbnd) # 加载 UZF 边界标志数组 irunbnd = np.loadtxt(fname_runbnd) # 加载径流边界标志数组 uzgag = [[106, 160, 121, 3], [1, 1, -122, 1]] # UZF 测站数据: [测站编号, 行, 列, 层] # finf = {0: 1.8250e-03, 1: 1.8250e-03} # 注释掉的入渗率数据，可能通过外部文件或 uzf_data 提供 uzf = flopy.modflow.ModflowUzf1( mf, nuztop=nuztop, # 非饱和带顶层标识 iuzfopt=iuzfopt, # 非饱和带计算选项 irun =irun , # 径流计算标志 iet =iet , # 蒸发计算标志 ipakcb=iuzfcb, # 输出控制标志 iuzfcb2=iuzfcb2, # 辅助输出控制标志 ntrail2=ntrail2, # 尾迹单元格数 nsets=nsets2, # 流量计算组数 surfdep=surfdep, # 地表凹陷深度 uzgag=uzgag, # UZF 测站数据 iuzfbnd=1, # UZF 边界标志 (这里直接设置为1，可能表示所有的iuzfbnd单元都参与UZF计算) irunbnd=0, # 径流边界标志 (这里直接设置为0，表示不使用额外的径流边界定义) vks=1.0e-6, # 垂向饱和导水率 eps=3.5, # 径流参数 thts=0.35, # 饱和含水量 thtr=0.15, # 残余含水量 thti=0.20, # 初始含水量 )实例化 MODFLOW-NWT 的排水 (DRN) 包
Python
fname = "elv.txt" fname_drnElv = data_path / folder_name / "drn_arrays" / fname pooch.retrieve( url=f"https://github.com/modflowpy/flopy/raw/develop/examples/data/{folder_name}/drn_arrays/{fname}", fname=fname, path=data_path / folder_name / "drn_arrays", known_hash=None, ) fname = "cond.txt" fname_drnCond = data_path / folder_name / "drn_arrays" / fname pooch.retrieve( url=f"https://github.com/modflowpy/flopy/raw/develop/examples/data/{folder_name}/drn_arrays/{fname}", fname=fname, path=data_path / folder_name / "drn_arrays", known_hash=None, ) drnElv = np.loadtxt(fname_drnElv) # 加载排水单元高程数据 drnCond = np.loadtxt(fname_drnCond) # 加载排水单元电导数据 # 将 NumPy 数组转换为 Pandas DataFrame，方便处理非零值 drnElv_lst = pd.DataFrame( { "lay": 1, # 排水单元位于第 1 层 (0-based 对应实际的第2层) "row": np.nonzero(drnElv)[0] + 1, # 获取非零元素的行索引，并转换为 1-based "col": np.nonzero(drnElv)[1] + 1, # 获取非零元素的列索引，并转换为 1-based "elv": drnElv[np.nonzero(drnElv)], # 获取非零元素的高程值 "cond": drnCond[np.nonzero(drnCond)], # 获取非零元素的电导值 }, columns=["lay", "row", "col", "elv", "cond"], # 定义列名 ) # 将 DataFrame 转换为列表的列表，以便 drn 构造函数使用 stress_period_data = drnElv_lst.values.tolist() # 创建一个字典，每个应力期有一个条目 stress_period_data = {0: stress_period_data, 1: stress_period_data} # 两个应力期使用相同的排水数据 drn = flopy.modflow.ModflowDrn(mf, ipakcb=ipakcb, stress_period_data=stress_period_data)实例化与质量传输路由 (LMT) 包的链接 (为 MODFLOW-NWT 生成链接文件)
Python
lmt = flopy.modflow.ModflowLmt( mf, output_file_name="NO3.ftl", # 输出流路径文件名称 output_file_header="extended", # 输出文件头格式为扩展格式 output_file_for t="for tted", # 输出文件格式为格式化文本 package_flows=["all"], # 将所有包的流量信息写入 linker 文件 )现在开始创建 MT3D-USGS 文件
Python
# 首先设置 MT3D-USGS 类，并传入 MODFLOW-NWT 对象以设置多个 BTN 数组 mt = flopy.mt3d.Mt3dms( modflowmodel=mf, # 关联的 MODFLOW 模型 modelname=modelname, # MT3D 模型名称 model_ws=modelpth, # MT3D 模型 职业空间 version="mt3d-usgs", # MT3D-USGS 版本 namefile_ext="mtnam", # MT3D 名称文件的扩展名 exe_name=mtexe, # MT3D-USGS 可执行文件名称 ftlfilename="no3.ftl", # MODFLOW 流量文件名称 (由 LMT 生成) ftlfree=True, # 流量文件是否为 free 格式 )实例化 MT3D-USGS 的基本传输 (BTN) 包
Python
ncomp = 1 # 组分数量 lunit = "FT" # 长度单位 (英尺) sconc = 0.0 # 初始浓度 (所有单元格的默认初始浓度) prsity = 0.3 # 孔隙度 cinact = -1.0 # 非活动单元格的浓度 (通常用于非活动单元格，例如 -1.0) thkmin = 0.000001 # 最小有效层厚度 nprs = -2 # 打印浓度 结局的频率。-2 表示在每个 时刻步结束时打印。 nprobs = 10 # 质量平衡计算的打印频率。 npr s = 10 # 全局质量平衡计算的打印频率。 dt0 = 0.1 # 初始传输 时刻步长 nstp = 1 # 每个流 时刻步内的传输 时刻步数 (这里每个流步只有一个传输步) mxstrn = 500 # 最大应变值 (用于自动调整 时刻步长) tt ult = 1.2 # 时刻步长乘数 (传输 时刻步长可以按此乘数增加) tt ax = 100 # 最大传输 时刻步长 # 这些观测点需要用 0-based 索引输入 obs = [[0, 104, 158], [1, 104, 158], [2, 104, 158]] # 观测点数据: [层, 行, 列] btn = flopy.mt3d.Mt3dBtn( mt, MFStyleArr=True, # 使用 MODFLOW 风格的数组输入 DRYCell=True, # 允许干燥单元格 lunit=lunit, # 长度单位 sconc=sconc, # 初始浓度 ncomp=ncomp, # 组分数量 prsity=prsity, # 孔隙度 cinact=cinact, # 非活动单元格浓度 obs=obs, # 观测点 thkmin=thkmin, # 最小厚度 nprs=nprs, # 打印频率 nprobs=nprobs, # 质量平衡打印频率 chk s=True, # 启用质量平衡检查 npr s=npr s, # 质量平衡打印频率 dt0=dt0, # 初始传输 时刻步长 nstp=nstp, # 每个流 时刻步内的传输 时刻步数 mxstrn=mxstrn, # 最大应变值 tt ult=tt ult, # 时刻步长乘数 tt ax=tt ax, # 最大传输 时刻步长 )实例化 MT3D-USGS 的平流 (ADV) 包
Python
mixelm = 0 # 平流项处理 技巧。0 表示上游加权有限差分法 (Upstream weighting) percel = 1.0000 # 允许在每个传输 时刻步内通过单元格的质量分数 (0.0 到 1.0) mxpart = 5000 # 粒子跟踪 技巧中的最大粒子数量 (如果 mixelm 为 3)。 nadvfd = 1 # 平流方案的选择标志。1 = 上游加权有限差分。 adv = flopy.mt3d.Mt3dAdv(mt, mixelm=mixelm, percel=percel, mxpart=mxpart, nadvfd=nadvfd)实例化 MT3D-USGS 的广义共轭梯度 (GCG) 求解器包
Python
mxiter = 1 # 最大外部迭代次数 (用于非线性 难题) iter1 = 50 # 内部迭代次数 isolve = 3 # 预处理共轭梯度法 (PCG) 求解器类型。3 表示直接稀疏求解器。 ncrs = 0 # 非交叉预条件器标志。0表示不使用。 accl = 1.000000 # 加速参数 (对于 PCG，通常设置为 1.0)。 cclose = 1.00e-06 # 停止迭代的收敛容差。 iprgcg = 5 # 打印 PCG 求解器诊断信息。 gcg = flopy.mt3d.Mt3dGcg( mt, mxiter=mxiter, iter1=iter1, isolve=isolve, ncrs=ncrs, accl=accl, cclose=cclose, iprgcg=iprgcg, )实例化 MT3D-USGS 的弥散 (DSP) 包
Python
al = 0.1 # 纵向弥散度 (单位：长度) trpt = 0.1 # 水平横向弥散度与 'AL' 的比率 (无量纲) trpv = 0.1 # 垂直横向弥散度与 'AL' 的比率 (无量纲) dmcoef = 1.0000e-10 # 分子扩散系数 (单位：长度^2/ 时刻) dsp = flopy.mt3d.Mt3dDsp(mt, al=al, trpt=trpt, trpv=trpv, dmcoef=dmcoef, multiDiff=True) # multiDiff=True: 允许使用每个单元格的弥散系数 (如果提供的话)实例化 MT3D-USGS 的源汇混合 (S ) 包
Python
# 没有与边界条件相关的用户指定浓度 mxss = 11199 # 源汇项的最大数量 s = flopy.mt3d.Mt3dS (mt, mxss=mxss)实例化 MT3D-USGS 的反应 (RCT) 包
Python
isothm = 0 # 吸附等温线类型。0表示无吸附。 ireact = 1 # 反应类型标志。1表示一阶衰减反应。 irctop = 2 # 反应选项。2表示使用单独的一阶衰减速率 rc1 和 rc2。 igetsc = 0 # 是否从文件读取初始吸附浓度。0表示不读取。 ireaction = 0 # 反应类型标志，通常在早期版本中使用，新版本中ireact更常用。 rc1 = 6.3258e-04 # 溶解相的一阶反应速率 (1/ 时刻) rc2 = 0.0 # 土壤层上的衰减率 (这里设置为0，表示无固相衰减) rct = flopy.mt3d.Mt3dRct( mt, isothm=isothm, ireact=ireact, igetsc=igetsc, rc1=rc1, rc2=rc2 )实例化 MT3D-USGS 的地表水流传输 (SFT) 包
Python
nsfinit = len(reach_data) # 地表水流包中包含的初始河段数量 mxsfbc = len(reach_data) # 地表水流包中的最大边界条件数量 (这里与河段数相同) icbcsf = 0 # 是否将地表水流传输的流量信息写入文件。0表示不写。 ioutobs = 92 # 输出观测值的间隔 (每 ioutobs 个 时刻步输出一次) isfsolv = 1 # SFT 求解器选项。1表示使用直接求解器。 wimp = 0.5 # 时刻加权因子 (介于 0 和 1 之间，0.5表示Crank-Nicolson) wups = 1.0 # 上游加权因子 (通常为 1.0) cclosesf = 1.0e-6 # SFT 求解器的收敛容差 mxitersf = 10 # SFT 求解器的最大迭代次数 crntsf = 1.0 # Courant 数限制 (用于自动调整 时刻步长) iprtxmd = 0 # 打印 SFT 诊断信息标志。0表示不打印。 coldsf = 0 # 初始浓度。0表示不指定。 dispsf = 0 # 弥散选项。0表示不启用弥散 (在 SFT 中)。 obs_sf = [225, 293, 326, 491, 614, 691, 8 , 1192, 1307] # SFT 观测点河段编号 (1-based) sf_stress_period_data = { 0: [0, 0, 0], # 应力期 0 的 SFT 边界条件数据 (这里为默认值) 1: [0, 0, 0], # 应力期 1 的 SFT 边界条件数据 2: [0, 0, 0], # 应力期 2 的 SFT 边界条件数据 (即使 MODFLOW 只有2个应力期，这里可能为预留) } gage_output = [None, None, "no3.sftobs"] # 测站输出文件列表。第三个元素用于 SFT 观测点输出。 sft = flopy.mt3d.Mt3dSft( mt, nsfinit=nsfinit, mxsfbc=mxsfbc, icbcsf=icbcsf, ioutobs=ioutobs, isfsolv=isfsolv, wimp=wimp, wups=wups, cclosesf=cclosesf, mxitersf=mxitersf, crntsf=crntsf, iprtxmd=iprtxmd, coldsf=coldsf, dispsf=dispsf, nobssf=len(obs_sf), # 观测点数量 obs_sf=obs_sf, # 观测点列表 sf_stress_period_data=sf_stress_period_data, # 应力期数据 filenames=gage_output, # 测站输出文件 )实例化 MT3D-USGS 的非饱和带传输 (UZT) 包
Python
mxuzcon = np.count_nonzero(irunbnd) # UZT 边界条件的最大数量 (基于径流边界的非零单元格数) icbcuz = 45 # UZT 输出单元格预算文件编号 iet = 0 # 蒸发计算标志。0表示不启用 UZT 中的蒸发。 wc = np.ones((nlay, nrow, ncol)) * 0.29 # 初始含水量分布 (通常用于UZT内部计算) sdh = np.ones((nlay, nrow, ncol)) # 地表滞留水深 (或地表排水高程，取决于具体应用) uzt = flopy.mt3d.Mt3dUzt( mt, mxuzcon=mxuzcon, # UZT 边界条件的最大数量 icbcuz=icbcuz, # 输出单元格预算文件编号 iet=iet, # 蒸发计算标志 iuzfbnd=iuzfbnd, # UZF 边界标志数组 (从 MODFLOW UZF 包继承) sdh=sdh, # 地表滞留水深 cuzinf=1.4158e-03, # 入渗水浓度 (或入渗流量) filenames="no3", # UZT 输出文件的前缀 )实例化 MT3D-USGS 的湖泊传输 (LKT) 包
Python
nlkinit = 1 # 初始湖泊数量 mxlkbc = 720 # 湖泊边界条件的最大数量 icbclk = 81 # 湖泊传输输出单元格预算文件编号 ietlak = 1 # 湖泊蒸发标志。1表示计算蒸发。 coldlak = 1 # 初始湖泊浓度标志。1表示所有湖泊初始浓度为 1.0 (或指定值)。 # 湖泊流量数据: {应力期索引: [[湖泊索引, 流量类型索引, 流量值]]} # 流量类型索引: 1=入流, 2=蒸发, 3=降雨, 4=径流, 5=取水, 6=地下水交换 (取决于MODFLOW-LAK的输出) # 这里 [[0, 1, 0.01667]] 表示湖泊 0，类型 1 (入流)，值为 0.01667 lkt_flux_data = { 0: [[0, 1, 0.01667]], # 应力期 0 的湖泊流量数据 1: [[0, 1, 0.02667]] # 应力期 1 的湖泊流量数据 (入流值不同) } lkt = flopy.mt3d.Mt3dLkt( mt, nlkinit=nlkinit, # 初始湖泊数量 mxlkbc=mxlkbc, # 最大湖泊边界条件数量 icbclk=icbclk, # 输出单元格预算文件编号 ietlak=ietlak, # 湖泊蒸发标志 coldlak=coldlak, # 初始湖泊浓度标志 lk_stress_period_data=lkt_flux_data, # 湖泊流量数据 )写入 MT3D-USGS 输入文件以进行检查和运行
Python
mf.write_input() # 写入 MODFLOW 输入文件 mt.write_input() # 写入 MT3D-USGS 输入文件Python
mf.run_model() # 运行 MODFLOW 模型 mt.run_model() # 运行 MT3D-USGS 模型使用 FloPy 简化 MT3DMS S 包的使用
这 一个多组分传输的演示。
Python
import os import sys from tempfile import TemporaryDirectory import numpy as np import flopy print(sys.version) print(f"numpy version: {np.__version__}") print(f"flopy version: {flopy.__version__}")输出:
3.12.10 | packaged by conda-forge | ( in, Apr 10 2025, 22:21:13) [GCC 13.3.0] numpy version: 2.2.5 flopy version: 3.9.3首先，我们将创建一个简单的模型结构
Python
nlay, nrow, ncol = 10, 10, 10 # 定义模型的层数、行数和列数 perlen = np.zeros((10), dtype=float) + 10 # 每个应力期的长度都设置为10个 时刻单位 nper = len(perlen) # 应力期总数 ibound = np.ones((nlay, nrow, ncol), dtype=int) # 定义活动区域，所有单元格都设置为活动单元 (1) botm = np.arange(-1, -11, -1) # 定义每层底部高程，从-1到-10 top = 0.0 # 定义模型顶部高程创建 MODFLOW 包
Python
# 临时目录 temp_dir = TemporaryDirectory() model_ws = temp_dir.name # 模型 职业空间路径 modelname = "s ex" # 模型名称 # 实例化 MODFLOW 对象 mf = flopy.modflow.Modflow(modelname, model_ws=model_ws) # 实例化离散化 (DIS) 包 dis = flopy.modflow.ModflowDis( mf, nlay=nlay, nrow=nrow, ncol=ncol, perlen=perlen, # 每个应力期的长度 nper=nper, # 应力期总数 botm=botm, # 各层底部高程 top=top, # 模型顶部高程 steady=False, # 所有应力期都设置为非稳定流 (瞬态流) ) # 实例化基本 (BAS) 包 # ibound: 定义活动、非活动和定水头单元格 # strt: 定义初始水头 bas = flopy.modflow.ModflowBas(mf, ibound=ibound, strt=top) # 实例化层属性流 (LPF) 包 lpf = flopy.modflow.ModflowLpf(mf, hk=100, vka=100, ss=0.00001, sy=0.1) # hk: 水平导水率 # vka: 垂直各向异性 (这里与垂直导水率相同) # ss: 比储存量 (针对承压层) # sy: 比出水度 (针对非承压层) # 实例化输出控制 (OC) 包 oc = flopy.modflow.ModflowOc(mf) # 使用默认输出控制设置 # 实例化预处理共轭梯度 (PCG) 求解器包 pcg = flopy.modflow.ModflowPcg(mf) # 使用默认求解器设置 # 实例化补给 (RCH) 包 rch = flopy.modflow.ModflowRch(mf) # 即使没有补给数据，也实例化以便 S 包识别我们将使用 MODFLOW 边界条件来跟踪 S 数据单元格位置。
获取一个字典 (itype)，其中包含每种边界类型的 S itype。
Python
itype = flopy.mt3d.Mt3dS .itype_dict() # 获取 MT3DMS S 包的 ityep 字典，它映射 MODFLOW 包类型到 S 源汇类型 print(itype) print(flopy.mt3d.Mt3dS .get_default_dtype()) # 打印 S 包默认的数据类型 s _data = {} # 初始化 S 数据字典输出:
{'CHD': 1, 'BAS6': 1, 'PBC': 1, 'WEL': 2, 'DRN': 3, 'RIV': 4, 'GHB': 5, 'MAS': 15, 'CC': -1} [('k', '<i8'), ('i', '<i8'), ('j', '<i8'), ('css', '<f4'), ('itype', '<i8')]注释：itype_dict 显示了不同 MODFLOW 边界条件类型对应的 S itype 值。例如，WEL (井) 对应 itype 2，GHB (一般水头边界) 对应 itype 5。s .get_default_dtype() 显示了 S 数据所需的基本列：层(k)、行(i)、列(j)、源汇浓度(css) 和 源汇类型(itype)。
添加一般水头边界 (ghb)。
一般水头边界的水头 (bhead) 在前 5 个应力期为 0.1，其中组分 1 (comp_1) 浓度为 1.0，组分 2 (comp_2) 浓度为 100.0。 接着 bhead 增加到 0.25，comp_1 浓度降低到 0.5，comp_2 浓度增加到 200.0。
Python
ghb_data = {} # 初始化 GHB 边界数据字典 print(flopy.modflow.ModflowGhb.get_default_dtype()) # 打印 GHB 包的默认数据类型 # 应力期 0 的 GHB 数据 # (层, 行, 列, 边界水头, 电导) ghb_data[0] = [(4, 4, 4, 0.1, 1.5)] # 在单元格 (4,4,4) 处设置一个 GHB # 应力期 0 的 S 数据 (对应 ghb_data[0] 中的 GHB) # (层, 行, 列, 源汇浓度 (通用), 源汇类型, 组分 1 浓度, 组分 2 浓度) # 注意：css 是通用浓度，但对于多组分模型，实际浓度通过后续的 cs (N) 参数传递 s _data[0] = [(4, 4, 4, 1.0, itype["GHB"], 1.0, 100.0)] # 应力期 5 的 GHB 数据 (从应力期 5 开始，边界条件发生变化) ghb_data[5] = [(4, 4, 4, 0.25, 1.5)] # 应力期 5 的 S 数据 (对应 ghb_data[5] 中的 GHB) s _data[5] = [(4, 4, 4, 0.5, itype["GHB"], 0.5, 200.0)] # 为所有层的第一列添加 GHB 边界，并设置浓度为 0 for k in range(nlay): for i in range(nrow): ghb_data[0].append((k, i, 0, 0.0, 100.0)) # 在第一列设置 GHB，水头为 0.0 s _data[0].append((k, i, 0, 0.0, itype["GHB"], 0.0, 0.0)) # 在第一列设置 S 浓度为 0 # 重新定义应力期 5 的 GHB 数据和 S 数据，覆盖之前为 (4,4,4) 设置的值 # 这是为了确保这个循环定义的边界也包含在应力期 5 的 S 数据中 ghb_data[5] = [(4, 4, 4, 0.25, 1.5)] # 再次添加单元格 (4,4,4) 的 GHB s _data[5] = [(4, 4, 4, 0.5, itype["GHB"], 0.5, 200.0)] # 再次添加单元格 (4,4,4) 的 S for k in range(nlay): for i in range(nrow): ghb_data[5].append((k, i, 0, -0.5, 100.0)) # 在应力期 5 的第一列设置 GHB，水头为 -0.5 s _data[5].append((k, i, 0, 0.0, itype["GHB"], 0.0, 0.0)) # 在应力期 5 的第一列设置 S 浓度为 0输出:
[('k', '<i8'), ('i', '<i8'), ('j', '<i8'), ('bhead', '<f4'), ('cond', '<f4')]注释：这里展示了 GHB 包数据结构。 由于 S 包的数据结构是字典，其键是应力期编号。每个键对应的值是列表，其中包含该应力期内每个源汇项的数据。对于多组分，MT3DMS S 包使用 cs (N) 字段来指定每个组分的浓度。
添加注入井。
所有应力期中，注入速率 (flux) 为 10.0，组分 1 浓度为 10.0，组分 2 浓度为 0.0。警告: 由于我们更改了应力期 6 的 S 数据（ 由于应力期 5 开始改变，会影响到 6），我们需要将井也添加到应力期 6 的 s _data 中（即 s _data[5]）。
Python
wel_data = {} # 初始化 WEL 井数据字典 print(flopy.modflow.ModflowWel.get_default_dtype()) # 打印 WEL 包的默认数据类型 # 应力期 0 的 WEL 数据 wel_data[0] = [(0, 4, 8, 10.0)] # 在单元格 (0,4,8) 处设置一个井，注入速率 10.0 # 应力期 0 的 S 数据，添加对应 WEL 的浓度 # 注意：css 是通用浓度，但对于多组分模型，实际浓度通过后续的 cs (N) 参数传递 s _data[0].append((0, 4, 8, 10.0, itype["WEL"], 10.0, 0.0)) # 应力期 5 的 S 数据，添加对应 WEL 的浓度（ 由于边界条件在应力期 5 变化， 因此需要明确添加） s _data[5].append((0, 4, 8, 10.0, itype["WEL"], 10.0, 0.0))输出:
[('k', '<i8'), ('i', '<i8'), ('j', '<i8'), ('flux', '<f4')]注释：这里展示了 WEL 包数据结构。对于 s _data，即使 WEL 的 flux 在应力期 5 之后没有变化， 由于其他边界条件在应力期 5 发生了变化，为了确保 S 文件正确生成，通常需要为发生变化的应力期重新完整定义 S 数据。
将 GHB 和 WEL 包添加到 mf MODFLOW 对象实例。
Python
ghb = flopy.modflow.ModflowGhb(mf, stress_period_data=ghb_data) # 创建 GHB 包实例 wel = flopy.modflow.ModflowWel(mf, stress_period_data=wel_data) # 创建 WEL 包实例创建 MT3DMS 包
Python
# 实例化 MT3DMS 主对象 mt = flopy.mt3d.Mt3dms(modflowmodel=mf, modelname=modelname, model_ws=model_ws) # 实例化基本传输 (BTN) 包 # sconc: 初始浓度 (对于多组分，这是第一个组分的默认初始浓度) # ncomp: 组分数量 (这里是 2 个组分) # sconc2: 第二个组分的初始浓度。如果有更多组分，可以定义 sconc3, sconc4 等。 btn = flopy.mt3d.Mt3dBtn(mt, sconc=0, ncomp=2, sconc2=50.0) # 实例化平流 (ADV) 包 adv = flopy.mt3d.Mt3dAdv(mt) # 实例化源汇混合 (S ) 包 # stress_period_data: 传入之前准备好的 S 边界条件数据字典 s = flopy.mt3d.Mt3dS (mt, stress_period_data=s _data) # 实例化广义共轭梯度 (GCG) 求解器包 gcg = flopy.mt3d.Mt3dGcg(mt)输出:
found 'rch' in modflow model, resetting crch to 0.0 S : setting crch for component 2 to zero. kwarg name crch2注释：这些警告来自 S 包，它检测到 MODFLOW 模型中存在 RCH 包，并自动将 RCH 的浓度 (crch) 设置为 0.0，特别是对于组分 2。这通常是 Flopy 内部处理不同包之间连接的机制。
让我们验证 stress_period_data 是否具有正确的数据类型
Python
print(s .stress_period_data.dtype) # 打印 S 包应力期数据的数据类型输出:
[('k', '<i8'), ('i', '<i8'), ('j', '<i8'), ('css', '<f4'), ('itype', '<i8'), ('cs (01)', '<f4'), ('cs (02)', '<f4')]注释：输出显示了 S 数据类型的详细结构。除了基本的 k, i, j, css, itype 之外，还出现了 cs (01) 和 cs (02) 字段，这正是 MT3DMS 处理多组分浓度的方式，证明 Flopy 正确地为 2 个组分配置了 S 数据。
创建 SEAWAT 包
Python
# 实例化 SEAWAT 对象，它结合了 MODFLOW 和 MT3DMS swt = flopy.seawat.Seawat( modflowmodel=mf, mt3dmodel=mt, modelname=modelname, namefile_ext="nam_swt", # SEAWAT 名称文件的扩展名 model_ws=model_ws, ) # 实例化可变密度流 (VDF) 包 (SEAWAT 特有) vdf = flopy.seawat.SeawatVdf(swt, mtdnconc=0, iwtable=0, indense=-1) # mtdnconc: 0 表示不使用浓度相关密度 # iwtable: 0 表示不使用密度表 # indense: -1 表示密度从 VDF 输入文件读取，而不是在 VDF 中计算写入输入文件
Python
mf.write_input() # 写入 MODFLOW 输入文件 mt.write_input() # 写入 MT3DMS 输入文件 swt.write_input() # 写入 SEAWAT 输入文件最后，修改 vdf 包以修复 indense。
Python
fname = f"{modelname}.vdf" # VDF 文件名 f = open(os.path.join(model_ws, fname)) lines = f.readlines() # 读取 VDF 文件的所有行 f.close() f = open(os.path.join(model_ws, fname), "w") # 以写入模式重新打开文件 for line in lines: f.write(line) # 写入原始行 for kper in range(nper): f.write("-1 ") # 对于每个应力期，写入 -1。这对应于 INDENSE 参数，指示读取外部密度 f.close()注释：这里对 .vdf 文件进行了手动修改。在 SEAWAT 中，当 indense = -1 时，它期望在每个应力期之后有一行 -1，表示密度来自 MT3DMS 的输出。Flopy 在生成时可能未完全符合这一要求，因此需要手动添加。
清理临时 职业空间
Python
try: # 忽略 Windows 上的 PermissionError temp_dir.cleanup() # 清理并删除临时目录 except: passSEAWAT 教程 1：Henry 盐水入侵 难题
在本教程中，我们将使用 FloPy 创建、运行和后处理 Henry 盐水入侵 难题，使用 SEAWAT 版本 4。
Python
# ## 入门 from pathlib import Path from tempfile import TemporaryDirectoryPython
import numpy as npPython
import flopyHenry 难题的输入变量
Python
Lx = 2.0 # 模型在 x 路线的长度 (米) Lz = 1.0 # 模型在 z 路线（垂向）的长度 (米) nlay = 50 # 模型层数 nrow = 1 # 模型行数 (2D 模拟， 因此只有 1 行) ncol = 100 # 模型列数 delr = Lx / ncol # 列宽 (米) delc = 1.0 # 行宽 (对于 2D 难题，此值通常不重要，设为 1.0) delv = Lz / nlay # 每层垂直厚度 (米) henry_top = 1.0 # 模型顶面高程 henry_botm = np.linspace(henry_top - delv, 0.0, nlay) # 各层底面高程，从顶面减去层厚向下均匀分布到 0.0 qinflow = 5.702 # 入流流量 (立方米/天) dmcoef = 0.57024 # 分子扩散系数 (平方米/天)；也可以尝试 1.62925 作为 Henry 难题的另一个案例 hk = 8 .0 # 水力传导系数 (米/天)创建基本的 MODFLOW 模型结构
Python
temp_dir = TemporaryDirectory() workspace = temp_dir.name # 模型 职业空间路径 (临时目录) name = "seawat_henry" # 模型名称 # 实例化 SEAWAT 主对象。SEAWAT 结合了 MODFLOW 和 MT3DMS。 swt = flopy.seawat.Seawat(name, exe_name="swtv4", model_ws=workspace) print(swt.namefile) # 打印 SEAWAT 的名称文件路径输出:
seawat_henry.nam保存单元格流量到单元号 53
Python
ipakcb = 53 # 输出控制标志，将单元格流量保存到单元号 53 的文件将 DIS 包添加到 MODFLOW 模型
Python
dis = flopy.modflow.ModflowDis( swt, nlay, # 层数 nrow, # 行数 ncol, # 列数 nper=1, # 应力期数 delr=delr, # 列宽 delc=delc, # 行宽 laycbd=0, # 无隔水层 top=henry_top, # 顶面高程 botm=henry_botm, # 各层底面高程 perlen=1.5, # 应力期长度 nstp=15, # 应力期内的 时刻步数 )Python
# BAS 包的变量 ibound = np.ones((nlay, nrow, ncol), dtype=np.int32) # 初始 ibound 数组，所有单元格设为活动单元 (1) ibound[:, :, -1] = -1 # 最右侧一列的单元格设为定水头边界 (-1)将 BAS 包添加到 MODFLOW 模型
Python
bas = flopy.modflow.ModflowBas(swt, ibound, 0) # 实例化 BAS 包，初始水头设为 0将 LPF 包添加到 MODFLOW 模型
Python
lpf = flopy.modflow.ModflowLpf(swt, hk=hk, vka=hk, ipakcb=ipakcb) # hk: 水平水力传导系数 # vka: 垂直水力传导系数 (这里与 hk 相同，表示各向同性) # ipakcb: 单元格预算输出标志将 PCG 包添加到 MODFLOW 模型
Python
pcg = flopy.modflow.ModflowPcg(swt, hclose=1.0e-8) # 实例化 PCG 求解器，水头收敛容差设为 1.0e-8将 OC 包添加到 MODFLOW 模型
Python
oc = flopy.modflow.ModflowOc( swt, stress_period_data={(0, 0): ["save head", "save budget"]}, # 在第一个应力期的第一个 时刻步保存水头和预算 compact=True, # 紧凑输出格式 )创建 WEL 和 S 数据
Python
itype = flopy.mt3d.Mt3dS .itype_dict() # 获取 S 包的源汇类型字典 wel_data = {} # 井数据字典 s _data = {} # 源汇混合数据字典 wel_sp1 = [] # 第一个应力期（索引 0）的井数据列表 s _sp1 = [] # 第一个应力期（索引 0）的 S 数据列表 for k in range(nlay): # 遍历所有层 # 井数据：在最左侧一列 (j=0) 添加注入井，流量平均分配到每层 wel_sp1.append([k, 0, 0, qinflow / nlay]) # S 数据： # 1. 对应井的浓度：最左侧一列 (j=0) 的浓度为 0.0 (淡水) s _sp1.append([k, 0, 0, 0.0, itype["WEL"]]) # 2. 定浓度边界：最右侧一列 (j=ncol-1) 的浓度为 35.0 (盐水) s _sp1.append([k, 0, ncol - 1, 35.0, itype["BAS6"]]) # 使用 BAS6 ityep 表示定浓度边界 wel_data[0] = wel_sp1 # 将井数据添加到字典，键为应力期索引 0 s _data[0] = s _sp1 # 将 S 数据添加到字典，键为应力期索引 0 wel = flopy.modflow.ModflowWel(swt, stress_period_data=wel_data, ipakcb=ipakcb) # 实例化 WEL 包，并指定 ipakcb 以便将井流量保存到预算文件创建基本的 MT3DMS 模型结构
Python
btn = flopy.mt3d.Mt3dBtn( swt, nprs=-5, # 浓度输出频率。-5 表示在每个应力期的结束 时刻步输出 prsity=0.35, # 孔隙度 sconc=35.0, # 初始浓度（这里设为盐水浓度，但在后续单元格会重新赋值） ifmtcn=0, # 浓度输出文件格式。0 表示二进制非格式化 chk s=False, # 不进行质量平衡检查（为简化起见） nprobs=10, # 质量平衡计算的打印频率 npr s=10, # 质量平衡计算 结局打印频率 dt0=0.001, # 初始传输 时刻步长 ) adv = flopy.mt3d.Mt3dAdv(swt, mixelm=0) # 实例化平流 (ADV) 包，mixelm=0 表示上游加权有限差分法 dsp = flopy.mt3d.Mt3dDsp(swt, al=0.0, trpt=1.0, trpv=1.0, dmcoef=dmcoef) # 实例化弥散 (DSP) 包 # al: 纵向弥散度 (设为 0.0，表示无机械弥散) # trpt: 水平横向弥散度与 al 的比率 (这里无意义， 由于 al=0) # trpv: 垂直横向弥散度与 al 的比率 (这里无意义， 由于 al=0) # dmcoef: 分子扩散系数 gcg = flopy.mt3d.Mt3dGcg(swt, iter1=500, mxiter=1, isolve=1, cclose=1e-7) # 实例化广义共轭梯度 (GCG) 求解器包 # iter1: 内部迭代次数 # mxiter: 外部迭代次数 # isolve: 求解器类型 (1 表示直接稀疏求解器) # cclose: 浓度收敛容差 s = flopy.mt3d.Mt3dS (swt, stress_period_data=s _data) # 实例化源汇混合 (S ) 包创建 SEAWAT 模型结构
Python
vdf = flopy.seawat.SeawatVdf( swt, iwtable=0, # 不使用密度表 densemin=0, # 最小密度（这里设为 0， 由于 denseref 和 denseslp 会定义） dense x=0, # 最大密度（这里设为 0） denseref=1000.0, # 参考密度 (淡水密度，公斤/立方米) denseslp=0.7143, # 密度斜率 (密度随浓度变化的斜率，公斤/立方米 per 单位浓度) firstdt=1e-3, # 初始密度相关的流动 时刻步长 )写入输入文件
Python
swt.write_input() # 将所有 SEAWAT 模型的输入文件写入到 职业空间运行模型
Python
success, buff = swt.run_model(silent=True, report=True) # 运行 SEAWAT 模型，静默模式并报告输出 assert success, "SEAWAT did not terminate nor lly." # 断言模型运行成功，否则报错后处理 结局
Python
import numpy as np # 再次导入 numpy，确保可用Python
import flopy.utils.binaryfile as bf # 导入用于读取 MODFLOW/MT3DMS 二进制文件的工具加载浓度数据
Python
ucnobj = bf.UcnFile(Path(workspace) / "MT3D001.UCN", model=swt) # 创建浓度文件对象 times = ucnobj.get_times() # 获取所有可用的 时刻步 concentration = ucnobj.get_data(totim=times[-1]) # 获取 最后一个 时刻步的浓度数据加载逐单元格流量数据
Python
cbbobj = bf.CellBudgetFile(Path(workspace) / f"{name}.cbc") # 创建单元格预算文件对象 times = cbbobj.get_times() # 获取所有可用的 时刻步 qx = cbbobj.get_data(text="flow right face", totim=times[-1])[0] # 获取 最后一个 时刻步的右侧流量 (x 路线) qy = np.zeros((nlay, nrow, ncol), dtype=float) # 对于 2D 模型，y 路线流量为 0 qz = cbbobj.get_data(text="flow lower face", totim=times[-1])[0] # 获取 最后一个 时刻步的下侧流量 (z 路线)创建包含浓度和流速向量的图
Python
import tplotlib.pyplot as plt # 导入 MatplotlibPython
fig = plt.figure(figsize=(12, 9)) # 创建一个图表 ax = fig.add_subplot(1, 1, 1, aspect="equal") # 添加一个子图，并设置纵横比相等 # 创建横截面绘图对象，这里是沿着第 0 行（ 由于只有 1 行） pmv = flopy.plot.PlotCrossSection(model=swt, ax=ax, line={"row": 0}) arr = pmv.plot_array(concentration) # 绘制浓度数组 pmv.plot_vector(qx, qy, -qz, color="white", kstep=3, hstep=3) # 绘制流速向量。注意 qz 前加负号， 由于 MODFLOW qz 通常是向下为正，但我们希望向量向上表示水流向上 # kstep 和 hstep 控制向量的绘制密度 plt.colorbar(arr, shrink=0.5, ax=ax) # 添加颜色条 ax.set_title("Simulated Concentrations") # 设置图 深入了解输出:
Text(0.5, 1.0, 'Simulated Concentrations')
加载水头数据
Python
headobj = bf.HeadFile(Path(workspace) / f"{name}.hds") # 创建水头文件对象 times = headobj.get_times() # 获取所有可用的 时刻步 head = headobj.get_data(totim=times[-1]) # 获取 最后一个 时刻步的水头数据创建包含水头等值线的图
Python
fig = plt.figure(figsize=(12, 9)) # 创建一个图表 ax = fig.add_subplot(1, 1, 1, aspect="equal") # 添加一个子图，并设置纵横比相等 pmv = flopy.plot.PlotCrossSection(model=swt, ax=ax, line={"row": 0}) # 创建横截面绘图对象 arr = pmv.plot_array(head) # 绘制水头数组 contours = pmv.contour_array(head, colors="white") # 绘制水头等值线，颜色为白色 ax.clabel(contours, fmt="%2.2f") # 在等值线上添加标签，格式为两位小数 plt.colorbar(arr, shrink=0.5, ax=ax) # 添加颜色条 ax.set_title("Simulated Heads") # 设置图 深入了解输出:
Text(0.5, 1.0, 'Simulated Heads')
Python
try: temp_dir.cleanup() # 尝试清理临时目录 except: # 避免 Windows 权限错误 pass