无极冀20井地震前后水流运动方向研究

董 博，纪春玲，章 阳，周安聘

（河北省地震局石家庄中心台，石家庄，050021）

0 引言

对承压含水层井孔而言，井水位对体应变潮汐的响应是井水位观测中普遍存在的一种现象，井水位潮汐中隐含有丰富的井-含水层系统信息[1]。深入分析井水位潮汐对推算含水层的各种参数，衡量井水位对应力应变响应量、了解含水层动态、解释井水位微动态等都具有重要的科学意义和现实必要性。

Bredehoeft以体应变固体潮理论为基础，结合水文地质学，从理论上分析了井-含水层系统对体应变潮汐的响应，并提出一种利用井水位潮汐信息来确定含水层的孔隙度、骨架体变模量和储水系数的方法[2]。之后，不同研究者对利用潮汐信息确定含水层参数的方法进行了补充[3-5]。其中，Hsieh等推导了井水位对含水层孔隙压力响应的振幅比和相位差数学表达式，并给出了利用井水位潮汐相位和振幅比来确定井水位渗透系数和储水系数的方法[4]；Elkhoury等利用Hsieh等提供的方法分析了加利福利亚地区两口井水位观测资料的潮汐振幅比和相位差，认为大震发生后的潮汐振幅比和相位变化是由于地震波增加含水层渗透系数所致[6]；Bower等运用孔弹性应力平衡理论推导出了含水层孔压与引潮位关系式，分析了井水位受裂隙产状及径向排水的变化特征[7]；刘春平等在此基础上，进一步对受垂向排水的井水位特征进行了研究，总结出井水位变化主要受裂隙产状、排水影响，并对变化特征进行归纳[8]。

此外，井水位潮汐是评价水位观测质量好坏的一个重要参数，对井水位的潮汐因子、振幅和相位等潮汐参数的研究可能获得震前异常信息[9-11]。前人对井水位潮汐因子、振幅和相位等潮汐参数动态变化的可能原因、影响因素讨论较少，且对地震前后井水流运动方向的研究略显不足。

本文以固体潮理论为研究基础，在确定了潮汐为该井水位变化主要影响因素的基础上，利用Baytap-G提供的潮汐分析方法[12],计算了无极冀20井水位2009—2019年M2和O1波的潮汐因子、相位及振幅值，深入分析了几次远场大震前后该井潮汐因子、振幅以及相位等潮汐参数的动态变化特征，并试探性地为潮汐参数的变化找到理论依据支持，对地震波能改变含水层储水系数、渗透系数等水文地质参数并使孔隙压力得到重新分布进行了验证，最后得出地震前后水流运动方向，为解释远距离地震触发和水文地质响应现象提供了实例。

1 观测井孔基本情况

无极冀20井观测站地处冀中拗陷，无极低凸起高点，衡水断裂北侧（图1）。该井位于河北省无极县境内，地理坐标为（38.2°N，114.9°E）。无极井于1989年开始观测，孔口标高44.5 m，井深2 984 m，属华北深井网。地下水类型为岩溶裂隙承压水，含水层为震旦系、寒武系灰岩、白云岩，观测层上覆泥岩、砂岩，下伏蚀变的白云岩，可能与含水段发生水力联系。井-含水层水流运动示意图见图2。无极井1986年开始进行静水位模拟观测，仪器为SW-40，2001年开始进行数字化观测，现使用的数字化观测设备为ZKGD3000NL型地下水数据监测系统。记录数据质量良好、潮汐清晰。

图1 井孔地质构造图

图2 井-含水层水流运动模型

2 井-含水层水流径向、垂向运动成因分析

河北局所管辖流体井口数量多，一般情况下研究井-含水层系统时通常会忽略含水层与隔水层之间的水流交换，但实际上几乎所有井所揭露的含水层都会与其上面覆盖的隔水层以及下面接触的隔水层都会发生微弱的水流交换。为更深入研究井-含水层水流运动情况，充分考虑含水层与隔水层之间的水流交换，建立井-含水层水流径向、垂向运动模型是十分必要的。

据现有研究，潮汐井水位振幅和相位受到裂隙产状和排水两类因素的影响。裂隙产状影响下，潮汐M2和O1波的水位-引潮高振幅比比值差异不大，但位相差符号相反[7]。排水影响分为径向和垂向排水两种情况。在多孔介质含水层中，不排水条件下含水层潮汐孔压-引潮高振幅比为常数E。在饱水多孔介质中，井水位与含水层孔压相等，水位-引潮高振幅比M=E（对含水层）或M′=E′（对隔水层）。

3 资料预处理

收集无极冀20井2009—2019年水位记录整点值，剔除因模拟水位换纸引起的错误数据，并采用克里金（kriging）插值[13]对数据进行插值补全。因井水位数据所包含的信息复杂，受环境、场地、认知因素影响较大，为了确定潮汐对井水位的影响程度，首先对井水位整点值数据（2015年8月）进行频谱分析（FFT变换），得到影响井水位的主要频率段（图3～4）。图中显示频率为0～0.5 h-1的信息振幅达0.6 kpa，经分析无极冀20井影响水位观测的因素，认为可能与井的固有频率及井口周边环境有关 ，有待进一步核实。

图3 2015年8月井水位整点值曲线

图4 井水位频谱分析

由上图分析可知，观测数据中明显含有日潮、半日潮、三分之一潮信息，且半日潮振幅最大、日潮次之、三分之一潮的振幅最小，这与固体潮理论是一致的。无极冀20井同震响应良好，许多学者已经做过研究[14-15]，不再叙述。

4 井水位潮汐特征分析

基于井水位潮汐变化，通过baytap-G程序、取计算窗长720 h、滑动步长168 h[1]计算无极冀20井2009—2019年水位整点值数据，得到M2和O1波的潮汐因子（各个谐波的观测振幅与理论振幅之比）、相位和振幅及其误差。由计算结果可知：M2波计算误差最小，其振幅误差主要分布在0～0.08 GPa之间，相位误差在0°～2.1°之间；O1波振幅误差主要分布在0～0.24 GPa，相位误差在1°～16°之间。

绘制M2波变化特征（图5）、M2波和O1波相位变化曲线（图6）。由图5可以看出，2009—2019年M2波潮汐因子、相位的动态变化可能与这期间中强地震的发生有关（地震已用箭头标出）。地震基本情况见表1，表中能量密度是根据前人推导的地震波能量密度公式计算得出[16-17]。

图6 M2波、O1波相位动态变化图

5 井-含水层水流方向改变分析与检验

5.1 井-含水层水流方向改变分析

实际情况下，井-含水层系统会受到排水情况以及裂隙产状共同影响，潮汐参数表现形式复杂[7-8]（表2）。

由图5结合表1～2可得出，2009—2019年无极冀20井的M2波潮汐因子、相位动态变化分为4个阶段。

第一阶段为2009—2011年。初潮汐因子、相位变化趋势不稳定，波动较大。由表2可以看出2009年9月30日苏门答腊7.7级地震、2010年4月7日北苏门答腊8.0级地震前后，M2波、O1波相位均为正值，M2波潮汐因子稍大于O1波潮汐因子，E'/E＜1，说明隔水层和含水层之间发生水流交换。由表1可知，这种情况井水位为垂向排水，水流运动方向为垂向。

表1 几次中强地震基本情况

表2 井水位变化特征成因分析

由图5～6可以看出，2009—2010年M2波、O1波相位值多次出现正负号的转变，且数值变化幅度较大。根据刘春平等[8]研究结果，这种情况的发生可能跟隔水层裂隙发育产状有关。2009年苏门答腊南部7.7级地震、2010年北苏门答腊8.0级地震，这两次地震M2波潮汐因子、相位均出现变化，且地震能量密度最大为2.38×10-4J/m3，与Wang[16]和Manga[17]给出的引起潮汐因子和相位发生变化的能量密度下限为10-3J/m3不符，此现象可能与裂隙发育有关。综合得出，此时间段内水流以垂向运行为主，受裂隙影响。

图5 M2波潮汐因子、相位、振幅随时间变化曲线

第二阶段为2011—2014年。2011年3月11日日本本州9.0级地震发生后M2波、O1波相位由正值变为负值，发生正负号的改变，且两波数值差异不大，说明井水位主要受径向排水。2011—2014年O1波潮汐因子均稍小于M2波潮汐因子，E'/E＜1说明该井水位变化主要受径向排水为主、垂向排水为辅，水流主要运动方向为径向，也说明震后含水层渗透系数得到改变并导致孔隙压力重新分布。

为了更进一步验证水流运动方向的改变，计算了3月11日日本本州9.0级地震的地震波能量密度为4.21×10-2J/m3。根据前人给出的研究结果：能引起含水层储水系数、渗透系数等水文地质参数改变的地震波能量密度最小数量级为10-3J/m3，这说明本次日本本州地震足以引起井-含水层水文地质参数的改变，从而影响水流的运动方向。

第三阶段2015年尼泊尔地震后至2018年下半年，2015年4月25日尼泊尔地震的发生使得M2波、O1波相位值均由负变为正，说明无极冀20井含水层与其相邻的隔水层之间发生了水流的交换，含水层在压应力的作用下向其上覆或者下伏隔水层排水。在E'/E＜1条件下，M2波潮汐因子均稍大于O1波潮汐因子，理论上排除了径向排水的可能，进一步证实了该井主要受垂向排水影响，水流运动方向为垂向。此现象维持至2018年8月19日斐济地震发生后，说明尼泊尔地震的发生使得无极冀20井由径向排水变为垂向排水。同样，计算此次地震波能量密度为1.399×10-3J/m3，符合前人研究结果。

由图5可知，尼泊尔地震发生后井水位M2波相位、振幅均变大，而O1波振幅变化不明显，表明该井井水位潮汐响应满足过渡区的特点。2016年所罗门群岛地震前后M2波和O1波振幅和相位没有明显变化，且与尼泊尔地震后基本一致。根据Doan等[18]建立的不同频率波振幅和相位值与渗透系数的模型，可知尼泊尔地震前无极冀20井满足排水条件（径向），震后该井满足不排水条件（径向），再次论证了尼泊尔地震后井水位受垂向排水影响，水流运动方向为垂向。

第四阶段，由2018年下半年开始，M2波、O1波相位调整至日本本州地震后的负值水平，潮汐因子、振幅变化也趋于平稳。斐济地震能量密度为2.229×10-5J/m3，能量水平不足以改变M2波的相位。通过计算2018年下半年全球MS≥5.0地震的能量密度均不足以引起M2波相位改变，因而初步推测，这可能是震后区域静态应力的调整过程。调整后两波相位均为负值，且差异不大，O1波潮汐因子稍小于M2波，说明应力场调整后的井水位主要受径向排水影响，E'/E＜1条件下垂向排水为辅，水流主要运动方向为径向。此现象一直持续到2019年末，根据M2波、O1波动态变化，可推断出这种情况将会维持一段时间，直至有新的大能量地震的发生。

绘制2009—2019年无极冀20井M2波和O1波的潮汐因子动态变化曲线（图7）。根据刘春平等[8]给出的理论，结合图6、表2可看出：2012—2014年M2波潮汐因子基本大于O1波，说明无极冀20井在该时间段属于表1中所列第1、2种情况；M2波相位基本小于O1波，因此该时间段内该井以径向排水为主，E'/E＜1条件下垂向排水为辅，水流运动方向主要为径向；2016年反之，为垂向排水，水流运动方向为垂向。与上述分析结果一致。

图7 M2波、O1波潮汐因子动态变化图

5.2 水流运动方向检验

由Hsieh等[4]建立的仅径向排水情况下井-含水层响应模型可知，M2波振幅系数与相位呈同向变化，而Doan等[18]得出仅垂向排水情况下井-含水层模型得到，M2波振幅系数与相位呈反向变化据此对上述分析结果进行检验。

绘制无极冀20井2012—2014年、2016年M2波潮汐因子和相位动态曲线图（为避免区域静态应力场调整干扰，本文只绘制了2016年潮汐参数曲线）（图8～9）。由图可知无极冀20井2012—2014年M2波潮汐因子和相位呈同向变化，主要为径向排水，水流运动方向主要为径向；2016年M2波潮汐因子和相位呈同反变化，为垂向排水，水流运行方向为垂向。进一步验证了上述结论。

图8 2012—2014年M2波潮汐因子和相位动态曲线

图9 2016年M2波潮汐因子和相位动态曲线

6 结语

本文通过对无极冀20井水位资料的潮汐参数分析，得到以下结论和初步认识。

1）通过对无极冀20井水位数据FFT变换可得出，该井水位变化主要受潮汐影响，观测质量较高。

2）无极冀20井潮汐因子和相位的动态变化可能与远距离强震有关。通过对无极冀20井2009—2019年水位整点值潮汐参数分析认为：①2009—2010年M2波和O1波相位值多次出现正负号的转变，且两波数值大小不定，表明此时间段内水流以垂向运行为主，受裂隙影响。②日本本州地震能量密度为4.21×10-2J/m3，震后水流运动方向由垂向变为径向，说明震后含水层渗透系数得到改变并导致孔隙压力重新分布；③尼泊尔地震能量密度为1.399×10-3J/m3，震后水流运动方向由径向变为垂向，含水层和隔水层之间发生水流交换，在压应力作用下含水层向隔水层排水；④2018年下半年井水位主要受径向排水影响，E'/E＜1条件下垂向排水为辅，水流主要运动方向为径向。

3）根据前人不同的研究方法、模型对所得结论进行检验，结果符合预期。

4）2009—2010年M2波和O1波相位值多次出现正负号的转变，根据前人研究结果，这可能与水流运动受裂隙影响有关，符合无极冀20井的地质构造条件。