基于井水位潮汐響應探討含水層的水力特征及其變化

趙可新，蘭雙雙，谷洪彪，梁文宇，喬鵬

（1.北京工業大學，北京 100124；2.防災科技學院，河北 三河 065201）

含水層的水力特征參數是表征地下水運動規律的重要指標，包括貯水系數、導水系數、滲透系數、越流系數等。當含水層所受應力狀態發生變化時，可能會引起其水力特征發生變化，進一步導致地下水流發生相應的變化［1-2］，從而威脅水庫壩基的穩定性以及地下水的安全性，甚至可能對深部地熱、石油、天然氣等重要地下資源的運移和儲存［3-4］產生一定的影響，故研究應力荷載作用下含水層水力特征的變化具有非常重要的實際工程意義。

含水層水力特征參數的傳統獲取方法主要有抽水試驗和注水試驗，但這些方法無法推求含水層滲透系數的實時連續變化，因此無法用來對外加荷載作用下含水層滲透性的變化進行分析，而含水層水頭對周期高頻性自然力（氣壓、固體潮等）的響應模型為該方面的研究提供了新的思路［5-6］，井水位的氣壓效應可以用來推求含水層的骨架材料特性［7-8］，井水位的潮汐響應可以用來計算封閉性承壓含水層的水平滲透性［9-11］。最新研究表明地震應力不僅可以改變地下水在井孔-含水層之間的水平流動狀態［12-13］，也可能通過裂隙的張開或閉合引起相鄰含水層垂向上的水力聯系［14］，Wang等［14］引入越流含水層模型，并通過越流系數的解析解分析了美國Oklahoma深井的滲漏現象，由此地震引起含水層垂向水力特征參數的改變引起學者們的高度重視。人們［15］通過川滇陜甘渝地區井水位的潮汐響應發現，大部分含水層并不是完全承壓的，在局部地區地下水存在著垂向運動，而以往的成果中對該方面的定量計算及討論涉及較少，地震引起含水層垂向水力特征參數改變的可參考研究成果也相對較少。

綜上，為了更加客觀地、定量地反映地震引起含水層水平及垂向水力特征的改變，本文在系統分析總結井水位潮汐響應與含水層水力特征理論關系的基礎上，以華鎣山斷裂帶附近大足井和北碚井為例，判斷井孔附近含水層系統中地下水流運動狀態，計算含水層水力特征參數，并以汶川地震為震例，討論地震應力對含水層水平及垂直水力特征及地下水流運動狀態的改變，探究地震引起兩口井水位的變化原因。該研究成果可拓展傳統地下水動力學的研究領域，進一步提高對外部荷載應力驅動含水層中地下水滲流運動的認識，為井孔-含水層系統對地震活動響應機制的研究提供參考依據。

1 理論基礎

在日、月引潮力作用下地球會發生周期性的壓縮與膨脹，造成含水層內部孔壓隨之增大與減小，引起觀測井與含水層之間產生壓力差并發生水流交換，最終導致井水位發生周期性波動，這種現象稱為井水位的固體潮效應。具體可分解為兩個過程：一是含水層內部水頭對潮汐體應變和潮汐垂向流的響應，前者是不排水孔彈性響應過程，后者是含水層垂向上的不均勻導致含水層水頭重新分布的過程，該過程只在垂向流存在的含水層中發生；二是井水位對含水層壓力水頭的響應，它主要受潮汐徑向流的影響，是一個典型的水動力學過程，在垂向流和徑向流中均存在此過程，井水位對壓力水頭的響應主要取決于含水層水動力學參數（如導水系數、貯水系數等）［16-17］?？梢?，井水位對固體潮的響應實質上反映了含水層的水力特征及地下水流運動狀態，故通過分析井水位的潮汐特征可以達到推求含水層水力參數、判定地下水流運動特征的目的，從而為解釋地震應力對井孔-含水層系統的影響提供依據。

1.1 計算貯水系數

貯水系數S和貯水率Ss是描述含水層釋水能力的參數，貯水率乘以含水層厚度等于貯水系數。

1.1.1 地球對潮汐的響應

在均質、彈性條件下，地球的潮汐水平應變［18］可以表示為

式中：εhh為潮汐水平應變；hn和l n為勒夫數；Wn為引潮位；a為地球半徑；n為含水層孔隙度。

垂向應變一般難以計算，利用泊松比v建立水平與垂向應變之間的關系，則體應變［19］表達式為

式中：εkk為潮汐體應變；εrr為潮汐垂向應變。

對于二階引潮位的變化（ΔW2），利用式（2）得到潮汐體應變［20］為

1.1.2 地下水對固體潮的響應

承壓含水層水位升降會引起含水層兩部分體積變化［21］：一部分是水位上升后井中水體的體積變化量dV1＝Ah；另一部分是水位上升導致的靜水壓力增加使得含水層內水體體積變化，即dV2＝－gρwh·（Vw／Kw）＝－gnρwh（V／Kw），則固體潮體應變引起含水層體積的變化可表示為

式中：ρw為水的密度，kg／m3；Kw為水的體積模量，Pa；V為含水層體積，m3。因為KwA＋gρwV，則承壓含水層水位變化為

假設含水層體積變化等于含水層孔隙度變化，由貯水率的定義可知含水層的貯水率［22］為

將式（6）代入式（5）可得井水位和貯水率的關系為

將式（3）代入上式可得到引潮位表示的貯水率［20］

式中：泊松比v＝0.25，h2＝0.606，l2＝0.084，地球半徑a＝6 371.393 km；水頭變化Δh是引潮位ΔW2引起的，W2＝g Kmbf（θ），Km是與地球質量、月球質量、地球與月球之間的距離、地球半徑有關的常數，約等于0.537 m［23］；b是與潮汐波分量周期相關的常數，對于M2波，b＝0.908；f（θ）是緯度θ的函數，f（θ）＝0.5cos2θ。

1.2 計算導水系數

在引潮力的作用下，含水層巖石產生形變，壓力水頭改變，驅動井與井周含水層之間水流交換，地下水運動以水平流為主，簡稱潮汐徑向流。值得注意的是，引潮力作用于含水層到井-含水層之間發生水流交換，存在著時間滯后，即相位差。Paul等［24］提出潮汐徑向流條件下壓力水頭振幅比A和相位差ηr與含水層參數的解析解為

其中

式中：ω為潮汐分波頻率，d－1；τ＝2π／ω是潮汐分波周期；Kei和Ker為零階Kelvin的函數；T為徑向導水系數，m2／d；rω為井濾水管半徑，m；rc為井筒（水位變動段）半徑，m；H為水井水位的變化量，m。

由式（9）、（10）可知，對于潮汐徑向流而言，振幅比A和相位差ηr是關于導水系數T和貯水系數S的函數，圖1顯示了相位差ηr和無量綱參數之間的關系，對于徑向流而言，相位差ηr＜0；當在一個可能的值域范圍內，相位差隨T或的增加而增大，故根據ηr的取值可估算參數，進而確定含水層的T值。

圖1 不同條件下，徑向流相位差（ηr）隨的變化［9］Fig.1 Under different conditions，variation of phase shift in radial flow with

1.3 計算越流系數

自然界中存在的井-含水層系統上、下巖層往往并不是絕對隔水的，其中一個或兩個可能是弱透水層：當這個含水層和相鄰含水層之間存在水頭差時，地下水便會發生越流運動；或當觀測井本身為非完整井時，觀測含水層與下伏含水層之間或許存在著水量交換。由此，Wang等［14］推導出越流含水層系統中潮汐力驅動地下水在含水層和井孔間運動的控制方程，并得出了振幅比和相位差的解析解

其中

在式（14）、（15）中，當K'＝0即沒有越流存在時，該方程與Hsieh的模型解析解一致，如圖2所示相位差是關于越流系數σ′＝K′／b′、貯水系數S以及導水系數T的函數：如果貯水系數S和導水系數T已知，當相位差為負值時，此時垂向上的水力聯系較少，含水層以水平徑向流為主；隨著相位差增大為正值，此時σ′的值變大，含水層垂向滲漏增強。因此，在S和T已知的含水層系統中，由相位差可以得到σ′的值。

圖2 越流含水層系統中相位差與越流系數的關系Fig.2 Relationship between the coefficient of leakage and phase shift in leaky aquifer system

2 觀測背景與數據

華鎣山斷裂帶位于中國西南地區，是一條右旋走滑型逆斷裂帶。該斷裂帶北起萬源，向南西經達川、榮昌至宜賓西南，全長約460 km［25］，是川東隔擋式褶皺帶和川中平緩構造的分界斷層。地震流體觀測井大足和北碚皆位于華鎣山斷裂帶附近10 km范圍內［圖3（a）］，隸屬四川盆地的重慶市，鉆孔揭露的地層主要為侏羅系中統沙溪廟組紫紅色泥巖和淺灰白色長石砂巖互層，即“紅層”。其中，長石砂巖為該區地下水賦存的主要含水層，單層厚度約10～50 m。含水砂巖上下均被相對隔水的泥巖所夾持，具有層間承壓性，且該類型含水層廣泛分布，在較大范圍內存在著水力聯系。多孔抽水試驗觀測數據［26］顯示：位于同一含水砂巖層中相距數百米或千米之外的鉆孔，一井抽水會引起另一井水位發生顯著下降。該區域紅層承壓水［27］主要在露頭區接受大氣降水的補給，其次是地表水的垂直入滲和部分越流補給，降水和地表水通過含水層暴露于地表發育的裂隙系統下滲。在重力作用下順層或沿砂巖裂隙發生徑流或垂向越流補給，當含水層被切割時，地下水以泉的形式排泄于地表或地表水體，也可以通過上下相鄰含水層進行越流排泄。

圖3 觀測井位置與柱狀Fig.3 Location of the observation wells and stratigraphy of the borehole

大足井位于華鎣山斷裂帶中部，西山背斜北西翼［16］。如圖3（b）所示，該井深108.7 m，46 m下設套管，直徑127 mm，46 m以下為裸孔，井徑150～120 mm。其中：地表出露地層主要是中侏羅統沙溪廟組泥巖夾砂巖；5.50～47.00 m為青灰色長石砂巖；層中32～35 m、38 m深度段集中發育傾角為30°～75°的裂隙；46.0～108.7 m巖性以泥巖為主。地下水類型屬孔隙裂隙承壓水，為非自流井。

北碚井位于華鎣山斷裂帶東側，觀音峽背斜東翼。如圖3（c）所示，井深105.36 m，Φ127的套管下至42.1 m處，套管外用水泥固井止水。井孔巖性自地表往下：0～3.54 m為覆土；3.54～29.90 m為紫色泥巖，巖層傾角35°，層面風化裂隙發育；29.90～70.24 m為主要觀測含水層，井孔巖性為中侏羅統沙溪廟組淺灰白色中-粗粒厚層狀長石砂巖，該段裂隙較發育，裂隙傾角55°～80°，且30.40、43.30 m處的裂隙被方解石半充填，38.40 m處的裂隙被黃鐵礦晶粒充填；70.24～105.36 m為紫紅色砂質泥巖，夾紫灰色粉砂巖條帶［16］。地下水類型屬裂隙承壓水，為非自流井。

北碚井和大足井水位觀測方式均為靜水位觀測，觀測儀器皆為LN-3A型水位儀，采樣頻率為每分鐘1次，觀測數據來自于中國地震前兆監測臺網。如圖4所示，北碚井和大足井年變化趨勢基本一致，受水文因素影響明顯，豐水期水位較高枯水期水位較低；井水位日動態變化穩定，正常變化幅度約6 cm，具有清晰的兩峰兩谷的潮汐變化形態。且北碚井和大足井水位的頻譜分析結果（圖5）表明：對兩口井水位影響較大的潮汐成分是半日波M2和S2，周期分別為0.517 5 d和0.500 0 d，故可以利用井水位的潮汐效應探討含水層的水力特征變化。

圖4 大足井和北碚井水位的年變化曲線Fig.4 Annual variation of the water level in wells Dazu and Beibei

圖5 大足井和北碚井水位頻譜分析結果曲線Fig.5 Results of spectrum analysis of the water level in wells Dazu and Beibei

3 含水層水力特征參數的計算

根據中國地震臺網可知2012年1月1日至5月1日川滇地區無Ms 5.0以上地震發生，且該時段屬于枯水期，井水位受降雨因素干擾少，故本次選取該時間段內大足井和北碚井水位小時值觀測數據進行潮汐分析，進而推求含水層的水力特征參數。潮汐波是由一系列不同周期的潮汐分量組成，其中M2波（半日波）水位數據振幅較大，受氣壓等因素干擾小，誤差相對較小，因此本次借助國際通用潮汐分析軟件Baytap-G程序進行井水位對M2波的潮汐分析。數據間隔設置為30 d，滑動窗長設置為15 d，分別提取兩口井水位對M2波的響應特征參數（A和η），計算結果見表1。

表1 大足井和北碚井對M2波的潮汐分析結果Tab.1 Results of tidal response of M2 wave in wells Dazu and Beibei

據表1可知大足井水位對M2波響應的相位差η＜0，表明該井附近含水層地下水流類型以水平徑向流為主。大足井管揭穿了整個含水層厚度，含水層分布穩定，地層傾角約為20°，地下水流運動方向受巖層傾角影響，接近水平流。利用大足井M2波的理論引潮高振幅和該井觀測段水頭響應振幅數據，據公式（8）推求該井的觀測含水段的貯水能力，結果見表2，大足井觀測含水層的貯水率為1.55×10－6m－1；由該井鉆孔資料將觀測含水層厚度取值為41.5 m，則含水層的貯水系數為6.45×10－5。由此利用公式（9）和（10）及相位差可計算觀測井附近含水層的導水系數見表3，導水系數平均值為1.74 m2／d，滲透系數（表3）為0.041 9 m／d，與前人［28］抽水試驗的結果0.053 4 m／d一致，表明利用井水位的潮汐響應計算含水層水力參數具有可行性。

表2 大足井和北碚井所處含水層貯水能力計算Tab.2 Calculation results of the storage capacity in the aquifer near wells Dazu and Beibei

表3 大足井所處含水層滲透系數Tab.3 The coefficient of permeability in the aquifer near well Dazu

北碚井的潮汐相位差η＞0（表1），說明北碚井附近含水層系統以垂向流為主、徑向流為輔。根據前人鉆孔資料，該井鉆孔揭露本區29.9～70.24 m的長石砂巖含水層，井深30.40 m和36.32 m處涌水，水頭高度分別為1.77 m和2.74 m，表明該含水層中地下水局部存在著自下而上的垂向流運動。原因為：一是據鉆孔資料知井觀測含水段裂隙較發育，裂隙傾角約55°～80°，判斷局部發育的裂隙為含水層內地下水垂向運動提供了優勢通道；二是井深30.40、43.30 m處的裂隙被方解石半充填，38.40 m處的裂隙被黃鐵礦晶粒充填，形成局部相對弱透水層與含水層互層，觀測段（42.1～70.24 m）可通過弱透水層與上覆含水層發生垂向水量交換。由于Hsieh模型適用于封閉性良好的承壓含水層徑向流（η＜0），Wang模型適用于存在垂向水量交換的越流含水層系統，故本次采用Wang模型計算北碚井附近含水層的水力特征參數，用以表征垂向上水量交換的能力。首先根據北碚井鉆孔資料將觀測含水層厚度M取值為40.34 m，利用潮汐M2波引潮高振幅和該井觀測段的井水位數據，據公式（8）估算北碚井含水段的貯水能力，結果見表2，求得觀測含水層的貯水率為1.41×10－6m－1，貯水系數為5.68×10－5。通過前人抽水試驗［28］可知該井觀測段含水層的水平滲透系數K＝0.014 9 m／d，則導水系數T＝KM＝0.60 m2／d，依據方程（14）和（15）可繪制出北碚井的越流系數和相位差的關系曲線（圖6）。由表4可知北碚井水位潮汐分析的相位差變化范圍約為25°～35°，則可估算越流系數的范圍為1.5×10－8～2.5×10－8s－1。

表4 北碚井所處含水層系統越流系數計算Tab.4 Calculation results of the coefficient of leakage in the aquifer near well Beibei

圖6 北碚井相位差與越流系數關系Fig.6 Relationship between the coefficient of leakage and phase shift in well Beibei

4 地震對含水層水力特征的影響

4.1 汶川地震引起的井水位變化

2008年5月12日14時28分，我國四川省汶川縣映秀鎮與臥龍鎮之間發生了Ms 8.0級地震，地震發生在青藏高原東緣的松潘-甘孜地塊與揚子地塊交界的龍門山主中央斷裂帶上，該地震造成了我國大部分井水位發生同震響應［29］，其中包括大足和北碚井。如圖7（a）所示，汶川地震發生時，大足井水位發生94 cm的階降變化，2天后呈現正常潮汐波形，但數月內未恢復到震前水平。北碚井水位出現明顯階降現象，幅度為93.6 cm，3個小時后水位重新達到新的平穩狀態，但數月內未恢復到震前狀態，見圖7（b）。

圖7 汶川地震引起井水位變化Fig.7 Variation of well water level caused by the Wenchuan earthquake

4.2 汶川地震引起含水層水力特征的變化

汶川地震發生時大足和北碚井水位皆出現了明顯的同震響應，利用井水位的潮汐效應定量分析地震前后含水層水力特征是否發生了變化。以大足井和北碚井2008年的水位小時值觀測數據為基礎，進行固體潮調和分析，為了準確地描述汶川地震前后含水層水力特征的變化，此次數據間隔設置為20 d，滑動窗長設置為10 d，計算地震前后大足和北碚井固體潮M2波的相位差。

如圖8（a）所示，汶川地震前大足井的相位差在0°上下浮動，地震發生后相位差最高達到11°，根據相位差和地下水流動類型的關系判斷，震后地下水發生了垂向流運動，表明含水層在垂向上存在著水量交換，由此引入越流系統模型來估計大足井附近含水層的垂向滲透能力，基于水平滲透系數以及含水層的貯水能力，根據公式（14）和（15），得到井水位固體潮響應相位差與含水層越流系數之間的關系（圖9）以及越流系數的計算結果（圖10），可見汶川地震發生前后大足井的水位相位差上升為正值，相應的越流系數增大，即垂向滲透性增強。

圖8 汶川地震前后相位差變化Fig.8 Variation of phase shift before and after the Wenchuan earthquake

圖9 大足井相位差與越流系數關系Fig.9 Relationship between the coefficient of leakage and phase shift in well Dazu

如圖8（b）所示，汶川地震的發生造成北碚井的相位差由震前平均值26°上升到67°，上升幅度為41°，變化量超過3倍均方差，表明含水層內部地下水運動仍以垂向流為主，且垂向流運動增強。由此引入越流含水層系統模型推求北碚井附近含水層越流系數的變化，以表征含水層垂向水量交換的能力。計算結果見圖10，汶川地震的發生導致含水層越流系數由震前平均值1.692×10－8s－1上升到2.541×10－6s－1，根據越流系數的定義可知，觀測含水層與相鄰含水層之間垂向水量交換增大?？梢?，北碚井和大足井的潮汐參數在汶川地震后均有所增大且大于0，表明兩口井附近含水層垂向流運動加強。引入越流系統模型推求含水層水力特征參數的變化，結果顯示兩口井所在含水層越流系數皆增大，定量說明了觀測含水層與相鄰含水層的垂向水流交換增強。

圖10 汶川地震前后含水層系統越流系數變化Fig.10 Variation of the coefficient of leakage in the aquifer system before and after the Wenchuan earthquake

4.3 含水層水力參數變化機制探討

地震的發生會造成區域構造應力（靜應力）的變化和地震波（動應力）的傳播，導致含水層形變甚至引起含水層水力特征（滲透性）的改變。對于靜應力引起的井水位變化的解釋，孔隙線彈性理論［30］和松散含水層固結理論［31］比較成熟。根據Biot［30］孔隙線彈性理論，震源破裂產生的同震靜應力引起地殼發生變形，導致含水層內部孔壓發生變化，故表現為壓縮區域井水位上升，膨脹區域井水位下降。固結理論認為，對于松散沉積物，剪應變引起顆粒運動進入到原有的孔隙中，導致孔隙體積減小，孔壓增大，井水位發生變化。對于動應力對含水層造成的影響，一般認為地震波作用于井-含水層系統，導致裂隙的疏通［32］（裂隙中膠體顆粒、液滴和氣泡的活化遷移）或堵塞，裂隙的張開或閉合［33］，甚至是新裂隙的形成［34］，引起含水層滲透性或層間水力聯系的改變，從而造成含水層內部孔隙壓力的重新調整與平衡，引起井水位的響應。對于地震引起井水位以及含水層水力特征變化的原因，可能并不是一種機制可以解釋的，還需要結合井孔周邊的水文地質條件以及井自身的狀況進行具體分析。

汶川地震后，北碚井和大足井的水位發生了明顯的階降變化且地下水垂向流增強。一方面，據Shi等［35］得到汶川地震同震破裂體應變空間分布與井水位同震響應形態對比（圖11），可知北碚井和大足井位于含水層膨脹區域，即汶川地震發生后，震源破裂產生的構造應力使北碚和大足井附近含水層發生了膨脹變形，孔壓降低，井水位明顯下降。另一方面，地震發生時，應力傳播至含水層，加快了地下水流動速度，可能致使堵塞在含水層局部裂隙中的鈣質、鐵質膠體或細小顆粒被清除，從而造成大足井和北碚井附近含水層中裂隙疏通，連通性增大，且因為裂隙傾角近垂直，所以大足井表現出垂向流的特征，北碚井垂向交換量增強。值得注意的是，該區域主要觀測含水層在較大范圍內存在著統一的水力聯系［26］，地震應力使得裂隙通道張開，各個井孔附近含水層間的水量交換加強，勢必造成在巖層露頭部位孔隙壓力的聚集和釋放，井水位迅速上升。如巴南井為該區一口自流井，兩口井距離約80 km，該區地勢較低，地下水沿砂巖裂隙通道出露地表，如圖12所示，汶川地震的發生導致該井水位發生幅度約為120 cm的階升現象，這也是區域地下水徑流運動狀態改變的結果。綜上，該區域上沙溪廟組砂巖含水層分布面積廣，且具有統一的水力聯系，大足井和北碚井皆地處徑流區，汶川地震的發生導致含水層發生膨脹變形和裂隙疏通，地震造成井水位下降及垂向水量交換增強，井水位變化形態基本一致，巴南井地勢較低且為自流井，孔隙壓力的聚集和釋放導致井水位上升。

圖11 汶川地震同震孔隙體應變［28］（膨脹為正）與井水位變化Fig.11 Co-seismic pore volume strain and well water level changes in the Wenchuan earthquake

圖12 巴南井水位變化曲線Fig.12 Variation of well water level of well Banan

5 結 論

本文系統分析了井水位潮汐響應與含水層水力特征的理論關系，以華鎣山斷裂帶上大足井和北碚井為例，基于井水位的潮汐分析計算了觀測段含水層水力特征參數（包括貯水系數、滲透系數、導水系數、越流系數），并探討了汶川地震對含水層水力特征及地下水流運動狀態的影響及機制。主要結論如下：

由潮汐分析知大足井附近含水層中地下水運動以水平徑向流為主，北碚井孔-含水層系統中裂隙較發育，地下水垂向流和徑向流并存；利用潮汐徑向流模型計算出大足井含水層滲透系數為0.041 9 m／d，與以往抽水試驗的結果（0.053 4 m／d）基本一致，表明利用井水位的潮汐響應計算含水層水力參數具有可行性，引入越流含水層系統模型計算得到北碚井越流系數范圍為1.5×10－8～2.5×10－8s－1。

汶川地震引起北碚井和大足井的水位均發生了階降變化，且含水層中地下水垂向流運動加強，原因是該區域侏羅系砂巖含水層在較大范圍內存在著統一的水力聯系，且局部裂隙發育，傾角較大，汶川地震的發生導致含水層發生膨脹變形和裂隙疏通，造成了徑流區北碚井和大足井水位下降及垂向水量交換增強，排泄區井水位上升。