基于数值模拟与含水层垂向应力反演的静乐井水位异常分析1

张淑亮 李宏伟 吕 芳 陈 慧

（山西省地震局，太原 030021）

引言

静乐井水位是山西乃至华北地区映震能力较强的地下流体观测井之一，井水位异常能客观反映区域应力场的变化，全球7 级以上地震均能记录到水震波。1983 年始测以来出现的几次高值异常与山西北部至晋冀蒙交界地区强震活动时间具有准同步性，在1989 年大同-阳高Ms6.1 级地震、1991 年忻州Ms5.1 级地震、1991 年大同-阳高Ms5.8 级地震、1998 年张北Ms6.2 级地震、1999 年大同-阳高Ms5.6 级等中强地震发生前有较好的异常显示。因此将静乐井水位高值异常作为晋冀蒙交界地区地震趋势预测的重要判据（张淑亮等，1997；车用太等，1999；王爱英等，1998）。 2017 年8 月静乐井水位又出现大于多年平均值的高值异常，到目前为止异常仍在持续。按照以往震例时间指示意义，目前已接近发震时间，未来晋冀蒙交界区震情形势如何发展？目前的异常是否真实反映了构造活动加剧所引起的区域应力场和深部地球物理场的改变？是否与中强震孕育过程有关？静乐井水位目前的异常特征不足以解析异常区的应力状态，从而也直接影响对晋冀蒙交界地区地震形势的认识和未来地震趋势的判定。

多年的观测表明，影响静乐井水位动态变化的主要干扰因素为短时强降雨与河流水体的荷载作用。无论是短时间的强降雨还是洪水期河流的荷载作用，均可在短时间内造成井水位的快速变化（张淑亮等，2005；车用太等，2004；鱼金子等，1994；张昭栋等，1989；张昭栋等，1990）。目前对这2 种干扰因素的识别与排除基本上以定性分析为主，定量分析方法较少。特别是对于同时受2 种因素影响的静乐井水位而言，在异常分析与判定工作中，仅采用与同期降雨量对比的方法分析井水位变化与降雨量的关系，而河流水位变化对井水位的影响分析工作至今尚未开展。为此，本文利用数值模拟与含水层垂向应力反演的方法，对降雨与河流荷载作用下的静乐井区应力进行定量计算，并从区域应力场的角度分析2 种方法计算结果的差异，以期为地震趋势判定提供较可靠的预测依据。

1 井区概况

静乐流体观测井位于吕梁山断块隆起区的静乐向斜，位于东碾河断裂带上，地处东碾河南岸，与东碾河相距约130m，其北部和南部均为山区（见图1）。静乐井深362.92m，井区为河谷地区，主要发育有下古生界灰岩、白云岩喀斯特含水层及第四系冲洪积砂砾石含水层，2 个含水层间有较强的水力联系。观测层岩性时代为中奥陶系灰岩裂隙溶洞水，其顶板埋深为46.3m，单位涌水量84.5L/（s·m），渗透系数为2.41—11.5m/d。观测层具有厚度大、透水性强的特点。含水层岩石溶洞发育，静乐井在230—290m 处通过破碎带，地下水由北西经井区向东南方向排泄，井区处于强径流区，地下水补给主要来自北、南山区基岩裸露区的大气降水补给和东碾河河水沿断裂破碎带入渗补给。

图1 井区地质构造 Fig. 1 Well area tectonic map

2 强降雨与河流荷载作用下静乐井水位动态特征

多年的观测表明，大于15mm 的集中降雨将引起静乐井水位的快速上升（张淑亮等，2005），可能由于降雨在地表形成荷载，使得含水层受力状态发生变化，引起水井水位上升。东碾河水位变化是影响静乐井水位变化的又一重要因素。当大面积洪水作用通过井区时，井水位迅速大幅上升，几乎没有时间滞后，洪水期水位上升，洪水过后水位下降，水位升降幅度取决于洪水深度。

图2 静乐井水位与降雨量和东碾河流量的对比 Fig. 2 Comparison of Jingle well level with rainfall and Dongnian river discharge

图2 所示为静乐井水位与井区降雨量和东碾河流量的对比。由图2 可知，2017 年8 月20 日左右单日降雨量为82.8mm 时，东碾河流量由原来的0.777m3/s 增至9.8m3/s，增加了9.023m3/s； 静乐井水位由9.074m 升至7.747m，水位上升幅度达1.327m。由此可知，降雨与河流荷载作用是影响静乐井水位动态变化的主要因素，但降雨结束后井水位并未迅速回落至上升前的水平。张淑亮等（1997）曾对静乐井水位高值变化与降雨的关系进行了较系统的研究，认为降雨量的增大虽可在短时间内造成井水位升降变化，但不能改变水位的年变动态，只有在区域应力场增强的背景条件下，降雨引起的附加力源才可能造成井水位年动态的异常变化，降雨荷载是区域应力场变化的诱发因素，也可能是在2017 年雨季结束后井水位仍处于高值的主要原因。

3 降雨与河流荷载作用下井区应力数值模拟分析

利用井区周围地壳结构、不同地层介质弹性参数建立三维地质模型，预测降雨与河流荷载作用下的井区应力。由于本文基于数值模拟与含水层垂向应力反演的方法进行定量计算，因此，依据静乐井井区地层岩性特征和井-含水层特征参数选取计算参数。因静乐井位于河谷地区，故计算时以变质岩系中大理岩力学性质指标经验数据为依据。

3.1 三维地质模型建立

根据静乐井区地层结构、河流位置及研究目标，设定有限元尺寸为1km×1km×360m，河流宽76m。采用SOLID185 单元，模型采用由面拉伸成体的建模方式，取井孔周围围岩弹性模量为5×1010MPa，泊松比为0.25，密度为2600kg/m3。

网格划分的目的是对模型实现离散化，并用适当数量的网格单元得到最精确的解。本文选用Workblench 中的自动网格划分工具，以满足模型要求。网格疏密程度直接影响计算结果的精度，但网格加密将增加CPU 计算时间，且需更大的存储空间。为保证模型网格密度和计算精度，本文选择的网格尺寸为1.5m，同时适当增加井区附近区域的网格密度。其他参数均采用默认值，可满足模型的需求。

载荷和约束是ANSYS 软件求解计算的边界条件，由所选单元的自由度形势定义。由于本文主要分析河流与降雨荷载的影响，因此选择垂向应力加载。为防止模型移动，模型底端与侧面采用固定约束，限制发生移动和转动。

考虑数值模拟结果的多样性及模型造成的误差，建模时选择不同的边界条件，并反复测试计算，以可反映降雨与最大河流荷载作用对井区的最大影响范围作为模型的边界条件。由于降雨与河流荷载的作用主要表现为垂向应力，因此模型初始地应力只考虑重力。因本文通过数值模拟与反演的方法分析静乐井水位在降雨与河流荷载作用下的垂向应力在量级上是否具有相似性，进而确定影响水位变化的主要因素，对计算精度的要求相对较低，故选取物理力学参数时以研究区内地层岩性的平均值作为模型计算参数。

3.2 结果分析

由于2017 年是近年来静乐井水位年变幅最大的年份，也是降雨与河流荷载影响最大的年份，故本文以2017 年8 月洪水期东碾河水位、静乐井区降雨量和静乐井水位为研究对象，计算降雨和河流荷载作用下的静乐井井区的垂向位移和垂向应力。

（1）单降雨荷载的影响

以2017 年8 月降雨量作为加载体，将降雨量等效为荷载，计算得到静乐井垂向位移和垂向应力如图3 所示。由图3 可知，降雨荷载对静乐井产生的垂向位移约为0.6mm，垂向应力约为12hPa。

（2）单河流荷载的影响

以洪水期东碾河水位最大变化量作为加载体，并将河水水位等效为荷载计算得到静乐井垂向位移和垂向应力如图3 所示。由图3 可知，洪水期河流荷载对静乐井产生的垂向位移约为0.3mm，垂向应力约为15hPa。

图3 单降雨荷载作用下静乐井垂向位移和垂向应力 Fig.3 Vertical displacement and vertical stress of Jingle well under single rainfall load

图4 单河流荷载作用下静乐井垂向位移和垂向应力 Fig.4 Vertical displacement and vertical stress of Jingle well under single river load

（3）降雨与河流荷载的综合影响

由于静乐井水位变化是在因降雨与河流荷载的共同作用，单河流和单降雨的影响仅能反映井水位的部分变化，因此，需考虑二者的共同作用。以2017 年8 月东碾河水位最大变化量和降雨量作为共同加载体，并将二者等效为荷载，计算得到静乐井垂向位移和垂向应力如图5 所示。由图5 可知，降雨和河流荷载共同作用对静乐井产生的垂向位移约为0.7mm，垂向应力约为30hPa，与单河流和单降雨影响的二者之和基本一致，表明降雨与河流荷载作用是影响静乐井水位变化的主要因素。

图5 河流与降雨荷载共同作用下静乐井垂向位移和垂向应力 Fig.5 Vertical displacement and vertical stress of Jingle well under combined action of river and rainfall loads

4 静乐井-含水层垂向应力反演

数值模拟结果表明，降雨和河流荷载是影响静乐井水位变化的主要因素。为进一步验证数值模拟结果能否真实反映井区实际应力变化，还需将其与实际水位观测值所反映的井区应力进行对比。在由井水位变化反演井区应力场变化方面，目前取得了一些进展，如黄辅琼等（2004）利用华北地区40 多个深井水位动态变化资料，探讨研究了华北地区现今构造应力场状态；孙小龙等（2011）运用小波分析法提取出能反映水位多年动态变化的趋势信息，反演出华北地区多年构造应力场的变化特征；丁风和等（2015）利用部分含水层介质参数、井水位变化量与含水层垂向应力变化量的关系式，定量分析了张渤带地区构造应力场的动态变化过程。

本文根据静乐井-含水层介质参数及井水位变化量等数据，反演井水位变化引起的含水层垂向应力。应力反演所用数据时段为2016—2017 年（根据东碾河水位与流量始测时间确定）。反演前首先对观测数据进行预处理（丁风和等，2015），去除长周期干扰因素、降雨、气压等年际变化对水位的影响；然后采用去趋势法和傅里叶周期法，剔除地下水开采和降水补给等长周期干扰对水位的影响。

4.1 反演所需各类参数计算

（1）气压系数与潮汐效率计算

利用卡尔曼滤波方法消除地球固体潮对所选时段水位、气压整点值数据的影响。利用滤波后的日均值数据，采用高阶差分和一元线性回归等方法获得静乐井气压系数。利用维尼迪科夫潮汐调和分析方法获得静乐井水位M2波潮汐效率。气压系数和潮汐效率计算结果如表1 所示。

表1 静乐井-含水层多种参数计算结果 Table 1 Calculation results of various parameters of Jingle well-aquifer

（2）静乐井-含水层介质参数计算

不排水状态下，井水位气压系数和潮汐效率可分别表示为（Bredehoeft，1967；李春洪等，1990；张昭栋等，1993）：

式中BP为井水位气压系数；Bg为M2波潮汐效率；n为含水层孔隙度；α为含水层固体骨架体积压缩系数；β为含水层内水的体积压缩系数；ρ为水的密度；g为重力加速度（ρg=0.098hPa/mm）。因此，n、β可依据式（3）与表1 中的潮汐效率和气压系数得到。根据式（1）或式（2）可求α。

4.2 静乐井-含水层垂向应力反演

利用表1 数据，根据水平层状含水层模式下，含水层介质参数（孔隙度、水和固体骨架体积压缩系数等）、井水位变化量与含水层垂向应力变化量的关系式（张昭栋等，1987）等来计算含水层垂向应力，含水层垂向应力计算如下：

式中 Δσz为含水层垂向应力变化量；E为含水层固体骨架的杨氏模量（E= 1/α）；△H为剔 除地下水开采、降雨和气压影响后的含水层应力变化引起的压力水头变化量，即井水位变化 量。当井-含水层系统所受应力增强，即 Δσz>0 时，井水位上升，水位埋深H变小，其变化量ΔH<0；当井-含水层系统所受应力减弱，即 Δσz<0 时，井水位下降，水位埋深H变 大，其变化量ΔH＞0。利用式（4）与表1 中有关参数计算静乐井2016—2017 年含水层垂向应力，结果如表1、图6 所示。可知，静乐井含水层垂向应力高值时段与降雨量高值时段、东碾河流量高值时段较吻合，再次表明降雨与河流荷载作用是造成含水层垂向应力高值变化的主要因素，井水位的大幅突升变化与垂向应力增大密切相关。

5 静乐井水位异常分析

5.1 模型预测值与井水位反演值对比分析

由数值模拟结果可知，无论单河流荷载还是单降雨荷载效应均对洪水期静乐井水位造成了一定影响，二者造成的垂向位移累计0.9mm，垂向应力累计27hPa。河流和降雨综合影响模型结果也得到类似结果（垂向位移累计0.7mm，垂向应力累计30hPa）。由含水层垂向应力反演结果可知，洪水期含水层的应力显著增强，由正常时段的4.2hPa 增至79hPa 左右。2 种方法计算结果均表明静乐井水位上升变化与河流荷载效应和降雨荷载效应密切相关。但2 种方法计算结果存在一定差异，降雨与河流荷载作用下井区垂向应力数值模拟结果小于含水层垂向应力反演值，相差39hPa。

图6 东碾河流量、静乐井-含水层垂向应力与降雨量的对比 Fig. 6 Flow rate of Dongnian river, vertical stress of Jingle well-aquifer and rainfall comparison chart

由于数值模拟计算得到的垂向应力是研究区自身重力、降雨与河流荷载作用下产生的垂向应力之和，而由井水位反演的垂向应力不仅包含上述因素产生的垂向应力，可能还含有一些因构造活动引起的垂向应力变化。事实上2 种方法的计算结果的确存在一定 差异。虽然2 种计算方法均受建模条件和参数选取的影响而存在计算误差，但根据以往的研究，误差值一般不会大于计算值，因此2 种计算结果的差异性是客观存在的。按照表1 中静乐井水位平均气压系数4.3521mm/hPa 可推算出2 种方法垂向应力计算结果差值引起的井水位变化幅度约为16.97cm，大于正常时段水位平均月变化幅度0.9cm 的水平。表明影响静乐井水位快速上升的因素除降雨与河流荷载作用外，可能还受其他因素的影响。

5.2 水位高值异常与构造活动的关系

为进一步探讨静乐井水位2017年8月以来的高值异常除受降雨与河流荷载影响外是否受构造活动的影响，本文从可反映区域应力状态变化的参数和测量值（缺震和b 值、地震矩释放、GPS 基线等）及同一构造区其他前兆测项变化特征与静乐井水位高值异常进行对比分析。

大量研究结果表明，b 值与其对应区域应力状态、地壳破裂强度有关（刘雁冰等，2017；朱艾斓等，2005；吴小平，1990）。1983 年以来山西北部至晋冀蒙交界区3 级以上地震缺震和b 值曲线有4 次异常过程（见图7），第1 次异常出现在1988 年，异常持续过程中在山西地震带出现一组中强地震活动，先后发生了大同-阳高6.1 级地震、侯马4.9 级和忻州5.1 级地震、大同-阳高5.8 级地震等，其中大同-阳高6.1 级地震被认为是华北地区在经历唐山大地震后进入新活跃幕的标志，是华北应力场趋于增强的结果；第2 次异常出现在1995 年初，异常持续过程中先后发生了1996 年5 月包头6.4 级地震、1998 年1 月张北6.2 级地震、1999年3 月张北5.6 级地震、1999 年11 月大同-阳高5.6 级地震等中强地震活动，既是1989 年开始的华北新一轮活跃幕的继续，又是华北应力场增强的结果；第3 次异常始于2007 年，这次异常过程中在山西北部至晋冀蒙交界区发生了几次4 级以上中等地震，部分学者认为，这组中等地震的发生是受汶川地震的影响，鄂尔多斯块体与华北平原块体相对挤压和扭错增强，导致山西带形变场与构造应力场由原来的构造拉张转为构造挤压，是汶川地震延迟触发的结果 （刘峡等，2013，朱艾斓等，2010；刘瑞春等，2014）；2017 年又出现同步异常变化，目前异常仍在发展，且缺震和b 值的几次异常与静乐井水位高值异常较为同步，因此，可认为静乐井水位高值异常反映的是华北区域应力场增强的过程，是山西地震带地震活动增强的重要判据（张淑亮等，1997）。2017 年静乐井水位高值异常也可能与区域应力场增强有一定关系。

图7 静乐井水位月均值与山西地震带3 级地震缺震和b 值的对比 Fig. 7 Comparison chart of monthly average water level of jingle well with magnitude 3 earthquake deficiency and b value in Shanxi seismic belt

由跨过静乐井的山西岢岚-山西太原GPS 基线时序曲线图8 可知，2016 年以来基线由原来的下降趋势（基线缩短）变为缓慢上升（基线拉长），应力状态在2016 年发生了变化，由静乐井所在区域地震矩加速释放时图9 可知，2017 年井区存在明显的加速释放低值区。GPS基线和地震矩释放均表明2017 年静乐井所在区域应力场在增强。

图8 山西岢岚-山西太原GPS 基线时序曲线 Fig.8 GPS baseline timing curve of Kelan, Shanxi-Taiyuan, Shanxi

静乐井构造上位于山西西部吕梁山隆起区，位于同一构造区距静乐井较近的前兆测点有宁武钻孔应变（60km）和神池钻孔应变（95km），在静乐井水位突升前后这2 个测点的应变观测值也出现大幅度反向变化，所以说静乐井水位高值变化不是独立事件（见图10），可能反映了所在构造区应力场的改变。

图9 静乐井所在区域地震矩加速释放时空图 Fig.9 Accelerated release diagram of seismic moment in jingle well area

综上所述，静乐井水位的高值异常除与降雨、河流荷载作用有关外，井区构造活动增强可能也是一种因素。

图10 静乐井水位与同一构造区其它前兆测项对比图 Fig.10 Comparison between Jingle well water level and other precursor measurements in the same structure area

6 结论与认识

通过数值模拟、井-含水层垂向应力反演、区域应力场变化特征等对静乐井水位2017 年以来出现的高值异常进行初步探讨，得出以下结论：

（1）根据井区周围地壳结构、不同地层介质弹性参数建立三维地质模型，并以2017 年洪水期东碾河水位最大变化量和降雨量作为共同加载体，模拟得到静乐井垂向位移约为0.7mm，垂向应力约为30hPa。

（2）利用静乐井含水层孔隙度、固体骨架体积压缩系数和水的体积压缩系数等参数，根据水平层状含水层模式下含水层介质参数、井水位变化量与含水层垂向应力变化量的关系式，得到静乐井-含水层垂向应力约为79hPa。

（3）数值模拟与井-含水层垂向应力反演结果均表明，静乐井水位高值异常与同时段降雨量增多、河流荷载效应增强关系密切，但2 种结果存在一定差异，含水层垂向应力反演值大于降雨与河流荷载效应，相差39hPa，这种差异可能反映静乐井水位除受降雨与河流荷载作用外，还可能受其他因素的影响。

（4）通过与同一构造区距静乐井较近的前兆测点宁武钻孔应变和神池钻孔应变对比分析可知，静乐井水位突升前后2 个测点应变观测值出现大幅度反向变化，表明静乐井水位高值变化不是独立事件，可能反映了所在构造区应力场发生了改变。

（5）可反映区域应力场变化特征的几种参数对2 种方法计算结果产生差异的原因表明，2017 年静乐井水位高值异常期间山西地震带3 级地震缺震和b 值、穿过静乐井的GPS 基线和地震矩释放均存在显著的异常，推测静乐井水位高值异常除受降雨与河流荷载作用的影响外，也可能受构造活动增强的影响。