基于隨機模擬的渾河沖洪積扇地區地下水壓采風險預報

蘇小四，杜守營，杜尚海，宋憲宗，邵廣凱，王 璜

1.吉林大學地下水資源與環境教育部重點實驗室，長春 130021

2.吉林大學水資源與環境研究所，長春 130021

0 引言

地下水數值模擬技術在地下水有關領域的廣泛使用，使其在地下水資源評價、開發利用、管理和規劃等方面起到了重要作用［1］。但地下水系統固有的隨機性以及研究者主觀認識的局限［2］，使數值模擬采用的參數帶有很大的不確定性，由此可導致模型計算結果的不確定性［3－4］。近年來考慮參數不確定性的地下水流數值模擬的研究，在國內外得到了較好的應用和發展［5］，且通過對相關水文地質參數的不確定性分析，可在建立地下水流數值模型的基礎上進行隨機預報，預報結果可為地下水資源管理提供更可靠依據。

由于長期過量開采地下水，渾河沖洪積扇地區已形成了大面積地下水位降落漏斗，截止到2010年的統計數據表明，僅沈陽城區33m地下水等水位線圈閉的漏斗區面積已達到150km2［6］。為有效地緩解城市地下水超采區地下水位持續下降問題，防止地下水環境的進一步惡化，地下水壓采成為解決該類問題的主要方案之一［7］，準確預測不同壓采方案下地下水水位的恢復情況對超采區地下水資源規劃和管理具有重要的實際意義［8－9］。筆者以渾河沖洪積扇地區為研究區，在水文地質參數不確定性分析和地下水流數值模型的基礎上，進行壓采條件下地下水水位恢復效果的隨機預報，并根據預報結果對地下水壓采方案進行風險評價。

1 研究區概況

研究區位于渾河沖洪積扇上，地勢由東北向西南逐漸降低，平均坡度7.5‰，海拔45～50m，如圖1所示。區內地貌以沖洪積平原為主，東北及東南部屬風積崗狀與波狀臺地，除東部丘陵區外，均被第四系松散堆積物覆蓋。多年平均降水量587.5 mm，多年平均蒸發量826.8mm。渾河是流經該區最大的河流，多年平均流量45.6m3/s，從渾河沖洪積扇的扇頂到扇緣，地表水與地下水轉換頻繁，在研究區內主要由地表水滲漏補給地下水，局部地段表現為地下水補給地表水。

區內地表為第四系黏性土和亞砂土，其下分別為：以第四系全新統砂礫石和卵石為主的孔隙潛水含水層（Q4）、第四系全新統砂卵石為主的孔隙微承壓含水層（Q3）、第四系中、下更新統砂卵石為主的孔隙局部承壓含水層（Q2＋1）及第三系的風化基巖裂隙水（圖2）。其中：孔隙潛水和孔隙微承壓含水層在區域內大面積連續分布，且二者水力聯系密切，因此將其進一步概化為潛水（微承壓）含水層；第四系孔隙潛水（微承壓水）和孔隙承壓水水頭近一致，兩者之間是由黏土及亞黏土構成不連續的弱透水層，該弱透水層在區內大部分地段存在缺失；承壓含水層底板為相對隔水的基巖；區域地下水流向從東北流向西南，在集中開采的中山、李官堡和競賽等供水水源地附近存在多處地下水位降落漏斗。

2 研究方法

2.1 水文地質概念模型及數值模型

本次利用地下水模擬系統（groundwater modeling system，GMS）軟件建立了研究區地下水流數值模型。在概化研究區水文地質概念模型的基礎上建立地下水流數值模型，具體過程如下。

1）含水層系統概化：對區內的鉆孔資料進行整理分析，查明上層孔隙潛水（微承壓水）與下層孔隙局部承壓含水層之間的亞黏土層不連續，且二者之間水頭較為一致，存在較為密切的水力聯系，因此本次將研究區第四系孔隙潛水（微承壓水）含水層與孔隙承壓含水層系統概化為一層。

圖1 研究區位置及概念模型示意圖Fig.1 Study area and conceptual model diagram

2）邊界條件概化：水平方向上，模擬區東北及西南部地貌類型屬風積崗狀臺地與風積波狀臺地，地下水位變化不大，側向補給量較穩定，研究中確定為二類邊界（圖1）；西部邊界因其延伸方向與地下水流向近似平行，且多年變化不大，研究中確定為零流量邊界；北部及東南部邊界附近雖然地下水動態變化較大，但豐富的地下水水位動態監測數據可以較為精確地刻畫邊界附近水位的變化規律，研究中確定為一類邊界；垂向上，含水層上部主要接受大氣降水補給，為水量交換邊界，下部因與相對隔水的黏土層、基巖接觸，所以定為隔水邊界。

3）源匯項概化：含水層的補給主要有降水入滲、河水入滲、地表灌溉回滲、渠系滲漏和地下水側向徑流等；地下水的排泄途徑主要以人工開采和潛水蒸發為主，其中開采井通過調整過濾器的位置控制開采層位及開采量，蒸發排泄量通過潛水位與蒸發高程和極限蒸發深度的關系自動計算；河流入滲補給量則通過“River”模塊處理成源匯項，根據渾河附近監測孔的地下水水位觀測值、河床沉積物滲透系數和厚度等參數對河流屬性賦值，由模型自動計算地下水與地表水之間的補排量［10］。

4）數值模型建立：研究區水力坡度平均小于1.0‰，地下水以水平運動為主，流速緩慢，滲流基本符合線性達西定律，且區內含水層的巖性和厚度均有不同程度變化，故模擬區地下水為非均質、各向同性的二維非穩定流，其數學表達式為

式中：K為潛水含水層滲透系數（m/d）；H為潛水水位（m）；z為含水層底板標高（m）；Kz為河床沉積物的垂向滲透系數（m/d）；dz為河床沉積物的厚度（m）；Hz為河流的水位（m）；Qr為補給強度（m/d）；Qe蒸發排泄強度（m/d）；Qi為開采強度（m/d）；μ為潛水含水層給水度；H0為初始水位（m）；H1為一類邊界點的水位（m）；q為二類邊界單寬流量（m2/d）；x、y為坐標（m）；D為計算區范圍；Γ1，Γ2分別為一、二類邊界。

2.2 參數靈敏度分析

靈敏度是衡量地下水流數值模型對參數變化敏感程度的指標，可以反映一種因子的變化對另一因子的影響程度［11］。限于模型概化及調參階段對水文地質條件認識的主觀局限性，模型識別后的率定參數具有較大的不確定性。通過靈敏度分析可對數值模型各參數對地下水位變化的影響程度排序，定性或定量地評價參數不確定性對模擬結果的影響［12－13］，識別出對系統總體輸出和動態影響較大的參數，進而確定隨機模型的變量參數。

靈敏度分析采用因子變換法，首先確定某一參數作為變量因子，同時保持其他參數不變，逐次運行模型后選定驗證時段末刻，分別輸出各個參數分區中典型觀測井的地下水水位值Hi，作為靈敏度分析的評價指標［3，13］。當作為變量因子的參數由初值αk變為（αk＋Δαk），相應地，因變量Hi（αk）變為Hi（αk＋Δαk）時，單指標的參數靈敏度系數計算公式為

式中：βi，k為水位H對第k個參數在第i個觀測點上的靈敏度系數；Hi為i點的水位；αk為第k個參數的初值。

為便于比較不同參數的靈敏度系數，需對式（2）求得的靈敏度系數進行無量綱處理：

2.3 地下水流預報的隨機模型

地下水流預報隨機模型的建立調用了GMS軟件中的“Stochastic Modeling”模塊，參數隨機化過程基于蒙特卡羅原理，其前提是假定滲透系數、給水度等隨機變量符合一定的概率分布規律，例如滲透系數一般被認為符合對數正態分布，而給水度服從均勻分布［14］。在設置好各個參數概率分布的條件下，首先為每個參數指定均值、標準差、最大值和最小值作為其概率分布的統計特征，然后對參數做多次抽樣試驗，生成多組輸入變量的隨機組合，對應每組隨機變量分別運行一次 MODFLOW，最后得到若干組水流模型的計算結果，由此獲得區內任意點計算水位的統計值［15］。

考慮到隨機因素的不確定性對預報結果可能產生波動性的影響，本次研究在參數靈敏度分析的基礎上，建立了研究區地下水位預報的隨機模型。預報結果不僅能夠對壓采方案下區域地下水位的動態變化趨勢進行預測，還可以根據典型觀測井水位的恢復效果對區域地下建筑工程的潛在風險進行評價。

3 結果與討論

3.1 地下水流數值模擬結果分析

本次研究將渾河沖洪積扇地區剖分成5 025個有效單元格，平均每個單元格面積為0.2km2。模擬期時間步長為1個月。其中：2009年4月25日－2010年4月25日為識別期，以30d為1個應力期，共12個應力期；2010年4月25日－2011年10月25日為驗證期，以30d為1個應力期，共18個應力期。

模型識別及驗證期典型觀測孔的計算水位與實測水位擬合效果如圖3所示。可以看出，計算水位與實測水位擬合較好，擬合結果的平均誤差為0.05 m，平均絕對誤差為0.27m，平均方差為0.31，擬合誤差小于0.50m的觀測井占其總數的84.50%。由此可見，所建模型基本能正確地反映區內地下水流系統的滲流規律，可在此基礎上對模型率定的水文地質參數進行敏感性分析。從計算結果還可以看出，現狀開采條件下，在中山、于洪、競賽、李官堡等水源地出現分散的地下水水位降落漏斗。

3.2 參數靈敏度分析

圖3 驗證期末刻（2011年10月25日）地下水流場擬合及水文地質參數分區圖Fig.3 Fitting curve of groundwater flow field and hydrological parameters zoning at the end of the verification period（Oct.25，2011）

結合研究區水文地質條件，本次研究選取滲透系數、給水度、降水入滲補給系數及河床沉積物的滲透系數4個參數進行靈敏度分析。考慮到區內含水層巖性以粗砂和礫石為主，滲透系數應為20～150 m/d，給水度應為0.15～0.35［16］；而驗證后的滲透系數為33～121m/d，給水度為0.11～0.29，水文地質參數分區情況如圖3所示。相對于滲透系數和給水度，河床沉積物的滲透系數和降水入滲補給系數無明顯的制約條件；因此在保證滲透系數和給水度取值范圍合理的前提下，確定各個參數的變化幅度在率定值的±30%之間。

根據研究區滲透系數分區和靈敏度計算方法，可以得到滲透系數的靈敏度計算結果（圖4）。從圖4可以看出，第5分區的滲透系數最為敏感，主要是因為該分區為渾河沖洪積扇的河漫灘，含水層巖性變異顯著，其他3個參數的計算結果也顯示第5分區的參數靈敏度系數比較大。因此本次研究對第5分區各個參數靈敏度分析結果做了對比分析，結果如圖5所示。可以看出：地下水水位對含水層滲透系數的變化最敏感；其次是給水度，而對河床沉積物滲透系數和降雨入滲補給系數的靈敏性較差；滲透系數和給水度在其率定值附近對稱性增加或減少時，靈敏度系數也呈現出對稱性的變化規律，即參數變化幅度增大時，靈敏度系數隨之增大。

圖4 各分區滲透系數靈敏度對比Fig.4 Permeability coefficient sensitivity of each partition

3.3 地下水流隨機預報結果分析

如果不考慮參數隨機性，則預報模型中水文地質參數是確定的，因此其模擬結果相當于隨機模擬的一個特例，不具有統計特征；而采用參數隨機模型進行地下水位預報，能更真實地逼近客觀的水文地質條件，統計結果能提高預測結果的精度，從而減少決策誤差。本次研究在參數靈敏度分析的基礎上，選取滲透系數和給水度作為隨機變量，并確定隨機變量的變化范圍，其中滲透系數的變化范圍為率定值的±30%之間，給水度的變化范圍為率定值的±20%之間，再為每個隨機變量指定初始值、最大值、最小值、均值及標準差，其中初始值取模型驗證階段得到的參數率定值，標準差取區間長度的20%，相關的隨機參數設定如表1所示。

圖5 第5區參數靈敏度分析結果Fig.5 Parameter sensitivity results of the fifth partition

表1 隨機參數的輸入設定Table 1 Random parameter settings

地下水流模型中一類邊界的地下水位是由大氣降水、蒸發等氣象因素及人類活動（主要是地下水開采）對影響半徑內各動態監測井的水位影響共同決定的，因此進行一類邊界條件水位預測時，分別考慮了上述2種影響因素后復合疊加處理［17］。根據Theis公式近似計算出設計開采量在一類邊界上各動態觀測孔處引起的地下水位降深S［18］，預報出各孔在人類活動影響后的末刻水位。大氣降水量采用歷史降水周期重現法進行預報［11］，在對模擬區降水統計資料分析后，設定2006－2022年的降水量為1964－1980年之間降水量的歷史重現，由此確定預報期內一類邊界在自然和人類活動共同影響下的水位以及研究區內的降水量［7］。蒸發以及河流水位均取多年平均值。預報期地下水的設計開采量是根據研究區未來地下水資源規劃方案、在現狀開采量的基礎上、調整開采量以每年5%的速度遞減、以點狀開采井的形式輸入到模型中的。

地下水流預報模型以驗證期末刻（2011年10月25日）的地下水位作為預報期的初始流場，預報時長為10a，在GMS軟件中更改模型的運行方式為“Stochastic Simulation”，設置模擬次數為200，則相應地生成200組隨機組合數并分別賦值到數值模型中。運行可得200組地下水隨機流場，即模型中任意節點在模型末刻均可計算出200個隨機水位。

考慮到本次地下水流預報是在限制地下水開采的條件下進行，因此選取以水源地為降落漏斗中心的觀測井，將預報末刻（2022年10月25日）作為時間節點，分別輸出中山水源地3＃觀測井和李官堡水源地9＃觀測井的隨機預報水位。通過SPSS 15.0軟件對輸出結果進行統計，可得到置信區間為95%的條件下地下水位的平均值、最大值、最小值以及標準偏差等信息（圖6）。從圖6可以看出，受最敏感的滲透系數影響，隨機模型的預報水位也基本符合正態分布規律。

選取預報期內每年10月25日作為時間節點，對第5分區部分水源地典型監測井的隨機預報結果進行統計分析，可得到預報期內平均地下水位年際動態變化趨勢（圖7）。從圖中可以看出：在地下水開采量以每年5%的速度壓采下，相比于預報初始（2011年10月25日），區內地下水水位得到了顯著恢復，平均上漲3.3m；地下水水位上升幅度最大的是中山水源地，說明地下水壓采可以緩解區域地下水水位下降及降落漏斗擴大等問題。

從壓采預報結果可以看出，雖然地下水壓采可以有效緩解地下水位持續下降帶來的環境問題，但當地下水水位恢復超過地下建筑工程安全設計水位時，地下建筑物的安全將會受到影響，因此需要對地下水壓采給周邊地下建筑帶來的潛在風險進行評價。本次研究中的風險性分析是指地下水開采決策過程中，在系統輸入有不確定因素的情況下，研究某一范圍內地下水位（H）在一時間段內超過某給定安全水位（Hi）的可能性，即當H≥Hi時的概率P［19－20］。通過對隨機預報結果進行統計，可以得到研究區任意一點在預報期內水位超過設計安全水位的概率，從而定量地評價給定的開采方案對地下水建筑安全的風險性。

以水位上升幅度最大的中山水源地為例，到預報末刻中山水源地3＃觀測井水位恢復情況如表2所示。可以看出：當中山水源地設計的安全水位是29.8m時，在預報期的第10年的滲水風險概率為1.5%，為小概率事件，發生的可能性較小；當該處設計的安全水位是28.8m時，在預報期的第6年（2018年）出現滲水風險的概率為6.5%，而第10年（2022年）發生滲水的風險概率為100.0%。由此可見，按照設計的壓采方案對研究區地下水進行開采，未來地下水位的恢復會不同程度地對附近的地下建筑構成滲水威脅，設計安全水位Hi越低，風險出現的時間越早。因此，考慮參數隨機性的地下水位風險預報結果能為區域地下水資源風險決策提供依據。

圖7 預報期內地下水水位動態變化趨勢Fig.7 Dynamic changing trends of groundwater level during the forecast period

表2 中山水源地預報期內不同安全水位對應的風險概率Table 2 Risk probability under different security groundwater levels during the forecast period of Zhongshan water source P/%

4 結論

1）參數靈敏度分析結果表明，地下水水位對含水層滲透系數的變化最敏感，其次是給水度，而對河床沉積物滲透系數和降雨入滲補給系數的靈敏性較差；滲透系數和給水度在其率定值附近增加或減少時，靈敏度系數隨之增加或減小。

2）研究表明，壓縮開采地下水資源能夠有效緩解地下水水位下降帶來的環境問題，地下水開采量以每年5%的速度壓采時，區內地下水水位平均上漲3.3m，但水位恢復的同時也可能誘發局部地下工程滲水，且地下建筑物的設計安全水位越低，滲水風險越大。

（

）：

［1］吳吉春，陸樂.地下水模擬不確定性分析［J］.南京大學學報：自然科學，2011，41（3）：227-234.

Wu Jichun，Lu Le. Uncertainty Analysis for Groundwater Modeling ［J］.Journal of Nanjing University：Natural Sciences，2011，41（3）：227-234.

［2］束龍倉，朱元生.地下水資源評價中的不確定性因素分析［J］.水文地質工程地質，2000（6）：6-8.

Shu Longcang， Zhu Yuansheng. Analysis of Uncertainty Factors in Evaluation of Groundwater Resources［J］.Hydrogeology ＆ Engineering Geology，2000（6）：6-8.

［3］劉佩貴，束龍倉.傍河水源地地下水水流數值模擬的不確定性［J］.吉林大學學報：地球科學版，2008，38（4）：639-643.

Liu Peigui，Shu Longcang.Uncertainty on Numerical Simulation of Groundwater Flow in the Riverside Well Field［J］.Journal of Jilin University：Earth Science Edition，2008，38（4）：639-643.

［4］黃修東，杜尚海，崔峻嶺，等.平谷地下儲水空間調蓄能力及不確定性計算［J］.工程勘察，2013（7）：35-39.

Huang Xiudong，Du Shanghai，Cui Junling，et al.Computation and Uncertainty Analysis of Groundwater Reservoir Regulation Capacity in Pingu Area［J］.Geotechnical Inverstigation ＆Surveying，2013（7）：35-39.

［5］杜守營，鹿帥，杜尚海.基于GMS的地下水流數值模擬及參數敏感性分析［J］.中國農村水利水電，2013（8）：77-80.

Du Shouying，Lu Shuai，Du Shanghai.Numerical Simulation of Groundwater and Sensitivity Analysis of Parameters Based on GMS［J］.China Rural Water and Hydropower，2013（8）：77-80.

［6］周浩，王殿武，孫才志，等.沈陽中心城區水源地漏斗恢復調蓄模擬研究［J］.水文，2011，31（5）：47-51.

Zhou Hao，Wang Dianwu，Sun Caizhi，et al.Numeric Simulation Research on the Recovery of Groundwater Funnel in Core Region of Shenyang City［J］.Journal of China Hydrology，2011，31（5）：47-51.

［7］杜守營，袁文真，杜尚海，等.壓采條件下地下水位恢復對沈陽市地下工程安全影響分析［J］.中國農村水利水電，2013（5）：20-24.

Du Shouying，Yuan Wenzhen，Du Shanghai，et al.Impacts of Water Table Recovery by Reducing Groundwater Abstraction on Underground Project Safety in Shenyang City［J］.China Rural Water and Hydropower，2013（5）：20-24.

［8］杜尚海，蘇小四，呂航.不同降水豐枯遭遇條件下滹沱河地下水庫人工補給效果研究［J］.吉林大學學報：地球科學版，2010，40（5）：1011-1018.

Du Shanghai，Su Xiaosi，Lü Hang.The Artificial Recharge Effects of Groundwater Reservoir Under Different Precipitation Plentiful-Scanty Encounter in Hutuo River［J］.Journal of Jilin University：Earth Science Edition，2010，40（5）：1011-1018.

［9］Du Shanghai，Su Xiaosi，Zhang Wenjing.Effective Storage Rates Analysis of Groundwater Reservoir with Surplus Local and Transferred Water used in Shijiazhuang City，China［J］.Water and Environment Journal，2013（27）：157-169.

［10］曹廣祝，強毅，李峰，等.李官堡水源地地下水流場數值模擬及預測［J］.中國水能及電氣化，2011（8）：7-13.

Cao Guangzhu，Qiang Yi，Li Feng，et al.Prediction and Numerical Simulation of Groundwater Flow Field of Liguanpu Wells［J］.China Water Power ＆Electrification，2011（8）：7-13.

［11］翟遠征，王金生，蘇小四，等.地下水數值模擬中的參數敏感性分析［J］.人民黃河，2010，32（12）：99-101.

Zhai Yuanzheng，Wang Jinsheng，Su Xiaosi，et al.Parameter Sensitivity Analysis in Numerical Simulation of Groundwater［J］.Yellow River，2010，32（12）：99-101.

［12］李文運，崔亞莉，蘇晨，等.天津市地下水流－地面沉降耦合模型［J］.吉林大學學報：地球科學版，2012，42（3）：806-813.

Li Wenyun，Cui Yali，Su Chen，et al.An Integrated Numerical Groundwater and Land Subsidence Model of Tianjin［J］.Journal of Jilin University：Earth Science Edition，2012，42（3）：806-813.

［13］束龍倉，王茂枚，劉瑞國，等.地下水數值模擬中的參數靈敏度分析［J］.河海大學學報：自然科學版，2007，35（5）：491-495.

Shu Longcang，Wang Maomei，Liu Ruiguo，et al.Sensitivity Analysis of Parameters in Numerical Simulation of Groundwater［J］.Journal of Hohai University：Natural Sciences，2007，35（5）：491-495.

［14］Freeze R A.A Stochastic Conceptual Analysis of Dimensional Groundwater Flow in Nonuiform Homogeneous Media ［J］. Water Resources Research，1975，11（5）：725-741.

［15］劉猛，束龍倉，劉波.地下水數值模擬中的參數隨機模擬［J］.水利水電科學進展，2005，25（6）：25-27.

Liu Meng，Shu Longcang，Liu Bo.Stochastic Parameter Simulation in Numerical Simulation of Groundwater ［J］. Advances in Science and Technology of Water Resources，2005，25（6）：25-27.

［16］王大純，張人權，史毅虹，等.水文地質學基礎［M］.北京：地質出版社，1995.

Wang Dachun，Zhang Renquan，Shi Yihong，et al.Hydrological Geology Basis［M］.Beijing：Geological Publishing House，1995.

［17］盧文喜.地下水運動數值模擬過程中邊界條件問題探討［J］.水利學報，2003（3）：33-35.

Lu Wenxi.Approach on Boundary Condition in Numerical Simulation of Groundwater Flows［J］.Journal of Hydraulic Engineering，2003（3）：33-35.

［18］薛禹群，朱學愚，吳吉春，等.地下水動力學［M］.北京：地質出版社，1997.

Xue Yuqun， Zhu Xueyu， Wu Jichun，et al.Groundwater Dynamics［M］.Beijing： Geological Publishing House，1997.

［19］崔亞莉，王婉麗，楊廣元.基于LHS方法的地下水隨機模擬及開采量風險性評價［J］.地學前緣，2010，17（6）：134-139.

Cui Yali， Wang Wanli， Yang Guangyuan.The Parameter Stochastic Simulation and Exploitation Risk Analysis Based on LHS Method［J］.Earth Science Frontiers，2010，17（6）：134-139.

［20］陸樂，吳吉春.地下水數值模擬不確定性的貝葉斯分析［J］.水利學報，2010，41（3）：264-271.

Lu Le， Wu Jichun. Bayesian Analysis of Uncertainties in Groundwater Numerical Simulation［J］.Journal of Hydraulic Engineering，2010，41（3）：264-271.