NUMERICAL SIMULATION OF POLLUTANT MIGRATION IN MULTI-LAYER MEDIA WATER-BEARING SYSTEM CONTAINING AQUITARD: A CASE STUDY OF TIEBUTIE POOL, DAQING CITY
-
摘要: 选取大庆市贴不贴泡区具有普遍意义的饱和-非饱和多层介质含水系统中石油类污染为研究对象.在分析水文地质条件的基础上, 利用实验研究获得含水系统中连接上下含水层的弱透水层的渗透规律, 然后运用数值模拟技术建立起此类含水系统中石油类污染质运移数值模拟模型.应用所建模型对污染质的污染趋势进行了预测, 并提出污染控制措施.Abstract: The typical oil pollution in the saturated-unsaturated multi-layer media water-bearing system characteristic of Tiebutie Pool in Daqing City is chosen as the object of the present investigation. On the base of the analysis of the hydrogeological data, the permeability law of the aquitard linking the upper and lower water-bearing beds in the water-bearing system has been achieved by an experimental research. Then the numerical simulation model of the oil pollutant migration in this kind of water-bearing system is constructed by means of the numerical simulation technology. The numerical model thus established can be used to forecast the pollution trend of the pollutants and to propose the corresponding measures to control the pollution. In summary, this research result is of wide application in China where the sustainable development strategy is strongly upheld.
-
Key words:
- aquitard /
- multilayer media /
- waterbearing system /
- numerical simulation /
- pollution
-
具有弱透水层的多层介质含水系统在自然界中分布广泛, 如我国的松嫩平原、松辽平原内都分布有这类地下水系统, 而一些大型油田, 如大庆油田、辽河油田等恰恰分布于其上, 这类地下水系统往往遭受来自油田开发区石油类污染质的污染.弱透水层既是下部承压含水层中优质地下水的补给通道, 又是污染通道.因此, 本文针对这一实际问题, 选取大庆市的一个典型的纳污湖泡——贴不贴泡区作为研究区, 首先认识弱透水层的渗透规律, 然后建立起这类含水系统中污染质运移数值模拟模型, 并利用该模型对研究区地下水的石油类污染状况进行预测, 从而为这类地下水系统的污染的综合治理提供科学依据和有效措施, 确保油田的可持续发展.
1. 研究区概况
研究区位于石油城市大庆市喇嘛甸镇西北, 为一小型闭流盆地, 面积约12.28 km2, 东面与向荣屯、良种场接壤, 西面与沈家北泡毗邻, 南面与三棱泡交界, 北侧与北十里相连(图 1).区内贴不贴泡呈NW—SE向分布, 纵向长约3 km, 横向宽约1 km, 四周被风沙、岗地围成盆状.地层自下而上依次是: 第三系泰康组和第四系白土山组砂砾岩和砂砾石承压含水层, 厚约80 m; 第四系荒山组饱水粘性土层(弱透水层), 厚约51 m; 第四系哈尔滨组的细砂、亚粘土、亚砂土层, 厚13~19 m.区内贴不贴泡已成为石油化工业的废水排放处, 石油类质量浓度在排污口处为30.60 mg/L, 湖泡中为0.6 mg/L (枯水期).其下部潜水被污染的范围是以湖泡为中心, 向四周呈椭圆形状扩散, 构成污染晕.其中心的污染质质量浓度为0.6 mg/L, 向外逐渐降至0.03 mg/L.
第四系荒山组饱水粘性土层是全区的弱透水层.它连接着潜水和承压水, 成为地下水和污染质运移的通道.因此, 在建立污染质运移数值模拟模型之前, 先要对弱透水层的渗透规律进行研究.
2. 弱透水层的渗透规律
将饱水的原状粘性土试样置于渗压容器中, 利用渗压仪进行不同水头下的渗透实验, 获得V-I关系曲线, 即OABC折线(图 2), 以及V-I关系方程[1, 2].
线段OA:
(1) 线段AB:
(2) 线段BC:
(3) 式中: I为水力梯度, I01为第二阶段开始时的初始水力梯度, 其值为0.291, I02为第三阶段开始时的初始水力梯度, 其值为0.440, V01为第二阶段开始时的初始渗透速度, 其值为0.938×10-4m/d, V02为第三阶段开始时的初始渗透速度, 其值为1.938×10-4m/d, K1为第一阶段的渗透系数, K2为第二阶段的渗透系数, K3为第三阶段的渗透系数.
可见, 在较小的水力梯度I (0 < I≤I01) 作用下, 只有重力水在大孔隙中渗透, 其渗透规律符合达西定律, 并可用线段OA来描述(图 2, 式(1)).此时, K1=3.223×10-4m/d.在I01 < I≤I02条件控制下, 大孔隙中的重力水运动被加强, 微孔隙中一部分具有较小抗剪强度的结合水开始移动.该阶段渗透系数变化相对较快, 但整个阶段相对较短.V和I的关系实际由曲线AB来描述, 由于曲线AB的弯曲度不大, 故用直线AB来近似代替曲线AB.因此达西定律能近似地适用于此阶段的渗透(图 2, 式(2)).此时, K2=6.711×10-4m/d.在较大的水力梯度I (I > I02) 作用下, 不仅大孔隙通道中的渗透增强, 而且微孔隙通道中有更多的结合水也已运动.尽管水力梯度很大, 而且具有小抗剪强度的结合水已经运动, 但是对于具有大抗剪强度的结合水来讲, 其运动还是很困难的.此阶段, 渗透系数的增加很缓慢, 而且有随水力梯度的增加趋于常数的趋势, V和I之间的关系可近似地由线段BC来描述, 故可用达西定律来描述该阶段的渗透规律(图 2, 式(3)).此时, K3=19.954×10-4m/d.
3. 含水系统中污染质运移数值模拟
3.1 含水系统概念模型
研究区内含水介质自下而上依次为: 砂砾岩和砂砾石含水层, 饱水粘性土层(弱透水层), 细砂、亚粘土、亚砂土层.全区所有侧向边界均为第二类边界, 其中承压含水层侧向边界为已知流量边界, 其余为隔水边界.顶部边界湖泡水体为第一类边界, 其余为入渗补给边界.承压含水层底板为底部隔水边界.地下水在多孔介质中的运动符合达西定律, 污染质浓度变化符合Fick定律, 无化学反应.地下水运动为非稳定流, 含水介质为非均质各向异性.
3.2 数学模型
根据上述含水系统的概化, 可建立起描述地下水含水系统的数学模型.该模型由两部分组成: 水分运移模型和污染质运移模型[3~5].
(1) 三维流水分运移数学模型.为了将饱和带和非饱和带进行统一描述, 采用压力水头为状态变量.数学模型为
(4a) (4b) (4c) (4d) 式中:
(5) (6) Ks为饱和渗透系数, 在饱和弱透水粘性土层中其值由式(1~3) 确定.ψ (x, y, z, t) 为压力水头, μ*为贮水率, θ为体积含水率, α, nl为经验系数, 由非饱和带中水分特征曲线确定, m=1-1/nl, 且0 < m < 1, G为补给项(降水入渗等), W为排泄项(蒸发排泄、抽水等), ψ0 (x, y, z) 为初始压力水头, ψ1 (x, y, z, t) 为第一类边界上的压力水头, f1 (x, y, z, t) 为第二类边界上的水分通量, n为边界外法线向量, Γ1, Γ2分别为第一、第二类边界.
(2) 三维流污染质运移数学模型.三维饱和—非饱和污染质运移的定解问题为
(7a) (7b) (7c) (7d) 式中: C (x, y, z, t) 为污染质质量浓度, Rd为阻滞因子, D为弥散系数, 由纵向弥散度αl和横向弥散度αt以及渗流速度确定, Ce为污染源的质量浓度, C0 (x, y, z) 为初始污染质质量浓度, C1 (x, y, z, t) 为第一类边界上的污染质质量浓度, f2 (x, y, z, t) 为第二类边界上的污染质弥散通量.
3.3 模型识别与验证
先将其进行平面二维三角形剖分, 并进行编号, 共剖分为137个结点、226个三角单元.垂向上分6层, 7个层面分别为地表面、潜水面(湖泡底部为地下5 m)、地下10 m、弱透水层顶板、弱透水层顶板以下8 m、弱透水层底板和承压含水层底板.采用三棱柱单元输入信息, 然后计算机程序自动将计算区剖分成959个结点、4 068个四面体单元, 边界结点除湖水面上的结点为第一类边界点外, 其余的均为第二类边界点.然后采用特征线法与有限单元法相结合来求解上述数学模型.计算时的初始流场和初始污染质浓度场由专门布设的17个测压管式观测孔和3个民井的实测值, 经泛Kriging空间插值获得.
由于计算区面积比较小, 地层岩性平面上变化不大, 故将计算区平面上作为一个参数区, 而在垂向上分出4个参数区, 即由上至下剖分成的第1层、第2-3层、第4-5层和第6层.参数初值参考已有水文地质资料和室内外实验测定值给出.
由于弱透水层的渗透系数是水力梯度的函数, 因此, 实际计算时需要利用上述实验研究得出的弱透水层的渗透规律来解决这一问题.
以1995年11月20日实测地下水水头及石油类质量浓度为初始水动力场和初始石油类污染质质量浓度场, 拟合1996年7月9日的实测地下水水头和石油类污染质质量浓度.经反复拟合计算, 计算水头与实测水头之差的绝对值小于0.5 m的观测孔数占观测孔总数的90%, 石油类污染质质量浓度计算值与实测值之差的绝对值小于0.1 mg/L的观测孔数占观测孔总数的85%.最后得到的模型参数见表 1.用1996年11月20日的实测数据进行了验证, 证实所建模型可靠.
表 1 模型参数Table Supplementary Table Parameters in the numerical model4. 污染质运移预测及污染控制
运用上述所建模型即可对大庆市贴不贴泡区的多层介质含水系统中的石油类污染质的运移做出预测.计算时, 源、汇项和边界均按多年平均来处理, 非饱和带的导水率由(6a) 式确定, 弱透水层的渗透系数由(1) ~ (3) 式确定.结果表明, 在已有条件均保持不变的情况下, 地下水中石油类污染质的质量浓度 变化不大, 10 a后污染质锋面未穿透弱透水层, 但已污染了弱透水层4 m多.在目的层大量开采造成承压水头下降10 m的情况下, 污染质扩散范围有所加大, 但速度比较慢, 10 a后污染质进入弱透水层近8 m.而当开采目的层减少开采量, 并将潜水和承压水之间的水头差保持在14.55 m, 也就是降低了弱透水层中的水力梯度, 结果污染质的下移基本被控制.由此可见, 弱透水层中非达西流的存在对控制污染质下移起到了决定性的作用.
由预测结果可知, 在现状条件下, 石油类污染质仍有污染下部承压含水层优质地下水的趋势.若加大下部承压含水层的开采强度, 则石油类污染质下移速度会加快.因此, 要控制污染质下移, 一是要减少污染质的排放量, 降低湖泡中石油类污染质的含量; 二是要节约用水, 减少承压水的开采量, 以减小潜水与承压水之间的水头差, 利用弱透水层的性质来阻止污染质下移.当污水排放量保持现状, 要使石油类污染质不下移, 则潜水和承压水之间的水头差不得大于14.55 m.
5. 结论
(1) 由饱和粘性土构成的弱透水层, 在自然界中是普遍存在的, 对其渗透规律的研究是建立多层介质含水系统数值模拟模型的基础; 因此, 该研究具有重要的理论意义. (2) 应用数值模拟方法和弱透水层的渗透规律, 在石油开发区建立起一个包括包气带、潜水、弱透水层和承压水在内的多层介质含水系统中石油类污染质运移的数值模拟模型, 是目前解决此类地质条件复杂地区地下水污染问题的最为有效、实用的方法.它可为该类地区地下水污染控制措施的制定提供可靠而精确的科学依据. (3) 根据模型预测结果, 提出大庆市贴不贴泡区地下水系统中石油类污染质下移控制主要措施: 减少污染质排放量, 控制承压水开采.
-
表 1 模型参数
Table 1. Parameters in the numerical model
-
[1] 张忠胤. 关于结合水动力学问题[M ]. 北京: 地质出版社, 1980. [2] 宿青山, 刘丽, 邱枚大, 等. 现代实验水文地质[M ]. 长春: 吉林科学技术出版社, 1991. [3] 陈崇希, 唐仲华. 地下水流动问题数值方法[M ]. 武汉: 中国地质大学出版社, 1990.157~182. [4] 孙讷正. 地下水流的数学模型和数值方法[M ]. 北京: 地质出版社, 1981.214- 241. [5] 孙讷正. 地下水污染: 数学模型和数值方法[M ]. 北京: 地质出版社, 1989.123~261. -