線彈性含水層井水位、孔壓對引潮位響應(yīng)的研究及其應(yīng)用

劉春平 唐彥東 廖 欣 萬 飛 石 云

(防災(zāi)科技學(xué)院，三河 065200)

線彈性含水層井水位、孔壓對引潮位響應(yīng)的研究及其應(yīng)用

劉春平 唐彥東 廖 欣 萬 飛 石 云

(防災(zāi)科技學(xué)院，三河 065200)

應(yīng)用線彈性介質(zhì)力均衡方程，研究不排水條件下飽水巖體在潮汐力作用下的體應(yīng)變，提出了承壓含水層孔壓對引潮高的線性響應(yīng)方程，并給出了該方程響應(yīng)系數(shù)(E)的物理意義。結(jié)合Hsieh等(1987)提出的井水位對孔壓響應(yīng)的振幅比(A)和位相差(α1)公式，進(jìn)一步推導(dǎo)出了井水位-引潮高振幅比M=EA和位相差α=α1+α2公式。M和α是基于實測井水位和理論引潮高數(shù)據(jù)計算的，在孔壓與引潮高位相差(α2)已知的條件下，由M和α，可依次計算含水層導(dǎo)水系數(shù)(T)、水位-孔壓振幅比(A)，以及孔壓-引潮高振幅比(E)。選擇川18、川06井作為實例，計算和分析了T、A、E和M的變化。

含水層 孔彈性介質(zhì) 引潮高 孔壓

0 引言

潮汐力、氣壓、土壤水分和雪荷載等高頻荷載的作用周期從日到季，在地質(zhì)意義上是非常短的。由于這些荷載的重復(fù)性和小量級，它們所引起的地球形變非常接近于線彈性，因此，這些荷載所導(dǎo)致的水壓變化可以由線性孔彈性理論精確描述(Wang，2000)。Biot(1962)構(gòu)建了孔壓和介質(zhì)形變耦合孔彈性理論，并廣泛用于解釋含水層孔壓與應(yīng)力-應(yīng)變關(guān)系(Bodvarsson，1970;Ge et al.，1992;Neuzil，2003)。潮汐力作為一個已知的小量級高頻荷載，所引起的井水位潮汐現(xiàn)象揭示了含水層孔壓與巖石形變的耦合關(guān)系，一直受到巖石力學(xué)和地震學(xué)界的關(guān)注。他們應(yīng)用Biot孔彈性理論建立了孔壓、井水位與應(yīng)力(或應(yīng)變)方程，研究了井水位對潮汐力的響應(yīng)函數(shù)(Bredehoeft，1967;van der Kamp et al.，1983; 張昭棟等，1991，1995)，并通過響應(yīng)(或固體潮)系數(shù)研究含水層介質(zhì)特征，估計水動力學(xué)參數(shù)(Hsieh et al.，1987;張昭棟等，1999;2002)。Hsieh等提出了水位和含水層孔壓振幅比(A)計算公式，從理論上建立了井水位-孔壓的響應(yīng)關(guān)系。本文運(yùn)用孔彈性應(yīng)力平衡理論推導(dǎo)出了含水層孔壓與引潮位關(guān)系式，結(jié)合Hsieh等(1987)給出的井水位與孔壓之間的關(guān)系，分析了井水位、含水層孔壓和引潮位三者之間的振幅比和位相差。

1 線彈性含水層孔壓隨引潮位變化方程

高頻荷載周期相對于排水或水力響應(yīng)發(fā)生的時間尺度很短，可應(yīng)用不排水條件的應(yīng)力-應(yīng)變和流體壓力式。

式(1)、(2)中σij是應(yīng)力;p是流體壓力;λ，μ為多孔介質(zhì)固體骨架的拉梅參數(shù);N為飽和地質(zhì)巖體的Biot模量，是與流/固組合的壓縮性有關(guān)的彈性參數(shù);β為Biot系數(shù)，又稱為有效應(yīng)力系數(shù); δij是 Kronecker函數(shù);εij是巖石應(yīng)變張量;I1= ε11+ ε22+ ε33是巖石體應(yīng)變。

如果不考慮巖體運(yùn)動加速度，在xi(i=1，2，3表示空間3個坐標(biāo)軸)方向上的應(yīng)力分量總和與引潮力相同方向上的分力是平衡的(方俊，1984)，即

式(3)中σij(j=1，2，3)表示xi方向上的應(yīng)力;ρ是巖體密度;ψn是n階引潮位。對于i=1，2，3，仿上式分別可以寫出3個方向上的力平衡方程。

將式(1)代入式(3)，并應(yīng)用應(yīng)變與位移關(guān)系式 εij=(?ui/?xj+ ?uj/?xi)/2(i≠ j) 和 εii=?ui/?xi可得到式(3)的體應(yīng)變和位移表達(dá)式

式(4)中 ?2= ?2/?x21+ ?2/?x22+ ?2/?x23是拉普拉斯微分算子。對于 i=1，2，3 分別得到 3 個坐標(biāo)方向上的位移方程。方程(4)是描述不排水條件下飽和巖石應(yīng)變(位移)與引潮位關(guān)系的偏微分方程;在不考慮孔壓(p=0)的條件下，巖石應(yīng)變(位移)與引潮位關(guān)系的方程(方俊，1984)為

比較可知，只要在式(4)中用λ'=λ+β2N代替式(5)中的λ，2個方程的形式便是一樣的，它們具有相同形式的解。對于二階引潮位，方程(5)關(guān)于體應(yīng)變的1個特解為I'1=-ρψ2/(λ'+2μ)，齊次方程( ?ψn/?xi=0)關(guān)于體應(yīng)變的通解為 I″1=3ρ(7λ'+6μ)ψ2/(λ'+2μ)(19λ'+14μ)(方俊，1984)。由潮汐應(yīng)力產(chǎn)生的總體應(yīng)變?yōu)?/p>

式(6)中Ku=λ+β2N+2μ/3是不排水條件下飽水巖石的體積模量。由于Ku＞λ+β2N-2μ/3，在式(6)中分母近似為20Ku，其相對誤差＜1/20=5%。所以，潮汐應(yīng)力作用下含水層的體應(yīng)變可近似為 I≈0.1ρψ2/Ku，代入式(2)得到引潮位作用下含水層孔壓 p=-0.1ρBψ2。其中B=βN/Ku是含水層巖石的Skempton系數(shù);-0.1ρψ2是潮汐應(yīng)力平均值。應(yīng)用孔壓與水柱高度換算關(guān)系p=ρwgH(ρw和g分別是水的密度和重力加速度)得到

式(7)中φ=-ψ2/g是引潮高(L);ρ'=ρ/ρw，孔彈性條件下ρ'的變化與巖石體應(yīng)變在同一個量級(10-6～10-8)上，變化很小可視為常數(shù)，一般為2.3～3.2(唐大雄等，1980);式(7)表明不排水條件下飽水巖石孔壓線性響應(yīng)引潮高的變化。巖石Skempton系數(shù)(B)在孔彈性變化范圍內(nèi)是常數(shù)。對于不同的巖石B值在0和1之間，堅硬致密巖石B=0，黏性土B趨于1。實驗獲得一些巖石的B值一般＞0.5(晏銳等，2008)，裂隙含水層巖石如砂巖、石灰?guī)r等的B值一般在0.5～0.9之間(Neuzil，2003)。根據(jù)前面ρ'和B的值域分析，孔壓對引潮高的響應(yīng)系數(shù)(E)一般位于0.115～0.288之間。

2 線彈性含水層井水位對引潮位的振幅比和位相差響應(yīng)

為消除其他高頻與低頻荷載對井水位的影響，利用潮汐分波的周期性特點(diǎn)，可以分解具有相同周期的潮汐井水位和孔壓。Hsieh等(1987)提出了在潮汐作用下，井水位對孔壓響應(yīng)的復(fù)振幅比為

振幅比和位相差分別為

式(8)—(10)中h0和H0分別是井水位和含水層孔壓的復(fù)振幅;A是水位-孔壓振幅比，值域為0＜A＜1;α1是水位-孔壓位相差(°);G和F的表達(dá)式為

式(11)中ω是潮汐分波頻率(1/T);Kei和Ker是零階Kelvin函數(shù);βw=(ωS/T)1/2rw;T和S分別是含水層導(dǎo)水系數(shù)和儲水系數(shù)。rw和rc分別是井濾水管(滲水段)半徑和井筒(水位變動段)半徑。因此，振幅比A和位相差α1也都是T和S的函數(shù)。Hsieh等(1987)根據(jù)(9)—(11)式計算了A和α1隨無量綱參數(shù)變化的標(biāo)準(zhǔn)曲線。該曲線反映了A和α1在S的變化區(qū)間(10-3～10-7)內(nèi)變化很小;但A隨T的增加逐步增加并趨于1，α1隨T的增加逐步減少并趨于零。

另一方面，根據(jù)潮汐的周期性，可由φ =φ0eiωt和H=H0eiωt分別描述引潮高變化和孔壓響應(yīng)。這里φ0和H0分別是引潮高和孔壓的復(fù)振幅，t是時間。由式(7)所描述的孔壓和引潮高關(guān)系適合于孔壓和引潮高周期性變化過程的任一點(diǎn)，因此亦有

從式(12)可知，孔壓-引潮高的振幅比為E，位相差為零。實際上由于含水層的不均質(zhì)性，以及裂隙產(chǎn)狀等因素的影響，孔壓-引潮高位相差一般不為零，記為α2(°)。若孔壓瞬時響應(yīng)體膨脹，孔壓與體膨脹位相是一致的，α2可由含水層體膨脹測量數(shù)據(jù)計算。

將式(12)代入式(8)，結(jié)合式(9)可得到井水位-引潮高振幅比和位相差為

式(13)、(14)中α是井水位-引潮位的位相差(°)，它等于井水位-孔壓(α1)和孔壓-引潮位(α2)位相差之和;M=EA是井水位-引潮高振幅比。當(dāng)導(dǎo)水系數(shù)(T)較大，A趨于1時，即不排水條件下，M=E是井水位-引潮高的最大振幅比。式(13)、(14)就是孔彈性含水層中井水位對引潮高振幅和位相的響應(yīng)方程。

利用Baytap分析軟件，逐月計算M2波井水位-引潮高振幅比(M)和位相差(α)。假設(shè)含水層與其他水體無水流交換，且M2波井水位-引潮高位相差是由井與含水層的水流交換產(chǎn)生的，即α2=0，α1=α。將α1值以及其他已知參數(shù)ω、rw和rc值代入式(10)、(11)，給定合理的S值，由Matlab軟件迭代計算導(dǎo)水系數(shù)(T)值。將S和T值代入式(11)計算G和F，并將G和F代入式(9)計算井水位-孔壓振幅比(A)，由式(13)M=EA計算孔壓-引潮高振幅比(E)。

3 計算實例分析

從四川省選擇滿足無垂向水流交換且位相差α2可忽略的川18和川06兩口監(jiān)測井。川06井表層2～3m為第四系黃色松散浮土層，向下至61m是裂隙比較發(fā)育的砂化白云石大理巖，它與表層浮土一起構(gòu)成潛水含水層。觀測含水帶位于251.28～283.27m深度，為大理巖破碎帶，滲透系數(shù)為0.0135m/d。觀測層上覆灰色和暗灰色輝綠巖(深61～152.24m)、蛇紋石化石英白云石大理巖(深152.24～251.28m)，厚度大、裂隙少且不聯(lián)通(無地下水)，是隔水頂板;觀測層下伏石英白云大理巖和鈣質(zhì)絹云母，厚度＞300m，沒有發(fā)現(xiàn)裂隙，是隔水底板。觀測層隔水頂、底板厚度大，裂隙不發(fā)育，與含水層無水力聯(lián)系，即無垂向水流。同理，根據(jù)川18井成井剖面，觀測含水帶分別是中厚層白云石大理巖(深364.66～371.35m)和斷層破碎帶(深443.29～460.54m)，上覆富水性極弱的黑云母化閃長巖，下伏蝕變的白云石大理巖，無地下水活動的痕跡。故川18井觀測含水帶與上覆和下伏地層之間也無水力聯(lián)系。

運(yùn)用川18和川06井逐時觀測井水位與理論引潮高數(shù)據(jù)，由Baytap軟件計算得到M2波井水位-引潮高逐月振幅比(M)和位相差(α)(圖1，3，4)。圖1是2井的井水位-引潮高位相差，負(fù)值表示位相滯后;2005年3月蘇門答臘8.7級地震、10月的南亞坎大陸7.6級地震和2007年6月云南普洱6.4級地震在2井的α曲線上均有明顯反映。如川18井受蘇門答臘地震的影響，位相差(α)急劇下降至0.5°后又迅速升高到38°，然后急劇下降并逐步趨于穩(wěn)定 (圖1)。由于川18井α的最小值趨于0，說明巖石膨脹(或孔壓)與引潮位位相差接近于0，即α2=0，α1=α。川06和川18井都位于四川省會理縣境內(nèi)，2井所處自然地理條件差異不大，也可認(rèn)為川06井的α2=0，α1=α。因此，圖1中2井的位相差主要是井與含水層水流交換造成的。此外，川18井除個別時段受地震影響位相差較大外，基本都在13°以下，川06井在蘇門答臘地震前的位相差是25°～30°，震后位相差穩(wěn)定在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階段是α升高，可能是上游振動產(chǎn)生的大量沉積物導(dǎo)致井周裂隙阻滯(T減小);第3階段α趨于穩(wěn)定，井周沉積物在一定的孔壓條件下趨于平衡，T也趨于穩(wěn)定。

將圖1中的α1(=α)值以及2井其他參數(shù)ω、rw和rc值代入式(10)、(11)，設(shè)定S=10-4，由Matlab軟件迭代計算2井的含水層導(dǎo)水系數(shù)(T)值(圖2)。其中川18井T在10-5量級內(nèi)變化，川06井T在10-6量級內(nèi)變化，這與2井抽水試驗計算的滲透系數(shù)分別為0.1457m/d和0.0135m/d，差1個數(shù)量級的結(jié)果一致。將含水層的S和T值代入到式(11)分別計算G和F，并代入到式(9)計算井水位-孔壓振幅比(A)，將A和M值代入式(13)計算孔壓-引潮高振幅比(或潮汐響應(yīng)系數(shù)，E)。這些計算結(jié)果分別如圖3，4所示。據(jù)圖3，4，結(jié)合圖2分析，可知井水位-孔壓振幅比(A)趨勢性地響應(yīng)導(dǎo)水系數(shù)(T)的變化;井水位-引潮高振幅比M(=AE)通過A的變化也趨勢性地響應(yīng)導(dǎo)水系數(shù)(T)的變化。理論上，在孔彈性條件下，孔壓-引潮高振幅比(E)應(yīng)為常數(shù)，但2井的E值都是變化的，也是趨勢性地響應(yīng)T的變化。這可能與在地震活動的影響下井周含水層裂隙沉積物的增減有關(guān)。從本文第2節(jié)對式(7)的分析可知，E主要依賴于B的變化。當(dāng)井周裂隙發(fā)生沉積，孔隙度減少，孔隙內(nèi)的壓力增加時，這部分內(nèi)壓是由水承擔(dān)的，即水承擔(dān)壓力的比例增加，B增加，E也增加;反之，則E減少。如圖4所示，2004年1—12月，T值小、A值小、M值小、E值大的現(xiàn)象，表明井周發(fā)生裂隙沉積導(dǎo)致E增加。2005年1月T值陡增，表明井周含水層裂隙疏通導(dǎo)致E值減小。因此，E值的變化是由過水?dāng)嗝娉练e物對孔壓擾動產(chǎn)生的，與T的變化仍有一定的關(guān)系。整體上，A和E隨T的變化分別呈正和負(fù)相關(guān)，故M(=EA)值比A值和E值相對穩(wěn)定。

圖2 川18和川06井含水層的逐月導(dǎo)水系數(shù)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 結(jié)語

本文針對線彈性含水層中孔壓與引潮高的線性響應(yīng)關(guān)系，清晰地定義了響應(yīng)系數(shù)(E)的物理意義和值域范圍。結(jié)合Hsieh等(1987)關(guān)于井水位-孔壓振幅比(A)和位相差(α1)的公式，進(jìn)一步建立了井水位-引潮高振幅比M=EA和位相差α=α1+α2關(guān)系式。應(yīng)用實測井水位和引潮位數(shù)據(jù)，計算了M2波井水位-引潮高振幅比(M)和位相差(α)。在已知體膨脹(孔壓)位相的條件下，基于M和α值，提出了估計含水層導(dǎo)水系數(shù)(T)、井水位-孔壓振幅比(A)和孔壓-引潮高振幅比(或孔壓潮汐響應(yīng)系數(shù)，E)等參數(shù)的方法流程。以川18和川06井為例，計算了逐月M和α值，分析了振幅比A、E和M隨T的變化。這些研究還將有助于進(jìn)一步探討井水位變化的形成機(jī)理，研究由于含水層(水動力和巖石力學(xué))特征變化導(dǎo)致井水位非線性響應(yīng)應(yīng)力變化的異常反映。

方俊.1984.固體潮［M］.北京：科學(xué)出版社.

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

唐大雄，孫愫文.1980.工程巖土學(xué)［M］.北京：地質(zhì)出版社.

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

晏銳，陳颙，高福旺，等.2008.從昌平井體應(yīng)變、水位對地震波的響應(yīng)特征求算含水層的Skempton常數(shù)［J］.地震學(xué)報，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.深井水位的固體潮效應(yīng)［J］.地震學(xué)報，13(1)：66—75.

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

張昭棟，鄭金涵，張廣城.1995.水井含水層系統(tǒng)的潮汐響應(yīng)函數(shù)［J］.西北地震學(xué)報，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.加卸載響應(yīng)比在體應(yīng)變固體潮中的應(yīng)用［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.地下水潮汐現(xiàn)象的物理機(jī)制和統(tǒng)一數(shù)學(xué)方程［J］.地震地質(zhì)，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，F(xiàn)arr 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改回。

地震行業(yè)科研專項(200808055，200808079)資助。

劉春平，男，1962年生，1983年本科畢業(yè)于華東地質(zhì)學(xué)院水文地質(zhì)工程地質(zhì)專業(yè)，1986年碩士畢業(yè)于中國建筑科學(xué)研究院勘察技術(shù)研究所巖土工程專業(yè)，1989年于國家地震局地質(zhì)研究所獲得構(gòu)造物理專業(yè)博士學(xué)位，教授，研究方向為地下水動力學(xué)和地震地下流體，電話：010-61596137，E-mail：lcp@fzxy.edu.cn。