基于耦合過程的抽水含水層壓縮線性預測方法研究

呂偉華，楊帥，辛文青,孫信，張繼周，薄冠中

(1. 南京林業大學土木工程學院，南京 210037；2. 水能資源利用關鍵技術湖南省重點實驗室，長沙 410014；3. 公路地質災害防治工程技術研究中心，銀川 750004；4. 中交一公局第二工程有限公司，江蘇 蘇州 215011；5. 江蘇交通控股有限公司，南京 210016)

地面沉降已成為我國各大城市的主要工程地質問題之一。其中，上海最大沉降中心已達2.6 m多，沉降超過2 m的城市還有天津、太原、西安等，蘇錫常地區最大沉降中心甚至達到2.8 m。國際上如墨西哥、日本和美國等國家的很多城市也都存在嚴重的地面沉降問題[1]。地下水的過量開采是引發地面沉降問題的最主要因素，而含水層的壓縮變形量不可忽略[2-3]。因此，如何科學評估含水層壓密變形是抽水地面沉降問題研究的重要課題之一。

開采地下水引起地下水位下降，地層的有效應力增加進而發生壓密變形，從而導致地面沉降發生[4-7]。但地質條件的復雜性和差異性、地下水位波動以及地表荷載作用使得要準確預測地面沉降非常困難。依靠抽水地面沉降實測數據統計規律進行預測分析是一種可靠的有效研究方法[8-9]，在基坑工程中也常根據降水實測統計數據資料對含水層的變形特性進行規律分析[10-13]，但這類方法往往需要長期觀測且易忽略土性本質；當采用不同概化模型數值計算分析進行地面沉降的預測[14-16]時，雖然參數選擇容易但合理性不易滿足，結果的可靠性難以令人信服。

本研究采用設計模型箱模擬含水層弱透水層交替地層系統，研究特定抽水控制條件下的含水層壓縮特性，通過水位降深與含水層變形的規律分析，建立基于水力坡降含水層壓密孔隙變化的耦合關系，針對地面沉降中含水層的變形規律進行趨勢預測。

1 模型試驗與結果

模型箱如圖1所示，側面為矩形，尺寸為60 cm×50 cm×100 cm，選用高強有機玻璃板和三角鋼條焊接加固，不考慮側向變形。槽內壁均勻涂抹凡士林作隔水邊界。含水層選用長江底清淤砂，弱透水層選用亞黏土，物理性質見表1，顆粒分析見圖2。如圖1所示，以PVC管在模型中心處模擬抽水井對底部含水層非完整井式抽水，底邊10 cm設進水口，紗布包扎以防堵塞，管壁與土層間隔套管允許自由滑動但不透水；通過槽側壁鉆孔，于含水層內部較低位置處設細水位管并定時開閥測水頭，在各層位交界處分別鋪設沉降板。模型填筑完成后沿抽水管慢速回注61 kg地下水，至填筑體頂面水體浸潤，覆蓋保護層以防水分散失。

圖1 模型槽試驗Fig. 1 Layout of the physical model of aquifer systems

表1 含水層系統填土物理性質參數Table 1 Physical property parameters of tested soil in aquifer system

圖2 砂土、黏土顆粒分析曲線Fig. 2 Particle size distribution of sand and clay

采用虹吸原理從含水層系統埋設的模擬抽水井中以定速/定量按一定間隔時間分級抽水，并將抽出的水密封保存，待整個系統最終沉降穩定后，按原方式全部注水回灌。完整的抽水注水時程曲線如圖3所示，為了讓試驗系統中的土層因釋水發生的緩慢變形得以充分發展，分級抽水的間隔時間較長。含水層的變形-水位降深-歷時曲線如圖4中黑色散點球所示，各坐標平面投影線為相對變化關系曲線。

圖3 抽水、注水時程曲線Fig. 3 Distribution of dewatering and recharging

圖4 抽水含水層變形-水位降深-歷時曲線Fig. 4 Diachronic curve of the pumped aquifer-water level dropdown

很明顯，抽水期結束后含水層的變形仍發生了緩慢增長。參考文獻[16]所用的統計方法，模型試驗所得的水位降深與抽水層的壓縮變形滿足雙線性關系，其斜率可以用來確定抽水含水層的儲水系數。因此，將抽水期的含水層變形進行擬合歸一得到如圖5所示的含水壓縮層累計變形時間半對數曲線。根據Cooper-Jacob半對數直線方法[17]有：

(1)

圖5 壓縮變形-時間半對數擬合曲線Fig. 5 Fitting curves between compression deformation versus time semi logarithmic

式中：s為某時段t所對應的水位降深，m；q為抽水井抽水速率，m3/d；T為含水層的導水系數，m2/d；r為與抽水井的距離，m；Ss為承壓含水層儲水系數，無量綱。

在短時間t1～t2之間水位降深Δs可表示為：

(2)

對于一維固結問題，含水層儲水系數即為：

(3)

式中：Δb為含水層的變形量，m；Δh為抽水層的水位降深，m。抽水含水層的儲水系數可表示為：

(4)

如圖6所示，基于抽水期含水層水位降深與壓縮變形統計后明顯滿足雙線性關系，即達到臨界水位降深時有效應力為含水層砂土層前期固結壓力。據前期固結壓力后的水位降深壓縮變形斜率即可得到抽水含水層儲水系數值。表2所示分別為采用Cooper-Jacob法、雙線法得到的抽水含水層儲水系數與薛禹群[18]、Marsily[19]所給參考值的比較，就系數量級而言，可認定本次模型試驗得到的抽水含水層儲水系數基本合理，模型試驗是有效的。

圖6 水位降深-壓縮變形曲線Fig. 6 Curves between water level dropdown and compression deformation

表2 含水層儲水系數計算對比Table 2 Calculation and comparison of water storage coefficient of pumped aquifers

2 考慮耦合過程的含水層變形

基于大量砂土滲透試驗的數據統計及規律分析，Lambe等[20]給出了初始孔隙比為e0的含水砂層滲透系數與孔隙比之間存在如下經驗關系：

(5)

式中：e為孔隙比，無量綱；k0為初始孔隙比為e0時的滲透系數，cm/s；k為孔隙比為e時的滲透系數，cm/s。

考慮一維流固耦合過程是基于滲流特性與固結條件關系的聯合作用，即可通過含水層儲水系數與砂土的壓縮性指標之間建立定量聯系實現。鑒于含水層的抽水與回灌從力學機制角度理解為對含水層的加載與卸載，其實質相當于在e-lgp曲線中水位降深引起含水層有效應力的增加從而導致抽水含水層砂顆粒間的壓密變形，而壓密后的含水砂層空隙比減少，反映了抽水水位降深的作用影響效果。令cc為超過前期固結壓力后的壓縮指數，建立水頭高度與孔隙比的關系式如下：

cc(lgh-lgh0)=e0-e

(6)

式中：h為對應孔隙比為e時的水位降深，m；h0為初始狀態孔隙比為e0時的水位降深初值，m。根據式(6)可以得到：

e=e0-cclg(h/h0)

(7)

將式(7)代入式(5)并乘以含水層厚度H可得到其導水系數表達式：

T(h)=Hk0[1-cclg(h/h0)/e0]3×

{1+e0/[e0-cclg(h/h0)]}

(8)

為實現對超過含水層前期固結壓力后的水位降深的簡化計算，將公式(8)基于數值計算方法就水位降深與導水系數間的定量關系進行線性等效，可得到如圖7所示的直線表達式：

T(h)=Ah+B

(9)

圖7 導水系數與水位降深關系簡化Fig. 7 The simplified relationship between water conductivity coefficient and water level dropdown

式中，A、B為考慮抽水含水層初始狀態、壓縮性、抽水水位降深關系的常系數，無量綱。

將式(9)代入式(1)，含水層抽水引起的滲流與豎向變形過程耦合，反映了孔隙的減小與滲透性減弱同步發展過程，得到抽水含水層計算的水位降深表達式為：

(10)

根據Newton迭代法，先構造一個迭代函數代替原方程進行求根計算，而后得到新的求解迭代格式，取某一個時刻t0的初始水位h0進行迭代試算，至任一時刻ti的水位降深為hi，直至迭代計算誤差限ε(如可取為0.001)滿足要求即止，滿足收斂性要求。

其中，某計算時刻的水位降深迭代解計算式合理數值解的表達式如下：

(11)

(12)

該方程對抽水含水層壓密變形計算是基于水位降深表達式(11)的迭代，即反映了簡單一維條件下水頭變化的滲流與固結變形的耦合過程。通過模型試驗實測數據分析所得的含水層物理力學和抽水指標參數(表3)，對該計算方法進行合理性驗證。對比結果如圖8所示，明顯可知考慮耦合過程的計算結果較非耦合值與實測結果更加接近。

表3 抽水含水層計算參數Table 3 Calculation parameters of pumped aquifer

圖8 含水層變形-水位降深-歷時曲線對比Fig. 8 Comparison between calculated results and monitored data of aquifer deformation-water level dropdown-diachronic curves

作為蘇錫常地區地下水主采層，第Ⅱ承壓水水位動態完全受開采控制。從開采層次看，區內地下水以第Ⅱ承壓水為主，據統計該層內開采井共2 978 眼，占總井數的70.28%，年開采量約2.2×108m3，占總開采量的65.42%，主要分布在滬寧線兩側，平均單井開采量約為1 000～2 500 m3/d。為進一步驗證上述方法合理性，選取常州清涼小學抽水水位降深、地面沉降歷時監測數據進行計算對比，所選計算參數如表4所示。

表4 原位含水層計算所取參數Table 4 Calculation parameters of pumped aquifer in situ

含水層水位降深和累計壓縮變形如圖9和圖10所示。由圖可見，發生的水位降深與累計壓縮變形的非耦合計算值與耦合計算值都較實測值偏大，但后者更接近實測結果，顯然更合理。與模型試驗相比，清涼小學的水位降深-壓縮變形的雙線性關系此處未見，僅呈現單向趨勢規律發展，原因是該處地下含水層系統已經歷長期抽水水位降深累積作用，即水位降深值早已超過前期固結壓力。綜上分析，本研究方法可為后續提出一種含水層壓縮量評估的預測方法提供基礎。

圖9 含水層水位降深歷時曲線計算對比Fig. 9 Comparison between water level drawdown-diachronic curve of pumped aquifer

圖10 含水層累計壓縮變形歷時曲線計算對比Fig. 10 Comparison between cumulative compression deformation-diachronic curve of pumped aquifer

3 地面沉降預測方法

抽水含水層的水位降深本質上可等效為對砂顆粒骨架的加載，含水砂層顆粒間孔隙被擠壓而體積縮小，整個土層發生壓縮變形，其滲透水能力將減弱。為考察水位降深對含水砂層滲透系數影響敏感性，設計一組理想計算模型試驗。基于公式(10)，令含水層中某一抽水井的抽水量q=45 m3/d，含水層平均厚度H=100 m，其儲水系數設為常數Ss=10-5，初始孔隙比為e0=0.756，當抽水時間為365 d時，滲透系數從8.8×10-4cm/s變化為1.2×10-5cm/s，導水系數從76 m2/d變化為1.0 m2/d，發生水位降深37.25 m。

距抽水井r=50 m處，水位降深、滲透系數隨抽水時間變化關系如圖11所示。以某一速率抽水，水位降深隨時間不斷增加，而降深較小的初段滲透系數減小顯著，當超過某水位降深后，其變化不是很明顯。圖12所示為理想抽水試驗模型的含水層壓縮變形、水位降深受抽水量影響的關系曲線。

圖11 理想計算模型變化曲線Fig. 11 Change curve of ideal calculation model and results

抽水是引起地面沉降的根本原因，因此基于沉降隨抽水量的發展變化規律，通過如圖13所示控制抽水速率和抽水時間從而限制總抽水量以達到控制地面沉降的目的。

圖13 抽水地面沉降預測方法Fig. 13 Prediction method of land subsidence

基于上述研究對常州清涼小學第Ⅲ承壓含水層進行實例分析，所選計算參數同上。基于實測值預測、非耦合計算預測與耦合計算預測結果3部分組成的對比情況見圖14，實線部分為已知的統計數據，虛線部分為預測。如果仍按2000年結束時的抽水速率，而不采取控制開采地下水量的措施，則至2010年左右，第Ⅲ承壓含水層的實際水位將降深至70 m，考慮耦合的計算預測為70.8 m，而不考慮耦合的計算預測水位則達到73.5 m；同理，至2010年，該含水層發生的累計壓縮變形根據實測數據預測為163.2 mm，根據耦合計算值預計將達到191 mm，而非耦合計算結果預測為223.9 mm。

圖14 實測數據與耦合、非耦合計算預測結果對比Fig. 14 Comparison between measured results and coupled and uncoupled calculation results

根據多年大量實測數據所做的預測理應更接近真實值，而考慮耦合過程得到的預測值比不考慮耦合情況下所得預測值顯然更合理。因此，在長期抽取利用地下水地區，對于缺乏足夠實測水位降深以及地面沉降資料的情況，可以借鑒采用本研究耦合計算結果的直線法預測控制沉降量，從而反推控制抽水量。

4 結 論

通過室內模型試驗研究了對含水層系統抽水水位降深、含水層壓密變形的時間變化規律，并進一步提出了抽水含水層孔隙變化影響滲透性的水位降深迭代數值解析方法，對抽水含水層壓縮量的預測方法進行了初步探討，得出以下結論：

1)模型試驗中抽水含水層的儲水系數數量級達到10-2，滿足含水層合理性賦存條件，抽水含水層的水位降深與壓縮變形之間滿足雙線性關系且確定轉折處為其前期固結壓力。

2)提出了一種通過Cooper-Jacob解代入滲透系數隨水位降深變化的迭代數值解法，并據此推求了一種簡化的流固耦合變形計算方法，通過模型試驗結果與現場實測數據對該方法的合理性進行了驗證。

3)含水層水位降深、累積壓縮變形隨時間變化存在明顯線性發展趨勢，提出了一種直線法預測含水層壓密變形方法，且將考慮耦合過程引入該線性預測算法的結果比非耦合算法更接近實測數據預測值。

4)對于缺乏足夠實測水位降深及地面沉降數據的情況，可以采用本研究所推薦的直線法預測壓縮量控制抽水量，這對防治由抽水地面沉降導致的地質災害評價具有參考意義。