运城东郭井水位异常与孔隙度变化的分析

姚林鹏，刘国俊，李 民，李晓锐，王洪峰

(1.山西省地震局夏县中心地震台,山西 夏县 044400；2.太原大陆裂谷动力学国家野外科学观测研究站,山西 太原 030025)

0 引言

由于地下水存在普遍性、流动性和难压缩性，当它处于一个封闭的承压系统时，井水位的动态变化能够客观、灵敏地反映地壳中应力应变状态的改变。前人已有的研究认为，承压井水位的动态能够反映含水层孔隙率的变化[1]，而研究含水层孔隙率的变化能够有助于分析井水位出现异常变化的成因[2]。

杨柳等，武晓军等利用华北地区及黑龙江地区承压井水位资料进行含水层参数的计算，来反映该地区应力应变状态[3-4]；方慧娜利用井水位气压效应反演含水层参数，对汶川地震前后四川、安徽地区的南溪、邛崃和巢湖井的井水位动态进行了分析，认为地震导致了含水层参数的变化或者对井孔周围地质构造产生了影响[5]；丁风和等研究了苏皖地区5口井含水层介质在孔隙度和渗透系数的变化情况，结合井-含水层系统，在受力与层储水水均衡两种不同状态下的微动态变化模型，分析了各井水位出现同步上升变化现象在水动力学上的成因机制[6]。

该文针对运城东郭井2013-2017年出现了多次较大幅度的水位异常，拟通过地下流体动力学和岩石力学等相关理论，对东郭井数字化水位、配套气压等资料进行处理，求解井水位动态变化中的潮汐因子与气压系数，进一步求解含水层孔隙度。在讨论含水层在受应力变化和水均衡，两种不同状态下的微动态变化后，梳理东郭井-含水层系统中孔隙度随时间的变化情况，与井水位的异常变化进行分阶段对比，明确2013-2017年东郭井几次大的水位异常变化，所反映的地下水动力或构造活动的变化信息，为晋陕豫交界区的区域震情研判和跟踪提供一定的科学依据。

1 观测井基础信息与资料选取及处理

1.1 观测井基础信息

东郭井是山西省地震局“十五”期间，自筹资金建成的一口地震地下流体观测井，位于汾渭断陷带的运城盆地，海拔高程400 m，构造上位于中条山断裂，呈弧形向东南的凸出点(见第34页图1)。井深614.71 m，观测含水层深度485 m，为第三系上新统松散岩，地下水类型为孔隙承压水，成井时水位埋深21.18 m，水面温度10 ℃。

1.2 资料选取及处理

选取2013-2017年，东郭井“十五”数字化水位变化数据，和同测点配备的气压变化整点值数据，进行如下整理。将水位和气压数据进行逐月检查、整理，缺数时，结合多次样条插值与一般多项式分段拟合结果进行替代，干扰及校测等产生的错误数据加以改正，保证数据的真实可靠；利用回归法获得仅余气压项和固体潮项的水位整点值资料。

2 东郭井水位变化特征

2013-2017年东郭井水位变化较大，变幅在40.9～53.1 m之间，共出现3次大的异常变化。分别为2014年9月突升异常、2015年9月突降异常、2017年2月突升异常。根据水位的变化动态，将5年间的变化分为4个阶段：第1阶段(2013年1月-2014年9月)，水位变化平稳，表现为“夏低冬高”的年变形态；第2阶段(2014年10月-2015年9月)，水位持续性下降，且下降速率较为平稳；第3阶段(2015年9月-2017年2月)，水位持续性上升，波动较大；第4阶段(2017年2月-12月)，水位持续性下降，下降速率较为平稳(见图2、表1)。

图1 山西中条山地区水文地质图Fig.1 Hydrogeological map of Zhongtiaoshan area, Shanxi Province

图2 东郭井水位年动态变化图Fig.2 Annual dynamic change of water level in Dongguo well

3 含水层参数的动态变化特征

3.1 气压系数和潮汐因子的计算

表1 东郭井水位不同阶段的变化特征Table 1 Variation characteristics of water level in Dongguo well at different stages

考虑把水位变化理解为趋势变化成分、随机变化成分、潮汐影响成分和气压影响成分等几个分项。利用卷积回归法进行潮汐影响成分的去除，然后结合一般多项式分段拟合结果,去除水位动态中随机变化成分和趋势变化成分，余下的水位变化数据将只包含气压影响成分分项。首先,将趋势变化表现明显的气压数据进行一阶拟合去趋势;然后,采用一元回归方法,利用仅剩余气压影响分项的水位数据和校正后的气压数据,求解气压影响系数Bp[7]。

计算潮汐影响因子的过程是先利用卷积回归法，进行井水位变化中气压影响成分的去除，再结合一般多项式分段拟合结果，去除水位动态中的趋势变化成分和随机变化成分，余下的水位变化数据只剩潮汐响应分项，使用Venedikov调和分析软件计算各波群的潮汐影响因子Bg[8]。

3.2 孔隙度计算

在假设含水层不排水的状态下，水位气压系数Bp和潮汐因子Bg可分别表示为[9-10]：

(1)

(2)

由式(1)、(2)可得：

(3)

式(1)、(2)中：α为含水层固体骨架的体积压缩系数；β为水的体积压缩系数，取值为4.69×10-10(m2/N)；n为含水层的孔隙度；ρ为水的密度；g为重力加速度；ρ和g相乘为水的重度，其值为0.098 hPa/mm。由式(1)、(2)得出的公式(3),消去了固体骨架的体积压缩系数α，在已知气压系数Bp与潮汐因子Bg，只有一个未知数n的情况下，通过将各月气压系数与潮汐因子进行代入，即可求得含水层的孔隙度n。

3.3 孔隙度变化特征分析

第35页图3的2013年-2017年东郭井孔隙度月值结果表明，第35页孔隙度出现3次大的变化，分别为2014年10月“突升”变化、2015年9月-12月“突升-突降”变化、2016年11月-2017年2月“突升-突降”变化。3次变化将整体的孔隙度数据分成4个阶段：第1阶段(2013年1月-2014年9月)，孔隙度在15%左右波动变化，从2014年4月开始缓慢下降；第2阶段(2014年10月-2015年8月)，孔隙度在25%波动，从2015年4月开始缓慢下降；第3阶段(2016年1月-10月)，孔隙度在15%左右平稳波动；第4阶段(2017年2月-12月)，孔隙度在20%左右波动变化。

图3 观测含水层孔隙度变化图Fig.3 Porosity change of observed aquifer

4 含水层孔隙度动态变化特征

在观测含水层为承压水条件下，井水位测区含水层因受力状态改变，而引起地下水位出现动态变化的过程中，当含水层受压变形时，孔隙度减小，孔隙压力增强，井水位呈现上升形态；反之，当含水层受张力发生变形，孔隙度增加，孔隙压力减弱，井水位呈现下降形态(见表2)。受力状态的改变和水均衡过程，均可以造成水位发生上升或下降形态的变化，要区分两种成因的不同，关键是了解在水位动态发生变化的过程中，孔隙度的变化情况。

表2 含水层受力状态变化引起的地下水动态变化关系表Table 2 Dynamic change of groundwater caused by change of stress state of aquifer

含水层介质由水均衡过程改变,而引起地下水位动态变化的过程中，当含水层补给量大于泄流量，补给增强时，孔隙度增大，孔隙压力增强，井水位呈上升变化形态；反之，补给量小于泄流量，水位呈下降变化形态(见表3)。

表3 含水层水均衡变化引起的地下水动态变化关系表Table 3 Dynamic change of groundwater caused by water balance change of aquifer

东郭井水位4个阶段的变化与孔隙度的变化具有准同步性，从各阶段的整体特征来看，第2阶段相较于第1阶段水位升高，孔隙度同步增大，反映出水位的上升变化与补给量多于泄流量、补给过程增强有关；第3阶段相较于第2阶段水位降低，孔隙度整体同步较小，反映出水位的低位变化与泄流量多于补给量、补给减弱有关；第4阶段相较于第3阶段水位升高，孔隙度同步增大，反映出水位的上升变化与补给量多于泄流量、补给增强有关(见图4)。

东郭井水位4个阶段之间的异常大幅变化，从孔隙度变化的角度来看，水位上升与下降和孔隙度增大与减小同步，与水均衡过程的孔隙度变化一致，与应力改变过的孔隙度变化不一致，反映出水位变化与区域含水层受力状态的改变相关性较弱，主要与含水层介质因水均衡状态过程发生变化有关。

图4 水位变化与含水层孔隙度变化对比图Fig.4 Comparison of changes in water level and porosity of aquifer

5 结论与讨论

针对东郭井水位2013年-2017年发生几次大的异常变化，通过Venedikov调和分析及气压影响系数、含水层参数(孔隙度)等的计算和分析，得出以下结论：

(1) 东郭井水位与孔隙度均反映出较明显的4个阶段性变化，各阶段的时间跨度及变化发生点具有准同步性。

(2) 引起东郭井水位3次异常变化的原因主要是水均衡的问题，即不同时间段内补给量与泄流量的大小问题，与区域含水层受力状态的改变相关性较弱。

(3) 通过计算含水层孔隙度的变化情况，能较好地应用于水位异常变化的成因分析。

该研究主要对东郭井水位2013-2017年出现的4个阶段变化进行分析，通过孔隙度整体的变化来讨论水位异常变化的因素，但在每次异常变化前后，孔隙度都会出现短期的大幅波动，反映出水均衡过程导致水位异常变化的复杂性，还需要进行更加深入的分析与讨论。