1976年唐山7.8级地震前地下水位变化与地震成核过程的关系

马玉川， 晏锐， 王广才， 余怀忠， 黎明晓， 丁志华， 张子广

1 中国地震台网中心， 北京 100045 2 中国地质大学(北京)水资源与环境学院， 北京 100083 3 河北省地震局， 石家庄 050021

0 引言

地震预测是尚未解决的世界性科学难题，而探索与地震有必然联系的地震前兆是亟需解决的科学问题之一.近几十年来，根据构造物理实验，学者们观测到断层在黏滑失稳前存在预滑和预滑区逐步扩展现象(Ohnaka and Kuwahara，1990；Ohnaka and Shen，1999；Mclaskey and Kilgore，2013)，以及断层在黏滑后期存在由应力积累为主转变为应力释放为主的不可逆变形的亚失稳态阶段(马瑾等，2012；马瑾和郭彦双，2014；马瑾，2016)，这些地震成核过程的研究为探索确定性地震前兆提供了实验基础，对开展地震成核理论研究和预测实践具有重要的现实指导意义.基于速率-状态依赖性摩擦本构关系(Dieterich，1978, 1979, 1981；Ruina，1983)，学者们对地震成核过程进行了数值模拟研究，结果表明在地震成核过程中存在成核尺度与加载速率关联、滑动速率持续增加、成核区不断扩展等现象(Dieterich, 1992；Kato et al., 1992；Kato and Hirasawa, 1996；He et al.，1998, 2021；何昌荣，2000； Ampuero and Rubin, 2008).在野外观测条件下，存在较为普遍的由地震成核前应力场时空演变所导致的大地震之前的小震时空迁移或错动加速现象(Chen et al., 2011; Kato et al., 2012；Kato and Nakagawa, 2014；马瑾和郭彦双，2014；Yagi et al., 2014； Zhuo et al., 2018；王凯英等，2021).由于野外构造的复杂性和野外观测的局限性，在野外利用其他地球物理监测手段(形变、电磁和流体等)辨认地震成核过程，特别是地震成核的后期阶段，仍然是需要进一步研究的科学问题.

1976年7月28日3时42分，在河北省唐山市丰南一带发生7.8级大地震，造成惨重的人员伤亡和经济损失.唐山地震发生时，京津唐及其邻近地区有400余口地下水位观测井，震中区有13口观测井，记录到了丰富的观测资料(唐山地震资料汇编组，1981).这些资料为探索野外地震成核过程提供了可能.本项研究对唐山地震前后发震断层附近的地下水位变化进行回溯分析，通过对比分析唐山地震前的地下水位变化与地震成核过程研究结果，以探索在野外观测中辨认地震成核过程.

1 区域地质构造

1976年唐山7.8级地震震中区位于唐山断陷中的开平向斜构造的西北翼，区内发育有北东走向的唐山断裂带(图1)，该断裂带主要由三条平行的断裂组成，自西向东分别为陡河断裂、长山—巍山断裂和唐山—古冶断裂(虢顺民等，1977；刘国栋等，1982；Liu et al.， 2013).陡河断裂全长约50 km，唐山市以南一段被开滦煤矿称为山Ⅰ断层，为正断层，断层面倾角50°～90°，断层总位移150～200 m，附近地层剧烈褶皱甚至倒转；长山—巍山断裂全长约20 km，被称为山Ⅱ和山Ⅲ断层，为逆断层，断层倾角随深度加深而变缓，500 m以上30°～60°，500 m以下小于15°，最大断距超过300 m；唐山—古冶断裂全长约30 km，被称为山Ⅳ和山Ⅴ断层，山Ⅳ断层为正断层，倾角60°，山Ⅴ断层为逆断层，倾角45°，垂直断距大于500 m(陈墨香，1981；刘国栋等，1982；谷世铭和周全度，1985).

图1 唐山7.8级地震震中区地质略图

震后的野外地质调查工作表明，唐山地震的发震断层与唐山断裂带中的唐山—古冶断裂有关.唐山地震形成的地表破裂带主要分布在唐山—古冶断裂内，北起唐山市胜利路南至安机寨附近，长度8～11 km，走向为北东向30°，具有右旋走滑特征，最大走滑量为2.3 m，多数地段垂直位移为0.2～0.7 m(虢顺民等，1977；杨理华，1982；杜春涛等，1985).近些年，有些学者基于地表塌陷和形变资料、钻孔剖面和探槽探测等工作，认为唐山地震的地表破裂由安机寨向南延伸至丰南西河以南，总长度超过47 km(邱泽华等，2005；江娃利，2006；Guo et al.，2017).

2 地下水位资料

唐山7.8级地震发生后，许多研究者对京津唐地区地下水位资料进行了系统的梳理.结果表明，唐山地震前存在显著的地下水位长期、中期、短期和临震变化(Wang et al.，1979；王雅灵，1980；吴锦秀等，1980；张洪波和李坤，1981；汪成民和尹伯忠，1982).研究者分析了这些变化和唐山地震的关系，对于1972年初至1975年春夏唐山和天津一带出现的水位长期下降变化，王雅灵(1980)认为天津地区封存咸水水位的变化是唐山地震的趋势异常，汪成民和尹伯忠(1982)认为不能完全归因于抽水和干旱，可能与唐山地震的孕育过程有关，而张洪波和李坤(1981)、贾化周等(1982)、车用太等(1994)认为是地下水过量开采所致.对于1975年春夏至1976年春唐山地区出现的中期水位下降速率变缓，以及唐山和天津一带1976年4、5月份出现的短期水位加速下降、1976年7月下旬出现的临震水位突变，无法用地下水的开采和补给解释，被认为是唐山地震的前兆异常(王雅灵，1980；张洪波和李坤，1981；贾化周等，1982；汪成民和尹伯忠，1982；车用太等，1994).

基于上述唐山7.8级地震前地下水位变化的相关研究，本研究利用文献(唐山地震资料汇编组，1981；杨文田和虞修，1981；张洪波和李坤，1981；汪成民和尹伯忠，1982；谷世铭和周全度，1985)中的地下水位资料，回溯了唐山地震前后的地下水位变化.因为1972年初至1975年春夏唐山和天津一带出现的地下水位长期下降变化与唐山地震的关系存在争论，所以重点收集了唐山地震前的中期、短期、临震和震后地下水位资料.由于地震成核过程重点在于震中区附近的变化，因而主要收集了唐山地震发震断层附近的地下水位资料.共收集到13口地下水位观测井，基本情况见表1、表2和图1.这些观测点分布在唐山地震地表破裂带以西至陡河断裂，距离震中(东经118.2°，北纬39.6°，深度12 km)3～18 km.

表1 唐山7.8级地震发震断层附近地下水位震前中短期变化和震后变化信息

表2 唐山7.8级地震前发震断层附近地下水位临震变化信息

3 唐山地震前后发震断层附近的地下水位变化

有7个测点记录到唐山7.8级地震前的中期和短期水位变化.中期变化出现在唐山地震前1年左右(1975年7月15日)，变化形态为同步的地下水位下降速率变缓，水位波动幅度不大(表1).短期变化出现在唐山地震前4个月左右(1976年4月5日)，变化形态为同步的地下水位加速下降，到地震发生前达到历史低值(表1，图2)，下降幅度超过1 m(唐山矿涌水量为4 m3·min-1).在这些地下水位出现中期和短期变化期间(1975—1976年)，唐山地区的地下水开采强度维持在50 m3·min-1左右(张洪波和李坤，1981).这些测点的基础资料比较完整，震前水位观测连续，观测层岩性以埋藏较深的奥陶系灰岩为主，山西水2井和丰南岳42井的水位观测为自记式水位仪，其他为每天一次的人工测量，观测精度为厘米级.

图2 唐山7.8级地震发生前后发震断层附近地下水位变化曲线

有6个测点记录到唐山地震前的临震变化.临震变化出现在唐山地震前16 h至几分钟，地下水位出现数米的突升或波动变化，有些甚至形成自喷和自流(表2)，最为显著的是丰南岳42井，在震前2 h，井水自埋深10 m以下喷出地表(由于没有收集到观测数据或曲线，这些信息来自汪成民和尹伯忠(1982)的文字描述).

有8个测点记录到唐山地震后的大幅上升，包含2口震前无变化的井孔(唐山水文站).震后水位上升的幅度为1.8～10 m(唐山矿涌水量在震后上升56.5 m3·min-1)，上升后维持在高值，1978年开始下降，1979年底恢复到震前水平(表1、图2).由于只收集到震后恢复观测后的数据，实际的变化幅度可能更大.例如山西水2井，该井的地下水位埋深在地震前为32 m左右，震前在水面以下2 m处安装了测温仪，震后测温仪被抛出井口(汪成民等，1988).对于同时记录到震后和临震变化的唐山矿和丰南岳42井，震后变化形态是临震变化的延续.

4 地下水位变化与地震成核研究结果的对比

基于构造物理实验和数值模研究，尤其在Dieterich(1978，1979，1981)和Ruina(1983)建立和完善了速率-状态依赖性摩擦本构关系之后，地震成核过程的相关研究取得了一系列重要的认识，包括发现了成核区的扩展现象、注意到了滑动速率的持续增加等(Dieterich, 1992； Kato et al., 1992；Kato and Hirasawa,1996；何昌荣，2000；Ampuero and Rubin, 2008)；然而，目前对地震成核过程各个阶段的划分、最终成核阶段的空间尺度大小和持续时间、以及最终成核阶段能否在野外条件下被检测到等问题仍然在研究当中(Zhang et al.，2020； Barbour and Beeler，2021；He et el.，2021).针对上述问题，基于花岗岩实验约束的速率-状态摩擦定律，He等(2021)对地震成核过程进行了新一轮的数值模拟研究，得到了一些新的认识，部分研究结果如图3所示.其中，地震成核过程包括应力积累、初始弱化和总体弱化三个阶段(图3a)；在地震成核后期的总体弱化阶段，会出现核心弱化区的现象，核心弱化区的边界具有随时间扩展并在后期收缩的特征(图3b)；核心弱化区的边界扩展开始时间、边界扩展速率、边界收缩时间受控于加载速率和断层尺度(图3c)；核心弱化区边界收缩的最后阶段，弱化区越来越小，弱化区内应力下降、滑动量加速增加，外围应力持续升高(图3d).为探寻唐山地震前地下水位是否是与野外地震成核过程有关，这里将He等(2021)的地震成核模拟结果与唐山地震前的地下水位观测结果进行对比.

图3 地震成核过程二维数值模拟研究结果(来自He et al.，2021)

根据He等(2021)的研究结果，地震成核各个阶段的形成、演化与断层加载速率(V0)和断层尺度(L)密切相关(图3).唐山地区的加载速率参考地壳形变和重复地震等资料进行确定，唐山地震的断层尺度参考地表破裂带和余震精定位结果来确定.对于唐山地震发生前震中区附近的地壳形变，张祖胜(1982)基于水准资料总结了唐山地震震中区附近的垂直形变特征：1954—1967年略有上升，平均变化速率约0.2 mm·a-1，1967—1969年上升，速率为8.8 mm·a-1，1969—1972年下降，速率为8.3 mm·a-1，1972—1975年略有上升，速率为2.7 mm·a-1；黄立人等(1988)根据唐山地区的三角测量结果给出唐山地区在1960—1971年期间相对周围块体发生北北东向运动，水平位移约153 mm，而该水平形变测量的精度较低，仅能够反映震前存在一定的水平运动.受观测数据的限制，目前无法获取唐山地震前较为精确的地壳运动量，所以主要采用唐山地震后的观测数据来估算加载速率.对于唐山地震发生后震中区附近的地壳形变，Li等(2007)依据精定位方法获取了重复地震，由此估算出唐山地震断层滑动速率为2.6 mm·a-1，与华北地区GPS测量结果相当(Wang et al., 2001)；武艳强等(2016)依据1999—2007年唐山地区的GPS资料得到地壳运动量小于3.6 mm·a-1；根据Wang和Shen(2020)给出的全国1991—2016年GPS解算结果，唐山地震震中区附近(东经117.2°—119.2°、北纬38.6°—40.6°)的位移量介于2.5～5 mm·a-1；上述震后结果基本在同一数量级，本文采用Li等(2007)基于重复地震所得的结果，将唐山地区的加载速率设定为2.6 mm·a-1.对于唐山地震的断层尺度，早期野外调查所得的唐山地震地表破裂带长度为8～11 km(虢顺民等，1977；杨理华，1982；杜春涛等，1985)；近些年有学者根据新的观测证据认为唐山地震地表破裂带长度超过47 km(邱泽华等，2005；江娃利，2006；Guo et al.，2017)；Li等(2007)对唐山地震震中区2001—2006年的小震进行了精定位，精定位后的地震分布在50 km范围内(图1)；参考这些研究结果，本文将唐山地震的断层尺度设定为50 km.

根据He等(2021)的模拟结果，如果将断层加载速率(V0)设定为2.6 mm·a-1、断层尺度(L)设定为50 km，由此得到核心弱化区边界扩展的开始时间 (Tex)为628天、核心弱化区边界开始收缩的时间为377天.对比唐山地震前地下水位的变化时间，唐山地震前地下水位中期变化开始时间(379天)与核心弱化区边界收缩的开始时间吻合.在核心弱化区继续收缩阶段，其外围应力持续升高，可能对孔隙水进行挤压，同时，断层滑移速率和滑移量出现加速(He et al., 2021)；在这样的情况下，根据Segall和Rice(1995)基于模拟实验得出的结果，如果考虑断层带的剪切变形可能引起的扩容作用，核心弱化区的断层带孔隙度就会增加；基于这些研究结果，唐山地震前地下水位短期变化(114天)可能与核心弱化区扩容有关.唐山地震临震水位变化的时间与模拟结果预期的成核过程最后阶段其核心弱化区周围应力快速上升的时间基本一致；成核区周围应力的快速上升，可以引起成核区附近水位的大幅上升，与临震阶段观测到的情况基本一致.

从空间分布的角度，根据He等(2021)对断层加载速率V0=2.6 mm·a-1、断层尺度L=50 km条件下的模拟结果，核心弱化区沿断层走向的长度大约为23 km；本研究中的水位观测井孔分布在沿发震断层方向上距离震中25 km范围之内；因而水位观测井分布在临近核心弱化区的区域.依据He等(2021)的研究结果，在地震成核过程中，核心弱化区在形成之后边界开始收缩，最后集中在空间尺度较小的成核区(图3)；水位观测井的空间分布，也存在一个类似的收缩过程，即中短期变化分布在发震断层外围，而临震变化主要集中在发震断层附近；2个同时记录到中短期和临震变化的水位测点(唐山矿和丰南岳42井)，分布在与唐山地震主破裂带较近的位置(图1).

因此，从时间变化和空间分布的对比结果来看，唐山地震前的地下水位变化与地震成核过程模拟结果具有一定的可比性，即中期水位下降速率变缓与地震成核后期阶段的核心弱化区收缩相关，而短期水位加速下降可能对应于核心弱化区断层带扩容，最后临震阶段的大幅上升与成核区周围应力快速上升基本吻合.

5 基于地震成核理论讨论地下水位变化的成因

由于唐山地震前的地下水位变化可以与地震成核过程模拟研究结果进行类比，因此，以下从地震成核过程中应力变化的影响出发，对唐山地震前地下水位变化的成因给出一种可能的解释.对于地下水位中期下降速率变缓，其可能的成因是由于核心弱化区的收缩，造成核心弱化区内部应力显著降低，而核心弱化区外围应力显著增强；因为地下水位观测井临近断层滑动面，经过一段时间的收缩，大部分观测井都会处于核心弱化区外围的应力增强区，所以在核心弱化区收缩时水位出现上升变化(下降背景下出现速率变缓).对于短期水位下降加速，其可能的成因是由于核心弱化区的扩容，形成扩容吸水现象，引起断层带周围地下水量减小；因为短期异常水位观测井临近强扩容区，所以出现水位下降现象(下降背景下出现速率加速).由于核心弱化区进一步收缩直至成核，成核区外部是应力增强区，在成核区周围的应力处于快速上升阶段，震源区水位会出现巨幅上升(或波动后上升).

许多研究者曾对唐山地震前地下水位异常的成因进行过讨论，例如：王雅灵(1980)认为是区域应力场作用下震源及附近地区应力集中与强化过程的反映；吴锦秀等(1980)将中长期、短期、临震、震时变化对应为震源体的弹性变形、塑形变形、微裂扩容和主破裂阶段；张洪波和李坤(1981)认为与震源应力场的强烈挤压和断裂面的急速张合有关；汪成民和尹伯忠(1982)认为可能反映了震中地区裂隙演变的过程；车用太等(1994)推测是深部物质上涌-微裂隙发育的结果.本文的解释与这些研究结果相比，同属于应力场变化引起介质演化进而导致地下水位变化的范畴，相对而言，地下水位中期和临震变化的机理与地震成核理论吻合的较好；然而，本文的解释和以往的研究都属于定性化的讨论，没有对地下水位的变化给出定量化的分析.

近年来，马瑾等(2012)，马瑾和郭彦双(2014)，马瑾(2016)研究团队提出了亚失稳理论，认为断层在黏滑失稳前存在一个亚失稳阶段，在这个阶段内，应力由积累为主转变为应力释放为主，是断层进入不可逆变形的关键阶段.亚失稳可细分为前期的准静态应变释放和后期的准动态应变释放两个阶段.依据马瑾等(2012)实验的温度场演化结果(应力变化的反映)，拐折断层的拐折点及附近区域在亚失稳阶段的应力变化明显(图4a).唐山地震前地下水位出现中期、短期、临震变化的井孔也分布在断层出现拐折的震中附近(图4b).因此，地下水位变化与拐折断层亚失稳实验结果在空间分布特征上有一定的相似性.

图4 亚失稳实验中的断层和温度测点、唐山7.8级地震的发震断层和附近的水位测点(图4a来自马瑾等(2012))

6 结论与讨论

根据1976年唐山7.8级地震前地下水位变化的相关研究，以及地震成核过程的最新研究成果，本文回溯了唐山7.8级地震发生前后发震断层附近的地下水位变化，对比分析了地下水位变化与地震成核研究结果，取得以下初步结论：

唐山7.8级地震前后地下水位存在显著的中期、短期、临震和震后变化，其中，中期地下水位下降速率变缓与地震成核后期阶段的核心弱化区收缩、短期地下水加速下降与核心弱化区扩容、临震水位突升或波动与成核区周围应力快速上升阶段在时间变化和空间分布特征上具有一定的相关性.

唐山7.8级地震的临震水位变化，可以用成核区外围应力的快速增加来解释，核心弱化区收缩引起的应力变化以及扩容吸水基本上能够解释中期和短期水位变化，因此最新的地震成核理论可以较为合理的解释唐山地震前地下水位变化的成因.据此认为唐山地震前地下水位变化可能反映的是野外观测到的地震成核后期阶段，在一定程度上为地震成核野外观测提供了实际震例.

虽然本研究仍是回溯性的个案研究，但研究结果可能反映了地震成核过程中的一种情况，对探索井水位地震前兆、地震成核的野外观测以及提高地震监测预测水平提供了可能.但由于仍是初步研究，对于唐山地震前地下水位变化与地震成核过程定量化的关系、唐山地震前地下水位长期变化与地震成核过程其他阶段的关系、唐山地震后地下水位巨幅上升的机理，目前没有给出科学的分析，有待进一步研究.

致谢在研究中得到蒋海昆研究员、车用太研究员、付虹研究员、李乐研究员、刘亢副研究员、姚琪副研究员的指教和帮助，审稿专家提供了宝贵的修改意见和建议，在此深表感谢.