通量上邊界與水頭上邊界方法的地下水流系統模擬對比

牛 宏，梁 杏，2，張人權

1.中國地質大學環境學院，武漢 430074

2.中國地質大學生物地質與環境地質國家重點實驗室，武漢 430074

0 引言

1940年，Hubbert通過勢的數學分析得出了具有上升和下降水流的河間地塊示意流網圖，打破了傳統認為的河間地塊地下水水平運動的定勢，這是地下水流系統理論出現的雛形［1］。之后，Tóth在加拿大阿爾伯達中部地區觀察到有別于Hubbert流網圖的地下水流現象，通過分析Hubbert等前人的成果，在均質各向同性介質與嚴格假定的邊界條件下，利用解析解繪制了潛水盆地中具有局部、中間及區域的嵌套式多級地下水流系統，創造性地建立了地下水流系統理論［2－7］，并且認為地形勢是控制盆地地下水流系統的關鍵因素。張人權等［8－9］指出，地下水流系統理論是構建人和自然協調的、良性運轉系統的理論基礎和有效工具，地下水流系統理論是當代水文地質學的核心概念框架。侯光才［10］、陳蓓蓓［11］應用地下水流系統理論分析了鄂爾多斯白堊系盆地地下水流系統及北京市地下水流系統的演化，對該區各級水流系統特征有了初步的認識。

劉宇、梁杏等［12－13］在進一步開展地下水流系統模式研究時發現，在Tóth地形勢控制的模型中，給定盆地上邊界水頭條件（簡稱水頭法），雖然容易獲得解析解，便于理想模型的精確求解，但當單獨改變盆地的其他因素（滲透性、盆地深度等）時，盆地補給和排泄也發生變化，水頭法不能很好地分析與理解單一因素對地下水流系統模式發育的影響；給定上邊界水頭固化了盆地的勢源與勢匯的位置及數目，這與實際盆地潛水面和勢匯的變化不相符合，也限制了地下水流系統模式的轉化研究。Haitjema和Mitchell-Bruker［14］對地形控制潛水地下水位提出了質疑，Mitchell-Bruker應用定流量上邊界的邊界元法，得出即使是復雜地形的盆地，改變入滲補給量，也可以出現與Tóth方法不同的水流模式，只有當入滲補給量對介質滲透性的比值較大時，地形才對地下水位起控制作用。為此，梁杏等［15－17］提出了通量上邊界研究方法（簡稱通量法），并開展了相關的物理模擬和數值模擬。

水頭法在進行水流系統研究中存在的問題究竟如何，與通量法偏差有多大，應用時限制條件如何，還沒有定量的對比模擬分析報道。筆者旨在用數值法進行水頭法和通量法的對比模擬，通過模擬結果的定量比較，具體指出水頭法模擬的問題與局限性。

1 Tóth典型地下水流系統求解簡介

Tóth分別采用一次函數和一次函數疊加正弦波函數，解析出簡單與復雜盆地地下水流系統模式［2］。其中，采用一次函數疊加正弦波函數的水頭法進行盆地二維流模型的數學方程如下：

式中：h為地下水水頭；x、z分別代表以谷底左下方為起點到盆地地下水流動區域中任何一點的水平和垂直距離；Kx、Kz分別為x和z方向的滲透系數；Lx表示x方向總距離；Lz表示盆地底部到盆地最低排泄谷底的深度；h0為上邊界定水頭值；tanα為盆地谷翼平均坡角正切值（盆地地形斜率）；a為波幅；b＝2π/λ。Tóth利用上述模型，取Lx＝6 096.00 m，Lz＝3 048.00m，a＝15.24m，tanα＝0.02，計算出存在5個高地勢和5個低地勢的波狀上邊界；繪制了盆地多級地下水流系統圖（圖1）；認為在地形復雜的均質各向同性潛水盆地中，地下水流呈現多級次特點，并分別定義為局部的（local system，LS）、中間的 （mittle system，MS）及區域的 （regional system，RS）地下水流系統。

圖1 地形復雜的均質各向同性盆地中的多級次水流系統Fig.1 Multilevel flow systems of homogeneous and isotropic basin in complex terrain

方程組（1）和圖1顯示，在Tóth的求解方法下，滲流區簡化為矩形區域，采用水頭法求解計算時潛水面、盆地勢源與勢匯的位置和數目不變；這樣求解與實際盆地條件不相符合。實際盆地中，補給區和勢匯位置已知時，潛水面（陡緩變化）受到補給量、盆地形態和介質等影響而變化，盆地滲流區僅與潛水面的變化保持一致。為克服解析解的局限性，筆者采用數值法，分別進行盆地水頭法和通量法的Tóth典型地下水流模型和變條件模型的改進模擬計算，并對比分析了改進前后的水流系統特征。

2 Tóth典型地下水流系統的數值模擬

2.1 水頭法的數值模擬

應用數值模擬軟件采用加拿大Waterloo水文地質公司在Modflow的基礎上應用現代可視化技術開發研制的Visual-Modflow，它是目前國際一致認可的三維地下水流和溶質運移模擬的標準可視化專業軟件系統［18］。

數值模擬方法對Tóth典型模式的滲流區域進行改進，其數學模型類同方程組（1），唯一的區別是模型上部滲流區與上邊界水頭變化一致，而非Tóth求解時的z＝Lz的水平線。數值模擬取Kx＝Kz＝0.10m/d，有效孔隙度0.25，給水度0.20，總孔隙度0.30，Lx＝6 096.00m，最低谷底高程zmin＝3 048.00m，x方向均分為66列，z方向均分為50層，對上部滲流區進行加密，第一層再加密均分為3層，共52層。模擬得出水流模式如圖2所示，采用水頭法的數值模擬結果與Tóth解析解的求解結果（圖1）基本一致，此時滲流區域最低排泄點為3 048.00m，潛水面最高水頭為3 170.00m。在此基礎上，筆者進行了通量法的模擬。

圖2 Tóth典型模式（滲流區修正）的水頭邊界數值模擬結果Fig.2 Given head boundary numerical simulation results of Tóth’s typical model with seepage zone modified

2.2 通量法的數值模擬

通量法的數值模擬是在水頭法基礎上，將Tóth典型模式（圖2）上邊界5個補給區（圖1中1、3、5、7、9勢源區）的補給量轉換成平均入滲補給量，在保持盆地補給條件一致的情況下，用通量上邊界代替定水頭上邊界進行模擬。其中：排泄區按照Tóth典型模型的5個排泄點對應設置為5個排泄溝，x方向均分為66列；對5個排泄溝進行加密，每格加密均分成4格，共81列；z方向均分為50層，對上部滲流區進行加密，第一層加密均分為3層，共52層。通量法的模擬結果如圖3所示，最低谷底排泄點水頭為3 048.00m，潛水面最高水頭為3 164.00 m，模擬的實際潛水面較水頭法和緩些。對比水頭法（圖2）與通量法（圖3）結果，盡管2種方法的水流模式相同，但通量法的潛水面是計算得出的，并非水頭法給定的規則“正弦波”。

圖3 通量上邊界條件的Tóth典型水流模式數值模擬結果Fig.3 Numerical simulation results of Tóth’s typical model with flux upper boundary

2.3 Tóth典型模式的數值模擬結果對比

用數值法進行Tóth典型模式對比模擬遵循的基本原則是保持盆地水均衡的一致性。Tóth典型模式滲流區修正前后，得出的水流模式相同，修正后盆地水均衡量略大于Tóth解析解（表1）；這是由于此時滲流區與潛水面保持一致，實際滲流區域就會有所增大，在上邊界定水頭求解時，根據達西定律水均衡量相應增大，圖2的補給強度較圖1條件增大了4.3%。表1中補給強度為各區補給量除以各區長度所得的單位長度補給量。

表1 Tóth典型模式的滲流區修正前后補給區補給強度對比Table 1 Recharge contrast of Tóth’s typical model with seepage zone modified before and after

Tóth典型模式的水頭法與通量法的數值模擬相比，在上邊界補給強度相同的條件下，兩者模擬的水流模式結果相同，即發育三級嵌套式水流系統（圖2與圖3對比）。并且通量法保持了5個補給區與水頭法的對應關系。利用通量法求解模擬時，潛水面不再固定，而是隨上邊界入滲條件而形成，模擬結果既符合實際，也能滿足改變盆地條件的水流系統模擬的要求。

3 Tóth模型變條件的通量上邊界對比模擬

3.1 變條件的水頭法模擬

Tóth在典型水流系統的基礎上，進行了改變盆地條件（表2）的模擬，如變斜率（A-1與 A-2）、變波幅（A-2與 A-3）以及改變盆地深度（A-3與 A-4）的模擬［2］。Tóth變條件的模擬結果如表2和圖4。

表2 Tóth定水頭方法變條件的地下水流系統模擬結果Table 2 Groundwater flow-systems simulation results of Tóth’s given head boundary with factor changed

根據變條件求解結果 Tóth［3］指出，地形起伏（對應地下水位波動）是形成嵌套式多級次水流模式的主要因素。當保持地下水位局部波幅不變，各谷底基線的區域性斜率增大時，多個排泄點相對高差增大，指向最低河谷的區域水流系統的范圍擴大，局部水流系統縮小退化為滯流帶或者消失（圖4a、b）；當區域性地形勢斜率不變，增大局部波幅時，局部水流系統將增大直至到達盆地底界，區域水流系統消失（圖4b、c）；其他條件不變，增加盆地深度（圖4c、d），水流下切深度增大，可形成中間和區域水流系統，呈現局部、中間及區域三級水流系統。

Tóth通過上述方法得出了盆地水流模式與“地形勢”和盆地深度的關系，從圖4和表3可知：1）水頭法的變條件模擬中，加大盆地上邊界斜率和地形起伏（波幅），實質都是在加大盆地的入滲補給，水流模式都會發生變化；2）定水頭條件，水流系統只能從局部系統轉化為多級系統；3）加大盆地深度，水流系統也能夠從局部系統轉化為多級系統，但上邊界入滲補給強度同時增加。也就是說，水頭法在改變其他任意參數時，存在盆地補給強度同時發生變化的問題。為克服這一模擬問題，下面利用通量法進行對比分析。

圖4 Tóth定水頭方法在變條件下地下水流數值模擬結果Fig.4 Groundwater flow-systems simulation results of Tóth＇s given head boundary with factor changed

3.2 變條件的通量法模擬

3.2.1 變斜率的模擬條件與結果

變斜率的通量法模擬采用圖4a計算的補給結果，在保持上邊界補給不變的條件下，設計3組逐漸增大的勢匯分布斜率，即調節5個勢匯分布高程，進行通量上邊界對比模擬，模擬結果如圖5和表3。

表3 變斜率的地下水流系統模擬比較Table 3 Simulation comparison of groundwater flow-system with the slope changed

2種方法僅從水流模式上看，都是由一級LS轉化成二級RS＋LS再轉化成三級RS＋MS＋LS，但水頭法在變斜率的同時，上邊界補給也發生了改變；而通量法在保持上邊界補給不變時，改變斜率即為潛水面發生變化，如果繼續加大勢匯分布斜率，部分高勢匯將失去排泄能力，如圖5c只有3個較低勢匯起作用。說明存在多個勢匯的盆地，除最低排泄區外，在沒有進行模擬計算前其他勢匯都只能稱為潛在勢匯，而采用定水頭上邊界求解時，固定了所有勢匯。

3.2.2 變波幅的模擬條件與結果

水頭法波幅加大，實質是上邊界補給加大，所以變波幅的通量法模擬參考圖4b條件采用3組變波幅的水頭法補給結果，進行通量法對比模擬，模擬結果如圖6和表4。

通量法與水頭法變波幅（圖6和表4）模擬結果表明：隨波幅加大，水流系統模式發生同規律的改變，但通量上邊界條件下，潛水面隨補給強度和勢匯位置自動形成；波幅減小到一定程度（相當于補給強度減小），部分高勢匯將失去排泄能力，如圖6a只有4個較低勢匯起作用，同樣也說明存在多個勢匯的盆地；除最低排泄區外，在沒有進行模擬計算前其他勢匯都只能稱為潛在勢匯；而采用水頭法求解時，固定了所有勢匯。

圖5 通量上邊界方法變斜率的地下水流系統模擬結果Fig.5 Groundwater flow-system simulation results of flux upper boundary with the slope changed

圖6 通量上邊界方法變波幅的地下水流系統模擬結果Fig.6 Groundwater flow-system simulation results of flux upper boundary with the amplitude changed

3.2.3 變深度的模擬條件與結果

變深度的通量法模擬采用圖4c計算的補給結果，保持上邊界補給不變，設計3組增大盆地深度，進行通量法對比模擬，模擬結果如圖7和表5。

從增大盆地深度模擬結果（圖7和表5）可以看出：2種方法模擬時，水流模式可能相同也可能不同；水頭法模擬，盆地深度增大的同時，盆地補給也增大；而通量法可以在保持上邊界補給不變條件下，單一因素地分析盆地深度變化對水流模式的影響。圖7結果表明，加大盆地深度，地下水流模式從簡單LS水流系統向三級RS＋MS＋LS水流系統方向發育，同時潛水面降低，水面坡度變緩。

表4 變波幅的地下水流系統模擬比較Table 4 Simulation comparison of groundwater flow-system with the amplitude changed

表5 變深度的地下水流系統模擬比較Table 5 Simulation comparison of groundwater flow-system with the depth changed

3.2.4 變滲透系數的模擬條件與結果

為了探討2種條件下滲透系數對水流系統模式的影響，筆者選取圖4b，保持上邊界補給不變，進行了一組改變盆地滲透系數的模擬，探討水頭法與通量法水流模式與滲透系數的關系。對比模擬結果如表6和圖8。

表6 變滲透系數的地下水流系統模擬比較

Table 6 Smulation comparison of groundwater flow-system with the hydraulic conductivity changed

注：E-2與表2中A-2的模擬條件相同。

圖7 通量上邊界變深度的地下水流系統模擬結果Fig.7 Groundwater flow-system simulation results of flux upper boundary with the depth changed

圖8 通量上邊界變滲透系數的地下水流系統模擬結果Fig.8 Groundwater flow-system simulation results of flux upper boundary with the hydraulic conductivity changed

從模擬結果（圖8）可以看出：盆地其他條件不變，不論滲透系數增大或減小，通量法的水流系統都會發生改變；而水頭法由于固定了上邊界水頭，根據達西滲流定律可知，模擬計算的補給強度隨滲透系數成比例增大（表6），兩者比值沒有變化，因此盆地水流模式不會發生改變［11］。所以，保持上邊界補給不變的通量法，揭示了水流系統變化特征，即隨滲透系數增大，潛水面將下降并趨于平緩，盆地部分較高勢匯將失去排泄能力（圖8c，a），即便是滲透系數極大的情況下，復雜盆地也將僅僅發育單級區域水流系統。

4 結論

1）Tóth利用規則的水頭上邊界條件，解析得出了典型地下水流模式。然而，水頭法模擬固化了盆地勢源與勢匯的位置和數目，在改變盆地條件的模擬中，盆地潛水面不能隨條件的改變而變化，這與實際不相符，也限制了地下水流模式的轉化研究；通量法能得出Tóth典型地下水流模式，且潛水面自動形成，在改變盆地條件模擬時，潛水面也隨條件而變化。

2）水頭法在改變盆地條件（上邊界水頭分布斜率、局部水頭變化波幅、盆地深度和滲透系數）的同時，上邊界補給也發生了改變，得出的地下水流模式并不是某種單一因素的影響結果；通量法能保持上邊界補給的一致，模擬分析單一因素對地下水流模式變化的影響。

3）基于通量法與水頭法在地下水流系統模擬中的優勢與不同，在進行盆地地下水流系統理論和實際研究時，應該綜合2種方法的特點，結合實際資料條件進行方法的選取與應用。

（

）：

［1］Hubbert M K.The Theory of Ground-Water Motion［J］.The Journal of Geology，1940，48（8）：785-944.

［2］Tóth J.A Theoretical Analysis of Groundwater Flow in Small Drainage Basins［J］.Journal of Geophysical Research，1963，68（16）：4795-4812.

［3］Tóth J.Mapping and Interpretation of Field Phenomena for Groundwater Reconnaissance in a Prairie Environment，Alberta，Canada［J］. Hydrological Sciences Journal，1966，11（2）：20-68.

［4］Tóth J.Groundwater Discharge：A Common Generator of Diverse Geologic and Morphologic Phenomena［J］.Hydrological Sciences Journal，1971，16（1）：7-24.

［5］Tóth J.Gravity-Induced Cross-Formation Water Flow：Possible Mechanism for Transport and Accumulation of Petroleum［J］.American Association of Petroleum Geologists Bulletin，1978，62（3）：567-568.

［6］Tóth J.Patterns of Dynamic Pressure Increment of Formation-Fluid Flow in Large Drainage Basins，Exemplified by the Red Earth Region，Alberta，Canada［J］.Bulletin of Canadian Petroleum Geology，1979，27（1）：63-86.

［7］Tóth J.Cross-Formational Gravity-Flow of Groundwater：A Mechanism of the Transport and Accumulation of Petroleum：The Generalized Hydraulic Theory of Petroleum Migration ［J］.Problems of Petroleum Migration，1980，29：121-167.

［8］張人權，梁杏，靳孟貴，等.當代水文地質學發展趨勢與對策［J］.水文地質工程地質，2005，32（1）：51-56.

Zhang Renquan，Liang Xing，Jin Menggui，et al.The Trends in Contemporary Hydrogeology ［J］.Hydrological and Engineering Geology，2005，32（1）：51-56.

［9］張人權，梁杏，靳孟貴，等.水文地質學基礎［M］.6版.北京：地質出版社，2011.

Zhang Renquan，Liang Xing，Jin Menggui，et al.Fundamentals of Hydrogeology［M］.6nd ed.Beijing：Geological Publishing House，2011.

［10］侯光才，林學鈺，蘇小四，等.鄂爾多斯白堊系盆地地下水系統研究［J］.吉林大學學報：地球科學版，2006，36（3）：391-398.

Hou Guangcai， Lin Xueyu， Su Xiaosi，et al.Groundwater System in Ordos Cretaceous Artisan Basin（CAB）［J］.Journal of Jilin University：Earth Science Edition，2006，36（3）：391-398.

［11］陳蓓蓓，宮輝力，李小娟.北京地下水系統演化與地面沉降過程［J］.吉林大學學報：地球科學版，2012，42（1）：373-379.

Chen Beibei，Gong Huili，Li Xiaojuan.Groundwater System Evolution and Land Subsidence Process in Beijing［J］.Journal of Jilin University：Earth Science Edition，2012，42（1）：373-379.

［12］劉宇，賈靜.基于 MATLAB小型潛水盆地地下水流動系統模擬研究［J］.地下水，2009（3）：1-3.

Liu Yu，Jia Jing.Study of Groundwater Flow Systems Simulation of Small Drainage Basin Based on Matlab［J］.Groundwater，2009（3）：1-3.

［13］梁杏，牛宏，張人權，等.盆地地下水流模式及其轉化與控制因素［J］.地球科學：中國地質大學學報，2012，37（2）：269-275.

Liang Xing，Niu Hong，Zhang Renquan，et al.Basinal Groundwater Flow Patterns and Their Transformation and Dominant Factors［J］.Earth Sciences：Journal of China University of Geosciences，2012，37（2）：269-275.

［14］Haitjema H M，MitchellB S.Are Water Tables a Subdued Replica of the Topography？［J］.Ground Water，2005，43（6）：781-786.

［15］劉彥，梁杏，權董杰，等.改變入滲強度的地下水流模式實驗［J］.地學前緣，2010，17（6）：111-116.

Liu Yan， Liang Xing， Quan Dongjie，et al.Experiments of Groundwater Flow Patterns Under Changes of Infiltration Intensity［J］.Earth Science Frontiers，2010，17（6）：111-116.

［16］權董杰.二維盆地地下水流模式與轉化規律分析［D］.武漢：中國地質大學，2011.

Quan Dongjie. Analysis of Two-Dimension Groundwater Flow Patterns and Transformation Rules［D］.Wuhan：China university of geosciences，2011.

［17］Liang Xing，Quan Dongjie，Jin Menggui，et al.Numerical Simulation of Groundwater Flow Patterns Using Flux as Upper Boundary［J］.Hydrological Processes，2012.

［18］楊春玲，邢世錄，韓愛中.地下水系統數值模擬的研究進展［J］.內蒙古科技與經濟，2007（21）：305-307.

Yang Chunling，Xing Shilu，Han Aizhong.The Research Progress of Groundwater System Numerical Simulation［J］.Inner Mongolia Technology and Economy.2007（21）：305-207.