# -*- coding: utf-8 -*- import flopy import numpy as np import matplotlib.pyplot as plt import os # 初始化模型域 model_name = "mf6_model" mf_exe = r"E:\modFlow\mf6.exe" model_dir = os.path.join(os.getcwd(), "modflow_example") os.makedirs(model_dir, exist_ok=True) sim = flopy.mf6.MFSimulation( sim_name=model_name, version="mf6", exe_name=mf_exe, sim_ws=model_dir ) # 时间配置（10年周期） tdis = flopy.mf6.ModflowTdis( sim, time_units="DAYS", nper=1, # 单应力期 perioddata=[(3650 * 10, 100, 1.0)] # 总时长3650天*10年，100个时间步 ) # 水流模型（GWF） gwf = flopy.mf6.ModflowGwf( sim, modelname=model_name+"_flow", save_flows=True, newtonoptions=["NEWTON"] ) # 空间离散（三维非结构化网格示例） nlay = 3 # 3个含水层 bot = [50.0, 30.0, 10.0] # 各层底板高程 dis = flopy.mf6.ModflowGwfdis( gwf, nlay=nlay, nrow=100, ncol=100, delr=10.0, # 10米网格分辨率 delc=10.0, top=60.0, botm=bot ) # 含水层参数（基于实验室测定与现场抽水试验） k11 = [1e-4, 5e-5, 1e-5] # 水平渗透系数(m/s) ss = [1e-5, 1e-5, 1e-5] # 储水率(1/m) sy = [0.2, 0.15, 0.1] # 给水度 npf = flopy.mf6.ModflowGwfnpf( gwf, icelltype=[1, 0, 0], # 1:可压缩非承压，0:刚性承压 k=k11, k33=[k*0.1 for k in k11] # 垂向各向异性系数 ) # 边界条件（定水头边界） chd = flopy.mf6.ModflowGwfchd( gwf, stress_period_data=[ [(0, 0, 0), 55.0], # 西边界水头55m [(0, 99, 99), 50.0] # 东边界水头50m ] ) # 溶质运移模型（GWT） gwt = flopy.mf6.ModflowGwt( sim, modelname=model_name+"_transport" ) # 溶质离散参数（基于ISO 5667-3:2018标准） dsp = flopy.mf6.ModflowGwtdsp( gwt, alh=10.0, # 纵向弥散度(m) ath1=1.0, # 横向弥散度 atv=0.1 # 垂向弥散度 ) # 运移过程配置（吸附+降解） mst = flopy.mf6.ModflowGwtmst( gwt, sorption="linear", bulk_density=1.8e3, # 介质密度(kg/m³) distcoef=0.5e-3, # 分配系数(m³/kg) decay=0.001, # 一级降解速率(1/day) decay_sorbed=True ) # 初始污染源（突发泄漏场景） source_loc = (8588952, 20002690,13 ), # 中间层中心位置 ssm = flopy.mf6.ModflowGwtssm( gwt, sources=[ [ # 第一层：所有应力期 ( # 第二层：单个应力期的源项 0, # 应力

<think>好的，我现在需要帮助用户解决如何使用FloPy构建包含水流和溶质运移的MODFLOW 6模型的问题。首先，我得回顾一下MODFLOW 6和FloPy的基本知识。MODFLOW 6是USGS开发的最新地下水流动模型，支持多种类型的模拟，包括地下水流动和溶质运移。FloPy是用于构建和运行MODFLOW模型的Python包，可以方便地进行前处理和后处理。 用户的问题涉及几个关键步骤：空间离散化、渗透系数定义、定水头边界条件，以及溶质弥散和降解参数的配置。我需要分步骤来解答每个部分，并确保使用正确的FloPy模块和类。 首先，空间离散化。MODFLOW 6使用结构化网格，通常通过DIS或DISV包来定义。FloPy中创建模型的第一步通常是建立模拟对象，然后添加时间步长和应力期，接着是空间离散化。比如，使用DIS包时，需要指定层数、每层的行数和列数，以及网格的尺寸。需要提醒用户注意网格的分辨率是否足够，以及是否会影响计算效率和精度。 接下来是渗透系数的定义。这属于地下水流动模型的参数，通常在NPF（Node Property Flow）包中设置。用户需要为每个模型层指定水平渗透系数（K）和垂向各向异性系数。在FloPy中，可以通过npf包的k属性来赋值，可能是一个标量或数组，取决于各层的均质性。 然后是定水头边界条件。这通常通过CHD（Constant Head）包来实现，用户需要指定哪些单元格是定水头，以及对应的水头值。需要指导用户如何创建CHD包，并正确设置应力期数据，包括开始和结束时间，以及各时间段的水头值。注意MODFLOW 6中应力期的管理方式，可能需要使用应力期输入的结构。 关于溶质运移部分，用户需要配置溶质弥散和降解参数。这涉及到GWT（GroundWater Transport）模型中的DSP（Dispersion）和IST（Mass Source Term）包。DSP包需要设置纵向弥散度、横向弥散度和垂向弥散度，可能还包括分子扩散系数。降解参数可能通过化学反应包来定义，比如一级衰减系数。需要确认用户是否需要进行反应运移模拟，或者只是简单的衰减模型。 此外，需要注意MODFLOW 6中水流模型（GWF）和运移模型（GWT）的耦合。用户需要分别建立这两个模型，然后通过交换（Exchange）将它们连接起来。在FloPy中，可能需要使用gwf与gwt模型的构造函数，并创建相应的交换对象。 可能遇到的常见问题包括：网格划分不正确导致模型不稳定，边界条件设置错误导致无法收敛，弥散参数单位不匹配，或者溶质运移模型未正确耦合。需要提醒用户检查输入数据的单位一致性，以及模型的时间步长设置是否合理。 最后，应该提供一个代码示例，展示如何将这些步骤整合起来。代码需要包括导入必要的模块、创建模拟对象、定义空间离散化、设置渗透系数、添加定水头边界、配置溶质运移参数，以及运行模拟。同时，代码中的注释应该清晰，帮助用户理解每一步的作用。 总结来说，我需要按照MODFLOW 6和FloPy的工作流程，分步骤解释每个关键配置，并提供示例代码。同时，生成的相关问题应该覆盖用户可能进一步感兴趣的方面，如参数敏感性分析、结果可视化、处理非均质性等。</think>以下是使用FloPy构建MODFLOW 6水流与溶质运移模型的详细步骤： ### 1. 空间离散化配置 使用$DIS$包定义结构化网格，假设建立3层模型，网格尺寸为$100 \times 100 \times 10$米： ```python from flopy.mf6 import MFSimulation sim = MFSimulation(sim_name="contaminant_transport") tdis = flopy.mf6.ModflowTdis(sim, time_units="DAYS") gwf = flopy.mf6.ModflowGwf(sim, modelname="gwf") dis = flopy.mf6.ModflowGwfdis( gwf, nlay=3, nrow=50, ncol=50, delr=100.0, delc=100.0, top=0.0, botm=[-10, -20, -30] ) ``` ### 2. 渗透系数定义 通过$NPF$包设置渗透系数场，假设各向同性介质： ```python npf = flopy.mf6.ModflowGwfnpf( gwf, icelltype=0, # 0表示饱和层 k=[1e-4, 5e-5, 2e-5], # 各层渗透系数(m/s) save_flows=True ) ``` ### 3. 定水头边界条件 使用$CHD$包设置左右边界为定水头： ```python chd_spd = [] for layer in range(3): for row in range(50): # 左边界 chd_spd.append([(layer, row, 0), 10.0]) # 水头10m # 右边界 chd_spd.append([(layer, row, 49), 8.0]) # 水头8m chd = flopy.mf6.ModflowGwfchd( gwf, stress_period_data=chd_spd, pname="CHD-1" ) ``` ### 4. 溶质运移模型配置 #### 4.1 创建GWT模型 ```python gwt = flopy.mf6.ModflowGwt(sim, modelname="gwt") # 耦合GWF-GWT gwfgwt = flopy.mf6.ModflowGwfgwt( sim, exgtype="GWF6-GWT6", exgmnamea=gwf.name, exgmnameb=gwt.name ) ``` #### 4.2 弥散参数设置 ```python dsp = flopy.mf6.ModflowGwtdsp( gwt, alh=10.0, # 纵向弥散度(m) ath1=1.0, # 横向弥散度(m) atv=0.1, # 垂向弥散度(m) diffc=1e-9 # 分子扩散系数(m²/s) ) ``` #### 4.3 污染物初始浓度 ```python ic = flopy.mf6.ModflowGwtic( gwt, strt=0.0, # 初始浓度全为0 filename="contam.ic" ) # 设置污染源区域 cnc_spd = {0: [[(0, 25, 25), 1000.0]]} # 顶层中心点初始浓度1000mg/L ssm = flopy.mf6.ModflowGwtssm( gwt, sources=cnc_spd ) ``` #### 4.4 一级降解反应 ```python ist = flopy.mf6.ModflowGwtist( gwt, decay_rate=0.001 # 降解速率常数(1/day) ) ``` ### 5. 完整模型执行流程 ```python # 设置求解器参数 ims = flopy.mf6.ModflowIms( sim, complexity="MODERATE" ) # 写入输入文件并运行 sim.write_simulation() success, _ = sim.run_simulation() if not success: raise Exception("模型运行失败") ```