Modflow6模型PEST自动调参

项目图示

视频教程

info

本页面中的短视频为简短操作介绍，获取完整EnviFusion教学课程，请访问 在线课程。

算例工程下载地址

1 算例简介 2 操作步骤 2.1 新建项目 2.2 创建MODFLOW模型 2.3 定义活动区域 2.4 定义网格高程 2.5 定义渗透系数 2.6 定义补给 2.7 定义定水头 2.8 定义抽水井 2.9 定义排水沟 2.10 解释运行MODFLOW 2.11 创建MT3DMS模型 2.12 解释运行MT3DMS模型 3 小结

算例简介

本算例先使用水平渗透初始值创建MODFLOW6水流模型，再基于10个观测井的观测水位值，使用PEST调参方法，对MODFLOW6管道模型中的水平渗透系数进行耦合迭代计算，以期得到适宜的参数值。

研究区域为一个中等大小盆地，该盆地面积有72.5平方公里，当地为半干旱气候，年平均降水量为0.381 m/yr，多数降水通过蒸散流失。含水层的补给最终排入到盆地中心的一条小溪流中。 这条溪流的流向向北，最终流入到高程为1000米的湖中。盆地中有十口观测井。

操作步骤

新建项目

1. 打开Envifusion软件，若软件界面仍有其他项目，点击菜单栏按钮关闭当前项目文件。

创建模型网格

首先，基于模拟区边界shp文件创建用于运行MODFLOW6模拟的网格，并进行模型初始化。

1. 点击开始菜单栏中的导入数据选项，或点击界面上方的按钮，打开算例文件夹中的Sources & Sinks_ploy.shp文件，在该图层属性面板中点击应用完成设置，生成Sources & Sinks_ploy.shp图层。

2. 选中Sources & Sinks_ploy.shp图层，点击计算->通用工具箱->曲面->曲面->曲面（均一平滑），生成曲面->曲面（均一平滑）1图层。在该图层属性面板数据子面板中，将目标边长设为1000，迭代次数设为2；在渲染子面板中，将渲染方式设为曲面网络，点击应用基于研究区边界生成三角网格。

   

   

3. 选中曲面->曲面（均一平滑）1图层，点击通用工具箱->散点->散点->曲面（泰森），生成散点->曲面（泰森）1图层。在该图层属性面板点击应用基于上一步中的三角网格顶点构建泰森多边形。

   

4. 选中散点->曲面（泰森）1图层，点击通用工具箱->属性->采样，以曲面->曲面（均一平滑）1图层为数据源，以散点->曲面（泰森）1图层为采样体，点击确定生成采样1图层。

   

5. 选中采样1图层，点击通用工具箱->网格->阈值裁切，生成阈值裁切1图层。在该图层属性面板数据子面板中，将标量数列设为有效采样，阈值范围设为1至1点击应用裁取有效采样网格，生成用于MODFLOW6模拟的多边形网格。

   

   

6. 点击计算->数值模拟->MODFLOW6水流->初始化ModFlow6模型，生成初始化ModFlow6模型1图层。将默认水平渗透系数设为8，点击应用完成设置。

   用户应自行确认模型单位的一致性，单位会在后续操作中进行设置。本算例长度单位为米，时间单位为天。

   

定义网格高程

导入包含模型网格顶底高程的csv文件定义网格高程。

1. 点击开始菜单栏中的导入数据选项，或点击界面上方的按钮，打开算例文件夹中top.csv文件，生成top.csv1图层，在该图层属性面板中点击应用完成导入。

2. 选择top.csv1图层，点击工具箱->散点->表格->散点，生成表格->散点1图层，在该图层属性面板数据子面板中，将XYZ坐标分别设为xyz列，点击应用完成设置。

   

3. 选初始化ModFlow6模型1图层，点击数值模拟-> MODFLOW6水流->MF6网格高程，选择表格->散点1图层为属性图层，点击确认生成MF6网格高程1图层，将其重命名为010-MF6网格高程1图层。

   在010-MF6网格高程1图层属性面板数据子面板中，将赋值模式设为散点Z坐标->单层，赋值对象设为模型顶板，勾选下方的上层优先选项，点击应用完成设置，为网格顶部设置高程。

   

4. 重复步骤1-3，导入包含网格底部高程的bot.csv文件，为模型网格赋值底部高程，生成020-MF6网格高程2图层。在该图层属性面板数据子面板中，将赋值模式设为散点Z坐标->单层，赋值对象设为单层底板，点击应用完成设置，为网格赋值底部高程。在3D模式下查看模型网格如下图所示：

   

定义渗透系数

本项目将基于水位观测值对单层模型的水平渗透系数K值分布进行PEST自动调参。在开始创建MODFLOW6模型时，先在研究区范围内确定一定数量的种子散点。准备的散点csv文件包含散点坐标、各点预设渗透系数参数、渗透系数调参范围等属性值。导入该文件后，先将种子散点中的预设K值分布映射到MODFLOW6模型，再耦合PEST进行参数调整。

1. 点击开始菜单栏中的导入数据选项，或点击界面上方的按钮，打开算例文件夹中K参数.csv文件，生成K参数.csv1图层，在该图层属性面板中点击应用完成设置。导入模型渗透系数文件。

   点击视图窗口上方的新建按钮，新建视图2表格窗口。在表格窗口中可看到导入的csv文件内容。该文件将作为PEST调参的入口文件，因此应至少包含散点XY坐标、current（模型中使用的当前参数值，在调参开始时等于初始值）、initial（参数初始值）、min（参数调参范围最小值）、max（参数调参范围最大值）、parTrans（调参拟合选项，取值为0表示该参数进行线性调参，大于0表示使用Log调参，小于0表示不参与调参）、type（参数分组标记，调整多类参数时使用，本例中不涉及）等属性列。

   

2. 选择K参数.csv1图层，点击工具箱->散点->表格->散点，生成表格->散点3图层。在该图层属性面板数据子面板中，为X、Y坐标分别选择x和y数列，勾选下方的2D散点和保留坐标选项，点击应用完成设置。

   

3. 选择020-MF6网格高程2图层，点击数值模拟-> MODFLOW6水流->MF6克里金赋值，选择表格->散点3图层为属性图层，点击确认生成MF6克里金赋值1图层，将其重命名为030-K-克里金赋值1图层。

   点击该图层属性面板顶部的显示/隐藏高级属性选项，在数据子面板中将属性类别选为水平渗透系数，层号为1，水平渗透系数选择current数列。在渲染子面板中将渲染方式改为曲面网络，点击应用预设模型水平渗透系数。

   

定义定水头边界

绘制折线定义模型定水头边界。

1. 在模型最北部绘制表示定水头边界的折线。点击计算->创建数据体->几何体->多边形，生成多边形1图层。按住Alt键点击鼠标左键在模型顶部绘制折线，穿过模型顶部一行网格，在该图层属性面板数据子面板中，点击坐标参数属性下删除按钮，删除（0，0，0）和（1，0，0）两点，点击应用完成绘制。

   

2. 选择030-K-克里金赋值1图层，点击数值模拟-> MODFLOW6水流->MF6线赋值，选择多边形1图层为属性图层，点击确认生成MF6线赋值1图层，将其重命名为040-定水头-线赋值1图层。

   在该图层数据子面板中将属性类别选为定水头边界，层号为1，固定水头选择指定值，指定值设为1000，在渲染子面板中将渲染方式改为曲面网络，点击应用设置模型定水头边界。

   

   

定义河流边界

导入包含河流边界条件的efd文件，设置河流边界。

1. 点击开始菜单栏中的导入数据选项，或点击界面上方的按钮，打开算例文件夹中河流.efd文件，生成河流.efd1图层，在该图层属性面板中点击应用完成设置。导入包含河流边界信息的折线，所包含的属性包括河流水位、河底高程、河流传导系数。

   

2. 选择040-定水头-线赋值1图层，点击数值模拟-> MODFLOW6水流->MF6线赋值，选择河流.efd1图层为属性图层，点击确认生成MF6线赋值2图层，将其重命名为050-河流-线赋值2图层。

   在该图层数据子面板中将属性类别选为河流，层号为1，河流水头、传导系数和河底高程分别选择响应属性列，在渲染子面板中将渲染方式改为曲面网络，点击应用设置模型河流边界。

   

   

定义补给边界

为模型赋值补给边界（垂向入渗边界）。

1. 选择050-河流-线赋值2图层，点击视图窗口上方的矩形选择网格按钮，选中视图窗口中所有网格，点击数值模拟-> MODFLOW6水流->MF6网格赋值，生成MF6网格赋值1图层，将其重命名为060-补给-网格赋值1图层。

   在该图层属性面板数据子面板中将属性类别选为垂向入渗边界，层号为1，垂向入渗率设为8e-5，点击应用完成设置，为模型设置垂向入渗（降水补给）边界。

   

定义观测井

导入包含观测井信息的csv文件设置观测井。

1. 点击开始菜单栏中的导入数据选项，或点击界面上方的按钮，打开算例文件夹中obswell.csv文件，生成obswell.csv1图层，在该图层属性面板中点击应用完成设置。导入包含观测井信息的文件。在表格窗口中可看到导入的csv文件内容，观测井文件至少要包含观测点XY坐标，观测水位三列属性。

   

2. 选择obswell.csv1图层，点击工具箱->散点->表格->散点，生成表格->散点4图层。在该图层属性面板数据子面板中，为X、Y坐标分别选择x和y数列，勾选下方的2D散点和保留坐标选项，点击应用完成设置。

   

3. 选择060-补给-网格赋值1图层，点击数值模拟-> MODFLOW6水流->MF6点赋值，选择表格->散点4图层为属性图层，点击确认生成MF6点赋值1图层，将其重命名为070-观测井-点赋值1图层。

   在该图层数据子面板中将属性类别选为观测站，测站类型设为水位观测井，在渲染子面板中将渲染方式改为曲面网络，点击应用设置模型观测井。

   

   

解释运行MODFLOW

解释运行创建好的水流模型并读取模型结果。

1. 选中070-观测井-点赋值1图层，点击数值模拟-> MODFLOW6水流->MF6解释运行，生成MF6解释运行1图层。在属性面板中，点击按钮选择模型文件生成的位置，将模型名保存为Pest。选择复杂度为线性，时间单位为天，长度单位为米，点击应用，系统将自动编译并运行MODFLOW6模型。

   

2. 点击开始->保存工程，将当前项目文件保存为pest调参.efp。重新点击EnviFusion图标打开一个新工程界面，将其保存为水头残差.efp，用于查看模型结果和观测井水头残差分布。

3. 在水头残差.efp新工程中，点击开始菜单栏中的导入数据选项，或点击界面上方的按钮，选择上一步中用户指定名称的Pest.mf6文件夹中的efModFlow6.nam6文件。项目浏览器中生成efModFlow6.nam6图层，点击应用导入MODFLOW6模型结果。

   在该图层属性面板渲染子面板中，将渲染方式设为曲面，调色板色彩设为散点ef水头，点击下方的选择预设调色板，选择RainbowRich调色板，将风格属性下的透明度调整为0.25；在视窗属性下，将背景设为纯色，RGB值分别为232，视图窗口如下图所示：

   

   

查看水头残差

导入观测井信息，基于观测水位和模型计算水位，计算并使用圆柱展示水头残差。

1. 在新工程水头残差.efp中导入obswell.csv文件，并使用表格->散点工具，生成表格->散点1图层。

2. 选中efModFlow6.nam6图层，点击通用工具箱->空间拓扑->散点穿刺，将efModFlow6.nam6图层设为待穿刺图层，将表格->散点1图层设为散点图层，点击确定生成散点穿刺1图层，在该图层属性面板数据子面板中，将穿刺模式设为内部单点，点击应用获得模型网格中观测井所在散点。

   

3. 选中散点穿刺1图层，点击通用工具箱->属性->采样，将efModFlow6.nam6图层设为数据源，将散点穿刺1图层设为采样体，点击确定生成采样1图层，在该图层属性面板数据子面板中，勾选保留散点数据选项，点击应用采样模型模拟水位数据。查看该图层信息面板，可看到该图层包含ef水头和观测水头两列属性值。

   

   

4. 选中采样1图层，点击通用工具箱->属性->计算器，生成计算器1图层。在该图层属性面板数据子面板中，将结果数列名称设为ef水头残差，点击计算器下方的标量按钮，将计算公式设为ef水头-观测水头，点击应用完成残差计算。

   

5. 选中计算器1图层，点击通用工具箱-> 散点->散点->线段（高级拉伸），生成散点->线段（高级拉伸）1图层。在属性面板数据子面板中，将线段长度设为ef水头残差，长度缩放系数设为100，点击应用完成设置。

6. 选中散点->线段（高级拉伸）1图层，点击通用工具箱-> 散点->散点->线段（高级拉伸），生成散点->线段（高级拉伸）1图层。在属性面板数据子面板中，将半径设为150；渲染子面板中，将色彩设为ef水头残差，点击应用可在视图中看到表示残差的圆柱。

   

7. 圆柱向上残差为正，模拟水头高于观测水头；圆柱向下残差为负，模拟水头低于观测水头。视图窗口如下图所示：可看出按照水平渗透系数预设值模拟所得水位结果与观测水位相比偏低。下一步对K值进行PEST自动调参。

   

设置PEST自动调参

在模型管道中加入PEST调参图层，准备调参运行文件，基于观测井水位数据调整水平渗透系数值。

1. 在pest调参.efp项目中，选择MF6解释运行1图层，点击数值模拟-> MODFLOW6水流->自动调参，选择K参数.csv图层为参数图层，MF6解释运行1图层为模型图层，点击确认生成自动调参1图层。

   

2. 在该图层属性面板数据子面板中，设置控制文件.pst的保存路径，本算例保存在新建的算例文件夹\PEST文件夹中；同时用户应自行准备观测文件：obs.csv文件，将鼠标置于观测文件选项处会弹出文件说明。

   将外部迭代设为5。Pest调参中存在内部迭代和外部迭代两种过程，每次外部迭代将包含多次内部迭代过程，其内部迭代次数由参与拟合的参数数量决定。外部迭代次数越多，优化拟合的程度越高，但同时也相应消耗计算资源。用户在优化过程中，可先尝试使用较低的外部迭代次数（如3-5次），后期根据拟合效果和计算时间灵活调整。

3. 自行准备的obs.csv文件格式如下图所示，包含两个必要字段列，model_obs（模型计算值）和real_obs（观测水位值，真值），如需设置不同观测井权重，可增加weight字段列。文件可包含多行数据，每行将各自指定一套观测文件映射规则（文件保存路径）。

   本例中的只有一组观测井，因此观测文件仅包含一行映射规则，由计算值文件（由模型运行生成）、观测值文件（即真值文件，用户自行准备）和权重文件（用户自行准备），分别对应model_obs、real_obs和weight字段列。

   

4. 当MODFLOW6模型中设置观测井时，会在输出文件夹（算例文件夹\Pest.mf6）中自动生成efModFlow6.obs.head.csv文件，我们使用此文件作为计算值文件，映射model_obs字段。

   真值文件和权重文件格式与计算值文件完全相同，仅存在数值差异。参考计算值文件fModFlow6.obs.head.csv，在同一文件夹下，创建efModFlow6.real.head.csv文件，根据观测井水位填写不同模型网格中水位观测值，映射real_obs字段；创建efModFlow6.weight.head.csv文件，设置观测井权重，映射weight字段。三个文件内容如下图所示。

   三个文件准备完成后，在上一步的算例文件夹\PEST\obs.csv文件中输入字段列的映射规则（文件路径）。

   注：本算例中，已在算例文件中为用户准备好需要创建的efModFlow6.real.head.csv、efModFlow6.weight.head.csv和\obs.csv文件，用户可将前两个文件放置于算例文件夹\Pest.mf6文件夹中；将obs.csv文件放置于算例文件夹\PEST文件夹中，并根据算例文件夹\Pest.mf6文件夹位置修改obs.csv文件中的映射规则。

   

   

   

5. 在自动调参1图层数据子面板中继续设置模型参数。选择参数文件.csv所在位置，即保存水平渗透系数预设值、调参初始值、阈值、参数分组等信息的K参数.csv文件所在位置。分别将参数值、参数分组、初始值、最小值、最大值和拟合选项设为current、type、initial、min、max和parTrans数列。点击应用进行调参。

   

6. 若已完成一次调参，修改参数后再次调参，应用按钮为灰色时，可将视图切换至表格视图，勾选自动调参图层左侧方框使其可见，启动自动调参。

7. 通过一定时间计算，完成调参，在表格视图查看自动调参1或K参数.csv1图层， current属性列即为调参之后各点的水平渗透系数。

   

8. 在水头残差.efp工程中，右击efModFlow6.nam6图层，重新导入文件，查看水头残差变化。可以看到与调参之前相比，多数观测井中的残差明显变小。

   

小结

完成本指南中的算例后，可以学到以下内容：

• 如何创建并运行地下水MODFLOW6水流模型；
• 如何准备MODFLOW6模型参数文件；
• 如何计算并展示水头残差；
• 如何准备PEST调参文件；
• 如何设置PEST模型参数。