沉井施工降水时地下水动态分析及对江堤影响的研究

张文杰 沉井施工降水时地下水动态分析及对江堤影响的研究 63 1 绪论1.1 地下水模拟的重要性地下水的动态分析实际上就是地下水系统的数值模拟。地下水是一种宝贵的自然资源，是维持人类社会可持续发展的经济和战略资源，是综合国力的有机组成部分1。当今世界所面临的“人口，资源，环境”三大问题，都直接或间接地与地下水有关。地下水含水层具有较大的调蓄能力，且水量稳定，不易受污染，水质较好，在河流、湖泊被大量污染的今天，地下水在我国城市生活和工农业用水中起着举足轻重的作用，我国约有三分之二的城市以地下水作为重要的供水水源。随着人口的增长、工农业生产及城市建设的迅速发展，地下水资源开发利用的规模越来越大，由此导致的诸如水质恶化、形成区域性的降落漏斗、地面沉降等一系列地质灾害和环境问题，给人类带来了不同程度的危害，因此查清地下水资源的形成条件、进行正确的总量评价、对研究区地下水系统进行模拟分析,并结合其时空分布特点制定合理的开采利用方案十分必要2。地下水系统3指在一定的水文地质条件下，在某一范围内形成的地下含水系统，水力联系密切并与相邻含水系统相对隔绝。广义上的地下水系统指包括地下含水系统和与之相关的社会、经济、环境要素的总体。实际的水文地质条件是复杂的：地形起伏多变,含水层厚度不均一，非均质，各向异性，结构复杂，多层不连续，源汇项的时空分布不均，水文地质参数的时空变异性明显等。数值模型对这些特征有相对合理的表征，是对真实地下水系统的仿真和模拟。地下水资源评价就是“对某系统地下水资源的数量质量时空分布特征和开发利用条件作出科学全面的分析和估计”45，通常地下水资源评价都是指地下水资源允许开采量的评价。1.2 地下水模拟的研究现状1.2.1 国外研究现状最早对地下水进行动态分析，采用的是比较直观也是最为简单的水均衡方法以及水文地质比拟方法 6。十九世纪中叶，达西(Henri Darcy)在总结前人实践的基础上，通过试验提出了水在孔隙介质中渗透的线性渗透定律，即达西定律；稍后，袭布衣(J.Dupuit)以达西定律为基础，研究了单向和平面径向稳定运动，奠定了地下水稳定流理论的基础。1905年，E.梅勒第一次用解析法论证了泉水流量的预测方法 7。1935年，泰斯(C.V.Theis)提出了地下水流向承压水井的非稳定流公式，开创了现代水文地质计算的历史。五十年代，随着深层承压水的开发利用，雅柯布(C.E.Jacob)、汉土什(M.S.Hautush)等人研究了有越流补给的情况，接着出现了考虑无压含水层迟后反应、非完整井等情况下的解析解。同时把稳定流计算中已经行之有效的叠加原理、映射法应用到非稳定流计算中来，以解决井群干扰和边界的影响以及抽水流量呈阶梯式变化等非稳定流动问题。五十年代后期，T.H.卡门斯基在解析法分析群孔潜水动态的基础上，系统地研究了存在降水入渗条件下的有限差分法，并用它来预测地下水动态的变化。六十年代后期随着计算技术的进步，数值模拟方法应用到地下水计算中来，理论和工程上分析地下水的能力都取得了突破性的进展，先后出现了二维流平面(剖面)模型、准三维流模型、三维流模型、耦合模型等 8-12。国外近年来出现的几个主要模型，例如R.Bravo(美国)等做的美国休斯敦模型，美国维吉尼亚州的滨海平原下水三维模型，确定了地下水的运动模式，主要的补给和排泄区域以及无约束系统中地下水的循环次数，澳大利亚Queensland大学系统研究了地下水动力学的数值解法和建立模型的方法步骤，在实践中模拟了Collie流域地下水及其与地表水的动态耦合问题；并模拟了海水入侵使海岸含水层污染物迁移的过程，将数值模型应用到含水量计算、开采预测、管理方案及含水层补救等，取得了显著成效。目前地下水系统数值模拟方法13主要有有限差分法(FDM)、有限单元法(FEM)、边界元法(BEM)和有限分析法(FAM)等。国外该领域的研究主要针对数值模拟法的薄弱环节，提出新的思维方法，采用新的数学工具，分析不同尺度下的变化情况，合理地描述地下水系统中大量的不确定性和模糊因素。根据Anderson等14提出的工作程序，要建立一个正确且有意义的地下水系统数值模型，应进行以下工作：确定模型目标，建立水文地质概念模型，建立数学模型，模型设计及模型求解，模型校正，校正灵敏度分析，模型验证和预报，预报灵敏度分析，模型设计与模型结果的给出，模型后续检查以及模型的再设计。针对Newman等的推测，Wood15提出了二维地下水运动有限元计算的时间步长条件。Kim16等对抽取地下水造成的noordbergum effect(reverse water level fluctuation)现象进行数值模拟，阐述了其机理性原因。Scheibe17等分析了在不同尺度下的地下水流及其运移行为。Ghassemi18指出三维模型可以详细说明含水层系统的三维边界条件以及抽水应力情况，而二维模型就不能恰当处理。Li Shu-guang等19指出数值模型还不能解决预报的不确定性因素问题，并开创性地提出一种随机地下水模型，可以解决均值分布和小尺度过程的不同尺度问题。Mehl20等提出二维局部网格细分法的有限差分地下水模型，提供了新的插值和误差分析的方法，模拟结果的可靠性得到了提高。1.2.2 国内研究现状我国地下水数值模拟起步较晚，开始于20世纪70年代，随着地下水数值模拟研究方面的发展，地下水数值模拟已成为研究地下水运动规律和定量评价地下水资源的主要手段，而且发展趋势已远远超出作为一种计算手段的原有范畴，成为模拟水文地质过程、发生、发展的演变规律的主要方法。现在我国已经建立了囊括国际地下水模拟中心(IGWMC) P.Vender Heijde分类中所有模型，即预报模型(包括水流模型、热量运移模型、形变模型、多目标模型)、管理模型和识别模型。研究范围涉及到饱和带、非饱和带与饱和-非饱和带，基本满足了我国国民经济建设发展的需要。随着非稳定流理论的发展，以及电子计算机的广泛应用，使得各种复杂条件下的地下水运动都可应用数值法求解。近年来，随着新技术、新方法的广泛应用，我国该领域科学家做了大量的工作，在建立地下水系统数值模拟模型中发现问题，在理论和方法上不断创新，通过数值模型理论与相关研究方向的理论结合，不断提高模拟结果的可靠性。卢文喜21对地下水运动数值模拟中的边界条件进行了分析，提出在模型预报前要考虑自然因素、人类活动因素及邻区水流条件因素产生的耦合效应问题，先对边界条件进行预报。薛禹群等22介绍了Ms-FEM(多尺度有限元法)的基本原理，并将其应用于非均质多孔介质中的流动问题，通过计算结果的比较得出多尺度有限元法比传统有限元法有效的结论。南京大学地球科学系23用Galerkin有限元法，建立了裂隙-岩溶泉的柳林泉平面二维区域地下水渗流模型，描述山西柳林泉域的地下水渗流，模拟效果良好，据此预报了柳林电厂水源地投入使用后对区域地下水流场的影响，以及对柳林泉的影响。合肥工业大学24基于徐州市张集水源地裂隙岩溶水贮存和运动特征，建立了三维有限元数学模型，利用等参有限元方法求解，并根据大型群孔干扰抽水试验资料，识别模型参数，进一步对裂隙岩溶水资源的开采方案进行了评价，为水源地开采提供了依据，同时对条件类似地区的地下水资源评价具有重要的参考价值。中国地质科学院水文地质环境地质研究所25对黄河三角洲利津地区，建立三维地下水数值模型，分析预测了开采浅层地下淡水引起的地下水动力场的变化和咸水入侵的可能性。中国地质大学26在分析水文地质条件和建立黄河下游(河南段)影响带数学模型的基础上，采用FEFLOW软件建立了三维地下水流模拟模型，详细描述了模型中截渗墙的处理，并用模型模拟了截渗墙建成后黄河侧渗量和地下水流场的变化。1.3 课题的提出、研究内容、方法及要达到的目的1.3.1 课题的提出本课题研究的是结合在建的泰州长江大桥特大型陆上沉井施工降水方案的研究而开展的。泰州长江大桥是世界上首座三塔两跨千米级悬索桥。泰州长江大桥东距江阴长江大桥57公里，西距润扬大桥66公里，是江苏省“五纵九横五联”高速公路网的重要组成部分。大桥工程全长62.088公里，由北接线、跨江主桥、夹江桥和南接线四部分组成。工程起于宁通高速公路宣堡枢纽，在泰州永安洲镇跨越长江，向西于扬中小泡沙跨越夹江，经镇江姚桥镇进入常州境内，止于沪宁高速公路汤庄枢纽。跨江主桥为主跨1080m2的三塔双跨特大跨径钢悬索桥，主桥及夹江桥全长9.726公里，桥面宽33米。主桥通航净空高度不小于50米，净宽不小于760米，能满足5万吨级巴拿马货轮的通航需要。夹江桥通航净高不小于18米，净宽不小于100米。1.3.2 研究内容南北锚碇是悬索桥的重要组成部分，锚碇结构采用沉井。根据实际的施工情况，本文只对北锚碇的施工降水方案进行分析研究。北锚碇设计尺寸分别为52.4m68.3m，基底标高-55.0m27。锚碇的剖面图如图1-1。图1-1 锚碇的剖面图针对沉井施工方案，需要在沉井施工过程中降水。根据地质资料，锚碇区全为松散沉积物覆盖，它由第四系和上第三系构成，土层富水性好，且与长江联系密切。沉井降水深度较深，会引起沉井周边地面沉陷、开裂，严重威胁周边已有建筑物的安全与稳定，尤其会对附近长江大堤的安全产生重大影响。江堤属一级堤防，安全事关重大，因此必须对沉井施工降水后江堤的土体变形，渗透变形，整体稳定性进行分析评估28。在地面沉降及江堤变形预测和控制研究中，都需要对地层的空间展布形态、结构特点、物理力学性质、地下水分布与抽取状况进行综合研究。因而，利用当代计算机信息化技术，建立地层可视化模型，直观再现锚碇区不同性质(地质成因，展布方式，物质成分，结构特征，物理力学性质，水力特征）地层的空间展布形态和地层间的组合关系，在本课题研究中是非常必要的。1.3.3 研究方法 目前最常用的数值模拟方法是有限单元法和有限差分法13。根据不同的方法，国外编制了不同的数值模拟软件，主要有：加拿大Waterloo水文地质公司的Visual MODFLOW、德国WASY水资源规划和系统研究所的有限单元软件FEFLOW、美国Brigham Young大学环境模型系统试验室的GMS等。29这些国际流行的地下水流模拟软件有其共同的特点，即：功能齐全的输入、输出数据手段；友好的人机界面，可视化效果好。本课题采用有限元软件GMS(Groundwater Modeling System)对泰州长江公路大桥陆上沉井施工降水的地下水流路径进行三维模拟计算，分析江堤的渗流稳定，江堤及地面的沉降变形。GMS 30是美国Brigham Young University 的环境模型研究实验室和美国军队排水工程试验工作站在综合MODFLOW、FEMWATER、MT3DMS、RT3D、SEAM3D、MODPATH、SEEP2D、NUFT、UTCHEM 等已有地下水模型的基础上开发的一个综合性的、用于地下水模拟的图形界面软件。其图形界面由下拉菜单、编辑条、常用模块、工具栏、快捷键和帮助条6 部分组成,使用起来非常便捷。目前，国内应用这个软件分析地下水运行系统并于实际工程的实现还很少见。以下是对本文所需的GMS模块的简介。MODFLOW31是美国地质调查局于80年代开发出的一套专门用于孔隙介质中地下水流动的三维有限差分数值模拟软件。MODFLOW 自从问世以来,由于其程序结构的模块化、离散方法的简单化和求解方法的多样化等优点32-34 ,已被广泛用来模拟井流、河流、排泄、蒸发和补给对非均质和复杂边界条件的水流系统的影响。MODPATH是确定给定时间内稳定或非稳定流中质点运移路径的三维示踪模型。它和MODFLOW一起使用,根据MODFLOW 计算出来的流场,MODPATH可以追踪一系列虚拟的粒子来模拟从用户指定地点溢出污染物的运动。这种追溯跟踪方法可以用来描述给定时间内井的截获区。SEEP2D是用来计算坝堤剖面渗漏的二维有限元稳定流模型。它可以用于模拟承压和无压流问题,也可以模拟饱和和非饱和带的水流;对无压流问题,模型可以只局限于饱和带。根据SEEP2D的结果可以作出完整的流网。PEST和UCODE是用于自动调参的两个模块。在自动进行参数估计时,交替运用PEST 或UCODE来调整选定的参数,并且重复用于MODFLOW、FEMWATER等的计算,直到计算结果和野外观测值相吻合。MAP可使用户快速地建立概念模型。在MAP模块下,以TIFF、JEPG等图件为底图,在图上确定表示源汇项、边界、含水层不同参数区域的点、曲线、多边形的空间位置,快速建立起概念模型。TINs即三角不规则网络(Triangulated Irregular Net2works),通常用来表示相邻地层的界面,多个TINs就可以被用来建立实体(Solid)模型或三维网格。实体是在不规则的三角形网络(TIN)建立完成后,通过一系列操作产生的实际地层的三维立体模型。Borehole Data(钻孔数据)用来管理样品和地层这两种格式的钻孔数据。样品数据用来作等值面和等值线，地层数据用来建立TIN、实体和三维有限元网格。1.3.4 技术路线地下水数值模型的建立，需要考虑地下水含水系统中含水层、弱透水层与隔水层在空间的组合分布，以及水文地质参数的时空变化性。根据本次研究的目标和要求，制定技术路线如图1-2所示。图1-2 论文技术路线1.3.5达到的目的本次模拟建立三维可视化地层，动态直观的显示水流运动状态和路径，为降水施工方案提供依据；地层的不均匀性导致了地面的不均匀沉降，模型分析过程中能预测周边地区和江堤的沉降情况，有利于控制地面沉降的发生和发展；江堤的渗流模拟将显示江堤剖面完整的流网及等势线，能比较准确地分析江堤的渗透稳定性。1.4 本章小结本章主要介绍了地下水模拟的重要性及国内外的研究现状，并根据研究区的实际情况提出了研究方法、技术路线、目的及论文的创新点。3 地下水系统数值模拟3.1 地下水系统概化模型对地下水水流进行数值模拟,首先要确定研究区范围。研究区范围是客观存在的,而客观存在实质上是一个连续的整体。只有准确认识研究区的内因(地质结构与水文地质条件)和外因(边界条件)后,才能将其作为一个模拟研究的整体。数值模拟的目的就是要再现客观存在,研究其具有的客观规律36。而实际的水文地质条件是十分复杂的，要想完善地建立描述计算区地下水系统的数值模型是困难的。因此，应根据水文地质条件和地下水资源评价的目的，对实际的水文地质条件进行简化。这一过程称为水文地质条件的概化，其原则为：根据评价的目的和要求，所概化的水文地质概念模型应反映地下水系统的主要特征；概念模型要简单明了；概念模型要能够被用于进一步的定量描述，以便于建立描述符合研究区地下水运动规律的微分方程的定解问题。水文地质条件的概化通常包含以下几个方面：计算区域几何形状的概化；含水性质的概化；边界性质的概化；参数性质(均质或非均质,各向同性或各向异性)的概化；地下水流状态(一维、二维或三维)的概化37。3.1.1 地下水系统含水层结构特征35 （1）第一含水层（潜水）含水层属全新统，岩性为灰黄、灰色亚粘土、粉砂、底部为淤泥质亚粘土，层厚为6.68.3m。透水性一般，富水性相对较差，主要接受大气降水与田间沟塘水补给，根据抽水资料，在抽水期间下层水位没有明显变化，说明第一层地下水与下层地下水没有明显的水力联系，越流影响不明显。（2）第二含水层（承压水）含水层属全新统，主要岩性为粉砂，局部地段为细砂，顶板标高4.015.76m，底板标高23.2123.66m，层厚17.90m19.20m。砂质纯，透水性较好，富水性强，通过两组抽水资料对比可以看出：靠近长江方向的试验井涌水量丰富。（3）第三含水层（相对弱含水层）含水层属全新统，主要岩性为亚粘土混粉砂或夹粉砂，局部地段为粉砂，顶板标高23.2123.66m，底板标高49.4651.91m，层厚27.80m28.70m。单井涌水量较小，受长江潮汐的影响不明显，抽水期间降落漏斗影响半径小，观测井水位降深较小，因此该含水层透水性较弱。3.1.2 地下水系统水量交换特征从空间上看，地下水流整体上以水平运动为主、垂向运动为辅，地下水系统符合质量守恒定律和能量守恒定律。含水层分布广、厚度大，在常温常压下地下水运动符合达西定律。考虑浅、深层之间的流量交换以及软件的特点，地下水运动可概化成空间三维流。地下水系统的垂向运动主要是层间的越流，地下水三维模型可以很好的解决越流问题。地下水系统的输入、输出随时间、空间变化，故地下水为非稳定流；参数随空间变化，体现了系统的非均质性，根据抽水试验资料，含水层有明显的方向性，所以含水层可以概化成非均质、各向异性。所以，根据研究区水文地质条件，通过研究区地下水补给和动态变化特征，将本研究区概化成非均质、各向异性、空间三维结构、非稳定地下水流系统。3.1.3 地下水系统边界条件边界条件的研究是建立水文地质立体概念模型的一个重要的方面。边界可分为内、外边界。外边界是指各含水层组的周围边界，内边界是指沟通多层含水层之间的研究区内边界。根据该区的水文地质条件，研究基本上构成了地下水补、径、排的完整系统。其主要边界有两个：一是侧向边界，南侧沿江一带，接受长江水补给强烈，在模型中可设其为一类已知水头边界，水位值为长江边界附近的水位；东侧和北侧地下水的等水位线基本上与边界垂直，所以将东侧和北侧概化为零通量边界(流线边界)。二是垂向边界，潜水含水层自由水面为系统的上边界，通过该边界，潜水与系统外发生垂向水量交换，如接受长江补给、河塘补给、大气降水入渗补给、蒸发排泄等，因研究区施工时间较短，未考虑大气降水入渗补给、蒸发排泄。3.2 地下水系统数学模型3.2.1 数学模型的表达要确定一个水文地质问题的数学模型，需要从概念模型出发，用一组数学关系式来刻画它的数量关系和空间形式，以反映所研究地质体的地质、水文地质条件和所要描述的地下水流动的基本特征，达到再现或复制一个实际地下水系统基本状态的目的。这一组数学关系式就是通常所说的数学模型3839。这一组数学关系式不仅要概括体系的已知部分(水头等),还容纳了未知部分(参数)。通过已知的数据来预测体系的属性,即数学模型的解。在大多数情况下,地下水的动态运移由三维偏微分方程描述。这些方程,附加边界条件和初始条件后就构成了数学模型36。对承压含水和非承压含水的水流运动，在不考虑水的密度变化的条件下，以水头作为基本变量，多孔介质中地下水流动的有限差分公式数学模型可以用下面的偏微分方程384041来表示： （3-1） （3-2） （3-3） （3-4）其中：渗透系数在方向上的分量，我们假定渗透系数的主轴方向与坐标轴的方向一致；水头；单位体积流量，代表流进汇或来自源的水源；孔隙介质的贮水率；时间；研究区各层初始水头值；研究区各层第一类边界上的实测水头值；研究区各层第二类边界上的单位面积流量；一类及二类边界；承压含水层边界上的内法线。公式（3-1）加上相应的初始条件和边界条件，便构成了一个描述地下水流动体系的数学模型，该数学模型的解就是一个描述水头值分布的代数表达式，在所定义的空间和时间范围内，所求得的水头h应满足边界条件和初始条件。3.2.2 数学模型的求解数学模型建立后，其求解包括正演问题和反演问题。如果给定了含水层的水文地质参数、边界条件和初始条件，求解含水层各个时刻的水头值，进行水头预报，这类问题称为正问题，即正演；如果根据动态观测资料或抽水试验资料反过来确定水文地质参数，这个过程就叫反演，即前者的逆问题或反求参数问题。为验证所建立的数学模型的正确性，求解时一般要先进行模型识别，即反演过程。通过反演验证模型后，再利用反演所求得的水文地质参数进行正演，求解所需要的项目，如水位预报、涌水量预测等。上述建立的数学模型很难采用解析法来进行求解，一般都要借助计算机采用数值法求解。本文运用GMS软件进行数值模拟，其数值方法采用的是有限差分法。有限差分法的基本思想4243是将研究区域用差分网格离散，对每个结点通过用差商近似控制方程中的微商，然后耦合初始条件及边界条件求封闭的线性代数方程组。在有限差分法求解过程中，连续的时间和空间被划分成为一系列离散的点，在这些点上，连续的偏导数也由水头差分公式来取代，将所有的未知点联合起来，这些有限差分式构成了一个线性方程组，然后对这个线性方程组进行联立求解，这样获得的解就是水头在各个离散点上的近似解。数值解虽然不能给出描述水头随时间和空间变化的代数表达式，但它可以用来解决大量的实际问题。3.3 地下水系统数值模型3.3.1 模型前期处理3.3.1.1 网格划分数值模拟中，研究区将由结点和有关的有限差分网格或有限个单元组成的离散域所代替。结点距的选取是网格设计中的一个关键步骤，在水头变化大的地方或水头分布需要了解得比较详细的部位，结点距要小一点，但相邻结点距或单元的大小不能相差悬殊，要逐渐过渡。观测井和抽水井应尽可能地放在结点上，这样提供的水头、流量比较精确，进行校核也比较方便 38。本研究区计算面积为181350，采用GMS软件进行自动三角形网格剖分，剖分单元13058，结点9338。由于在河流和抽水井处水头变化比较大，因此对网格进行了加密。网格划分见图3-3。3.3.1.2 时间离散时间步长的选择和空间步长的选择一样是模型设计的一个关键步骤。在抽水试验时，水头下降速度不均匀，经历着一个由迅速下降至逐渐缓慢的过程，因此，抽水初期用较小的时间步长，以后随着水头下降速度的减缓，逐步加大时间步长。时间步长的增大一般按下式进行：，de Marsily44认为最好的选择是应用,本文所用的时间步长的增大系数为。3.3.1.3 水文地质参数和初值的确定水文地质参数是表征含水介质储水能力、释水能力和地下水运动能力的指标。分析、研究、求取水文地质参数，是地下水资源计算与评价的重要环节之一，参数选取的准确与否，直接影响地下水资源计算与评价的精度。水文地质参数常通过野外试验、实验室测试及根据地下水动态观测资料采用有关理论公式计算求取，或采取数值法反演求参等。本次模拟采用GMS中的LPF（Layer Propery Flow）子程序包，需要输入每一含水层的水平向渗透系数和垂直向渗透系数。在非稳定流计算中还需要输入贮水率和给水度。414546 水平和垂直向渗透系数的初步确定采用研究区已有的钻孔地质资料和抽水试验资料得到的实测参数值，然后基于这些离散点上的值采用特定方法的空间区域化方法形成研究区的初始参数值，模型的最终参数值由模型调试确定。（1）贮水率弹性贮存或弹性释放是含水层的重要水文地质特征。为了从数量上评价，需要表明含水层中贮存水量的变化和承压水头相应变化之间的关系，通常采用贮水率或贮水系数来表明这种关系。贮水率表示当含水层水头降低一个单位时，由于含水层内骨架的压缩和水的膨胀而从单位体积含水层中释放出来的总水量。当承压含水层的厚度固定不变时，有，它所释放的水量完全是弹性水量。潜水贮水系数所释放的水量包括两部分：一部分是含水层由于压力变化所释放的弹性水量；另一部分是水头变化一个单位时所疏干含水层的重力水量，这一部分水量正好等于含水层的给水度，由于潜水含水层的弹性变形很小，可以用给水度近似代替贮水系数。贮水系数是没有量纲的参数，其确定方法是根据非稳定流抽水试验资料，利用配线法（标准曲线对比法）的降深时间配线法推求41。由式： （3-5） （3-6）式中：抽水井流量； 抽水井井壁降深；抽水井半径；抽水时间； 含水层的贮水系数；导水系数；井函数； 对式（3-5）、（3-6）取对数，得 （3-7） （3-8）从解析几何中知道，如果有，那么曲线和的形状是相同的，只是曲线相对于曲线在横坐标上位移了，在纵坐标上位移了。因此，如果在双对数纸上做的关系曲线，在另一张模数相同的双对数纸上做关系曲线，显然，这两条曲线的形状也是相同的，只是坐标位移了一个距离而已。降深时间配线法的方法如下：根据数值表，在双对数纸上作一张的关系曲线，称之为标准曲线。在另一张模数相同的透明双对数纸上，根据抽水资料作关系曲线（即降深时间关系曲线），称之为实际资料曲线。把实际资料曲线重叠在标准曲线上，保持坐标轴平行，移动到二根曲线完全重合为止。在重合的双对数纸上任选一点作为配合点，读出改点在两张对数纸上相应的坐标，和的值，再把值代入公式（3-5）、（3-6）中，即可解出和来。具体计算见表3-1。表3-1 贮水率计算表井号配比点(h)(m)S1-110.08331.358.11.776.25.082.9310-31.57710-420.2782.53034.3632.1113.920055.30S1-210.08330.281.50.42.8820.2220.7412.8831.6672292.93.02S4-110.20.44100224.761.3820.60.930041.35341.1200051.38S4-210.11670.0830.42.5920.19440.21512.4731.390.423022.41S5-110.10.24100062.940.0954.110-41.45110-5220.32100040.953100.56500070.95S5-2110.2850050.158220.29100050.1523100.34500060.156（2）给水度给水度是表征潜水含水层给水能力或储水能力的一个指标4647。它与岩性有关，随排水时间、潜水埋深、水位变化幅度及水质的变化而变化。本次模拟的给水度根据研究区含水层岩性及颗粒组成，取给水度为0.06。（3）初始条件初始条件为某一个选定的初始时刻渗流区内的水头分布。初始时刻一般理解为按照需要任意选定的某一瞬间，但不一定是实际开始抽水时刻。本文模拟的初始水头是利用这一时刻各观测孔的实测水头用插值的方法得到的。在动态稳定流条件下，流入系统的水流等于流出的水流，水头在空间上随地而异。Anderson和Woessner认为“选用识别过的模型所得稳态水头解作为初始条件是一种标准做法”。48因此，在具体模拟的时候选择了一个水位相对稳定的时期作为初始时刻，输入边界条件，求出模型的稳态解，这个稳态解即可作为下一步非稳定流模拟的初始条件。在建立拟合期初始流场时采用现场抽水试验期2006年12月13日8:002006年12月17日20:00的各观测水位，用线性插值的方式来求出同一层面上各单元中心点的初始水位值，然后将模拟初始时刻水位值赋给模型。在模型预报中，选择了一个水位相对稳定的时期作为初始时刻，输入边界条件，运用所建模型求它在很短时间内的稳态解作为下一步非稳定流模拟的初始条件。3.3.1.4 源汇项的处理1.补给项本研究区地下水的侧向补给主要来源于河流渗漏补给和研究区四周边界补给，垂向补给主要为含水层越流补给，调用了GMS中的河流子程序包。2.排泄项地下水的排泄项主要以抽水井的形式来表示，调用了GMS中的水井子程序包，在模型中用结点表示，并给以一定的流量，负值表示抽水。3.3.2 建立三维可视化地层模型信息系统针对研究区含水层的埋藏深度、厚度、岩性，构建GMS地下水数值模拟模型。建模过程如下:（1）根据地质钻孔资料，由钻孔坐标确定每一钻孔在计算区域的坐标位置,生成Borehole柱状图。（2）将地面和江堤的实际位置和尺寸导入模型。（3）根据TINS模块，用Borehole中的Inverse Distance Weighted 差值方法生成Solids地层，见图3-1、3-2。（4）利用网格法将实体模型转化为网格模型。该步实现了有限差分方程中的计算系数转换。网格模型指的是地下水的建模过程和模型识别验证都是直接在三维网格上操作。网格法是一种最简单、最基本的建模方法,对于水文地质条件比较简单、研究区面积比较小的模型,网格法是一种有效的方法。网格法建模是基于对网格的直接操作,因此参数分区的赋值、边界条件的指定和调参都非常繁琐,其工作量相当可观而且枯燥容易出错，在本次建模过程中采用网格法与概念模型联合建模。概念模型建模的思想就是把三维网格的属性用coverage图层相应实体的属性表示,模型赋值和调参直接操作coverage上的实体。一个coverage可以包含点(point,代表曲线的顶点；node,代表曲线的端点/结点)、弧(arc,代表线性实体,如河流、沟渠)和多边形(polygon,代表面状实体,如渗透系数分区)。这些实体都有唯一标识号,能赋予属性49。本次模型的计算网格模型见图3-3。（5）在三维网格下设置Modflow的初始条件, 给定井的抽水量、河流、边界水位变化情况,选择计算包及迭代计算方法。在MAP模块下,将MAP文件中的数据转化到Modflow中,该步实现了有限差分方程中的源汇项的转换。(6)检查所建模型及定解条件的给定是否正确, 不断调整、改正后运行Modflow进行计算。图3-1 研究区实体模型图3-2 研究区实体地层剖面图3-3 研究区网格模型3.3.3 模型识别与检验通过对本研究区水文地质的合理概化，建立了水文地质的概念模型，根据这个物理的概念模型又建立与之对应的地下水渗流数学模型，经过这一系列的高度概化和对应模型的转换，尽管每个阶段都有理论依据，但转换后的最终模型是否能全面、客观地表征研究区的具体水文地质条件和特征呢？这就需要进行模型识别验证。模型的识别与检验过程是整个模拟中极为重要的一步工作50。模型的识别是为了证明所建立的模型确实有再现现场实测水头和水流场的能力，模型所反映的这些场符合实际情况51。为此必须对模型进行检验，即把模型预测的结果和通过某种试验或对含水层施加某种影响后所得到的实际观测曲线或一个地区地下水动态长期观测资料进行比较，看两者是否一致。若不一致，就要对模型进行校正，通常要进行反复地修改参数和调整某些源汇项才能达到较为理想的拟合结果。只有经过识别校正后的模型才能初步说它确实能代表所研究的地质体，是实际水流系统的复制品了，才可以用于进一步的计算或预测。找出模型参数完成识别的基本方法有两大类：（1）试估-校正法，人工调整参数；（2）自动调整参数。本模型的识别与检验过程采用试估-校正法5253。该方法预先给定一组参数的估计值，输入模型中进行计算，比较计算结果与实测结果的误差，如果误差没有达到精度要求，则对刚才输入的一组参数做出调整，直到达到一定的精度时停止计算，这时所用的一组参数值被认为是符合实际的参数值。本文数值模拟结果的可靠性是通过比较计算的地下水水位与实测地下水水位的拟合程度。模型的识别和验证主要遵循以下原则54-56:模拟的地下水流场要与实际地下水流场基本一致，即要求地下水模拟等值线与实测地下水位等值线形状相似；模拟地下水的动态过程要与实测的动态过程基本相似，即要求模拟与实际地下水位过程线形状相似；从均衡的角度出发，模拟的地下水均衡变化与实际要基本相符；识别的水文地质参数要符合实际水文地质条件。根据以上四个原则，对研究区地下水系统进行了识别和验证。通过反复调整参数，识别水文地质条件，确定了模型结构和参数。本次研究，由于缺乏足够的、可以控制全区的长期水位观测资料，无法通过地下水长期动态资料来识别模型，所以采用现场抽水试验期间的水位动态观测资料，将识别时段定为抽水试验期。把各观测孔在抽水试验期间的水位观测资料代入模型中，观察各水位观测孔的计算水位和观测水位拟合情况，作为调整参数的依据。GMS提供了一种比较直观的校正标志功能，借助此功能可以在模型的平面视图中直接查看拟合效果，如图3-4。该标志定义如下：当模型在观测井处求得的误差值（计算值与观测值之差）在设定的允许范围内时，该标志显示为绿色短条，随着误差的增大，其颜色将按照绿色黄色红色渐变，长度也变长。该标志还可以对每一口观测井赋权重，以体现各个观测井对拟合效果的贡献是有差别的。利用校正标志，可以直观地看到各观测井处水位拟合的一个总体情况，方便了调参工作的进行。图3-4 校正标志校正标志只能用来粗略的了解研究区的水位拟合情况，要想定量知道拟合效果，还要通过计算观测水位对比图。本文对各含水层的水平向渗透系数和水平向渗透系数与垂直向渗透系数的比值作为需要调整的参数。在模型识别过程中，根据设定的容许误差，遵循模型识别的原则，根据水位实测值对不同参数进行反复调节和识别，给出了两个参数的范围值，水位识别的结果见图3-5图3-12。从图中可以看出，经过识别的参数计算模拟的水位线与实测水位线趋势相同，误差都在容许范围内 。图3-5 第一含水层1号观测井图3-6 第一含水层2号观测井图3-7 第二含水层1号观测井图3-8 第二含水层2号观测井图3-9 第二含水层3号观测井图3-10 第二含水层4号观测井图3-11 第三含水层1号观测井图3-12 第三含水层2号观测井3.3.4 模型识别评估模型识别是通过实测水头和模拟水头等值线图的对比，提供了一种看得见的定量方法，同时给出了识别误差空间分布的大致信息。但根据现场资料绘制的等值线图包含了绘制引起的误差，不宜作为识别的唯一证据，因此需要对识别结果进行定性和定量的评估。实测水头和模拟水头的散点图是另一种显示拟合程度的方法。各点对直线的偏离应是随机分布的。用列表法表示水头观测值、模拟值以及它们间的误差值和某种误差值的平均是展示识别结果的另一种常用方法。误差值的平均就作为识别过程中的平均误差。识别的目的就是使这个误差最小。一般采用均方根（）误差38来表示水头观测值和模拟值间误差平方和的均值：式中：为识别值的数目。根据模型识别的结果，对各参数组合进行识别评估，各含水层的见表3-2表3-4。从表中所列的，选择最小误差的参数组合作为各含水层的参数值，具体值见表3-5。表3-2 潜水含水层模型参数识别评估汇总表Kh=1.0,Kh/Kv=10Kh=0.9,Kh/Kv=8Kh=0.8,Kh/Kv=6Kh=0.8,Kh/Kv=7Kh=0.7,Kh/Kv=5Kh=0.7,Kh/Kv=6Kh=0.6,Kh/Kv=4Kh=0.5,Kh/Kv=3S3-10.1080.10430.10630.0910.09610.0770.0880.0851S3-20.0420.04460.05040.0420.05040.04040.0530.0594Total0.150.14890.15670.1330.14650.11740.1410.1445表3-3 相对弱含水层模型参数识别评估汇总表Kh=6.8,Kh/Kv=24Kh=6.9,Kh/Kv=24Kh=7.0,Kh/Kv=24Kh=7.1,Kh/Kv=24Kh=6.8,Kh/Kv=22Kh=6.9,Kh/Kv=22Kh=6.7,Kh/Kv=20Kh=6.8,Kh/Kv=20Kh=6.6,Kh/Kv=18S4-10.31780.28110.24080.20740.28460.2490.29330.25690.2983S4-20.33620.35070.37550.38990.37150.3860.38340.39750.3987Total0.6540.63180.61630.59730.65610.6350.67670.65440.697表3-4 承压含水层模型参数识别评估汇总表Kh=21.4,Kh/Kv=6Kh=21.6,Kh/Kv=7Kh=21.6,Kh/Kv=10Kh=21.8,Kh/Kv=8Kh=21.8,Kh/Kv=10Kh=21.8,Kh/Kv=11Kh=22,Kh/Kv=8Kh=22,Kh/Kv=10Kh=22,Kh/Kv=12Kh=22.2,Kh/Kv=16Kh=22.2,Kh/Kv=18Kh=22.4,Kh/Kv=22Kh=22.6,Kh/Kv=22S1-10.25360.18570.21120.12610.13960.14660.08060.08660.09380.08230.08310.09560.1496S1-20.03740.07620.060.12330.11050.10580.17580.16340.1530.19180.18190.21750.2668Total0.2910.26190.27120.24940.25010.25240.25640.250.24680.27410.2650.31310.4164Kh=23,Kh/Kv=8Kh=23.2,Kh/Kv=8Kh=23.2,Kh/Kv=10Kh=23.4,Kh/Kv=8Kh=23.4,Kh/Kv=14Kh=23.6,Kh/Kv=8Kh=23.6,Kh/Kv=10Kh=23.6,Kh/Kv=14Kh=23.8,Kh/Kv=10Kh=23.8,Kh/Kv=14Kh=23.8,Kh/Kv=20Kh=24,Kh/Kv=14Kh=24,Kh/Kv=18Kh=24.2,Kh/Kv=16Kh=24.2,Kh/Kv=18S4-10.15260.1530.15030.170.150.19970.18550.16370.21950.1910.16160.22630.19920.2510.2364S4-20.140.1130.1240.10050.1240.10.09850.1060.1020.09780.1140.1010.09740.10680.1014Total0.29260.2660.27430.27050.2740.29970.2840.26970.32150.28880.27560.32730.29660.35780.3378表3-5 各含水层参数值含水层名称（）潜水含水层0.76-0.06承压含水层22.581.57710-4-相对弱含水层7.0241.45110-5-3.3.5 模型检验由于识别中存在不确定性，通过识别模型得到的一组参数有可能并不精确代表现场的实际数值，为进一步验证所建立的数学模型和模型识别后确定的水文地质参数的可靠性，需要利用已有的另外时段的地下水位动态观测资料对数学模型进行检验。模拟值和现场实测值都应在预先设定的容许误差范围内一致，容许误差应小于拟合计算期间水位变化值的10%，在水位变化值较小（小于5m）的情况下，水位拟合误差一般应小于0.5m。38如所有结果显示两者一致，说明模型可靠，具有实用性，可用于预报。如不一致，则必须对识别所得的上一组参数进行修正，重新进行识别、检验，直至一组新的参数在识别和检验阶段都显示模拟值和实测值两者在预先设定的容许误差范围内一致。本研究区由于缺少长期观测资料，检验通过另外三组独立的稳定流数据来完成。图3-13 第一含水层计算水位与实测水位过程线对比图3-14 第二含水层计算水位与实测水位过程线对比1图3-15 第二含水层计算水位与实测水位过程线对比2图3-16 第三含水层计算水位与实测水位过程线对比通过对地下水计算水位（流场）与实测水位（流场）拟合误差进行统计，模型的识别与检验阶段全部观测孔的水位拟合情况较好，模型的计算值与观测井反映的实测值基本具有相同的变化趋势。拟合误差小于0.5m，能够较好反映出该点水位动态趋势，如图3-13图3-16。说明含水层结构、边界条件概化、水文地质参数的选取及源汇项的确定都是合理的，所建立的数学模型能较真实地刻画研究区地下水系统特征，仿真性强，可以利用建立的数值模型进行地下水资源的评价。经分析，产生误差的主要来源是各源汇项的统计误差、野外观测的水文地质参数的误差，还有就是水文地质边界的确定，也常常带有很大的人为因素和不确定因素。另外，由于地层资料的精度问题导致地层模拟的误差。还有一点需要说明的是，模型算法采用的是迭代求解，通过迭代法所得到的解，仅仅是差分方程的近似解。其精度受很多因素的影响，如选定的收敛指标以及所用的迭代方法本身。即使每个时间段所得到的解是精确的，也仅仅是相对于在该时间段内所建立的差分方程组而言。对于偏微分方程公式来说，它的数值解也是近似解而已。因为通过有限差分法所得到的解与对应的解析解相比，总带有截断误差。一般来说，这种截断误差会随网格间距的增加而增加。即使是获得了对基本偏微分方程的精确解，对于野外条件而言，这种解仍然是近似解而已。 3.3.6 模型预报经过模型的识别和检验，说明所建立的模型可以作为真实地下水系统的代表，可以用其对未来地下水的运动状态进行预报。所谓预报，就是将规划的地下水开采量和经预测得到的地下水各种源汇项代入到已经过识别和检验的数学模型之中，通过模型的运转计算，计算出沉井不同降水深度下的地下水位及流场的响应状况。3.3.6.1 一类边界的预测一类边界处的水位受自然因素控制，在预报阶段，在模型预报之前先对边界条件做出预报。由于在桥位处没有水文测站，缺少水文资料，因此根据邻近测站水文资料推算了多年平均水位和五年一遇水位，分别为3.0m和5.5m。由于泰州大桥处在感潮河段，所以以上水位都为高潮位，计算中间水位由程序自己插值完成。3.3.6.2 模型中源汇项的预报在模型预报之前，应先对模型中的源汇项做出预报，其中主要是对抽水量的预报。根据沉井施工方案，抽水井等间距的分布在沉井四周，每边5口抽水井，共16口。抽水井与沉井间距为10m，抽水井的布置见图3-17。图3-17 研究区施工降水井布置图3.3.6.3 预报时段的确定根据沉井施工方案，施工降水总共分为四个时段，每个时段历时1个月，对预报时段的划分见表3-6。表3-6 预报时段划分预报时段经历的时间2008-3-12008-4-12008-4-12008-5-12008-5-12008-6-12008-6-12008-6-303.3.6.4 模型预报结果根据施工计划，沉井中心处在四个时段降水深度分别为8m、13m、18m、23m，经过模型计算预报，各抽水井的涌水量见表3-7。表3-7 抽水井的涌水量（单位：）抽水井标号第时段涌水量第时段涌水量第时段涌水量第时段涌水量12504506808502250450680850325045068085042504506808505250450680850625045068085071704506808508170230350720917023035050010170230350500111702303505001217023035050013170230350500141702303507201517045068085016250450680850总涌水量32805660857011590根据南北锚碇水文地质勘察报告的北锚碇基坑涌水量的估算值，当基坑水位降至-15m时，基坑内水主要来源于第二层含水层，基坑涌水量由公式（3-9）计算得： （3-9）式中：基坑涌水量，； 基坑水位降深，； 潜水含水层厚度，； 渗透系数， 基坑影响半径，； 基坑等效半径，。当基坑水位降至-23m时，基坑内水主要来源于第二层含水层以及第三层上层，基坑涌水量由公式（3-9）计算得：将模型预报的涌水量与水位降深作拟合，得出当基坑水位降至-15m时，涌水量，与理论公式计算结果的误差为3.66%；当基坑水位降至-23m时，两种涌水量计算结果的误差为20.55%。预报模型考虑了地层分布的不均匀性和渗透系数在不同方向的差异性，并且考虑的边界条件更接近实际情况，经数值模拟预报的基坑涌水量和一般理论公式计算的基坑涌水量是吻合的，误差也是在可以接受的范围内。这进一步验证了模型的可靠性和准确性，为施工的顺利进展提供了可靠的依据。施工降水过程中近长江侧水头剖面、上游方向水头剖面的等水头分布图和平面等水头分布图如图3-18图3-29。图3-18 降水8m时近长江侧水头剖面