魏海濱，谷洪彪，孔慧敏，3，遲寶明，3

(1.防災(zāi)科技學(xué)院，河北 三河 065201；2.河北省地震動力重點實驗室，河北 三河 065201；3.中國地震局工程力學(xué)研究所，黑龍江 哈爾濱 150006)

0 引言

井水位的同震變化反映了脆性地殼受到地震的影響，導(dǎo)致含水層介質(zhì)的應(yīng)力狀態(tài)發(fā)生改變，使井水位呈現(xiàn)出階躍、振蕩或持續(xù)性變化。地震時井水位表現(xiàn)出的不同響應(yīng)方式表明其物理機(jī)制復(fù)雜且影響因素較多，如井震距、含水層特性及井孔結(jié)構(gòu)等。靜態(tài)應(yīng)變理論及含水層固結(jié)理論對井水位同震階躍變化進(jìn)行了合理解釋(Wakita，1975；Roeloffs，1996；Quilty，Roeloffs，1997；Grecksch，Roth，1999；Jónsson，2003；Itaba，2008；Wang，2009；Wang，Manga，2009)；地震波通過時造成含水層中裂隙的開合、裂隙中膠體顆粒的堵塞與疏通或是孔隙中氣泡的遷移等導(dǎo)致含水層的滲透性發(fā)生變化，可能是造成井水位同震振蕩與持續(xù)性變化的原因(Roeloffs，Evelyn，1998；Elkhoury，2011；Manga，2012；史浙明，2015)。地震造成的靜態(tài)應(yīng)力與峰值動態(tài)應(yīng)力均會對井水位變幅產(chǎn)生影響，Ge和Stover(2000)認(rèn)為與孔壓擴(kuò)散時間尺度相比地震波引起的動態(tài)應(yīng)力會在含水層中迅速衰減，且在地震波通過后含水層一般不會產(chǎn)生永久應(yīng)變，故本文忽略動態(tài)應(yīng)力的影響，主要關(guān)注同震靜態(tài)應(yīng)力引起的含水層與井孔水動力響應(yīng)過程。

利用靜態(tài)應(yīng)變理論研究近場井水位的同震響應(yīng)，涉及不同的時間與空間尺度，是一個復(fù)雜的三維流固耦合過程，該過程利用解析法求解具有極大的挑戰(zhàn)性。數(shù)值模擬是目前研究地震水動力學(xué)問題的重要手段，數(shù)值模擬過程中利用多種軟件的耦合計算來描述同震水位響應(yīng)過程已有一些成功案例(Ge，Stover，2000；李悅，2011；Nespoli，2016；Zhang，2020；Lei，2021)。然而，近場同震水位響應(yīng)過程物理機(jī)制復(fù)雜、影響因素眾多，利用數(shù)值模擬方法模擬近場同震水位響應(yīng)仍存在一些問題。首先，這些數(shù)值模擬模型均未考慮井-含水層系統(tǒng)內(nèi)地下水由孔隙或裂隙介質(zhì)中的流動變?yōu)榫字屑兞黧w域中的流動這一過程，其次已有研究僅使用簡單的孔壓擴(kuò)散方程對研究區(qū)震后孔壓恢復(fù)過程進(jìn)行模擬，未考慮流固耦合過程。

針對上述問題，本文基于Okada位錯理論計算2014年魯?shù)?.5地震同震靜態(tài)應(yīng)變場分布，建立應(yīng)力場-滲流場耦合控制方程組，以魯?shù)榈卣鹜痨o態(tài)應(yīng)變場為初始條件，模擬魯?shù)榈貐^(qū)(200×200)km孔壓的響應(yīng)及恢復(fù)過程，最后利用達(dá)西定律與Navier-Stokes方程將含水層內(nèi)孔壓擾動值轉(zhuǎn)化為會澤井孔內(nèi)水位變幅值。本文在一定程度上還原會澤井水位對魯?shù)榈卣鸬耐痦憫?yīng)過程，分析近場地震活動與井水位變化之間的關(guān)系，為其它近場地震引發(fā)的井水位變化研究提供思路借鑒，也為今后探索更為復(fù)雜的流-固-熱-化學(xué)多物理場耦合井水位同震響應(yīng)的模擬方法提供參考。

1 魯?shù)榈卣鸺皶删桓艣r

1.1 魯?shù)榈卣鸢l(fā)震構(gòu)造

2014年8月3日16時30分云南魯?shù)辇堫^山鎮(zhèn)發(fā)生6.5強(qiáng)震，震中位置(27.11°N，103.35°E)，震源深度12 km。本次地震發(fā)生在巴顏喀拉塊體、川滇塊體和華南塊體交匯區(qū)南部，馬邊斷裂南段、昭通斷裂、蓮峰斷裂之間。由于本次地震造成的斷層破裂長度有限且未達(dá)到地表，因此破裂面方位不易識別。震后許多學(xué)者及機(jī)構(gòu)通過地震烈度分布和地震破裂過程反演對魯?shù)榈卣鸬陌l(fā)震斷層進(jìn)行了研究，得出魯?shù)榈卣鸢l(fā)震構(gòu)造為一條新發(fā)現(xiàn)的近北西走向的斷層——包谷垴—小河斷裂(房立華等，2014；李西等，2014；徐錫偉等，2014；李艷娥等，2015)。但又有學(xué)者通過視震源時間函數(shù)分析、余震分布及震源機(jī)制解研究提出魯?shù)榈卣鹂赡苁怯山鼥|西向和近南北向的共軛斷層先后破裂導(dǎo)致(張勇等，2015)。

本文研究參照中國地震局對于魯?shù)榈卣鸬某醪椒治鰣蟾?Cheng，2015)，魯?shù)榈卣馂橐淮胃邇A角左旋走滑地震，發(fā)震斷層為走向340°的包谷垴—小河斷裂(圖1)。

圖1 魯?shù)榈卣鹫鹬屑捌涓浇貐^(qū)發(fā)震構(gòu)造圖

1.2 會澤井水位觀測

云南會澤井(滇01號井)為靜水位觀測井，地理坐標(biāo)為(26.52°N，103.15°E)。井孔地處川滇地塊東南部，位于小江、蓮峰、則木河3大斷裂的交匯部位，構(gòu)造位置特殊，井孔周邊水文地質(zhì)條件如圖2a所示。會澤井受白霧街不對稱雙曲弧形構(gòu)造控制，張性裂隙發(fā)育，為地下水的富集提供了有利條件，井孔周邊致密玄武巖柱狀節(jié)理發(fā)育，易于形成地下水的運(yùn)移通道，加上兩側(cè)巖溶水的側(cè)向補(bǔ)給以及大氣降水的天然補(bǔ)給，使得井孔所處含水層具有較好的補(bǔ)、徑、排條件，能夠形成一定規(guī)模的含水系統(tǒng)。會澤井深103.15 m，井孔半徑54 mm，套管深度87.8 m，水位埋深約30 m，濾水管為34.06～87.80 m，觀測段為34.06～103.15 m，觀測層巖性主要為第四系風(fēng)化玄武巖(圖2b)，地下水類型為裂隙承壓水，水溫約16 ℃(Sun，2019)。

會澤井自2012年開始數(shù)字化觀測，水位對其周邊幾次地震均有響應(yīng)，映震性較好(圖3a)。多年來，水位整體呈緩慢降低型，水位下降速率約為0.3 m/a，大氣降水對水位影響較小，氣壓效應(yīng)不明顯，有固體潮效應(yīng)。在2014年魯?shù)?.5地震時，距離震中約70 km的會澤井水位有同震響應(yīng)，水位在地震發(fā)生時瞬時上升0.33 m，之后開始下降，并在震后50 d內(nèi)逐步恢復(fù)至震前水平(圖3b)。

圖3 2012—2014年會澤井水位變化(a)和2014年魯?shù)镸S6.5地震會澤井水位的同震響應(yīng)(b)

2 魯?shù)榈卣鹜痨o態(tài)應(yīng)變場

假設(shè)地震造成井水位的變化是由地震靜態(tài)應(yīng)變場導(dǎo)致基質(zhì)骨架變形產(chǎn)生的應(yīng)力，從巖石向孔隙流體傳遞引起的。本文采用美國地質(zhì)調(diào)查局(USGS)開發(fā)的Coulomb 3.3軟件，利用Okada(1992)有限矩形源模型對魯?shù)榈卣鹪斐傻耐痨o態(tài)應(yīng)變場進(jìn)行計算，參數(shù)選取見表1，得到2014年魯?shù)?.5地震同震靜態(tài)應(yīng)變場等值線圖(圖4)。由圖4可見，紅色區(qū)域應(yīng)變值為正值，屬于拉張區(qū)域，藍(lán)色區(qū)域應(yīng)變值為負(fù)值，屬于壓縮區(qū)域。應(yīng)變拉張與壓縮區(qū)域呈現(xiàn)四象限分布，極值分布在斷裂的兩側(cè)，量級為10，遠(yuǎn)離斷層應(yīng)變逐漸減小。

表1 同震靜態(tài)應(yīng)變場計算參數(shù)

圖4 2014年魯?shù)镸S6.5地震同震靜態(tài)應(yīng)變場分布(斷裂名稱同圖1)

3 井水位同震響應(yīng)數(shù)值模擬

地震造成地殼產(chǎn)生形變，使得含水層骨架應(yīng)力向孔隙中的流體傳遞，導(dǎo)致含水層中孔壓受到擾動并在壓力梯度的驅(qū)動下擴(kuò)散恢復(fù)，造成地下水從井孔流入或流出，最終反映為井水位的變化。水位同震響應(yīng)主要包括地震造成的靜態(tài)應(yīng)變場分布、含水層孔壓擾動及擴(kuò)散、壓力擾動下井-含水層系統(tǒng)水位響應(yīng)3個部分。

3.1 含水層孔壓擾動及擴(kuò)散數(shù)學(xué)模型

含水層中孔壓受地震靜態(tài)應(yīng)力的影響，在壓力梯度的驅(qū)動下擴(kuò)散恢復(fù)，這個過程中含水層的彈性變形和孔隙水壓的擴(kuò)散是一個流固耦合的過程，因此，需要建立含水層應(yīng)變與孔壓之間的耦合數(shù)學(xué)模型來定量解釋二者之間的關(guān)系。本文對含水層孔壓擾動及擴(kuò)散控制方程的建立和求解作如下假設(shè)：①含水層為均質(zhì)各向同性、連續(xù)的承壓含水層，頂、底板水平；②多孔介質(zhì)骨架為各向同性彈性體，骨架變形為小變形，符合線性孔彈性理論；③地下水流符合達(dá)西定律；④地下水流為等溫滲流；⑤滲透率各向同性；⑥單元體內(nèi)流體無外部源匯項。

本文只關(guān)注靜態(tài)應(yīng)力對含水層孔壓的擾動，忽略加速度的影響，在含水層中選取典型單元體并進(jìn)行受力分析。利用靜力平衡方程、廣義胡克定律、有效應(yīng)力原理、達(dá)西定律與質(zhì)量守恒定律，建立含水層中任一單元體在受到靜態(tài)應(yīng)力作用下應(yīng)變與孔壓的耦合控制方程組，具體推導(dǎo)過程見魏海濱(2021)的研究：

(1)

圖5 含水層孔壓擾動及擴(kuò)散模型范圍及震中、井孔分布圖

3.2 壓力擾動下井-含水層系統(tǒng)水位響應(yīng)數(shù)學(xué)模型

在井-含水層系統(tǒng)中，當(dāng)?shù)叵滤趬毫μ荻鹊尿?qū)動下由含水層向井中流動時，水流方向由縱向變?yōu)榇瓜?，且井中為純流體域。因此，地下水在井中的運(yùn)動有別于水流在含水層中的運(yùn)動，在井-含水層系統(tǒng)中孔隙水壓擾動從壓力擾動源到井中水位變幅，應(yīng)分為含水層中壓力擴(kuò)散與井中水柱運(yùn)動兩部分。為了刻畫在靜態(tài)應(yīng)力的影響下孔壓在含水層中擴(kuò)散并在井中以水柱的形式運(yùn)動這一過程，本文建立了以井為中心的二維軸對稱的小尺度模型(圖6a)。

該模型包含承壓含水層與井兩部分，井孔半徑54 mm，含水層概化為厚100 m、半徑100 m。在側(cè)邊界施加地震靜態(tài)應(yīng)力導(dǎo)致的孔壓變化值，邊界處孔壓變化值的作用使得含水層內(nèi)產(chǎn)生孔壓梯度。在孔壓梯度的驅(qū)動下，含水層中的水流向井或井水中流向含水層運(yùn)動(井孔只在側(cè)壁進(jìn)水流速為)(圖6b)。當(dāng)?shù)叵滤诔袎汉畬又羞\(yùn)動時，假設(shè)含水層均質(zhì)各向同性，頂、底板水平，水流符合達(dá)西定律：

(2)

式中：為貯水率；為徑向滲透系數(shù)。

當(dāng)?shù)叵滤M(jìn)入井孔后，水流由徑向流動變?yōu)榇瓜蛄鲃樱揖诪榧兞黧w域，故采用Navier-Stokes方程描述井孔中水流的運(yùn)動：

(3)

式中：為井孔中水流速度；g為重力加速度；為孔壓；和分別為水的密度和動力粘滯系數(shù)。

在初始時刻，含水層中各點孔壓為靜水壓力，即=g；邊界條件為：=g+Δ。

(a)壓力擾動下含水層-井系統(tǒng)水位響應(yīng)二維軸對稱模型圖

(b)模擬水位同震變化過程示意圖

3.3 數(shù)值模擬過程

將魯?shù)榈卣痨o態(tài)應(yīng)變場轉(zhuǎn)化為含水層孔壓擾動值，以此作為初始條件，利用COMSOL Multiphysics中的PDE模塊，在200 km×200 km×0.5 km的研究區(qū)范圍內(nèi)對建立的含水層孔壓擾動及擴(kuò)散數(shù)學(xué)模型進(jìn)行求解，獲得研究區(qū)0.1 km深度處含水層孔壓演化規(guī)律。最后，在COMSOL中利用內(nèi)置的達(dá)西定律、層流模塊，建立以井孔為中心、半徑100 m的小尺度井-含水層系統(tǒng)模型，以靜水壓強(qiáng)為初始條件，利用含水層孔壓擾動及擴(kuò)散模型中會澤井孔所處剖分網(wǎng)格的平均孔壓擾動值作為邊界條件，對會澤井水位的變化趨勢進(jìn)行模擬。

3.3.1 含水層孔壓擾動及擴(kuò)散數(shù)值模擬

基于研究區(qū)概念模型(圖5)構(gòu)建三維幾何模型用于數(shù)值模擬，將式(1)表示為廣義矩陣形式，輸入COMSOL Multiphysics的PDE模塊：

(4)

其中，

研究區(qū)淺層含水層巖性主要為二疊系峨眉山組風(fēng)化玄武巖，本文利用等效連續(xù)介質(zhì)理論，將裂隙介質(zhì)等效為多孔介質(zhì)，參照水文地質(zhì)手冊經(jīng)驗值對模型進(jìn)行賦值，參數(shù)設(shè)置見表2。結(jié)合模型初始條件和邊界條件，確定可能影響含水層孔壓擴(kuò)散的因素，設(shè)計不同情景模擬魯?shù)榈卣鸷蠛畬涌讐簲U(kuò)散的過程，模擬情景設(shè)置見表3。采用自由四面體網(wǎng)格剖分計算網(wǎng)格，考慮計算機(jī)的算力限制及計算效率，在保證模型計算精度的前提下共剖分5 417 817個網(wǎng)格，模擬時間步長為1 d，模擬時段為研究區(qū)0.1 km深度孔壓擾動值消散99%以上所需時長。

表2 數(shù)值模擬參數(shù)

表3 不同模擬情景參數(shù)

以魯?shù)榈卣鹜痨o態(tài)應(yīng)變場作為初始條件，對不同模擬情景下研究區(qū)含水層孔壓擾動及擴(kuò)散過程進(jìn)行模擬。初始時刻，研究區(qū)孔壓分布與同震靜態(tài)應(yīng)變場分布一致，但符號相反，最大正、負(fù)超孔隙水壓變化為分別3.52×10Pa、-3.74×10Pa。研究區(qū)內(nèi)距離斷層兩端較遠(yuǎn)的區(qū)域產(chǎn)生的同震孔壓擾動值在±(10～10)Pa，相較斷層兩端附近的孔壓擾動值小2～3個數(shù)量級。

情景1、2模擬不同滲透系數(shù)條件下含水層孔壓演化規(guī)律，如圖7所示。模擬結(jié)果表明，斷層兩端的孔壓梯度決定了孔壓的擴(kuò)散模式。地下水沿著斷層?xùn)|西兩側(cè)的壓力最大梯度流動。情景1中含水層孔壓經(jīng)過50 d的演化，孔壓變化的最大值已降低了99%以上，孔壓基本恢復(fù)至震前水平。情景2中魯?shù)榈卣鹪斐傻目讐簲_動值在低滲透系數(shù)模型中經(jīng)過500 d消散99%以上，基本恢復(fù)至震前水平。

圖7 在情景1(a)和情景2(b)下模擬不同時間研究區(qū)0.1 km深度孔壓分布

通過比較情景1、2模擬結(jié)果可知，兩個模型的孔壓消散模式類似，但當(dāng)滲透系數(shù)降低一個量級時，較低的滲透系數(shù)使孔壓消散的時間增加了約10倍。因此，滲透系數(shù)是地震后孔壓擴(kuò)散時長的主控因素，研究區(qū)地質(zhì)介質(zhì)類型是決定地震后孔壓擴(kuò)散速率的主要變量。

情景3、4分別將孔隙度設(shè)置為0.2、0.6，其它條件與情景1保持一致，模擬不同孔隙度對研究區(qū)孔壓擴(kuò)散規(guī)律的影響。模擬結(jié)果顯示(圖8)，不同孔隙度并不會影響孔壓消散模式，在相同時間內(nèi)不同的孔隙度只會對孔壓的消散速率有所影響但影響有限，孔隙度越大孔壓消散越慢，不同孔隙度下含水層孔壓基本都在50 d內(nèi)恢復(fù)至震前水平。

情景5、6分別將楊氏模量設(shè)置為8×10和8×10，其它條件與情景1保持一致，模擬不同風(fēng)化程度的玄武巖含水層對研究區(qū)孔壓擴(kuò)散規(guī)律的影響。模擬結(jié)果顯示(圖9)，楊氏模量會顯著改變研究區(qū)孔壓擴(kuò)散的時間，在其余相同的條件下，楊氏模量越大，孔壓擴(kuò)散時間越短，風(fēng)化嚴(yán)重的玄武巖孔壓會在350 d內(nèi)恢復(fù)至震前水平，而未風(fēng)化的玄武巖孔壓會在20 d內(nèi)恢復(fù)至震前水平。

實際統(tǒng)計結(jié)果顯示，研究區(qū)內(nèi)多口監(jiān)測井水位基本均在50 d內(nèi)恢復(fù)至震前穩(wěn)定值，因此，情景1、3、4模擬得到的孔壓恢復(fù)時間較為吻合。由模擬結(jié)果可知，含水層滲透系數(shù)與楊氏模量會顯著改變孔壓擴(kuò)散時長，而孔隙度對含水層孔壓擴(kuò)散時長影響較小。

圖8 在情景3(a)和情景4(b)下模擬不同時間研究區(qū)0.1 km深度孔壓分布

圖9 在情景5(a)和情景6(b)下模擬不同時間研究區(qū)0.1 km深度孔壓分布

3.3.2 壓力擾動下含水層-井系統(tǒng)水位響應(yīng)數(shù)值模擬

獲得研究區(qū)含水層內(nèi)孔壓演化過程后，利用COMSOL Multiphysics中內(nèi)置的達(dá)西定律及層流模塊，以靜水壓強(qiáng)為初始條件，選擇情景1中會澤井孔處網(wǎng)格單元含水層平均孔壓演化曲線，作為壓力擾動下含水層-井系統(tǒng)水位響應(yīng)模型的邊界條件，模擬會澤井水位對魯?shù)榈卣鹜痦憫?yīng)過程及震后恢復(fù)過程，結(jié)果如圖10所示。

圖10 會澤井同震水位變化模擬及實測對比圖

模擬結(jié)果顯示，會澤井水位先呈階躍上升形式，之后在50 d內(nèi)恢復(fù)至震前水平，井水位最大振幅為0.45 m。模擬與實測水位變化趨勢一致，但在水位最大振幅的幅值上有所差異，實測水位最大振幅為0.33 m，較模擬值偏小，可能由于水位同震響應(yīng)實際上升過程中有井損的產(chǎn)生所導(dǎo)致。震后井水位實測值下降速率較緩，并且在下降期間還存在幾次小幅度的上升，查詢當(dāng)?shù)貧庀筚Y料，發(fā)現(xiàn)研究區(qū)在8月和9月降雨較多，均勻的降水使井水位恢復(fù)曲線較為平緩，但幾次暴雨可能使實測水位小幅上升，但也不排除地震或其它構(gòu)造應(yīng)力的影響。本文對井水位恢復(fù)過程的模擬，未考慮降雨及其他構(gòu)造應(yīng)力的影響，因此，模擬與實測水位恢復(fù)曲線有所偏差。

3.3.3 模擬計算中存在的不足

首先，本文雖計算出地震造成的同震靜態(tài)應(yīng)變場，建立含水層孔壓擾動及擴(kuò)散、壓力擾動下含水層-井系統(tǒng)水位響應(yīng)2個耦合的數(shù)值模型，模擬了含水層孔壓擾動、井水位響應(yīng)整個過程，但并未建立表示這一過程的全耦合模型，因此，2個模型耦合項的模擬結(jié)果即含水層內(nèi)孔壓的模擬結(jié)果會對最終結(jié)果產(chǎn)生影響。其次，在含水層孔壓擾動及擴(kuò)散的模擬過程中，未考慮隨深度增加溫度變化對模擬結(jié)果的影響，且未考慮地質(zhì)體的非均質(zhì)性對于孔壓擴(kuò)散的影響，這兩個因素可能會使最終結(jié)果產(chǎn)生誤差。同時，本文未考慮降雨、井孔泄壓、井損以及其它構(gòu)造應(yīng)力對于水位變化的影響，從而導(dǎo)致井水位恢復(fù)過程的模擬與實測結(jié)果有一定偏差。最后，本文未考慮動態(tài)應(yīng)力造成的滲透性變化導(dǎo)致的水位變化，在今后研究中應(yīng)考慮該影響。

4 結(jié)論

本文以會澤井水位對2014年魯?shù)?.5地震同震響應(yīng)為例，利用數(shù)值模擬方法對近場地震水位同震階變的過程及影響因素進(jìn)行分析。建立2個耦合模型模擬計算會澤井水位對魯?shù)榈卣鸬耐痦憫?yīng)過程，結(jié)論如下：

(1)魯?shù)榈卣鹪斐裳芯繀^(qū)內(nèi)地殼壓縮區(qū)與膨脹區(qū)呈四象限分布，壓縮區(qū)與膨脹區(qū)的分界沿斷層走向延伸，極值點分布在斷層北段兩側(cè)，最大壓縮應(yīng)變值為6.56×10、最大膨脹應(yīng)變值為6.99×10。

(2)同震靜態(tài)應(yīng)變場的產(chǎn)生導(dǎo)致研究區(qū)內(nèi)含水層孔壓擾動。本文設(shè)計6種模擬情景對研究區(qū)內(nèi)孔壓的響應(yīng)與擴(kuò)散進(jìn)行模擬，結(jié)果表明，影響含水層孔壓擾動與擴(kuò)散的參數(shù)為滲透系數(shù)與楊氏模量，滲透系數(shù)越大，含水層孔壓擴(kuò)散越快，楊氏模量越小，含水層孔壓擴(kuò)散越慢，同時，研究區(qū)內(nèi)多口井水位的實際恢復(fù)時間，與各情景模擬孔壓恢復(fù)時間進(jìn)行對比，可知研究區(qū)實際水文地質(zhì)參數(shù)與情景1、3、4所設(shè)置參數(shù)相近。

(3)利用壓力擾動下含水層-井系統(tǒng)水位響應(yīng)模型，選取情景1設(shè)置模型參數(shù)，模擬會澤井孔水位同震響應(yīng)及恢復(fù)過程，得出會澤井水位在同震階段階升0.45 m，之后在50 d內(nèi)恢復(fù)至震前水位，模擬結(jié)果與實測值在趨勢上相符，但由于本模型在研究區(qū)內(nèi)只有一口井，未考慮存在多口井泄壓、井損效應(yīng)及降雨對地下水的補(bǔ)給，因此，對水位的變化細(xì)節(jié)的刻畫還不夠精確。