董 博，紀(jì)春玲，章 陽，周安聘

（河北省地震局石家莊中心臺(tái)，石家莊，050021）

0 引言

對(duì)承壓含水層井孔而言，井水位對(duì)體應(yīng)變潮汐的響應(yīng)是井水位觀測(cè)中普遍存在的一種現(xiàn)象，井水位潮汐中隱含有豐富的井-含水層系統(tǒng)信息[1]。深入分析井水位潮汐對(duì)推算含水層的各種參數(shù)，衡量井水位對(duì)應(yīng)力應(yīng)變響應(yīng)量、了解含水層動(dòng)態(tài)、解釋井水位微動(dòng)態(tài)等都具有重要的科學(xué)意義和現(xiàn)實(shí)必要性。

Bredehoeft以體應(yīng)變固體潮理論為基礎(chǔ)，結(jié)合水文地質(zhì)學(xué)，從理論上分析了井-含水層系統(tǒng)對(duì)體應(yīng)變潮汐的響應(yīng)，并提出一種利用井水位潮汐信息來確定含水層的孔隙度、骨架體變模量和儲(chǔ)水系數(shù)的方法[2]。之后，不同研究者對(duì)利用潮汐信息確定含水層參數(shù)的方法進(jìn)行了補(bǔ)充[3-5]。其中，Hsieh等推導(dǎo)了井水位對(duì)含水層孔隙壓力響應(yīng)的振幅比和相位差數(shù)學(xué)表達(dá)式，并給出了利用井水位潮汐相位和振幅比來確定井水位滲透系數(shù)和儲(chǔ)水系數(shù)的方法[4]；Elkhoury等利用Hsieh等提供的方法分析了加利福利亞地區(qū)兩口井水位觀測(cè)資料的潮汐振幅比和相位差，認(rèn)為大震發(fā)生后的潮汐振幅比和相位變化是由于地震波增加含水層滲透系數(shù)所致[6]；Bower等運(yùn)用孔彈性應(yīng)力平衡理論推導(dǎo)出了含水層孔壓與引潮位關(guān)系式，分析了井水位受裂隙產(chǎn)狀及徑向排水的變化特征[7]；劉春平等在此基礎(chǔ)上，進(jìn)一步對(duì)受垂向排水的井水位特征進(jìn)行了研究，總結(jié)出井水位變化主要受裂隙產(chǎn)狀、排水影響，并對(duì)變化特征進(jìn)行歸納[8]。

此外，井水位潮汐是評(píng)價(jià)水位觀測(cè)質(zhì)量好壞的一個(gè)重要參數(shù)，對(duì)井水位的潮汐因子、振幅和相位等潮汐參數(shù)的研究可能獲得震前異常信息[9-11]。前人對(duì)井水位潮汐因子、振幅和相位等潮汐參數(shù)動(dòng)態(tài)變化的可能原因、影響因素討論較少，且對(duì)地震前后井水流運(yùn)動(dòng)方向的研究略顯不足。

本文以固體潮理論為研究基礎(chǔ)，在確定了潮汐為該井水位變化主要影響因素的基礎(chǔ)上，利用Baytap-G提供的潮汐分析方法[12],計(jì)算了無極冀20井水位2009—2019年M2和O1波的潮汐因子、相位及振幅值，深入分析了幾次遠(yuǎn)場(chǎng)大震前后該井潮汐因子、振幅以及相位等潮汐參數(shù)的動(dòng)態(tài)變化特征，并試探性地為潮汐參數(shù)的變化找到理論依據(jù)支持，對(duì)地震波能改變含水層儲(chǔ)水系數(shù)、滲透系數(shù)等水文地質(zhì)參數(shù)并使孔隙壓力得到重新分布進(jìn)行了驗(yàn)證，最后得出地震前后水流運(yùn)動(dòng)方向，為解釋遠(yuǎn)距離地震觸發(fā)和水文地質(zhì)響應(yīng)現(xiàn)象提供了實(shí)例。

1 觀測(cè)井孔基本情況

無極冀20井觀測(cè)站地處冀中拗陷，無極低凸起高點(diǎn)，衡水?dāng)嗔驯眰?cè)（圖1）。該井位于河北省無極縣境內(nèi)，地理坐標(biāo)為（38.2°N，114.9°E）。無極井于1989年開始觀測(cè)，孔口標(biāo)高44.5 m，井深2 984 m，屬華北深井網(wǎng)。地下水類型為巖溶裂隙承壓水，含水層為震旦系、寒武系灰?guī)r、白云巖，觀測(cè)層上覆泥巖、砂巖，下伏蝕變的白云巖，可能與含水段發(fā)生水力聯(lián)系。井-含水層水流運(yùn)動(dòng)示意圖見圖2。無極井1986年開始進(jìn)行靜水位模擬觀測(cè)，儀器為SW-40，2001年開始進(jìn)行數(shù)字化觀測(cè)，現(xiàn)使用的數(shù)字化觀測(cè)設(shè)備為ZKGD3000NL型地下水?dāng)?shù)據(jù)監(jiān)測(cè)系統(tǒng)。記錄數(shù)據(jù)質(zhì)量良好、潮汐清晰。

圖1 井孔地質(zhì)構(gòu)造圖

圖2 井-含水層水流運(yùn)動(dòng)模型

2 井-含水層水流徑向、垂向運(yùn)動(dòng)成因分析

河北局所管轄流體井口數(shù)量多，一般情況下研究井-含水層系統(tǒng)時(shí)通常會(huì)忽略含水層與隔水層之間的水流交換，但實(shí)際上幾乎所有井所揭露的含水層都會(huì)與其上面覆蓋的隔水層以及下面接觸的隔水層都會(huì)發(fā)生微弱的水流交換。為更深入研究井-含水層水流運(yùn)動(dòng)情況，充分考慮含水層與隔水層之間的水流交換，建立井-含水層水流徑向、垂向運(yùn)動(dòng)模型是十分必要的。

據(jù)現(xiàn)有研究，潮汐井水位振幅和相位受到裂隙產(chǎn)狀和排水兩類因素的影響。裂隙產(chǎn)狀影響下，潮汐M2和O1波的水位-引潮高振幅比比值差異不大，但位相差符號(hào)相反[7]。排水影響分為徑向和垂向排水兩種情況。在多孔介質(zhì)含水層中，不排水條件下含水層潮汐孔壓-引潮高振幅比為常數(shù)E。在飽水多孔介質(zhì)中，井水位與含水層孔壓相等，水位-引潮高振幅比M=E（對(duì)含水層）或M′=E′（對(duì)隔水層）。

3 資料預(yù)處理

收集無極冀20井2009—2019年水位記錄整點(diǎn)值，剔除因模擬水位換紙引起的錯(cuò)誤數(shù)據(jù)，并采用克里金（kriging）插值[13]對(duì)數(shù)據(jù)進(jìn)行插值補(bǔ)全。因井水位數(shù)據(jù)所包含的信息復(fù)雜，受環(huán)境、場(chǎng)地、認(rèn)知因素影響較大，為了確定潮汐對(duì)井水位的影響程度，首先對(duì)井水位整點(diǎn)值數(shù)據(jù)（2015年8月）進(jìn)行頻譜分析（FFT變換），得到影響井水位的主要頻率段（圖3～4）。圖中顯示頻率為0～0.5 h-1的信息振幅達(dá)0.6 kpa，經(jīng)分析無極冀20井影響水位觀測(cè)的因素，認(rèn)為可能與井的固有頻率及井口周邊環(huán)境有關(guān) ，有待進(jìn)一步核實(shí)。

圖3 2015年8月井水位整點(diǎn)值曲線

圖4 井水位頻譜分析

由上圖分析可知，觀測(cè)數(shù)據(jù)中明顯含有日潮、半日潮、三分之一潮信息，且半日潮振幅最大、日潮次之、三分之一潮的振幅最小，這與固體潮理論是一致的。無極冀20井同震響應(yīng)良好，許多學(xué)者已經(jīng)做過研究[14-15]，不再敘述。

4 井水位潮汐特征分析

基于井水位潮汐變化，通過baytap-G程序、取計(jì)算窗長720 h、滑動(dòng)步長168 h[1]計(jì)算無極冀20井2009—2019年水位整點(diǎn)值數(shù)據(jù)，得到M2和O1波的潮汐因子（各個(gè)諧波的觀測(cè)振幅與理論振幅之比）、相位和振幅及其誤差。由計(jì)算結(jié)果可知：M2波計(jì)算誤差最小，其振幅誤差主要分布在0～0.08 GPa之間，相位誤差在0°～2.1°之間；O1波振幅誤差主要分布在0～0.24 GPa，相位誤差在1°～16°之間。

繪制M2波變化特征（圖5）、M2波和O1波相位變化曲線（圖6）。由圖5可以看出，2009—2019年M2波潮汐因子、相位的動(dòng)態(tài)變化可能與這期間中強(qiáng)地震的發(fā)生有關(guān)（地震已用箭頭標(biāo)出）。地震基本情況見表1，表中能量密度是根據(jù)前人推導(dǎo)的地震波能量密度公式計(jì)算得出[16-17]。

圖6 M2波、O1波相位動(dòng)態(tài)變化圖

5 井-含水層水流方向改變分析與檢驗(yàn)

5.1 井-含水層水流方向改變分析

實(shí)際情況下，井-含水層系統(tǒng)會(huì)受到排水情況以及裂隙產(chǎn)狀共同影響，潮汐參數(shù)表現(xiàn)形式復(fù)雜[7-8]（表2）。

由圖5結(jié)合表1～2可得出，2009—2019年無極冀20井的M2波潮汐因子、相位動(dòng)態(tài)變化分為4個(gè)階段。

第一階段為2009—2011年。初潮汐因子、相位變化趨勢(shì)不穩(wěn)定，波動(dòng)較大。由表2可以看出2009年9月30日蘇門答臘7.7級(jí)地震、2010年4月7日北蘇門答臘8.0級(jí)地震前后，M2波、O1波相位均為正值，M2波潮汐因子稍大于O1波潮汐因子，E'/E＜1，說明隔水層和含水層之間發(fā)生水流交換。由表1可知，這種情況井水位為垂向排水，水流運(yùn)動(dòng)方向?yàn)榇瓜颉?/p>

表1 幾次中強(qiáng)地震基本情況

表2 井水位變化特征成因分析

由圖5～6可以看出，2009—2010年M2波、O1波相位值多次出現(xiàn)正負(fù)號(hào)的轉(zhuǎn)變，且數(shù)值變化幅度較大。根據(jù)劉春平等[8]研究結(jié)果，這種情況的發(fā)生可能跟隔水層裂隙發(fā)育產(chǎn)狀有關(guān)。2009年蘇門答臘南部7.7級(jí)地震、2010年北蘇門答臘8.0級(jí)地震，這兩次地震M2波潮汐因子、相位均出現(xiàn)變化，且地震能量密度最大為2.38×10-4J/m3，與Wang[16]和Manga[17]給出的引起潮汐因子和相位發(fā)生變化的能量密度下限為10-3J/m3不符，此現(xiàn)象可能與裂隙發(fā)育有關(guān)。綜合得出，此時(shí)間段內(nèi)水流以垂向運(yùn)行為主，受裂隙影響。

圖5 M2波潮汐因子、相位、振幅隨時(shí)間變化曲線

第二階段為2011—2014年。2011年3月11日日本本州9.0級(jí)地震發(fā)生后M2波、O1波相位由正值變?yōu)樨?fù)值，發(fā)生正負(fù)號(hào)的改變，且兩波數(shù)值差異不大，說明井水位主要受徑向排水。2011—2014年O1波潮汐因子均稍小于M2波潮汐因子，E'/E＜1說明該井水位變化主要受徑向排水為主、垂向排水為輔，水流主要運(yùn)動(dòng)方向?yàn)閺较?，也說明震后含水層滲透系數(shù)得到改變并導(dǎo)致孔隙壓力重新分布。

為了更進(jìn)一步驗(yàn)證水流運(yùn)動(dòng)方向的改變，計(jì)算了3月11日日本本州9.0級(jí)地震的地震波能量密度為4.21×10-2J/m3。根據(jù)前人給出的研究結(jié)果：能引起含水層儲(chǔ)水系數(shù)、滲透系數(shù)等水文地質(zhì)參數(shù)改變的地震波能量密度最小數(shù)量級(jí)為10-3J/m3，這說明本次日本本州地震足以引起井-含水層水文地質(zhì)參數(shù)的改變，從而影響水流的運(yùn)動(dòng)方向。

第三階段2015年尼泊爾地震后至2018年下半年，2015年4月25日尼泊爾地震的發(fā)生使得M2波、O1波相位值均由負(fù)變?yōu)檎?，說明無極冀20井含水層與其相鄰的隔水層之間發(fā)生了水流的交換，含水層在壓應(yīng)力的作用下向其上覆或者下伏隔水層排水。在E'/E＜1條件下，M2波潮汐因子均稍大于O1波潮汐因子，理論上排除了徑向排水的可能，進(jìn)一步證實(shí)了該井主要受垂向排水影響，水流運(yùn)動(dòng)方向?yàn)榇瓜?。此現(xiàn)象維持至2018年8月19日斐濟(jì)地震發(fā)生后，說明尼泊爾地震的發(fā)生使得無極冀20井由徑向排水變?yōu)榇瓜蚺潘?。同樣，?jì)算此次地震波能量密度為1.399×10-3J/m3，符合前人研究結(jié)果。

由圖5可知，尼泊爾地震發(fā)生后井水位M2波相位、振幅均變大，而O1波振幅變化不明顯，表明該井井水位潮汐響應(yīng)滿足過渡區(qū)的特點(diǎn)。2016年所羅門群島地震前后M2波和O1波振幅和相位沒有明顯變化，且與尼泊爾地震后基本一致。根據(jù)Doan等[18]建立的不同頻率波振幅和相位值與滲透系數(shù)的模型，可知尼泊爾地震前無極冀20井滿足排水條件（徑向），震后該井滿足不排水條件（徑向），再次論證了尼泊爾地震后井水位受垂向排水影響，水流運(yùn)動(dòng)方向?yàn)榇瓜颉?/p>

第四階段，由2018年下半年開始，M2波、O1波相位調(diào)整至日本本州地震后的負(fù)值水平，潮汐因子、振幅變化也趨于平穩(wěn)。斐濟(jì)地震能量密度為2.229×10-5J/m3，能量水平不足以改變M2波的相位。通過計(jì)算2018年下半年全球MS≥5.0地震的能量密度均不足以引起M2波相位改變，因而初步推測(cè)，這可能是震后區(qū)域靜態(tài)應(yīng)力的調(diào)整過程。調(diào)整后兩波相位均為負(fù)值，且差異不大，O1波潮汐因子稍小于M2波，說明應(yīng)力場(chǎng)調(diào)整后的井水位主要受徑向排水影響，E'/E＜1條件下垂向排水為輔，水流主要運(yùn)動(dòng)方向?yàn)閺较?。此現(xiàn)象一直持續(xù)到2019年末，根據(jù)M2波、O1波動(dòng)態(tài)變化，可推斷出這種情況將會(huì)維持一段時(shí)間，直至有新的大能量地震的發(fā)生。

繪制2009—2019年無極冀20井M2波和O1波的潮汐因子動(dòng)態(tài)變化曲線（圖7）。根據(jù)劉春平等[8]給出的理論，結(jié)合圖6、表2可看出：2012—2014年M2波潮汐因子基本大于O1波，說明無極冀20井在該時(shí)間段屬于表1中所列第1、2種情況；M2波相位基本小于O1波，因此該時(shí)間段內(nèi)該井以徑向排水為主，E'/E＜1條件下垂向排水為輔，水流運(yùn)動(dòng)方向主要為徑向；2016年反之，為垂向排水，水流運(yùn)動(dòng)方向?yàn)榇瓜颉Ｅc上述分析結(jié)果一致。

圖7 M2波、O1波潮汐因子動(dòng)態(tài)變化圖

5.2 水流運(yùn)動(dòng)方向檢驗(yàn)

由Hsieh等[4]建立的僅徑向排水情況下井-含水層響應(yīng)模型可知，M2波振幅系數(shù)與相位呈同向變化，而Doan等[18]得出僅垂向排水情況下井-含水層模型得到，M2波振幅系數(shù)與相位呈反向變化據(jù)此對(duì)上述分析結(jié)果進(jìn)行檢驗(yàn)。

繪制無極冀20井2012—2014年、2016年M2波潮汐因子和相位動(dòng)態(tài)曲線圖（為避免區(qū)域靜態(tài)應(yīng)力場(chǎng)調(diào)整干擾，本文只繪制了2016年潮汐參數(shù)曲線）（圖8～9）。由圖可知無極冀20井2012—2014年M2波潮汐因子和相位呈同向變化，主要為徑向排水，水流運(yùn)動(dòng)方向主要為徑向；2016年M2波潮汐因子和相位呈同反變化，為垂向排水，水流運(yùn)行方向?yàn)榇瓜?。進(jìn)一步驗(yàn)證了上述結(jié)論。

圖8 2012—2014年M2波潮汐因子和相位動(dòng)態(tài)曲線

圖9 2016年M2波潮汐因子和相位動(dòng)態(tài)曲線

6 結(jié)語

本文通過對(duì)無極冀20井水位資料的潮汐參數(shù)分析，得到以下結(jié)論和初步認(rèn)識(shí)。

1）通過對(duì)無極冀20井水位數(shù)據(jù)FFT變換可得出，該井水位變化主要受潮汐影響，觀測(cè)質(zhì)量較高。

2）無極冀20井潮汐因子和相位的動(dòng)態(tài)變化可能與遠(yuǎn)距離強(qiáng)震有關(guān)。通過對(duì)無極冀20井2009—2019年水位整點(diǎn)值潮汐參數(shù)分析認(rèn)為：①2009—2010年M2波和O1波相位值多次出現(xiàn)正負(fù)號(hào)的轉(zhuǎn)變，且兩波數(shù)值大小不定，表明此時(shí)間段內(nèi)水流以垂向運(yùn)行為主，受裂隙影響。②日本本州地震能量密度為4.21×10-2J/m3，震后水流運(yùn)動(dòng)方向由垂向變?yōu)閺较颍f明震后含水層滲透系數(shù)得到改變并導(dǎo)致孔隙壓力重新分布；③尼泊爾地震能量密度為1.399×10-3J/m3，震后水流運(yùn)動(dòng)方向由徑向變?yōu)榇瓜?，含水層和隔水層之間發(fā)生水流交換，在壓應(yīng)力作用下含水層向隔水層排水；④2018年下半年井水位主要受徑向排水影響，E'/E＜1條件下垂向排水為輔，水流主要運(yùn)動(dòng)方向?yàn)閺较颉?/p>

3）根據(jù)前人不同的研究方法、模型對(duì)所得結(jié)論進(jìn)行檢驗(yàn)，結(jié)果符合預(yù)期。

4）2009—2010年M2波和O1波相位值多次出現(xiàn)正負(fù)號(hào)的轉(zhuǎn)變，根據(jù)前人研究結(jié)果，這可能與水流運(yùn)動(dòng)受裂隙影響有關(guān)，符合無極冀20井的地質(zhì)構(gòu)造條件。