王西寶 程樹岐
(中國山東 276400 山東省地震局沂水地震臺)
地下水位觀測值是一個包含氣壓、固體潮、降雨、地應力場等多種因素影響值的復合參數(shù)(車用太等,2006)。地下水位主要受構造應力如地震、斷層蠕動和構造板塊相互作用以及非構造應力如大氣壓力、降雨載荷、潮汐應力等的影響。在這些影響因素中,氣壓和潮汐為長期加載于井孔—含水層的2種自然載荷,降雨和斷層活動則具有較大的隨機性。張昭棟等(2002)研究地下水潮汐現(xiàn)象的物理機制時發(fā)現(xiàn),固體潮、氣壓和地表水體負荷潮汐三者雖然對承壓含水層作用的機理不同,但它們對承壓含水層的應力作用是一致的,在總結水位對體應變潮汐、大氣潮汐與荷載潮異同的基礎上,張昭棟等(2002)將固體潮、氣壓、地表水體負荷潮汐對承壓含水層作用的偏微分方程進行了歸納總結,給出了井水位潮汐的統(tǒng)一數(shù)學物理方程和解析解表達式;Elkhoury等(2006)利用Hsieh等(1987)的方法,分析了美國加利福尼亞地區(qū)2口井水位觀測資料的潮汐振幅比和相位差,認為大震發(fā)生后潮汐振幅比和相位變化是地震波增加了含水層滲透系數(shù)所致;對井水位潮汐振幅和相位進行動態(tài)跟蹤還可能捕捉到地震前兆信息(張昭棟等,1994)。大量觀測和研究表明(陳建民等,1994;張昭棟等,1997,1999a,1999b;廖欣,2010),大震前井水位固體潮可發(fā)生變化。
本文以魯14井水位觀測資料為研究對象,利用Baytap-G潮汐分析軟件(Tamura et al,1991;單維鋒等,2017)計算其潮汐響應特征參數(shù),分析遠場大震前后潮汐因子和相位差的變化特征,以期深入了解井水位微動態(tài)變化特征,為該井水位觀測資料的動態(tài)評價和效能評估提供參考。
魯14井(編號37038)位于山東省莒南縣石蓮子鎮(zhèn)小官莊(35.35°N,118.66°E),當?shù)睾0胃叨?5 m,屬丘陵地區(qū),構造上處于沂沭斷裂帶的昌邑—大店斷裂上。該井水泥固封井口,井深320.20 m。觀測井水位高出地面約6.5 m,套管高7.362 m。套管深度 20 m,0—20 m井徑150 mm,20—100 m井徑110 mm,100—324 m井徑91 mm。地層巖性自上而下分別為,0—1 m第四系覆蓋土層;1—251m 白堊系砂巖、頁巖、礫巖互層;251.9—320.2 m石炭系灰?guī)r夾火成巖侵入體。其中,觀測段為251.9—320.2 m,含水帶類型為裂隙承壓水(圖1)。

圖1 地質構造和井孔柱狀圖Fig.1 Geological structure and wellbore histogram
井水位受氣壓、固體潮、斷層活動、降雨、地應力等多種因素的影響。使用國際上通用的Baytap-G程序提供的方法進行計算分析。該方法采用迭代法估計潮汐參數(shù)、漂移和氣象等時間序列的回歸系數(shù),并且采用貝葉斯信息原理來判斷最優(yōu)參數(shù)估計,從而提高了計算精度(Tamura et al,1991)。其基本原理是將井水位時間序列進行如下分解

式中,M為潮汐分波的個數(shù);Amm、ωmn、φmn分別為第m波群中第n個潮波分量的理論振幅、角頻率和初始相位;δm、Δφm分別為第m波群待估算的振幅因子(潮汐因子)和相位滯后(相位差),其中,當Δφm為正時表示相位超前,Δφm為負時表示相位滯后;Dr(t)為長趨勢項;bk為延遲響應因子;xi為附加信號(外部加載),如氣壓或水文觀測值;h為階變值;tz為階變函數(shù),階變之前其值為0,階變之后其值為1;ε(t)為觀測誤差。右邊第1項為潮汐變化項;第2項為長期趨勢項;第3項是相關聯(lián)數(shù)據(jù)項(本文指氣壓);第4項為階梯變化項;第5項為噪聲項(亦稱殘差項)。
選取2008—2017年魯14井井水位觀測資料整點值進行潮汐參數(shù)的計算,選擇魯14井2008年4月1日至6月30日、2011年2月1日至4月30日井水位觀測數(shù)據(jù)進行各因素分離分析。地震目錄使用李盛樂等開發(fā)的EQDown程序下載,選取國內(nèi)外6級以上并引起潮汐參數(shù)變化的地震(表1)。計算前,對觀測資料中的錯誤數(shù)據(jù)進行剔除;對于觀測中由儀器等引起的短時間缺數(shù)問題,采用mapsis內(nèi)置的插值法進行內(nèi)插處理(蔣駿等,2000;陸遠忠等,2002)。

表1 所選地震的基本信息Table 1 Basic information of selected earthquakes
魯14井為靜水位觀測,觀測所用儀器為中國地震局地震預測研究所研制生產(chǎn)的LN-3A型數(shù)字水位記錄儀,采樣間隔為1 min,經(jīng)過預處理后儀器可產(chǎn)出水位整點值和日均值。據(jù)晁洪太(2008),該井水位的固體潮響應明顯,受固體潮影響的水位最大差值幅度可達23 cm。自2014年4月起魯14井水位下降明顯,同年5月相關專家到莒南魯14井進行水位異常調查,最后確定該井水位下降速率加快的主要原因是受干旱影響,與構造應力調整間的關系不明顯。圖2為魯14井2008年4月1日至6月30日、2011年2月1日至4月30日和2015年3月1日至5月31日井水位觀測數(shù)據(jù)與氣壓及部分遠震的對應曲線。由圖2(a)、2(b)可見,該井水位與氣壓間基本成負相關,而在2015年這種相關性被打破[圖2(c)],隨著水位不斷下降,這種相關性變得不太明顯,這或許與2013年4月后出現(xiàn)淤堵且沒有疏通有關。

圖2 魯14井水位觀測數(shù)據(jù)與氣壓(a)2008年汶川MS 8.0地震;(b)2011年日本MS 9.0地震;(c)2015年尼泊爾MS 8.1地震Fig.2 The curve of water level and air pressure observations
潮汐的強能量主要集中在幾個頻段,如O1、K1等日波頻段,S2、M2等半日波頻段。其中,M2波能量最強,通常具有較大的振幅和較小的計算誤差,且受氣壓和溫度的影響較小(來貴娟,2014)。采用Baytap-G程序分析井水位固體潮效應,該方法可以有效剔除除降水、抽水等非周期或長周期因素的影響,非周期性的偶然干擾也不會影響M2波相位差和振幅的變化趨勢,并且這些外界干擾不會影響關于區(qū)域應力變化的基本結論(唐彥東,2015)。因此,主要分析M2波相位差和潮汐因子的變化趨勢來研究魯14井潮汐參數(shù)變化特征。
為提取魯14井水位潮汐因子和相位差的動態(tài)變化信息,取計算窗長為720 h(30天),滑動步長為360 h(15天),用Baytap-G程序計算M2波的潮汐因子、相位差及計算誤差,并繪制各參數(shù)隨時間變化圖(圖3),部分數(shù)據(jù)見表2。

表2 魯14井水位相位差、振幅、潮汐因子Table 2 Tidal factor and phase difference of the water level of Lu 14 well

圖3 魯14井水位(a)和M2波相位差、潮汐因子(b)隨時間的變化序號1—8表示的地震與表1中地震相對應,序號a表示井出現(xiàn)淤堵Fig.3 The variation of water level,tidal factor,and M2 phase difference
由圖3可見,潮汐因子和相位差表現(xiàn)為同向變化且相位差均小于0°,表明潮汐地下水流類型以徑向流為主。圖4為潮汐因子與相位差之間的關系。由圖4可知,潮汐因子和相位差間呈正相關,表明二者同向變化,這進一步佐證了魯14井潮 汐地下水流類型以徑向流為主。由表2可見,2008年5月12日汶川MS8.0地震后,井水位潮汐因子和相位差同步下降,之后同向緩慢上升,但并沒有恢復到震前水平。直到2010年2月26日琉球群島MS7.0地震后,潮汐因子和相位差才達到比汶川地震前還高的水平,之后約有1年時間基本維持不變。2011年3月11日日本MS9.0地震發(fā)生后,潮汐因子和相位差顯著增大,潮汐因子由0.218 7增大到0.281 2,增加了29%,相位差由-50.771°上升到-34.442°,上升了14.329°。這3次地震導致潮汐因子和相位差同步變化幅度(階升或階降)均較大,并且二者均沒有恢復到震前水平,這可能說明含水層滲透性在一定程度上可能存在永久性改變。

圖4 潮汐因子與相位差之間的關系Fig.4 Relationship between Tidal factor and phase difference
根據(jù)Hsieh等(1987)提出的徑向流條件下井水位對含水層壓力水頭響應模型及其推導出的徑向流條件下井水位—含水層孔壓的振幅比和相位差公式,在徑向流條件下,隨著導水系數(shù)的增加(或減少),潮汐因子和相位差同步增大(或減小)(劉春平,2017)。在上述第1—6次地震后(第7、8兩次地震由于2013年后魯14井出現(xiàn)淤堵且沒有疏通,故不作分析),潮汐因子和相位差均出現(xiàn)不同程度的變化,其中,第1、4、6次地震發(fā)生后潮汐因子和相位差均出現(xiàn)不同程度階降,說明震后導水系數(shù)減弱,由于導水系數(shù)等于滲透系數(shù)與含水層厚度的乘積,而在地震過程中含水層厚度基本保持不變,因此這說明含水層滲透性降低;第2、3、5次地震震后潮汐因子和相位差均同向階升,尤其在第2、3次地震發(fā)生后階升明顯,這表明此時導水系數(shù)增大,即含水層滲透性增強。
經(jīng)過梳理發(fā)現(xiàn),影響魯14井水位潮汐參數(shù)的地震均為震中距大于1 000 km的遠場地震。由于遠場地震引發(fā)的同震體應變過小,無法引起潮汐參數(shù)較大幅度的變化,因此認為,可能是地震波作用使含水層滲透性發(fā)生了變化從而引起潮汐參數(shù)的改變(Elkhoury et al,2006)。地震波的傳播過程即能量釋放的過程,研究表明,地震波能量密度與震中距、震級間存在定量關系(Wang,2007)。對于同一口觀測井,震中距越小,震級越大,則地震波能量密度越大,所能引起的井水位變化越明顯。從引起含水層滲透性變化的幾次地震震級、震中距與能量密度之間的關系來看(圖5),對于魯14井而言,地震波能量密度達到1.00×10-4J/m3時即引起含水層滲透性的變化。Wang 等(2010)研究認為,地震波引起含水層滲透系數(shù)變化的能量密度下限約為1.00×10-4J/m3,當?shù)卣鸩芰棵芏却笥?.00×10-3J/m3時,地震波對裂隙的剪切作用會更強,從而能更顯著地引起含水層滲透系數(shù)的變化。本文研究結果與上述結果基本一致。從圖5可知,引起魯14井滲透性變化的地震波能量密度為1.00×10-4—1.00×10-1J/m3,而引起含水層滲透性增大或減小的地震波能量密度并沒有一定的規(guī)律可循,并非地震波能量密度大(或小)就一定能引起含水層滲透性增大(或減小)。如表1 所示,在幾次地震中第1次地震的地震波能量密度較大,但卻沒有引起含水層滲透性增加;而第5次地震的地震波能量密度(0.47×10-3J/m3)小于地震4的(0.68×10-4J/m3)卻引起含水層滲透性增加,這可能說明地震波能量密度并非是引起含水層滲透性增大或減小的主要因素,含水層滲透性或許還受地震波振幅及低頻面波衰減速度的影響,同時也表明,地震引起的含水層介質滲透性參數(shù)變化還存在各向異性。

圖5 地震波能量密度分布與震級、震中距之間的關系Fig.5 Distribution of earthquake-induced permeability changes as functions of earthquake magnitude and distance
氣壓和固體潮對地下水位產(chǎn)生的影響量可利用Baytap-G程序計算地下水位受大氣壓力、固體潮影響的分量得到。利用Baytap-G程序對魯14井數(shù)據(jù)進行處理,井水位觀測值趨勢項、氣壓對水位的影響量、固體潮對水位的影響量及噪聲對水位的影響量如圖6所示。由圖6可見,氣壓對魯14井影響較小,從2次地震前后的情況來看,氣壓影響分量分別小于1.0 cm、4.5 cm;而潮汐影響分量分別小于9.0 cm、11.3 cm;而從噪聲對水位的影響量來看,2008年4月1日至6月30日噪聲對水位的影響量一般小于0.1 cm,2011年2月1日至4月30日噪聲對水位的影響量小于0.35 cm,正常情況下,若沒有地震活動的影響,噪聲對水位的影響量曲線近似為直線,僅顯示觀測誤差。而在地震發(fā)生時,噪聲對水位的影響量并沒有明顯異常波動,2次地震前后噪聲對水位影響量的不同表現(xiàn)可能說明不同時期系統(tǒng)觀測誤差是不同的。

圖6 水位時間序列分析結果(a)汶川MS 8.0地震;(b)日本MS 9.0地震Fig.5 Analysis results of water level time series
固體潮和氣壓是持續(xù)作用于井-含水層系統(tǒng)的2個最重要的非構造應力,通過對氣壓、固體潮影響量的分離可知,對魯14井來說,固體潮對該井的影響量更大一些。
通過對魯14井水位的分析,得到以下初步認識。
(1)魯14井水位與氣壓間基本成負相關。自2013年4月出現(xiàn)淤堵后這種相關性逐漸消失,這可能與淤堵且沒有疏通有關。
(2)魯14井潮汐地下水流類型以徑向流為主。由潮汐因子和相位差隨時間的變化特征可知,受固體潮影響的潮汐參數(shù)變化相對穩(wěn)定。井水位潮汐因子和相位差受遠場大震的影響較顯著。
(3)地震波不僅可引起魯14井含水層滲透性增強,也可引起滲透性減小。但是,地震波能量密度并非引起含水層滲透性增大或減小的主要因素,含水層滲透性或許還受到地震波振幅及低頻面波衰減速度的影響。地震波能量密度與震級、震中距之間的關系表明,引起該井滲透性變化的能量密度范圍為1.00×10-4—1.00×10-1J/m3。
(4)對氣壓、固體潮影響量的分析表明,固體潮對魯14井的影響量更大一些。
以上僅是一些初步認識,關于地震波對魯14井含水層滲透性影響的研究還需積累更多資料和震例來進行分析。