（水文学及水资源专业论文）基于gis的威海市海水入侵可视化数值模拟系统研究.pdf

摘要 地下水不仅是水文循环的一个重要组成部分，也是人类不可缺少的一种自然 资源，其水量稳定，水质良好，在社会发展和人民生活中的地位越来越重要了。 随着人类社会的不断发展，人1 ：3 的迅速增长，人们对地下水的需求量日益增大， 开采强度也在不断增加，地下水的超量开采使地下水位急剧下降，导致沿海地区 发生了大面积的海水入侵，成为制约社会经济发展的重要因素。 本文选用美国地质调查局所开发的饱和带非饱和带可变密度的地下水流运 动和溶质或能量运移模型( s u t r a 模型) ，从源程序入手，剖析其程序设计原理、 方法、结构和流程，详细探讨了s u t r a 模型的特点，并根据实际情况对其加以 改进。 出于s u t r a 三维模型参数和数据的复杂性，s u t r a 模型的前处理和后处 理过程中涉及到数量巨大的数据处理和分析工作。本文利用地理信息系统处理海 量数据和空间分析的功能，开发出g i s 与s u t r a 模型有机融合在一起的集成系 统一海水入侵可视化数值模拟系统。 模型改进及系统构建包括：( 1 ) 改进s u t r a 模型的降水入渗补给、潜水蒸 散发消耗和机井抽水量等源汇项计算程序，使之可以有效反映气候变化及地下水 开采对海水入侵的影响；( 2 ) 利用g i s 空间分析、数据管理和可视化表达功能， 实现了s u t r a 模型的可视化构建、含水层空间和模拟时间的可视化离散、参数 的可视化赋值、模型的可视化模拟及其模拟结果的可视化表达等功能，为海水入 侵研究提供了有效的模拟和分析工具。 关键词：地理信息系统，海水入侵，地下水数值模型，s u t r a 模型 a b s t r a c t g r o u n d w a t e rp l a y sa ni m p o r t a n tr o l ei nw a t e rc y c l ea n di so n eo ft h e i n d i s p e n s a b l en a t u r er e s o u r c e s i ti sr e l a t i v e l ys t e a d yi nw a t e rq u a n t i t ya n dc l e a ni n w a t e rq u a l i t y g r o u n d w a t e ri sm o r ea n dm o r ei m p o r t a n tf o rt h es o c i a ld e v e l o p m e n t a n dl i v i n g w i t ht h er a p i dd e v e l o p m e n to fe c o n o m i c sa n di n c r e a s eo ft h ep o p u l a t i o n ， i n c r e a s eo fg r o u n d w a t e rw i t h d r a w a lr e s u l t si nd e c r e a s eo fg r o u n d w a t e rl e v e l ，a n dt h u s l c a d st oe n v i r o n m e n t a ld i s a s t e r , s u c ha ss e a w a t e ri n t m s i o ni nt h ec o a s ta r e a s n l i sh a s b e c o m ea ni n f l u e n c ef a c t o ri nl i m i t i n gs o c i a la n de c o n o m i cd e v e l o p m e n t i nt h i sp a p e r , av a r i a b l ed e n s i t yg r o u n d w a t e rf l o wa n ds o l u t et r a n s p o r tm o d e l ( s u t r a ) w a ss e l e c t e dt os i m u l a t et h es e a w a t e ri n t m s i o n t h ep r i n c i p l e ，m e t h o d ， m o d e ls t r u c t u r ea n dc o m p u t a t i o np r o c e s s e so ft h es u t r am o d e lw e r ei n v e s t i g a t e do n t h eb a s i so ft h ef o r t r a np r o g r a m sp u b l i s h e db yu n i t e ds t a t eg e o l o g i c a ls u r v e y 1 1 1 e m o d e ls t r u c t u r ew a sm o d i f i e dt oe n h a n c ei t sm o d e l i n gc a p a c i t y t h e3 dm o d e lc o m p l e x i t ya r i s e sf r o mm el a r g en u m b e ro fm o d e lp a r a m e t e r sa n d i n p u ta n do u t p u td a t a w jd e v e l o p e dav i s u a l i z e dn u m e r i c a lm o d e l i n gs y s t e mw h i c h i n t e g r a t e sg i sa n ds u t r at of e a s i b l yd e a l i n gw i t hs p a t i a la n a l y s i s a n dd a t a m a n a g e m e n t t h ep r i m a r yi m p r o v e m e n to f t h em o d e la n ds y s t e md e v e l o p m e n tf e a t u r e si n c l u d e ( 1 ) t om o d i f rg r o u n d w a t e rs o u r c ea n ds i n ks u b r o u t i n e sf o rc o m p u t i n gp r e c i p i t a t i o n r e c h a r g e e v a p o t r a n s p i r a t i o nl o s sa n dg r o u n d w a t e rw i t h d r a w a l t h i s m o d i f i c a t i o n e n a b l e st h em o d e lt os i m u l a t ei n f l u e n c e so fc l i m a t ea n dp u m p i n gw e l l so ns e a w a t e r i n t r u s i o n ；( 2 ) t od e v e l o pv i s u a l i z e dm o d e l i n gs y s t e mo nt h eb a s i so fs p a t i a la n a l y s i s f u n c t i o n ，d a t am a n a g e m e n ta n dv i s u a le x p r e s s i o no fg i s t h i sv i s u a l i z e ds y s t e mh a s c a p a c i t yi nv i s u a l i z e dm o d e lc o n s t r u c t i o n ，i n c l u d i n gd a t ai n p u t ，m o d e ls i m u l a t i o na n d o u t p u te x p r e s s i o n t h ed e v e l o p e ds y s t e mp r o v i d e sa ne m c i e n t t o o lf o rs i m u l a t i o na n d a n a l y s i so f t h es e a w a t e ri n t r u s i o n k e yw o r d s ：g i s ，s e a w a t e ri n t r u s i o n ，n u m e r i c a lm o d e l i n g ，s u t r a 本研究得到以下项目资助 教育部新世纪优秀人才支持计划( n c e t - 0 4 0 4 9 2 ) 学位论文独创性声明： 本人所呈交的学位论文是我个人在导师指导下进行的研究工作 及取得的研究成果。尽我所知，除了文中特别加以标注和致谢的地方 外，论文中不包含其他人已经发表或撰写过的研究成果。与我一同工 作的同事对本研究所做的任何贡献均已在论文中作了明确的说明并 表示了谢意。如不实，本人负全部责任。 论文作者( 签名) ： 学位论文使用授权说明 7 年哆月l 日 ，f 河海大学、中国科学技术信息研究所、国家图书馆、中国学术期 刊( 光盘版) 电子杂志社有权保留本人所送交学位论文的复印件或电 子文档，可以采用影印、缩印或其他复制手段保存论文。本人电子文 档的内容和纸质论文的内容相一致。除在保密期内的保密论文外，允 许论文被查阅和借阅。论文全部或部分内容的公布( 包括刊登) 授权河 海大学研究生院办理。 论文作者( 签名) ：p 7 年日 河海大学硕士学位论文第一章绪论 1 1 研究背景 第一章绪论 随着社会经济建设和城市化的快速发展以及人口的迅速增长，人类对地下 淡水的需求量越来越大，开采强度也在不断加强，地下水的超量开采使地下水 水位急剧下降，导致沿海很多地区发生大面积的海水入侵成为制约社会经济发 展的重要因素。 威海市近几年由于工农业生产发展对地下水的需求量不断增加加之遭遇 1 9 9 9 - 2 0 0 1 年严重干旱年份，地下水超采严重，地下水位急剧下降，破坏了内 陆淡水与海水的固有水动力平衡，引起海水入侵。据统计，威海市的三市一区 均存在不同程度的海水入侵现象，全市已有2 0 3 ，4 k m 2 遭受不同程度的浸染，1 4 万亩农耕地处于侵染区，7 8 眼机井报废。由于地下水咸化，粮食平均每年减产 4 2 万吨，工业产值每年损失2 亿元。侵染严重地区，人畜吃水发生了困难，据 不完全统计，侵染区目前已有5 ，8 万人吃水发生困难。为了解决以上问题和困 难，威海市急需合理开发利用地下水资源，在最大限度地满足工农业用水和生 活用水的同时，遏制海水入侵。为此，在现有地下水实时动态监测管理系统的 基础上，需进一步研究威海市海水入侵规律、海水入侵与地下水开发利用的关 系，建立海水入侵预测报警系统、灾害评估系统和地下水资源开发利用管理决 策系统以及地下水资源开发利用管理决策系统，为地下水资源合理开发利用和 沿海生态环境保护提供决策依据。 1 2 海水入侵机理 天然条件下，海岸带含水层中的淡水和咸水f n 维持着一种平衡关系，它们 之间有一个界面，界面以上的淡水流向海岸。抽取淡水后，会引起潜水位( 或 测压水头) 的下降。如果淡水抽水量超过了它的补给量，使得海岸带附近的潜 水位( 或测压水头) 下降到出现淡水体的水头低于附近海水楔形体的水头时， 界面就要向陆地推进，直至形成新的平衡，这种现象称为海水入侵n 可见海 河海大学硕士学位论文第一章绪论 水入侵是在天然和人类活动复合作用下的地下咸水和淡水共同运动的流体动力 学过程。 从物理机制上分析可知，海水入侵主要有两个驱动力 2 1 ： ( 1 ) 咸淡水密度差引起咸水向陆地渗流：由于水的压力p = p g i l ，式中g 为 重力加速度，对同一地理位置g 不变，因此，即使咸淡水界面两侧相距很近的 两点上咸水和淡水具有相同的压力水头h ( 取决于研究点在水面以下的深度) ， 但由于咸淡水具有不同的密度p ，故这两个点上的淡水和咸水却具有不同的静 水压力p ，这个压力差驱使咸水向内陆渗流： ( 2 ) 咸淡水中盐的浓度差驱使咸水向内陆扩散：浓度差直接引起咸水向陆 地扩散，由于介质透水性的不均一性，在扩散过程中发生力学弥散，扩散和力 学弥散合称为弥散。正是这个弥散作用造成了咸淡水之间的过渡带和入侵海水 向大海的回流。 事实上，天然条件下海水入侵之所以不会无限发展，是因为除上述两个动 力外，还有一个阻力，即咸淡水水位高程差引起的淡水向海洋渗流。一般情况 下，淡水位高于海水位，淡水向海洋渗流对海水入侵起到抗衡和抵御的作用， 使入侵和反入侵因子间达到平衡。海水入侵发展全过程的本质，就是这个平衡 的破坏和重建过程。 另外，由于淡水和咸水很容易混合，它们之间的接触带由于水动力弥散作 用常形成一个由淡水、低矿化水逐渐变为高矿化水、咸水的过渡带。过渡带的 厚度和形状由多种因素影响着，如岩性、地质构造、水动力特征、含水层的补 排关系等，当过渡带厚度远小于含水层厚度时，可以把过渡带近似看作是突变 界面，即通常所说的突变界面模型。此模型是一种非常理想的概化模型，其求 解只能获得一种近似结果；反之，当过渡带厚度与含水层厚度比较不可以忽略 时，这种模型称为过渡带模型，这种模型只能用数值法求解，但它可以更好地 揭示海水入侵的入侵机理和入侵规律。由于本文研究区的含水层厚度不大，过 渡带厚度与含水层厚度相比不可忽略，因此本文采用过渡带模型模拟海水入侵， 再现研究区的海水入侵过程。 河海大学硕士学位论文 第一章绪论 1 3 国内外研究现状 全球范围内海水入侵的普遍性，已引起国际社会的共同关注，有关国家对 此积极开展研究和治理工作。国外对海水入侵的研究已有一百多年的历史，最早 提出滨海地区咸淡水界面理论的是荷兰人b a d o n g h y b e n p l 和德i 雪人h e r z b e 唱【3 】， 他们分别于1 8 8 8 年和1 9 0 1 年独立提出了著名的计算咸淡水交界面的g h y b e n h e r z b e r g 关系式，直到现在，仍被沿用着。1 9 6 4 年b e a r 和d a g a n 探讨了 g h y b e n h e r z b e 唱关系式的正确性，发现了在等厚度的水平承压含水层中，确定 了稳定流动下界面尾端的深度并由此确定入侵楔形体长度等【4 】f 5 t 。近二十多年以 来，计算机技术的飞跃发展，出现了一大批海水入侵数值模拟模型，如 3 d f e m f a t 2 d f e m f a t 、c o d e s a 3 d 、d s t r a m ? f a s t - c ( 2 d 3 d ) ，f e f l o w 、 h s t 3 d 、s h a r p 、s u t r a 等。其中饱和带非饱和带变密度地下水流及溶质或能 量运移模型( am o d e ls a t u r a t e d - u n s a t u r a t e dv a r i a b l ed e n s i t yg r o u n d - w a t e rf l o w w i t hs o l v eo re n e r g yt r a n s p o r t ，简称s u t r a 模型) 是国外在海水入侵研究中应 用最广泛的模型口1 ，具有源代码公开、模块化设计、代码适用性强、数值解稳定、 精确清晰、便于用户修改等优点。 国外对海水入侵进行了大量的研究f 3 j ，其中使用s u t r a 模型对海水入侵的研 究主要有： 1 9 9 9 年，b a t a i e a s h t i a n i 等对s u t r a 进行了改进，研究边界条件周期变化对 非饱和带的影响【6 1 ； 2 0 0 0 年，h o w a r dw r e e v e s 等通过修改s u t r a 模型研究总应力变化对盐咸滩 地下水的影响【7 1 ； 2 0 0 0 年，s m e h t a 等对美国南部高原盐碱化进行了研究【8 i 。 我国近年来针对沿海地区海水入侵机制也进行了一系列研究： 1 9 9 1 年，薛禹群等在烟台莱州市建立了第一个符合当地实际情况的三维海水 入侵数学模型，全面研究了当地的海水入侵的现状、规律、水化学特征等1 9 】； 2 0 0 1 年，成建梅、陈崇希等人建立了三维变密度对流弥散水质数学模型来研 究山东省烟台夹河中、下游地区咸淡水界面的运移规律1 0 】； 2 0 0 4 年，丁铃使用s u t r a 二维模型对大连市周水子地区进行海水入侵模 拟，并探讨了一些防治措施【” 。 河海大学硕士学位论文第一章绪论 1 4 研究内容和目的 本文选用威海市高技术产业开发区为研究区，由于该地区含水层较薄， 所以使用过渡带模型对其进行分析研究，在总结前人研究成果的基础上，选 、用美国地质调查局开发的饱和带非饱和带变密度地下水溶质或能量运移模型 ( 简称s u t r a 模型) 为海水入侵数值计算模型。为了处理s u t r a 模型的海 量数据，本文利用g i s 的空间分析功能和处理海量数据的能力，在v i s u a lb a s i c 平台上，使用a r c o b j e c t s 组件开发了基于g i s 的海水入侵可视化数值模拟系 统。具体研究内容如下： ( 1 ) 在研究海水入侵机理的基础上，详细研究s u t r a 模型的机理和输入 输出文件的格式，以及文件中所涉及到的各个量的含义；在此基础上结合威海 市海水入侵的特点对s u t r a 模型进行改进，使之可以方便地利用观测资料( 降 水量、蒸散发量、抽水量等) ； ( 2 ) 利用g i s 的空间拓扑分析和数据管理功能及可视化表达等功能，在可 视化编程语言v i s u a lb a s i c 平台下，对s u t r a 模型进行集成，开发出基于g i s 的威海市海水入侵可视化数值模拟系统； ( 3 ) 将该数值模拟系统在威海市高技术产业开发区进行应用。 本文研究的目的是： ( 1 ) 对s u t r a 模型的改进为海水入侵数值模拟提供理论基础； ( 2 ) 实现海水入侵模型( s u t r a 模型) 的一体化管理和可视化表达，为 海水入侵机理研究和水资源管理提供了技术支持； ( 3 ) 利用g i s 的空间分析技术，使s u t r a 模型的前处理、后处理更加高 效方便。 河海人学硕士学位论文 第二章研究区概况 2 1 社会经济概况 第二章研究区概况 威海市地处山东半岛最东端，位于北纬3 6 0 4 2 至3 7 。3 5 ，东经1 2 1 。4 2 至 1 2 2 。4 2 ( 如图2 1 ) ，东、南、北濒黄海，西与烟台接壤，辖文登、荣成、乳山 三市和环翠区，总面积为5 4 8 5 3 k m 2 ，海岸线长9 8 5 9 k m ，占全省的三分之一。 2 0 0 0 年末，全市总人1 口2 4 6 9 5 万人，其中农业人口1 5 6 5 0 万人，非农业人口9 0 4 5 万人，人口平均密度为4 5 0 人k m 2 。威海市是全国第一个卫生城市，被联合国评 为全球1 0 0 个最适合人类居住的城市之一，经济比较发达，其中机械电子、纺织 服装、交通运输、海洋捕捞和海洋养殖业尤为发达。2 0 0 0 年，全市国内生产总 值5 6 0 9 亿元，工农业总产值1 3 0 6 5 亿元【1 2 1 。 图2 1 研究区位置示意图 河海大学硕士学位论文 第二章研究区概况 2 2 气象 本区位于北温带季风气候区内，具有明显的海洋性气候特征。四季变化和 季节进退都比较明显，由于三面环海，受海洋调节，与同纬度内陆地区相比， 具有雨水适中、空气湿润、气候温和的特点，但四季差异显著，春季受西南大 风影响，气温回升快，空气干燥，蒸发量大，多春旱发生。夏季受东南季风控 制，降水适中，降水量占全年总降水量的6 0 左右，秋季受蒙古高压影响，夏 季风南退，降水减少，冬季受极地大陆气团所控制，冷空气活动频繁，受渤海 暖洋面影响，经常出现冷流降雪天气。威海市全市多年平均降水量为7 6 0 m m ， 降水的显著特点为：一是降水主要集中在汛期，汛期多年平均雨量5 4 5 m m ，占 年降水量的7 2 ；二是年际变化大，丰枯差异悬殊，最大为1 9 5 9 年1 1 2 6 4 m m ， 最小为1 9 9 9 年3 7 4 o m m ，相差7 5 2 4 m m 1 1 2 1 0 该区季节风比较明显，冬季风速最大，春季次之，夏秋季最小。冬春季多 为北风和西北风，风向频率1 0 一4 0 ，最大风速1 2 级，夏秋季以南风和东南风 为主，风速较小，频率l o 1 5 1 1 2 】。 2 3 水资源及其开发利用状况 由于该区的降水特点使得水资源的开发难度加大，利用率降低。据最新水 资源评价，全市多年平均地表水资源量为1 4 2 6 亿m 3 ，地下水资源量为4 7 5 亿 m 3 ，扣除两者重复计算量2 7 3 亿m 3 ，全市多年平均水资源总量为1 6 2 8 亿m 3 。 地下水的来源主要靠降水补给，还受降水量和降水强度的影响。地下水按含水 层的性质划分主要有三种类型：一是河谷平原和滨海地区多为砂层含水层，埋 深一般1 - l o m ，分布均匀，连续呈层状，水量较丰富，是地下水富集区，是农 业灌溉的重要水源；二是丘陵山区的变质岩内裂隙水，地下水埋深一般为1 0 一 3 0 m ，水量较小，且受降水影响较大；三是基岩断裂带的深层水，地下水埋深 3 0 1 4 0 m ，主要分布于较大的地质构造的断层中。自1 9 9 9 2 0 0 1 年的严重干旱， 使威海市的工农业及城市生活用水告急，几座向市区供水的水库几乎出现无水 可供的局面。于是作为储备水源的地下水出现超采、滥采的现象，地下水位急 剧下降。原先向市区供水的杨家滩和西捞台两个水源地被迫关闭【l “。 河海大学硕士学位论文第二章研究区概况 2 4 地貌及地质构造 威海市地处胶东地区的东北部，是一长期隆起地带，属低山丘陵地形，海 拔高度不大，除少数山峰海拔5 0 0 米以上外，大部分为2 0 0 3 0 0 米的波状丘陵 如图2 2 。区内出露的地层以下远古带胶东群的各类变质岩为主。自上远古代至 晚第三纪，一直处于隆起上升状态，遭受风化侵蚀，没有接受沉积，直至新生 代第四纪中更新世开始有残坡积、洪积冲积、海积等堆积层。它们的分布与厚 度显然受古地理条件的控制【l “。 图2 1 2 研究区三维地形图 区内第四系地层主要为中上更新世的残积坡积层、洪积冲积层和全新世的 海相沉积层。残积坡积层、洪积冲积层二者连续过渡，分布于山区河流两岸、 山区盆地、山坡及山麓地带。厚度较火，1 1 0 米不等；海相沉积层主要分布 于滨海平原地带。厚度不一，一般存2 0 米左右，由砂石、淤泥质砂土等组成。 乳山一威海复背斜，为胶东地区古老的构造形式，是一个规模较大的强烈 河海大学硕士学位论文 第二章研究区概况 构造带，对胶东地区东部的构造具有骨架定型作用。威海处于复背斜的北端， 由于多次受岩浆活动的影响，其褶皱形态受到严重破坏。威海北部为一单斜构 造层，岩层走向一般在3 1 0 0 - - 3 3 0 0 ，倾向4 0 。- - 6 0 。，倾角5 0 。- - 6 5 。，局部产状 稍有变化。 在研究区内还发育着一系列的北西向、北东向、近东西向和近南北向的压 性或压扭性断裂和积压破碎带。区域岩浆活动以市区望岛村南中生代燕山期粗 粒花岗岩为代表，属文登岩体的一部分，呈岩柱或岩枝产出，与胶东群变质岩 系侵入接触，在变质岩系中多见片麻理、片理方向侵入的中基性岩脉及后期石 英岩脉、长英岩脉、伟晶岩脉及惶斑岩脉穿插。 2 5 水文地质概况 地下水类型主要是第四系孔隙潜水与微承压水，主要含水层为第四系海相 的粉细砂、中砂；陆相坡洪积的粗砂、粉细砂等，补给来源主要为大气降水、 地下径流，排泄方式主要有大气蒸散发及地下径流。 2 6 地层组成 研究区内出露的地层由上而下依次为： 1 第四系人工堆积层( q 4 叫) 素填土：褐色黄褐色，松散，稍湿，主要由砂及粘性土组成，局部含有 碎石土或少量生活垃圾，顶部有的为植物生长土。几乎全区发育。层厚为1 0 3 o m 。 2 第四系植物层( q 4 p d ) 表土：褐色，稍湿湿，松散，为植物生长土，砂及粘性土组成，富含植 物根系，厚度0 3 1 1 0 m 3 第四系风积层( q 4 e 0 1 ) 中砂：黄色，湿饱和、松散、颗粒均匀，局部含有少量粘性土。厚度约 为1 4 2 1 m 4 第四系沼泽相沉积层( q 4 “) 含淤泥质土粉砂：灰灰褐色，饱和，松散，为植物生长十，含较多植物 河海大学硕士学位论文 第二章研究区概况 根系，有臭味，局部夹淤泥质土薄层。层厚约为o 9 0 1 9 0 m 。 5 第四系新近冲洪积层( q 4 “卜p 1 ) 粉砂：黄褐色，湿，松散，颗粒均匀，局部夹软塑可塑状态粉土薄层， 该层层厚约为0 8 0 2 4 0 m 粉质粘土：褐色，软塑可塑状态，土质均一，含少量小砾石及铁锰氧化 物，局部夹粉细砂及粉土薄层，层厚约o 4 0 2 o o m 6 第四系海相堆积层( q 4 “) 粉细砂：灰灰褐色，松散稍密，湿饱和，颗粒较均匀，分选一般， 含半腐烂植物根系及贝壳碎片，具腥臭味，局部夹黄褐色中粗砂薄层，层厚 1 1 0 2 3 0 m 中砂：灰黄黄褐色，稍密中密，饱和，主为长石、石英颗粒，颗粒不 均匀，分选性较差，含少量贝壳碎片及长英质小砾石，局部夹灰色粗砂及粉细 砂薄层。层厚1 , 9 0 3 7 0 m 细砂：灰色，饱和，稍密中密，颗粒较均匀分选性良好，含较多贝壳碎 片，具腥臭味。层厚1 2 4 3 m 粉质粘土：灰灰褐色，软塑状态，土质较均一，含沙量较高。层厚o 8 0 2 o o m 淤泥质土：灰黑色，软塑，局部可塑，土质细腻，含腐烂植物碎屑，有臭 味，局部含较多细砂颗粒。层厚o 6 0 1 6 0 m 7 第四系海陆交互相沉积层( q 4 “) 中砂：黄褐色，饱和，中密，颗粒均匀，混粘性土。最大层厚3 7 0 m 粉质粘土：灰灰褐黄褐色，可塑，局部软塑，上部土质较均一，下部 含砂颗粒较多。层厚2 7 0 4 7 0 m 8 第四系陆相坡洪积层( q 4 d l + p 1 ) 粗砂：黄褐色，中密，饱和，颗粒较均匀，分选一般，含长英质小砾石， 局部夹含粘性土粉砂小薄层。层厚o 6 0 2 1 0 m 粉细砂：黄褐色，中密密实，顶部松散稍密，饱和，颗粒均匀，分选 性良好，局部含较多粘性土，夹含粘性土粉砂小薄层。层厚4 4 0 4 9 0 m 含粘性土细砂：黄白色，饱和，稍密，颗粒较均匀，含粘性土 河海大学硕士学位论文 第二章研究区概况 粉质粘土：灰绿黄褐色，可塑硬塑状态，土质均一，含少量小砾石及 铁锰氧化物，局部含砂量高，局部含较多长英质碎砾石，夹中细砂及碎石薄层。 层厚0 6 1 2 9 m 粉砂：灰黄色，饱和，稍密中密，颗粒均匀，层厚0 7 0 2 1 0 m 碎石：黄色，饱和，密实，碎石成分主要为强中风化斜长片麻岩及长英 质碎块，砂及粘性土填充。最大厚度为o 8 0 m 含粘性土粉砂：黄褐色，中密，饱和，颗粒较均匀，分选性好，含少量长 英质小砾石。层厚1 7 0 m 中粗砂：黄褐色，饱和，中密密实，颗粒极不均匀，分选性较差，含碎 石块，最大粒径为6 c m ，底部为碎石薄层，碎石成分主要为长英岩块和斜长片 麻岩。层厚2 0 0 4 o o m 9 下元古界胶东群变质岩系( p 。1 ) 强风化斜长片麻岩：黄褐色，主要矿物成分：长石、石英及少量黑云母， 中细粒变晶结构，片麻状构造，强风化程度。层厚：1 0 0 2 5 0 m 河海大学硕士学位论文第三章s u t r a 模型及其改进 3 1 概述 第三章s u t r a 模型及其改进 s u t r a 模型【1 4 1 ”1 是美国地质调查局公开发表、c l i f f o r di v o s s 等人开发的 用于饱和带非饱和带可变密度的地下水流和溶质或能量运移的计算机程序，可 进行二维或三维海水入侵模拟计算，其结构高度通用和模块化，允许用户直接 添加新的模拟方法或过程。此外，该模型采用f o r t r a n9 0 的编程语言，强调程 序的清晰性和模块性，这些非程序的效率，便于用户根据自己的需要对程序进 行修改。 s u t r a 模型采用有限单元和有限差分混合法离散控制方程，模拟饱和带、 非饱和带变密度地下水流动和地下水中溶质或能量运移两个过程，可以模拟化 学物质运移过程、吸附、生成以及降解等诸过程，应用在地下水污染物运移和 含水层恢复设计，以及垃圾填埋场周围变密度淋滤液的运动模拟、抽水井附近 含水层中海水入侵的剖面模拟( 咸淡水的突变界面或咸淡水之间存在扩散带) 。 s u t r a 模型的能量运移模拟程序可用于模拟含水层热能分布、地下热传导、含 水层热能存储系统、地热库、含水层热污染以及天然水文地质对流系统u 5 i t l 6 】。 s u t r a 模型的输入文件采用a s c i i 码文件，可以使用任何文本编辑器加 以生成和修改，文件内部还可以添加注释，便于用户解读、理解和修改；对于 输入的数据，s u t r a 模型具有检查功能，在模拟过程中，还可以存储模拟的中 间结果，用于在模拟过程中的任意中断时间段的重新模拟；s u t r a 的输出结构 既包括任意时刻的流体速度等矢量数据，也包括流体压力、溶质浓度或含水层 能量、饱和度和流体质量和溶质质量或能量均衡，以及系统内各观测点的计算 结果。s u t r a 模型中的网格建立非常灵活，适用于任意几何形状的模拟区域。 在二维笛卡儿系统或柱坐标系统中，采用四边形有限元剖分；在三维系统中， 采用六面体有限元剖分。含水层和流体的物理参数，如渗透系数，可以为各向 同性的，也可以是各项异性的，在整个模拟系统中其大小和方向均有变化。边 界条件、源汇项可以是固定值，也可以是时间相关的。 河海大学硕士学位论文第三章s u t r a 模型及其改进 3 2 数学模型 对于考虑过渡带的海水入侵问题，s u t r a 模型采用两个方程来描述【”1 ： 一是地下水流动力方程，该方程用来描述溶质浓度的不断变化，引起密度的不 断变化，从而影响流体压力的变化；二是溶质质量平衡方程，该方程是用来描 述地下水中的溶质随地下水的流动的运移过程。 基于压力的地下水流动力方程为： ( & 峨+ 印i a s w j 、百 + ( 瓯塑o c ) 丝a t 一里- ( 警 匦一p 纠= g cs ， 式中：s 。【1 为指饱和率；【l 】为孔隙度：p m ( l s 2 ) 】为压力：e l m 。m 】为溶质质 量分数( 溶质质量与流体总质量之比) ；t s 】为模# a h , y f 日q ；“l 2 】为渗透率张量； k ， 1 】为相对渗透率； m l s 】为水动力粘滞系数；g l s 2 】为重力加速度张量； o 。 m ( l 3 s ) 为水质量源汇项；川m ，l 3 为流体密度，可以表达成： p 兰风+ 罢( c 0 l 式中：c o m s ，m 】为参考流体溶质质量分数；p o m l 3 】为参考流体密度，当c = c o 时；罢为密度随浓度的变化量，一般假设为固定值，如淡水和盐水在2 0 。c d l 混合，c o = o ，p o9 9 8 2 k g m 3 】时，罢z7 0 0 峙m 3 ；s o p m ( l s 2 ) 】为压力贮 0 l 水率， s 。= ( 1 一b + 妒 ( 3 3 ) 式中：口 m ( 1 s 2 ) 】_ 1 为多孔介质压缩系数； m ( i s 2 ) 】。为流体压缩系数； s u t r a 模型中使用通用的达西定律来描述多孔介质中地下水流的运动： 一 戋卜哩) 4 ， 式中：! l ，s 】为平均流体速度。 溶质质量平衡方程： 河海大学硕士学位论文第三章s u t r a 模型及其改进 堕兰趔；地二兰超生】_ 西o t( 3 5 ) 一里。p 雯) + 里。p ( d ，+ 里) v c ) + 岱。+ ( 1 一占h + g c 式中：p 。【m 。乓】为固体颗粒的密度；c m 。，m g 】为吸附在固体颗粒上的比浓 度；d 。【l 2 s 】为多孔介质中分子扩散系数：“1 为单位张量矩阵( 对角线值为l ， 其余为零；二维空间中为2 x 2 矩阵；三维空间中为3 x 3 矩阵) ；d l 2 s 为弥散 系数张量( 二维空间中为2 2 矩阵；三维空间中为3 x 3 矩阵) ；l ， m ：m s 】为溶 质生产反应中流体的溶质源汇项；一【m ；, q v l o s 】为吸附物质本身生产反应中的吸 附质量源汇项。 3 3 数值解法 s u t r a 模型采用有限元和有限差分混合算法，在加权剩余方法体系中应用 求解。当对研究区域进行适当的时空离散化后，这种方法非常准确。标准有限 元逼近只用于描述流体质量、溶质质量和能量的平衡方程条件下，其它所有非 通量条件都使用积分有限差分方法的有限元网格形式近似表示。这种方法道理 上浅显易懂，是一种最简单、最经济的方法，它既保留了有限元方法的数学精 确性和几何灵活性，还具有有限差分方法的有效性。本节主要介绍和分析 s u t r a 模型的数值解法，为了更好地理解数值解法，先解释一下对控制方程离 散时所要用到的一些概念。 3 3 1 结点、计算单元和均衡单元 s u t r a 模型中把模拟区分成结点和计算单元，所有相邻的计算单元的一部 分组成的均衡单元，结点是与均衡单元相对应的，二维或三维中每个计算单元都 分成四个( 如图3 1 ) 或八个部分( 如图3 2 ) 。结点的参数赋值( 如s u t r a 模型 中的压力、饱和度和浓度或温度) 是通过结点的均衡单元与参数分区叠加分析 进行参数赋值。 在二维模型中，计算单元是四边形( 如图3 1 ) ；三维模型中，计算单元是 广义的六面体( 如图3 2 ) 。因为三维六面体的各个面可以是曲面，所以计算单 河海大学硕士学位论文第三章s u t r a 模型及其改进 元不是严格的六面体。这里需要注意的是所有计算单元组成的雅克比矩阵的行 列式不为零。2 】，所以计算单元不可能存在两个相邻结点组成的线或三个相邻结 点组成的面可以置换计算单元内部的线或面【1 7 1 。 计算单元 均衡单元 结点 图3 1 s u t r a 二维网格中的结点、计算单元和均衡单元【1 7 l 计算单元 均衡单元 结点 图3 2s u t r a 三维网格中的结点、计算单元和均衡单元【1 7 】 3 3 2 数值解法 s u t r a 对地下水流动力方程的求解是利用有限元和有限差分混合的方法， 结合剩余加权法，把地下水流动力方程离散成下面方程【1 5 ： 芸 ( 等一心毛pn 。， 叫剁州+ ( 筹卜时( 射 式中：滓【1 ，n - n ；n n 是结点总数；当i ：j 时，民2 0 ；当i 2 j 时，6 u 5 1 母卜，+ 哮 ， ( 3 7 ) 河海大学硕士学位论文第三章s u t r a 模型及其改进 c p , = ( 岱。詈) ， ( 3 8 ) = 圳( ( 丝) ) ( 臀卅，卜蝴 ， = 川陋) ( k ) ) ) 乳一 n 盱m ( ( 酬铷甄卜蛐 t ， 。= “水丝) ) ( ) m 耕卿舳 其中b g ，y ) 兰n nb ，谚b ，y ) ，b 。是结点i 的网格厚度，一k l 是渗透率张量，它 的单元离散后表示为( ( 丝) ) q 【m s 是均衡单元i 的总质量源汇项，包括井流量、潜水蒸散发量、降水入渗 量等：p 。 m l 2 黼i 上确定压力值：v 。， l s 】是均衡单元i 上基于确 定压力源的压力传导系数：一【l 3 】是均衡单元i 的体积；其余符号表示同上。 3 4 程序设计结构及其流程图 s u t r a 具有程序设计结构模块化、从上而下易于代码解读、容易逻辑跟踪 和最终修改的特点 1 5 1 。 s u t r a 中，每个子程序清晰明了，为用户编制了可以用来修改的子程序， 以便对随时间变化的边界条件和源汇项进行控制。此外s u t r a 还提供了非饱 和带水流模拟的子程序，用户可以根据研究区的不同，编制或修改适合研究区 实际情况的非饱和带水流模拟的子程序。 s u t r a 的输入输出也是模块化设计的，其输入文件包括：文件分配文件( 即 f i l ea s s i g n m e n t sf i l e ，定义了整个模拟过程中需要使用的文件类型，及相应的 文件代码和文件名) ；主输入文件( 即m a i ni n p u tf i l e ，包含了所有模拟过程中 河海大学硕士学位论文第三章s u t r a 模型及其改进 描述模型结构的信息、模拟过程中的各种参数、对模型的时空离散化等) ；初始 条件文件( 即i n i t i a lc o n d i t i o n sf i l e ，包含了每个结点的初始压力值和浓度值( 或 能量值) 。输出文件包括：主输出列表文件( 即m a i no u t p u tl i s t i n gf i l e ，列出 用户需要的一系列信息，如：输入参数、流体和溶质质量均衡计算值等) ；重启 条件文件( 即r e s t a r tc o n d i t i o n sf i l e ，用来存储最后时间段的计算结果，当模拟 热启动时可以作为初始条件使用) ；结点结果文件( n o d e w i s er e s u l t sf i l e ，用来 存储用户指定的时间段序列中每个结点的坐标值、压力、浓度和饱和度的计算 结果) ：计算单元结果文件( e l e m e n t w i s er e s u l t sf i l e ，用来存储用户指定的时 间段序列中每个计算单元的形心坐标值、流速各分量的计算结果) ：观测点结果 文件( 即o b s e r v a t i o n sf i l e ，用来输出用户指定的含有观测资料的结点的计算结 果) ；模拟概要文件( 即s i m u l a t i o ns u m m a r yf i l e ，用来输出模拟过程中的收敛 信息和出错信息，默认文件名是“s u t r a s m y ”) 。 程序的主要逻辑流程图如图3 3 。s u t r am a i n 控制着模拟过程的开始和 结束，读取主输入文件中的1 - 8 个数据集( 部分通过子程序i n d a t 0 来读取) ， 然后使用输入数据来计算各种数组的维数，并为其分配存储空间，通过 c o m m o n 把这些信息传给子程序s u t r a 。子程序s u t r a 通过调用i n d a t l 、 s o u r c e 、b o u n d 和c o n n e c 等子程序来初始化一些常量和读取主输入文件 的数据。通过子程序b a n w i d 来计算带宽，用p t r s e t 来设置储存矩阵的指 针，然后直接调用i n d a t l 来读取初始条件文件的数据，在调用o u t l s t 2 ( 对 二维模型) 或o u t l s t 3 ( 三维模型) 来向“1 s t ”文件打印初始条件。 如果要计算稳定流的压力，子程序会自动确定循环参数，并调用“z e r o ” 来初始化数组。对于非稳定流的压力求解，要求用户设定时间步长循环参数， 确定每个时间步长中应该求解压力还是浓度( 温度) ，还是两者都求，具体判断 是由n p c y c 和n u c y c 两个变量以及子程序s u t r a 中的开关变量m l 决定的。 加 压力和浓度( 温度) 同时求解 仅求解压力 ( 3 1 3 ) 仅求解浓度( 温度) 稳定流与非稳定流的控制变量为i s s f l o ，具体如下 m“坦， l lm 河海大学硕士学位论文第三章s u t r a 模型及其改进 1 0 i s s f l o ： 1 1 2 假定为非稳定流 压力时间步长之前，假定为稳定流( 3 1 4 ) 压力时间步长开始之后，假定为稳定流 这里要注意的是当时间段i t 为0 时，求解稳定压力，而对第一个运移时间 步长是单个递增的。 河海大学硕士学位论文第三章s u t r a 模型及其改进 ： 幢 p t r s e t - t 1 1 d i o u t l s t 2 z e r o t r j s e t t o d ) o u t l s t 3 t 二堡亟阐e * 圆圈 _ o u t r s t 图3 3s u t r a 逻辑流程图( 据j b e a r 修改) ，其中f x o p e n 和i n f s s 是改进 后添加的子程序 子程序s u t r a 在新的时问段的末尾时刻增加一个模拟块t s e c ，在下一个 河海大学硕士学位论文第三章s u t r a 模型及其改进 时间步长开始时，将当前的向量转化成下一个时间步长的初始向量。如果存在 随时间变化的源汇项或边界条件，就会调用b c t i m e 子程序进行处理；若要考 虑吸附作用，则要调用子程序a d s o r b 进行相应的计算。当进行逐个计算单元 计算时，就要调用e l e m n 2 ( 对二维模型) 或e l e m n 3 ( 对三维模型) 构造矩 阵方程。子程序n o d a l 用来为全局矩阵进行结点或均衡单元计算。子程序b c 是根据边界条件的