牛 宏，梁 杏，2，張人權(quán)

1.中國(guó)地質(zhì)大學(xué)環(huán)境學(xué)院，武漢 430074

2.中國(guó)地質(zhì)大學(xué)生物地質(zhì)與環(huán)境地質(zhì)國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢 430074

0 引言

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

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

水頭法在進(jìn)行水流系統(tǒng)研究中存在的問題究竟如何，與通量法偏差有多大，應(yīng)用時(shí)限制條件如何，還沒有定量的對(duì)比模擬分析報(bào)道。筆者旨在用數(shù)值法進(jìn)行水頭法和通量法的對(duì)比模擬，通過模擬結(jié)果的定量比較，具體指出水頭法模擬的問題與局限性。

1 Tóth典型地下水流系統(tǒng)求解簡(jiǎn)介

Tóth分別采用一次函數(shù)和一次函數(shù)疊加正弦波函數(shù)，解析出簡(jiǎn)單與復(fù)雜盆地地下水流系統(tǒng)模式［2］。其中，采用一次函數(shù)疊加正弦波函數(shù)的水頭法進(jìn)行盆地二維流模型的數(shù)學(xué)方程如下：

式中：h為地下水水頭；x、z分別代表以谷底左下方為起點(diǎn)到盆地地下水流動(dòng)區(qū)域中任何一點(diǎn)的水平和垂直距離；Kx、Kz分別為x和z方向的滲透系數(shù)；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，計(jì)算出存在5個(gè)高地勢(shì)和5個(gè)低地勢(shì)的波狀上邊界；繪制了盆地多級(jí)地下水流系統(tǒng)圖（圖1）；認(rèn)為在地形復(fù)雜的均質(zhì)各向同性潛水盆地中，地下水流呈現(xiàn)多級(jí)次特點(diǎn)，并分別定義為局部的（local system，LS）、中間的 （mittle system，MS）及區(qū)域的 （regional system，RS）地下水流系統(tǒng)。

圖1 地形復(fù)雜的均質(zhì)各向同性盆地中的多級(jí)次水流系統(tǒng)Fig.1 Multilevel flow systems of homogeneous and isotropic basin in complex terrain

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

2 Tóth典型地下水流系統(tǒng)的數(shù)值模擬

2.1 水頭法的數(shù)值模擬

應(yīng)用數(shù)值模擬軟件采用加拿大Waterloo水文地質(zhì)公司在Modflow的基礎(chǔ)上應(yīng)用現(xiàn)代可視化技術(shù)開發(fā)研制的Visual-Modflow，它是目前國(guó)際一致認(rèn)可的三維地下水流和溶質(zhì)運(yùn)移模擬的標(biāo)準(zhǔn)可視化專業(yè)軟件系統(tǒng)［18］。

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

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

2.2 通量法的數(shù)值模擬

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

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

2.3 Tóth典型模式的數(shù)值模擬結(jié)果對(duì)比

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

表1 Tóth典型模式的滲流區(qū)修正前后補(bǔ)給區(qū)補(bǔ)給強(qiáng)度對(duì)比Table 1 Recharge contrast of Tóth’s typical model with seepage zone modified before and after

Tóth典型模式的水頭法與通量法的數(shù)值模擬相比，在上邊界補(bǔ)給強(qiáng)度相同的條件下，兩者模擬的水流模式結(jié)果相同，即發(fā)育三級(jí)嵌套式水流系統(tǒng)（圖2與圖3對(duì)比）。并且通量法保持了5個(gè)補(bǔ)給區(qū)與水頭法的對(duì)應(yīng)關(guān)系。利用通量法求解模擬時(shí)，潛水面不再固定，而是隨上邊界入滲條件而形成，模擬結(jié)果既符合實(shí)際，也能滿足改變盆地條件的水流系統(tǒng)模擬的要求。

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

3.1 變條件的水頭法模擬

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

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

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

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

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

3.2 變條件的通量法模擬

3.2.1 變斜率的模擬條件與結(jié)果

變斜率的通量法模擬采用圖4a計(jì)算的補(bǔ)給結(jié)果，在保持上邊界補(bǔ)給不變的條件下，設(shè)計(jì)3組逐漸增大的勢(shì)匯分布斜率，即調(diào)節(jié)5個(gè)勢(shì)匯分布高程，進(jìn)行通量上邊界對(duì)比模擬，模擬結(jié)果如圖5和表3。

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

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

3.2.2 變波幅的模擬條件與結(jié)果

水頭法波幅加大，實(shí)質(zhì)是上邊界補(bǔ)給加大，所以變波幅的通量法模擬參考圖4b條件采用3組變波幅的水頭法補(bǔ)給結(jié)果，進(jìn)行通量法對(duì)比模擬，模擬結(jié)果如圖6和表4。

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

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

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

3.2.3 變深度的模擬條件與結(jié)果

變深度的通量法模擬采用圖4c計(jì)算的補(bǔ)給結(jié)果，保持上邊界補(bǔ)給不變，設(shè)計(jì)3組增大盆地深度，進(jìn)行通量法對(duì)比模擬，模擬結(jié)果如圖7和表5。

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

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

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

3.2.4 變滲透系數(shù)的模擬條件與結(jié)果

為了探討2種條件下滲透系數(shù)對(duì)水流系統(tǒng)模式的影響，筆者選取圖4b，保持上邊界補(bǔ)給不變，進(jìn)行了一組改變盆地滲透系數(shù)的模擬，探討水頭法與通量法水流模式與滲透系數(shù)的關(guān)系。對(duì)比模擬結(jié)果如表6和圖8。

表6 變滲透系數(shù)的地下水流系統(tǒng)模擬比較

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

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

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

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

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

4 結(jié)論

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

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

3）基于通量法與水頭法在地下水流系統(tǒng)模擬中的優(yōu)勢(shì)與不同，在進(jìn)行盆地地下水流系統(tǒng)理論和實(shí)際研究時(shí)，應(yīng)該綜合2種方法的特點(diǎn)，結(jié)合實(shí)際資料條件進(jìn)行方法的選取與應(yīng)用。

（

）：

［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］張人權(quán)，梁杏，靳孟貴，等.當(dāng)代水文地質(zhì)學(xué)發(fā)展趨勢(shì)與對(duì)策［J］.水文地質(zhì)工程地質(zhì)，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］張人權(quán)，梁杏，靳孟貴，等.水文地質(zhì)學(xué)基礎(chǔ)［M］.6版.北京：地質(zhì)出版社，2011.

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

［10］侯光才，林學(xué)鈺，蘇小四，等.鄂爾多斯白堊系盆地地下水系統(tǒng)研究［J］.吉林大學(xué)學(xué)報(bào)：地球科學(xué)版，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］陳蓓蓓，宮輝力，李小娟.北京地下水系統(tǒng)演化與地面沉降過程［J］.吉林大學(xué)學(xué)報(bào)：地球科學(xué)版，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小型潛水盆地地下水流動(dòng)系統(tǒng)模擬研究［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］梁杏，牛宏，張人權(quán)，等.盆地地下水流模式及其轉(zhuǎn)化與控制因素［J］.地球科學(xué)：中國(guó)地質(zhì)大學(xué)學(xué)報(bào)，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］劉彥，梁杏，權(quán)董杰，等.改變?nèi)霛B強(qiáng)度的地下水流模式實(shí)驗(yàn)［J］.地學(xué)前緣，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］權(quán)董杰.二維盆地地下水流模式與轉(zhuǎn)化規(guī)律分析［D］.武漢：中國(guó)地質(zhì)大學(xué)，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］楊春玲，邢世錄，韓愛中.地下水系統(tǒng)數(shù)值模擬的研究進(jìn)展［J］.內(nèi)蒙古科技與經(jīng)濟(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.