2019, Vol. 62
2. 中国科学院大学, 北京 100049;
3. 中国科学院地球科学研究院, 北京 100029;
4. 中国石油大学(华东)地球科学与技术学院, 山东青岛 266580;
5. 中国地震局地球物理研究所, 北京 100081;
6. 中国地震台网中心, 北京 100045
2. University of Chinese Academy of Sciences, Beijing 100049, China;
3. Institutions of Earth Science, Chinese Academy of Sciences, Beijing 100029, China;
4. China University of Petroleum(East China), School of Geosciences, Qingdao Shandong 266580, China;
5. Institute of Geophysics, China Earthquake Administration, Beijing 100081, China;
6. China Earthquake Networks Center, Beijing 100045, China
前人研究中多种类型的地震水力响应被记录(刘耀炜和施锦, 2000; 黄辅琼, 2008; Wang and Manga, 2010; Yan et al., 2016), 很多发生在距离破裂断层较远的区域，这些区域的静态应力变化相对较小(Sil and Freymueller, 2006; Zhang et al., 2017).研究表明地震波的传播能够导致远场区域明显的水力变化(Geballe et al., 2011; Shi and Wang, 2014; Zhang et al., 2015).
岩石的非均匀结构会引起实际储层中不相混融流体的非均匀分布，当地震波通过这种斑块饱和介质传播时，该非均匀分布会产生孔隙压梯度(Ba et al., 2017).孔隙压梯度会驱使其中的流体流动(Ba et al., 2015), 这一现象在低渗透率的砂泥岩中尤为明显(Ba et al., 2016).因此，渗透率会在很大程度上决定流体的运移和流体流动的强度.
前人提出了多种机理来解释同震流体的变化.“地震波传播引起的动态应力导致裂缝中淤积颗粒清空致使渗透率增大”理论被很多研究认可，并且被广泛应用于解释很多远场范围内的同震水位变化(Brodsky et al., 2003; Elkhoury et al., 2006; Wang and Chia, 2008; Wang and Manga, 2010; Shi and Wang, 2014; Shi et al., 2014; Yan et al., 2016).Manga等(2012)得出结论：自然界中应变幅度只要达到10-6就有可能会导致渗透率增大.多孔介质孔隙通道中的毛细淤积被清空也是一种主要的孔隙尺度的解释机理，通过此种作用，弹性波的传播会导致自然渗透率增大(Beresnev et al., 2011).实验室研究表明渗透率的增加及其恢复能够用水流的快速冲刷以及在孔喉或者裂缝间隙的细砂颗粒的逐渐重新淤积来解释(Liu and Manga, 2009; Candela et al., 2014).
Yan等(2014)研究了2011年MW9.1日本东北地震引起的中国区域的地下流体变化，发现地震引起的渗透率的短期变化有可能发生在43%呈现持续性同震水位变化的井孔中，但是在所有观测井孔中发生的几率小于15%.此外此研究还发现在某些观测井中同震相位差下降了，并且统计分析指出同震水位变化与除了潮汐因子外的其他参数(如：井孔震中距、方位角、井深或者含水层岩性)之间并没有明显的关联性.这些研究结果表明远场地震引起的地下流体变化机理是非常复杂的.Wang等(2016)研究发现大地震作用能够打破隔水层从而引起垂向渗透率的增大.Shi和Wang(2016)的研究也表明地震能够重新打开现存垂直向裂隙(之后会恢复)使封闭含水层变为半封闭的状态.上述两个研究中都存在正相位差，表明地震后可能发生了垂直向排水.
左家庄和宝坻是中国华北地区的两口相邻的井(相距约~50 km).两个远场区域(震中距>1000 km)的大地震—2008MW7.9汶川地震和2011 MW9.1日本东北地震引起了这两口井中截然不同的同震水位响应.左家庄井中的同震水位变化呈现非常大的幅度(~2 m)，然而宝坻井中的同震水位变化幅度却非常小(~0.05 m).本研究中，我们采用不同频率域水位的气压响应分析、水位对潮汐应变的响应分析以及井孔测井资料等来解释这两口井同震响应的不一致性.通过测井分析发现宝坻井的观测含水层中存在页岩，并且计算分析发现此井受裂隙的影响很大.页岩的脆性和多裂隙以及各向异性特性可能是导致观测含水层含有页岩的宝坻井受裂隙影响明显的主要原因.计算证实此井存在垂直向渗流，这很可能是由裂隙影响导致宝坻井封闭性不够好引起的.并且气压分析也表明此井伴随较差的井孔贮藏效应，这也可能是由于封闭性不好导致的.综合上述分析，我们将这两口井划分为两种模型—1.水平流动模型; 2.水平流动+垂直流动的混合流动模型.由于裂隙影响，宝坻井(观测含水层含有页岩)所在处的含水层介质与外界的水力沟通性在震前就较强(震前渗透率就比较大)，所以该井观测含水层与外界的孔隙压差别较小，导致该井同震渗透率变化较小甚至没有变化，这是宝坻井同震水位变化幅度总是非常微小的原因.
1 选台原则及基础观测资料研究中我们收集了左家庄(ZJZ)和宝坻(BD)两口井的水位、固体潮(计算理论固体潮数据)和气压观测数据.如图 1所示，三角形为左家庄和宝坻两口水位观测井，左家庄属于北京市，宝坻属于天津市，两者相距仅为50 km.相对于图中2008 MW7.9汶川地震和2011 MW9.1日本东北地震两个远场(震中距>1000 km)大地震而言，可以认为这两口井是相邻的观测井.两口井的同震水位变化幅度差异特别大，且具有历史传承性(见表 1；附图 1).每次左家庄同震水位下降约为2 m，幅度比较大；而每次宝坻同震水位下降幅度仅为0.05 m左右，幅度非常微小.
|
图 1 左家庄(ZJZ)和宝坻(BD)井的地理位置(红色三角形)，以及2008年汶川地震和2011年日本东北地震的震源位置图.“西瓜皮”表示地震的震源机制解(从全球CMT目录中下载).黑线表示中国境内的浅层断层(Deng et al., 2006) Fig. 1 Map showing the geographic locations of the Zuojiazhuang (ZJZ) well and Baodi (BD) well (red filled triangles), and the epicenters of the 2008 MW7.9 Wenchuan earthquake and the 2011 MW9.1 Tohoku earthquake.'Beach ball' shows the lower hemisphere projection of the focal mechanism (from Global CMT) for each earthquake. Black lines indicate the location of surface faults on the Chinese mainland (Deng et al., 2006) |
|
附图 1 两口井在2008汶川地震(a)和2011日本东北地震(b)前后的原始观测水位曲线, 水位中都有清晰可见的潮汐响应记录(承压井).水位记录数值表示井口到水面的距离.虚线表示地震发生的时刻(北京时间) Appendix Fig. 1 Original water-level records of the two wells for the 2008 MW7.9 Wenchuan (a) and the 2011 MW9.0 Tohoku (b) earthquakes, both with resolvable tidal responses. The water level is the distance between the wellhead and the water surface inside the well. The dotted lines indicate the start time of the earthquake (Beijing time) |
|
表 1 左家庄和宝坻井的基础信息 Table 1 Basic information about the Zuojiazhuang and Baodi wells |
由于远场区域静态体应变的变化非常小(Kilb et al., 2002; Zhang and Huang, 2011；Shi et al., 2013)，并且不管处于膨胀区还是压缩区(附图 2)，两口井中的同震水位总是下降，所以此区域的同震水位变化不可能是由静态体应变的变化导致的.二井的同震水力响应机理有待探究.我们选取了2008 MW7.9汶川地震和2011MW9.1日本东北地震两个地震前后的数据进行分析计算，因为此时间段的水位数据并没有受到抽水或者其他影响因素的干扰.需要有足够长时间段的连续的水位记录才能够获得相对稳定的水位与理论固体潮应变之间的相位差信息.前人研究发现，距离海岸几十千米远的井中的水位都会受到较明显的海潮的影响(Beaumont and Berger, 1975)，为了避免海潮的干扰，我们选择研究的这两口井离开海岸的距离大于100 km.表 1概括了井孔的基础信息，包括井深、井径、含水层岩性以及地质结构等.此二井安装的是LN-3A型数字水位记录仪，其观测精度≤0.2% F.S., 采样率为1 sample/min, 分辨率是1 mm.气压观测用的是RTP-1型数字气压记录仪，采样率为1 min.理论固体潮的数据是用Mapseis软件(Lu et al., 2002)计算得到的.
|
附图 2 基于Okada弹性半空间位错模型计算的2008 MW7.9汶川地震和2011 MW9.0日本东北地震引起的同震静态体应变的变化(Okada, 1992).虚线表示压缩区，实线表示膨胀区.图中五角星为震中位置，三角形为左家庄和宝坻水位观测台站位置 Appendix Fig. 2 Spatial distribution of the static volumetric strain changes for the 2008 MW7.9 Wenchuan (a) and the 2011 MW9.0 Tohoku (b) earthquakes, which were calculated using the elastic half-space dislocation model (Okada, 1992). The curved solid lines indicate inflation, whereas the curved dashed lines represent compression. The stars are the epicenters of the earthquakes, and the triangles represent the stations |
通过库伦破裂应力的计算(附图 2)发现，两口井不论是处于压缩区还是膨胀区，同震水位都呈下降，并且远场区域左家庄和宝坻所处位置的同震静态应变变化数值特别微小(表 2)，因此静态应力完全不可能引起如此明显的同震水位变化.
|
图 2 水位和理论固体潮中提取M2波计算所得的相位差和幅度比.幅度比为水位对固体潮应变潮汐M2波的幅度比. 虚线为(a)2008 MW7.9汶川地震的发震时刻, (b) 2011 MW9.1日本东北地震的发震时刻(北京时间) Fig. 2 Amplitude ratios and phase responses over time for the two wells at the frequency of the M2 wave. The amplitude response is the amplitude ratio of water level over the Earth tides. The dotted lines show the start time of the 2008 MW7.9 Wenchuan earthquake and the 2011 MW9.1 Tohoku earthquake (Beijing time) |
|
|
表 2 两个大地震前后水位对潮汐M2波响应的幅度比和相位差 Table 2 Water level amplitude and phase responses to the M2 wave of tidal strain before and after the two large earthquakes |
如附图 3所示，左家庄井深2605 m，观测含水层岩性m，观测含水层主要为页岩.
|
附图 3 两口观测井的井孔柱状测井图 图片来源于中国地震监测志.红色箭头标注部分为井孔观测层开口筛管段的深度. Appendix Fig. 3 Lithologic logs of the two wells. Both records are from the seismic monitoring records of China The red arrows and numbers on the right hand side show the depth of the screened sections. |
水位、气压和潮汐应变三者是相互关联的，天体引潮力的作用导致气压潮加载、海潮加载和固体潮加载，三者共同作用，导致日常水位的周期性变化.潮汐和气压以不同的方式和机理加载在井-含水层系统上，这两种输入信号在某些频段是互相紧密关联的.
我们采用普通的互相关函数，从频率域分析了左家庄和宝坻两口井记录的水位、潮汐和气压三者之间的相干性，计算了相干系数γxy2(附图 4).
|
附图 4 宝坻(a)和左家庄(b)井水位、气压、潮汐应变三者两两之间的互相关函数 Appendix Fig. 4 Ordinary coherence functions among the water level, barometric pressure, and tidal strain for the (a) Baodi and (b) Zuojiazhuang wells |
相干函数的定义:
|
(1) |
式中Gxx(ω)和Gyy(ω)是两个输入信号的功率谱，而Gxy(ω)是两者间的交叉谱.
计算过程中，采用的数据时间段为2010年12月1日至2011年3月10日，数据连续、稳定，并且不受地震或降雨等的影响(附图 1).相干系数计算过程中，选取的汉明窗的窗长和步长分别为30天和15天.与Lai等(2013)的研究类似，我们假定水位受到固体潮影响较大的频段(0.5~8 cpd)为中间频段，f< 0.5 cpd和f>8 cpd分别为低频段和高频段.
如附图 4所示，低频部分(0.1cpd<f < 0.5cpd)左家庄和宝坻井的水位与气压的相干系数都高于0.9，表明水位与气压在周期为数天到几十天的低频部分相关性较好.水位与固体潮在中间频段(~1 cpd/~2cpd)有较为一致的波成分，在日波(O1波:频率0.9295368 cpd)和半日波(M2波:频率1.9322736 cpd)频段水位与潮汐的相干系数接近为1，表明水位和固体潮在O1波和M2波部分的相关性较好.由于气压在M2波和O1波频段的能量很弱，在此二波频段气压与水位的相干系数下降非常明显.在高频部分，水位中气压的信噪比很低，气压信号相对较弱.两口井中水位与气压的相干系数都随着频率的增大而降低.由于M2波的幅度大，几乎不受气压影响，通常是最佳潮汐分析成分的选择.
2.2 潮汐响应分析：相位差和幅度比计算潮汐分析中, M2波的幅度和相位比O1波的更稳定、可信(Bower, 1983).因此，与很多前人的研究类似(晏锐, 2008; Elkhoury et al., 2006; Zhang et al., 2009; Zhang and Huang, 2011; Lai et al., 2013; Yan et al., 2014; Shi et al., 2015; Liao et al., 2015)，本研究也从水位和理论固体潮中提取M2波(周期=745.2 min；角频率=1.4053×10-4 rad/s)来计算两者之间的相位差和幅度比.我们采用互相关分析的方法计算水位和理论固体潮之间的相位差(Zhang et al., 2015, 2017).每15天计算一个数值，每次滑动3天(图 2).
附图 3所示，左家庄是较为普通的井，其主要观测含水层的岩性为白云岩.计算所得两个大地震之前此井的相位差都为负值(图 2).
然而，水位和潮汐应变M2波的正相位差(水位变化导致潮汐体应变的变化)却普遍出现在宝坻井中(图 2).计算还发现，常规岩性(白云岩)观测层的左家庄井，同震相位差明显上升，表明同震渗透率有大幅度的增加(Elkhoury et al., 2006)(图 2).但是，观测含水层主要为页岩的宝坻井同震相位差却并未增大(图 2；表 2).为了选择合理的模型进行分析计算，首先我们需要判断每口井的裂隙的影响并进行储水效应分析.
2.3 裂隙影响分析为了分析是否存在裂隙的影响，我们计算了水位与理论固体潮中提取的M2波的相位差和提取的O1波的相位差(图 3).观测含水层为页岩的宝坻井中，M2波的相位差和O1波的相位差符号相反，因此，正如Bower(1983)研究指出：不排水条件下的此井可看作受到平面裂隙的方向性影响.平面裂隙对井中潮汐的方向性效应影响也将引起水位与理论固体潮中所提取的潮汐分波间的正相位差(0°~180°)(Bower, 1983).因此，观测含水层岩性中含有页岩的宝坻井中正的相位差很有可能是由介质中的裂隙引起的.然而对于观测含水层中不含有页岩的左家庄井，计算所得的水位与理论固体潮中提取的M2波的相位差和提取的O1波的相位差在震前具有相同的符号，并且地震前两者都为负值，表明存在分别由摩擦效应以及井孔贮藏效应导致的排水效应(Bower, 1983).
|
图 3 两口井中水位对潮汐的相位响应分析，水位和理论固体潮中提取M2波计算得到的相位差(红色点)；水位和理论固体潮中提取O1波计算得到的相位差(黑色点).虚线为(a)2008 MW7.9汶川地震的发震时刻; (b) 2011 MW9.1日本东北地震的发震时刻(北京时间) Fig. 3 The phase responses over time for the two wells at the frequency of the M2 wave (red dots) and O1 wave (black dots). The vertical dashed lines show the start time of the (a) 2008 MW7.9 Wenchuan earthquake and (b) 2011 MW9.1 Tohoku earthquake |
为了判断井孔贮藏效应是否同时存在于左家庄和宝坻两口井中，我们计算了水位和气压间的气压系数和相位差(图 4).
|
图 4 交叉谱估算得到的低频-中频-高频部分水位对气压的响应系数和相位差(黑点) (a)宝坻井; (b)左家庄井.负的相位差表示水位滞后于气压 Fig. 4 Barometric efficiencies and phase shifts (black dots) in the low-, intermediate-, and high-frequency bands by cross-spectra estimation, for the Baodi (a) and Zuojiazhuang (b) stations. Negative phase shifts mean water level lags behind the barometric pressure |
在中间频段(0.5~8 cpd)，潮汐的影响是很显著的，因此会导致气压响应不稳定.在频率相对较高的频段(f>6 cpd)，水位的信噪比较低，水位中的气压信号较弱，因此高频段的相位差变化非常杂乱无章.根据常识，在受到噪声或干扰时，相比幅度比(气压系数)相位差可能变化地更快更明显，因此我们可以忽略两口井在较高频段f>6 cpd的相位差分析.
如图 4所示，左家庄的气压系数随着频率从低频段过度到中频段而下降，可能是由井孔贮藏效应导致的(Rojstaczer, 1988a).因此我们能够认为左家庄井存在井孔贮藏效应.相位差在低频段较为稳定并且接近于0度，随着频率从低频段过度到中间频段而下降，这也证明了井孔贮藏效应很可能存在.然而，由于宝坻井低频段和中频段的气压响应差异不是那么明显，宝坻井的井孔贮藏效应不是很明显.因此，总体来说两口井都存在储水效应，但页岩井(宝坻井)中的储水效应不如常规的白云岩(砂岩)井(左家庄井)明显.
详细分析计算方法如下：
水位对气压和潮汐加载响应的传递函数能够通过解下述矩阵计算而获得(Bendat and Piersol, 1986; Rojstaczer, 1988a, b):
|
(2) |
BB和TT分别表示气压和潮汐的功率谱，BT是气压和潮汐之间的交叉谱，TB是BT的复数共轭；BW和TW分别表示气压和水位以及潮汐和水位之间的交叉谱；HB和HT分别为水位对气压和潮汐响应的传递函数.
由方程(2)可以推导得到
|
(3) |
|
(4) |
在分析传递函数过程中，我们采用了类似Lai等(2013)的数据处理方法.在中间频段，根据方程(3)和(4)我们分别计算了水位对气压和潮汐响应的传递函数.在低频和高频部分，潮汐的影响很小，我们根据方程(3)仅计算了水位对气压响应的传递函数.计算过程中，采用的数据时间段为2010年12月1日至2011年3月10日，数据连续、稳定，并且不受地震或降雨等的影响(附图 1).
在低频频段(f < 0.5 cpd)，我们对预处理过的水位和气压数据进行滤取周期1 h至12 d的滤波处理.之后每216 min(约45.5 d)记录一次功率谱和交叉谱的计算结果，汉明窗的窗长和步长分别为~11.4 d和5.7 d.在中间频段(0.5~8 cpd)，滤取的频段为3 min~3 d，每215min(约22.8d)作一次记录，汉明窗的窗长和步长分别为~2.8d和1.4d.在高频段部分(f>8 cpd)，滤取的频段也是3 min~3 d.由于高频段水位响应的信噪比比较低，气压信号相对较弱，我们将数据叠加，每214 min(约11.4 d)作一次记录，汉明窗的窗长和步长分别为~8.5 h和4.3 h.两口井的气压系数和相位差响应如图 4所示.
3 目前存在的两种理想模型 3.1 水平流模型目前为止，Hsieh等(1987)提出的水平流动模型是潮汐响应研究中的标准模型.此模型假设有一均匀、各向同性、水平层状无限延伸、封闭的理想井-含水层系统.在潮汐应力的作用下，含水层发生变形从而压缩孔隙中的流体进入井孔，因而产生水位变化.由于井孔储水效应的存在，使得含水层的水流入井孔需要一定的时间，因此导致含水层孔隙压变化与井孔水位之间存在时间滞后(相位滞后).假设含水层压力水头以大小hf作周期性波动，井孔水位变化记做x，Hsieh等(1987)通过求解地下水扩散方程得出井水位与波动频率之间的表达式：
|
(5) |
|
(6) |
hf是含水层中压力水头的震荡(复数)；ho压力水头震荡的复数；x是水位相对静态位置的位移(复数)；xo是水位震荡幅度的复数；i=(-1)1/2；t表示时间；ω=2π/τ是震荡的角频率; τ是震荡的周期.
幅度比响应A定义为水头和压力水头震荡的幅度比：
|
(7) |
相位差定义为
|
(8) |
公式中arg(z)是复数z的幅角.相位差能够被理解为2πtp/τ, 式中tp是水头震荡达到顶峰到压力头震荡达到顶峰的滞后时间.
井水位的响应依赖于流体流过多孔介质的过程，因此对含水层介质的渗透系数和储水系数都比较敏感.很多情况下其余与井水位响应相关的因素还包括：井的几何结构，震动的周期以及惯性效应.对于长周期的固体潮加载作用惯性效应是可以忽略的(Cooper et al., 1965).另外两个因素是相互独立且较好地被约束的.因此，井水位对于长周期的潮汐振荡加载的幅度比和相位差响应为(Hsieh et al., 1987)
|
(9) |
|
(10) |
式中
|
(11) |
T是导水系数，S是储水系数，Ker、Kei是零阶开尔文函数，rw是井半径，rc是套管内径，ω是潮汐角频率.
此模型伴随有明显的负相位差(-90°~0°)，并且水位对潮汐的幅度和相位响应可以分别用来监测含水层的储水率和渗透率(e.g.Hsieh et al., 1987; Elkhoury et al., 2006; Doan et al., 2006).这种情况下，垂向的排水效应被忽略，相位差一直为负值.
3.2 垂直流模型Wang(2000)研究指出，当存在垂向水流时，水位潮汐响应可用施加于半空间表面的周期性加载来模拟.潜水面的边界条件是排水的，无限深处则为不排水的边界条件.在这些条件下水位对潮汐应力的响应为(Wang, 2000)
|
(12) |
|
(13) |
上式中, 
|
图 5 垂直向排水引起的压力弥散至水表面对应的幅度比和相位差，计算基于潮汐M2波(角频率ω=2pi/(12.42×3600)=1.4053×10-4 rad/s)进行 Fig. 5 Amplitude and phase of due to the diffusion of the pressure to the surface as the vertical drainage, with the M2 wave (angular frequency ω=2 pi/(12.42×3600)=1.4053×10-4 rad/s) |
本研究涉及两种不同的水力响应模型：1.流体水平流动模型；2.流体水平流动+垂直流动混合流动模型.
4.1 水平流动模型——左家庄井(观测含水层无页岩)按照上述分析，左家庄井有明显的负相位差(小于-1°), 此外，还具有显著的井孔贮藏效应, 并且震前该井受裂隙的影响不明显，这些都与Hsieh等(1987)提出的流体水平流动模型相符, 此模型也是目前潮汐响应研究中的标准模型.此模型假设观测含水层为无限空间的单一、均匀、各向同性的封闭含水层.水位和潮汐应变之间的相位差被认为是由井孔贮藏效应导致而不是由含水层中裂隙的排列方向导致(Hsieh et al., 1987).显著的负相位差(-90°~ 0°)(流体流动滞后于潮汐应变的作用)出现在此模型中.相位差的上升(如, -5°至-1°)代表着渗透率的增大(Hsieh et al., 1987; Elkhoury et al., 2006).
计算显示，左家庄的导水系数T处于10-6至10-5m2·s-1范围(图 6).在此范围内，水位对潮汐M2波响应的相位差与导水系数将有明显的同步变化趋势(附图 5; Doan et al., 2006).同震相位差明显上升，与之对应的渗透率也发生了显著增大(图 7a1).
|
附图 5 水位对潮汐响应的幅度比和相位差 此计算是在储水系数S=10-4, rw=rc=10 cm的井中进行的(Doan et al., 2006). Appendix Fig. 5 Predicted amplitude and phase (lags negative) of the response of a reservoir to tides Here a storativity of 10-4 and a well where rw=rc=10 cm was used (Doan et al., 2006). |
宝坻井震前即出现较多正的相位差，可能表示从观测含水层到水面发生了垂直向的流体流动(水力弥散)(Roeloffs, 1996; Lai et al., 2013; Zhang et al., 2016; Shi and Wang, 2016; Wang et al., 2016).因此，Hsieh等(1987)的水平层状流动模型不能够完全适用于宝坻井.该井相位差的误差是3.745°，在扣除误差后(+3.745°/-3.745°)，此井绝大部分相位差仍然在正值和负值之间交替变化，表明宝坻井可以被归为水平流+垂直流的混合流动模型.相位差序列可能是井孔贮藏效应(导致水位M2波相位滞后于潮汐体应变M2波相位)以及垂直向排水效应(导致水位M2波相位超前于潮汐体应变M2波相位)共同作用的结果(Lai et al., 2013; Zhang et al., 2016).正的相位差可能是由页岩的复杂裂隙或者各向异性导致的，表明含水层介质中可能发生了垂直向排水.目前只有纯水平向流体流动模型(第3.1节)和纯垂直向流体流动模型(第3.2节)两种极端的物理模型可用来计算对应的水平向/垂直向渗透率.我们尝试将扣除误差后仍为负相位差的数值对应的相位差和幅度比的数据应用于Hsieh等(1987)的纯水平流动模型来计算水平向的渗透率，但是由于存在正相位差(可能存在垂向流动)，含水层的封闭性较差，所以相位差或许比封闭性好的状态下的数值要高很多，因此计算后得不到合理的渗透率和储水率的数值.所以Hsieh等(1987)的水平层状流动模型对于宝坻井并不完全适用.因此，我们借鉴Elkhoury等(2006)的做法，将所有相位差-15度(去除垂直向渗流效应)后再进行计算，得到了宝坻井水平向的渗透率和储水率(图 7a2、7b2)，此计算忽略了垂直向的渗漏，所得数值与真实数值之间必定是有差别的，但仍然可以在一定程度上反映宝坻井观测含水层的渗透率和储水率.
|
图 6 左家庄井在2008年MW7.9汶川地震前后的导水系数，垂向虚线表示地震的发生时刻(北京时间) Fig. 6 Transmissivity over time from the Zuojiazhuang well before and after the 2008 Wenchuan earthquake. The vertical dashed lines indicate the start time of the earthquake (Beijing time) |
与此同时，我们还用扣除误差后仍为-1°至45°的相位差的数值对应的相位差和幅度比的数据应用于纯垂直流模型(Wang, 2000)来计算垂直向的渗透率(算法见第3.2节)(图 7a3、b3).由计算结果可以确定宝坻井存在垂直向渗流, 该井的封闭性不够好(或许由页岩的脆性、多裂隙影响以及各向异性导致).封闭性不够好也会导致较差的井孔贮藏效应，这个与气压响应反应出宝坻井的井孔贮藏效应不够好是一致的.
或许是由于观测含水层中页岩裂隙的影响，导致宝坻井震前渗透率就比较大(图 7; 表 3)，该井所在位置的观测含水层介质与外界的水力沟通性在震前就较强，所以该井观测含水层介质与外界的孔隙压差别较小，较小的压力梯度导致孔隙中同震流体冲刷颗粒的速度较小，从而导致同震渗透率上升较小甚至没有变化(图 7; 表 3)，这些因素导致宝坻井同震水位变化幅度也较小.此外，宝坻井的封闭性较差，也会加剧同震时候流体的变化幅度较小.
|
图 7 根据均匀水平层状流动模型计算得到的左家庄井的储水率和渗透率(a1)(b1)；根据均匀水平层状流动模型和垂直流模型计算得到的宝坻井的储水率和渗透率(a2)(a3)/(b2)(b3).虚线表示地震的发生时刻(2008汶川地震(A)；2011日本东北地震(B)).空心圆圈表示不适用此模型的数据点，即超出模型阈值范围的数据点 Fig. 7 Specific storage and permeability of well Zuojiazhuang obtained from the horizontal flow model(a1)(b1); and Specific storage and permeability of well Baodi obtained from the horizontal and vertical flow models (a2)(a3)/(b2)(b3). The vertical dashed lines indicate the start time of the (A) 2008 Wenchuan and (B) 2011 Tohoku earthquakes. The open circles are the data invalid for the model and indicating values exceeding the bound. |
|
|
表 3 根据水平层状流动模型计算的两个大地震前后两口井的渗透率和储水率 Table 3 Permeability and specific storage of those two wells before and after the two large earthquakes calculated from the horizontal flow model |
由图 7可见，宝坻井垂直向渗流的渗透率数值比水平向的渗透率数值还要大，这可能与地应力作用导致垂直向的裂隙尺度较大造成的，也有可能存在断层等影响因素，此问题有待下一步结合实地探测地应力等的数据后进行综合论证.
Hsieh等(1987)提出的水平层状流动模型假设有一均匀、各向同性、水平层状无限延伸、封闭的理想井-含水层系统.但实际情况下，含水层介质是非均匀的，并且有可能存在各向异性(尤其是页岩)，封闭性不够好等情况.离开这些理想假设条件的话，基于Hsieh等(1987)模型计算得到的结论就不成立了.
左家庄井中几乎所有的水位与理论固体潮M2波之间的相位差都为负值，除了2011年日本MW9.1东北地震后出现正相位差，表明2011年日本东北地震导致的巨大能量的地震波振动有可能导致了含水层介质性质发生了永久性改变(Zhang et al., 2016).今后的研究中我们将进一步论证此问题.
两口井的井深不同(附图 3)，宝坻井的井深要比左家庄的浅很多.当含水层封闭性很好，Hsieh等(1987)的模型能够完全适用时，井水位的变化与井深没有关系(Hsieh et al., 1987; Doan et al., 2006).而当含水层封闭性较差，观测含水层中的水呈现完全垂直流动时，水位与潮汐应变之间的相位差都为正，此时水位变化与井管开口段的深度有关系(Doan et al., 2006).本研究中，宝坻井水位与潮汐应变的相位差既有正值又有负值，我们认为水平向流体流动和垂直向流体流动可能是同时发生的，因此可以认为观测含水层呈半封闭状态.因此我们可以忽略井孔深度对水位变化的影响.
我们考虑过加入分析其他更多的大地震期间的数据来阐明机理，例如2012 MW8.6苏门答腊地震以及2015MW7.9尼泊尔地震等.但是由于这些地震期间宝坻井的水位观测受到很多干扰并且有较多数据缺失，导致不得不放弃这些数据分析.
近场同震流体变化主要由静态应力变化导致的地层的压缩或膨胀引起而不是由渗透率的变化引起(Wang and Chia, 2008; Zhang and Huang, 2011; Shi et al., 2013).由于近场同震流体变化的机理与中-远场同震流体变化的机理截然不同，本研究中我们并未对近场地震或者区域地震作研究.
数值模拟方法有可能对精确机理的探究有帮助，以后将进一步研究探讨.未来的工作中亟待建立完善的水平流动+垂直流动的混合流模型，并进行模型测试、应用.
6 结论不管是国内还是国外，前人对地震地下流体机理的研究基本都是基于砂岩等常规岩性的观测含水层进行的，以页岩为观测含水层进行同震流体和气压、潮汐响应机理研究目前尚是空白.页岩具有低孔-低渗、各向异性、脆性并且含有大量复杂裂隙等特性.井孔打钻时，大多常规井的观测含水层一般都会避开页岩层，因此此类观测井数量非常有限.笔者基于宝坻和左家庄两口相邻的井做了一些探索性的对比研究.对于2008 MW7.9汶川地震和2011 MW9.1日本东北地震来说，这两口井都处于远场区域.我们运用了流体对潮汐以及气压的响应来分析解释观测含水层含有页岩的井孔的一系列特殊的同震变化现象.研究发现：1)页岩的脆性和多裂隙以及各向异性特性可能导致观测含水层含有页岩的宝坻井受裂隙影响明显，导致其封闭性不够好(伴随较差的井孔贮藏效应)，存在垂直向渗流；2)由于裂隙影响，宝坻井(观测含水层含有页岩的井)所在位置的孔隙压与外界的水力沟通性在震前就较强(震前渗透率就比较大)，所以该井观测含水层与外界的孔隙压差别较小，导致同震渗透率上升较小甚至没有变化，最终导致该井同震水位变化幅度总是非常微小.3)通过多方面分析计算论证，本研究所涉及的两口相邻井可分别归为两个模型—水平流动模型以及水平流+垂直流的混合流动模型.
致谢 我们感谢晏锐研究员、王其允教授、史浙明副教授以及张文辉博士给予的帮助和讨论.
Ba J, Carcione J M, Sun W T. 2015. Seismic attenuation due to heterogeneities of rock fabric and fluid distribution. Geophysical Journal International, 202(3): 1843-1847. DOI:10.1093/gji/ggv255 |
Ba J, Zhao J G, Carcione J M, et al. 2016. Compressional wave dispersion due to rock matrix stiffening by clay squirt flow. Geophysical Research Letters, 43(12): 6186-6195. DOI:10.1002/2016GL069312 |
Ba J, Xu W H, Fu L Y, et al. 2017. Rock anelasticity due to patchy saturation and fabric heterogeneity:A double double-porosity model of wave propagation. Journal of Geophysical Research:Solid Earth, 122(3): 1949-1976. |
Beaumont C, Berger J. 1975. An analysis of tidal strain observations from the United States of America:Ⅰ. The laterally homogeneous tide. Bulletin of the Seismological Society of America, 66: 1613-1629. |
Bendat J S, Piersol A G. 1986. Random Data: Analysis and Measurement Procedures. New York: John Wiley, 640.
|
Beresnev I, Gaul W, Vigil R D. 2011. Direct pore-level observation of permeability increase in two-phase flow by shaking. Geophysical Research Letters, 38: L20302. DOI:10.1029/2011GL048840 |
Bower D R. 1983. Bedrock fracture parameters from the interpretation of well tides. Journal of Geophysical Research, 88(B6): 5025-5035. DOI:10.1029/JB088iB06p05025 |
Brodsky E E, Roeloffs E, Woodcock D, et al. 2003. A mechanism for sustained groundwater pressure changes induced by distant earthquakes. Journal of Geophysical Research:Solid Earth, 108(B8): 2390. DOI:10.1029/2002JB002321 |
Candela T, Brodsky E E, Marone C, et al. 2014. Laboratory evidence for particle mobilization as a mechanism for permeability enhancement via dynamic stressing. Earth and Planetary Science Letters, 392: 279-291. DOI:10.1016/j.epsl.2014.02.025 |
Cooper H H Jr, Bredhoeft J D, Papdopulos I S, et al. 1965. The response of well-aquifer systems to seismic waves. Journal of Geophysical Research, 70(16): 3915-3926. DOI:10.1029/JZ070i016p03915 |
Deng Q D, Zhang P Z, Ran Y K. 2006. Distribution of Active Faults in China (1:4000000). Beijing: Science Press.
|
Doan M L, Brodsky E E, Prioul R, et al. 2006. Tidal Analysis of Borehole Pressure:A Tutorial. California: University of California.
|
Elkhoury J E, Brodsky E E, Agnew D C. 2006. Seismic waves increase permeability. Nature, 441(7097): 1135-1138. DOI:10.1038/nature04798 |
Geballe Z M, Wang C Y, Manga M. 2011. A permeability-change model for water-level changes triggered by teleseismic waves. Geofluids, 11(3): 302-308. DOI:10.1111/gfl.2011.11.issue-3 |
Hsieh P A, Bredehoeft J D, Farr J M. 1987. Determination of aquifer transmissivity from earth tide analysis. Water Resources Research, 23(10): 1824-1832. DOI:10.1029/WR023i010p01824 |
Huang FQ. 2008. Response of wells in groundwater monitoring network in Chinese Mainland to recent large earthquakes[Ph (in Chinese). Beijing: Institute of Geophysics, China Earthquake Administration.
|
Kilb D, Gomberg J, Bodin P. 2002. Aftershock triggering by complete Coulomb stress changes. Journal of Geophysical Research, 107(B4): ESE 2-1-ESE 2-14. DOI:10.1029/2001JB000202 |
Lai G J, Ge H K, Wang W W. 2013. Transfer functions of the well-aquifer systems response to atmospheric loading and Earth tide from low to high-frequency band. Journal of Geophysical Research, 118(5): 1904-1924. |
Liao X, Wang C Y, Liu C P. 2015. Disruption of groundwater systems by earthquakes. Geophysical Research Letters, 42(22): 9758-9763. DOI:10.1002/2015GL066394 |
Liu W Q, Manga M. 2009. Changes in permeability caused by dynamic stresses in fractured sandstone. Geophysical Research Letters, 36: L20307. DOI:10.1029/2009GL039852 |
Liu YW, Shi J. 2000. Information characteristics of ground fluid precursor on continental strong earthquakes. Acta Seismologica Sinica (in Chinese), 22(1): 102-107. |
Lu Y Z, Li S L, Deng Z H, et al. 2002. Seismology Analysis and Prediction System Based on GIS (Mapsis Software). Chengdu: Chengdu Map Press: 232.
|
Manga M, Beresnev I, Brodsky E E, et al. 2012. Changes in permeability caused by transient stresses:Field observations, experiments, and mechanisms. Reviews of Geophysics, 50(2): RG2004. DOI:10.1029/2011RG000382 |
Okada Y. 1992. Internal deformation due to shear and tensile faults in a half-space. Bulletin of the Seismological Society of America, 82(2): 1018-1040. |
Roeloffs E. 1996. Poroelastic techniques in the study of earthquake-related hydrologic phenomena. Advances in Geophysics, 37: 135-195. DOI:10.1016/S0065-2687(08)60270-8 |
Rojstaczer S. 1988a. Determination of fluid flow properties from the response of water levels in wells to atmospheric loading. Water Resources Research, 24(11): 1927-1938. |
Rojstaczer S. 1988b. Intermediate period response of water levels in wells to crustal strain:sensitivity and noise level. Journal of Geophysical Research, 93(B11): 13619-13634. DOI:10.1029/JB093iB11p13619 |
Shi Z M, Wang G C, Liu C L. 2013. Co-seismic groundwater level changes induced by the May 12, 2008 Wenchuan earthquake in the near field. Pure and Applied Geophysics, 170(11): 1773-1783. DOI:10.1007/s00024-012-0606-1 |
Shi Z M, Wang G C. 2014. Hydrological response to multiple large distant earthquakes in the Mile well, China. Journal of Geophysical Research:Earth Surface, 119(11): 2448-2459. DOI:10.1002/2014JF003184 |
Shi Z M, Wang G C, Wang C Y, et al. 2014. Comparison of hydrological responses to the Wenchuan and Lushan earthquakes. Earth and Planetary Science Letters, 391: 193-200. DOI:10.1016/j.epsl.2014.01.048 |
Shi Z M, Wang G C, Manga M, et al. 2015. Mechanism of co-seismic water level change following four great earthquakes-insights from co-seismic responses throughout the Chinese mainland. Earth and Planetary Science Letters, 430: 66-74. DOI:10.1016/j.epsl.2015.08.012 |
Shi Z M, Wang G C. 2016. Aquifers switched from confined to semiconfined by earthquakes. Geophysical Research Letters, 43(21): 11166-11172. DOI:10.1002/2016GL070937 |
Sil S, Freymueller J T. 2006. Well water level changes in Fairbanks, Alaska, due to the great Sumatra-Andaman earthquake. Earth, Planets and Space, 58(2): 181-184. |
Wang C Y, Chia Y. 2008. Mechanism of water level changes during earthquakes:Near field versus intermediate field. Geophysical Research Letters, 35(12): L12402. DOI:10.1029/2008GL034227 |
Wang C Y, Manga M. 2010. Earthquakes and Water:Lecture Notes in Earth Sciences. Berlin: Springer: 18-24.
|
Wang C Y, Liao X, Wang L P, et al. 2016. Large earthquakes create vertical permeability by breaching aquitards. Water Resources Research, 52(8): 5923-5937. DOI:10.1002/2016WR018893 |
Wang H F. 2000. Theory of Linear Poroelasticity with Applications to Geomechanics and Hydrogeology. Princeton: Princeton University Press..
|
Yan R, Woith H, Wang R J. 2014. Groundwater level changes induced by the 2011 Tohoku earthquake in China mainland. Geophysical Journal International, 199(1): 533-548. DOI:10.1093/gji/ggu196 |
Yan R, Wang G C, Shi Z M. 2016. Sensitivity of hydraulic properties to dynamic strain within a fault damage zone. Journal of Hydrology, 543: 721-728. DOI:10.1016/j.jhydrol.2016.10.043 |
Zhang Y, Huang F Q, Lai G J. 2009. Research on Skempton's coefficient B based on the observation of groundwater of Changping station. Earthquake Science, 22(6): 631-638. DOI:10.1007/s11589-009-0631-z |
Zhang Y, Huang F Q. 2011. Mechanism of different coseismic water-level changes in wells with similar epicentral distances of intermediate field. Bulletin of the Seismological Society of America, 101(4): 1531-1541. DOI:10.1785/0120100104 |
Zhang Y, Fu L Y, Huang F Q, et al. 2015. Coseismic water-level changes in a well induced by teleseismic waves from three large earthquakes. Tectonophysics, 651-652: 232-241. DOI:10.1016/j.tecto.2015.02.027 |
Zhang Y, Fu L Y, Ma Y C, et al. 2016. Different hydraulic responses to the 2008 Wenchuan and 2011 Tohoku earthquakes in two adjacent far-field wells:the effect of shales on aquifer lithology. Earth, Planets and Space, 68(1): 178. DOI:10.1186/s40623-016-0555-5 |
Zhang Y, Wang C Y, Fu L Y, et al. 2017. Mechanism of the coseismic change of volumetric strain in the far field of earthquakes. Bulletin of the Seismological Society of America, 107(1): 475-481. DOI:10.1785/0120160253 |
黄辅琼. 2008.中国大陆地震地下水观测井对大地震的响应[博士论文].北京: 中国地震局地球物理研究所.
|
刘耀炜, 施锦. 2000. 强震地下流体前兆信息特征. 地震学报, 22(1): 102-107. DOI:10.3321/j.issn:0253-3782.2000.01.014 |
晏锐. 2008.影响井水位变化的几种因素研究[硕士论文].北京: 中国地震局地震预测研究所.
|