基于數值模擬與含水層垂向應力反演的靜樂井水位異常分析1

張淑亮 李宏偉 呂 芳 陳 慧

（山西省地震局，太原 030021）

引言

靜樂井水位是山西乃至華北地區映震能力較強的地下流體觀測井之一，井水位異常能客觀反映區域應力場的變化，全球7 級以上地震均能記錄到水震波。1983 年始測以來出現的幾次高值異常與山西北部至晉冀蒙交界地區強震活動時間具有準同步性，在1989 年大同-陽高Ms6.1 級地震、1991 年忻州Ms5.1 級地震、1991 年大同-陽高Ms5.8 級地震、1998 年張北Ms6.2 級地震、1999 年大同-陽高Ms5.6 級等中強地震發生前有較好的異常顯示。因此將靜樂井水位高值異常作為晉冀蒙交界地區地震趨勢預測的重要判據（張淑亮等，1997；車用太等，1999；王愛英等，1998）。 2017 年8 月靜樂井水位又出現大于多年平均值的高值異常，到目前為止異常仍在持續。按照以往震例時間指示意義，目前已接近發震時間，未來晉冀蒙交界區震情形勢如何發展？目前的異常是否真實反映了構造活動加劇所引起的區域應力場和深部地球物理場的改變？是否與中強震孕育過程有關？靜樂井水位目前的異常特征不足以解析異常區的應力狀態，從而也直接影響對晉冀蒙交界地區地震形勢的認識和未來地震趨勢的判定。

多年的觀測表明，影響靜樂井水位動態變化的主要干擾因素為短時強降雨與河流水體的荷載作用。無論是短時間的強降雨還是洪水期河流的荷載作用，均可在短時間內造成井水位的快速變化（張淑亮等，2005；車用太等，2004；魚金子等，1994；張昭棟等，1989；張昭棟等，1990）。目前對這2 種干擾因素的識別與排除基本上以定性分析為主，定量分析方法較少。特別是對于同時受2 種因素影響的靜樂井水位而言，在異常分析與判定工作中，僅采用與同期降雨量對比的方法分析井水位變化與降雨量的關系，而河流水位變化對井水位的影響分析工作至今尚未開展。為此，本文利用數值模擬與含水層垂向應力反演的方法，對降雨與河流荷載作用下的靜樂井區應力進行定量計算，并從區域應力場的角度分析2 種方法計算結果的差異，以期為地震趨勢判定提供較可靠的預測依據。

1 井區概況

靜樂流體觀測井位于呂梁山斷塊隆起區的靜樂向斜，位于東碾河斷裂帶上，地處東碾河南岸，與東碾河相距約130m，其北部和南部均為山區（見圖1）。靜樂井深362.92m，井區為河谷地區，主要發育有下古生界灰巖、白云巖喀斯特含水層及第四系沖洪積砂礫石含水層，2 個含水層間有較強的水力聯系。觀測層巖性時代為中奧陶系灰巖裂隙溶洞水，其頂板埋深為46.3m，單位涌水量84.5L/（s·m），滲透系數為2.41—11.5m/d。觀測層具有厚度大、透水性強的特點。含水層巖石溶洞發育，靜樂井在230—290m 處通過破碎帶，地下水由北西經井區向東南方向排泄，井區處于強徑流區，地下水補給主要來自北、南山區基巖裸露區的大氣降水補給和東碾河河水沿斷裂破碎帶入滲補給。

圖1 井區地質構造 Fig. 1 Well area tectonic map

2 強降雨與河流荷載作用下靜樂井水位動態特征

多年的觀測表明，大于15mm 的集中降雨將引起靜樂井水位的快速上升（張淑亮等，2005），可能由于降雨在地表形成荷載，使得含水層受力狀態發生變化，引起水井水位上升。東碾河水位變化是影響靜樂井水位變化的又一重要因素。當大面積洪水作用通過井區時，井水位迅速大幅上升，幾乎沒有時間滯后，洪水期水位上升，洪水過后水位下降，水位升降幅度取決于洪水深度。

圖2 靜樂井水位與降雨量和東碾河流量的對比 Fig. 2 Comparison of Jingle well level with rainfall and Dongnian river discharge

圖2 所示為靜樂井水位與井區降雨量和東碾河流量的對比。由圖2 可知，2017 年8 月20 日左右單日降雨量為82.8mm 時，東碾河流量由原來的0.777m3/s 增至9.8m3/s，增加了9.023m3/s； 靜樂井水位由9.074m 升至7.747m，水位上升幅度達1.327m。由此可知，降雨與河流荷載作用是影響靜樂井水位動態變化的主要因素，但降雨結束后井水位并未迅速回落至上升前的水平。張淑亮等（1997）曾對靜樂井水位高值變化與降雨的關系進行了較系統的研究，認為降雨量的增大雖可在短時間內造成井水位升降變化，但不能改變水位的年變動態，只有在區域應力場增強的背景條件下，降雨引起的附加力源才可能造成井水位年動態的異常變化，降雨荷載是區域應力場變化的誘發因素，也可能是在2017 年雨季結束后井水位仍處于高值的主要原因。

3 降雨與河流荷載作用下井區應力數值模擬分析

利用井區周圍地殼結構、不同地層介質彈性參數建立三維地質模型，預測降雨與河流荷載作用下的井區應力。由于本文基于數值模擬與含水層垂向應力反演的方法進行定量計算，因此，依據靜樂井井區地層巖性特征和井-含水層特征參數選取計算參數。因靜樂井位于河谷地區，故計算時以變質巖系中大理巖力學性質指標經驗數據為依據。

3.1 三維地質模型建立

根據靜樂井區地層結構、河流位置及研究目標，設定有限元尺寸為1km×1km×360m，河流寬76m。采用SOLID185 單元，模型采用由面拉伸成體的建模方式，取井孔周圍圍巖彈性模量為5×1010MPa，泊松比為0.25，密度為2600kg/m3。

網格劃分的目的是對模型實現離散化，并用適當數量的網格單元得到最精確的解。本文選用Workblench 中的自動網格劃分工具，以滿足模型要求。網格疏密程度直接影響計算結果的精度，但網格加密將增加CPU 計算時間，且需更大的存儲空間。為保證模型網格密度和計算精度，本文選擇的網格尺寸為1.5m，同時適當增加井區附近區域的網格密度。其他參數均采用默認值，可滿足模型的需求。

載荷和約束是ANSYS 軟件求解計算的邊界條件，由所選單元的自由度形勢定義。由于本文主要分析河流與降雨荷載的影響，因此選擇垂向應力加載。為防止模型移動，模型底端與側面采用固定約束，限制發生移動和轉動。

考慮數值模擬結果的多樣性及模型造成的誤差，建模時選擇不同的邊界條件，并反復測試計算，以可反映降雨與最大河流荷載作用對井區的最大影響范圍作為模型的邊界條件。由于降雨與河流荷載的作用主要表現為垂向應力，因此模型初始地應力只考慮重力。因本文通過數值模擬與反演的方法分析靜樂井水位在降雨與河流荷載作用下的垂向應力在量級上是否具有相似性，進而確定影響水位變化的主要因素，對計算精度的要求相對較低，故選取物理力學參數時以研究區內地層巖性的平均值作為模型計算參數。

3.2 結果分析

由于2017 年是近年來靜樂井水位年變幅最大的年份，也是降雨與河流荷載影響最大的年份，故本文以2017 年8 月洪水期東碾河水位、靜樂井區降雨量和靜樂井水位為研究對象，計算降雨和河流荷載作用下的靜樂井井區的垂向位移和垂向應力。

（1）單降雨荷載的影響

以2017 年8 月降雨量作為加載體，將降雨量等效為荷載，計算得到靜樂井垂向位移和垂向應力如圖3 所示。由圖3 可知，降雨荷載對靜樂井產生的垂向位移約為0.6mm，垂向應力約為12hPa。

（2）單河流荷載的影響

以洪水期東碾河水位最大變化量作為加載體，并將河水水位等效為荷載計算得到靜樂井垂向位移和垂向應力如圖3 所示。由圖3 可知，洪水期河流荷載對靜樂井產生的垂向位移約為0.3mm，垂向應力約為15hPa。

圖3 單降雨荷載作用下靜樂井垂向位移和垂向應力 Fig.3 Vertical displacement and vertical stress of Jingle well under single rainfall load

圖4 單河流荷載作用下靜樂井垂向位移和垂向應力 Fig.4 Vertical displacement and vertical stress of Jingle well under single river load

（3）降雨與河流荷載的綜合影響

由于靜樂井水位變化是在因降雨與河流荷載的共同作用，單河流和單降雨的影響僅能反映井水位的部分變化，因此，需考慮二者的共同作用。以2017 年8 月東碾河水位最大變化量和降雨量作為共同加載體，并將二者等效為荷載，計算得到靜樂井垂向位移和垂向應力如圖5 所示。由圖5 可知，降雨和河流荷載共同作用對靜樂井產生的垂向位移約為0.7mm，垂向應力約為30hPa，與單河流和單降雨影響的二者之和基本一致，表明降雨與河流荷載作用是影響靜樂井水位變化的主要因素。

圖5 河流與降雨荷載共同作用下靜樂井垂向位移和垂向應力 Fig.5 Vertical displacement and vertical stress of Jingle well under combined action of river and rainfall loads

4 靜樂井-含水層垂向應力反演

數值模擬結果表明，降雨和河流荷載是影響靜樂井水位變化的主要因素。為進一步驗證數值模擬結果能否真實反映井區實際應力變化，還需將其與實際水位觀測值所反映的井區應力進行對比。在由井水位變化反演井區應力場變化方面，目前取得了一些進展，如黃輔瓊等（2004）利用華北地區40 多個深井水位動態變化資料，探討研究了華北地區現今構造應力場狀態；孫小龍等（2011）運用小波分析法提取出能反映水位多年動態變化的趨勢信息，反演出華北地區多年構造應力場的變化特征；丁風和等（2015）利用部分含水層介質參數、井水位變化量與含水層垂向應力變化量的關系式，定量分析了張渤帶地區構造應力場的動態變化過程。

本文根據靜樂井-含水層介質參數及井水位變化量等數據，反演井水位變化引起的含水層垂向應力。應力反演所用數據時段為2016—2017 年（根據東碾河水位與流量始測時間確定）。反演前首先對觀測數據進行預處理（丁風和等，2015），去除長周期干擾因素、降雨、氣壓等年際變化對水位的影響；然后采用去趨勢法和傅里葉周期法，剔除地下水開采和降水補給等長周期干擾對水位的影響。

4.1 反演所需各類參數計算

（1）氣壓系數與潮汐效率計算

利用卡爾曼濾波方法消除地球固體潮對所選時段水位、氣壓整點值數據的影響。利用濾波后的日均值數據，采用高階差分和一元線性回歸等方法獲得靜樂井氣壓系數。利用維尼迪科夫潮汐調和分析方法獲得靜樂井水位M2波潮汐效率。氣壓系數和潮汐效率計算結果如表1 所示。

表1 靜樂井-含水層多種參數計算結果 Table 1 Calculation results of various parameters of Jingle well-aquifer

（2）靜樂井-含水層介質參數計算

不排水狀態下，井水位氣壓系數和潮汐效率可分別表示為（Bredehoeft，1967；李春洪等，1990；張昭棟等，1993）：

式中BP為井水位氣壓系數；Bg為M2波潮汐效率；n為含水層孔隙度；α為含水層固體骨架體積壓縮系數；β為含水層內水的體積壓縮系數；ρ為水的密度；g為重力加速度（ρg=0.098hPa/mm）。因此，n、β可依據式（3）與表1 中的潮汐效率和氣壓系數得到。根據式（1）或式（2）可求α。

4.2 靜樂井-含水層垂向應力反演

利用表1 數據，根據水平層狀含水層模式下，含水層介質參數（孔隙度、水和固體骨架體積壓縮系數等）、井水位變化量與含水層垂向應力變化量的關系式（張昭棟等，1987）等來計算含水層垂向應力，含水層垂向應力計算如下：

式中 Δσz為含水層垂向應力變化量；E為含水層固體骨架的楊氏模量（E= 1/α）；△H為剔 除地下水開采、降雨和氣壓影響后的含水層應力變化引起的壓力水頭變化量，即井水位變化 量。當井-含水層系統所受應力增強，即 Δσz>0 時，井水位上升，水位埋深H變小，其變化量ΔH<0；當井-含水層系統所受應力減弱，即 Δσz<0 時，井水位下降，水位埋深H變 大，其變化量ΔH＞0。利用式（4）與表1 中有關參數計算靜樂井2016—2017 年含水層垂向應力，結果如表1、圖6 所示。可知，靜樂井含水層垂向應力高值時段與降雨量高值時段、東碾河流量高值時段較吻合，再次表明降雨與河流荷載作用是造成含水層垂向應力高值變化的主要因素，井水位的大幅突升變化與垂向應力增大密切相關。

5 靜樂井水位異常分析

5.1 模型預測值與井水位反演值對比分析

由數值模擬結果可知，無論單河流荷載還是單降雨荷載效應均對洪水期靜樂井水位造成了一定影響，二者造成的垂向位移累計0.9mm，垂向應力累計27hPa。河流和降雨綜合影響模型結果也得到類似結果（垂向位移累計0.7mm，垂向應力累計30hPa）。由含水層垂向應力反演結果可知，洪水期含水層的應力顯著增強，由正常時段的4.2hPa 增至79hPa 左右。2 種方法計算結果均表明靜樂井水位上升變化與河流荷載效應和降雨荷載效應密切相關。但2 種方法計算結果存在一定差異，降雨與河流荷載作用下井區垂向應力數值模擬結果小于含水層垂向應力反演值，相差39hPa。

圖6 東碾河流量、靜樂井-含水層垂向應力與降雨量的對比 Fig. 6 Flow rate of Dongnian river, vertical stress of Jingle well-aquifer and rainfall comparison chart

由于數值模擬計算得到的垂向應力是研究區自身重力、降雨與河流荷載作用下產生的垂向應力之和，而由井水位反演的垂向應力不僅包含上述因素產生的垂向應力，可能還含有一些因構造活動引起的垂向應力變化。事實上2 種方法的計算結果的確存在一定 差異。雖然2 種計算方法均受建模條件和參數選取的影響而存在計算誤差，但根據以往的研究，誤差值一般不會大于計算值，因此2 種計算結果的差異性是客觀存在的。按照表1 中靜樂井水位平均氣壓系數4.3521mm/hPa 可推算出2 種方法垂向應力計算結果差值引起的井水位變化幅度約為16.97cm，大于正常時段水位平均月變化幅度0.9cm 的水平。表明影響靜樂井水位快速上升的因素除降雨與河流荷載作用外，可能還受其他因素的影響。

5.2 水位高值異常與構造活動的關系

為進一步探討靜樂井水位2017年8月以來的高值異常除受降雨與河流荷載影響外是否受構造活動的影響，本文從可反映區域應力狀態變化的參數和測量值（缺震和b 值、地震矩釋放、GPS 基線等）及同一構造區其他前兆測項變化特征與靜樂井水位高值異常進行對比分析。

大量研究結果表明，b 值與其對應區域應力狀態、地殼破裂強度有關（劉雁冰等，2017；朱艾斕等，2005；吳小平，1990）。1983 年以來山西北部至晉冀蒙交界區3 級以上地震缺震和b 值曲線有4 次異常過程（見圖7），第1 次異常出現在1988 年，異常持續過程中在山西地震帶出現一組中強地震活動，先后發生了大同-陽高6.1 級地震、侯馬4.9 級和忻州5.1 級地震、大同-陽高5.8 級地震等，其中大同-陽高6.1 級地震被認為是華北地區在經歷唐山大地震后進入新活躍幕的標志，是華北應力場趨于增強的結果；第2 次異常出現在1995 年初，異常持續過程中先后發生了1996 年5 月包頭6.4 級地震、1998 年1 月張北6.2 級地震、1999年3 月張北5.6 級地震、1999 年11 月大同-陽高5.6 級地震等中強地震活動，既是1989 年開始的華北新一輪活躍幕的繼續，又是華北應力場增強的結果；第3 次異常始于2007 年，這次異常過程中在山西北部至晉冀蒙交界區發生了幾次4 級以上中等地震，部分學者認為，這組中等地震的發生是受汶川地震的影響，鄂爾多斯塊體與華北平原塊體相對擠壓和扭錯增強，導致山西帶形變場與構造應力場由原來的構造拉張轉為構造擠壓，是汶川地震延遲觸發的結果 （劉峽等，2013，朱艾斕等，2010；劉瑞春等，2014）；2017 年又出現同步異常變化，目前異常仍在發展，且缺震和b 值的幾次異常與靜樂井水位高值異常較為同步，因此，可認為靜樂井水位高值異常反映的是華北區域應力場增強的過程，是山西地震帶地震活動增強的重要判據（張淑亮等，1997）。2017 年靜樂井水位高值異常也可能與區域應力場增強有一定關系。

圖7 靜樂井水位月均值與山西地震帶3 級地震缺震和b 值的對比 Fig. 7 Comparison chart of monthly average water level of jingle well with magnitude 3 earthquake deficiency and b value in Shanxi seismic belt

由跨過靜樂井的山西岢嵐-山西太原GPS 基線時序曲線圖8 可知，2016 年以來基線由原來的下降趨勢（基線縮短）變為緩慢上升（基線拉長），應力狀態在2016 年發生了變化，由靜樂井所在區域地震矩加速釋放時圖9 可知，2017 年井區存在明顯的加速釋放低值區。GPS基線和地震矩釋放均表明2017 年靜樂井所在區域應力場在增強。

圖8 山西岢嵐-山西太原GPS 基線時序曲線 Fig.8 GPS baseline timing curve of Kelan, Shanxi-Taiyuan, Shanxi

靜樂井構造上位于山西西部呂梁山隆起區，位于同一構造區距靜樂井較近的前兆測點有寧武鉆孔應變（60km）和神池鉆孔應變（95km），在靜樂井水位突升前后這2 個測點的應變觀測值也出現大幅度反向變化，所以說靜樂井水位高值變化不是獨立事件（見圖10），可能反映了所在構造區應力場的改變。

圖9 靜樂井所在區域地震矩加速釋放時空圖 Fig.9 Accelerated release diagram of seismic moment in jingle well area

綜上所述，靜樂井水位的高值異常除與降雨、河流荷載作用有關外，井區構造活動增強可能也是一種因素。

圖10 靜樂井水位與同一構造區其它前兆測項對比圖 Fig.10 Comparison between Jingle well water level and other precursor measurements in the same structure area

6 結論與認識

通過數值模擬、井-含水層垂向應力反演、區域應力場變化特征等對靜樂井水位2017 年以來出現的高值異常進行初步探討，得出以下結論：

（1）根據井區周圍地殼結構、不同地層介質彈性參數建立三維地質模型，并以2017 年洪水期東碾河水位最大變化量和降雨量作為共同加載體，模擬得到靜樂井垂向位移約為0.7mm，垂向應力約為30hPa。

（2）利用靜樂井含水層孔隙度、固體骨架體積壓縮系數和水的體積壓縮系數等參數，根據水平層狀含水層模式下含水層介質參數、井水位變化量與含水層垂向應力變化量的關系式，得到靜樂井-含水層垂向應力約為79hPa。

（3）數值模擬與井-含水層垂向應力反演結果均表明，靜樂井水位高值異常與同時段降雨量增多、河流荷載效應增強關系密切，但2 種結果存在一定差異，含水層垂向應力反演值大于降雨與河流荷載效應，相差39hPa，這種差異可能反映靜樂井水位除受降雨與河流荷載作用外，還可能受其他因素的影響。

（4）通過與同一構造區距靜樂井較近的前兆測點寧武鉆孔應變和神池鉆孔應變對比分析可知，靜樂井水位突升前后2 個測點應變觀測值出現大幅度反向變化，表明靜樂井水位高值變化不是獨立事件，可能反映了所在構造區應力場發生了改變。

（5）可反映區域應力場變化特征的幾種參數對2 種方法計算結果產生差異的原因表明，2017 年靜樂井水位高值異常期間山西地震帶3 級地震缺震和b 值、穿過靜樂井的GPS 基線和地震矩釋放均存在顯著的異常，推測靜樂井水位高值異常除受降雨與河流荷載作用的影響外，也可能受構造活動增強的影響。