華北平原地下水位驅動下的地面沉降現狀與研究展望

郭海朋，李文鵬，王麗亞，陳 曄，臧西勝，王云龍，朱菊艷，卞躍躍

（1.京津冀平原地下水與地面沉降野外科學觀測研究站，北京 100081；2.中國地質環境監測院，北京100081；3.北京市水文地質工程地質大隊，北京 100195；4.中國地質大學（北京）水資源與環境學院，北京 100083；5.中國地質博物館，北京 100034）

地面沉降是世界各地普遍發生的環境地質問題，而地下水動力與土層變形的互饋關系研究是地面沉降區地下水資源量評價的前提，是區域地面沉降預測預警的基礎，也是環境地質學科發展的需求。目前國際上關于地下水動力與土層變形互饋機制的研究，主要集中在地下水位驅動下黏性土和砂土變形規律研究，時序沉降數據用于計算區域水文地質參數和校正模型參數，以及土-水耦合模型研究三個方面。華北平原是我國地面沉降影響面積最大的區域，地面沉降對京津冀協同發展、雄安新區、北京城市副中心等國家重大戰略區規劃建設和南水北調、高速鐵路等國家重大工程區安全運營存在重要影響。地下水位變化與土層變形互饋關系研究，對華北平原地面沉降監測網的優化、監測數據的定量化分析、地面沉降模型仿真性和預測預警準確性的提高等方面具有科學意義，可為有效緩解該地區地下水利用引起的地面沉降災害提供技術支撐。

1 國內外研究現狀

1.1 地下水位驅動下黏土和砂土變形規律研究

近幾十年來，含水層系統中地下水位變化影響下的土層變形規律受到國內外許多學者的高度關注。例如，He等[1]實施了物理模型試驗，研究地下水補給和排泄引起的砂層和黏土層變形特征及孔隙水壓力變化規律，提出了三種有利于減緩地面沉降發展的地下水開發模式。郭海朋等[2]初步總結了天津平原地下水水位變化模式，分析了不同地下水位變化模式下的土層變形特征。總體而言，該領域研究成果主要集中在土層尤其是黏性土層的流變變形特征研究和土層變形應力應變關系定量化研究等方面。

1.1.1 地下水位驅動下土層變形流變特征研究

隨著地面沉降問題研究的深入，國內外越來越多研究者逐漸認識到抽水誘發下土層的壓縮量中包含了一定比例的流變變形量，不僅黏性土可能發生蠕變[3-5]，砂礫含水層也可能發生蠕變[6-7]。過去在研究地面沉降時，通常將含水砂土層作為線彈性體處理，或者認為砂土層的變形對地面總沉降的貢獻可以忽略不計[8]。試驗結果表明在單調荷載下，砂土的蠕變變形可以達到總變形的10%[9]。分層標監測數據也表明，含水砂層的變形不僅包括彈性變形，而且包括塑性變形和蠕變變形[2,10]。研究人員根據長三角地區大量地面沉降和水位觀測資料，結合室內土工試驗，揭示了砂土層在不同的應力條件下，有的表現為彈性變形，有的則表現為以塑性變形為主并包含有蠕變的非線性變形特征[10-12]。李玉岐等[13]利用自行設計的砂土排灌水試驗裝置，分析地下水排灌引起的砂樣宏觀豎向變形及細觀移動，結果表明：在砂樣排灌水的初期階段，砂樣的結構發生了明顯重組，砂樣沉降不僅發生在排水時，而且在回灌時繼續增加，砂樣產生了較大的、不可恢復的塑性和黏性變形；砂樣結構達到相對穩定后，排水時產生的豎向變形變小，而回灌時砂樣發生較大的回彈。郭海朋等[2]通過土體高壓固結試驗表明冀中坳陷內武清凹陷地區及滄縣隆起地區中晚更新世地層黏性土體在不同荷載條件下，蠕變特征明顯。

1.1.2 地下水位驅動下土層變形本構關系研究

Terzaghi固結理論是簡單的彈性變形模型，由于機理簡單和控制方程的線性特征，被廣泛應用于土層壓縮分析。考慮到過于復雜的黏彈塑性應力-應變模型很難直接用于區域地面沉降模型中，Corapcioglu等[14]基于劃分主、次固結的黏土固結理論，提出了一維線性黏彈本構關系（Merchant模型）。吳林高等[15]在分析抽灌水引起地面沉降的力學機制基礎上，得出了廣義Kelvin本構模型。以上兩個模型的適用性都缺乏試驗數據驗證。Ye等[16]通過分析長三角地區土層變形特征，采用修正的Merchant模型刻畫土層的變形規律，在此基礎上建立了三維變參數水流和垂向一維沉降組成的部分耦合模型，對長三角區域的地面沉降進行了模擬。如果能獲取土層變形的野外監測和室內試驗數據，可能會建立更加貼近實際的土層變形本構模型并用于地面沉降模擬和預測[17]。Tsai等[18]研究提出了一種黏彈塑性模型，該模型由一個黏彈性模型和一個黏塑性模型組成，并由單刀雙擲開關控制黏塑性模型的使用，經驗證，該模型計算結果具有較小的相對誤差。

1.2 時序沉降數據用于計算區域水文地質參數和校正模型參數

利用分層標時序土層變形和地下水位數據反演計算區域水文地質參數已經為大家所熟悉，應力-應變圖解法是常用的手段之一。該方法最早由Riley[19]提出，根據應力應變關系曲線可求出彈性貯水率和非彈性貯水率，Cleveland等[20]提出的圖解法可進一步確定垂向滲透系數。葉淑君等[21]結合上海沉降和水位觀測資料，求出了上海地區弱透水層的彈性貯水率、非彈性貯水率和垂向滲透系數，為上海地面沉降的數值模擬提供了較好的初始參數值。Zhang等[22]利用北京平原某分層標監測數據計算了彈性貯水率和非彈性貯水率，并與美國加利福尼亞州圣荷昆谷（San Joaquin Valley）的對應參數值進行了相互比較和評價。

區域地面沉降監測技術在過去的幾十年經歷了較快的發展，其主要監測手段從經典的水準測量發展至日趨成熟的GPS技術與合成孔徑雷達干涉測量技術（Interferometric Synthetic Aperture Radar，InSAR），獲取了大量的區域地表形變數據[23-25]。水文地質參數是地面沉降數值模擬過程中必不可少的信息，但由于含水層的非均質等復雜性質，許多水文地質參數的空間分布在現有技術條件下無法快速準確獲取，因而利用一些易獲取且精度較好的監測數據（如地下水位等）推估難以獲取的水文地質參數分布成為了地下水領域的研究熱點[26]。當含水層系統相對簡單，不存在地下水位變化和地面沉降之間明顯滯后的情況下，區域地面沉降和地下水位監測數據可以用來直接反演計算彈性貯水率、非彈性貯水率等參數，進而可以補齊時間或空間上地下水位觀測數據的缺失[27-29]。Zhuang等[30]探索了利用地質統計方法反演求取弱透水層滲透系數、彈性和非彈性貯水率的方法，并成功用于常州某地弱透水層水文地質參數的計算。Zhuang等[31]基于歐拉公式，獲得了地下水位變化驅動下超固結弱透水層彈性變形的解析解，并用來評估上海某35.54 m厚弱透水層的垂向滲透系數和彈性儲（釋）水率，結果表明該解析方法可以定量解釋地面沉降滯后特征并有效評估超固結弱透水層的水力參數。

高分辨率水文地質參數是地面沉降模型實現準確預測的基礎，而InSAR、水準、GPS等多種手段獲取的高時空分辨率地面沉降監測數據使有效校正模型的水文地質參數成為可能。UCODE[32]是一種基于梯度的自動反演算法，研究人員已經開展了利用InSAR解譯計算的地面沉降數據和地下水位監測數據并基于UCODE進行模型水文地質參數校正的研究[27,33-34]。Zhang等[34]基于UCODE并利用InSAR獲取的高時空分辨率地面沉降數據和地下水位數據校正了含水層導水系數、彈性和非彈性貯水率、斷層導水系數等水文地質參數，結果表明聯合利用地面沉降和地下水位監測數據校正的模型比僅用地下水位數據校正的模型更符合實際。但是，基于梯度的反演方法主要缺點是計算量大，不能快捷地評估參數，而且預測具有不確定性。

1.3 土-水耦合模型研究

抽排水引起的地面沉降問題實際上就是一個滲流場與應力場耦合的問題，土-水耦合模型按滲流場與應力場結合方式分為不耦合的兩步計算模型、部分耦合模型和完全耦合模型。不耦合的兩步計算模型分為完全獨立的兩步，先計算孔隙水壓力，再計算變形，水流及沉降方程中各參數不隨時間變化。兩步計算模型最初在研究威尼斯的由多層含水層與弱透水層組成的地下水系統因抽水引起的地面沉降問題時提出[35]。Shearer[8]以太沙基一維固結理論為基礎，在Modflow程序的基礎上開發了考慮黏性土層孔隙水壓力變化滯后于含水層水位變化的夾層排水模塊（IDP），通過沉降計算獲取弱透水層的釋水，實現了大厚度含水層釋水壓縮的模擬。部分耦合模型中孔隙水壓力和變形分步計算，但兩者之間相互影響。含水層水頭下降會導致相鄰弱透水層中的地下水產生滲流，弱透水層的變形具有明顯的非線性特征；同時，土層的孔隙比、壓縮性和透水性等參數也隨著土體變形量的變化而變化。國內外地面沉降模型大多基于土-水部分耦合模型[36-39]。完全耦合模型基于Biot固結理論，孔隙水壓力和土層變形同時計算，在理論上這種模型最符合沉降物理機制，但是物理場控制方程更為復雜、計算量大、需要參數多，很難在實際問題中應用。近年來研究人員在完全耦合模型的數值計算方法方面進行了探索[40-41]，并將完全耦合模型成功應用于意大利威尼斯、美國拉斯維加斯和中國上海、滄州、德州等地的地面沉降機理研究及模擬預測[42-46]。駱勇等[47]分別采用基于滲流-沉降分步計算的土力學經驗公式、滲流-沉降部分耦合的GMS軟件中SUB模型和基于滲流-沉降完全耦合的COMSOL Multiphysics模型，對疏排水引起的地面沉降進行了模擬計算，研究結果表明，完全耦合模型計算結果更為合理、更符合實際沉降特征。Wang等[48]探索研究了隨機非均質多孔介質中地下水流動和土層變形的動態相互作用。Pham等[49]基于多孔介質彈性力學理論，開發了程序模塊SUB+，可以實現水-應力的完全耦合模擬。TOUGH2是一套功能強大、應用廣泛的模擬孔隙或裂隙介質中多相流的系列程序，在水氣兩相流數值模擬中得到了廣泛應用[50-51]。研究人員以TOUGH2為基礎進行二次開發[52]或者將其與其它地質力學模擬程序結合[53]，實現了地下水流與地質力學的耦合模擬，顯示出該程序在土-水耦合模擬方面有廣泛的應用前景。

2 華北平原地面沉降現狀

我國地面沉降主要發生在華北平原、長江三角洲、汾渭盆地、淮河平原、珠江三角洲、東北平原以及山區斷陷盆地等地區，共有22個省（區、市）超過100個地級以上城市發生地面沉降。其中，華北平原是目前我國地面沉降影響面積最大、沉降速率最快的區域[54]。華北平原人均水資源量296 m3，為全國人均水資源的1/8，與干旱缺水國家以色列的人均水資源量相當，屬于水資源極度短缺地區。因地表水嚴重短缺，生產生活供水主要依賴地下水[55]，2000年以來，地下水占總供水量的比例達到64%。南水北調工程一定程度上緩解了供水緊張局面，但依然不能滿足社會經濟快速發展需求。地下水長期超采形成多個地下水降落漏斗，并引發嚴重的地面沉降，累計沉降量較大地區與深層地下水漏斗區的分布基本吻合，天津、滄州、衡水、德州地面沉降已經連成一片（圖1）。依據累計地面沉降量和地面沉降速率兩個指標（表1），對京津冀平原地面沉降現狀發育程度進行評價，結果表明地面沉降發育程度強的區域主要分布在北京平原東部、天津平原中南部、河北平原滄州、衡水、廊坊、邢臺、邯鄲等地，總面積達1.4×104km2（圖2）。

圖1 京津冀平原主要淺層和深層地下水漏斗以及累計沉降量較大地區（大于50 cm）分布圖Fig.1 Distribution of shallow and deep groundwater cones and the area with the cumulative subsidence greater than 50 cm in Beijing—Tianjin—Hebei Plain

表1 地面沉降現狀發育程度評價標準Table 1 Evaluation standard of the development of land subsidence

圖2 京津冀平原地面沉降現狀發育程度分區圖Fig.2 Zoning map of land subsidence development in Beijing—Tianjin—Hebei Plain

分層標監測數據顯示，北京市100 m以下深部地層沉降貢獻率為82%，其中天竺、張家灣等地深部地層貢獻率超過90%，淺部部分地層出現回彈。天津市300 m以下深部地層貢獻率占65%，其中漢沽營城、西青鄭莊子等地深部地層貢獻率超過70%，淺部部分地層出現回彈。河北省實施深層地下水禁限采的衡水、廊坊等地區，深部地層對地面沉降的貢獻率逐漸減少，其中衡水市區150 m以深地層的沉降貢獻率由2009—2015年的66%降至2019年的34%。必須高度重視地下水尤其是深層地下水超采問題，開展地面沉降防控，減少經濟損失。

近年來，南水北調工程運行為地下水超采治理提供了置換水源，減少了地下水開采量，工程沿線大中型城市區地下水位下降趨勢得到遏制，地面沉降防控取得了初步成效，天津、滄州、衡水等重點城市主城區地面沉降得到有效控制，北京、天津地面沉降嚴重區面積（本文地面沉降嚴重區指年沉降量大于50 mm的地區）總體呈現減少趨勢（圖3）。南水北調進京5年以來，北京地下水位整體處于恢復上升狀態，尤其是100 m以淺的含水層，水位回升幅度較大，與2014年相比，平原區地下水位平均上升了3 m，地下水降落漏斗面積減少了38%（圖4）。地下水位恢復導致地面沉降減緩，地面沉降嚴重區面積減少了71%。滄州市自2005年開始實施關停單位自備井的禁采措施，市區年沉降量已由60～80 mm降至目前的10～15 mm。滄州市區分層標組于2010年6月開始監測，共由五座分層標組成，分別監測第一壓縮層（5～69 m）、第二壓縮層（69～196 m）、第三壓縮層（196～253 m）、第四壓縮層（253～375 m）。2014年第一至四壓縮層的變形量分別為?2.4，?5.4，?4.6，?8.9 mm（負號表示壓縮，正號表示回彈，下同），2019年地面沉降進一步減緩（圖5），四個壓縮層的變形量分別為+0.6，?2.7，?2.2，+4.2 mm。2013—2018年，河北衡水市主城區年平均沉降量由15 mm減至10 mm。衡水分層標主要監測第一壓縮層（41～150 m）、第二壓縮層（150～267 m）、第三壓縮層（267～401 m），其中第一壓縮層對應淺層含水層組，第二、三壓縮層對應深層含水層組。2009年10月—2015年8月，第一至三壓縮層的平均變形速率為?19.3，?31.3，?6.1 mm/a，對地面沉降量的貢獻率分別為34%、55%、11%。衡水地區從2014年開始進行地下水超采區綜合治理，大幅壓減地下水開采量，衡水市城區禁采深層地下水，深層地下水水位下降減緩甚至回升，有效遏制了地面沉降的快速發展。2019年第一至三壓縮層的變形量減至?10.9，?3.4，?2.1 mm，對地面沉降量的貢獻率為66.2%、20.8%、13%，主要壓縮層位從之前的以第三、四含水層組為主逐步轉變為以一、二含水層組為主。

圖3 北京、天津和河北平原地面沉降嚴重區面積占各省（市）平原區面積比例（2012—2019年）Fig.3 The ratio between the severe subsidence area and the total area for Beijing,Tianjin and Hebei plain（from 2012 to 2019）

圖4 北京平原地面沉降嚴重區、地下水漏斗面積占平原區總面積比例與平均水位埋深（2012—2019年）Fig.4 The ratio of area of severe subsidence and groundwater depression cones to the total plain area and the average waterlevel depth for Beijing plain（from 2012 to 2019）

圖5 滄州市區各壓縮層累計變形隨時間變化曲線（2019年）Fig.5 The temporal variation of the deformation of different compression layers in the downtown of Cangzhou（2019）

雖然華北平原地面沉降出現減緩態勢，北京、天津、河北滄州等治理區地面沉降速率大幅下降，但華北平原尤其是河北平原地面沉降總體上仍然處于較快發展階段，地面沉降防控形勢依然嚴峻。為滿足農業用水需求，山前平原超采淺層地下水，中東部平原的衡水、滄州等地超采深層地下水，農業區超采是當前華北平原地面沉降嚴重的主要影響因素之一。2019年河北全省平均降水量437.6 mm，比2018年減少71.6 mm，為確保農業生產，加大了灌溉期地下水開采力度，使農灌區地下水位普遍下降。據中國地質調查局地下水位統測數據，2019—2020年，華北平原淺層和深層地下水位平均分別下降0.6 m和1.0 m，水位下降區主要分布在農業灌區，加劇了地面沉降發展。目前，華北平原地面沉降嚴重區面積超過全國沉降嚴重區總面積的80%，沉降速率較快地區主要分布在北京朝陽與通州的交界處，天津武清、靜海、西青、北辰，河北廊坊、滄州、保定、衡水、邯鄲等地。河北雄安新區、北京城市副中心等地存在地面沉降嚴重區，京津城際、京滬高鐵等高速鐵路穿越地面沉降嚴重區，應予以高度重視。

3 研究展望

地下水位變化驅動下的土層變形特征及其機制研究仍屬國內外地面沉降研究的前沿課題。如前所述，土層變形規律、水文地質參數反演校正和土水耦合模型應用等三個研究方向聯成一個統一的整體，近年來在基礎理論上有了較大的發展，在技術方法上更加注重利用InSAR等新技術、新方法獲取的數據開展研究，為華北平原地面沉降防治研究提供了很好的經驗借鑒。同時，華北平原具有獨特的水文地質條件，是我國含水層結構最為復雜的大型平原，受長期地下水超采等人類活動影響，形成了世界上面積最大的地下水降落漏斗和最嚴重的地面沉降區。不同地區地質條件的復雜性和時空上的差異性導致許多研究結論可能無法在該地區適用，需要不斷提高地下水與地面沉降監測技術水平與監測精度，因地制宜地開展相關研究工作。在總結地下水位變化模式和土層變形規律的基礎上，研究地面沉降差異性特征及變形成因機理，將大幅提升華北平原地面沉降機理的認識水平，極大促進環境地質學科發展，為地面沉降防災減災提供重要技術支撐。

（1）深化地面沉降機理和預測預警研究。華北平原構建了地面沉降監測網，監測面積達9×104km2，覆蓋了京津冀整個平原區；其中，建成基巖標14個，分層標31組，水準監測點3 677個，GPS監測點272個，地下水監測井5 000余眼，InSAR實現了京津冀平原區全覆蓋連續動態監測。由于沉積過程、沉積相、地層巖性及應力加載過程條件（地下水位變化模式）等不同，華北平原不同埋深土層在地下水位變化影響下的變形特征存在較大差異。多年積累的地下水位、土層變形等信息對于地面沉降預測預警來說無疑十分寶貴，可是迄今對這些數據缺乏科學化、規律化的分析和提煉，沒有對不同地下水位變化模式下土層變形規律這一關鍵問題進行系統的研究。如何聯合利用室內土工試驗和基于現場調查監測數據的數理分析揭示土層變形應力應變規律？如何量化這一規律從而形成簡單實用的土層變形本構關系，甄別誘發地面沉降快速發展的地下水位及其下降速率等地面沉降控制指標？如何將土層變形本構模型與地下水流模型耦合用于地面沉降監測預警？對這些科學問題的解決，直接關系到對華北平原地面沉降形成機理的認識，也是區域地面沉降預測預警的迫切需求，有助于因地制宜地提出有針對性的地面沉降防控措施。圖6概括了地下水位變化驅動下土層變形規律的研究思路。

（3）提升對華北平原地下水資源屬性的認識。地下水動力與土層變形的互饋關系研究是地面沉降區地下水資源量評價的前提條件。以華北中東部平原為例，深層地下水補給條件較差，可更新能力低，在相同地下水開采強度下，更容易產生地面沉降。近40年以來，深層地下水位最大累計降幅已經超過100 m，目前天津、滄州、衡水、德州深層地下水漏斗已經連成一片，面積超過2.0×104km2。華北平原深層地下水開采量中含水層和弱透水土層壓密釋水量所占比例達40.1%～43.8%，滄州地區高達57.6%[56]。正確評價地面沉降區地下水資源量，需要充分利用地面沉降觀測數據，科學確定沉降層位、沉降量與地下水開采量的定量關系。

圖6 地下水位變化驅動下土層變形規律研究思路Fig.6 Research ideas on groundwater-derived deformation law of soil layers

（4）加強地熱開發與地面沉降關系研究。京津冀地區地熱資源開發始于20世紀70年代初，21世紀以來開發利用規模不斷擴大，目前是我國地熱資源開發利用程度最高的地區，主要開發新近系明化鎮組、館陶組、古近系東營組3個孔隙熱儲層和寒武-奧陶系、薊縣系2個巖溶熱儲層。地熱水開發利用對地面沉降的影響目前尚無定論，主要原因是缺少觀測數據和研究成果支持。深入研究地熱水開發引起的水位變化與土層變形的關系，提出合理的水位控制指標，可以為實現地熱資源持續利用和地面沉降有效控制提供理論依據和數據支撐。