线弹性含水层井水位、孔压对引潮位响应的研究及其应用

刘春平 唐彦东 廖 欣 万 飞 石 云

(防灾科技学院，三河 065200)

线弹性含水层井水位、孔压对引潮位响应的研究及其应用

刘春平 唐彦东 廖 欣 万 飞 石 云

(防灾科技学院，三河 065200)

应用线弹性介质力均衡方程，研究不排水条件下饱水岩体在潮汐力作用下的体应变，提出了承压含水层孔压对引潮高的线性响应方程，并给出了该方程响应系数(E)的物理意义。结合Hsieh等(1987)提出的井水位对孔压响应的振幅比(A)和位相差(α1)公式，进一步推导出了井水位-引潮高振幅比M=EA和位相差α=α1+α2公式。M和α是基于实测井水位和理论引潮高数据计算的，在孔压与引潮高位相差(α2)已知的条件下，由M和α，可依次计算含水层导水系数(T)、水位-孔压振幅比(A)，以及孔压-引潮高振幅比(E)。选择川18、川06井作为实例，计算和分析了T、A、E和M的变化。

含水层 孔弹性介质 引潮高 孔压

0 引言

潮汐力、气压、土壤水分和雪荷载等高频荷载的作用周期从日到季，在地质意义上是非常短的。由于这些荷载的重复性和小量级，它们所引起的地球形变非常接近于线弹性，因此，这些荷载所导致的水压变化可以由线性孔弹性理论精确描述(Wang，2000)。Biot(1962)构建了孔压和介质形变耦合孔弹性理论，并广泛用于解释含水层孔压与应力-应变关系(Bodvarsson，1970;Ge et al.，1992;Neuzil，2003)。潮汐力作为一个已知的小量级高频荷载，所引起的井水位潮汐现象揭示了含水层孔压与岩石形变的耦合关系，一直受到岩石力学和地震学界的关注。他们应用Biot孔弹性理论建立了孔压、井水位与应力(或应变)方程，研究了井水位对潮汐力的响应函数(Bredehoeft，1967;van der Kamp et al.，1983; 张昭栋等，1991，1995)，并通过响应(或固体潮)系数研究含水层介质特征，估计水动力学参数(Hsieh et al.，1987;张昭栋等，1999;2002)。Hsieh等提出了水位和含水层孔压振幅比(A)计算公式，从理论上建立了井水位-孔压的响应关系。本文运用孔弹性应力平衡理论推导出了含水层孔压与引潮位关系式，结合Hsieh等(1987)给出的井水位与孔压之间的关系，分析了井水位、含水层孔压和引潮位三者之间的振幅比和位相差。

1 线弹性含水层孔压随引潮位变化方程

高频荷载周期相对于排水或水力响应发生的时间尺度很短，可应用不排水条件的应力-应变和流体压力式。

式(1)、(2)中σij是应力;p是流体压力;λ，μ为多孔介质固体骨架的拉梅参数;N为饱和地质岩体的Biot模量，是与流/固组合的压缩性有关的弹性参数;β为Biot系数，又称为有效应力系数; δij是 Kronecker函数;εij是岩石应变张量;I1= ε11+ ε22+ ε33是岩石体应变。

如果不考虑岩体运动加速度，在xi(i=1，2，3表示空间3个坐标轴)方向上的应力分量总和与引潮力相同方向上的分力是平衡的(方俊，1984)，即

式(3)中σij(j=1，2，3)表示xi方向上的应力;ρ是岩体密度;ψn是n阶引潮位。对于i=1，2，3，仿上式分别可以写出3个方向上的力平衡方程。

将式(1)代入式(3)，并应用应变与位移关系式 εij=(∂ui/∂xj+ ∂uj/∂xi)/2(i≠ j) 和 εii=∂ui/∂xi可得到式(3)的体应变和位移表达式

式(4)中 ∇2= ∂2/∂x21+ ∂2/∂x22+ ∂2/∂x23是拉普拉斯微分算子。对于 i=1，2，3 分别得到 3 个坐标方向上的位移方程。方程(4)是描述不排水条件下饱和岩石应变(位移)与引潮位关系的偏微分方程;在不考虑孔压(p=0)的条件下，岩石应变(位移)与引潮位关系的方程(方俊，1984)为

比较可知，只要在式(4)中用λ'=λ+β2N代替式(5)中的λ，2个方程的形式便是一样的，它们具有相同形式的解。对于二阶引潮位，方程(5)关于体应变的1个特解为I'1=-ρψ2/(λ'+2μ)，齐次方程( ∂ψn/∂xi=0)关于体应变的通解为 I″1=3ρ(7λ'+6μ)ψ2/(λ'+2μ)(19λ'+14μ)(方俊，1984)。由潮汐应力产生的总体应变为

式(6)中Ku=λ+β2N+2μ/3是不排水条件下饱水岩石的体积模量。由于Ku＞λ+β2N-2μ/3，在式(6)中分母近似为20Ku，其相对误差＜1/20=5%。所以，潮汐应力作用下含水层的体应变可近似为 I≈0.1ρψ2/Ku，代入式(2)得到引潮位作用下含水层孔压 p=-0.1ρBψ2。其中B=βN/Ku是含水层岩石的Skempton系数;-0.1ρψ2是潮汐应力平均值。应用孔压与水柱高度换算关系p=ρwgH(ρw和g分别是水的密度和重力加速度)得到

式(7)中φ=-ψ2/g是引潮高(L);ρ'=ρ/ρw，孔弹性条件下ρ'的变化与岩石体应变在同一个量级(10-6～10-8)上，变化很小可视为常数，一般为2.3～3.2(唐大雄等，1980);式(7)表明不排水条件下饱水岩石孔压线性响应引潮高的变化。岩石Skempton系数(B)在孔弹性变化范围内是常数。对于不同的岩石B值在0和1之间，坚硬致密岩石B=0，黏性土B趋于1。实验获得一些岩石的B值一般＞0.5(晏锐等，2008)，裂隙含水层岩石如砂岩、石灰岩等的B值一般在0.5～0.9之间(Neuzil，2003)。根据前面ρ'和B的值域分析，孔压对引潮高的响应系数(E)一般位于0.115～0.288之间。

2 线弹性含水层井水位对引潮位的振幅比和位相差响应

为消除其他高频与低频荷载对井水位的影响，利用潮汐分波的周期性特点，可以分解具有相同周期的潮汐井水位和孔压。Hsieh等(1987)提出了在潮汐作用下，井水位对孔压响应的复振幅比为

振幅比和位相差分别为

式(8)—(10)中h0和H0分别是井水位和含水层孔压的复振幅;A是水位-孔压振幅比，值域为0＜A＜1;α1是水位-孔压位相差(°);G和F的表达式为

式(11)中ω是潮汐分波频率(1/T);Kei和Ker是零阶Kelvin函数;βw=(ωS/T)1/2rw;T和S分别是含水层导水系数和储水系数。rw和rc分别是井滤水管(渗水段)半径和井筒(水位变动段)半径。因此，振幅比A和位相差α1也都是T和S的函数。Hsieh等(1987)根据(9)—(11)式计算了A和α1随无量纲参数变化的标准曲线。该曲线反映了A和α1在S的变化区间(10-3～10-7)内变化很小;但A随T的增加逐步增加并趋于1，α1随T的增加逐步减少并趋于零。

另一方面，根据潮汐的周期性，可由φ =φ0eiωt和H=H0eiωt分别描述引潮高变化和孔压响应。这里φ0和H0分别是引潮高和孔压的复振幅，t是时间。由式(7)所描述的孔压和引潮高关系适合于孔压和引潮高周期性变化过程的任一点，因此亦有

从式(12)可知，孔压-引潮高的振幅比为E，位相差为零。实际上由于含水层的不均质性，以及裂隙产状等因素的影响，孔压-引潮高位相差一般不为零，记为α2(°)。若孔压瞬时响应体膨胀，孔压与体膨胀位相是一致的，α2可由含水层体膨胀测量数据计算。

将式(12)代入式(8)，结合式(9)可得到井水位-引潮高振幅比和位相差为

式(13)、(14)中α是井水位-引潮位的位相差(°)，它等于井水位-孔压(α1)和孔压-引潮位(α2)位相差之和;M=EA是井水位-引潮高振幅比。当导水系数(T)较大，A趋于1时，即不排水条件下，M=E是井水位-引潮高的最大振幅比。式(13)、(14)就是孔弹性含水层中井水位对引潮高振幅和位相的响应方程。

利用Baytap分析软件，逐月计算M2波井水位-引潮高振幅比(M)和位相差(α)。假设含水层与其他水体无水流交换，且M2波井水位-引潮高位相差是由井与含水层的水流交换产生的，即α2=0，α1=α。将α1值以及其他已知参数ω、rw和rc值代入式(10)、(11)，给定合理的S值，由Matlab软件迭代计算导水系数(T)值。将S和T值代入式(11)计算G和F，并将G和F代入式(9)计算井水位-孔压振幅比(A)，由式(13)M=EA计算孔压-引潮高振幅比(E)。

3 计算实例分析

从四川省选择满足无垂向水流交换且位相差α2可忽略的川18和川06两口监测井。川06井表层2～3m为第四系黄色松散浮土层，向下至61m是裂隙比较发育的砂化白云石大理岩，它与表层浮土一起构成潜水含水层。观测含水带位于251.28～283.27m深度，为大理岩破碎带，渗透系数为0.0135m/d。观测层上覆灰色和暗灰色辉绿岩(深61～152.24m)、蛇纹石化石英白云石大理岩(深152.24～251.28m)，厚度大、裂隙少且不联通(无地下水)，是隔水顶板;观测层下伏石英白云大理岩和钙质绢云母，厚度＞300m，没有发现裂隙，是隔水底板。观测层隔水顶、底板厚度大，裂隙不发育，与含水层无水力联系，即无垂向水流。同理，根据川18井成井剖面，观测含水带分别是中厚层白云石大理岩(深364.66～371.35m)和断层破碎带(深443.29～460.54m)，上覆富水性极弱的黑云母化闪长岩，下伏蚀变的白云石大理岩，无地下水活动的痕迹。故川18井观测含水带与上覆和下伏地层之间也无水力联系。

运用川18和川06井逐时观测井水位与理论引潮高数据，由Baytap软件计算得到M2波井水位-引潮高逐月振幅比(M)和位相差(α)(图1，3，4)。图1是2井的井水位-引潮高位相差，负值表示位相滞后;2005年3月苏门答腊8.7级地震、10月的南亚坎大陆7.6级地震和2007年6月云南普洱6.4级地震在2井的α曲线上均有明显反映。如川18井受苏门答腊地震的影响，位相差(α)急剧下降至0.5°后又迅速升高到38°，然后急剧下降并逐步趋于稳定 (图1)。由于川18井α的最小值趋于0，说明岩石膨胀(或孔压)与引潮位位相差接近于0，即α2=0，α1=α。川06和川18井都位于四川省会理县境内，2井所处自然地理条件差异不大，也可认为川06井的α2=0，α1=α。因此，图1中2井的位相差主要是井与含水层水流交换造成的。此外，川18井除个别时段受地震影响位相差较大外，基本都在13°以下，川06井在苏门答腊地震前的位相差是25°～30°，震后位相差稳定在15°左右，且所有位相差都在12°以上 (图1)。

图1 川18和川06井水位-引潮高位相差Fig.1 Phase shift of water level to tide generating height from the Chuan 18 and 06 wells.

从图1分析，每次震后α的变化都可分为3个阶段：第1阶段使α值急剧下降，反映了井周含水层的迅速疏通(T增加);第2阶段是α升高，可能是上游振动产生的大量沉积物导致井周裂隙阻滞(T减小);第3阶段α趋于稳定，井周沉积物在一定的孔压条件下趋于平衡，T也趋于稳定。

将图1中的α1(=α)值以及2井其他参数ω、rw和rc值代入式(10)、(11)，设定S=10-4，由Matlab软件迭代计算2井的含水层导水系数(T)值(图2)。其中川18井T在10-5量级内变化，川06井T在10-6量级内变化，这与2井抽水试验计算的渗透系数分别为0.1457m/d和0.0135m/d，差1个数量级的结果一致。将含水层的S和T值代入到式(11)分别计算G和F，并代入到式(9)计算井水位-孔压振幅比(A)，将A和M值代入式(13)计算孔压-引潮高振幅比(或潮汐响应系数，E)。这些计算结果分别如图3，4所示。据图3，4，结合图2分析，可知井水位-孔压振幅比(A)趋势性地响应导水系数(T)的变化;井水位-引潮高振幅比M(=AE)通过A的变化也趋势性地响应导水系数(T)的变化。理论上，在孔弹性条件下，孔压-引潮高振幅比(E)应为常数，但2井的E值都是变化的，也是趋势性地响应T的变化。这可能与在地震活动的影响下井周含水层裂隙沉积物的增减有关。从本文第2节对式(7)的分析可知，E主要依赖于B的变化。当井周裂隙发生沉积，孔隙度减少，孔隙内的压力增加时，这部分内压是由水承担的，即水承担压力的比例增加，B增加，E也增加;反之，则E减少。如图4所示，2004年1—12月，T值小、A值小、M值小、E值大的现象，表明井周发生裂隙沉积导致E增加。2005年1月T值陡增，表明井周含水层裂隙疏通导致E值减小。因此，E值的变化是由过水断面沉积物对孔压扰动产生的，与T的变化仍有一定的关系。整体上，A和E随T的变化分别呈正和负相关，故M(=EA)值比A值和E值相对稳定。

图2 川18和川06井含水层的逐月导水系数Fig.2 Monthly transmissivity of the aquifer of Chun 18 and 06 well.

图3 川18井孔压-引潮高、井水位-孔压和井水位-引潮高振幅比Fig.3 Amplitude ratio ofwater level to pore pressure(A)，pore pressure and water level to tidal height(E and M)from well Chuan 18.

图4 川06井孔压-引潮高、井水位-孔压和井水位-引潮高振幅比Fig.4 Amplitude ratio of water level to pore pressure(A)，pore pressure and water level to tide generating height(E and M)from the well Chuan 06.

4 结语

本文针对线弹性含水层中孔压与引潮高的线性响应关系，清晰地定义了响应系数(E)的物理意义和值域范围。结合Hsieh等(1987)关于井水位-孔压振幅比(A)和位相差(α1)的公式，进一步建立了井水位-引潮高振幅比M=EA和位相差α=α1+α2关系式。应用实测井水位和引潮位数据，计算了M2波井水位-引潮高振幅比(M)和位相差(α)。在已知体膨胀(孔压)位相的条件下，基于M和α值，提出了估计含水层导水系数(T)、井水位-孔压振幅比(A)和孔压-引潮高振幅比(或孔压潮汐响应系数，E)等参数的方法流程。以川18和川06井为例，计算了逐月M和α值，分析了振幅比A、E和M随T的变化。这些研究还将有助于进一步探讨井水位变化的形成机理，研究由于含水层(水动力和岩石力学)特征变化导致井水位非线性响应应力变化的异常反映。

方俊.1984.固体潮［M］.北京：科学出版社.

FANG Jun.1984.Earth Tides［M］.Science Press，Beijing(in Chinese).

唐大雄，孙愫文.1980.工程岩土学［M］.北京：地质出版社.

TANG Da-xiong，SUN Su-wen.1980.Science of Engineering Rock and Soil［M］.Geological Publishing House，Beijing(in Chinese).

晏锐，陈颙，高福旺，等.2008.从昌平井体应变、水位对地震波的响应特征求算含水层的Skempton常数［J］.地震学报，30(2)：144—151.

YAN Rui，CHEN Yong，GAO Fu-wang，etal.2008.Calculating B value ofaquifer from volume strain and water level response to seismic waves at Changping seismic station ［J］.Acta Seismologica Sinica，30(2)：144—151(in Chinese).

张昭栋，郑金涵，冯初刚.1991.深井水位的固体潮效应［J］.地震学报，13(1)：66—75.

ZHANG Zhao-dong，ZHENG Jin-han，FENG Chu-gang.1991.Effect of earth tide on deep well water level［J］.Acta Seismologica Sinica，13(1)：66—75(in Chinese).

张昭栋，郑金涵，张广城.1995.水井含水层系统的潮汐响应函数［J］.西北地震学报，17(3)：66—71.

ZHANG Zhao-dong，ZHENG Jin-han，ZHANG Guang-cheng.1995.Response functions of well aquifer system to tide［J］.Northwestern Seismological Journal，17(3)：66—71(in Chinese).

张昭栋，王秀芹，董守德.1999.加卸载响应比在体应变固体潮中的应用［J］.地震，19(3)：217—222.

ZHANG Zhao-dong，WANG Xiu-qin，DONG Shou-de.1999.Application of response ratio of load and unload to bulk strain earth tide［J］.Earthquake，19(3)：217—222(in Chinese).

张昭栋，郑金涵，耿杰，等.2002.地下水潮汐现象的物理机制和统一数学方程［J］.地震地质，24(2)：208—214.

ZHANG Zhao-dong，ZHENG Jin-han，GENG Jie，et al.2002.Physicalmechanism and unitarymathematical equation for tidal phenomena of ground water［J］.Seismology and Geology，24(2)：208—214(in Chinese).

Biot M A.1962.Mechanics of deformation and acoustic propagation in porousmedia［J］.JAppl Phys，33(4)：1482—1498.

Bodvarsson G.1970.Confined fluids as strainmeters［J］.JGeophys Res，75(14)：2711—2718.

Bredehoeft JD.1967.Response ofwell-aquifer systems to earth tides［J］.JGeophys Res，72(12)：3075—3087.

Ge S，Garven G.1992.Hydromechanicalmodeling of tectonically driven groundwater flow with application to the Arkoma foreland basin［J］.JGeophy Res，97(B6)：9119—9144.

Hsieh P A，Bredehoeft JD，Farr JM.1987.Determination of aquifer transmissivity from earth tide analysis［J］.Water Resources Research，23(10)：1824—1832.

Neuzil C E.2003.Hydromechanical coupling in geologic processes［J］.Hydrogeology Journal，11：41—83.

van der Kamp G，Gale JE.1983.Theory of earth tide and barometric effects in porous formations with compressible grains［J］.Water Resources Research，19(2)：538—544.

Wang H F.2000 Theory of Linear Poroelasticity［M］.Princeton University Press，Princeton，New Jersey.

WELL WATER LEVEL CHANGEW ITH TIDE GENERATING HEIGHT FOR LINEAR POROELASTIC AQUIFER AND THEIR APPLICATION

LIU Chun-ping TANG Yan-dong LIAO Xin WAN Fei SHIYun
(Institute of Disaster Prevention，China Earthquake Administration，Sanhe 065200，China)

Volume strain for the saturated rock in the undrained condition under tide force is studied in this paper with mechanical balance equation of linear elastic medium.The equation relating to pore pressure in confined aquifer responding to tide generating height(TGH)is proposed，and the coefficient(E)in this equation is defined with a clear physical meaning.In combination with the formula proposed by Hsieh(1987)for the amplitude ratio(A)and phase shift(α1)of well water level response to pore pressure，the formula is derived for the amplitude ratio M(=EA)and phase shift α(= α1+α2)of well water level response to TGH.M andα can be calculated by measured water level and theoretical earth tide data.Assuming the phase shift(α2)of pore pressure to TGH approximates zero，the transmissivity(T)of the aquifer，the amplitude ratio(A)and response coefficient(E)can be in turn determined by M and α.As an example，Chuan 18and 06 well data are used to calculate M and α，and to estimate T，A and E，and the changes of A，E and M with T are analyzed.

aquifer，poreelastic medium，tide generating height，pore pressure

P312.4

A

0253-4967(2011)01-0133-08

10.3969/j.issn.0253-4967.2011.01.013

2010-10-10收稿，2011-03-04改回。

地震行业科研专项(200808055，200808079)资助。

刘春平，男，1962年生，1983年本科毕业于华东地质学院水文地质工程地质专业，1986年硕士毕业于中国建筑科学研究院勘察技术研究所岩土工程专业，1989年于国家地震局地质研究所获得构造物理专业博士学位，教授，研究方向为地下水动力学和地震地下流体，电话：010-61596137，E-mail：lcp@fzxy.edu.cn。