天津市平原区地下水开采现状与治理对策

天津市平原区地下水开采现状与治理对策

0地面沉降耦合模型天津大量地下水资源，水流集中在漏水中心，形成了城市、塘沽区、汉沽区和海口河下游等沉降中心。地下水开采对地面沉降的影响相当严重，必须尽快解决。因此，我们必须研究地下水流-地面沉降耦合模型，更好地分析地下水资源的组成，预测地下水流场和地面沉降的变化趋势，为社会、经济和环境的可持续发展提供重要参考。基于理论法的地下水流-地面沉降数值模型包括水流模型与沉降模型2个部分,按结合方式分为两步走计算模型、部分耦合模型和完全耦合模型。地面沉降模型的计算最常用的是一维太沙基原理,同时还有Boit固结理论,考虑流变特性的流固耦合计算理论,以及考虑黏弹性、黏弹塑性变形的计算方法等。美国地调局在Modflow基础上,开发出ISB程序包,用于模拟水头变化而造成的地层压缩过程。1993年,崔小东与英国水文地质学家Kiching等人考虑了黏土层中的沉降滞后效应,共同开发了IDP程序包。在ISB的基础上,美国地调局最新研究的Sub程序包,能够模拟随时间变化的排水量以及夹层(弱透水层)滞后排水和压缩的过程,从而更加准确地确定地面沉降量。由于大区域地下水流-地面沉降耦合模型所涉及的因素复杂、所需区域地面沉降相关资料难以获得,国内外建立此类模型甚少。天津市早期有关地面沉降模型研究的成果主要有郭占荣建立的准三维水流和一维地层压缩耦合模型;崔小东利用IDP程序包,将其与Modflow耦合,对天津汉沽含水层系统进行了模拟;徐鸣建立了滨海新区地下水-地面沉降耦合模型。上述研究过程中没有充分考虑三维渗流特点以及黏性土层的厚度变化特点等因素,有一定的弊端。本次研究从Leake提出的“夹层”理论出发,运用Modflow2005以及Sub程序包,考虑到渗流方向、夹层厚度变化等特点,通过各个时段夹层边界的水头值将地面沉降模型耦合于水流模型中,建立了天津市平原区三维渗流与一维地面沉降耦合模型。1地下水开采量分析天津市平原区东南临渤海,北依燕山山脉,西北部为北京市和河北省,面积约为1.06×104km2。研究区概况见图1。区内主要含水层为第四系,自上而下可分为全新统天津组(Qh4)、晚更新统塘沽组(Qp3)、中更新统佟楼组(Qp2)、早更新统杨柳青组(Qp1),新近系明化镇组上段含水层相对发育,在本区也作为深层承压水加以开发利用。自山麓至渤海之滨的广大平原区大致以宝坻断裂为界,北面分布有第四系孔隙地下水与下伏基岩地下水,南面分布有巨厚新生界,赋存孔隙地下水。本区地下水多年处于超量开采状态。据统计,近20多年来地下水年开采量为(6.0~11.5)×108m3,地下水位呈持续下降趋势,产生了严重的地面沉降问题。宝坻断裂以南的广大平原区均有不同程度的下沉,截至到2002年,地面沉降面积达到8798.12km2,其中累计沉降量超过1m的面积达到4080.48km2,沉降区域已与邻近的河北省地面沉降区连成一片,构成华北平原地面沉降区的一部分。2地下水流与地面沉降的结合模型2.1理论模型及边界条件根据研究区水文地质结构特征,结合地下水开发利用状况,将研究区概化为6个含水层组:第I层为潜水含水层组;第Ⅱ—Ⅵ层为承压含水层组。考虑到第Ⅶ含水层组(明化镇组)与其上的含水层组的水量交换,将其作为底部越流层。在边界处理上,除了第I含水层沿着渤海湾边界为给定水头边界(水头值为0m)及第Ⅰ、Ⅱ含水层北部是已知流量边界外,其余研究区边界(各层)为通用水头边界(general-headboundary,GHB)。水文地质参数是反映源汇项分布以及水文地质结构的重要数据,主要有大气降水入渗系数(α)、含水层的渗透系数(K)、潜水含水层给水度(μ)及各承压含水层的弹性储水率(Ss)。水文地质分区主要参考岩性、富水性分区,各区的水文地质参数初值主要参考相关的水文地质报告(1),在综合分析水文地质条件的基础上,给出各水文地质参数分区和初值,为模型调试的依据。其中,各层骨架弹性储水率一为6×10-61/m、骨架非弹性储水率统一为6×10-41/m,模型参数及边界条件见表1。地下水开采是研究区地下水的主要排泄项之一,主要是农业、生活、工业以及林牧副渔业用水。此外,大气降水入渗、农田灌溉入渗、蒸发、侧向流入流出、压缩释水也是研究区的主要源汇项。概念模型示意图见图2。由于潜水含水层补给丰富,地下水恢复较快,水位年内变差较小,因此,在本次研究中不考虑潜水层的压缩释水。而研究区深层地下水主要依靠越流补给,在人工开采条件下,水位变化幅度较大,且呈现多年下降趋势,因此将第Ⅱ-Ⅵ含水层作为压缩层。本研究仅考虑由于抽取地下水引起的土体变形。根据有效应力原理,土体变形可用有效应力原理来解释,分为砂性土的弹性变形和黏性土的塑性变形。根据Biot理论,当黏性土厚度较大时,黏性土的释水和变形存在着明显的滞后现象,本研究将黏性土层≤1.5m的作为非滞后黏性土层,将>1.5m的作为滞后性黏性土层。非滞后黏性土层与含水层砂层统一概化为非滞后压缩层,而在同一个含水层中的滞后性黏性土层作为滞后性压缩层,采用双面排水扩散方程处理。2.2数学模型2.2.1非均质流系统研究区水位随时间变化,水文地质参数随空间变化,参数在水平方向上无明显的方向性,因此地下水模型可以概化为三维、非均质水平各向同性的非稳定地下水流系统。在笛卡尔坐标系下可描述为式中:h为水头(m);K为含水介质的水平渗透系数(m/d);Kz为含水介质垂向渗透系数(m/d);ε为含水层的源汇项(1/d);Ss为自由面以下含水层储水率(1/m);t为时间(d)。2.2.2各含水层的水土流失地面沉降模型中,由水位下降引起的含水层土体变形量Δb为其中:式中:Ssk为骨架储水率(1/m);Sske为弹性骨架储水率(1/m);b为压缩层的初始厚度(m);Δh为水头变化值(m);Sskv为非弹性骨架储水率(1/m);h′为当前时段计算的水头(m);h′min为和前期固结应力对应的前期最低水头(m)。非滞后夹层对于厚度较薄(≤1.5m)的黏性土夹层,地面沉降不考虑滞后作用。此时,从非滞后压缩层中释放的水量为式中:q1为单位体积内弱透水层释放水量(1/d);γ为黏性土厚度占含水层厚度的比例。滞后夹层对于厚度较大(>1.5m)的滞后性黏性土夹层,其固结过程概化为一个双面排水过程,可用下式来表示:式中:z为垂向距离(m);S′s为黏性土层储水率(1/m);K′z为黏性土层垂向渗透系数(m/d)。在同一个模拟层中,将具有相同垂向渗透系数和相同弹性、非弹性储水率的滞后夹层组合成一个滞后夹层系统。等价夹层和含水层通过2个界面(上部边界和下部边界)交换的水量为则通过夹层(部分)系统的所有夹层交换的水的体积为单元格j中每一个滞后夹层的压缩量计算公式为(6)—(8)式中:Qj为第j个节点处通过两界面交换的水量(m3/d);Δxj、Δyj分别为x、y方向单元格的大小(m);Δz为夹层厚度(m);hjm为m时段末时刻含水层中单元格j的水头,即夹层上下边界处的水头(m);Qjtotal为通过夹层系统交换的总水量(m3/d);nequiv为等效比例因子;Δbjm为第m个应力期网格点j处的沉降量(m);S′mski、S′mskN分别为第m个应力期网格点i、N处的储水率(1/m);h′im、h′Nm分别为第m个应力期网格点i、N处的水头(m);h′im-1、h′Nm-1分别为第m-1个应力期网格点i、N处的水头(m);H′im、H′im-1分别为第m、第m-1个应力期网格点i处的预固结水头(m)。将上述方程组合,便可得到地下水-地面沉降耦合模型,通过各个时刻夹层边界的水头,将地面沉降模型耦合于水流方程中,夹层中释放出来的水量最终叠加到水流方程中的源汇项中进行计算。3地下水流梯度模型3.1计算结果与分析研究区面积约1.1×104km2。用Modflow2005程序对地下水流-地面沉降进行模拟,采用有限差分法,进行矩形剖分,单元长500m,宽500m,共剖分了7层,各层346行248列,单元总数为600656个,其中活动单元格276336个。由于大区域模型的复杂性,模型在最初一段应力期是一个自我调整过程,计算结果与实际观测数据之间误差较大。因此,长时间序列对模型的精度起着至关重要的作用,本次模拟选定1998年1月到2008年12月为模型识别验证期,以一个月为一个应力期,共132个应力期,每个应力期分为2个时间步长。3.2数据库的建立以1997年12月底地下水水位资料为基础,采用克里金插值法获取第Ⅰ—Ⅵ含水层的初始水位。为了模型的拟合需要,将初始沉降量设为0m,模型中含水层最低水位通过观测孔的历史水位数据得到。水文地质参数的识别是模型主要任务之一,将各项参数分区初值输入到模型中,通过拟合各个含水层流场和典型孔的动态过程曲线,识别含水层的参数,最后确定各参数分区值。研究区源汇项可以分为点、线、面3种要素。本次研究利用GIS软件的空间分析功能,统一采用数据库管理模式存储随时间变化的大量数据,最后利用程序写入到相对应的Modflow格式文件中。将源汇项以乡镇为单位,概化为井的形式,分配到研究区的各个网格点中。根据钻孔资料统计的含水层的砂层、黏性土层层数以及对应每一层的厚度,分别求解出每一个观测孔的滞后、非滞后压缩层相对应的参数,然后利用地质统计学中克里金插值方法得到每一个模型网格的参数数据,最后形成Sub文件。3.3模型识别验证3.3.1地面沉降模型拟合本次模型识别与检验采用试估-校正法,即反复地调整参数和某些源汇项,达到较为理想的拟合结果。本次侧重于对深层地下水研究,深层地下水水位观测孔数较多,第Ⅱ—Ⅵ含水层组参与拟合的水位长观孔分别有136、87、53、42四个,地面沉降分层标共有6组参与拟合,可进行区域流场及历时曲线的拟合(图3—5),从总均衡上控制承压含水层的区域流场。计算出来的地面沉降中心基本上与实际吻合,沉降标处的沉降过程线拟合良好。产生误差的主要来源是各源汇项的统计误差、地层结构概化误差、水文地质参数空间分布的差异性以及边界条件的处理等。此外,一个重要因素是模拟范围大,源汇项众多,存在很多概化问题;特别是本次开采量以乡镇作为统计单元,即开采量平均分配到每一个乡镇的网格中,不可避免地导致一定的误差。模型网格大小为500m×500m,由网格中心水位表征250000m2区域上水位变化情况,当进行单点动态曲线拟合、水位波动较大时,仍会有较大的误差。本次地面沉降模型是在拟合校正后的地下水流模拟基础上进行的,由于市区流场复杂,拟合效果不是很理想,市区附近水位出现明显的下降,直接导致了市区大直沽计算的地面沉降与实际有所偏差;由此可知,只有水流模拟保持了较高的精度,才能保证地面沉降拟合的精度。3.3.2地下水的补给、径流、边界进出情况通过建立1998—2008天津市平原区地下水流-地面沉降耦合模型,得到天津市平原区地下水均衡项,以多年平均的均衡项进行分析,如表2所示。在本次研究中,模型共分为7个含水层组,深层地下水主要考虑第Ⅱ—Ⅵ含水层组,第Ⅶ含水层组只作为越流层考虑,因此,表1中的均衡情况没有将第Ⅶ含水层组计算在内,仅将其与第Ⅵ含水层组发生的越流交换量统计到深层含水层组均衡中。从表2可知:天津市多年平均地下水均衡中,总补给量为16.96×108m3/a,主要的补给资源来自大气降水入渗(占补给量的54.01%)、农业灌溉入渗(占14.15%)、侧向边界流入(13.97%),此外还有来自第Ⅶ含水层组少量越流补给(仅占补给量的0.65%)。总排泄量为18.33×108m3/a,其中:人工开采为8.76×108m3/a,占总排泄量的47.79%;蒸发为9.00×108m3/a,占总排泄量的49.10%。总补排差为-1.37×108m3/a,地下水系统处于负均衡状态。在深层地下水补给量中,越流补给为3.41×108m3/a(包括来自浅层和第Ⅶ含水层组越流补给),占深层含水层补给量41.84%;压缩释水占32.15%;侧向边界流入占24.17%。由此可以看出,在多年开采条件下,越流补给、压缩释水、侧向边界流入已经成为深层地下水补给的主要方式。深层地下水主要消耗于人工开采,其开采量占总开采量的89.61%。对于浅层水而言,水位保持在一定的范围内,地下水由蓟县山前接受补给,向中部平原、滨海平原径流。蓟县、宝坻、武清处水位较高,水流通过其边界侧向流出一部分到河北和北京相邻区域。对于深层水,在南部边界附近,水位高于周边,部分水侧向流到河北省;冬季地下水开采量小,水位稍有回升,亦使得小部分水侧向流入相邻区域。4天津市平原区地面沉降趋势预测南水北调实施后,将大幅减小深层地下水开采量,以减缓地面沉降下降速率,达到地下水涵养与保护的目的。本次研究以2010年为起始年,以南水北调规划开采量为基础,对天津市平原区地面沉降发展趋势进行了预测。预测结果(图6)显示:在对天津市区实施全部停采措施后,随着地下水位回升,市区大直沽地面沉降有所恢复(图中负值代表地面沉降恢复值);在汉沽D1点,预测初期仍有地面沉降发生,2012年沉降量达到最大,为20mm