建立MODFLOW模型有两种方法:网格方法或概念模型方法。网格方法是直接利用3D网格,应用surces/sinks和其它模型参数,基于单元的形式建立模型。概念模型方法利用Map模块的GIS工具来建立概念模型。在概念模型里,sources/sinks,含水层参数,模型边界等都可以在概念模型里设定。当模型建立完成后,将生成网格并把概念模型转化为网格模型,所有的单元属性都自动设定。
- 问题描述
本教程使用的模型是模拟垃圾堆放区的溶质运移问题。研究区域情况如下图所示:
模型北部为含水层露头,有两条河流汇聚,含水层分两层,在模型中处理为潜水含水层和承压含水层。
北部边界处理为0流量边界,其它边界处理为定水头边界 (据河流)。假设系统的补给源主要为降雨补给。区域内一些小溪有时会干涸,但由于受地下水的补给有时会有流水。将以排水的方式来处理这些小溪。另外在模型中还有两个生产井。
- 启动GMS
- 本例用到的模块
Grid,Geostatistics,Map,MODFLOW。
- 导入背景图
打开tutfiles\modfmap\start.gpr。此工程里存储了模型的底图。
- 保存工程
另存为easttex。
- 定义单位
点击Edit | Units命令,定义length为ft,time为d,其它不变,点击OK。
- 定义边界
- 新建层
- 切换到Map模块,右击Map Data文件夹,选择New Conceptual Model,取名为East Texas,模型选择MODFLOW,点击OK。
- 右击East Texas,选择New Coverage,取名为Boundary,设置Default elevation为700,Default layer range设为1到2,点击OK。
- 建立arc
点击Create Arc工具,以下图所示位置点击模型的西边。南部边界以河流为准,北部边界为基岩为准。双击到起点位置结束。(注:取消上一个点按backspace,取消整条线按ESC)
- 建立Local Source/Sink层:这一层定义区域边界,local sources/sinks(包括:井,河流,排水沟和常水头边界),依据概念模型,属性值可以赋给feature objects。在建立feature objects前,先在Coverage Setup对话框里更改属性。
右击Boundary层,选择Duplicate命令,把层的名字改为Sources/Sinks
右击Sourcds/sinks层,点击Coverage Setup
在Preset选择Sources/sink Coverage选项。这将打开一系列标准的source/sink类型的属性。
在Sources/Sinks/BCs清单里,关掉以下本教程不需要的项:specified flow,general head,river,seepage face,barrier.
确保选中了Use to define model boundary(active area)
点击OK。
- 建立定水头边界线
在南部和东部边界建立 定水头边界,在此之前,需要把先前建立的边界线分成3条线,一条在顶部将定义为0流量边界,其它两条将定义为两条河。通过在线上选择一个或多个顶点,将其转化为节点可分割线。
- 点击select vertices工具
- 选择下图所示的两个顶点。顶点1在两条河流的交汇处,顶点2在东部边界河流的出口处。先点击一个顶点,然后按住shift再点击第二个顶点可同时选中两个顶点。
- 点击feature objects | vertices <> nodes命令。
现在已经定义了三条线了,下面将定义两个河上的线为定水头边界:
- 点击select arc工具
- 选择南部和东部河流上的边界线(按shift同时选择两条)
- 点击properties
- 在all行选择spec.head类型
- 点击OK
- 点击模型其它部分取消取边界经的选择状态,此时河流的两条线将变颜色。
下面将在边界线的终点节点上定义水头。沿着定水头边界线上的水头假设为沿着边界线线性变化。
- 点击select points/nodes工具
- 双击南部边界线的西(左)部的节点
- 在head-stage里输入常数697
- 点击OK
- 以同样的方式,设置交汇点处为685,右边顶部节点为703.
- 定义排水边界线。在小溪河床处定义drains。
- 点击create arc工具
- 如下图所示建立三条弧线(从南部边界向北依次点):
现在把线定义为drains交设置它的conductance和高程
- 点击select arcs
- 选中所有的drains线
- 点击properties
- 在all行,设置type为drain
- 输入传导系数为6000.
- 为每条线把from layer和to layer都设为1,表示所有的drains都只在第1层。
- 点击OK
高程是定义到两个节点上的,其间的高程通过插值得到
- 点击select points/node
- 双击下图所示节点2,注意此节点有2个 属性,这是因为它附加到了2个不同类型的弧线上
- 在drain属性的bot.elev输入696,不要更改spec.head属性里的内容,点击OK
- 重复以上步骤,为其它节点 输入高程值。
- 建立polygons
在local sources/sinks层,模型的整个区域都必须由非重叠的多边形覆盖。这将定义网格的活动区域。大多数情况下,所有的多边形都是变水头多边形(默认),但是,其它的多边形可能用到,例如,模拟一个湖,将会用到常水头多边形。定义多边形最简单的方法是首先建立图层里用到的所有的线,然后点击build polygons命令。这个命令搜索所有的闭合的弧线来建立多边形。这些多边形默认为变水头,但可以通过选择多边形后在它的属性里更改为其它类型。现在所有的线都已经建立好,开始建立多边形。此模型所有的多边形都是变水头多边形。
- 点击feature object | build polygons
尽管建立了多边形,但在显示上没有任何显示。可以通过在display属性里设置打开polygons(fill)选项来查看多边形。
- 建立井
在local sources/sinks层最后需要定义的是井。井被定义 为点要素,这里将建立两个井。
- 点击create point
- 在图上大概位置处单击建立井,当新的点选中后,在GMS窗口的顶部输入X,Y坐标(2741,4673)。
- 点击属性,选择类型为well。
- 在flow(pumping)rate处输入定流量-24100
- 把from layer和to layer属性改为1,表示这个井只在第1层。点击OK
- 以同样的方式建立第二个井,坐标(10557,3290),流量-100000.但这个井的from layer和to layer改为2,表示这口井仅在第2层。
网格重定义:在井的附件,为了更精确的模拟水流,需要进行重新定义(通常是加密),可以通过GMS设置重定义数据来完成。
- 点击select points/node工具
- 选择两个井
- 点击properties按键
- 找到refine列,在all行,打开选项,更改base size为75,bias为1.1,max size为500
- 点击OK即可
- 描绘补给区
建立模型的下一是建立一个定义了补给区域的层。 假设在模型的补给是一致的,除了垃圾场以外。在垃圾场处的补给由于 垃圾掩埋系统的作用会减少。
- 复制边界
- 右击boundary层,点击duplicate,取名为recharge
- 右击recharge层,点击coverage setup
- 在areal properties清单里,打开recharge rate属性,点击OK
- 建立垃圾场边界
接下来,建立描述垃圾场边界的弧线。首先在垃圾场处建立一个矩形的闭合圈。然后编辑节点跟顶点,使定义的边界与垃圾场一致。
- 点击create arc工具
- 建立如下图所示的矩形多边形。此时不用考虑精确的坐标。
- 点击select vertices
- 围绕垃圾场拖出一个矩形框,以选择所有的顶点
- 点击feature objects | vertices<>nodes
- 点击select points/nodes
- 选择矩形框的其中一个顶点
- 此时可以在编辑窗口输入精确的坐标
- 重复以上步骤以设置其它点
- 建立多边形
- 点击feature objects|build polygons
- 设置补给值
对垃圾场多边形输入一个值,其它面输入另一个值
- 点击select polygons
- 双击垃圾场多边形
- 设置recharge rate为0.0002
- 点击OK
- 双击其它多边形,设置recharge rate为0.0228,点击OK
- 定义渗透系数
为每个层都输入渗透系数值。通常要定义多个多边形以输入不同的参数。为了简便,此教程为每个层输入一个常量
- 复制边界
- 右击boundary层,点击duplicate,取名为Layer1
- 右击Layer1层,点击coverage setup
- 在areal properties清单里,打开Horizontal K和Vertical anis属性,更改default layer range为1 to 1,点击OK
- 右击Layer1层,点击duplicate,取名为Layer2
- 右击Layer2层,点击coverage setup
- 更改default layer range为2 to 2,点击OK
- 顶层
输入顶层的K值
- 在data tree里选择layer1层
- 点击feature objects|build polygons
- 点击select polygons,双击多边形
- 把horizonal K改为16,vertical anis改为4,点击OK
- 底层
输入底层的K值
- 在data tree里选择layer2层
- 点击feature objects|build polygons
- 点击select polygons,双击多边形
- 把horizonal K改为32,vertical anis改为4,点击OK
- 此时就完成了定义概念模型的层,在 建立网格之前,先使sources/sinks层成为活动层,在data tree里选择sources/sinks层。
- 定位Grid Frame
现在层都已经完成,将建立网格。首先是确定网格的位置与方向,利用Grid Frame来做。
- 点击feature objects|grid frame
- 点击new frame
- 把Z坐标改为550,dimension改为200.这为MODFLOW层高设置一个初始值,稍后将插值层高。
- 点击OK退出Grid frame。
- 建立Grid
- 点击feature objects|map->3D grid。注意网格通过grid frame的数据限制了大小,如果grid frame不存在,网格默认为模型外边界5%处左右。注意在X,Y区域的单元数不能更改,这是由于行列数以及单元边界的位置都由在井处重定义的点数据来控制。
- 在Z-Dimension更改Number cells为2
- 点击OK
- 定义活动/非活动区域
现在已经生成了网格,下一步是定义 活动和非活动区域。这可以利用local sources/sinks层的信息来自动完成。
- 选择一个多边形
- 点击properties
- 确认layer assignment为1到2,点击OK
- 点击feature objects|activate cells in coverage(s)
- 初始MODFLOW数据
转化概念模型为基于网格的数值模型,在这之前,需要初始MODFLOW数据。
- 切换到3D Grid模块
- 点击MODFLOW|New simulation,点击OK
- 插值层高
在保存模型并运行MODFLOW之前的最后一步为定义层高和初始水头。由于使用的是LPF包,不管层的类型,都需要定义顶底板高程。两层模型需要定义第一层顶板高程、底板高程以及第二层底板高程。定义含水层高程的最简单方法为导入一系列定义高程的散点数据,并插值高程到含水层。有时,可以使用一系列散点数据来完成,有时,要使用两个系列的散点系列:一个用于定义地表,另一个用于定义第1层和第二层的底板。地表点通常通过地图数字化后得到,层高可以通过钻孔数据得到。
- 导入地表散点数据
散点数据已经被读取进来,因为它已经包括在工程文件之中。这些点来自于导入的text文件,在2D Geostatistices教程里有描述。这些散点是隐藏的,先要 取消隐藏。
- 切换到2D scatter point模块
- 在数据树,选中terrain和elevs两个系列
- 选中terrain使其成为活动数据
- 插值水头和高程
- 点击interpolation|interpolate->MODFLOW Layers。弹出设置窗口。
- 使ground_elev和Starting Heads为高亮状态,点击Map按钮
- 使ground_elev数据系列和Top Elevations layer1为高亮状态,点击map按钮
- 点击OK执行插值
- 插值含水层高程
- 选择elevs散点系列
- 点击interpolation|interpolate->MODFLOW Layers。此时GMS为自动的把bottom elevations layer 1和bottom elevations layer 2数据系列进行匹配。
- 点击OK
- 调整显示
- 在data tree里不选散点数据
- 切换到map模块
- 点击display options
- 关掉grid frame
- 点击OK
- 查看模型剖面
为了检查插值结果,我们查看剖面
- 切换到3D Grid模块
- 选择居中的某个单元
- 点击view j axis按钮
为了更好的查看剖面,重新设置Z轴比例
- 点击display | settings
- 把Z magnification设为5
- 点击OK
注意此时剖面的右边有尖灭部分,这将在运行模型的时候进行修正
- 修复高程数据
在Model checker有许多修复含水层数据的工具
- 点击MODFLOW|Check simulation
- 点击Run Check
- 点击Fix layer errors
可以看到在layer2许多错误提示,将用truncate to bedrock选项进行修正
- 选择Truncate to bedrock
- 点击fix affected layers
- 点击OK退出对话框,点击Done退出Model checker
注意层错误已经修正,查看含水层的修正方法是在plan视图里查看
- 点击view k axis切换到Plan view
- 在mini-grid显示里,点击向下按钮查看第二层,此时北部区域有非活动单元
- 点击向上按钮切换回第1层
- 转换概念模型
- 切换到map模块
- 点击feature objects | map -> modflow/modpath
- 选中all applicable coverages,点击OK
- 检查模拟
现在,已经完成定义MODFLOW的数据,准备运行模拟,先运行Model checker以检查是否还有错误存在。
- 切换到3D Grid模块
- 点击Modflow | check simulation
- 点击Run check,这里应该没有任何错误。
- 点击Done退出
- 保存工程
- 运行MODFLOW
- 点击MODFLOW | Run MODFLOW,此时MODFLOW启动了,并且Model Wrapper显示出来
- 当求解完成,点击close
- 显示水头等值线
当完成后会显示一系列的等值线,但为了更好的显示,可以更改线的颜色
- 点击data | contour options
- 在contour的line style按钮处点击向下按键
- 选择dark blue色
- 点击OK退出
要显示第二层的等值线
- 在mini-grid处点击向下按钮
- 在查看完等值线后,点击向上按钮返回到第一层
- 在Side View里查看潜水面
- 点击select cell
- 在模型的右边的井的附件选择一个单元
- 点击view j axis
- 点击view k axis
- 查看流量平衡
MODFLOW求解方案包括水头文件和单元的流量文件(CCF)。GMS可以利用CCF文件来显示流量平衡值。例如,我们可能想了解是否有水从排水沟里排出,可以通过点击drain弧线来实现。
- 切换到map
- 点击select arcs
- 点击最右边的drain arc
可以在窗口底部看到通过此线的流量,正面我们查看流到河流的水量
- 点击其中一个定水头边界线
- 按住shift时点击每一个定水头边界,此时可以看到选中边界线上的总流量。以此方式可以查看通过多边形的总流量,例如看补给多边形的补给量。要查看一系列单元的流量可以按以下步骤:
- 切换到3D grid
- 通过拖到一个矩形来选择一组单元
- 点击data | flow budget
这个对话框显示了选中单元详细的流量平衡数据
- 点击Done退出
- 点击模型的边界外以取消对单元格的选择
- 结论
- 可以导入背景图来辅助建立概念模型
- 定义一个模型边界层,然后再复制它成为别的层是一个好方法
- 可以使用coverage setup对话框来设置关联到点、线、面的属性
- 一些线的属性,例如水头,其给定值不是选择线来设置的,而是线的终端两个节点,其中沿着线中间的属性是通过线性插值得到
- 通过map -> modflow/modpath命令把概念模型转化为网格