地面沉降对含水层参数及给水能力的影响研究

刘 蓉，曹国亮，赵 勇，陆垂裕，孙青言，严聆嘉，彭 鹏

(1.中国水利水电科学研究院流域水循环模拟与调控国家重点实验室，北京 100038；2.河北省水文水资源勘测局，河北 石家庄 050031)

地面沉降已经成为一个全球性的地质灾害问题，导致地面沉降的因素很多，过度开采地下水是造成地面沉降的主要原因，地面沉降的沉降速率和沉降量与地下水位的下降速率和下降量密切相关[1]。目前我国形成的面积较大的地下水漏斗区多达56个，漏斗中心地下水埋深最大达百余米，并且各个区域的沉降中心仍在不断发展，沉降空间分布范围有逐渐连成一片的趋势[2]。由于深层地下水开采量很大一部分来自弱透水层的压密释水，这部分变形为不可恢复的非弹性变形，因此弱透水层压密释水量是一次性消耗量[3-5]。

土层压缩变形会导致土层相关水文地质参数改变，其中释水系数和渗透系数是关键影响参数，现有的地面沉降模拟中多将参数设置为固定值[6]，造成很大的误差[7]。随土层的压缩变形，土壤孔隙度逐渐减小，渗透系数、储水系数也随之减小，不同时期不同深度土层的渗透系数和储水系数都不同[8-9]，因此考虑时变性参数的含水层系统压缩模拟已成为地面沉降研究的关键和趋势。地面沉降过程中水文地质参数的变化规律是水文地质研究的一个难点，现有少量研究也基本是通过室内岩土压缩实验分析参数变化规律[10-15]。除了使地质参数发生改变外，土层压缩变形带来的更为严重的影响是对含水层系统给水能力和储水能力的影响，弱透水层非弹性压缩变形造成的含水层系统储水能力的永久性减小对地下水资源保障能力的影响尤为重要。

本文以含水层压缩过程中的物理机制为依据，建立一维压缩模型，并以沧州地区为研究区进行沉降计算，旨在通过数学手段研究地下水开采量、沉降量、地质参数变化以及给水能力变化之间的关系。浅层的越流补给量是地下水开采资源的重要组成部分[16]，但关于越流补给量所占的比例却一直存在争议，本文同时探讨了越流补给规律。

1 研究方法

1.1 水文地质概念模型

沧州是我国华北平原最典型的沉降漏斗区，沉降中心地区的累计沉降量已达2 600 mm[17]，近几年由于限制地下水的开采，沉降速率得以减缓但仍在持续发生。为了研究地面沉降对于含水层系统地质参数变化和给水能力变化的影响，本文将沧州地面沉降区的地层结构[17]槪化为一个多层含水层系统(图1)，并对地面沉降概念模型作如下假设：①承压含水层系统始终为饱和系统；②承压含水层地下水开采量主要由承压含水层弹性释水、弱透水层非弹性释水以及潜水含水层向承压含水层的越流补给量构成，不考虑承压含水层的侧向渗流；③潜水含水层水位保持不变；④土层的形变为均匀形变；⑤不考虑弱透水层压缩变形的时间滞后性；⑥不考虑岩土固体部分和水本身的压缩，土层的压缩变形全为多孔介质孔隙度的减小。

图1 多层含水层系统地层结构及弱透水层压缩变形概念模型示意图Fig.1 Schematic diagram of the structure of themulti-layer aquifer system and the conceptualmodel of the aquitard compression deformation

根据太沙基有效应力原理，作用于含水层系统上的上覆岩土体、地表建筑物、大气压力等荷载形成的总压应力σ由有效应力σ′和孔隙水应力μ与之动态平衡。当地下水开采后，地下水位下降导致孔隙水应力μ减小，总压应力σ保持不变，有效应力σ′增大，从而压缩多孔介质，引起土层压缩变形，地面发生沉降。

深层地下水的开采量主要由以下几部分构成：承压含水层的弹性释水量、弱透水层的非弹性释水量、越流补给量、侧向补给量。对于沧州地区，张兆吉等[18]根据水均衡计算其深层地下水资源由8.9%的弹性释水量、41.1%的黏土释水量构成、46.7%的越流补给量和3.3%的侧向补给量构成；郭永海等[19]计算的沧州深层地下水开采量由24.4%的砂层压密释水量、36.45%的黏土层压密释水量、36.19%的越流补给量和3%的侧向补给量构成。以上研究说明，研究区的侧向补给量非常小，本次研究不予与考虑。因此在本模型中，地下水开采量由承压含水层的弹性释水量、弱透水层的非弹性释水量以及潜水层越流补给量三部分构成。随着地下水位下降，含水层系统的给水能力受到影响，地下水可开采量、含水层的弹性释水量和非弹性释水量、越流补给量均发生变化。

1.2 含水层系统压缩变形数学模型

为进一步研究深层地下水开采过程中地面沉降对于含水层系统地质参数变化和给水能力变化的影响，本次研究在提出的概念模型基础上，结合已有的理论知识和经验公式[20-24]，建立一个一维数值模型对研究区进行模拟计算。

承压含水层地下水开采量：

V=Ve+Vv+Vl

(1)

式中：V——地下水开采量模数(单位面积的开采量)/m；

Ve——承压含水层的弹性释水量模数/m；

Vv——弱透水层的非弹性释水量模数/m；

Vl——来自潜水层的越流补给量模数/m。

有效应力：

σ′=σ-μ

(2)

假定总应力σ不变，所以有效应力的变化Δσ′即为孔隙水压力的变化：

Δσ′=γwΔh

(3)

式中：γw——水的容重；

Δh——承压水头的变化。

承压含水层弹性释水量：

Ve=SkeΔh

(4)

Ske=Sskeb

(5)

(6)

(7)

(8)

式中：Ske——弹性储水系数；

Sske——弹性储水率；

b——承压含水层厚度/m；

Cr——弹性压缩指数；

e——承压含水层砂性土的孔隙比；

Δb——承压含水层厚度变化/m；

b0——初始承压含水层厚度/m；

Δe——孔隙比变化；

e0——初始孔隙比。

弱透水层非弹性释水量：

Vv=SkvΔh

(9)

Skv=Sskvm

(10)

(11)

式中：Skv——黏土层的非弹性储水系数；

Sskv——非弹性储水率；

m——弱透水层厚度/m；

Cc——非弹性压缩指数；

e——黏性土的孔隙比。

弱透水层厚度变化Δm、黏性土的孔隙比变化Δe的计算公式与承压含水层相同。

潜水层的越流补给量：

(12)

(13)

式中：K——渗透系数；

Δt——越流时间/d；

Hu——潜水水位；

Hc——承压水位；

C——与土体渗透系数相关的常数；

n——经验指数的参数，一般取5；

e——弱透水层黏性土的孔隙比。

黏土层越流计算假设越流和水位变化同步，实际由于黏土层的压缩滞后作用，存在越流滞后的现象，当黏土层相邻含水层出现水头差后，造成黏土层的压密释水，当相邻含水层的水头变化逐渐波及整个黏土层后，才开始发生越流，越流发生时间要滞后于相邻含水层中水头变化发生的时间。根据华北地区黏土层厚度、渗透系数和储水系数计算的黏土层越流滞后时间为3～5 a[25]，相比本次模拟的时间尺度来说较短，对计算过程中的水量平衡影响不大。

1.3 含水层系统压缩模拟计算

将承压水头的变化率Δh作为控制条件，采用迭代法进行计算。在初始状态下，给定一个承压水头变化值，带入模型计算出开采量和变化后的各项地质参数的值，然后将这些变化后的参数作为下一次模拟计算的初始条件进行计算，如此迭代循环计算直至计算结束。

模型参数分为两种：固定参数，见表1；迭代变化参数，其初始值见表2。所有参数的选取均结合研究区实际情况和经验值进行设置[20-22,26]。

表1 模拟过程中固定不变的参数

表2 模拟过程中变化参数及其初始值

沧州漏斗区地下水位埋深变化趋势如图2所示[26-27]，深层地下水位埋深累计下降约100 m，浅层地下水位埋深的变幅不大，基本维持在10 m左右。参考沧州漏斗区地下水位变化情况，为了更清晰剥离分析水文地质参数变化对含水层系统变形及开采量的影响，此次模拟将地下水位持续下降阶段按匀速下降处理，沧州市沉降中心地下漏斗形成过程中地下水位的平均下降速率为2.8 m/a，将地下水位的计算变化步长Δh定为0.1 m，对应时间Δt为13 d，模拟结束时刻将累计地下水位下降值设置为200 m；根据式(13)反推得到与土体渗透系数相关的常数C；有效应力σ′由总应力(主要是含水层的上覆土体荷载)减去孔隙水压力得到，假设总应力基本保持不变，σ′的变化等于孔隙水应力的变化；假设初始条件下承压水位与潜水位相同，其水头差Hu-Hc初始值为0，因为潜水位的变幅很小，所以假定潜水位不变，随着承压水位的降低，承压水位与潜水位间水力梯度逐渐变大。

图2 沧州漏斗区深层地下水(第III含水组)埋深年际变化趋势图Fig.2 Interannual change trend of the depth to the deepgroundwater in the groundwater level depressioncone area near Cangzhouu数据来源：根据韩占涛等[27]和张廉钧[28]研究结果所绘。

2 结果与分析

随着承压水位逐渐下降，含水层系统的多孔介质受到压缩，土层发生压缩变形，孔隙比减小，储水系数、渗透系数等相关水文地质参数发生变化，含水层系统的给水能力和储水能力也变化。

2.1 含水层系统压缩变形分析

图3表示的是承压水位逐渐下降的过程中，弱透水层和承压含水层的厚度变化情况。弱透水层的压缩变形量明显大于承压含水层，承压水头下降100 m时，弱透水层的变形量为4.1 m，承压含水层为1.4 m；假使承压水头累计下降达到200 m时，弱透水层的变形量为6 m，承压含水层为2 m。深层地下水的开采造成的地面沉降主要是由弱透水层的压缩变形造成的，且这部分压缩变形为非弹性压缩，不可恢复，是一种永久性的地质灾害。压缩释出的是弱透水层内地质历史时间尺度上形成的地下水，地下水资源储存空间减小的体积，对于深层地下水资源同样是一种永久性损害。承压水头匀速下降时，弱透水层和承压含水层的变形速率逐渐减小。

图3 弱透水层和承压含水组的土层厚度的变化Fig.3 Changes in thickness of the aquitardand confined aquifer

2.2 含水层系统水文地质参数变化分析

图4表示的是随着地下水位的下降，含水层系统相关水文地质参数的变化过程。其中，图4(a)表示的是弱透水层和承压含水层介质孔隙比的变化规律。孔隙比的变化规律与土层厚度变化规律基本相同，都是逐渐减小且减小的速度逐渐变慢，弱透水层的孔隙比减少量是承压含水层的3倍。图4(b)表示的是弱透水层的非弹性储水系数和承压含水层的弹性储水系数随水头下降的变化，二者都逐渐减小，且减小的速度越来越慢。从图中可看出地下水位下降的前期，储水系数的变化非常明显，当地下水位累计降深超过70 m后，二者的变化速率开始变得缓慢并开始逐渐趋于稳定。另外，当累计水头下降100 m时，弹性、非弹性储水系数都减少了约62%；当累计水头下降200 m时，弹性、非弹性储水系数都减少了约77%，即结束地面沉降发生后含水层系统储水系数仅为地面沉降发生前的20%～40%。含水层系统压缩，储水系数减小趋势显著，在地面沉降地区用固定不变的水文地质参数进行地下水资源量计算是极其不准确的。储水系数也表征含水层系统的给水能力，随着地下水位的不断下降，整个含水层系统的给水能力受到严重威胁，假设对深层地下水持续进行开采，地下水位继续下降，土层继续压缩，储水系数不断减小，最终可能导致土层压缩到极限，不再有非弹性压缩释水。图4(c)表示的黏土层的渗透系数随水头下降的变化，也是逐渐减小且减小的速度越来越慢。由式(13)可知，黏土层的渗透系数由介质的孔隙比决定，并且成正相关，但渗透系数的变化程度远大于孔隙比，当水头累计下降100m时，渗透系数减少37%，而弱透水层孔隙比仅减小了9.8%。

图4 弱透水层和承压含水层水力参数的变化Fig. 4 Changes in hydraulic parameters of the aquitard and confined aquifer

2.3 给水能力的变化分析

图5表示的是承压含水层的弹性释水量、弱透水层的非弹性释水量、浅层越流补给量以及总开采量在单位水头降深下的变化和累计变化，反映出随地下水位下降，含水层系统给水能力的变化规律。当水头累计下降100 m时，单位水头降深可开采的弹性、非弹性释水量减少62%，越流量增加将近3倍；当水头累计下降200 m时，单位水头降深可开采的弹性、非弹性释水量减少77%，越流量增加5倍多。随着土层的压缩变形，弹性释水量和非弹性释水量在大幅度减小但减小的速度越来越慢，这与弹性、非弹性储水系数的变化相吻合，都反映出弱透水层和承压含水层的给水能力随承压水位下降逐渐变小。而越流补给量却随水位的下降呈明显增大趋势，与弱透水层渗透系数的变化规律相反，这是因为渗透系数虽然在逐渐减小，但潜水层与承压含水层之间的水头差却在不断增大，所以越流补给量呈缓慢增长的趋势。由于弱透水层、承压含水层给水能力逐渐减小且减小的速度越来越缓慢并且趋于稳定，浅层越流补给量不断增加且后期增幅较小，所以整个含水层系统深层地下水的给水能力逐渐减小，前期减小速度较快后期逐渐趋于稳定。

图5 可开采量的变化Fig.5 Changes in the amount that can be exploited

图6表示的是单位水头降深可开采量中的弹性释水、非弹性释水、越流补给三部分的占比和累计可开采量中的各部分占比。图6(a)可看出，随着水位的不断下降，地下水位每下降1 m可开采量中，弹性释水量、非弹性释水量所占比例持续减小，累计水头降深达100 m时，弹性释水量由24.2%下降到12%，非弹性释水量由75.7%下降到38%，均下降了将近一半；而单位水头降深下的越流补给量所占比例却在不断增加，从原来的接近0增长到50%。图6(b)可看出随着地下水位的下降，弹性、非弹性释水量所占累计开采量的百分比在逐渐减小，而越流补给量所占的百分比在逐渐增大。当水头累计下降100 m时，弱透水层的累计非弹性压密释水量约占总开采量的55.9%，与前人[15]研究得到的57.6%相差不大。

对于深层地下水开采造成的沉降漏斗区，越流补给量已经成为深层地下水开采的重要组成部分。目前关于越流补给量的计算方法主要有三种：张兆吉等[17]根据地下水均衡测算出沧州平均每年的越流补给量能达到47.6%；石建省等[29]根据沉降体积占深层开采量的比例反推出京津冀德平原累计越流量与侧向补给量所占的比例为55%～58%；王莹等[30]根据黏土层中氟离子浓度的变化推算出越流补给量占深层地下水开采量的36%。本次计算得到的沧州目前越流补给量占深层地下水开采量的49%，未来随着承压水头的下降，越流补给量所占的比例还会继续增加。随压缩变形对弹性、非弹性释水量的影响，承压含水层的地下水开采量必须通过越流量的增加而维持。根据以上分析，深层地下水的开采主要由黏土层的非弹性释水构成，但随着地下水位的不断下降，土层逐渐被压实，这种压密释水逐渐减小，更多的开采量由来自浅层的越流补给量替代，越流补给量对于深层地下水开采量的贡献越来越大。

图6 各部分可开采量的占比Fig.6 Percentage of the recoverable amount of each part

3 讨论

储水系数是地下水流模型计算中的关键水文地质参数，其大小动态变化，且随着沉降的发展减小幅度很大。对于地下水资源量计算，使用固定不变的储水系数计算地下水资源量存在很大误差，容易忽略地下水资源量的亏空。对沉降漏斗区开展地质调查，确定新的地质参数，并用动态发展的计算方法计算地下水资源量对于准确评价目前的地下水资源量来说至关重要。

从给水能力的变化来看，随着地下水位的下降，整个含水层系统的给水能力逐渐变弱，虽然后期有微小的增长，但这部分增长的开采量由越流补给量贡献。对于本研究建立的理想的一维沉降模型来说，确保浅层水量充足、潜水位与承压水位差越来越大才能使越流补给量一直增大，而现实中含水层的越流补给情况其实非常复杂，但由于深层地下水的过度开采造成的对浅层地下水的袭夺确实客观存在，其对整个含水层调蓄能力的影响以及伴生的咸水下移等一系列负面影响必须在未来地下水资源的管理中被考虑。

弱透水层和承压含水层的压密释水量随地下水位的下降越来越小，这对于地下水资源管理非常重要，在计算未来地下水可开采量时，需要考虑这种变化的影响。更为重要的是，由于弱透水层的压缩变形非弹性不可逆，所以其储水系数的减小为永久性减小，未来即使对地下水进行回灌和补给，使地下水位上升到初始水平，但弱透水层的压缩变形却得不到回弹，其给水能力并不能恢复。所以一旦地面沉降发生，即使对地下水进行补给，含水层系统储存的水资源量并不能完全恢复，这对地下水库容也是一种永久性的损失，含水层系统所能储蓄的水资源量变少，对于地区的用水安全和地下水资源的保障能力影响巨大。

沧州深层地下水位累计下降约100 m，根据本模型计算出累计压缩变形量约为4 m，这与实际观测到的2.6 m[16]有一定差距。分析原因主要有以下几方面：(1)因为本文仅做了一个简单的一维沉降模型，而实际的沉降区情况特别复杂，现实中监测到的沉降量为开采量作用到一个区域上的结果，本文是开采量作用到一个点造成的形变量，必然比实际测量数据大很多；(2)很多计算所需的参数并无准确的观测值，所以大多根据经验确定，会造成很大的误差；(3)本次研究设定的条件较为理想化，不考虑补给、沉降时间滞后性等因素，也是造成误差的来源之一。本文的目的并不是将模拟的结果与实际情况做比较，而是通过数学手段对沉降发生过程中的水文地质参数的变化规律做计算。虽然仅基于沧州地区的具体情况进行了模拟研究，但通过模拟计算揭示的地面沉降造成的承压含水层给水能力降低应该是地面沉降严重地区普遍存在的现象，这些地区虽然地面沉降造成的地质灾害问题已经引起高度重视，但伴随其发生的水文地质参数变化造成的承压含水层系统给水能力的永久性降低也迫切地需要研究者和管理者重视。

4 结论

(1)地面沉降过程中，随土层孔隙度减小，弱透水层储水系数和渗透系数均减小，以沧州地区为例，弱透水层储水系数最高可减小77%、渗透系数最高可减小50%。

(2)弱透水层的非弹性压缩释水是深层地下水开采量的主要构成来源，其对开采量的贡献量随弱透水层压缩而逐渐减小，模拟发现，承压水位累计下降达200 m时，非弹性释水量可减小77%；承压含水层的弹性释水量占开采量比例较小，呈逐渐减小的趋势，本次模拟中最多也可减小77%；潜水层越流补给量随潜水位与承压水位差增大而增大，最高可增大5倍。整个含水层系统深层地下水的给水能力逐渐减小。

(3)弱透水层的压缩变形不可恢复，造成的含水层系统地下水储量是永久性损失，即使通过深层地下水的禁采或人工补给等手段，承压含水层水位能够有一定水平恢复，但损失的弱透水层水储量并不能得到补充，储水系数减小也降低了含水层系统的调节能力。地面沉降导致的整个深层含水层系统的给水能力降低，对地下水安全和供水保障能力有很大的影响，假设发生极端的自然灾害，沉降区的深层地下水能否还能应对应急用水需求值得深入研究。

(4)实际的含水层系统非常复杂，本文为深层地下水开采导致的地面沉降造成的水文地质参数和含水层系统给水能力的变化规律提供了一种新的研究手段和思路。本次研究对沉降的时间滞后性、区域差异性等并未考虑，今后将从这几方面继续进行深入研究。