MODFLOW-SFR地表地下水交互模拟
- 项目图示
- SFR图层管道
- MODFLOW图层管道
算例简介
本算例研究区域地图如下图所示，区域内有多条区段（Segment）组成的河网，在对该区域地下水进行数值模拟时，既需要考虑地表水与地下水之间的补排交互过程，也需要逐段考察河道中的地表水如何向下游汇集。前者与River程序包类似，地表水与地下水之间的补排关系由它们的水头值决定；而后者需要考虑地表河网的水量分配、河流流经过程中的地表径流、降水补给与蒸发、河道形状与河床岩性、地表水与地下水交互等因素的影响。因此研究复杂河网和地下的交互作用需要使用MODFLOW-SFR程序包。
在SFR程序包中，需要创建河流区段（Segment）来代表河流的空间位置，根据河流的流向和水量分配，从上游向下游为每个区段和河段依次进行编号，构建河网拓扑结构。随后使用模型网格切割区段，在每个单独的网格中得到河段（Reach）。
本算例将首先创建单层网格的MODFLOW水流模型，再根据地表河网的分布，构建河网拓扑结构，为区段赋值河床参数、区段水力参数等属性，完成SFR河网程序包创建。将带有SFR程序包信息的图层链接入已有的MODFLOW水流模型中，解释运行模型展示地表与地下水交互模拟得到的地下水水位分布。
创建MODFLOW水流模型
新建项目
- 打开Envifusion软件，若软件界面仍有其他项目，点击菜单栏
按钮关闭当前项目文件。
创建MODFLOW模型
首先，创建用于运行MODFLOW2005模拟的网格。
点击创建数据体-> MODFLOW->创建MODFLOW2005模型，创建流场模拟（MODFLOW2005）图层。将模型网格划分设置为(10，15，1)，即XYZ方向的网格数分别为10、15和1。模型X边界分别设置为(0，50000)，Y边界分别设置为(0，75000)，Z边界分别设置为(0，1)。
勾选水平渗透系数选项，并将默认水平渗透系数设为10；勾选非稳定流选项，并将默认给水度和默认贮水系数设为0.15和1e-6，点击应用完成设置。
用户应自行确认模型单位的一致性，单位会在后续操作中进行设置。本算例长度单位为英尺，时间单位为秒。此时系统将产生初始MODFLOW模拟图层。
定义活动区域
导入场地地图，根据地图设置模型活动区域。
点击界面上方的
按钮，打开并导入算例文件夹中的basemap2.jpg文件，选择jpge读取器，点击应用生成basemap2.jpg图层。选中basemap2.jpg图层，点击通用通用工具箱->空间拓扑->空间变换生成空间变换1图层，在该图层中属性面板数据子面板中，将X和Y方向的缩放设为36.36和41，点击应用完成设置。使底图和创建的模型网格大小一致。
点击视图窗口上方的俯视视角：Z轴向外
选项，将视图变为俯视图。点击创建数据体->几何体->多边形，生成多边形1图层，将其重命名为iBound图层。在该图层属性面板数据子面板中，勾选显示手柄、封闭折线选项，将下方的添加类型设置为多边形。此时根据basemap2.jpg图层中的研究区域边界创建弧段。
按住键盘上的Alt+Ctrl键，依次点击鼠标左键创建弧段。创建弧段前，数据子面板弧点坐标栏中有两个原始点（0，0，0）（1，0，0）。可以选中原始弧点，按住鼠标左键拖动弧点改变弧点位置；也可创建完所需弧段后，在坐标栏选中原始弧点，点击
按钮删除原始弧点。弧段创建完毕在属性面板中点击应用创建表示研究区域的多边形。在创建弧段过程中，为保证最终生成的多边形可以覆盖活动区域，尽量在Z轴向外的视图下创建弧段，必要时可以点击视图窗口上方的3D
按钮，将视图先转换为2D视图。
在图层面板中选中iBound图层，点击通用工具箱->属性->计算器，生成计算器1图层。在该图层属性面板数据子面板中，将结果数列名称设为iBound，将数值设为0，点击应用完成设置。
选择流场模拟（MODFLOW2005）图层，点击数值模拟-> MODFLOW2005-> MF2005多边形赋值，选择计算器1图层为属性图层，点击确认生成MF2005多边形赋值1图层，将其改名为010iBound- MF2005多边形赋值图层。
在010iBound- MF2005多边形赋值图层属性面板中将属性类别设为有效网格，层号为1，iBound选择iBound数列，勾选下方的反转选项；在渲染子面板中将渲染方式设为曲面网络，色彩设置为efIBound，点击应用完成设置。视图窗口如下图所示：红色部分为进行数值模拟的活动区域。
定义网格高程
导入顶底板高程文件，定义网格高程。
点击界面上方的
钮，打开并导入算例文件夹中的top.csv文件，点击应用生成top.csv图层。选中top.csv图层，点击通用工具箱->散点->表格->散点生成表格->散点1图层，在该图层中属性面板数据子面板中，为X、Y和Z分别选择x、y和z数列，点击应用完成设置。
选择010iBound-多边形赋值图层，点击数值模拟-> MODFLOW2005-> MF2005网格高程，选择表格->散点1图层为属性图层，点击确认生成MF2005网格高程1图层，将其改名为020 MF2005网格高程1图层。
在020 MF2005网格高程1图层属性面板中将赋值对象设为模型顶板，勾选下方的上层优先选项，点击应用为网格赋值顶部高程。
重复步骤1-3，导入bottom.csv文件，并将其转换为散点数据，生成表格->散点2图层。使用该图层为020网格高程1图层赋值底部高程，使用模拟-> MODFLOW2005->网格高程工具，生成 030网格高程2图层。
在030网格高程2图层属性面板数据子面板中将赋值对象设为单层底板，层号为1，取消勾选下方的上层优先选项；在渲染子面板中将渲染方式设为曲面，点击应用完成设置。网格如下图所示：
定义渗透系数
导入包含渗透系数分布的表格，定义模型渗透系数。
点击界面上方的
按钮，打开并导入算例文件夹中的K.csv文件，点击应用生成K.csv图层。选中K.csv图层，点击通用工具箱->散点->表格->散点生成表格->散点3图层，在该图层中属性面板数据子面板中，为X、Y分别选择x、y数列，勾选下方的2D散点选项，点击应用完成设置。
选择030网格高程2图层，点击数值模拟-> MODFLOW2005->MF2005层赋值，选择表格->散点3图层为属性图层，点击确认生成MF2005层赋值1图层，将其改名为040水平渗透系数- MF2005层赋值1图层。
在040水平渗透系数- MF2005层赋值1图层属性面板数据子面板中将属性类别设为水平渗透系数，层号设置为1，水平渗透系数选择K数列；在渲染子面板中，将渲染方式设为曲面，色彩设为ef水平渗透系数，点击应用完成设置。模型水平渗透系数分布如下图所示：
定义应力期
导入包含应力期信息的表格，定义模型应力期。
点击界面上方的
按钮，打开并导入算例文件夹中的stress.csv文件，点击应用生成stress.csv图层。选择040水平渗透系数-层赋值1图层，点击数值模拟-> MODFLOW2005->应力期，选择stress.csv图层为属性图层，点击确认生成应力期1图层，将其改名为050应力期图层。
在050应力期图层属性面板数据子面板中将应力期长度选择Length数列，当期步数选择Tstep数列，稳定流选择Steady数列，点击应用完成应力期设置。
定义给水度、垂向补给、蒸散与初始水头
分别导入包含模型给水度、垂向补给、蒸散和初始水头信息的表格，定义模型参数。
重复2.5节中的操作，从算例文件夹中分别导入包含模型给水度、垂向补给、蒸散和初始水头信息的表格：sy.csv、recharge.csv、et.csv和startH.csv。点击应用生成sy.csv、recharge.csv、et.csv和startH.csv图层。
选中sy.csv图层，点击通用工具箱->散点->表格->散点生成表格->散点4图层，在该图层中属性面板数据子面板中，为X、Y分别选择x、y数列，勾选下方的2D散点选项，点击应用完成设置。
选择050应力期图层，点击数值模拟-> MODFLOW2005->MF2005层赋值，选择表格->散点4图层为属性图层，点击确认生成MF2005层赋值1图层，将其改名为060给水度- MF2005层赋值图层。
在060给水度- MF2005层赋值图层属性面板数据子面板中将属性类别设为给水度（潜水），层号设置为1，给水度（潜水）选择sy数列，点击应用完成设置。
操作相同步骤，生成070垂向入渗- MF2005层赋值图层、080蒸散- MF2005层赋值图层、090初始水头- MF2005层赋值图层。各图层属性面板数据子面板设置如下图所示。
在各图层属性面板渲染子面板中，将渲染方式设为曲面，色彩设为不同图层对应的属性，点击应用可查看不同属性在网格中的分布。
点击开始->保存工程，将上述结果保存为flow.efp文件。
创建SFR模块
定义河网弧段
导入场地地图，根据河网分布绘制弧段。
点击EnviFusion图标，新建EnviFusion工程文件，将新工程保存为SFR.efp文件。
点击击界面上方的
按钮，打开并导入算例文件夹中的basemap2.jpg文件，点击应用生成basemap2.jpg图层。选中basemap2.jpg图层，点击通用工具箱->空间拓扑->空间变换生成空间变换1图层，在该图层中属性面板数据子面板中，将X和Y方向的缩放设为36.36和41，点击应用完成设置。使底图和创建的模型网格大小一致。
点击视图窗口上方的俯视视角 :Z轴向外
选项，将视图变为俯视图。点击创建数据体->几何体->多边形，生成多边形1图层，将其重命名为A图层。在该图层属性面板数据子面板中，勾选显示手柄选项，将下方的添加类型设置为折线。此时根据basemap2.jpg图层中的左上角的区段（Segment）位置创建弧段。具体创建方法可参考2.3节步骤3，创建完成后点击应用完成设置。
在图层面板中选中A图层，点击通用工具箱->图层->添加多边形，生成添加多边形1图层，将其重命名为B图层。在该图层属性面板数据子面板中勾选显示手柄选项，按住键盘Alt+Ctrl键，按照下图所示创建B区段弧段，弧段绘制完成后点击应用完成设置。
重复步骤5操作，依次创建C、D、E、F、G和H图层。创建好的弧段方向为终止点指向起始点，与上图中方向相反，需要在下一步中进行操作调整区段方向。
创建SFR模块
创建SFR模块，定义河网的拓扑结构和渗流参数。
点击界面上方的
按钮，打开并导入算例文件夹中的基本参数.csv文件，点击应用生成基本参数.csv图层。点击视图窗口上方的新建按钮，点击表格，创建视图2表格窗口，在图层面板中选中基本参数.csv图层，查看表格数据。
选择H图层，点击数值模拟-> MODFLOW-SFR>创建SFR2模块，选择H图层为一维河流图层，选择基本参数.csv图层为河流参数图层，点击确认生成创建SFR2模块1图层，将其改名为010创建SFR2模块图层。
在010创建SFR2模块图层属性面板数据子面板中，分别为各参数选择相对应的数列，河网的区段ID是区段唯一标识；上下游ID定义了河网的拓扑关系；水流方向是一列0/1数列，用来逆转河水流动方向；引水工程也是一列0/1数列，用来标识从上游河道截取部分河水的引水工程；ICALC参数定义了SFR模块中使用河网河段中心流量求取河水深度的方法，有0-4共五种选项。
点击界面上方的
按钮，打开并导入算例文件夹中的河床参数.csv文件，点击应用生成河床参数.csv图层。该图层包含上下游河床的渗透系数、河床厚度和河床高程信息。选择010创建SFR2模块图层，点击数值模拟-> MODFLOW-SFR->河床渗流参数，选择010创建SFR2图层为一维河流图层，选择河床参数.csv图层为河流参数图层，点击确认生成河床渗流参数1图层，将其改名为020河床渗流参数图层。
在该图层属性面板数据子面板中，分别为各参数选择相对应的数列，点击应用完成设置。
定义河段水力参数
由于不同区段（Segment）根据流量求取河流深度的方法不同，ICALC取值不同。0代表人工指定河流深度；1代表使用曼宁公式（矩形渠道）求取河流深度；2代表使用曼宁公式（自定义剖面）求取河流深度；3代表使用经验公式（指数）求取河流深度；4代表使用水文过程线求取河流深度。不同求取方法需要输入的区段水力参数不同，为不同ICALC区段输入水力参数。
点击界面上方的
按钮，打开并导入算例文件夹中的ICALC0.csv文件，点击应用生成ICALC0.csv图层。选择020河床渗流参数图层，点击数值模拟-> MODFLOW-SFR->区段水力参数，选择020河床渗流参数图层为一维河流图层，选择ICALC0.csv图层为水力参数图层，点击确认生成区段水力参数1图层，将其改名为030区段水力参数0图层。
在该图层属性面板数据子面板中，选择参数资源文件保存路径，并将其保存为para.sfr2.efs；应力期号设为1-2；ICALC设置为0：人工指定；分别为各参数选择相对应的数列，点击应用完成设置。
重复步骤1-2，分别导入ICALC1.csv、ICALC2.csv、ICALC3.csv、ICALC4.csv文件，生成ICALC1.csv图层、ICALC2.csv图层、ICALC3.csv图层、ICALC4.csv图层。
基于030区段水力参数0图层，点击数值模拟-> MODFLOW-SFR->区段水力参数，生成 040区段水力参数1图层、050区段水力参数2图层、060区段水力参数3图层、070区段水力参数4图层。
在040区段水力参数1图层属性面板数据子面板中，分别为各参数选择相对应的数列，点击应用完成设置。
在050区段水力参数2图层属性面板数据子面板中，分别为各参数选择相对应的数列。此外还需要在附加参数中为每一条ICALC值为2的区段定义剖面的形状，每个剖面需要8个点进行定义。此处有两条区段的ICALC值为2，所以共需要输入16行数据。因为剖面涉及两个维度：距离和高程，因此仅需要为附加参数中的第一列和第二列赋值，系统将忽略第三列的输入。两组剖面的节点输入如下所示：
| 区段G | 区段H |  | ||
|---|---|---|---|---|
| 距离 | 高程 | 距离 | 高程 | |
| 0 | 20 | 0 | 25 | |
| 10 | 13 | 10 | 17 | |
| 80 | 10 | 80 | 13 | |
| 100 | 2 | 100 | 4 | |
| 150 | 0 | 150 | 0 | |
| 170 | 10 | 170 | 10 | |
| 240 | 13 | 240 | 16 | |
| 250 | 20 | 250 | 20 |
在060区段水力参数3图层属性面板数据子面板中，分别为各参数选择相对应的数列，点击应用完成设置
在070区段水力参数4图层属性面板数据子面板中，分别为各参数选择相对应的数列。还需要在附加参数中定义水文过程线信息，附加参数的1-3列分别为水文过程线的流量、深度和宽度值。
| 流量 | 深度 | 宽度 |  |
|---|---|---|---|
| 0.5 | 0.25 | 3 | |
| 1 | 0.4 | 3.5 | |
| 2 | 0.55 | 4.2 | |
| 4 | 0.7 | 5.3 | |
| 7 | 0.8 | 7 | |
| 10 | 0.9 | 8.5 | |
| 20 | 1.1 | 12 | |
| 30 | 1.25 | 14 | |
| 50 | 1.4 | 17 | |
| 75 | 1.7 | 20 | |
| 100 | 2.6 | 22 |
在该图层属性面板渲染子面板中，将渲染方式设为曲面，色彩设置为ef区段入口；在风格属性中，将线宽设为5，并勾选下方的将线渲染为圆柱选项，点击应用查看河网中区段分布。
- 在图层面板中选择070区段水力参数4图层，点击开始->导出图层，或点击保存
按钮保存图层，将该图层保存为sfrLayer.efd文件，在保存选项对话框中不勾选保存所有时间步选项，点击确定保存图层。
解释运行MODFLOW
SFR链接MODFLOW
导入包含SFR模块的图层，链接到MODFLOW程序中。
打开flow.efp工程文件，点击界面上方的
按钮，打开并导入上一步中生成的sfrLayer.efd文件，点击应用生成sfrLayer.efd图层。在图层面板中选中sfrLayer.efd图层，点击数值模拟-> MODFLOW-SFR->链接MODFLOW，选择sfrLayer.efd图层为一维河流图层，选择090初始水头-层赋值图层为模型工程图层，点击确认生成链接MODFLOW1图层，将其改名为100链接MODFLOW图层。
在该图层属性面板数据子面板中，选择SFR2模块文件保存路径，并将其保存为model.sfr，层位设为1，流量单位设置为立方英尺每天，平滑倍数设为12，点击应用完成设置。
解释运行MODFLOW
解释运行创建好的水流模型。
选中100链接MODFLOW图层，点击数值模拟-> MODFLOW2005->MF2005解释运行，生成解释运行MODFLOW1图层。在属性面板中，点击
按钮选择模型文件生成的位置，将模型名保存为flow。选择解算器为GMG，时间单位为秒，长度单位为英尺，保存应力期为1，点击应用，系统将自动编译并运行MODFLOW模型。
点击界面上方的
按钮，选择上一步中用户指定名称的MODFLOW文件夹中的efModFlow2005.nam文件。图层面板中生成efModFlow2005.nam图层，点击应用导入MODFLOW模型结果。其中hed文件为水头文件，lmt6文件为MODFLOW生成的链接文件，其中记载了边界信息、参数信息、以及逐网格的水力梯度信息，cbc文件为水均衡文件。
由于非活动区域水位值较低，需要对该区域进行裁切。选中efModFlow2005.nam图层，点击通用工具箱->空间拓扑->等维裁切，生成等维裁切1图层。在该图层属性面板数据子面板中，将切块类型设为数据系列，将数据系列设为ef水头，裁切值设为0；在渲染子面板中，将渲染方式设为曲面，色彩设为ef水头，点击下方的重置颜色阈值到数值边界
按钮，点击应用完成设置。模拟区水位如下图所示
在图层面板中打开sfrLayer.efd图层、100链接MODFLOW图层和等维裁切1图层显示。将100链接MODFLOW图层属性面板渲染子面板中的渲染方式改为网络，在坐标变换属性中将Z方向缩放设为10；将等维裁切1图层属性面板渲染子面板中的渲染方式改为矢量流线，在坐标变换属性中将Z方向缩放设为10，地下水流线如下图所示：
小结
完成本算例后，可以学到以下内容：
- 如何创建并运行地下水MODFLOW水流模型；
- 如何设置应力期和模型单位；
- 如何基于研究区底图创建用于SFR的河网拓扑结构；
- 如何为区段（Segment）赋值河流参数；
- 如何为区段（Segment）赋值区段水力参数；
- 如何将SFR模块链接MODFLOW并解释运行。