云南會(huì)澤井水位對(duì)2014年魯?shù)镸S6.5地震同震響應(yīng)過(guò)程模擬*

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

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

0 引言

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

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

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

1 魯?shù)榈卣鸺皶?huì)澤井水位概況

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

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

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

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

1.2 會(huì)澤井水位觀測(cè)

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

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

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

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

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

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

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

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

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

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

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

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

(1)

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

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

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

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

(2)

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

當(dāng)?shù)叵滤M(jìn)入井孔后，水流由徑向流動(dòng)變?yōu)榇瓜蛄鲃?dòng)，且井孔為純流體域，故采用Navier-Stokes方程描述井孔中水流的運(yùn)動(dòng)：

(3)

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

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

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

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

3.3 數(shù)值模擬過(guò)程

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

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

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

(4)

其中，

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖10 會(huì)澤井同震水位變化模擬及實(shí)測(cè)對(duì)比圖

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

3.3.3 模擬計(jì)算中存在的不足

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

4 結(jié)論

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

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

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

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