（水文学及水资源专业论文）秦皇岛地区海水入侵的二维数值模拟.pdf

摘要 随着沿海地区社会经济的发展和人口数量的增加，地下淡水资源对人类的生存和发展 越来越重要。然而由于各种自然因素和人类活动的影响，滨海含水层中海水入侵现象不断加 剧。奏皇岛洋、戴河地区由于其特殊的地理位置，海水入侵问题更是备受关注。 本文建立了一个研究海水入侵问题的二维平面非稳定流数学模型，以秦皇岛洋、戴河平 原区为例，糨合溶质运移方程和水流方程，模拟出了地表河流与地下水系统的相互作用。模 型两部、东部雨北部以缸山和丘陵为边界，南部以海岸线为边界，采用三结点三角形单元对 计算区进行剖分，前沿洋、戴河两岸和抽水井附近加密了网格，然后采_ l = j 不等距时间步长的 上游加权迎风有限元法对方程组进行求解。 文中首先利用a r c v i e w 软件建立了模型计算区内各类水文地质参数的地理信息系统。该 系统不仅便于本模型的各类数据输入，而且便于今后区内已有数据的空间检索。 其次模型通过研究区内1 9 9 8 年5 月2 0 日到2 0 0 3 年1 2 月3 1 日的观测井实测数据校正 和检验了模型的各类参数和边界条件，然后对区内5 年之后海水入侵进行了预测。结果发现 按照2 5 保证率的降雨量使得海水入侵的面积增加2 4 9 4 k m 2 ，按照5 0 保证率的降雨量使 得海水入侵的面积增加2 5 0 5k m 2 ，按照7 5 保证率的降雨量使得海水入侵的面积增加2 5 5 8 k m 2 。 在该模型的基础上仍以1 9 9 8 年5 月2 0 日到2 0 0 3 年1 2 月3 1 日为模拟时段，对三种防 治海水入侵的t 程方案进行了数值模拟，结果发现在洋河口修建防潮坝使得海水入侵的面积 减少3 ，6 3k m 2 沿海岸线兴建地下水坝可使海水入侵的面积减少3 8 6k m 2 ，通过沿岸布置注 水井更可使海水入侵豹面积减少4 1 8k m 2 ，将三个防治工程结合起来可使海水入侵的面积减 少3 2 4k m 2 。 最后鉴于实际地质参数的空间变异性，本文以t h e i s 模型为简单的算例，结合有限单元 法初步探讨了遗传算法在水文地质参数反演中的适用性，并进一步研究了遗传算法中各种参 数取值的变化对于反演结果的影响。灵敏度分析的最终结果表明，实测数据的测罱误差对于 反演参数的精度影响晟大。 关键词：海水入侵，数值模拟有限单元法，遗传算法 a b s t r a e l a b s t r a c t w i t ht h ed e v e l o p m e n to fe c o n o m i ca n dt h eg r o w t ho fp o p u l a t i o n ，g r o u n d w a t e r b e c o m e sa l li n c r e a s i n gi m p o r t a n tr e s o u r c et h a tl i m i t sh u m a ns u r v i v ea n dd e v e l o p m e n t d u et ot h ei m p a c to fn a t u r a lf a c t o ra n dh u m a na c t i v i t i e s ，h o w e v e r , s e a w a t e ri n t r u s i o n h a si n c e s s a n t l yd e t e r i o r a t e di nt h ec o a s t a la q u i f e r t h er e g i o nl o c a t i n gi nt h ep l a i n b e t w e e nt h ey a n gr i v e ra n dt h ed a ir i v e ro fq i n h u a n g d a oh a sa t t a c h e dm o r ea n dm o r e a t t e n t i o n t a k et h i sr e g i o nf o ra ne x a m p l e ，t h et e x ti n c o r p o r a t e st h ee f f e c to fd e n s i t y - d e p e n d e n t f l o wt oc o n s t r u c tat w o d i m e n s i o n a lu n s t e a d ys t a t es a l tw a t e ri n t r u s i o nm o d e lw h i c hi s a b l et os i m u l a t et h ei n t e r a c t i o nb e t w e e n 也er i v e ra n da q u i f e r f i r s t l yt h i sm o d e lt a k e i n t oa c c o u n tt h en a t u r a lm o u n t a i nb o r d e r sa se a s t e r n ，n o r t h e r na n dw e s t e r nb o u n d a r y s e c o n d l yt h em o d e li sd i s c r e t i z e db yl i n e a rt r i a n g l ee l e m e n ta n dr e f i n e da r o u n dt h e p u m p i n gw e l l a n dt h er i v e rb a n k t h i r d l yt h em o d e lu s et h es t r e a m l i n eu p w i n d g a l e r k i nf i n i t ee l e m e n ta p p r o a c ht os o l v et h ep a r t i a ld i f f e r e n t i a le q u a t i o n st h a tt i m e a r ed i s c r e t i z e di n t ov a r i a b l et i m es t e p s a r c v i e ws o f t w a r ei sc o n s i d e r e dt oc o n s t r u c tas u i to fg i sw h i c hi n c l u d ea l lk i n d so f h y d r o g e o l o g yp a r a m e t e rn e e d e df o rt h i sm o d e l t h i ss u i to fg i s i sn o t o n l y c o n v e n i e n tt oi n p u td a t a b u ta l s ot os c a nt h e s ee x i s t i n gd a t a 1 1 1 i sm o d e 】c a l i b r a t e sa n dv e r i f i c a t e st h e s ep a r a m e t e ra n db o u n d a r yc o n d i t i o n sf r o m t h e s eo b s e r v a t i o n v a l u e sw h i c hs p a nt i m ef r o mm a y2 0i n1 9 9 8t od e c3 1i n 2 0 0 3 n l e nt h ea r e ao fs a l tw a t e ri n t r u s i o ni sp r e d i c t e di nt h ef u t u r e5y e a r s t h e r e s u i t sd e m o n s t r a t et h a tt h ea r e ae x p a n d2 4 9 4 k i n 2u n d e rt h e2 5 r a i np r o b a b i l i t y c o m p a r e dt o 也ea r e ai n2 0 0 3 a 1 s ot h ea r e ae x p a n d s2 5 0 5 k i n zu n d e rt h e5 0 r a i n p r o b a b i l i t ya n d2 5 5 8 k r n ou n d e rt h e7 5 r a i np m b a b i l i t y b a s e do nt h i sv e r i f i c a t e dm o d e la n dt a k et h et i m es p a nf r o mm a y2 0i n1 9 9 8t od e c 31i n2 0 0 3a st h es i m u l a t i o np e r i o d ，n u m e r i c a ls i m u l a t i o no ft h r e et y p e so fs e a w a t e r i n t r u s i o np r o t e c t i o np r o j e c t i o nw e l ef m i s h e d t h er e s u l t sd e m o n s t r a t et h a tb u i l d i n g t h ed a mi ny a n gr i v e rm o u t hd e c r e a s et h ea r e ao fs a l tw a t e ri n t r u s i o nb y3 6 3 k m 2 l j c o m p a r e dt ot h ea r e aw i t h o u tp r o t e c t i o np r o j e c ti n2 0 0 3 a l o n gt h es h o r eb u i l d i n gt h e u n d e r g r o u n dw a l lm a k et h e 船o f s a l tw a t e ri n t r u s i o ns h r h l k3 8 6k m 2a n dl a y i n go u t i n j e c t i o nw e l lm a k et h e 掀s h r i k4 1 8i o n 2 f i n a l l yi ti sd i f f i c u l t yt oo b t a i nad e t e r m i n a t ev a l u ef o rt h es p a t i a lv a r i a b i l i t i e so ft h e r e a lh y d r o g e o l o g i cp a r a m e t e r t a k et h et h e i sm o d e la sas i m p l ee x a m p l e ，t h i st e x t p r i m a r yi n v e s t i g a t et h e 印p l i c a b i l i t yo fg e n e t i ca l g o r i t h mt o i n v e r s em o d e l i n go f h y d r o g e o l o g i cp a r a m e t e r a f t e r w a r dt h es e n s i t i v i t ya n a l y s i so fk i n d so fp a r a m e t e r u s e di ng e n e t i ca l g o r i t h mi sc o n d u c t e d f i n a l l yr e s u l t si n d i c a t et h a tt h em e a s u r e m e n t c n o ri st h em o s ta f f e c t i v ef a c t o rt ot h ep r e c i s i o no fi n v e r s em o d e l i n go fp a r a m e t e r k e yw o r d 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 ls i m u l a t i o n ，f i n i t e e l e m e n tm e t h o d ， g e n e t i ca l g o r i t h m 1 1 1 学位论文版权使用授权书 本人完全了解同济大学关于收集、保存、使用学位论文的规定， 同意如下各项内容：按照学校要求提交学位论文的印刷本和电子版 本；学校有权保存学位论文的印刷本和电子版，并采用影印、缩印、 扫描、数字化或其它手段保存论文；学校有权提供目录检索以及提供 本学位论文全文或者部分的阅览服务：学校有权按有关规定向国家有 关部门或者机构送交论文的复印侔和电子版；在不以赢利为目的的前 提下，学校可以适当复制论文的部分或全部内容用于学术活动。 学位论文作者签名：槐镟 川年 月，l 日 本授权书。 指导教师签名：刍l ；着蔓学位论文作者签名：钝瞧 同济大学学位论文原创性声明 本人郑重声明：所呈交的学位论文，是本人在导师指导下，进行 研究工作所取得的成果。除文中已经注明引用的内容外，本学位论文 的研究成果不包含任何他人创作的、已公开发表或者没有公开发表的 作品的内容。对本论文所涉及的研究工作做出贡献的其他个人和集 体，均已在文中以明确方式标明。本学位论文原创性声明的法律责任 由本人承担。 签名：毒也使 妒f 年，月，日 第1 章绪论 1 1 前言 第1 章绪论 海水入侵是指由于各种自然园素和人类活动的影响，滨海地区的地下水含水层中原有的 咸、淡水体之间动态平衡遭到破坏而引起海水进一步侵入到淡水含水层系统，直至两者达到 新的平衡状态这一环境水文地质现象。海水入侵是现代社会所特有的环境问题，它使得我国 部分沿海地区的环境遭受到严重破坏，阻碍了当地工农业和城市经济的发展，海水入侵问题 已经引起国际社会的广泛关注。 我国自7 0 年代以来陆续出现零星的海水和卤水入侵，进入8 0 年代后随着沿海地区的改 革开放，工农业经济建设的迅速发展人类活动地增加使得海水入侵的范屡逐渐扩大。目前 情况较严重的有山东的莱州、龙1 2 、寿光、昌邑，河北的秦皇岛，辽宁的大连，浙江的宁波 以及广西的北海等地。 长期以来，由于秦皇岛沿海经济的高速发展、人口的急剧增加，地下水的开采量逐年增 大，使得地下水降落漏斗不断变深变太，研究区内原有的咸淡水平衡被彻底打破，该问题的 直接后果是导致了成淡水界面持续的向内陆方向纵深。伴随着海水入侵的发生，秦皇岛洋、 戴河平原区的水质开始恶化土壤盐碱化，大量的抽水井报废，取水成本逐年上升。同时海 水入侵还进一步引发了一系列的环境问题，如铁路路基、建筑物、地下管道等下沉开裂，河 道堤防和行洪能力出现危机，大量的工业冷却设备被腐蚀。农田的土壤盐碱化程度加深等。 查清海水入侵的发展过程和演变趋势，科学合理的提出海水入侵的防治方案已经成为滨海地 区人们生存和发展的迫切要求。 1 2 国内外海水入侵相关问题的研究现状 由于海水入侵给滨海地区人类的生产生活造成了严重的影响，因此无论是在国内还是在 国外海水入侵问题都引起人们的广泛关注，例如自1 9 6 8 年起国际上每隔两年定期召开一次 的海水入侵会议( s e aw a t v ri n t r u s i o nm e e t i n g ) 主要关注于海水入侵的实际监测、数值模拟 和地下水资源管理这几个方面，其中数值模拟在解决海水入侵问题中有着独特的优势。一般 地模拟海水入侵海题的数学模型可分为两类。第一类为突变界面模型，即假设淡水和海水是 互不溶混流体，两者之间存在一个突变界面，流体动力弥散作用忽略不计。而且裘布依假定 近似成立，即无论在海水区还是在淡水区，水头随深度的变化可忽略不计，海水和淡水基本 是水平流动，同时含水层的水力传导系数和比贮水系数都不随深度变化，含水层的底板不透 水。 然而实际上海水和淡水是可以相互溶混的，因此他们之问必然存在一个由于流体动力弥 散作用而形成的过渡带，通过这个过渡带，水中的盐分从海水的浓度变到淡水的浓度。一些 第1 章绪论 研究结果表明过渡带的变动范围很大，在某些地区过渡带只有十几米宽，和含水层的整个厚 度相比可以忽略不计，此时可近似认为海水与淡水之间存在一个突变界面；但是在另外些 地区，过渡带的宽度可达到百米以上，不存在突变界面。因此产生了第二类模拟海水入侵问 题的数值模型溶质对流弥敖模型，该模型认为海水和淡水是互溶的，由于水动力的弥散 作用，二者之间存在一个浓度由淡水变化到海水的过渡带。 1 2 1 国外海水入侵相关问题的研究现状 早在1 8 世纪末】9 世纪初b a d o n - o h y b e n 1 8 8 8 】和h e 咖f 1 9 0 1 】最先提出突变界面假设， 并按照静水压强分布求解咸淡水界面问题。随后p i n d e r 和p a g e 1 9 7 6 对存在的突变界面情形 给予了完整的数学描述，b e a r 1 9 7 9 1 列出了许多应用突变界面假设的解析解。然而由于实际 地质条件的复杂性和方程的高度非线性，使得求取解析解遇到了困难，于是数值模拟成为海 水入侵研究的主流。c o l l i n s 等 1 9 7 9 1 及j a m e s 等 1 9 8 0 j 分别对承压含水层的突交界面问题进 行了研究，r a l p h 等 1 9 6 8 1 、p a n i g r a p h i 1 9 8 8 等分别对潜水突变界面问题进行了研究。r u b i n 和p i n d e r 1 9 7 7 、d i e r s c h 【1 9 8 4 1 、d i e s c h 和n i l l e r t 1 9 9 0 、r e i l l y 和g o o d m a n 1 9 8 7 】、 h o l z b e c h e r 1 9 9 5 对抽水井所引起的海水入侵升锥问题进行了研究。m e r c e r 1 9 8 0 】、h e r b e r t 等 1 9 8 8 、d o r g a r t e n 和t s a n g 1 9 9 1 】、o l d e n b u r g 和p r u e s s 1 9 9 5 、o l d e n b u r g 等 1 9 9 6 1 、j o h n s 和r i v e r a 1 9 9 6 、k o l d i t z 等 1 9 9 s 、k o n i k o w 等 1 9 9 7 】模拟了滨海含水层污灌所引起的咸淡 水界面问题，h o l z b e e h e r 1 9 9 8 和y o u n e s 等 1 9 9 9 对这个问题给出了很好的综述性讨论。l i n 和c h e n g 1 9 8 1 利用边界积分法求解了移动的突变界面问题，f r i n d 1 9 9 2 研究了可变密度地 下水流的长期瞬态模拟问题，j o s s e n l i n 和j o n g 等【1 9 8 6 对突变界面的横向弥散问题进行了研 究，b e a r 1 9 9 9 、c h e n g 和q u a z a r 1 9 9 9 分别对突变界面的模型概化及其解析解做了较为全 面的回顾和阐述。b a k k e r 2 0 0 3 1 基于d u p u i t 近似给出了一个模拟海水入侵的三维突变界面模 型的新表达式。 如前所述，突变界面模型以咸淡水互不溶混为前提假设显然是符合实际的。因而自6 0 年代起学界开始转入对溶质对流一弥散模型的研究。h e n r y 9 6 q 以承压、稳定流、均质、 简化边界条件为基础，首次求得一个与海岸线正交的垂向剖面上盐分浓度的解析解，h 蜘l y 模型也因此成为过渡带模型数值实验的一个基准( 成建梅，1 9 9 9 ) 。p i n d e r 和c o o p e r 1 9 7 0 1 则首次将可溶混成淡水流理论应用于沿海实际问题。l e e 等【1 9 7 4 】对二维剖面的h e n r y 问题 做了进一步研究，s e g o l 、p i n d e r 和g - r a y 1 9 7 5 用有限元法求解了非稳定流状态下盐分浓度 在剖面上二维分布，之后s e g o l 和p i n d c r 1 9 7 6 1 3 h j f 究了海水入侵的非稳定流数值模拟问趣。 h e i n f i c h 等 1 9 7 7 提出了用逆风式有限元法求解对流占优的海水入侵问题。v o l k e r 和 r u s l l t f 1 9 8 2 】对海水入侵问题中某些参数的重要性做了评估，并对突变界面模型和弥散模型 进行了对比。h u y a k o m 等 1 9 8 7 以参考水头和浓度为自变量建立了剖面二维有限元模型，证 明可推广于三维的情况，因而又进行了三维海水入侵数学模型的有限元计算。d i e r s c h 1 9 8 8 】 对由于密度差引起的海水洄流问题进行了有限元计算p u t t i 和p a n i e o n i 1 9 9 5 r r j - 比了海水入 侵问题中耦合数学模型的两种数值迭代解法，s o r e k 等 2 0 0 1 1 给出了模拟可变密度地下水流的 二维区域模型。近期在多孔介质中可变密度水流和溶质运移耦合的理论和应用方面， s i m m o n s 等【1 9 9 9 】、b e a r 等 1 9 9 9 、d i e r s c h 和k o l d i t z 2 0 0 2 给出了很好的综述和评论。 2 第1 章绪论 以上的研究都属于理论上的方法性探讨，算例都很理想化，应用于野外的实际算例比较 少见。a n d e r s o n 等 1 9 8 8 将- - 维模型计算区简化理想的正杰面体，芳将边界条件简化处理。 g u p t a 和y a p a 1 9 8 2 在研究曼谷附近海水入侵时，采用的对流一弥散模型忽略了浓度变化对 流速的影响。自9 0 年代以来艟着计算机技术的极大发展和计算机的迅速普及，模拟海水入 侵的野外三维大尺度数值模型得到了广泛的应用。g h a s s e m i 等 1 9 9 6 对比了海水入侵问题的 2 维和3 维数值模拟结果，b i b a r 等 1 9 9 7 考虑了在浓度对密度影响下的巴塞罗那某海滨海水 入侵问题，n i s h i k a w a 1 9 9 7 综合评测了当时可用于海水入侵问题数值模拟的各种软件。b u e s 和o l t e a n 2 0 0 0 对比了分别利用混台杂交有限元法和不连续有限元法数值模拟海水入侵问题 的结果，g u r a n a s e n 等 2 0 0 0 数值模拟了f l o r i d a 州区域尺度下的海水入侵问题，e s s i n k 2 0 0 1 】 利用改编的特征法m o c d e n s 3 d 程序对荷兰一滨海含水层海平面上升所引起的海水入侵问 题进行了数值研究，a l i e w i 等 2 0 0 1 、s h e r i f 和h a m z a 2 0 0 1 数值模拟了抽水情况下海水运 移问题，p a n i c o n i 等 2 0 0 1 模拟和分析了突尼斯一滨海含水层的海水入侵问题。e s s i n k 2 0 0 2 】 利用m o c d e n s 3 d 程序模拟了意大利一滨海含水层在设置了虚拟的生态环保型地下阻水坝 后的海水入侵问题。 随着人们认识和实践水平的提高，人们逐渐发现数值模拟结果的精确度主要取决于对模 拟地区各类水文地质参数值的了解程度，然而现实世界中此类物理参数永远无法精确测得。 这就需要人们寻找一种即高效又精确的最优化算法，以便从部分实测值反演求得各类模型参 数值，遗传算法正是其中的一种方法。由于遗传算法的多项优点，其被广泛地应用于水文地 质学科的各个领域，如地下水管理问题( m c k i r m e y 和l i n , 1 9 9 4 ) 、多目标地下水污染问题 ( r i t z e l 等，1 9 9 4 ) 、多目标地下水监测问题( c i e n i a w s k i 等，1 9 9 5 ) 、抽水井的晟优开采问题 ( j o h n s o n 和r o g e r s ，1 9 9 5 ) 、水文地质参数识别问题( p r a s a d 穆r a s t o g i ，2 0 0 1 ) 等。 1 2 2 国内海水入侵相关问题的研究现状 韩再生 1 9 8 s 对秦皇岛市北戴河地区海水入侵及其防治方案进行了数值模拟。范家爵 f 1 9 8 8 忽略密度对水位的影响建立了一个二维平面有限差分模型，对大连市大魏家一水源 地海水入侵问题进行了数值模拟。薛禹群等 1 9 9 1 ，1 9 9 2 ，1 9 9 3 ，1 9 9 7 构建了海水入侵数学 模型并对龙口市、菜州市、招远县海水入侵问题进行了数值模拟。尹泽生等 1 9 9 2 对莱州市 滨海区域海水入侵的综合防治开展了研究，创造性的提出了海水入侵通道类型的划分原则。 艾康洪 1 9 9 3 考虑密度对水位的影响，采用上游加权有限元法建立了一个= 维剖面的数 值模型，研究了广西漫尾岛咸淡水过渡带的运移规律。吴吉春【1 9 9 q 首次建立了反映含水 层中水一盐阳离子交换作用的海水入侵数学模型，并将该模型成功的应用于山东龙口市海水 入侵的数值模拟。吴吉春等 1 9 9 6 提出一种改进特征有限元法，并证明该法能有效的求解描 述海水入侵问题方程组中的非线性问题。 周训等 1 9 9 7 对广西北海市海水入侵的分布地点、出现时间、变化范围和主要监测井的 氯离子含量进行了分析，提出了防治的对策。张永祥等 1 9 9 7 对莱州湾南部海水入侵的水文 地质和水化学特征进行了调查研究。庄振业等 1 9 9 9 针对山东莱州湾沿岸海水入侵灾害的发 展进程展开了研究，首次给出了海水入侵灾害的发展模式，并根据灾害发展的不同阶段提出 了相应的防、减灾措藏。 3 第1 章绪论 寇英杰和王颖 2 0 0 0 综台定性的分析了辽宁省海水入侵的分布、成因和防治措施。马凤 山和蔡祖煌 ：z o o o 在综合国内外当前控制海水入侵灾害技术的基础上，提出了适台于我国国 情的海水入侵综合防治应用技术。王国利和周惠成 2 0 0 0 以大连市海水入侵为例，分别建立 了沿海边界节点的单位脉冲一海水入侵响应矩阵和模型内部节点的单位脉冲一水位降深响 应矩阵，并使之与地下水流模型相结合，构造了一个地下水系统优化管理模型。 陈鸿汉等 2 0 0 0 ，2 0 0 2 1 分析了潍河下游海水入侵水文地质水化学特征，并进行了动态的 三维数值模拟。成建梅等 2 0 0 q 建立了一个考虑变密度的海水入侵数学模型，得出了山东省 烟台市夹河中、下游地区咸淡水界面的运移规律。袁益让等 2 0 m 提出了一种迎风分数步差 分格式，并将其应用于山东省莱州湾地区海水入侵防治工程后效的数值模拟。武强等 2 0 0 l 】 首次应用g i s 对秦皇岛市洋戴河地区海水入侵进行了数值模拟。c h e n 和c h e n g 2 0 0 1 利用三 维数学模型对多层含水层可变密度的海水入侵问题进行了数值模拟，刘艾礼等 2 0 0 2 给出了 秦皇岛市沿海地区海水入侵灾害的最新监测结果，并对该地区海水入侵今后的进一步研究提 出了思路。 此外，有关遗传算法在水文地质学中的应用研究方面，邵景力等 1 9 9 8 首次将遗传算法 应用于求解地下水资源管理模型。李竞生和姚磊华 2 0 0 l ，2 0 0 2 结合有限单元法，将遗传算法 应用于水文地质参数的反演识别。刘强等 2 0 0 2 ) 针对传统求解地下水资源管理模型的缺点， 采用遗传算法来求解，茹将其应用于松原市区的地下水资源管理。 1 3 本文的主要工作 客观地说，很多人在秦皇岛洋、戴河平原区海水入侵的数值模拟问题上作出了杰出的 工作，但是仍存在或多或少的不足。本文充分借鉴前人研究成果，从模型的概化到防治方案 的数值模拟，开展了以下较有创新性的工作： i ) 建立了一个考虑密度对水位影响的二维平面潜水非稳定流模型； 2 ) 充分扩展了模拟的区域。概化出一个相对独立的地下水系统，既能模拟出留守营、 樊各庄地区地下水漏斗的演交过程，又能反映滨海地区成淡水界面的运移规律； 3 ) 在数值模型中首次考虑了地表水体一河流与地下水系统的相互影响和作用； 4 ) 利用g i s 地理信息系统便捷的数据管理和直观的空间查询特点，构建了一套包括本 模型所需的所有输入数据的地理信息系统； 5 ) 针对各类水文地质参数值难以确定的问题，本文以一个简单的算例证明了遗传算法 在水文地质参数反演中具有广阔的应用前景。 4 第2 章海水入侵数学模型的理论研究 第2 章海水入侵数学模型的理论研究 2 1 地下水状态方程 众所周知，无论是地下水运动方程和连续方程的推导计算，还是地下水溶质运移方程的 推导计算，都离不开对于整个地下水系统中各类状态变量的求导。为了以后推导各类方程的 叙述简便，以下将依次说明流体密度p 7 、孔隙率、含水层厚度b 和动力粘滞系数7 这 几个状态变量的本构关系。 2 1 1 流体密度p 7 的本构方程 一般来说，流体密度p 7 被视作水力水头 、浓度c 和温度丁的函数： p 7 = p 7 ( ，c ，r ) ( 2 1 ) 对p 7 求全导数可得： 如，：型砌+ 芝犯+ 望刃 a犯a t 犯+ 防铡肿 、- - - - u - - - 一 一声 = 乃7 砌+ 瓦詈雨p 7 犯一劫7 刀 ( 2 - 2 ) 其中，为水压作用下的流体体积压缩系数，口为通过最大浓度e 规范化的流体密度差 率，可近似定义为云= 印r 7 p 7 = l d 7 ( e ) 一彳】“。万为热作用下流体的体积膨胀系 数。假设，、口和都是常数，则对( 2 - 2 ) 式定积分可得流体密度本构方程为： p ，：彳p 卅一卜面哥。 耐【1 + 砸也) + j ( c , - l c o ) ( c c o 卜- ( 卜t o ) j 眈- 3 ) 其中彳为考虑参考水头、参考浓度c o 和参考温度t o 的参考密度。 为了应用的简便，b a x t e r w a l l a c e 给出了对应于最大浓度e 的流体密度p 7 经验式： p j zp 七0 7 cs(2-4) 由此我们可得密度差率苫的经验估计式为： 第2 章海水入侵数学模型的理论研究 苫0 7 c , ( 2 - 5 )口) “ d i e r s c h 【2 0 0 2 】证明密度差率云与最大浓度e 之间的线性近似足以满足一般应用的精度 要求。 2 1 2 孔隙率占的本构方程 改变多孔介质的孔隙率，必然引起土体骨架的压缩变形。而由土力学知识我们知道， 土体骨架的压缩变形使得土体排出一部分地下水，致使土体中的孔隙水压力减小而土体有效 压力增加。自然的，孔隙率只是水力水头h 的函数，而与浓度c 和温度r 无关： s=占(_lz)(2-6) 其中y 代表土体骨架的压缩系数- 捌珊= y ( 1 - c ) 砌 协， 2 1 3 含水层厚度b 的本构方程 对于潜水含水层，含水层厚度口只是水力水头h 的函数，而与浓度c 和温度丁无关： b = b ( 办) ( 2 8 ) 对( 2 - 8 ) 式求全导数得： d b ：旦竺d h ：劬( 2 - 9 ) a 矗 。 其中f 为一个无量纲的抑止因子( 0sf 1 ) ，它的物理意义指的是当含水层处于灌水 时期，每升高单位水头所能引起水面升高的高度，既含水层中可填充自由水韵孔隙空间所折 算出的高度。反之，当含水层处于抽水时期，每降低单位水头所能引起水面降低的高度，既 含水层中可抽出自由水的孔隙空间所折算出的高度。 2 1 4 流体动力粘滞系数，的本构方程 一般的。流体动力粘滞系数口7 只是浓度c 和温度r 的函数，而与水力水头h 无关： 7 = o ( c ，乃 ( 2 1 0 ) h a s 湖i z a d c h ( 1 9 8 8 ) 提出了考虑高浓度梯度影响的j 与浓度c 的经验公式： f ( c ) = f ( 1 + 1 g s a 7 4 1 矿2 + 4 4 5 巧3 ) ( 2 1 1 ) 其中g = c fd j 。 进一步的，p i n d e r ( 1 9 7 7 ) 提出了考虑温度影响的u 7 与温度丁的经验公式： 6 | i 舶 缸瓦 霎| 虹 导全求式击 对 第2 章海水入侵数学模型的理论研究 l ；! ! ：! 旦箜二旦：坚! ! 垄： u j l l ( 2 1 2 ) 其忙等a 综合考虑浓度c 和温度丁的影响，推得相对于参考粘滞系数面( c o ，t o ) 的经验公式： ! ：! 垒：! 竺! ! = ! ：坚! ! 笙： 7 ( c ，t ) 1 + 1 8 5 a r - 4 1 n r 2 + 4 4 5 盯 进而可得相对于参考浓度c 。和参考温度瓦的粘滞性函数无； 2 志 ， 1 + 1 8 s 搿f c 喝) 一4 1 吨喝) + 4 4 5 口蟊喝) 2 i i 西云i 面i 磊广 2 2 地下水水流方程和运动方程 14 - 0 7 0 6 3 - 一o 0 4 8 3 2 - 。 x - 。_ 。_ _ 。_ 。_ _ _ _ 。_ 。一 1 + o 7 0 6 3 f f t r 喝) 一o 0 4 8 3 2 e ( r = r o ) ( 2 1 3 ) ( 2 ，1 4 ) 传统的三维非均质各向异性饱和含水层中的地下水水流方程和运动方程为： 瓯警+ 譬呜 ( 2 _ 1 5 ) ”8 ta x 9 吼f ：一k , f p f j g _ a l l = - k f 掣 ( 2 - 1 6 ) 吼一广i”i 屹。 其中h = p 7 ( 户7 9 ) + b 为该点处的莘际水头，p 7 为该点处的流体压强，p 7 为该 点处流体的实际密度，7 为该点处流体的实际动力粘滞系数，g 为重力加速度，k ：为该 点处的渗透系数张量，世，为该点处相对于实际的流体密度和实际的流体动力粘滞系数的水 力传导率张量，s 。为比贮水系数或称贮水率或压缩系数。q 。为单位体积多孔介质中流体的 源汇流量，g ? 为地下水d c y 流速，t 为时问，z 。为笛卡儿坐标，f ，歹= 1 , 2 , 3 遵循e i n s t e i n 求和约定。 但是方程( 2 1 5 ) 、( 2 1 6 ) 不能反应流体密度变化和动力粘滞系数变化对水流运动所造 成的影响。众所周知，在极高浓度海水入侵的含水层，尤其是在咸淡水过渡带附近，地下水 的浓度变化很大即浓度梯度很大，溶液浓度的变化必然引起溶液密度的变化。在这些高浓度 梯度的地下水中，压强相同的矗，其水头不一定相等，水头梯度的负值( 一v h ) 已经不再 反映地下水的流向。浓度变化造成垂向密度有差异，在重力作用下会产生水体间的自然对流。 为了考虑海水入侵含水层中流体密度变化对水流运动的影响，必须重新推导地下水水流方程 和地下水运动方程 7 第2 章海水入侵数学模型的理论研究 d i 目 s c h ( 1 9 8 5 ，2 0 0 2 ) 推导出地下水系统守恒方程的最一般式； 詈c p 厶+ 毒c 删+ 善= 协 该式可分别适用于质量守恒、动量守恒和能量守恒，；f ，、，和s 在各类守恒方程中所 对应的物理量如表2 - 1 所示。其中p 为广义密度，e 为实际流速而非d a r c y 流速。 表2 - 1 微观守恒方程物理量表 通过一系列简化，d i e r s c h ( j 9 3 5 ) 导出各类物理量各自独立的夺恒方裎： 流体质量守恒方程 鲁( 矽) + 砉( 妒，v j ) = 句d 7 睇f 流体动量守恒方程 ”芳降如，卜 2 2 1 地下水流运动方程 ( 2 一1 8 ) ( 2 1 9 ) 地下水流运动方程本质上就是地下水动量守恒方程，即d a r c y 方程。比较( 2 - 1 6 ) 式和 ( 2 1 9 ) 式。初看之下似乎两者并不相同。实际上两者完全一致，只是( 2 1 6 ) 式中想求的 未知量是各点处的实际水头日，而( 2 1 9 ) 式中待求未知量为各点处所对应的压强p 7 ，只 需将= p ( p 7 9 ) + 毛代入( 2 1 9 ) 式即可得到( 2 一t 6 ) 式。但是我们前面已经提到， 对于极高浓度海水既卤水所入侵的含水层，已经不能单纯用所求得的实际水头目来判断地 下水的流向，压强p ，自然也不能用于判断。因为此时的地下水流流向已经不仅仅由水头梯 度这一个驱动力决定，它还受到浓度变化引发的密度梯度这个驱动力的影响。如此一来就产 生一个问题，究竟如何考虑这个由浓度梯度引发的密度梯度驱动力对地下水流速v ，的影 响? h 瑟跚啦划e h ( 1 9 8 6 b ) 推导了高浓度梯度对流体运动方程的影响： 酽一等( 等+ p s g e 3 一薏一髫考 c z 其中吼为d a m y 流速，口= 以p 为流体中溶质口的质量百分数，为流体中溶质 口的密度( 数值上等于溶质口的浓度) ，彤和万分别是流体中高浓度梯度和高温度梯度 8 第2 章海水入侵数学模型的理论研究 的交叉作用系数，r 代表流体中温度，e 3 为与重力方向相反的单位矢量，即右手螺旋坐标 系中大拇指= 方向的单位矢量。 对于海水入侵的地下水流等变密度的水流运动，实际的水力水头片= p 7 ，( p g ) + 黾 不仅仅随着测点的压强和高程变化而变化，而且还随着测点的实际密度p 7 变化而变化。因 此，对于有着相同的压强值和高程值不同的密度值的两个测点，实际的水力水头目是不同 的。为了使计算得的等水位线图能够反映地下水的真实流向，我们假设测压管中的流体密度 和实测点的流体密度不同，所有点的测压管中的流体密度取个共同的参考密度“，既淡 水密度。则得参考水头h 如图2 - 1 所示： h = 与+ 墨( 2 - 2 1 ) 蒯 将式( 2 2 1 ) 代入式( 2 2 0 ) 得： 铲一竽浯+ 筝寸彤考一形善 = 一号笋等 若+ 丛p o ：白 一珥等一形考p p 1 a x | j。1 。融j a x | = 喝厶( 考十乞笋屯 一彤等一彤考 c m z ， 其中k 。为相对于参考密度岛和参考动力粘滞系数j 的水力传导率，不同于式( 2 - 1 6 ) 中的k “，兀为粘滞性函数，反映了流体中实际的浓度、温度变化对于动力粘滞系数的影 响，其余符号同式( 2 - 2 0 ) 。 图2 - 1 参考水头和实际水头对比圈 9 p ，。 。f + x 3 p g 第2 章海水入侵数学模型的理论研究 由式( 2 - 2 2 ) 可以看出，d 孤靶流速由四部分驱动力构成。其中第一项表示参考水头h 的 梯度对1 ) a u r c y 流速的影响，第二项反映的是由于浓度梯度所产生的密度差所表现出重力差异 引起的自然对流，第三项单独反映浓度梯度对于d a r c y 流速的影响，第四项反映温度梯度对 于d 哪流速的影响。 式( 2 - 2 2 ) 考虑了极高浓度梯度、温度梯度对于地下水运动方程的影响，如卤水入侵、 地热入侵等。但是对于一般海水入侵的含水层中，地下水中溶质浓度一般小于3 s g l ，地下 水密度小于1 0 2 5 9 l ，地下水中浓度梯度比较小，因而密度的变化幅度也较小( 一般小于5 ) ，同时海水入侵含水层中地下水的温度变化也不大。因此，我们可以忽略式( 2 2 2 ) 中 第三、四项对地下水流速的共轭影响。而且相对低浓度的海水对于地下水流动的阻滞切应力 较小，因而对于动力粘滞系数影响较小，故在海水入侵问题中可近似认为，= 1 。通过以 上一系列简化后可得： 铲喝浯+ 华岛 式( 2 - 2 3 ) 适用于反映地下水浓度接近或低于海水， 下水在含水层中的运动。 ( 1 - 2 3 ) 密度和动力粘滞系数变化较小的地 由于本文采用的二维数值模型，因此对式( 2 - 2 3 ) 作垂商平均后，即可得二维地下水运 动方程： 一q 。- - - b 巧睁筝白 q 睁筝勺 其中弓为地下水传导系数张量- 2 2 2 地下水水流方程 ( 2 2 4 ) 地下水水流方程实质上就是地下水的水流质量守恒方程，即式( 2 1 8 ) 。同理，由于本 文采用二维数值模型，故对其进行二维的垂向平均，可得： 昙( 妒7 占) + 未( 妒，b q ) = 印7 嘭 ( ：甾) 为简化起见令万= 出彰，i = 幽，。，可得： o ( e p r b ) + 毒( 户7 两= 户7 万 进一步展开( 2 - 2 6 ) 式得， j 0 ( 2 - 2 6 ) 星兰皇塑查堡墼兰蔓型塑里堕堡壅一 举+ 牮十筚+ 拳+ 擘万 ( 2 _ 2 7 ) 对于i 由密度p 的本构方程( 2 - 3 ) 得： 曲- e - ；却巾警+ 南詈一万别 对于2 由孔隙率e 的本构方程( 2 7 ) 得： 加詈酬t 叫詈 对于3 由含水层厚度b 的本构方程( 2 - 9 ) 得