大區(qū)域含水層中局部復雜水流過程的多尺度模擬

曾季才，查元源，楊金忠

（武漢大學 水資源與水電工程科學國家重點實驗室，湖北 武漢 430072）

1 研究背景

長期以來，大區(qū)域地下水流運動數(shù)值模擬受到計算精度與成本的相互制約。在區(qū)域尺度問題的模擬中，往往通過尺度提升來將點尺度實測信息高度平滑地離散到大時間和空間網(wǎng)格上，進而與區(qū)域尺度模型進行尺度匹配［1-3］，本質上放棄了對局部小尺度問題的精細刻畫。而保留全區(qū)域大尺度過程描述，并將研究重點放在具有小尺度水流運動特征的局部區(qū)域，符合許多實際問題的求解需求，因此產(chǎn)生了一系列在區(qū)域尺度和局部尺度上實現(xiàn)多尺度水流運動模擬的數(shù)值模型。這些模型往往借助非均勻網(wǎng)格離散方法來求解地下水方程，例如基于有限單元法的地下水模型求解軟件（如FEFLOW［4］和S3D［5］），它們主要采用直接局部加密方法對全區(qū)域內不同分辨率的網(wǎng)格節(jié)點進行統(tǒng)一求解，容易在小尺度模擬范圍外增加大量不必要的漸變節(jié)點，同時其計算效率和精度易受非均勻網(wǎng)格剖分質量的顯著影響，從而增大實際應用的難度。此外在溶質運移等非對稱矩陣問題的求解中，這類方法的網(wǎng)格加密速度和加密極限受算法收斂性和跨尺度網(wǎng)格混疊誤差的影響［6］；與此同時，基于有限差分法的地下水模型求解軟件（如MODFLOW-USG［7］），則因為非結構化的網(wǎng)格容易產(chǎn)生不規(guī)則的稀疏帶狀矩陣以及條件對角占優(yōu)問題，增大了應用于大區(qū)域尺度時的矩陣求解難度；此外，基于有限差分局部網(wǎng)格加密的多網(wǎng)格耦合模型（如MODTMR［8-9］和MODFLOW LGR［10-11］），主要通過連接不同分辨率的獨立網(wǎng)格節(jié)點來求解不同尺度的地下水流運動問題，這類模型方法在處理小尺度復雜邊界條件及源匯項等問題上的靈活性及刻畫精度受到有限差分網(wǎng)格幾何形態(tài)的限制［11］，目前尚不能實現(xiàn)大區(qū)域含水層中的多時空尺度水流過程的高效模擬［8-10］，且用不同網(wǎng)格之間進行迭代耦合（如MODFLOW-LGR［10］）來求解高度非均質含水層等存在復雜水力條件的問題時，其計算效率和精度并不能占絕對優(yōu)勢［12］。

本文提出一種多時空尺度耦合模型，它采用尺度分離策略，在獲取區(qū)域大尺度模型粗網(wǎng)格解的同時，精細刻畫局部小尺度水流運動過程，將不同尺度模型的時間和空間信息進行關聯(lián)和傳遞。即在大尺度含水層構建粗網(wǎng)格母模型，在局部小尺度構建三維精細網(wǎng)格子模型，采用自由非匹配網(wǎng)格及自適應時間步長加密的非迭代耦合方案，進行多時空尺度的子母模型耦合，從而實現(xiàn)對復雜小尺度地下水流運動過程更高效和更精細的模擬。本文模型將不同尺度模型的參數(shù)輸入、邊界條件、初始條件和源匯項等信息進行尺度分離，通過模型邊界信息傳遞來耦合不同時空尺度的模型。其中，大區(qū)域尺度問題求解采用有限差分單元網(wǎng)格地下水模型MODFLOW 2005［13］，該模型對大區(qū)域地下水模擬問題有較高的穩(wěn)定性和求解效率［14］并已有較大量的代碼維護和改進［15-17］；而局部小尺度問題求解則采用控制體積有限單元模型S3D［5］，該模型能較好地適應小尺度高密度不規(guī)則網(wǎng)格的求解，并能適應局部大水力梯度、復雜邊界條件、非均質傾斜含水層等特殊問題［5］。區(qū)別于傳統(tǒng)多網(wǎng)格耦合模擬軟件（如MODFLOW-LGR［10］及MODTMR［8］），本文模型為獨立開發(fā)，并側重體現(xiàn)小尺度模型在求解局部復雜水流運動過程時所具有的優(yōu)勢。在后期的研究中，該模型已衍生出飽和三維溶質運移模型［5］、擬三維飽和非飽和水流運動及溶質運移模型［18］，以及三維飽和非飽和水流運動及溶質運移模型［19］，對于本文模型的持續(xù)開發(fā)具有較好的延續(xù)性及擴展性。

2 三維地下水多尺度耦合模型

2.1 多尺度模型耦合框架 本文模型采用的地下水流運動控制方程及邊界條件可以歸納如下：

本文模型對不同尺度的概念模型進行不同方法的數(shù)值離散，其中區(qū)域尺度模型為有限差分網(wǎng)格［13］，而局部尺度模型則為控制體積有限單元網(wǎng)格［5］，如圖1所示。在空間上，兩種空間尺度的網(wǎng)格可以通過各自的網(wǎng)格密度、分層密度和邊界形狀進行獨立創(chuàng)建。對于本文區(qū)域大尺度模型而言，其繼承了MODFLOW 2005軟件的所有子程序包，能有效處理區(qū)域尺度地下水流的入滲、蒸發(fā)、排水、抽水等問題［13］，而本文局部尺度模型S3D［5］則側重于刻畫高密度網(wǎng)格上具有小尺度特征的地下水流運動過程。模型為非迭代耦合，基本框架分為4部分（見圖2）：（1）求解區(qū)域尺度粗網(wǎng)格模型，并得到全局節(jié)點水頭；（2）通過三維水頭插值得到局部尺度模型在大時間尺度節(jié)點上的邊界節(jié)點水頭；（3）通過自適應時間步長調整及二階時間插值得到任意時間節(jié)點上的小尺度模型邊界水頭；（4）求解小尺度模型，得到局部小尺度模型高精度解。相比MODFLOW-LGR［10］的迭代耦合方案，本文模型的非迭代耦合方案雖然損失了一定的計算精度，但由于不存在小尺度解的反饋機制，在大區(qū)域復雜水流過程及非均質問題中具有更高的計算效率［12］。

2.2 多空間尺度信息傳遞 在多空間尺度耦合問題中，本文模型采用三維Darcy水頭插值方法來實現(xiàn)區(qū)域尺度模型水頭信息向局部尺度模型邊界節(jié)點的傳遞，該方法基于二維平面水頭Darcy-Planar插值公式［10］，主要通過地下水流運動的水頭和通量本構關系——Darcy定律，來求解任意小尺度模型邊界節(jié)點相對于大尺度模型節(jié)點水頭所存在的水頭差，從而得到其插值水頭［20-21］。相比傳統(tǒng)線性空間插值，三維Darcy水頭插值方法代碼實現(xiàn)簡單，且更多地考慮了節(jié)點水力參數(shù)的空間變異性，因而更適用于有限差分模型中任意空間位置節(jié)點的水頭插值。設任意小尺度模型邊界節(jié)點為C，其最臨近的大尺度模型節(jié)點為P0。P0在j=xyz三個軸向相鄰節(jié)點之間的水流通量qpj滿足達西定律：

圖1 多子域自由網(wǎng)格匹配的多尺度耦合模型示意圖

圖2 本文模型流程圖

式中：Cj為P0在j=x，y，z三個軸向上與相鄰節(jié)點之間的導水度（定義及公式見文獻［13］），m2/d；H0為大區(qū)域尺度模型節(jié)點P0的水頭解，m；H0′為j軸向上與P0相鄰的大區(qū)域尺度模型節(jié)點Pi（i=1，2，…，6）的水頭解，m。同時，節(jié)點P0控制范圍內的j方向流量qc j可以表示為：

式中：Kj為節(jié)點P0在j=x，y，z三個軸向上的水力傳導度，m/d；Aj為節(jié)點P0在j=x，y，z三個軸向上的過水橫截面積，m2；δlj為節(jié)點C與P0在j=x，y，z方向上的坐標差值，m；δhj為節(jié)點C與P0在j=x，y，z方向上的水頭差值。

本文插值方法假設相鄰母網(wǎng)格之間的水流通量與單個母網(wǎng)格內部的水流通量相等，即qpj=qcj，則由式（2）和式（3）可以得到

將j=x，y，z三個軸向上的δhj累加到P0節(jié)點的水頭H0，得到節(jié)點C的插值水頭hC：

2.3 多時間尺度信息傳遞 在多時間尺度耦合問題中，本文模型采用自適應時間步方法［22］獲取小尺度模型的時間步長，用以獲取自適應時間步長加密條件下的小尺度模型邊界信息。該方法主要依據(jù)當前求解結果的水頭變化值來調整子模型時間步長：

式中：Δt′為上一次完成求解所采納的時間步長，d；Δhe為控制參數(shù)，用于限定某個時間步長內子模型節(jié)點所容許的水頭差，m；為當前子模型解h1與上一時刻有效解h0的最大絕對差值，m；Δt為更新的時間步長，d。式可用于評估當前所使用時間步長Δt′的合理性，當ε≥ 1時認為時間步長Δt′合理，則Δt用于估算下一步時間步長；當ε＜1時，認為Δt′不合理，則采納Δt重新求解子模型。

求解子模型前，通過二階插值方法得到任意t時刻的小尺度模型一類邊界條件（圖3），該二階曲線擬合公式如下：

圖3 基于邊界水頭插值的多時間和多空間尺度模型耦合

3 模型驗證及分析

為驗證本文模型在多空間尺度模擬、多時間尺度模擬和復雜小尺度水流運動過程模擬上具有的優(yōu)勢，設計了均質含水層中的2個數(shù)值驗證算例，包括（1）多空間尺度大水力梯度模擬，和（2）多時空尺度復雜局部水流運動過程模擬。本文模型的驗證算例建立在一個長×寬×高為1 000 m×1 000 m×50 m的均質潛水含水層內（水力傳導度ks=10 m/d，給水度Sy=0.068 6，彈性釋水系數(shù)μs=1.0×10-5m-1），為避免相同邊界條件下不同網(wǎng)格密度對模型解的影響［23］，將四周和底板設為隔水邊界，含水層初始水頭為42 m。算例在子域1內（見圖4中矩形點劃線方框范圍）以I=300 mm/d進行持續(xù)均勻灌溉。在含水層獲取了具有全區(qū)域大尺度補給系數(shù)隨機場（均值為=0.167，相關長度為100 m，標準差為0.01）。

圖4 多空間尺度模型參數(shù)（補給系數(shù)為例）獲取示意圖

3.1 多空間尺度大水力梯度模擬 本算例布設了兩個分別位于（x，y）=（500 m，500 m）和（550 m，500 m），濾水段均為zbottom～ztop=20～30 m的獨立抽水井，各以定流量500 m3/d和250 m3/d進行為期60 d的持續(xù)抽水。采用克里金插值方法得到如圖4子域1內的復雜不規(guī)則邊界形態(tài)的補給系數(shù)分區(qū)。算例將全區(qū)域劃分為均勻有限差分粗網(wǎng)格（節(jié)點間距Δx=Δy=22.222 m），將局部尺度上的子域1采用高密度非結構化有限單元三角網(wǎng)（節(jié)點間距Δl≈2.5 m），等效的局部網(wǎng)格水平向加密比為8.9∶1，垂向加密比為1∶1。算例以通用地下水數(shù)值模擬軟件MODFLOW［13］在高密度網(wǎng)格（節(jié)點間距Δx=Δy=2.5 m）和最小時間步長（Δt=0.1 d）下的計算結果作為高精度參照解。

圖5給出了子域1在持續(xù)60 d大流量抽水條件下的末時刻水頭剖面。其中，本文模型小尺度解在抽水井附近節(jié)點出現(xiàn)極大水頭誤差（0.094 m），子域1內平均水頭絕對誤差為0.003 9 m，其中92%的節(jié)點水頭絕對誤差低于0.005 m；同時，本文模型大尺度解與精細網(wǎng)格參照解的絕對水頭誤差在抽水井附近達到最大（約1.74 m）。圖5中虛線代表了粗網(wǎng)格大尺度模型解，兩個抽水井產(chǎn)生的漏斗降深、形狀和方位均發(fā)生失真，這是由于粗網(wǎng)格區(qū)域尺度模型解難以用于刻畫局部尺度細節(jié)，其在小尺度區(qū)域的指示精度不再作為參考。由于小尺度模型解的誤差較小，可以認為這粗網(wǎng)格處理方式在局部大水力梯度問題中對全局大尺度解的精度影響不大。

對于大流量持續(xù)抽水問題而言，末時刻（t=60 d）誤差將累積到最大值。統(tǒng)計顯示，該時刻小尺度模型邊界插值節(jié)點水頭誤差均值為0.001 8 m，極大值為0.008 m，且90%的插值點水頭誤差低于0.003 m，這些誤差主要來自大尺度模型解的自身誤差，并通過時間和空間的水頭插值向小尺度模型節(jié)點傳遞并累積。小尺度模型解所獲取的誤差傳遞相比大尺度模型解在子域1范圍內的巨大誤差僅占比0.46%，因而可以忽略。同時從圖5可以看出，子域1范圍內的粗網(wǎng)格節(jié)點源匯項（井2）輸入存在不可避免的空間錯位，這種錯位對子域1外的大區(qū)域尺度解影響較小，但在子域1范圍內影響較大，證明了對于局部尺度問題采用小尺度模型求解的必要性。對比子域1內小尺度模型解與MODFLOW高密度時空網(wǎng)格下的參照解等值線圖（見圖6），可見本文模型在大區(qū)域尺度的局部大水力梯度問題的模擬中能實現(xiàn)較高精度。

圖5 對比子域1內本文模型大尺度和小尺度解與參照解（y=500 m處的x軸向剖面）

圖6 對比末時刻子域1內本文模型小尺度解（虛線）與參照解（實線）

3.2 多時空尺度局部復雜水流過程模擬 在3.1節(jié)中驗證算例的基礎上，針對復雜時變抽水井問題，在井1處設置了如圖7所示帶有隨機擾動（6個應力期內的抽水量均值為500、400、600、400、550和350 m3/d，標準差分別為28.9和10.7 m3/d）的抽水量時間序列，各應力期持續(xù)10 d不間斷抽水；井2處抽水量設置為井1的0.5倍。其中，本文模型在局部尺度上的模擬則采用自適應時間步長，為子域1中井1和井2施加的隨機擾動來表征實際抽水量的劇烈波動。此外，在子域1內，仍然以算例1中相同的入滲補給系數(shù)隨機場進行相同灌溉量下的人工補給。在如圖8所示的子域2內構建變密度非結構化Delaunay三角網(wǎng)。模型垂向模擬范圍位于局部含水層（z=20～50 m高程范圍內），并根據(jù)源匯項求解需要將子域2剖分為8層（見圖8），該子網(wǎng)格無需與母網(wǎng)格進行嚴格的空間匹配，底部節(jié)點高程無需與區(qū)域尺度模型底板高程匹配。算例在子域2內設置如圖8所示6個定流量（Q=－50 m3/d）抽水井（編號為3，4，…，8），濾水段高程均為zbottom～ztop=30～40 m。

圖7 子域1內井1和井2抽水流量均值和隨機擾動實值時間序列

圖8 子域2模型非匹配網(wǎng)格示意圖

算例在井1和井2上施加劇烈擾動的變流量抽水井，并設置觀測點A1（x=500 m，y=500 m，z=25 m）和A2（x=525 m，y=500 m，z=25 m），得到如圖9（a）和（b）中觀測水頭隨時間的變化。圖9中，觀測點的小尺度模型解均能較好地與參照解進行匹配，并能精細刻畫局部水頭變化過程；而點A1的大尺度模型解則由于施加了大時間尺度應力項而產(chǎn)生顯著誤差（圖9（a）），點A2的大尺度模型解則由于稀疏網(wǎng)格解自身誤差而與參照解之間存在大幅度偏離（圖9（b））。針對本文小尺度模型在子域1內的求解精度，統(tǒng)計z=25 m水平剖面上所有節(jié)點的絕對水頭誤差隨時間的變化，得到小尺度模型水頭誤差極大值保持在0.026 m以下（位于抽水井附近），而模擬期內所有節(jié)點平均絕對水頭誤差為0.005 6 m，大尺度模型“傳遞誤差”均值為0.004 2 m，可以認為本文模型對于多時間尺度、大水力梯度問題的模擬具有較高的計算精度。統(tǒng)計子域1的小尺度模型邊界插值節(jié)點的絕對水頭誤差，均低于1.0×10-4m，且均值約為2.0×10-5m，該均值約占粗網(wǎng)格大尺度模型“傳遞誤差”均值的1.1%，可見本文時間和空間尺度水頭插值方法能獲取較高精度的小尺度模型邊界水頭。

圖9 對比在觀測點A1和A2處的本文模型解粗網(wǎng)格大尺度解、高密度網(wǎng)格小尺度解與參照解隨時間的變化

圖10（a）細虛線為全區(qū)域內粗網(wǎng)格大尺度模型解；實線為高密度網(wǎng)格參照解；粗點劃線為高密度網(wǎng)格小尺度模型解。統(tǒng)計得到粗網(wǎng)格模型解在全區(qū)域的絕對水頭誤差均值為0.003 1 m，可見對于存在局部時變大水力梯度的復雜流場之外的大尺度區(qū)域，本文模型粗網(wǎng)格解能滿足大部分范圍內的精度，而局部區(qū)域內（子域1和子域2）則需要采用小尺度模型進行刻畫。在復雜多時空尺度、多子域、非匹配自由網(wǎng)格模擬條件下，圖10（b）對比了本文模型在子域2中的小尺度解與參照解，其絕對水頭誤差均值為0.004 m，極大值為0.028 m，其中95%的子網(wǎng)格節(jié)點絕對水頭誤差低于0.01 m，對于大水力梯度條件下的區(qū)域尺度模擬而言計算精度較高。統(tǒng)計末時刻子域1內本文模型小尺度解的絕對誤差，當單獨運行子域1模型時，其均值為0.005 6 m；當同時運行子域1和2模型時，其均值為0.006 3 m。子域1內的誤差出現(xiàn)明顯增長，主要是由于子域1臨近的子域2區(qū)域出現(xiàn)持續(xù)大水力梯度干擾，對母網(wǎng)格在子域1邊界節(jié)點處的求解精度有微弱影響，因此在實際應用中，應盡量采用單個高密度子網(wǎng)格求解（覆蓋）所有臨近大水力梯度問題。

圖10 對比本文模型解（等高線）與參照解

表1 本文模型計算成本及誤差統(tǒng)計（對比參照模型）

3.3 優(yōu)勢與不足 本文模型繼承了MODFLOW［13］作為區(qū)域尺度地下水模擬軟件所具有的先天優(yōu)勢，將其應用于粗網(wǎng)格大區(qū)域尺度模型求解，能獲得較高的計算效率和精度；而本文小尺度模型S3D則作為一種能精細刻畫局部水流運動過程的小尺度模型被提出，其非結構化網(wǎng)格能極大程度上滿足小尺度空間離散的需求，此外自適應時間步長方法能高效處理小尺度大水力梯度及劇烈時變源匯項。與此同時，由于S3D模型尚處于開發(fā)階段，計算效率有較大的優(yōu)化和提升的空間?？紤]到不同數(shù)值算法在數(shù)據(jù)封裝和矩陣求解方法上的差異，本文耦合模型的計算效率難以簡單通過計算耗時進行評價。表1對比了各算例中本文模型與參照模型的時間步數(shù)、節(jié)點總數(shù)、計算耗時和平均絕對水頭誤差的統(tǒng)計結果。相比傳統(tǒng)單一尺度高密度網(wǎng)格模型，本文模型能將矩陣規(guī)模和求解工作量減少90%以上，并在一定程度上降低計算耗時，同時在局部尺度復雜水流運動過程的模擬中保持較高的求解精度。

4 結論

針對傳統(tǒng)單一尺度模型在模擬大區(qū)域含水層中的小尺度復雜水流運動過程上的不足，本文建立了一種能進行多時空尺度非匹配自由網(wǎng)格模擬的多網(wǎng)格耦合模型。通過一系列數(shù)值試驗，驗證了該模型在帶有小尺度特征的大水力梯度問題的模擬上具有良好的適應性和計算精度。相比傳統(tǒng)單一尺度模型方法，本文模型在區(qū)域大尺度和局部小尺度水流問題的求解上進行了空間和時間尺度的分離，有利于不同尺度模型的穩(wěn)定性和收斂性。而非匹配自由時空網(wǎng)格耦合方法的應用，則極大地保證了本文模型對局部尺度問題中復雜邊界條件及源匯項和非均質含水層參數(shù)的高效數(shù)值離散，且子域模型網(wǎng)格可自由劃分及加密，能較大程度上降低網(wǎng)格節(jié)點數(shù)量并提升局部問題的求解精度。

為了更準確地刻畫大區(qū)域含水層中的小尺度復雜水流運動過程，本文模型需要對模型中不同水力參數(shù)、邊界條件、初始條件和源匯項進行尺度匹配，進一步的研究可將升尺度和降尺度方法聯(lián)合運用到不同尺度模型之間的信息交互中，以提高模型在長時間序列模擬中的穩(wěn)定性和可靠性。此外，由于小尺度模型中需要考慮非達西流的近似解析解［24］，目前S3D模型尚不能精確模擬抽水井附近小范圍內的井流。通過本文模型來精細刻畫局部尺度復雜水流運動過程，對于模擬小尺度復雜的溶質運移、飽和-非飽和水量交換、作物根系吸水等過程具有重要意義。

參 考 文 獻：

［1］ 張祥偉，竹內邦良.大區(qū)域地下水模擬的理論和方法［J］.水利學報，2004（6）：7-13.

［2］ DURLOFSKY L J，JONES R C，MILLIKEN W J.A nonuniform coarsening approach for the scale-up of displace?ment processes in heterogeneous porous media［J］.Advances in Water Resources，1997，20（5）：335-347.

［3］ KARIMI-FARD M，DURLOFSKY L J.Accurate resolution of near-well effects in upscaled models using flowbased unstructured local grid refinement［J］.Spe Journal，2012，17（4）：1084-1095.

［4］ TREFRY M G ，MUFFELS C.Feflow：A finite-element ground water flow and transport modeling tool［J］.Ground Water，2007，45（5）：525-528.

［5］ LIN L，YANG J，ZHANG B，et al.A simplified numerical model of 3-D groundwater and solute transport at large scale area［J］.Journal of Hydrodynamics，2010，22（3）：319-328.

［6］ ODMAN M T，RUSSELL A G.On local finite-element refinements in multiscale air-quality modeling［J］.Envi?ronmental Software，1994，9（1）：61-66.

［7］ PANDAY S，LANGEVIN C D，NISWONGER R G，et al.MODFLOW-USG version 1：An unstructured grid ver?sion of MODFLOW for simulating groundwater flow and tightly coupled processes using a control volume finite-dif?ference formulation［R］.US Geological Survey，2013.

［8］ LEAKE S A，CLAAR D V.Procedures and computer programs for telescopic mesh refinement using MODFLOW［R］.Center for Integrated Data Analytics Wisconsin Science Center，1999.

［9］ BOOTH C J，BREUERE P.Using MODFLOW with TMR to Model Hydrologic Effects and Recovery in the Shal?low Aquifer System Above Longwall Coal Mining［C］//10th IMWA Congress：Mine Water and the Environment.Karlocy Vary，Czech Republic，2008.

［10］ MEHL S W，Hill M C.MODFLOW-LGR—Documentation of Ghost Node Local Grid Refinement（LGR2）for Multiple Areas and the Boundary Flow and Head（BFH2）Package［R］.Section A Ground Water in Book Model?ing Techniques，2013.

［11］ 楊楊，邵景力等，基于LGR的不規(guī)則MODFLOW模型嵌套技術研究［J］.水文地質工程地質，2015（1）：22-28.

［12］ MEHL S，HILL M C.Development and evaluation of a local grid refinement method for block-centered finite-dif?ference groundwater models using shared nodes［J］.Advances in Water Resources，2002，25（5）：497-511.

［13］ HARBAUGH A W.MODFLOW-2005，the US Geological Survey modular ground-water model：The ground-wa?ter flow process.2005：US Department of the Interior［R］.US Geological Survey Reston，VA，USA，2005.

［14］ 陳皓銳，高占義，王少麗，等.基于MODFLOW的潛水位對氣候變化和人類活動改變的響應［J］.水利學報，2012，43（3）：344-353.

［15］ VERMEULEN P T M，HEEMINK A W，ESTROET C B M T.Reduction of large-scale numerical ground water flow models［C］//Computational Methods in Water Resources Proceedings.Elsevier Science Bv：Amsterdam.2002，397-404.

［16］ 劉路廣，崔遠來.灌區(qū)地表水-地下水耦合模型的構建［J］.水利學報，2012，43（7）：826-833.

［17］ JEPPESEN J，CHRISTENSEN S.A MODFLOW infiltration device package for simulating storm water infiltration［J］.Groundwater，2015，53（4）：542-549.

［18］ ZHU Y，SHI L，LIN L，et al.A fully coupled numerical modeling for regional unsaturated-saturated water flow［J］.Journal of Hydrology，2012，475：188-203.

［19］ ZHA Y，SHI L，YE M，et al.A generalized Ross method for two-and three-dimensional variably saturated flow［J］.Advances in Water Resources，2013，54：67-77.

［20］ WASSERMAN M.Local grid refinement for three-dimensional simulators［C］//SPE Symposium on Reservoir Sim?ulation.Society of Petroleum Engineers.1987.

［21］ HAEFNER F，Boy S.Fast transport simulation with an adaptive grid refinement［J］.Ground Water，2003，41（2）：273-279.

［22］ FORSYTH P，SAMMON P.Practical considerations for adaptive implicit methods in reservoir simulation［J］.Journal of Computational Physics.1986，62：265-81.

［23］ MEHL S，Hill M C.Grid-size dependence of Cauchy boundary conditions used to simulate stream-aquifer inter?actions［J］.Advances in Water Resources，2010，33（4）：430-442.

［24］ 文章，黃冠華，劉壯添，等.越流含水層中抽水井附近非達西流兩區(qū)模型近似解析解［J］.地球科學：中國地質大學學報，2011，36（6）：1165-1172.