承压含水层应力场与井水位变化的数值模拟

谢雨晴 赵永红†

ACTA Scientiarum Naturalium Universitatis Pekinensis - - 北京大学学报 -

谢雨晴 赵永红

北京大学地球与空间科学学院, 北京 100871; † 通信作者, E-mail: zhaoyh@pku.edu.cn

摘要 将含水层看成孔隙弹性介质, 建立二维承压含水层理想模型, 计算周围应力场变化引起的井水位变化的分布和量级, 讨论井水位突变的消散规律。根据消散过程, 用反褶积的方法, 从井水位数据计算周围应力应变场的变化, 并用数值方法和会理川06井的固体潮数据验证这种方法的可行性。应用数值模拟结果探讨同震井水位阶变的机制, 认为简单的孔隙弹性模型不能解释同震井水位变化特征。关键词 井水位; 渗流; 应力场变化中图分类号 P553

为质量迁移和能量传递两种。固体应力 渗流问题是由应力 应变场与渗流场的耦合引起的。大地震的发生伴随着区域应力场的调整, 引起深部承压含水层体应变的变化[3], 导致孔隙压力发生改变。孔隙压力的变化传到地表后表现为承压井水位的涨落。地应力的增加产生压缩区和张性区, 压缩区水位会逐渐抬升, 张性区水位下降[4]。当地层发生倾斜时, 地下水在重力场中取得新的平衡, 也会发生流动。由于地下水的流向、流速和流量发生变化, 水中的化学元素组成也会发生变化。Roeloffs[5]将岩石看成孔隙弹性介质, 研究承

[6]压井水位对应力的响应。Rice 等 利用孔隙弹性(poroelasticity)模型, 将弹性力学场与渗流场完全耦合。他们假设岩体为充满饱和水的多孔介质, 当岩石变形时, 产生的体应变会对渗流场产生影响,而渗流场的水头梯度作为体积力直接作用在应力 应变场上, 方程为

其中, u为孔隙介质的位移矢量, σ 为应力张量, p 是渗流场中的孔隙压力, 为流体密度, evol表示体积应变, U 是渗流场中的流速, αb 为 Biot-willis 耦合系数, κ 为渗透率, μ 为动力黏滞系数, D为高程。机制是渗流‒损伤模型。震源体的应力积累达到岩除应力场与渗流场的耦合作用外, 另一类作用

会进入微破裂‒膨胀(扩容)阶段[1],石的屈服强度后,这时震源体产生的非弹性变形使裂隙增生、微破裂出现和串通以及裂隙开闭, 使孔隙介质的渗透性发生变化, 通过物质上的联系, 将深部的应力状态传到地表, 造成井水位的涨落。

1.2 井水位变化与应力场灵敏度研究

‒应变大小。车用太等[7]对井水位对应力 1 mm 井水位对应的含水层应力井的水位对各类应力‒应变响应的量级(表位于内蒙古的丰镇井进行了较详细的研究, 得出该1)。其中潮汐效应由水位潮汐记录与重力固体潮理论值动态对比得出, 地表水体的计算是将一次洪水看成无限 长、300 m 宽、1.5 m 深水体作用在地表, 通过土力学中矩形荷载下地基附加垂直应力公式得到。列车荷载的作用在丰镇井上是幅度为 1~6 mm 的脉冲记录, 每个脉冲对应一列火车通过, 列车通过后水位可恢复原始水平, 也有少数列车没有响应。这些体应变估算值与井水位变化量的观测值之比, 就是相应的灵敏度。结果表明, 井水位对低频大面积作用力的响应灵敏度大, 对高频局部作用力的响应灵敏度小。

水位灵敏度还与井的孔径[4]和井所处位置的水文地质结构有关, 水位观测段的直径每缩小 10 mm,振荡信息的放大率提高 15%; 当含水层很薄, 封闭性很好时, 对应力场的变化会有更灵敏的响应。

Rice 等[6]认为, 在孔隙弹性模型中, 当渗流速度远小于应力变化速度时, 流体的流动可以忽略,相当于土力学中的不排水情况。由式(3)可以得到应力变化造成井水位的瞬时变化[5]:

∆p = −(2GB/3)[(1 + νu)/(1 − 2νu)]∆ε , (5) ∆p 表示含水层中水压的变化值, ∆ε 表示含水层体应变的增量(拉应变为正), νu表示不排水泊松比, G为岩石剪切模量, B 为 Skempton 孔压系数[8]。

[3]张昭栋等 认为, 如果忽略含水层骨架材料本身的体积变化, 根据弹性理论和地下流体动力学理论, 利用井水位的气压效应, 可以得到含水层垂直向应力变化和压力水头与含水层参数之间的定量关系:

g  v H , (6) 1 n  /(   n) V其中, ρ 为含水层内水密度(kg/m3), n为含水层介质孔隙度, αv和 β 分别为含水层骨架垂直方向压缩系数和水的体压缩系数(m2/n)。利用固体潮效应得到的含水层体应力变化和压力水头与含水层参数之间的定量关系[3]为

α为含水层固体骨架的体压缩系数。

利用这些公式可以计算含水层应力变化, 但由于一般含水层参数未知, 因此通常利用其他方法对井水位与应力之间的比例系数进行估计。例如可以用井水位固体潮效应反演含水层的应力变化(m),通常只反演一个系数, 公式[3]如下:

∆Hw= kw∆m , (8)其中, kw 是一个因井而异的常数, 可以通过固体潮理论计算, 即井水位潮汐幅度与当地平均潮汐力主

[3]应力之比。张昭栋等 用这个理论研究日本秋田Ms 7.7 级地震在中国大陆地壳薄层的同震响应, 计算得出的应力变化呈交替现象, 得到的应力变化在87~9600 Pa 范围内。

刘澜波等[9]对 3口井两年的水位观测数据进行分析, 推出体膨胀引起的压力扰动为

P =−Ea∆α,

(9)

Ea为孔隙介质的体压缩模量, ∆α为孔隙介质的体膨胀, P表示孔隙压力。通过观测资料可以求出含水层体积压缩模量、垂直压缩模量和平均孔隙度。

尹京苑等[10]用有限差分法对 1995 年云南孟连Ms 7.3 级地震和 1996 年丽江 Ms 7.0 级地震震前保山井水位的异常进行模拟, 并用井水位资料反演得到应力场的变化。以上方法或给出弹性均匀介质中井水位变化与应力变化关系的解析解, 或直接给出用实际数据拟合出的经验公式。在研究某些问题的时候, 这些方法是有效的, 但当含水层结构和地应力分布较复杂时, 需要进行应力场的数值模拟。应力场的数值模

拟基于物理规律, 可以通过反演等手段确定参数。相对于解析法, 可以研究更复杂的问题; 相对于经验公式, 能明了数据背后的物理过程, 验证物理机制的正确性, 可以利用井水位数据研究应力变化的情况, 进而研究地壳运动规律。

本研究将含水层看成孔隙弹性介质, 建立二维承压含水层理想模型, 用该模型计算周围应力场变化引起的井水位变化的分布和量级, 讨论井水位突变的消散规律和用反褶积计算含水层受力变化的方法, 为更复杂的井水位与地壳应力应变关系的正演和反演研究奠定基础。

2 理想模型计算

为了研究含水层受应力场作用产生的井水位变化, 本文建立二维承压含水层模型, 在所建模型中使用均匀、各向同性孔隙弹性介质和理想含水层几何模型。

2.1 模型、边界条件及参数

本研究使用 comsol multiphysics 中的 poroelasticity模块进行模拟计算。几何模型为一个长 10 km, 深 1 km 的理想承压含水层二维模型, A 和 B为两个测点, a 为一条测线, 共 218 个单元(图 1)。

渗流场中, 上下边界为不透水边界, 左右边界为 3000 m 的恒定水头。弹性力学场中, 左右边界施加的压力为 f0 (使含水层压缩为正, 拉伸为负)。上下边界为滚轴约束, 模拟过程中最初设 f0 为零,耦合场在所设置的水头边界条件和重力作用下达到平衡。以平衡状态为初始状态, 在含水层左右两边进行加载。本研究分别计算不同 f0作用下的结果。模型中使用的参数如表 2 所示, 其中孔隙度、泊松比和 Biot-willis 耦合系数 αb值来自 Berea 砂岩[11]。渗透率 κ 数据来自 COMSOL 说明文件① , 取值与Berea 砂岩相近。

2.2 加载瞬间井水位变化

理论上, 井水位会在加载瞬间产生突变, 突变的大小与加载瞬体应变的变化量成正比。根据[6] Rice 等 的解析公式(式(1)和(5))以及广义胡克定律, 可以得到 ∆H 和 ∆f0理论上的关系为 ∆H (2 GB / 3)[(1  u)/(1 2  )] u (10)。∆ f0  g使用 2.1 节的边界条件、参数和几何模型进行模拟。模拟结果显示, 井水位在加载瞬间产生突变,图 2 为测点 A 加载 f0瞬间井水位变化量, 其中斜线表示用解析公式(10)得到的结果。将模型中所取参数值代入式(10), 得到∆H  4.524 105 。

∆f0 用模拟出的数据拟合得到的表达式为

∆H = 3.851×10−5 ∆f0 (f0>0)。 (11)从模拟得到的曲线可以看出, 井水位瞬时变化与荷载大小呈正比关系, 说明加载瞬间孔隙弹性介质的力学性质是线性的, 即所引起的体应变与加载力呈正比。这是因为水在岩石中的渗流速度通常很慢, 水压力变化来不及扩散, 因此在加载瞬间, 水和岩石骨架组成的复合介质是弹性的。

当 ∆f0=105 Pa 时, 得到点 A处的体应变变化为−4.74×10−6, 而该点的井水位变化 ∆H 为 3.851 m,因此可以得到, 对于当前几何模型和参数, 井水位与体应变的大致关系为

∆H = −8.1×105∆ε 。(12)经过模拟计算, 可以得到水位变化与体应变的关系。对于真实含水层, 可用实测井水位数据对该模型进行调试。例如, 地球不同位置固体潮引起的体应变可以通过解析公式得到, 结合固体潮引起的井水位变化测量值, 可以得到反映真实介质的含水层模型, 用来反演地壳运动引起的体应变的变化。

2.3 井水位异常消散过程

模拟结果表明, 在井水位发生突变之后, 经过一定时间, 会回到初始状态的平衡位置。图 3(a)和(b)分别为 f0=105 Pa 时 t=2×104 s 和 t=2×105 s 的水头分布与流速矢量图, 水头沿 y 轴方向均匀分布,并且水从中间向两端挤出。这是因为在压力的作用下, 中部的岩石体积压缩造成孔隙压力上升, 而左右两端受恒定孔隙压力边界条件的影响保持不变。通过比较图 3(a)和(b)中前后两个时间点的水头分布, 发现水头梯度最初集中在两端, 因此两端的水压首先释放, 之后压力梯度逐渐向中间扩散。

图 4 为点 A处水头随时间的变化曲线, 正值表示压应力, 负值表示拉应力。可以发现, 在加载瞬间, 3种情况下的水头相对于初始水头 1000 m都有一个瞬时的变化, 其中拉应力下水头下降, 压应力下水头上升。这个变化值在 103 s 时间内基本上不变, 并在 107 s左右回到平衡值。

从图 3 看出, 渗流速度方向基本上是沿 x 轴,所以这里只讨论流速的 x 分量。图 5 为 3 种加载力下点 B 处流速随时间的变化, 速度正值表示往 x 正方向流。从图 5 可以看出, 施加压力时, 水先从中部往两侧流, 然后达到平衡值; 施加拉力时, 水先

从两侧往中部流, 然后达到平衡值。并且, 加载后渗流速度从某一值开始上升, 在 104 s 时达到最大,之后再逐渐变成零。

图 6 为 t=2×104 s 时 x 方向渗流速度沿测线 a的变化。从图 6 看出, 速度都是中心对称的, 也就是两侧的速度是反向的。受到压力时, 流体往外侧 移动; 受到拉力时, 两侧流体往中心移动。并且, 两侧速度较大, 中间速度较小, 中点速度为零。从水头、速度等物理量的计算结果可以看出,在加载瞬间, 边界上压力的改变会在瞬间造成含水层体应变的突变, 使井水位产生阶变。水在岩石孔隙和裂隙系统中的流动需要一定的时间, 在加载瞬间孔隙压力来不及扩散, 此时模型相当于土力学中的不排水条件。由于模型左右边界是恒定水头, 在水压力梯度的作用下, 最终水压梯度可以完全消散,此时相当于排水条件。

2.4 考虑水位异常消散的加载过程反褶积计算

从 2.3 节的结果可以推测, 分析水位变化时,在一定情况下需要考虑水位异常的消散过程。晏锐[12]等 提出体应变信号可以看成是水位信号与系统[13]函数的褶积。史浙明等 利用褶积回归方法去除[14]气压效应和降水等噪音。张昭栋等 提出用褶积滤波方法处理井水位对固体潮响应的滞后效应, 并通过鲁 29 井的现场试验, 表明井水位对信息响应存在“记忆”滞后效应, 滞后程度与含水层导水系数

有关。他们随后提出可以用水平层状承压含水层解释滞后效应, 利用响应时间进行泰勒级数展开的方法去除这种影响, 处理固体潮、气压等噪音的滞后

[15]问题 。图 7 为井水位变压试验结果。在试验中人为改变井口气压, 使其产生阶变, 同时记录井孔水位的相应变化, 直到含水层孔隙压力达到平衡并稳定一段时间。图 7 中两条曲线分别为气压和记录到的井水位随时间的变化。

本文采用 2.1 节的几何模型、边界条件及参数模拟水位消散过程, 研究井水位的记忆效应, 并用褶积方法对其效果进行预测, 最后用反褶积方法对井水位进行校正, 从而得到真实的应力变化。当应力场的变化速度与水位异常消散的时间相 当时, 需要考虑水位异常消散过程。从理论、实测数据和上述模拟中可以看出, 由于孔隙弹性介质的本构关系是线性的, 可以采取褶积方法, 从应力变化和恒定应力造成的水位突变及消散过程得到水位变化的整个过程; 也可以用反褶积方法, 从水位变化得到应力的变化过程。假设 H0为 ∆f0=1 Pa 荷载作用下的水位突变和消散曲线, f0(t)为荷载随时间的变化函数, 则水位的函数可以用下式计算: H( t )  f ( )  H (t  )d 。 (13)

0 0得到水位变化曲线后, 可以通过反褶积运算得到加载过程。

图 8(a)为 f0(t)函数, 取 f0(t)=5000 sin(10−7πt)。图 8(b)为 f0=1 Pa 时, 利用 2.1 节中有限元理想模型‒时间曲算出的 H0(t)曲线。图 8(c)中虚线是用有限元理想模型, 在图 8(a)的加载条件下算出的水位线; 实线是利用 f0(t)和算出的 H0(t), 用褶积公式(式(13))计算得到。图 8(d)中实线是用图 8(c)实线的结果, 根据图 8(b)的曲线反褶积运算得到的曲线(即加载方式); 虚线是用图 8(c)中虚线的结果(即有限元模拟出的结果), 用 H0(t)反褶积运算得出。‒时间曲线,

从图 8(c)可以看出, 无论是有限元计算结果, ‒时间曲线有很大不同。第一还是由褶积公式得到的井水位 曲线形状都与图 8(a)中荷载个峰值出现的时间相对于应力变化的峰值提前很多, 最小值出现的时间也不一样。比较图 8(a)和(c)

中 t=107 s 和 t=2×107 s 两个点, 它们边界上的荷载是一样的, 而井水位值却分别是最大值和最小值。这是由于影响水位变化的因素是体应变的变化速度, 而不是体应变的大小。从式(3)和(4)可以看出,弹性力学场和渗流场通过体应变与孔隙水压对时间的导数联系在一起, 同时受水位消散的影响。对式(13)求导可以得到H ( t )  f ( )  H 0( t  )d , (14)  0因此水位的上升和下降取决于加载历史和水位异常消散的速度。图 8(c)中井水位最小值没有出现在应力变化率最低时, 而是滞后, 这是因为渗流场边界条件是含水层两侧边界上的水头恒定不变, 而水位变化的测点在含水层中心, 此边界对水头变化的影响体现为消散曲线 H0(t), 其形状和特征与加载点至观测点的距离和渗流速度有关。渗流场的响应有记忆效应, 因此井水位的变化与含水层上的加载历史有关, 在相位上会滞后, 体现为式(14)中等式右边的项, 其中 ∆f0′(τ)和 H0(t − τ)分别代表加载历史和含水层的井水位异常消散特性。可以认为, 当 t= 2×107 s 时, 受历史加载过程和向两侧边界消散的影响, 渗流场仍然是从中心流向两侧, 孔隙水压仍在减小。当应力变化的时间与井水位异常消散的时间相当时, 有必要考虑井水位异常消散过程, 而不是直接用体应变量与水位变化量的正比关系得到体应变量。

比较图 8(c)和(d)中实线和虚线可以得出, 用褶积方法研究井水位变化是比较符合实际的。先求出典型的压应力和拉应力的水位异常消散曲线, 利用这条曲线做反褶积滤波, 将水位变化曲线转化为真实的应力变化曲线, 对理解地震前、同震和震后的地应力变化规律有重要意义。需要说明的是, 本研究中用到的有限元模型及边界条件是极其简化的,与实际井水位对应力变化的响应过程可能存在差距, 这里有限元模型的作用在于说明褶积方法的可行性, 并讨论相位滞后的可能原因。

2.5 实例研究

为说明反褶积方法处理实际数据的可行性, 本研究对会理川 06 井 2004 年 1 月 1 日至 2 月 1 日的固体潮记录进行处理。该井位于四川会理县中厂镇(102.06°E, 26.31°N), 属于国家级井网。1991 年1月1日起, 使用 SW40-1 型日记水位计观测静水位, 清晰地记录了固体潮和水震波。该井深度为

600.26 m, 观测层位深 251~283 m, 含水层为石英白云石大理岩断层磨碎带, 地下水类型为裂隙承压水, 渗透系数为 0.0135 m/d [16]。研究时段内理论潮汐体应变值与井水位记录值如图9所示, 数据分辨率为每小时一个值。

从图 9 可以看出, 井水位与体应变的信号之间存在一个很小的相位差, 与 2.4 节中有限元模型模拟结果相似。将理论潮汐体应变值作为输入信号,实测井水位记录当成输出信号, 采用最小二乘反褶积法对水位异常的消散过程曲线进行反演[17], 得到图 10 的固体潮井水位消散曲线(或响应曲线)。

从图 10 可以看出, 对于量级为 10−9 体应变的井水位响应在接近初始时刻时达到最大, 约为 5× 10−4 m, 之后逐渐下降, 并在约 10 h 的时候下降到最小值 0, 之后一直稳定在 0 附近, 这些特性与图8(b)中有限元模型模拟出来的水位异常消散规律相似, 它比实测数据得到的消散曲线更为平滑。可能的原因有实测数据中分辨率不足(每小时一个数据)以及理论潮汐体应变值和实测值的误差等。

图 11为消散曲线与理论潮汐体应变值褶积所

得的水位变化曲线与实测值对比。可以看出, 图9显示的固体潮体应变值与实测井水位的相位差得到较好的修正。可以认为, 反褶积计算所得的水位异常消散曲线可以反映造成该相位差的本质, 即含水层在排水边界条件下井水位响应的滞后现象。通过分析理论潮汐体应变值与井水位实测值的相位差和振幅比, 可以研究含水层的特性。廖欣[16]通过比较 2004 年苏门答蜡 Ms 8.9 地震前后初始相位与振幅比, 发现这两个量在地震前后发生阶变,这可以解释为地震后含水层渗透性增大。与直接解释相位差、振幅比等方法相比, 本文所用的反褶积方法从实测数据出发, 把可能涉及的关于含水层的几何、物理特性等概括为一条消散曲线, 既能在物理上解释相位差, 在数据处理时能消除相位差, 同时也能提取出含水层的物理特性。结合有限元正演模型, 可以进一步研究该含水层的响应特性。还可以利用反演出的井水位响应曲线, 通过褶积计算研究其他应力作用下含水层体应变的规律(如利用同震井水位得到含水层的同震体应变等)。这一方法需要更多的数据进行验证。

3 结论与讨论

本文通过讨论井水位变化的物理机制, 将含水层看成孔隙弹性介质, 建立了二维承压含水层模型,计算出应力场变化引起的井水位变化的量级, 并与解析解进行比较, 给出井水位瞬间变化量与应力变化的关系。当应变变化的周期较长时, 研究应力场变化引起的井水位变化需要考虑井水位异常的消散规律。本文给出计算井水位变化与水位异常消散曲线的反褶积方法, 用数值方法和会理川 06 井的固 体潮数据验证了这种方法的可行性。通过结合孔隙弹性有限元模型和实测数据, 研究了井水位消散过模型与数字信号处理研究井水位对应力‒应变场变程及井水位响应的相位滞后问题, 为以后结合物理化的响应以及地震同震和整个地震周期中井水位变化及含水层体应变变化奠定了基础。本研究中的有限元计算模型为二维理想含水层模型, 今后可以结合水文地质资料建立三维几何模型, 更加真实地模拟特定含水层中特定井对应力应变的响应。汶川地震的发震断层龙门山断裂带是一个带左旋走滑分量的逆冲断层。汶川地震的震源深度为15 km, 最大应力降为 53 MPA, 平均应力降为 18 MPA[18]。用孔隙弹性理论和位错模型计算得到在远场地表附近引起的体应变在 10−8 量级[19], 用式(10)估算得到井水位的变化约为 8 cm。因此, 在绝大多数区域, 很难用孔隙弹性理论解释同震井水位阶变

[13]的量级。史浙明等 利用固体潮效应计算孝义井等 5 口井的汶川地震同震体应变, 其中介休井、孝义井和巢湖井体应变变化在 107~105 量级, 位错理论得到的这几个井的体应变与实测体应变相差约两个量级。Brodsky 等[20]和车用太等[21]发现距离震源很远的井水位也可能有很大量级的响应, 这与体应变随距离快速衰减的理论相矛盾。对此, Brodsky等[20]解释为, 地震波经过含水层时引起的流体快速流动移开了一些沉积物构成的障碍, 导致渗透率提‒含水层系

[21]升。车用太等 则认为地下水位震前异常的模式有可能不只是岩土力学过程, 还包括井统的水动力学过程。通过本文的模拟计算, 我们认为上述差别的可能原因是, 这种估算没有考虑地震发生时渗流通道的改变对井水位变化的影响, 也没有考虑介质参数的不均一性, 因此简单的孔隙弹性模型解释不了同震井水位变化的量级。

需要说明的是, 本研究选取的几何模型和边界条件都比较简单, 为了更真实地理解地震与井水的关系, 需要更多的实验数据、实测数据以及更准确的模型对理论进行验证。

致谢 本文研究和撰写过程得到中国科学院大学王世民教授的建议和帮助, 防灾科技学院廖欣老师提供数据支持, 在此一并表示感谢。

参考文献

277‒286 [1] 王铁城, 鄂秀满. 中国地震地下流体监测系统的现状与展望. 中国地震, 1994, 10(3):

171‒179 [2] 车用太, 鱼金子. 地下水动态映震机制的试验与观测研究. 地震研究, 1992, 15(2): [3] 张昭栋, 赵淑平, 董传富, 等. 井水位阶变与含水

222‒229层所受体应力之间的定量关系. 地球物理学报, 1994, 37(增刊 2): [4] 陆明勇, 黄辅琼, 周峥嵘, 等. 地壳形变与地下水

98‒104异常关系研究进展. 大地测量与地球动力学, 2004, 24(3):

177‒209 [5] Roeloffs E A. Hydrologic precursors to earthquakes: a review. Pure Appl Geophys, 1988, 126: [6] Rice J R, Cleary M P. Some basic stress-diffusion solutions for fluid-saturated elastic porous media with

227‒241 compressible constituents. Rev Geophys Space Phys,

井水位对地壳应力‒ 1976, 14:

113‒ [7] 车用太, 刘五洲, 鱼金子, 等.应变响应灵敏度的研究. 地震, 2003, 23(3): 120

143‒147 [8] Skempton A W. The pore pressure coefficients A and B. Geotechnique, 1954, 4:

7‒12 [9] 刘澜波, 郑香媛. 井水固体潮分析结果及其在地震预报中的应用. 地震, 1985(1):

397‒401 [10] 尹京苑, 赵利飞. 保山井水位异常的数值模拟. 西北地震学报, 2000, 22(4): [11] Detournay E, Cheng H D A. Fundamentals of Poroelasticity // Hudson J A. Comprehensive rock

24‒25 engineering: principles, practice & projects. Oxford: Pergamon Press,1993: [12] 晏锐, 高福旺, 陈颢. 由井 303‒反演含水层体应变. 中国地震, 2007, 23(3): 309 [13] 史浙明, 王广才, 刘国春. 基于汶川地震同震地下

215‒223水位变化反演含水层体应变. 地震学报, 2012, 34(2):

23‒29 [14] 张昭栋, 王立忠, 高玉斌, 等. 用褶积滤波处理井水位对固体潮响应的滞后. 地震, 1993(4):

21‒27 [15] 张昭栋, 张华, 吴子泉. 井水位的“记忆”滞后效应.地震, 1998, 18(1): [16] 廖欣. 承压井水位潮汐异常机理研究: 以会理川-06 36‒37井和川-18 井为例[D]. 长沙: 湖南师范大学, 2010: [17] 李学聪, 刘伊克, 常旭, 等. 均衡多道 1 范数匹配

963‒973多次波衰减的方法与应用研究. 地球物理学报, 2010, 4(53): [18] 张勇, 冯万鹏, 许力生, 等. 2008 年汶川大地震的1186‒1194时空破裂过程. 中国科学: 地球科学, 2008, 38(10): [19] Shi Z, Wang G, Manga M, et al. Continental-scale

310‒320 water-level response to a large earthquake. Geofluids, 2015, 15: [20] Brodsky E E, Roeloffs E, Woodcock D, et al. A mechanism for sustained groundwater pressure changes induced by distant earthquakes. Journal of 503‒518 Geophysical Research: Solid Earth, 2003, 108(B8):

1‒7 [21] 车用太, 鱼金子, 杨会年. 试论地下水位震前异常的来源、机制与模式. 地震, 1988(4):

图 1理想承压含水层计算几何模型Fig. 1 Ideal model of a 2D confined aquifer

图 2井水位瞬时变化与应力变化的关系Fig. 2 Relation between transient variation of groundwater and the variation of stress

图 5 3种加载力下点 B处流速随时间的变化Fig. 5 Variation of Darcy’s velocity with time of point B under loading force of different value

图 4 Fig. 4 3种加载力下点 A处水头随时间的变化Variation of hydraulic head with time of point A under loading force of different value

图 6 3种加载力下测线 a上的速度分布Fig. 6 Distribution of Darcy’s velocity along line an under loading force of different value

图 3两个时刻的水头分布与速度场Fig. 3 Distribution of hydraulic head and Darcy’s velocity field of two epochs

(c)中虚线为有限元计算结果, (d)中虚线是将(c)中有限元计算结果经过反褶积计算, 去除消散过程影响后的结果图 8荷载随时间的变化(a)、水位异常消散曲线(b)、井水位变化曲线(c)和反褶积所得荷载变化曲线(d) Fig. 8 Variation of loading force with time (a), dissipation curve (b), variation of groundwater with time (c), and variation of loading force from deconvolution (d)

图 7 井水位变压实验曲线[15] Fig. 7 Variation of groundwater level with time in variable pressure test[15]

图 10 反褶积得到的水位异常消散曲线Fig. 10 Dissipation curve from deconvolution

图 9会理川 06井理论潮汐体应变值与井水位记录值Fig. 9 Theoretical volume strain of earth tide and recorded water level Huili Sichuan-06 well

图 11 消散曲线与理论潮汐体应变值得到的水位变化曲线与实测值的对比Fig. 11 Comparison between water level calculated from dissipation curve and theoretical volume strain of earth tide and recorded water level

Newspapers in Chinese (Simplified)

Newspapers from China

© PressReader. All rights reserved.