腾冲井同震水位动态孔弹性响应的渗透率反演_高楚林.pdf

2023 年第 38 卷 第1期2023，38(1):0012-0019地球物理学进展Progress in Geophysicshttp:/wwwprogeophyscnISSN 1004-2903CN 11-2982/P高楚林，张煜，张双喜，等 2023 腾冲井同震水位动态孔弹性响应的渗透率反演 地球物理学进展，38(1):0012-0019，doi:10 6038/pg2023GG0178GAO ChuLin，ZHANG Yu，ZHANG ShuangXi，et al2023 Permeability inversion of coseismic water level dynamic poroelastic response in Tengchong wellProgress in Geophysics(in Chinese)，38(1):0012-0019，doi:106038/pg2023GG0178腾冲井同震水位动态孔弹性响应的渗透率反演Permeability inversion of coseismic water level dynamic poroelastic response inTengchong well高楚林1，张煜1，2*，张双喜1，2，黄旭东1，姚永康1GAO ChuLin1，ZHANG Yu1，2*，ZHANG ShuangXi1，2，HUANG XuDong1，YAO YongKang1收稿日期2022-06-07;修回日期2022-11-15投稿网址http:/www progeophys cn基金项 目国 家 自 然 科 学 基 金 项 目(41774119，41874169，41774056，42074176，U1939204)、中 央 高 校 基 础 研 究 专 项 基 金 项 目(2042017kf0228)和武汉大学实验技术专项基金项目(WHU-2019SYJS-12)联合资助第一作者简介高楚林，女，1997 年生，硕士研究生，主要从事井水位同震响应机制方面研究 E-mail:clgao whu edu cn*通讯作者张煜，男，1982 年生，副教授，主要从事地震波传播与成像研究 E-mail:yuzhang sgg whu edu cn1 武汉大学测绘学院，武汉4300792 武汉大学地球空间环境与大地测量教育部重点实验室，武汉4300791 School of Geodesy and Geomatics，Wuhan University，Wuhan 430079，China2 Key Laboratory of Geospace Environment and Geodesy，Ministry of Education，Wuhan University，Wuhan 430079，China摘要获得地下含水层的水文参数对于适当评价地下水资源至关重要 除了实验室方法和现场测井实验，目前多利用水位对外部周期性载荷的响应函数进行含水层参数的反演 其中，基于似稳态孔弹性理论的方法有很多，但这些方法通常忽略多孔含水层中可能存在的频率依赖的固流耦合 本文提出基于动态孔弹性理论的水文同震响应模型，推演了频率依赖的含水层孔压和围压的关系，构建动态响应函数 根据云南省腾冲井观测到的近场同震水位响应并基于动态模型反演了井含水层渗透率的长期演化 计算结果显示动态孔弹性模型能很好地解释水位的同震响应并且适当反演含水层参数 近场地震会导致井含水层的短期渗透率大幅度增加并在不久后恢复，这可能与地震触发应力导致的瞬时固体形变、孔隙疏通和颗粒物的重分布有关关键词井水位;同震响应;动态孔弹性理论中图分类号P315文献标识码Adoi:10 6038/pg2023GG0178AbstractEstimating the hydralic properties of subsurfaceaquifersiscrucialtoproperlyevaluategroundwaterresources Apart from laboratory methods and field loggingexperiments，theinversionofaquiferparametersiscurrently mostly performed using the response function ofwater level to external periodic loadings Among them，therearemanymethodsbasedonquasistaticporoelasticity，but these approaches have ignored thepossible frequency-dependent solid/fluid coupling in theporous aquiferIn this paper，a hydrological coseismicresponse model based on dynamic poroelastic theory isfirstlyproposedtoderivethefrequency-dependentrelationship between aquifer pore pressure and stress Byconstructing a dynamic response function，the long-termevolution of well aquifer permeability was inferred from thenear-field coseismic water level response observed in theTengchong well in Yunnan Province The results show thatthe dynamic poroelastic theorycanwellexplainthecoseismic water level response and appropriately invert theaquifer parameters Near-field seismic can cause a largeshort-term increase in the permeability of the well aquiferand recover shortly afterwards，which may be related to thetransient solid deformation，pore opening，and particledistribution due to the near-field static stressKeywordsWater level;Coseismic response;Dynamicporoelastic theory2023，38(1)高楚林，等:腾冲井同震水位动态孔弹性响应的渗透率反演(www progeophys cn)0引言地下水分布广泛、易流动、不可压缩、对地下介质中的应力和固体变形非常敏感(Sun et al，2015)当含水层位于一个紧密的承压系统时，地下水可作为一种高度敏感的应变仪，放大井含水层系统的应力扰动，响应多种外部载荷，例如地球潮汐(Hsieh et al，1988)、大气负载(Van Der Kamp andGale，1983)、断裂活动(Zhang and Huang，2011)、地震波传播(Kitagawa et al，2011)其中，地震孕育过程中的应力积累，断层位错引起的静应力，地震波传播触发的动应力都极易导致孔隙压力和岩体水文性质的变化地震监测中同震井水位响应特征逐步受到重视(崔瑾等，2019;苏鹤军等，2020;向阳等，2020)实践中，地震触发地下水位变化在很多方面有重要意义:例如指示地壳变形(Johnston et al，1995)、油井的供水和生产(oberts et al，2003)、土壤液化(Wang，2007)、余震分布(Bosl and Nur，2002)、地震触发和火山爆发(Prejean，2004)，以及提供水文地质 地 震 前 兆 的 证 据(Matsumoto and oeloffs，2003)不同近震触发地下水位响应机制各异，多数记录显示阶梯状的同震井水位变化(Chia et al，2004;Wang and Chia，2008;胡小静等，2020;陈玮等，2020)基于似稳态孔弹性理论，有学者提出这可能由不排水扩张和沉积物的固结所致，并与响应地震应变的单向传递及饱和水的扩散有关(Zhangand Huang，2011)，应力在“不排水”条件下忽略固流的互耦合从而触发孔隙流体流动(Wang et al，2018)而在实际观测中，“不排水”条件通常难以持续(Zhang et al，2015，2017)，尤其对于像地震波这样的高频响应(Zhang et al，2009)，固流之间应存在相互作用本文因此利用动态孔弹性理论充分考虑同震井水位响应形态及固流耦合的频率依赖，研究了云南省腾冲井由腾冲 MW5 0 级地震触发的井水位响应，观测数据集包括流体观测井水位和大气压力 基于潮汐分析结果，改写了频率依赖的孔压和围压比值(特征 B()系数)，估算了井含水层的渗透率演化过程，并同 1998 年测井实验得到的渗透率值进行比较，分析了此方法在机理解释上解决的问题和同震响应对井含水层渗透率的影响1观测和数据腾冲地震台流体观测井(以下简称“腾冲井”)位于 NE 向大盈江断裂附近，在腾冲县城东南郊腾越镇董库村地震台内，海报高度 1620 m，地理坐标为 25 03N、98 5E 腾冲井井深 120 m，非自流井，属静水位观测，由于受断裂影响而构造发育，裂隙率达 12%，观测水层为第四系英安岩、安山岩及熔结凝灰岩沙层，地下水温度为 21 5，井孔构造图详见图 1a 水位观测仪器为 LN-3 型数字水位仪，测量精度为 0 2%F S，稳定度为 0 25%F S，分辨率为 1 mm 水位测量系统由压力传感器和数据采集器组成，采样率为 1 min 大气压力观测来自同一观测站的 TP-1 雨量气温气压观测仪，采样率同为1 min为了研究地震触发井水位响应，选择了距离腾冲井 30 41 km 发生于 2011 年 6 月 20 日 18 时 16分，震中25N，98 8E，深度15 7 km 的腾冲 MW5 0级 地 震(震 源 参 数 来 自 GCMT:https:/wwwglobalcmt org/CMTsearch html)，以下简称 6 20 地震 该地震震中位于腾冲小断裂附近，和腾冲井的空间位置详见图 1b图 2a 为 6 20 地震触发的腾冲井水位同震响应，表现为阶梯状变化伴随脉冲振荡(纵坐标为静水位测量值，即井口到水平面距离)，且在震后 1 天左右出现大幅度振荡 图 2 对比了水位和气压在6 20 地震前后的长期变化趋势，呈现出相似的周期性振荡 2011 年 2 月至 2011 年 4 月，水位出现季节性下降，对应于腾冲地区的旱季，后水位因为雨季的到来降雨的增多逐渐上升至 9 月进入平稳期 6 20地震触发的水位异常振荡位于雨季的长趋势上升中，因此，水位的同震及震后响应清晰区别于背景场，且同震振荡振幅约为日常潮汐振幅的 4 倍2频率域的含水层参数估算2 1潮汐分析为了确定水位和气压的频谱响应，采取小波变换进行时频分析 图 3a、b 为水位、大气压力降采到每小时后得到的小波功率谱，尽管水位和气压在256 h 及以上周期有一些模糊响应，但在 M2 潮汐频段(周期 T=745 2 min=12 42 h)均呈现出显著异于背景场且持续的响应，在图 3a 的其他频段，水位也呈现出明显的同震响应(白色虚线为发震时刻)图 3c、d 为水位和气压的交叉小波变换及小波相干，31地球物理学进展www progeophys cn2023，38(1)图 1腾冲流体观测井的井孔构造图(a)和腾冲井及 6 20 地震空间分布图(b)图(b)中的红色曲线表示断裂，蓝色曲线表示河流Fig 1Aquifer lithology and wellbore structure of Tengchong(TC)well(a)andlocations of TC well and epicenters of 6 20 earthquake(b)The red curves represent faults and the blue curves represent rivers in fig(b)图 2腾冲台记录的 2011 年 2 月 19 日至 2011 年10 月 16 日水位(a)、大气压力(b)时序变化红色箭头为地震发震时刻(2011-06-20，18:16)，右上角图为井水位响应于 6 20 地震的 24 h 时序Fig 2Water-level(a)and barometric-pressure(b)variations time series in the TC well during fromFebruary 19，2011 to October 16，2011The red arrow shows the time of the earthquake(2011-06-20，18:16)The top right graph shows the 24-hour time series ofwell water levels in response to the 6 20 earthquake可以看到 M2 频段为两个时序持续公共高功率的区域且高度相关，整个数据集在 8 16 h 的波段内(主要围绕半日波段)相干系数均大于 0 9 而因 M2 潮汐频段的高相关性和日间温度的低污染，气压和潮汐导致的水位响应常在此频段被分析(Doan et al，2006;ojstaczer，1988)在一个流体观测井研究系统中，气压通常会同时加载到地球表面和井孔，为垂向单向影响，地下水则会由于压差流进流出 当潮汐体积应变加载到含水层时，水位会响应于其导致的含水层固体形变(压缩和扩张)(ojstaczer，1988)由于井含水层限制，水位与潮汐间会存在相位滞后，且与含水层水文参数(渗透率等)密切相关，因此利用相移特征反演水文参数的研究很多(Zhang et al，2019;向阳等，2017)图 4 描述了简化的观测井筒圆柱模型的响应进程，采用了三维柱坐标系基于小波结果和前人经验，利用 Baytap-G 软件(Tamura et al，1991)选择 30 天移动窗口 4 天步长在 M2 频段获得了水位对气压的幅值比和相移(图5)，该软件使用了贝叶斯模型程序进行时间序列的分析以及 Akaike 贝叶斯信息准则(ABIC)寻找拟合最好的模型图 5 同震期间水位的幅值比和相移变化显著，震后幅值放大因子保持稳定，相移则在地震应力扰动下突变后不久恢复，这和之前一些研究类似(Elkhoury et al，2006;Yan et al，2016)这表明地震波的动态应力扰动改变了井含水层渗透性，造成 M2 频段的一段时间内固流耦合关系改变，通过动态响应关系并利用振幅和相位响应即能推算出渗412023，38(1)高楚林，等:腾冲井同震水位动态孔弹性响应的渗透率反演(www progeophys cn)图 3腾冲台记录的 2011 年 2 月 19 日至 2011 年 10 月 16 日水位(a)、大气压力(b)小波变换、两者的交叉小波变换(c)和小波相干(d)其中白色虚线为 620 腾冲地震的发震时刻Fig 3The wavelet transform of water level(a)，atmospheric pressure(b)，cross wavelet transform of both(c)andwavelet coherence(d)recorded at Tengchong station from February 19，2011 to October 16，2011The white dotted line is the time of 620 Tengchong earthquake图 4响应于大气压力的圆柱井含水层系统简化模型Fig 4A cylinder well-aquifer response tobarometric pressure透率的变化2 2动态孔压响应估算基于似稳态孔弹性理论，井中水位受气压触发的关系可以写成(Sil and Freymueller，2006):B=3pa1 vu1+v()u，(1)其中 a 是气压改变，p 是孔压改变(p=fgh)B 是多孔介质下的不排水 Skempton 系数，定义为触图 5利用 Baytap-G 软件(Tamura et al，1991)、30 天窗口和 4 天步长在 M2 频段得到的 2011 年 2 月 19 日至 2011 年 10 月 16 日水位对大气压力的幅值比和相移红色箭头所指数据点计算窗口为 2011 年 5 月 22 日至 2011 年 6 月 20 日Fig 5Amplitude amplification factor and phase shift ofwater level to barometric-pressure at M2 frequencywith a moving time window of 30 days and 4-daysteps from February 19，2011 to October 16，2011 byBaytap-G software(Tamura et al，1991)The data point calculation window indicated by the redarrow is from May 22，2011 to June 20，201151地球物理学进展www progeophys cn2023，38(1)发孔压 p 和含水层围压 kk/3 的比(式(2)反映固体应变对气压对水位的单向作用(图 4)vu是不排水 泊 松 比，表 示 触 发 过 程 中 流 体 质 量 不 变(oeloffs，1996;Wang，2000)式(2)为:B=pkk/3(2)如前文分析，不排水条件的假设，忽略了流体-固体的相互耦合作用 事实上，水位和气压之间的相移(图 5)并不包含在不排水的假设中，需要建立不包含固相的流体扩散模型解释 但若考虑动态孔弹性模型，可以将固体应变和流体流动统一到一个物理框架中，这时使用动态孔弹性理论修改公式(1)得到动态 B()系数的关系式:B()=h()a()1 vu1+v()u3fg(3)修正的频率域 B()系数，是一个包含幅值和相位的复数 h()和 a()同样指频率域的水位和气压，其对应了包含固流耦合的孔压和围压关系:B(，k)=p()kk()/3，(4)考虑井筒结构，其中 kk=zz+rr(方向的扭转力不会影响井筒体积改变)根据动态孔弹性方程(Biot，1962)，可以推导动态孔弹性方程中频率-波数域触发孔隙压力(p)与、Z 方向应力扰动(zz和 rr)在柱坐标系下的一般解图 5 中正相移表明了腾冲井垂直流动扩散行为，因此选择边界条件 z=0 时(0)=0，从而得到对应的频率依赖的 B()系数:B()=(+1)Me 1z(+1)Me 2z23(+)+(+1)M e 1z+23(+)+(+2)M(+1)(+2)e 2z，(5)其中 a=1 Km/Ks(Biot-Wills 系数)、1/M=()/Ks+/Kf，Km、Ks和 Kf分别为固体框架、固体颗粒和液态的体积模量，和 为固体框架的拉梅系数，i、i和 3为:i=2 k2Pi，i=1，2i=(+2)(2i 2)+2 f2(m2 ib)f2，i=1，23=f ib+m，(6)其中代表瑞雷波的复波数，可用 ks/0 92 近似=(1 )s+f为含有饱和流体的体积密度，为孔隙度，s和 f分别为固体颗粒和液体的密度 阻力-阻尼耦合系数 b=/，为动态流体黏度，为宏观渗透率2 3渗透率反演动态 B()系数受多参数控制，基于腾冲井的测井报告，相对于井含水层渗透率的变化，弹性和密度系数可以认为是稳定的，即 Ks=50 GPs，s=2700kg/m3(谢仁海等，2007)，f=977 882 kg/m3，g=9 8 m/s2，u=0 29 取 Kf=2 18 GPa，0=0 974 103Ns/m2，=0 12 根据 SCA 自相容模型(Mavko et al，1998)，孔隙纵横比 0 06，可得到固体框架体积模量和第二拉梅常数 Km=7 6674 GPa，m=6 1286 GPa 尽管在这里，B()系数与泊松比u也有关，但该参数在同震期间通常假定不变 另外，Smit 等人(2002)发现 B 系数和泊松比 u的变化幅度分别为 6 3%和 0 3%，前者远大于后者，因此 u的变化可以忽略不计选择渗透率 为 B()系数拟合的变量，B()系数幅值和相位的频散特征见图 6，其中 u和 l分别为任取渗透率的上下限图7 为基于公式(3)获得的 M2 频段 B()系数进一步扩展拟合动态 B()系数解析式(公式(5)从而反演得到的腾冲井渗透率 长期演化结果，将腾冲井 1998 年测井实验获得的渗透率值 0=6 9164 1013m2与模型反演的渗透率进行比较，可以发现2011 年2 月19 日至2011 年10 月16 日间反演的大部分渗透率略小于测井实验值，但也基本在 0一个数量级波动，6 20 地震发生后，伴随相移的增大，渗透率也突然增大，这可能与近场地震应力导致的固体变形和裂隙的打开以及颗粒的重分布有关 该渗透率异常约两个月恢复，表明应力的扰动带来的是弹性响应而不是塑性的改变3机理讨论地下水对应力和地壳的固体变形非常敏感 一些研究记录了井中水位对含水层压头的放大效果(Cooper et al，1965;Liu et al，1989)在基于似稳612023，38(1)高楚林，等:腾冲井同震水位动态孔弹性响应的渗透率反演(www progeophys cn)图 6动态围压与孔压比例系数 B()的频散颜色标代表不同的渗透率，其他参数按照理论值定义，(a)和(b)分别为幅值和相位随频率的变化，黑色虚线为 M2 频段位置Fig 6Dispersion of dynamic surrounding pressure and pore pressure scaling factor B()The color scale represents different permeability，other parameters are defined according to theoretical values，(a)and(b)are respectivelythe variation of amplitude and phase with frequency，black dotted line is M2 frequency band图 7基于动态 B()系数拟合的含水层渗透率时间演化Fig 7Time evolution of aquifer permeability based ondynamic B()-coefficient fitting态孔弹性理论的水位潮汐分析中，通常假设固体应变在“不排水”条件下触发流体流动(Zhang andHuang，2011)，水位和应力的相移被解释成纯流体不同方向的扩散，忽略固流耦合的频率依赖(Cooperet al，1965;ojstaczer，1988)而对于不排水条件，当潮汐应变的波长远大于含水层大小时，的确可以假定含水层系统不排水(Zhang and Huang，2011)但是观测同样已经证明该假设存在局限性(Kitagawa et al，2011;Zhang etal，2015)，因为含水层系统既不是完全封闭的也不是完全开放的 在低频的加载下，封闭含水层可能在其边界进行流体运动;高频加载下，一个开放的含水层可能展现一些封闭性质(Wang et al，2018)，从动态 B()系数的频散曲线特征可以看出(图6)，相比于似稳态孔弹性理论，动态孔弹性理论可以更好的解释该现象图 5 中水位对于气压的幅值放大因子和相移也清晰的表明了含水层流体流动的频率依赖特征(Elkhoury et al，2006;Wang and Chia，2008)这一现象同样可以用考虑固流耦合和地震波引起应力扰动诱导多孔流动的动态孔隙弹性模型更好的解释，其中进一步还可以引入多孔通道中的毛细管截留(Beresnev et al，2011)当引入动态孔弹性理论解释动态触发时，驱动流体的孔压不再是似稳态理论中的 Skempton 系数(B 系数)受围压单向传递，而会反过来作用于固体骨架，形成动态耦合过程，机理上可能能更完备地解释包括地震波触发的裂隙打开闭合等动态过程，一定程度上解释水位同震响应形态根据前人对传统垂直扩散模型(Hsieh et al，1988)的研究结果(Shi et al，2015;Wang et al，2018)，图 5 中，震后幅值放大因子的减小和相移的增大可能意味着含水层和隔水层的渗透率都有所增加(Zhang et al，2019)，这和通过动态孔弹性模型反演的结果相匹配(图 7)另外，通过将动态 B()系数反演的渗透率值与测井实验值相比，反演结果处在合理区间，而渗透率在地震发生之后明显的增加和不久后的恢复同样值得关注，这表明地震会导致井含水层渗透率可逆的改变4结论腾冲井提供了水位和大气压力的观测，计算两者 M2 频段的幅值因子和相移，识别到了近场地震71地球物理学进展www progeophys cn2023，38(1)导致的震后幅值因子的减小和正相移的增大，理论上对应于井含水层渗透率的增大，结合 M2 频段幅值因子和相移的结果通过使用动态孔弹性理论，可以反演长时间序列的井含水层渗透率的演化，反演的结果准确识别出了因地震扰动而带来的渗透率增大，与潮汐分析结论匹配 在动态孔弹性理论模型中，可以重写幅值放大因子和相移，从而估算频率依赖的孔压和围压的动态比值 B()系数，这解决了传统似稳态模型中固流耦合忽略相移的问题 另外，若进一步利用同震时间序列，结合动态 B()系数和由地震波导致体应力扰动，通过该模型也可以预测整个频带里的同震水位或直接约束动态 B()系数从而反演井含水层渗透率 最后，本实验主要研究了腾冲井的单个近场同震触发，未来仍然需要更多数据和研究对该模型进行验证和多情况分析致谢感谢审稿专家提出的修改意见和编辑部的大力支持!eferencesBeresnev I，Gaul W，Vigil D 2011 Direct pore-level observation ofpermeability increase in two-phase flow by shaking Geophys esLett，38(20):L20302，doi:10 1029/2011GL048840Biot M A 1962 Mechanics of deformation and acoustic propagation inporous media J Appl Phys，33(4):1482-1498，doi:101063/1 1728759Bosl W J，Nur A 2002 Aftershocks and pore fluid diffusion followingthe 1992 Landers earthquake J Geophys es:Solid Earth，107(12):2366，doi:10 1029/2001JB000155Chen W，Zhang L F，Zhang W H，et al 2020 Analysis on comparisonsof co-seismic response of underground water along the Xiannshanfault zone Journal of Geodesy and Geodynamics(in Chinese)，40(1):104-110Chia Y，Wang Y S，Chiu J J，et al 2004 Changes of GroundwaterLevel due to the 1999 Chi-Chi Earthquake in the Choshui iverAlluvial Fan in Taiwan Bull Seismol Soc Am，91(5):1062-1068，doi:101785/0120000726Cooper H H，Bredehoeft J D，Papadopulos I S，et al 1965Theresponse of well-aquifer systems to seismic wavesJGeophyses，70(16):3915-3926，doi:10 1029/JZ070i016p03915Cui J，Ding F H，Zeng X W，et al 2019 Study on coseismic responsecharacteristics and mechanism of well water level observation inNingxia Progress in Geophysics(in Chinese)，34(4):1281-1287，doi:106038/pg2019CC0209Doan M L，Brodsky E E，Prioul，et al 2006Tidal analysis ofboreholepressure-atutorialCambridgeMA:Schlumbergeresearch eportElkhoury J E，Brodsky E E，Agnew D C 2006 Seismic waves increasepermeability Nature，441(7097):1135-1138，doi:10 1038/nature04798Hsieh P A，Bredehoeft J D，ojstaczer S A 1988 esponse of wellaquifer systems to Earth tides:Problem revisited Water esourcesesearch，24(3):468-472，doi:10 1029/W024i003p00468Hu X J，Fu H，Li T，et al 2020 Study on Co-seismic esponse andPrediction Significance of Groundwater Level in the Dazhai WellJournal of Seismological esearch(in Chinese)，43(2):340-347Johnston M J S，Hill D P，Linde A T，et al 1995 Transient deformationduring triggered seismicity from the 28 June 1992 MW=7 3 Landersearthquake at Long Valley volcanic caldera，CaliforniaBullSeismol Soc Am，85(3):787-795Kitagawa Y，ItabaS，MatsumotoN，etal2011Frequencycharacteristics of the response of water pressure in a closed well tovolumetric strain in the high-frequency domain J Geophys es:Solid Earth，116(B8):B08301，doi:101029/2010JB007794Liu L B，oeloffs E，Zheng X Y 1989 Seismically induced water levelfluctuations in the Wali Well，Beijing，China J Geophys es:SolidEarth，94(B7):9453-9462，doi:101029/JB094iB07p09453Matsumoto N，oeloffs E A 2003 Hydrological response to earthquakesin the Haibara Well，central Japan II:possible mechanism inferredfrom time-varying hydraulic properties Geophys J Int，155(3):899-913Mavko G，Mukerji T，Dvorkin J 1998 The ock Physics HandbookCambridge:Cambridge University PressPrejean S G 2004 emotely triggered seismicity on the United StatesWest Coast following the MW7 9 Denali fault earthquake BullSeismolSocAm，94(6B):S348-S359，doi:10 1785/0120040610oberts P M，Esipov I B，Majer E L 2003 Elastic wave stimulation ofoil reservoirs:Promising EO technology?The Leading Edge，22(5):448-453oeloffs E 1996Poroelastic techniques in the study of earthquake-related hydrologic phenomena Adv Geophys，37:135-195，doi:101016/S0065-26870860270-8ojstaczer S 1988 Intermediate period response of water levels in wellsto crustal strain:Sensitivity and noise level J Geophys es:SolidEarth，93(B11):13619-13634，doi:101029/JB093iB11p13619Shi Z，Wang G，Manga M，et al 2015 Continental-scale water-levelresponse to a large earthquake Geofluids，15(1-2):310-320，doi:101111/gfl12099Sil S，Freymueller J T 2006 Well water level changes in Fairbanks，Alaska，due to the great Sumatra-Andaman earthquakeEarth，Planets and Space，58(2):181-184，doi:101186/BF03353376Smit T H，Huyghe J M，Cowin S C 2002 Estimation of the poroelasticparameters of cortical bone J Biomech，35(6):829-835，doi:101016/S0021-9290(02)00021-0Su H J，Cao L