2简介SUB(thesubsidenceandaquifer-systemcompactionpackage)是由美国地质调查局开发的用来模拟含水层压缩和地面沉降的程序包。Sub程序包模拟包括弹性和非弹性夹层的压缩,也能够模拟夹层的滞后和非滞后排水。2.1概要下面是操作步骤:1、把sub程序包加入到网格法模拟模型中去2、创建一个简单的概念模型去说明sub程序模块是怎样被概化并且插入到modflow中去的3问题描述我们将要解决的问题将在图1中说明,模型基于美国地质调查局模型,描述如下:羚羊峡谷(antelope)是莫哈维沙漠(mojave)西部的地形封闭的盆地,距离洛杉矶大约50英里。羚羊峡谷地下水流域面积大约940平方英里,北部被低矮山丘和断层分开,1972年以前,地下水占总供水源的90%多;自1972年以来,供应量占50-90%。包括正在迅速的发展的城市lancaster和palmdale在内,在羚羊峡谷大规模抽取地下水。本区地下水系统由上层、中层和下层含水层组成。含水层由砾石、砂、淤泥和粘土冲积沉积和粘土及淤泥质粘土湖相沉积,为非固结沉积。研究区之前的地下水开发,补给来自周围山区的径流入渗。地下水从补给区流向沙漠盆地通过蒸发和泉流出露进行排泄,部分地下水流水平向阻隔,例如断层,已经在地下水流区域确定。地下水位下降由于地下水开发已经减少了自然因素的排泄,而农村和城市的抽水已经成为了该区地下水系统的主要排泄因素。回水入渗已经成为了该区地下水系统重要的补给来源。该模型被离散为网格模型,为43行,60列,3层。每层对应相应含水层。模拟期从1915年-1995年80年,并且第一年为稳定状态。4开始如果必要,打开gms。如果gms已经运行,选择File|New命令确保项目设置是默认状态5读取工程首先,我们读取工程:1、首先File|Open2、浏览所建工程文件夹3、打开start.gpr.你会看到图1的modflow模型6用网格法加入sub程序包6.1重命名并保存模型我们将要做一下变更,接下来是重命名并且保存模型1、选择file→saveas2、确保一直在工程文件夹中3、对该工程重命名为avgrid.gpr.4、保存6.2激活sub程序包我们需要打开sub程序包1、选择MODFLOW|GlobalOptions打开MODFLOWGlobal/BasicPackage对话框2、选择Packages按钮打开MODFLOWPackages对话框3、在Optionalpackages.下勾选SUB-Subsidence程序包4、确认两次退出对话框6.3定义非滞后夹层总结:定义夹层的前期固结水头、弹性与非弹性储水系数,其中前期固结水头为含水层中初始水头(也即是那个2d散点)首先我们先为第一和第二模型层加入非滞后夹层1、选择MODFLOW|OptionalPackages|SUB-Subsidence打开MODFLOWSUBPackage对话框这个步骤是打开SUBpackage对话框。对话框左面Options标签包括程序包常用值,滞后夹层区域值在右面。Delayinterbed和no-delayinterbed被用来创建夹层和设置夹层相关数据序列2、选择NoDelayInterbeds标签3、通过点击两次在No-delayinterbedlayers表格底部的按钮插入两个新的夹层系统,。4、为第第一夹层赋值1,第二层夹层赋值2.选择一个夹层系统,其相应的序列值就被显示出来。5、点击夹层系统1,它的值就显示在了序列编辑界面6、选择2Ddataset→array按钮设置前期固结水头(使用初始水头生成的网格中的最大值和最小值)7、选择preconsolidatedhead数据集然后选择ok按钮退出对话框8、在View/Edit弹出菜单项选择(sfe)elasticskeletalstoragecoeff(夹层的弹性储存系数)(本项之所以会出现弹性与非弹性的设置都是遍布整个网格的数值,是因为这些是设定值,意思是说,假如夹层发生了弹性形变就用sfe计算,发生非弹性就用sfv计算)。9、选择Constant-Array按钮,设置序列数值为1.5e-410、在View/Edit弹出菜单项选择(sfv)inelasticskeletalstoragecoef11、选择Constant-Array按钮设置,设置序列数值为8.0e-312、同步骤5,点击夹层系统2,它的值就显示在了序列编辑界面13、仿照步骤5-11进行设置,同样选择PreconsolidatedHead数据集,(Sfe)Elasticskeletalstoragecoef值为9.0e-5,(Sfv)Inelasticskeletalstoragecoef值为5.0e-3。6.4定义滞后夹层下面为第一二层定义滞后夹层1、选择DelayInterbeds标签2、通过点击两次在Delayinterbedlayers表格底部的按钮插入两个新的夹层系统,。3、输入以下序列值。用2DDataSet-Array按钮对Dstart(夹层初始水头)和DHC(初始预固结水头)进行赋值(DZ等效厚度,夹层号NZ)SystemlayerRNB(等效夹层数)Dstar(夹层初始水头)DHC(初始预固结水头)DZ(等效厚度)NZ(夹层号)111.0Starting1Preconsolidatedhead5.51221.0Starting2Preconsolidatedhead4.724、在对话框上部选择Options标签5、更改NumberofMaterialZones(NMZ)为26、输入以下滞后夹层的地层属性值IDVerticalK垂直渗透系数Elasticspec.storage弹性单位储水量inelasticspec.storage非弹性单位储水量11.0e-65.0e-66.0e-420.5e-65.0e-66.0e-46.5输出垂直位移下面开始垂向位移数据的生成,相当于第一层的沉降生成。这些垂向位移在MODFLOWsolution里将会作为数据集显示。1、选择SUBOutputOptions按钮,打开MODFLOWSUBPackageOutputOptions对话框。2、选择PopulateTimeSteps(填入时间步长)按钮3、更改弹出菜单项为Specifiedoutputlasttimestepeachstressperiod(额定输出每个应力期的时间步长)并且选择ok关闭对话框4、在本单元格里,在allrows.(所有行)勾选Savevert.disp(Ifl8)5、选择ok两次退出MODFLOWSUBPackageOutputOptions和theMODFLOWSUBPackage7运行modflow现在,我们保存文件然后与运行MODFLOW。1、选择保存按钮2、选择MODFLOW|RunModflow命令3、当MODFLOW运行完,选择close按钮4、选择保存按钮保存工程8检查结果现在我们将会看到更接近的计算结果。首先我们会看到水均衡整合到SUBpackage8.1水均衡1、在工程列表选择Head数据集,然后,在TimeStepWindow选择最新的时间步长2、选择MODFLOW|FlowBudget命令打开FlowBudget对话框SUBpackage水均衡值包括INST.IBSTORAGE和DELAYIBSTORAGE。其值显示如下表。我们将会使用这些值去对照后面的我们用概念法建立的同样的模型。TypeFlowInFlowOutINST.IBSTORAGE298,000-6,400DELAYIBSTORAGE150,000-1,2103、选择ok关闭对话框8.2查看垂直位移1、在工程列表里点击avgrid(MODFLOW)里面的VerticalDisplacement数据集2、如果有必要在TimeStep窗口里滚动鼠标并且点击最后的时间步长模型里面1486单元的垂向位移不同于附近的区域高达6.8ft3、在ProjectExplorer,在VerticalDisplacement和DrawDown之间来回变更数据集并且观察这两个数据集的等值线的相似性。8.3创建垂直位移图像下面我们生成一个图像来显示单个单元格的垂向位移。1、确保VerticalDisplacement被选中2、选择Grid|FindCell并且在CellID里输入1486。3、Ok退出对话框,此时这个单元格将会在GraphicsWindow被选中4、点击。或者选择Display|PlotWizard5、在PlotType列表里选择ActiveDataSetTimeSeries并且点击完成按钮生成的图像如图2所示。在VerticalDisplacement和DrawDown之间切换数据集展示这两个之间的关系6、保存工程9创建概念模型下面我们用概念模型法把同样的夹层系统加入到初始模型。9.1重命名并且保存模型1、选择File|New关掉网格模型2、打开start.gpr文件3、另存为avconc.gpr并保存9.2创建概念模型1、在工程列表空白处右键并且选择New|ConceptualModel2、在ConceptualModelProperties对话框,重命名为AnetelopeValley.3、把FlowPackage项变更为BCF.4、点击ok9.3创建第一个图层1、在AntelopeValley概念模型上右键并且在弹出菜单上选择NewCoverage2、在新图层对话框中重命名为layer1.3、在面域属性表,勾选如下选项:•SUBDelayInterbed•SUBNon-delayInterbed4、ok并退出9.4创建多边形1、选择ayer1图层使其处于激活状态2、选择创建弧3、创建一个围住网格模型的弧。如图3那个外框4、选择FeatureObjects|BuildPolygons。9.5设置第一层多边形的属性1、变换选择工具2、双击任意多边形3、在属性对话框,设置如下的值,其他的属性设为默认。SUBSfe(elast.skel.st.coef,ND)0.00015SUBSfv(inelast.skel.st.coef,ND)0.008SUBRNB(nequiv,D)1.0SUBDZ(bequivequiv.thick.,D)5.5SUBVerticalk(D)1.0e-006SUBElasticspec.storage(D)5.0e-006SUBInelasticspec.storage(D)0.00064、ok退出对话框9.6创建第二层1、在layer1图层右键在弹出的菜单上点击Duplicate2、重命名为layer2,把Defaultlayerrang中的1to1t变为2to2.3、ok退出对话框9.7设置第二层多边形的属性1、确保layer2图层是激活的2、变换选择工具3、双击任意多边形4、在属性对话框,设置如下的值,其他的属性设为默认。SUBSfe(elast.skel.st.coef,ND)0.00009SUBSfv(inelast.skel.st.coef,ND)0.005SUBRNB(nequiv,D)1.0SUBDZ(bequivequiv.thick.,D)4.7SUBVerticalk(D)5.0e-007SUBElasticspec.storage(D)5.0e-006SUBInelasticspec.storage(D)0.00065、Ok退出对话框10map→modflow概念模型被建立,所以我们可以把它绘制到MODFLOW网格里1、选择FeatureObjects|Map→MODFLOW2、Ok11sub程序包数组值让我们看一下这些在modflow里被绘制入概念模型里面的SUB程序包的数据1、选择MODFLOW|OptionalPackages|SUB–Subsidence命令。注意到两个地层已经在Delayinterbedmaterialzoneproperties(