楊 海，姜月華，周權(quán)平，楊 輝，劉 林

(中國地質(zhì)調(diào)查局南京地質(zhì)調(diào)查中心，江蘇 南京 210016)

土體中的多維飽和-非飽和流模型被廣泛用于入滲過程中的水流運動模擬[1]和土體中的滲流問題[2-3]等，是非常重要的土體水分定量化計算工具[4-5]。早期建立的多維飽和-非飽和流方程稱為改進的Richards方程，是基于單一水頭(h)的方程[6]。該形式在非飽和流計算中有其局限性[7-8]，因滲流參數(shù)在時段內(nèi)受毛管水頭變化影響劇烈，易造成較大的質(zhì)量守恒誤差和對下滲深度的錯誤估算[9-10]。而利用水頭(h)和土壤含水量(θ)構(gòu)建的混合式方程在求解中可有效保證質(zhì)量守恒[9]，因此成為最通用的多維飽和-非飽和流方程。在眾多求解該偏微分方程的方法中，有限差分法因其簡易方便而被廣泛使用。

在以往的飽和-非飽和流模型中，每個離散單元需申請維度內(nèi)最大可能連接單元數(shù)的系數(shù)存儲空間[11](如三維中為6個)，而后根據(jù)單元與相鄰單元間的連接情況計算各系數(shù)。當(dāng)相鄰單元為無效單元或相互間無水力聯(lián)系時，系數(shù)為0。因此，該方法中難免會存儲部分零系數(shù)，且零系數(shù)的規(guī)模隨著邊界單元、內(nèi)部水力聯(lián)系復(fù)雜度的增加而增加，這將極大影響模型在大規(guī)模網(wǎng)格單元中的計算效率。此外，當(dāng)模型從高維度降至較低維度，計算時若維持最大可能連接單元數(shù)不變，則將產(chǎn)生更多零系數(shù)；但若減少最大可能連接單元數(shù)，則需修改局部代碼。顯然，這種建模求解方式并不高效、通用。因此，需要尋找合理的系數(shù)矩陣優(yōu)化方法，在保證模擬精度的前提下，提高模型在大規(guī)模網(wǎng)格單元中的計算效率及通用性。

空間鏈接器(spatial linker)是指多維空間中有效計算單元間的虛擬鏈接，被用于記錄兩相鄰有效單元間的連接關(guān)系?？臻g鏈接器式建模思路即預(yù)先分類提取不同方向(維度)上有效計算單元間的空間鏈接器，如此可避免產(chǎn)生零系數(shù)，節(jié)省數(shù)組存儲空間，優(yōu)化計算效率。該建模思想最早被用于大型河網(wǎng)水流計算中[12]，一維河道中的兩控制節(jié)點利用鏈接器進行連接，而后又被推廣至二維水流模擬中[13]。飽和-非飽和流數(shù)值模型中也出現(xiàn)過類似的思路，例如TOUGH2軟件[14]中對無流量邊界設(shè)置無流量連接，但尚未采用該通用模式進行整體性建模。陳景波等[15]曾在三維飽和-非飽和流建模中應(yīng)用此建模思路，但模擬精度一般，且無法在復(fù)雜邊界條件下應(yīng)用，限制條件較多。因此，尚未有基于該建模思路對飽和-非飽和流模型開展通用化模式的探討與研究。

本文針對多維飽和-非飽和流方程求解系數(shù)矩陣中常存在零系數(shù)的問題，提出在模型構(gòu)建中預(yù)先記錄有效單元間的連接關(guān)系(如連接方向、前后單元序號等)，避免計算和存儲零元素，實現(xiàn)有效連接的壓縮存儲；并結(jié)合矩陣標(biāo)識法，建立模擬精度較高的空間鏈接器式多維通用飽和-非飽和流模型。該項工作也旨在構(gòu)建一套具有自主知識產(chǎn)權(quán)、可求解復(fù)雜多維飽和-非飽和流問題的程序。

1 模型構(gòu)建

1.1 空間鏈接器及矩陣標(biāo)識法

飽和-非飽和流方程的詳細(xì)差分離散步驟可參考Dogan一文中[11]。所有非邊界單元的方程組最終可寫成如下矩陣表達(dá)式：

AH=R

(1)

式中：A——系數(shù)矩陣，具有主對角對稱、高稀疏性特征；

H——當(dāng)前時間步中的未知總水頭矩陣，為單列矩陣；

R——常數(shù)項矩陣，為單行矩陣。

假設(shè)在n個有效計算單元中有m個鏈接，則方程中的非零系數(shù)個數(shù)為n+2m，其中n對應(yīng)主元系數(shù)Ai,i的個數(shù)，2m則對應(yīng)非主元系數(shù)的個數(shù)。鏈接器兩端節(jié)點所記錄的兩非主元系數(shù)值相等且在矩陣中呈對角對稱，即Ai,j=Aj,i。將系數(shù)矩陣A中的非零系數(shù)逐行存儲于單行矩陣D中。定義空間鏈接器的數(shù)據(jù)存儲結(jié)構(gòu)為{istart, iend, d, a, b}，其中istart為首節(jié)點序號，iend為末節(jié)點序號，d為方向代號(0、1、2分別代表X方向、Y方向和Z方向，各方向均以靠近零點處為首節(jié)點)，a為首節(jié)點為主元時末節(jié)點系數(shù)Ai,j在矩陣D中的序號，b為末節(jié)點為主元時首節(jié)點系數(shù)Aj,i在矩陣D中的序號。多維中的鏈接器如圖1中所示。

圖1 空間鏈接器示意圖Fig.1 Schematic diagram of spatial linkers

由空間鏈接器存儲結(jié)構(gòu)可知，建立空間鏈接器信息過程中，需對各計算單元進行編碼。為適用于不規(guī)則邊界區(qū)域和微地形環(huán)境，需進行兩套單元的編碼，即空間網(wǎng)格單元(包括有效單元和無效單元)及有效計算單元，且為便于地表邊界條件設(shè)置，又需將有效計算單元分為地表單元和內(nèi)部單元兩類。以三維為例，離散網(wǎng)格單元按照自上而下、自前至后、從左往右的順序進行編碼，選取圖2(b)中的DEM數(shù)據(jù)進行三維網(wǎng)格單元剖分，垂向離散尺寸為0.01 m，對單元及鏈接器進行編碼，得到結(jié)果如圖2(c)所示。4層剖分單元總數(shù)為24，但有效單元僅有12個，地表網(wǎng)格按順序的存儲編號為{1,3,4,7,11}，鏈接器共16個，其中x方向有6個([1]—[6])，y方向有3個([7]—[9])，z方向有7個([10] —[16])，依次記錄空間鏈接器的數(shù)據(jù)存儲信息。

圖2 空間鏈接器與編碼Fig.2 Spatial linkers and encoding

以空間鏈接器提前記錄兩相鄰有效單元間的連接關(guān)系，實際上避免了不透水邊界單元及無水力聯(lián)系單元間的判斷設(shè)置，僅需對非零系數(shù)進行計算和存儲。設(shè)計多維方向代號可便于在多維度下進行拓展，模型更為通用。值得注意的是，空間鏈接器式建模方法需要在前期提取鏈接器信息，這將耗費一定時間，但隨著單元數(shù)量的增多及邊界條件復(fù)雜程度的增加，這部分耗時所占整體計算時間的比例將越來越少。

為更有效存儲矩陣中的非零系數(shù)、提高后期求解效率，本文引入矩陣標(biāo)識法，將原大規(guī)模系數(shù)矩陣寫成簡化的計算公式，其思想是將D中存儲的主元、非主元系數(shù)分開記錄。首先按行將系數(shù)矩陣A中的非零元素所在的列序號存于單行矩陣b數(shù)組中，然后將每行第一個非零元素在D中的序號存于單行矩陣J中，最后將每行主元系數(shù)在D中的序號存于主元標(biāo)識k中。原始矩陣可改寫成：

(2)

每個主元的水頭值可使用非主元水頭值表示：

(3)

1.2 邊界條件及模型參數(shù)

模型中考慮了三類邊界條件，包括水頭邊界條件(Dirichlet)、流量邊界條件(Neumann)及滲流面邊界條件(Seepage face)，其中滲流面邊界條件是一種混合型邊界條件，每一時間步長中滲流面的實際位置需要通過迭代計算得到[16]。

求解飽和-非飽和流方程時，還需建立壓力水頭-土壤含水量(h)及壓力水頭-非飽和導(dǎo)水率K(h)之間的本構(gòu)關(guān)系。壓力水頭與土壤含水量之間的關(guān)系被稱為土壤水分特征曲線，選用經(jīng)典的VG模型[17]對該曲線進行擬合，公式如下：

(4)

式中：m——經(jīng)驗形態(tài)參數(shù)，m=1-1 /n(n>1)；

θs——飽和含水率/(cm3·cm-3)；

θr——殘余含水率/(cm3·cm-3)；

α——進氣值的倒數(shù)/cm-1；

n——多孔介質(zhì)的均一性。

非飽和導(dǎo)水率K(h)的確定采用van Genuchten基于Mualem模型提出的改進公式[17]：

(5)

式中：Ks——飽和導(dǎo)水率/(cm·s-1)。

1.3 矩陣求解

模型中使用Intel?MKL函數(shù)庫中提供的并行稀疏直接求解器(Parallel Sparse Direct Solver，PARDISO) 對生成的矩陣進行求解。該求解器是目前最快的線性稀疏矩陣求解方法之一，對于大規(guī)模計算問題，表現(xiàn)出極高的計算效率和并行性，也在FEFLOW[18]中被應(yīng)用。求解中的收斂判斷選用混合指標(biāo)，即飽和帶以壓力水頭的變化值進行判斷，非飽和帶則以土壤含水量的變化值進行判斷。該模型在VS2008平臺利用C++語言開發(fā)。計算機配置為Intel(R) Core(TM) i7-6500U(2.50 GHz)處理器、16 G內(nèi)存。

1.4 誤差分析

選取平均絕對誤差(MAE)和均方根誤差(RMSE)這兩個指標(biāo)統(tǒng)計不同模擬結(jié)果與實測數(shù)據(jù)之間的偏差，公式如下所示：

(6)

(7)

式中：Si——模擬值；

Oi——實測值；

N——樣本數(shù)。

2 經(jīng)典案例模擬驗證

為檢驗?zāi)Ｐ驮诓煌瑮l件下的適用性及模擬精度，本文選取多維經(jīng)典案例進行模擬。

2.1 一維非飽和入滲案例

選取文獻(xiàn)[19]中的計算案例驗證模型準(zhǔn)確性。實驗土柱的長、寬均為50 cm，高為200 cm，數(shù)值模擬中垂向離散成40個網(wǎng)格單元(z=5 cm)。地表輸入雨強為50 mm/d，歷時10 d，時間步長為0.1 d，如圖3(a)所示。土柱中填充均質(zhì)土壤，飽和導(dǎo)水率Ks=10 cm/d，孔隙度n=0.45(近似為飽和含水量)。土壤中的本構(gòu)關(guān)系為：

(8)

式中：h——壓力水頭/cm；

ha——進氣壓/cm；

Sw——飽和度；

Swr——凋萎飽和度；

krw——相對導(dǎo)水系數(shù)。

已知Swr=θr/n=0.333，ha=0，代入式(8)，則Sw=1+0.667h/100，krw=(Sw-0.333)/0.667。初始時刻底板處為地下水位(壓力水頭為0)，地表處壓力水頭為-90 cm (Sw=0.3997，krw=0.1)，其余地方均為-97 cm(Sw=0.353，krw=0.03)。計算后輸出不同時刻沿深度的壓力水頭分布，與原文中的結(jié)果及UNSAT2模型[20]計算結(jié)果進行對比(圖3b)。

圖3 一維非飽和土柱邊界條件示意圖(a)和壓力水頭剖面分布(b)Fig.3 Boundary condition of the (a) 1D unsaturated soil column and (b) pressure head profiles for the unsaturated flow verification case

由圖3可知，本文計算的不同時刻土柱垂向壓力水頭分布及變化過程與原文模型結(jié)果非常接近。以原文中的模擬結(jié)果為實測值，對比本文與UNSAT2的模擬結(jié)果誤差(表1)。由表中可知，本文計算的各時刻、各深度處壓力水頭誤差指標(biāo)均小于UNSAT2。因此，本文提出的多維通用模型比UNSAT2模型更適用于一維非飽和帶中土壤水分運動模擬。

表1 模擬誤差統(tǒng)計

2.2 二維及三維地下水位局部起漲案例

二維地下水位局部起漲案例選用Vauclin等[21]開展的實驗，該案例被廣泛用于驗證飽和-非飽和地下水流模型的計算精度[22]。實驗裝置為填滿砂土的600 cm×200 cm土槽，以槽底為零基面，初始水位65 cm。頂部中間100 cm寬度的土壤表面設(shè)置恒定補給流量R=3.55 m/d，其余區(qū)域封閉。因?qū)嶒瀰^(qū)域的幾何對稱性，僅需對一半?yún)^(qū)域(300 cm×200 cm)進行模擬，地表0～50 cm區(qū)間有穩(wěn)定補給流量。離散網(wǎng)格單元中Δx=10 cm，Δz=5 cm，初始總水頭H=h+z=65 cm，右邊界單元為定水頭邊界，總水頭恒為65 cm，底部為不透水邊界。模擬時長8 h，時間步長Δt=1 min。已知儲水率Ss=0，飽和導(dǎo)水率Ks=8.40 m/d，本文沿用Clement等[23]使用的VG模型參數(shù)，分別為θr=0.01 cm3/cm3，θs=0.30 cm3/cm3，α=3.3 m-1，n=4.1。因初始時刻表面土壤的含水量極低，前期計算時段內(nèi)土壤含水量變化極快，因此需迭代150～250次才能得到較穩(wěn)定、準(zhǔn)確的數(shù)值解，水量平衡誤差可控制在5%以內(nèi)。模擬結(jié)果如圖4～5所示。

圖4 不同時刻地下水位模擬值與原文中的結(jié)果對比Fig.4 Comparison of observed and simulated water table elevations at different times

圖5 不同深度的壓力水頭變化Fig.5 Changes in hydraulic heads with time at the different depths

由圖4中可知，本文模型在不同時刻對比原文中的地下水位起漲過程擬合較好，補給區(qū)范圍內(nèi)的地下水位抬升最多，8 h后漲幅達(dá)55 cm，離補給區(qū)越遠(yuǎn)，地下水位漲幅越小。將本文模型同另外三種模型(Hydrus-2D[24]、SU3D[11]、VSF[25]模型)的模擬結(jié)果同Vauclin文中實測值比較，結(jié)果如表2所示，可知本文模型整體模擬精度與Hydrus-2D、SU3D、VSF模型基本一致。本案例在Hydrus-2D中的計算耗時為1.7 s，而在本文模型中的計算耗時為3.5 s。圖5顯示不同埋深處的壓力水頭變化，其數(shù)值變化過程與原文中結(jié)果基本一致。X=5 cm剖面處于補給區(qū)正下方，隨著水量逐漸向下滲透，土壤中的壓力水頭沿埋深依次響應(yīng)；而在X=165 cm剖面，土壤水更多通過飽和側(cè)向流補給，因此越靠近地下水位，壓力水頭值越早發(fā)生變化，位置1處離地下水位最遠(yuǎn)，壓力水頭始終不變，說明此處沒有受到水分的側(cè)向補給。

表2 模擬誤差統(tǒng)計

將上述案例中的y軸向外延拓，即可升至三維空間進行模擬，令y軸方向上的區(qū)域長為300 cm，離散網(wǎng)格單元中Δy=10 cm，表面水量補給區(qū)域范圍為50 cm×50 cm，保持其余條件、參數(shù)不變，最終的區(qū)域概況如圖6(a)所示。選取y=0剖面，將不同時刻地下水面線輸出，并與Dogan[26]博士論文中相同案例的計算結(jié)果進行比較(圖6b)，并將8 h補給后的三維地下水面位置進行展示(圖6c)。由圖5(b)可知，穩(wěn)定補給8 h后，y=0剖面補給區(qū)正下方的地下水面較之初始時刻上漲約15 cm，僅為二維案例中對應(yīng)位置處漲幅的四分之一，可見擴展至三維后，補給區(qū)的水量向更多周邊非飽和區(qū)域補給。與Dogan利用SU3D模型的結(jié)果相比，本文在0～100 cm范圍內(nèi)的地下水位模擬值略低于SU3D模型，最大差值約1 cm，但隨著計算時間的增加，至8 h時兩者的計算值已基本一致。本文t=8 h時的地下水面光滑平穩(wěn)，較好地刻畫了區(qū)域水頭梯度及補給方向，較之陳景波等[15]模擬所得有明顯突變的三維地下水面要合理很多。這些結(jié)果表明，本文模型的模擬精度與SU3D模型相當(dāng)，優(yōu)于陳景波等[15]提出的三維飽和-非飽和模型。

圖6 三維模擬區(qū)域概況(a); y=0剖面處隨時間變化的地下水面(b); 8 h補給后的三維地下水面(c)Fig.6 Description of (a) the 3D study area, (b) water table elevations at y=0 of various time and (c) 3D water table surface at the end of the rainfall

2.3 二維滲流面排水案例

為驗證滲流面這一特殊邊界條件的適用性，本文選取Vauclin等[27]于1975年設(shè)計的排水實驗進行模擬，實驗裝置與2.2節(jié)二維地下水起漲案例中的一致。土槽中的初始地下水位為145 cm，右側(cè)邊界水位瞬間降至75 cm，隨后下緣總水頭穩(wěn)定在75 cm。土柱中的地下水經(jīng)右側(cè)滲流面排水，整體地下水位逐漸下降至75 cm處。左側(cè)及下部邊界均設(shè)為不透水邊界，右側(cè)滲流面范圍隨地下水位下降而變化。模擬區(qū)域為300 cm×200 cm，離散網(wǎng)格單元中Δx=10 cm、Δz=5 cm，時間步長設(shè)為0.01 h。

使用原排水試驗?zāi)M中的Haverkamp關(guān)系式[28]及VG模型兩套土壤本構(gòu)模型進行結(jié)果比較。Haverkamp關(guān)系式中的土壤水分特征曲線及相對飽和導(dǎo)水率系數(shù)的表達(dá)式如下所示：

(9)

經(jīng)對比多篇參考文獻(xiàn)及模型結(jié)果驗證，確定其中參數(shù)：A=4.0×104，B=2.90，θs=0.30，θr=0.0，C=4.0×105，D=4.5，飽和導(dǎo)水率Ks=9.6 m/d。VG模型中的參數(shù)與 2.2節(jié)一致。加入Hydrus-2D進行模擬對比，計算中的土壤本構(gòu)關(guān)系選擇VG模型。多個時刻的結(jié)果對比如圖7所示。

圖7 不同時刻模擬的地下水位和實測地下水位的對比Fig.7 Measured and simulated water table positions at different times

由圖7可知，使用Haverkamp關(guān)系式(曾用于Vauclin 1975年模擬試驗)的地下水位模擬結(jié)果略好于VG模型，VG模型在0.5 h后的計算水位均高于實測值，在5 h后與實測值及原文模型結(jié)果一致。由表3可知，不同時刻本文模型的模擬精度與Hydrus-2D相當(dāng)。本案例在Hydrus-2D中的計算耗時為2.1 s，在本文模型中的計算耗時為4.5 s。本文模型可較準(zhǔn)確定位滲流面的邊界位置，可有效用于土體滲流問題的模擬。同時也發(fā)現(xiàn)，所選擇的土壤本構(gòu)模型及參數(shù)對土壤中模擬變量的分布影響極大。

表3 模擬誤差統(tǒng)計

3 田間降雨入滲過程驗證

3.1 研究區(qū)概況及參數(shù)測定

選取常州金壇區(qū)朱林鎮(zhèn)紅旗圩村一片稻麥輪作田進行采樣、監(jiān)測，并開展田間入滲模擬。區(qū)域地下水位埋深較淺，年內(nèi)平均埋深僅55 cm。土壤為脫潛型水稻土，剖面呈現(xiàn)分層特征：0～13 cm為耕作層(A)，13～27 cm為犁底層(Ap)，27～54 cm為滲育層(P)，54～75 cm為潴育層(W)，75 cm以下為潛育層(G)。其中0～60 cm深度主要為粉質(zhì)黏壤土，下層60～140 cm 為粉質(zhì)黏土。區(qū)域內(nèi)設(shè)立土壤水分監(jiān)測剖面，在10 cm、20 cm、40 cm、60 cm和80 cm處埋設(shè)土壤水分探頭，各層利用環(huán)刀取土樣測定干容重，0～10 cm、10～20 cm埋深處土壤干容重分別為1.20 g/cm3和1.23 g/cm3，20～40 cm、40～60 cm、60～80 cm埋深處土壤干容重均值在1.46 g/cm3左右，顯然表層土壤較之下層更為疏松。區(qū)域內(nèi)還建有氣象站與地下水位觀測井，詳見圖8。

圖8 實驗區(qū)布局概化圖Fig.8 Layout of the experimental area

田間各層土壤的水力參數(shù)、滲透性能各不相同，需要確定不同埋深處土壤水分運動參數(shù)的量級范圍。本文在田間3個剖面(0～10 cm、10～20 cm、20～40 cm、40～60 cm、60～80 cm)5個埋深范圍內(nèi)采土樣測定參數(shù)。對比采樣層和土壤發(fā)生層埋深范圍發(fā)現(xiàn)(表4)，第1、2、3采樣層范圍與對應(yīng)的土壤發(fā)生層范圍重合較多，此三層測量參數(shù)可代表土壤發(fā)生層內(nèi)參數(shù)；第4、5采樣層與對應(yīng)的土壤發(fā)生層范圍雖然錯位較多，但仍有局部交集，采樣層內(nèi)實測參數(shù)值可近似代表對應(yīng)發(fā)生層中的參數(shù)。利用壓力膜儀法測定土壤水分特征曲線，定水頭法測定土壤飽和導(dǎo)水率。以土壤吸力為 0 m 所對應(yīng)的含水量為飽和含水量；吸力為150 m所對應(yīng)的含水量為凋萎含水量。利用RETC軟件估算VG模型中的參數(shù)，測得的各參數(shù)如表4所示。由表可知，VG模型中的參數(shù)變異性不一：各土層的θs均一性相對較好；θr和n在20 cm以上均一性較好，而在20 cm以下變異性明顯大；α值變化最為敏感。各層飽和導(dǎo)水率的變異性也較大，并隨埋深逐漸增大，至60 cm以下逐漸減小，這可能與土壤質(zhì)地的變化有關(guān)。

3.2 參數(shù)敏感性分析

由式(4)、(5)中可知，土壤水分特征曲線受VG模型參數(shù)影響[29]，非飽和導(dǎo)水率則受α、n及飽和導(dǎo)水率影響[30]。因這兩個公式極其關(guān)鍵且均呈現(xiàn)高度非線性，為便于后期案例中的參數(shù)率定，本文選用擾動分析法[31-34]對公式各參數(shù)進行敏感性分析，即在不改變其他參數(shù)的前提下，將單一參數(shù)增加或減少一定比例。選取一組原始參數(shù)值：θr=0.09，θs=0.465，α=0.01，n=1.25，Ks=0.09，單位與表4中一致。對土壤水分特征曲線公式各參數(shù)加入±5%的擾動(-5%即原參數(shù)值乘以0.95的系數(shù))，對非飽和導(dǎo)水率公式各參數(shù)加入±5%、±10%的擾動，對比結(jié)果如圖9所示。由圖可知，n的擾動對土壤水分特征曲線的影響最大，其次是θs和θr，α影響最小；n的擾動也對非飽和導(dǎo)水率影響最大，其次是α，Ks的影響最小。因此，需要優(yōu)先對n進行率定調(diào)試。

3.3 參數(shù)率定及模擬結(jié)果對比

根據(jù)土壤發(fā)生層特征將土體分5層進行模擬，分別是：0～13 cm，13～27 cm，27～54 cm，54～75 cm，75～150 cm。利用相同機理的Hydrus-1D軟件手動率定土壤水運動參數(shù)，忽略地下水出流，參數(shù)l統(tǒng)一設(shè)為0.5。選取表層土壤含水量變化較大的20160402次降雨模擬一維下滲過程中不同埋深處的土壤含水量。結(jié)合測量的參數(shù)范圍及相似質(zhì)地土壤中的參數(shù)量級，最終得到率定后的各層參數(shù)，見表5。

表4 實驗區(qū)測定的土壤水力參數(shù)

圖9 參數(shù)敏感性分析Fig.9 Parameter sensitivity analysis注：(a)、(b)分別為VG模型參數(shù)加入±5%擾動對土壤水分特征曲線的影響; (c)、(d)分別為α、n、Ks三個參數(shù)加入±5%、±10%擾動對非飽和導(dǎo)水率的影響。

表5 率定的土壤水力參數(shù)

率定的VG模型參數(shù)大多在測定參數(shù)范圍內(nèi)，僅第1、2層的θr、α、n超出參數(shù)范圍較多，這可能與田間三處采樣剖面表層土壤的水力特性差異有關(guān)，表層土壤長期受人工耕作影響，疏松度不一，土壤水分特征曲線變異性較大，因此實測得到的參數(shù)范圍也受到較大影響。率定后的Ks隨埋深逐漸減小，均在測量值變化范圍內(nèi)。此外，次降雨中的表層土壤最大含水量與雨強呈正相關(guān)關(guān)系，更深處土壤層為地下水位頻繁變動區(qū)，該范圍內(nèi)土壤最大含水量受地下水埋深影響。具體而言，當(dāng)?shù)叵滤怀^某位置時，該位置由非飽和態(tài)變?yōu)轱柡蛻B(tài)，最大含水量也有所升高，這可能與非飽和帶土壤中滯留的空氣有關(guān)。因此，在后期計算中需要根據(jù)次降雨實測的最大含水量調(diào)整模型中的飽和含水量，變動范圍在0.04 cm3/cm3以內(nèi)。

為驗證各層率定參數(shù)的準(zhǔn)確性，對各層單一參數(shù)加入相同比例擾動，對比土壤含水量模擬結(jié)果。因θs受實測最大含水量制約，θr僅對土壤水分特征曲線有所影響，因此選取n、α、Ks這三個參數(shù)進行擾動?？紤]三個參數(shù)的敏感性差異，對n設(shè)計±5%的擾動，α設(shè)計±5%、±10%的擾動，Ks設(shè)計±10%、±20%的擾動，結(jié)果如圖10所示。由圖可知,n值的極小變化都對入滲過程影響極大，α在擾動達(dá)到10%時對20 cm和40 cm處的土壤水分變化過程造成一定影響，而設(shè)計最大擾動變化的Ks對含水量的變化過程影響最小?？傮w而言，所率定的各層土壤水力參數(shù)已達(dá)到較高的模擬精度。

圖10 土壤水力參數(shù)擾動對各層土壤水分模擬結(jié)果對比Fig.10 Comparison of simulated soil moistures of 5 layers using soil hydraulic parameter sets with some additional disturbance注：(a)為α、n擾動的影響; (b) 為Ks擾動的影響。

3.4 降雨入滲過程模擬

利用率定后的參數(shù)在本文構(gòu)建的模型中計算，除率定的20160402次降雨外，又選取20160406、20160420次降雨進行模擬。各案例中根據(jù)地表植被茂密度，扣除降雨前期一定的植被截留量，忽略降雨入滲過程中的蒸發(fā)，得到結(jié)果如圖11所示，其中圖11(a)為表層40 cm以上(包括40 cm)的土壤含水量變化，圖11(b)則為40 cm以下的變化。Hydrus-1D軟件常用于模擬降雨入滲補給過程[35-36]，此案例中使用該軟件計算結(jié)果和計算效率進行比較。經(jīng)比較，本文模型模擬結(jié)果與Hydrus-1D軟件中的結(jié)果十分相似，僅在60 cm、80 cm處略有差異。Hydrus-1D中三場次降雨過程計算耗時在0.4～0.8 s之間，本文模型中耗時約0.8～1.8 s，計算效率略低于Hydrus-1D。

圖11 三場次降雨模擬的土壤含水量和實測值的對比Fig.11 Comparison of measured and simulated soil moistures of three rainfall events

由圖11可知，整體而言，土壤含水量模擬值與實測值擬合較好，各深度處土壤含水量的絕對誤差均小于0.03 cm3/cm3，可滿足野外模擬的精度需求。在雨強相對較大的20160402和20160406次降雨中，表層40 cm以上模擬值的響應(yīng)時間有明顯延遲，這很可能與表層土壤中較多的大孔隙路徑有關(guān)，在雨強較大時，部分水量通過大孔隙快速入滲。在初始地下水埋深相對較淺的20160420次降雨中，下層土壤水分的模擬也出現(xiàn)一定的延遲現(xiàn)象，這也暗示著大孔隙流可能會造成地下水位的快速響應(yīng)，致使毛管邊緣帶中的土壤水分響應(yīng)也較快，但變幅遠(yuǎn)小于表層土壤。此外，利用模型統(tǒng)計了3場次降雨末期(計算時間小于32 h)各埋深層中土壤水分增量，發(fā)現(xiàn)次降雨計算時段末期，超過80%的入滲水量仍滯蓄在表層40 cm的土壤中，其中0～10 cm土層中的土壤水分增量約為次降雨總?cè)霛B量的40%～50%，10～20 cm土層中約占總?cè)霛B量的20%～30%，20～40 cm土層則在20%左右，而轉(zhuǎn)化為潛水的水量約占總?cè)霛B量的4%～12%，20160402次降雨潛水補給比例最小，20160406次降雨補給比例最大。

4 結(jié)論

(1)空間鏈接器式多維通用飽和-非飽和流模型利用多維方向上的空間鏈接器提前記錄兩相鄰有效單元間的連接關(guān)系，避免了不透水邊界單元及無水力聯(lián)系單元間的判斷設(shè)置，實現(xiàn)了系數(shù)矩陣的非零化計算和存儲，并引入矩陣標(biāo)識法以簡化系數(shù)矩陣。

(2)選取地下水位起漲、滲流面排水等四個經(jīng)典飽和-非飽和流案例驗證模型模擬精度及效率，結(jié)果表明：該模型在多維度、不同邊界條件(包括定流量邊界、滲流面等)下的模擬精度與Hydrus、VSF等成熟軟件相當(dāng)，計算效率略低于Hydrus軟件。

(3)對田間三場長歷時小雨入滲過程進行一維概化模擬發(fā)現(xiàn)，VG模型中的參數(shù)n在率定過程中敏感性最強，各層土壤含水量模擬值絕對誤差均小于0.03 cm3/cm3，結(jié)果與Hydrus軟件基本一致。因模型中尚未考慮土壤中大孔隙流影響，各層土壤水分響應(yīng)時間滯后于實測過程。三場次降雨計算時段末期，超過80%的入滲水量仍滯蓄在表層40 cm的土壤中，僅有4%～12%的水量已轉(zhuǎn)化為潛水。

本文模型被證明是求解多維飽和-非飽和流問題的高效工具，有望成為傳統(tǒng)多維飽和-非飽和流模型的重要補充。該通用化建模思路也非常適合將地表水、地下水模塊進行耦合。然而，記錄空間鏈接器信息時需循環(huán)遍歷所有離散單元，該部分前處理過程不可避免地增加了一定存儲空間和計算耗時。此外，模型中僅支持固定時間步長求解計算，這將增加簡單邊界條件下的計算耗時，未來仍需加入變時間步長計算方案，優(yōu)化計算效率。