Visual_modflow实例分析

实例：一个由两种不同地层构成的含水系统，北部和南部为无流量边界，东部和西部为河流，与地下水有充分的水力联系，可视为定水头边界，西边界与东边界的水头分别为9m，8m。

Fig. 8-1 Configuration of the sample problem

含水系统为非承压各向同性，第一层、第二层的水力传导系数分别为0.0001m/s和0.0005m/s，两层的纵向传导系数假设为横向的10%，有效孔隙率为25%，地面标高为10m。第一层、第二层的厚度分别为4m、6m。含水层接受面状补给量为8×10^－9mm。在第一层中靠近西部边界有一污染源。现在我们的任务是，如何用一个位于东部边界附近的完整井来隔离该污染源。

那么，在该区域内，需要建立一个数学模型来逐步计算抽水井所需的抽水量。只有当抽水量足够大时，被污染的地下水才能被抽水井全部截获。我们将使用PMWIN来建立数学模型，使用PMPATH来计算抽水井的截获带。基于所计算的地下水流场，使用MT3D和MOC3D来模拟污染源的运移情况。并演示如何使用PEST和UCODE来校正水流模型。最后创建一个动态演示方式，演示地下水被污染的发展变化过程。

为了示范运移模型，我们假想污染物溶解于地下水中的速度为1×10^-4μg/s/m²，含水层纵向弥散和橫向弥散系数分别为10m和1m，延迟系数为2，初始浓度、分子扩散系数及衰减速率均为0。我们将计算三年之后污染浓度分布情况，显示在两个含水层中点[x,y]=[290,310],[390,310]的浓度－时间突变曲线。

稳定流模拟

稳定流模拟中的六个主要步骤：

1. 建立模拟模型

2. 模型赋值

3. 执行水流模拟

4. 模拟结果检验

5. 子区水均衡计算

6. 输出结果。

步骤一：建立模型

1、打开File菜单，选择NewModel。

NewModel对话框弹出。选择一个保存模型数据的文件夹，如：C:\PM5DATA\SAMPLE，键入文件名：SAMPLE作为示例模型。模型的文件扩展名必须是“.PM5”。在Window95/98/NT下，文件名的有效使用字符数为120。最好在一个单独的文件夹中保存每一个模型及输出数据。可以同时运行多个模型（多任务处理）。

2、点击OK。

PMWIN只需几秒钟即可建立一个新模型，模型的文件名显示在标题栏中。

步骤二：模型赋值

包括：创建模型网格、定义边界条件、模型单元赋值。

在整个模型中，PMWIN要求单位使用统一。例如，若以m表示长度单位，s作为时间单位，水力传导系数必须以m/s表示，抽水量必须以m³/s表示，弥散系数必须以m表示。

在MODFLOW中，一个含水系统可以由一个包含一系列节点的离散区域或一些关联的有限差分块取代。图8-2，显示了一个包含网格和节点的含水系统的空间离散化情况，在每一节点均可计算出水头值。节点网格形成了数值模型的框架。一个地质单元可由一个或多个模型层表示。每一计算单元格的厚度、长度、宽度均可改变。计算单元所处的位置可以用列、行、层来表示，PMWIN使用索引符号[J,I,K]来指定计算单元的位置。如：位于第2列、第6行、第1层的计算单元记为[2,6,1]。

Fig. 8-2 Spatialdisretization of an aquifer system and the cell indices

一、 创建模型网格

1、打开Grid菜单，选择Mesh Size。

Model Dimension对话框弹出。

2、输入含水层的层数3，行数30、列数30，行高、列宽20。

第一含水层和第二含水层分别由一个模型层和两个模型层来代替。

3、点击OK。PMWIN的界面发生变化，此时显示了模型的整体网格。

PMWIN允许用户改变和旋转模型网格。用户可以随意改变模型网格行、列宽度，也可以增加/删除行和列。在该实例中，不需更改模型的网格。

4、从File菜单中选择Leave Editor或点击Leave Editor按钮退出。

二、定义含水层类型

1、打开Grid菜单，选择Layer Type。

Layer Options对话框弹出。

2、点击标签为Type列中的某一单元格，单元格中将出现一个带下箭头的按钮，点击该按钮，出现一个可供选择含水层类型的列表。

3、为第一层含水层选取1：Unconfined，其它含水层选取0：Confined，然后点击OK，关闭对话框。

三、定义边界条件

在模型中，边界上的各个计算单元均用不同的指定代码表示不同的边界类型：（1）“1”表示变水头边界，（2）“－1”表示定水头边界，（3）“0”表示无流量边界。

1、从Grid菜单中选择Boundry Condition→IBOUND(Modflow)。

PMWIN的Data Editor界面弹出，界面显示一个模型网格平面图，单元格指针位于[1,1,1]，即第一层的左上角。当前格的数值显示在状态栏底部。IBOUND默认值为1。单元格指针可以使用键盘的箭头键移动或使用鼠标点击你所需要的位置。使用PgUp、PgDn，或在工具栏的编辑框中输入含水层的层号后按Enter键，可以从一层移到另一层。

2、击鼠标右键，PMWIN显示一个Cell Value对话框。

3、在对话框中输入－1，点击OK。

计算单元[1,1,1]被定义为定水头计算单元。

4、点击Duplication按钮，打开复制开关。

如果Duplication键下陷，表示复制开关处于开启状态。移动网格指针，可将当前计算单元的数值复制到单元格指针经过的所有计算单元，再次点击Duplication键，关闭复制开关。

5、从左上角[1,1,1]移动单元格指针至左下角[1,30,1]，-1值被复制到模型西部边界的所有计算单元中。

6、移动单元格指针到右上角[30,1,1]。

7、再将单元格指针从右上角[30,1,1]移至右下角[30,30,1]，-1值被复制到模型东部边界的所有计算单元中。

8、点击Layer Copy按钮，打开层复制开关。

该按钮下陷，表示层复制开关处于开启状态。此时，从一层移至另一层，当前层的所有值将被复制到另一层中。再次点击该按钮，将关闭层复制开关。

9、按PgDn键，从第一层移至第二层、第三层。

第一层的值被复制到第二层、第三层。

10、点击Leave Editor按钮或从File菜单中选择Leave Editor退出定义边界条件界面。

四、模型几何参数赋值

（一）模型层的顶板高度赋值

1、打开Grid菜单，选择Top of Layer(Top)。

PMWIN显示模型网格。

2、打开Value菜单选择Reset Matrix或按Ctrl+R键。

Reset Matrix对话框出现。

3、在对话框中输入10，点击OK。

第一层的顶板高度定义为10。

4、按PgDn键移至第二层。

5、重复步骤2、3，输入第二层、第三层顶板高度6、3。

6、从File菜单中选择Leave Editor或点击Leave Editor按钮退出。

（二）模型层的底板高度赋值

1、从Grid菜单中选择Bottom of Layer(BOT)。

2、重复上述步骤，输入第一层、第二层、第三层的底板标高值6、3、0。

3、点击Leave Editor按钮退出。

五、时间及空间参数赋值

在本实例中，时间参数包括时间单位、应力期数、时段、水质运移时段。空间参数包括初始水头、水平及垂直方向的水力传导系数、有效孔隙率。

（一）时间参数赋值

1、打开Parameters菜单，选择Time。

Time Parameters对话框出现。在MODFLOW中，模拟时间被划分为几个应力期，也就是说，在外部应力作用下，时间间隔为常数，应力期又划分为时段。在大多数水质运移模型中，每一水量模型的时段又被分成为更小的运移时段，时段的长度并不与稳定流模型对应。然而，当我们使用MT3D和MOC3D来运行污染运移模拟模型时，实际的时段长度需在表中给定。

2、输入第一应力期的长度9.46728E+07（秒）。

3、点击OK接受其它默认值。

（二）初始水头赋值

1、打开Parameters菜单，选择Initial Hydrolic Head。

PMWIN显示模型网格。

2、从Value菜单中选择Reset Matrix（或按Ctrl+R键），在对话框中输入8，点击OK。

3、将网格指针移至左上角[1,1,1]。

4、点鼠标右键，在对话框Cell Value中输入9。

5、点击Duplication按钮，打开复制开关。

6、从左上角向下拖动鼠标，当前计算单元中的值9将被复制到鼠标所经过的所有计算单元中。

7、点击Layer Copy 按钮，打开层复制开关。

8、按PgDn键，移至第二层、第三层。

9、点击Leave Editor按钮，退出。

（三）导水系数赋值

1、从Parameters菜单中选择Vertical Hydraulic Conductivity。

2、打开Value菜单，选择Reset Matrix（或按Ctrl+R），在对话框中输入0.00001，点击OK。

3、按PgDn键，移至第二层。

4、重复步骤2、3，分别在第二层、第三层中输入0.00005。

5、点Leave Editor按钮退出。

（四）有效孔隙率赋值

1、从Parameters菜单中选Effective Porocity。

输入有效孔隙率值，并保存。模型的默认值为0.25。

2、退出。

（五）面状补给量的赋值

1、打开Models菜单，选择MODFLOW→Recharge。

2、从Value菜单中选择Reset Matrix…，在对话框的Recharge Flux[L/T]编辑栏中输入8E-9。

3、退出。

（六）定义抽水井位置及抽水量赋值

在MODFLOW中，一个注水井或抽水井由一个节点（或计算单元）取代，用户可用一个节点定义一个注水井或一个抽水井，并且假设井穿透整个含水层。MODFLOW可以模拟穿透多个含水层的抽水井、或抽水井在每一层中有不同的抽水量。多层井的总抽水量等于各单层抽水量之和。每一层的抽水量可以根据下式计算：

T_k：含水层的导水系数；∑T：被井穿透所有含水层的导水系数之和。但是，当第一层含水层为非承压时，如果不假设一个饱水厚度，我们无法确切知道井所在位置含水层的饱水厚度和导水系数，上述公式无法应用。模拟一个多层井的另一种可能是：对井所在的每一个计算单元设一个很大的垂向导水系数（或垂向渗漏），如1m/s，总抽水量赋予井的最底部计算单元。为了演示，在井中的其它计算单元中赋一个非常小的抽水量（1×10^-10m³/s），这样，来自每个被穿透含水层的水量将被MODFLOW计算出，并且该值可通过Water Budget Calculator来获得。

由于我们不能确定截获所有污染物所需的抽水量，我们假设一抽水量0.0012m³/s来进行调试。

1. 打开Models菜单，选择MODFLOW→Well。

2. 将网格指针移到[25,15,1]。

3. 点鼠标右键，输入-1E-10，点击OK。负值表示抽水井。

4. 按PgDn键，移到第二层。

5. 点鼠标右键，输入-1E-10，点击OK。

6. 按PgDn键，移到第三层。

7. 点鼠标右键，输入-0.0012，点击OK。

8. 退出。

步骤三：执行水流模拟

1、打开Models菜单，选择MODFLOW→Run.。Run Modflow对话框弹出。

2、点击OK，开始水流模拟。

在运行MODFLOW之前，PMWIN将利用用户自定义的数据为MODFLOW（及可选MODPATH）产生一个输入文件，该文件列在Run Modflow对话框的列表中。如果Generate标记为 ，表示只产生了一个输入文件。可以点击该按钮锁定Generate标记。一般情况下，不需更改该标记。

步骤四：模型结果检验

在水流模型模拟中，MODFLOW在列表文件path\OUTPUT.DAT中有一个详细运行记录。如果成功地完成了一个水流模拟，MODFLOW将以多种非格式化文件（二进制）保存模拟结果。如下表。在运行MODFLOW之前，用户可以控制这些非格式化（二进制）文件的输出。具体操作：从Models菜单中选择Modflow→Output Control。若Interbed Storage Package被激活，则产生path\INTERBED.DAT输出文件。

为了检验模拟结果的质量，MODFLOW在每个时段结尾为整个模型作了水量均衡计算，并保存在OUTPUT.DAT文件中。一个水均衡提供了一个数学解的所有可解性迹象。

MODFLOW的输出文件

文件	内容
Path\OUTPUT.DAT	详细运行记录和模拟报告
Path\HEADS.DAT	水头
Path\DDOWN.DAT	降深、初始水头与计算水头之差
Path\BUGET.DAT	各计算单元的流量项
Path\INTERBET.DAT	每一层中整个含水层的沉降、压密及压密之前的水位
Path\MT3D.FLO	与MT3D/MT3DMS的界面文件，这个文件由MT3D/MT3DMS提供的软件包LKMT创建。

Path为保存模型数据的文件夹。

步骤五：子区水均衡计算

为了计算方便，每个计算单元的流量项被保存在文件path\BUDGET.DAT中，这些单一计算单元的流量项实际上是一个单元流向一单元的流量，有四种类型：（1）一个单元向一个单元的应力流量；或受外部应力影响，流入或流出一个单元的流量，如抽水井或面状补给；（2）计算单元间的储存项：在各计算单元中，储存量的增加和减少的速率。（3）计算单元间的定水头流量项：流向或来自各个定水头计算单元的净流量。（4）计算单元间内部的流量项：通过各计算单元截面的流量，也即相邻计算单元间的流量。Water Budget Calculator利用逐个计算单元的流量来计算整个模型的水量均衡、子区及子区间的水量均衡。

子区水均衡计算

1．打开Tool菜单，选择Water Budget。

Water Budget对话框弹出。在稳定流计算中，不需改变Time的设置。

2．点击Zones。

一个Zone代表一个水均衡计算子区，由0－50间的一个代号（区号）表示，每个区的区号必须赋予每一个计算单元。区号0表示该计算单元与其它区无关。

3． Value菜单打开ResetMatrix…对话框，输入1，点击OK。

4． PgDn键，移到第二层。

5． 从Value菜单中，选择Reset Matrix…，输入2，点击OK。

6． 点击退出按钮，退出。

7． 点击OK，退出Water Budget对话框。

PMWIN在文件path\WATERBDG.DAT中计算并保存流量项。流量单位为[L³/T]。对每一时段、每一层的每一个区分别计算。流入区内的流量记为IN，子区之间的流量记在Flow Matrix中，通过子区边界的水平方向流量由HORIZEXCHANGE给出。EXCHANGE（UPPER）给出了来自（IN）或流向（OUT）上部相邻层的流量，EXCHANGE（LOWER）给出了来自（IN）或流向（OUT）下部相邻层的流量，如LAYER＝1，ZONE＝1的EXCHANGE（LOWER）（表示第一层、第一区与下部相邻层的流量交替项），从第一层到第二层为流量为2.587256E-03m³/s，计算结果中的PERCENT DISCREPANCY可由下式计算出：

在该实例中，整个模型中每层、每区的流入、流出量的百分比误差很小，这意味着模型方程有正确的解。

为了计算流向抽水井的确切流量，重复上述步骤计算水量均衡。这次，只需对计算单元[25,15,1]定义为1区，计算单元[25,15,2]定义为2区，计算单元[25,15,3]定义为3区，所有其它计算单元均定义为0区，水量计算结果列在表2.4中。抽水井在第一层中的抽水量为7.7992809E-05m³/s，第二层为5.603538E-04m³/s，第三层为5.5766129E-04m³/s，几乎所有抽水量均来自第二层，与实际含水层结构吻合。

步骤六：结果输出

除水量均衡计算之外，PMWIN提供了多种模拟结果检验、创建图形输出的可能。流线及速度矢量可以通过PMPATH显示。利用ResultExtractor，可将任一时段任一层的模拟结果从非格式化（二进制）结果文件中读取，或保存在ASCII Matrix文件中。一个ASCIIMatrix文件含有一个模型层的计算值（每计算单元一个值）。PMWIN能将ASCII Matrix文件装入模型网格中。以下给出操作步骤：

1． 用Results Extractor读取或保存计算的水头。

2． 于第一层的计算水头生成等水位线图。

3．使用PMPATH计算水流运行轨迹及抽水井的截获带。

（一）读取及保存计算的水头

1、从Tool菜单中选取Results Extractor。

ResultsExtractor对话框弹出。对话框中含有四个选项标签，分别为MODFLOW、MOC3D、MT3D、MT3DMS。在MODFLOW工作表中，可从Result Type下拉列表中选择结果类型，也可以指定要读取数据结果的层、应力期、时间步长。该电子表格显示了一系列的行和列，一行与一列的交汇处为一个单元格。其中的每一单元格与模型层中计算单元一一对应。下列步骤2、6以三个ASCIIMatrix文件为每一层保存了水头值。

2、点击Result Type的下拉箭头按钮，选择Hydraulic Head。

3、在Layer编辑框中输入1。

该实例中（稳定流模拟只有一个应力期一个时段）应力期数及时段数均为1。

4、点击Read按钮。

第一层、第一应力期，时段1的水头值将被读入电子表中。点击滚动条，可游览整体结果。

5、点击Save按钮。

Save Matrix对话框出现。在Save as Type选项中，可根据用户需要将数据结果保存为ASCII格式或SUFFER格式。给定文件名H1.DAT，保存文件，点击OK。

6、重复步骤3、4、5，分别保存第二层、第三层的水头计算结果H2.DAT，H3.DAT。

7、点击Close，关闭对话框。

（二）生成计算水位等值线图

1、从Tools菜单中选取Presentation。

在Presentation中定义的数据并不能被PMWIN使用，但可以用Presentation来临时保存数据或显示模拟结果的图形。

2、从Value菜单中选择Matrix…（Ctrl+B）。

Brows Matrix对话框弹出。点击Load按钮，可将ASCIIMtrix文件装入表中，或点击Save按钮，将表中数据保存为一个ASCIIMatrix文件。还可通过菜单Value中的Result Extractor读取水位数据，或使用Apply按钮在Presentation数据表中添加数据。

3、点击Load按钮。

Load Matrix对话框出现。

4、点击Open按钮，选取由Result Extractor保存下的文件H1.DAT，点击OK，H1.DAT被装入。

5、点击OK，关闭Browse Matrix对话框。

6、打开Option菜单，选择Environment（Ctrl+E）。

Environment Option对话框出现。对话框中包括三个选项列表。Appearance和 Coordinate System用来修改模型网格的外观及位置，Contours用来生成等值线图。

7、点击Contour标签，选中Visible复称选框，点击Restore Default按钮。

点Restore Default按钮，PMWIN默认等值线数为11，并使用当前层中的最大和最小值定义等值线的最大、最小值。如果FillContours被选取，点击列标签Fill，可自选颜色填充等值线图。用Label Format，可定义合适的数据格式。

注意：当退出编辑时，PMWIN将清除Visible复选框。

8、点击OK，退出Environment Options对话框。

9、从File菜单中选择Save Plot As…或Print Plot，保存或打印图形文件。

10、按PaDn移至第二层，重复步骤2－9，装入文件H2.DAT，显示和保存图形文件。

11、点击Leave Editor按钮，然后点击Yes，退出并保存Presentation结果。

通过以上步骤，可以用输入的数据、各种模拟结果、或以ASCIIMatrix格式保存的数据生成等值线图。如：可利用初始水位生成等水位线图，或使用ResultExtractor读取污染浓度分布数据显示等值线图，也可用Field Interpolator或Generator生成实测等值线图。

（三）绘制抽水井截获带

1、打开Models菜单，选择PMPATH(Pathlines and Contours)。

PMWIN调用水平运移模型PMPATH，并将当前模型自动装入PWPATH中，PMWIN使用“网格指针”来指定所显示平面图的行和列。按住Ctrl键并点击鼠标左键，可将网格指针移动到你所需要的位置。

注意：如果你随后在PMWIN中修改或计算了一个模型，你必须再次装入所修改的模型，以确保所做的修改被PMWIN确认。点击Open aModel按钮，从Open model对话框中选取带有扩展名“.PM5”的模型文件，装入文件。

2、计算抽水井的截获带

（1）点击SetParticle按钮。

（2）移动鼠标指针到模型区域，鼠标指针变为十字线。

（3）将十字线移到抽水井左上角。

（4）按住鼠标左键向右下方拖拽十字线，拉出一个小窗口将整个抽水井包括进来。

（5）释放鼠标左键。

弹出Add New Particles对话框。在编辑框中输入质点数，点击Properties标签，然后点击颜色按钮，为新质点选择合适的颜色，点击OK。

（6）按PgDn键，移到下一层，重复（3）、（4）、（5），为第二层、第三层设置抽水井周围的质点，可选用不同颜色表示。

（7）点击＜按钮，开始追溯质点的运移轨迹。

PMPATH计算并显示流程线的方向及抽水井的截获带。

为了详细地观看平面流程线的运移情况，可从Options菜单中选择Environment…，打开一个Enviroment Options对话框，在Cross Setions表中的Exaggeration编辑框中输入一个放大倍数值。剖面图中，一些流线的末端上行至地表，表示地下水接受补给，这是三维流模型与二维流模型的主要不同之处。在二维流模型中，由于不存在垂向流量（或始终为0），所以，在有降水补给的地区，也不会有流线翘向地表。

PMWIN还可生成抽水井截获带的动态演示图。具体操作：打开Option菜单，选择Particle Tracking(Time)，弹出Particle Tracking(Time)Properties对话框，然后，点击＜按钮。从这里可以看出，由于第一层比其它层的传导系数小，第一层的截获带也比其它层更小。