谢雨晴, 赵永红
北京大学地球与空间科学学院, 北京100871
Yuqing XIE, Yonghong ZHAO
中图分类号:
P553
通讯作者:
Corresponding authors:
收稿日期:
2016-04-11
修回日期:
2016-05-13
网络出版日期:
2017-07-13
版权声明:
2017 《北京大学学报(自然科学版)》编辑部 《北京大学学报(自然科学版)》编辑部 所有
基金资助:
展开
摘要 将含水层看成孔隙弹性介质, 建立二维承压含水层理想模型, 计算周围应力场变化引起的井水位变化的分布和量级, 讨论井水位突变的消散规律。根据消散过程, 用反褶积的方法, 从井水位数据计算周围应力应变场的变化, 并用数值方法和会理川06井的固体潮数据验证这种方法的可行性。应用数值模拟结果探讨同震井水位阶变的机制, 认为简单的孔隙弹性模型不能解释同震井水位变化特征。
关键词：
Abstract The aquifer can be treated as poroelastic media in current research. An ideal 2D confined aquifer model was founded to calculate the response of the groundwater level to the variation of stress field. Then the dissipation law of abnormity of groundwater level was discussed. A deconvolution method to invert stress-strain field from groundwater level data when considering the dissipation process was also discussed. Results from forward calculation and earth tide data of Huili Chuan-06 well were analyzed to confirm the deconvolution method. At last, the mechanism of co-seismic step variation in groundwater level was analyzed. It was concluded that a simple poroelastic model cannot explain the amplitude of variation.
Keywords：
井水位反映井水所在含水层对应位置和深度的静水压强, 并受井所在地区气压的影响。这个关系可以用下式表示:
其中, Hw表示水井水头(m), H表示承压含水层内井所在位置的压力水头(m), Pa表示大气压强(Pa), ρ表示流体密度, g表示重力加速度。
由于地下流体具有传递压力的能力, 因此有从地下深处(目前最深的井深度可以超过3000 m)携带岩体的物理状态和介质的连通性等信息的可能性。当地震孕育和发生时, 地下流体会产生异常现象, 并且也参与地震孕育过程[1]。
车用太等[2]把震前地下水位动态异常信息分为4个环节: 1) 地壳应力活动与含水层应力状态的变化; 2) 含水层变形破坏与孔隙压力的变化; 3) 孔隙压变化信息在井‒含水层系统中的传递; 4) 孔隙压力变化信息在井孔地下水动态中的表现。他们认为异常信号的传递至少有 3 个途径, 首先分为固体介质传递和液体介质传递两大类, 液体介质传递又分为质量迁移和能量传递两种。
固体应力‒渗流问题是由应力‒应变场与渗流场的耦合引起的。大地震的发生伴随着区域应力场的调整, 引起深部承压含水层体应变的变化[3], 导致孔隙压力发生改变。孔隙压力的变化传到地表后表现为承压井水位的涨落。地应力的增加产生压缩区和张性区, 压缩区水位会逐渐抬升, 张性区水位下降[4]。当地层发生倾斜时, 地下水在重力场中取得新的平衡, 也会发生流动。由于地下水的流向、流速和流量发生变化, 水中的化学元素组成也会发生变化。
Roeloffs[5]将岩石看成孔隙弹性介质, 研究承 压井水位对应力的响应。Rice 等[6]利用孔隙弹性(poroelasticity)模型, 将弹性力学场与渗流场完全耦合。他们假设岩体为充满饱和水的多孔介质, 当岩石变形时, 产生的体应变会对渗流场产生影响, 而渗流场的水头梯度作为体积力直接作用在应力‒应变场上, 方程为
其中, u为孔隙介质的位移矢量, σ为应力张量, p是渗流场中的孔隙压力, ρ为流体密度, evol表示体积应变, U 是渗流场中的流速, αB 为 Biot-Willis 耦合系数, κ为渗透率, μ为动力黏滞系数, D为高程。
除应力场与渗流场的耦合作用外, 另一类作用机制是渗流‒损伤模型。震源体的应力积累达到岩石的屈服强度后, 会进入微破裂‒膨胀(扩容)阶段[1], 这时震源体产生的非弹性变形使裂隙增生、微破裂出现和串通以及裂隙开闭, 使孔隙介质的渗透性发生变化, 通过物质上的联系, 将深部的应力‒应变状态传到地表, 造成井水位的涨落。
井水位对应力‒应变响应的灵敏度指 1 mm 井水位对应的含水层应力‒应变大小。车用太等[7]对位于内蒙古的丰镇井进行了较详细的研究, 得出该井的水位对各类应力‒应变响应的量级(表 1)。其中潮汐效应由水位潮汐记录与重力固体潮理论值动态对比得出, 地表水体的计算是将一次洪水看成无限长、300 m 宽、1.5 m 深水体作用在地表, 通过土力学中矩形荷载下地基附加垂直应力公式得到。列车荷载的作用在丰镇井上是幅度为 1~6 mm 的脉冲记录, 每个脉冲对应一列火车通过, 列车通过后水位可恢复原始水平, 也有少数列车没有响应。这些体应变估算值与井水位变化量的观测值之比, 就是相应的灵敏度。结果表明, 井水位对低频大面积作用力的响应灵敏度大, 对高频局部作用力的响应灵敏度小。
表1 井水位对各类地壳应力‒应变扰动的响应[
Table 1 Response of water level to different disturbance of stress and strain in the earth’s crust[
水位灵敏度还与井的孔径[4]和井所处位置的水文地质结构有关, 水位观测段的直径每缩小 10 mm, 振荡信息的放大率提高 15%; 当含水层很薄, 封闭性很好时, 对应力场的变化会有更灵敏的响应。
Rice 等[6]认为, 在孔隙弹性模型中, 当渗流速度远小于应力变化速度时, 流体的流动可以忽略, 相当于土力学中的不排水情况。由式(3)可以得到应力变化造成井水位的瞬时变化[5]:
Δp = -(2GB/3)[(1 + νu)/(1 - 2νu)]Δε , (5)
Δp 表示含水层中水压的变化值, Δε 表示含水层体应变的增量(拉应变为正), νu 表示不排水泊松比, G为岩石剪切模量, B 为Skempton孔压系数[8]。
张昭栋等[3]认为, 如果忽略含水层骨架材料本身的体积变化, 根据弹性理论和地下流体动力学理论, 利用井水位的气压效应, 可以得到含水层垂直向应力变化和压力水头与含水层参数之间的定量 关系:
其中, ρ为含水层内水密度(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级地震震前保山井水位的异常进行模拟, 并用井水位资料反演得到应力场的变化。
以上方法或给出弹性均匀介质中井水位变化与应力变化关系的解析解, 或直接给出用实际数据拟合出的经验公式。在研究某些问题的时候, 这些方法是有效的, 但当含水层结构和地应力分布较复杂时, 需要进行应力场的数值模拟。应力场的数值模拟基于物理规律, 可以通过反演等手段确定参数。相对于解析法, 可以研究更复杂的问题; 相对于经验公式, 能明了数据背后的物理过程, 验证物理机制的正确性, 可以利用井水位数据研究应力变化的情况, 进而研究地壳运动规律。
本研究将含水层看成孔隙弹性介质, 建立二维承压含水层理想模型, 用该模型计算周围应力场变化引起的井水位变化的分布和量级, 讨论井水位突变的消散规律和用反褶积计算含水层受力变化的方法, 为更复杂的井水位与地壳应力应变关系的正演和反演研究奠定基础。
为了研究含水层受应力场作用产生的井水位变化, 本文建立二维承压含水层模型, 在所建模型中使用均匀、各向同性孔隙弹性介质和理想含水层几何模型。
本研究使用comsol multiphysics中的poroela-sticity 模块进行模拟计算。几何模型为一个长 10 km, 深 1 km 的理想承压含水层二维模型, A 和 B 为两个测点, a 为一条测线, 共 218 个单元(图 1)。
渗流场中, 上下边界为不透水边界, 左右边界为 3000 m 的恒定水头。弹性力学场中, 左右边界施加的压力为f0 (使含水层压缩为正, 拉伸为负)。上下边界为滚轴约束, 模拟过程中最初设f0为零, 耦合场在所设置的水头边界条件和重力作用下达到平衡。以平衡状态为初始状态, 在含水层左右两边进行加载。本研究分别计算不同f0作用下的结果。模型中使用的参数如表2所示, 其中孔隙度、泊松比和Biot-Willis耦合系数αB值来自Berea 砂岩[11]。渗透率 κ 数据来自 COMSOL 说明文件①（① COMSOL Multiphysics, Co. Ltd., Solved with COMSOL Multiphysics 5.1, 2015）, 取值与Berea砂岩相近。
表2 理想模型中所取参数
Table 2 Parameters used in the ideal model
理论上, 井水位会在加载瞬间产生突变, 突变的大小与加载瞬体应变的变化量成正比。根据Rice 等[6]的解析公式(式(1)和(5))以及广义胡克定律, 可以得到ΔH和Δf0理论上的关系为
使用2.1节的边界条件、参数和几何模型进行模拟。模拟结果显示, 井水位在加载瞬间产生突变, 图2为测点A加载f0瞬间井水位变化量, 其中斜线表示用解析公式(10)得到的结果。
图2 井水位瞬时变化与应力变化的关系 Fig. 2 Relation between transient variation of groundwater and the variation of stress
将模型中所取参数值代入式(10), 得到
用模拟出的数据拟合得到的表达式为
Δ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)
经过模拟计算, 可以得到水位变化与体应变的关系。对于真实含水层, 可用实测井水位数据对该模型进行调试。例如, 地球不同位置固体潮引起的体应变可以通过解析公式得到, 结合固体潮引起的井水位变化测量值, 可以得到反映真实介质的含水层模型, 用来反演地壳运动引起的体应变的变化。
模拟结果表明, 在井水位发生突变之后, 经过一定时间, 会回到初始状态的平衡位置。图 3(a)和(b)分别为f0=105 Pa时 t=2×104 s 和 t=2×105 s 的水头分布与流速矢量图, 水头沿 y 轴方向均匀分布, 并且水从中间向两端挤出。这是因为在压力的作用下, 中部的岩石体积压缩造成孔隙压力上升, 而左右两端受恒定孔隙压力边界条件的影响保持不变。通过比较图 3(a)和(b)中前后两个时间点的水头分布, 发现水头梯度最初集中在两端, 因此两端的水压首先释放, 之后压力梯度逐渐向中间扩散。
图3 两个时刻的水头分布与速度场 Fig. 3 Distribution of hydraulic head and Darcy’s velocity field of two epochs 图 4 为点 A 处水头随时间的变化曲线, 正值表示压应力, 负值表示拉应力。可以发现, 在加载瞬间, 3 种情况下的水头相对于初始水头 1000 m 都有一个瞬时的变化, 其中拉应力下水头下降, 压应力下水头上升。这个变化值在 103 s 时间内基本上不变, 并在107 s左右回到平衡值。
图4 3种加载力下点A处水头随时间的变化 Fig. 4 Variation of hydraulic head with time of point A under loading force of different value
从图 3 看出, 渗流速度方向基本上是沿 x 轴, 所以这里只讨论流速的x分量。图5为3种加载力下点B处流速随时间的变化, 速度正值表示往x正方向流。从图 5 可以看出, 施加压力时, 水先从中部往两侧流, 然后达到平衡值; 施加拉力时, 水先从两侧往中部流, 然后达到平衡值。并且, 加载后渗流速度从某一值开始上升, 在 104 s 时达到最大, 之后再逐渐变成零。
图5 3种加载力下点B处流速随时间的变化 Fig. 5 Variation of Darcy’s velocity with time of point B under loading force of different value
图6为t=2×104 s时x方向渗流速度沿测线a的变化。从图 6 看出, 速度都是中心对称的, 也就是两侧的速度是反向的。受到压力时, 流体往外侧移动; 受到拉力时, 两侧流体往中心移动。并且, 两侧速度较大, 中间速度较小, 中点速度为零。
图6 3种加载力下测线a上的速度分布 Fig. 6 Distribution of Darcy’s velocity along line an under loading force of different value
从水头、速度等物理量的计算结果可以看出, 在加载瞬间, 边界上压力的改变会在瞬间造成含水层体应变的突变, 使井水位产生阶变。水在岩石孔隙和裂隙系统中的流动需要一定的时间, 在加载瞬间孔隙压力来不及扩散, 此时模型相当于土力学中的不排水条件。由于模型左右边界是恒定水头, 在水压力梯度的作用下, 最终水压梯度可以完全消散, 此时相当于排水条件。
从 2.3 节的结果可以推测, 分析水位变化时, 在一定情况下需要考虑水位异常的消散过程。晏锐等[12]提出体应变信号可以看成是水位信号与系统函数的褶积。史浙明等[13]利用褶积回归方法去除气压效应和降水等噪音。张昭栋等[14]提出用褶积滤波方法处理井水位对固体潮响应的滞后效应, 并通过鲁 29 井的现场试验, 表明井水位对信息响应存在“记忆”滞后效应, 滞后程度与含水层导水系数有关。他们随后提出可以用水平层状承压含水层解释滞后效应, 利用响应时间进行泰勒级数展开的方法去除这种影响, 处理固体潮、气压等噪音的滞后问题[15]。图 7 为井水位变压试验结果。在试验中人为改变井口气压, 使其产生阶变, 同时记录井孔水位的相应变化, 直到含水层孔隙压力达到平衡并稳定一段时间。图7中两条曲线分别为气压和记录到的井水位随时间的变化。
图7 井水位变压实验曲线[ Fig. 7 Variation of groundwater level with time in variable pressure test[
本文采用 2.1 节的几何模型、边界条件及参数模拟水位消散过程, 研究井水位的记忆效应, 并用褶积方法对其效果进行预测, 最后用反褶积方法对井水位进行校正, 从而得到真实的应力变化。
当应力场的变化速度与水位异常消散的时间相当时, 需要考虑水位异常消散过程。从理论、实测数据和上述模拟中可以看出, 由于孔隙弹性介质的本构关系是线性的, 可以采取褶积方法, 从应力变化和恒定应力造成的水位突变及消散过程得到水位变化的整个过程; 也可以用反褶积方法, 从水位变化得到应力的变化过程。假设H0为Δf0=1 Pa荷载作用下的水位突变和消散曲线, f0(t)为荷载随时间的变化函数, 则水位的函数可以用下式计算:
得到水位变化曲线后, 可以通过反褶积运算得到加载过程。
图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 荷载随时间的变化(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)
从图 8(c)可以看出, 无论是有限元计算结果, 还是由褶积公式得到的井水位‒时间曲线, 曲线形状都与图 8(a)中荷载‒时间曲线有很大不同。第一个峰值出现的时间相对于应力变化的峰值提前很多, 最小值出现的时间也不一样。比较图 8(a)和(c)中t=107 s和t=2×107 s两个点, 它们边界上的荷载是一样的, 而井水位值却分别是最大值和最小值。这是由于影响水位变化的因素是体应变的变化速度, 而不是体应变的大小。从式(3)和(4)可以看出, 弹性力学场和渗流场通过体应变与孔隙水压对时间的导数联系在一起, 同时受水位消散的影响。对式(13)求导可以得到
因此水位的上升和下降取决于加载历史和水位异常消散的速度。图 8(c)中井水位最小值没有出现在应力变化率最低时, 而是滞后, 这是因为渗流场边界条件是含水层两侧边界上的水头恒定不变, 而水位变化的测点在含水层中心, 此边界对水头变化的影响体现为消散曲线H0(t), 其形状和特征与加载点至观测点的距离和渗流速度有关。渗流场的响应有记忆效应, 因此井水位的变化与含水层上的加载历史有关, 在相位上会滞后, 体现为式(14)中等式右边的项, 其中 Δf0′(τ)和 H0(t - τ)分别代表加载历史和含水层的井水位异常消散特性。可以认为, 当 t= 2×107 s 时, 受历史加载过程和向两侧边界消散的影响, 渗流场仍然是从中心流向两侧, 孔隙水压仍在减小。当应力变化的时间与井水位异常消散的时间相当时, 有必要考虑井水位异常消散过程, 而不是直接用体应变量与水位变化量的正比关系得到体应变量。
比较图 8(c)和(d)中实线和虚线可以得出, 用褶积方法研究井水位变化是比较符合实际的。先求出典型的压应力和拉应力的水位异常消散曲线, 利用这条曲线做反褶积滤波, 将水位变化曲线转化为真实的应力变化曲线, 对理解地震前、同震和震后的地应力变化规律有重要意义。需要说明的是, 本研究中用到的有限元模型及边界条件是极其简化的, 与实际井水位对应力变化的响应过程可能存在差距, 这里有限元模型的作用在于说明褶积方法的可行性, 并讨论相位滞后的可能原因。
为说明反褶积方法处理实际数据的可行性, 本研究对会理川 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 会理川06井理论潮汐体应变值与井水位记录值 Fig. 9 Theoretical volume strain of earth tide and recorded water level Huili Sichuan-06 well
从图 9 可以看出, 井水位与体应变的信号之间存在一个很小的相位差, 与 2.4 节中有限元模型模拟结果相似。将理论潮汐体应变值作为输入信号, 实测井水位记录当成输出信号, 采用最小二乘反褶积法对水位异常的消散过程曲线进行反演[17], 得到图10的固体潮井水位消散曲线(或响应曲线)。
从图 10 可以看出, 对于量级为 10-9体应变的井水位响应在接近初始时刻时达到最大, 约为 5×10-4 m, 之后逐渐下降, 并在约 10 h 的时候下降到最小值 0, 之后一直稳定在 0 附近, 这些特性与图8(b)中有限元模型模拟出来的水位异常消散规律相似, 它比实测数据得到的消散曲线更为平滑。可能的原因有实测数据中分辨率不足(每小时一个数据)以及理论潮汐体应变值和实测值的误差等。
图 11 为消散曲线与理论潮汐体应变值褶积所得的水位变化曲线与实测值对比。可以看出,
图11 消散曲线与理论潮汐体应变值得到的水位变化曲线与实测值的对比 Fig. 11 Comparison between water level calculated from dissipation curve and theoretical volume strain of earth tide and recorded water level
通过分析理论潮汐体应变值与井水位实测值的相位差和振幅比, 可以研究含水层的特性。廖欣[16]通过比较2004年苏门答蜡Ms 8.9地震前后初始相位与振幅比, 发现这两个量在地震前后发生阶变, 这可以解释为地震后含水层渗透性增大。与直接解释相位差、振幅比等方法相比, 本文所用的反褶积方法从实测数据出发, 把可能涉及的关于含水层的几何、物理特性等概括为一条消散曲线, 既能在物理上解释相位差, 在数据处理时能消除相位差, 同时也能提取出含水层的物理特性。结合有限元正演模型, 可以进一步研究该含水层的响应特性。还可以利用反演出的井水位响应曲线, 通过褶积计算研究其他应力作用下含水层体应变的规律(如利用同震井水位得到含水层的同震体应变等)。这一方法需要更多的数据进行验证。
本文通过讨论井水位变化的物理机制, 将含水层看成孔隙弹性介质, 建立了二维承压含水层模型, 计算出应力场变化引起的井水位变化的量级, 并与解析解进行比较, 给出井水位瞬间变化量与应力变化的关系。当应变变化的周期较长时, 研究应力场变化引起的井水位变化需要考虑井水位异常的消散规律。本文给出计算井水位变化与水位异常消散曲线的反褶积方法, 用数值方法和会理川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]则认为地下水位震前异常的模式有可能不只是岩土力学过程, 还包括井‒含水层系统的水动力学过程。通过本文的模拟计算, 我们认为上述差别的可能原因是, 这种估算没有考虑地震发生时渗流通道的改变对井水位变化的影响, 也没有考虑介质参数的不均一性, 因此简单的孔隙弹性模型解释不了同震井水位变化的量级。
需要说明的是, 本研究选取的几何模型和边界条件都比较简单, 为了更真实地理解地震与井水的关系, 需要更多的实验数据、实测数据以及更准确的模型对理论进行验证。
致谢 本文研究和撰写过程得到中国科学院大学王世民教授的建议和帮助, 防灾科技学院廖欣老师提供数据支持, 在此一并表示感谢。
The authors have declared that no competing interests exist.
中国地震地下流体监测系统的现状与展望 地下水动态映震机制的试验与观测研究 井水位阶变与含水层所受体应力之间的定量关系 地壳形变与地下水异常关系研究进展 Hydrologic precursors to earthquakes: a review Some basic stress-diffusion solutions for fluid-saturated elastic porous media with compressible constituents 井水位对地壳应力‒应变响应灵敏度的研究 The pore pressure coefficients and B 井水固体潮分析结果及其在地震预报中的应用 保山井水位异常的数值模拟 由井‒含水层系统的水位动态反演含水层体应变 基于汶川地震同震地下水位变化反演含水层体应变 用褶积滤波处理井水位对固体潮响应的滞后 井水位的“记忆”滞后效应 承压井水位潮汐异常机理研究: 以会理川-06井和川-18井为例[D] 均衡多道1范数匹配多次波衰减的方法与应用研究 年汶川大地震的时空破裂过程 Continental-scale water-level response to a large earthquake A mechanism for sustained groundwater pressure changes induced by distant earthquakes 试论地下水位震前异常的来源、机制与模式
/
1 井水位变化与应力场关系
1.1 物理模型
1.2 井水位变化与应力场灵敏度研究
应力‒应变来源 体应变/mm-1 潮汐效应 6.6×10-9 地表水体荷载效应 9.6×10-7 列车荷载效应 1.0×10-2
错误！未定义书签。
1 COMSOL Multiphysics, Co. Ltd., Solved with COMSOL Multiphysics 5.1, 2015
2 理想模型计算
2.1 模型、边界条件及参数
参数 取值 Biot-Willis耦合系数αB 0.79 渗透率κ/m2 2.894×10-11 动力黏滞系数μ/Pa · s-1 10-3 泊松比 0.2 水密度/kg · m-3 103 岩石密度/kg · m-3 2750 杨氏模量/Pa 1.44×1010 孔隙度 0.19
2.2 加载瞬间井水位变化
2.3 井水位异常消散过程
2.4 考虑水位异常消散的加载过程反褶积计算
2.5 实例研究
3 结论与讨论
参考文献
[1]
[2]
[3]
[4]
[5]
[6]
[7]
[8]
[9]
[10]
[11]
[12]
[13]
[14]
[15]
[16]
[17]
[18]
[19]
[20]
[21]
〈
〉
