空間鏈接器式多維通用飽和-非飽和流模型研究

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

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

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

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

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

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

1 模型構(gòu)建

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

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

AH=R

(1)

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

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

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

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

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

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

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

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

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

(2)

每個(gè)主元的水頭值可使用非主元水頭值表示：

(3)

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

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

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

(4)

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

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

θr——?dú)堄嗪?(cm3·cm-3)；

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

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

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

(5)

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

1.3 矩陣求解

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

1.4 誤差分析

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

(6)

(7)

式中：Si——模擬值；

Oi——實(shí)測(cè)值；

N——樣本數(shù)。

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

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

2.1 一維非飽和入滲案例

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

(8)

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

ha——進(jìn)氣壓/cm；

Sw——飽和度；

Swr——凋萎飽和度；

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

已知Swr=θr/n=0.333，ha=0，代入式(8)，則Sw=1+0.667h/100，krw=(Sw-0.333)/0.667。初始時(shí)刻底板處為地下水位(壓力水頭為0)，地表處壓力水頭為-90 cm (Sw=0.3997，krw=0.1)，其余地方均為-97 cm(Sw=0.353，krw=0.03)。計(jì)算后輸出不同時(shí)刻沿深度的壓力水頭分布，與原文中的結(jié)果及UNSAT2模型[20]計(jì)算結(jié)果進(jìn)行對(duì)比(圖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可知，本文計(jì)算的不同時(shí)刻土柱垂向壓力水頭分布及變化過(guò)程與原文模型結(jié)果非常接近。以原文中的模擬結(jié)果為實(shí)測(cè)值，對(duì)比本文與UNSAT2的模擬結(jié)果誤差(表1)。由表中可知，本文計(jì)算的各時(shí)刻、各深度處壓力水頭誤差指標(biāo)均小于UNSAT2。因此，本文提出的多維通用模型比UNSAT2模型更適用于一維非飽和帶中土壤水分運(yùn)動(dòng)模擬。

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

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

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

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

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

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

圖6 三維模擬區(qū)域概況(a); y=0剖面處隨時(shí)間變化的地下水面(b); 8 h補(bǔ)給后的三維地下水面(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 二維滲流面排水案例

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

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

(9)

經(jīng)對(duì)比多篇參考文獻(xiàn)及模型結(jié)果驗(yàn)證，確定其中參數(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進(jìn)行模擬對(duì)比，計(jì)算中的土壤本構(gòu)關(guān)系選擇VG模型。多個(gè)時(shí)刻的結(jié)果對(duì)比如圖7所示。

圖7 不同時(shí)刻模擬的地下水位和實(shí)測(cè)地下水位的對(duì)比Fig.7 Measured and simulated water table positions at different times

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

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

3 田間降雨入滲過(guò)程驗(yàn)證

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

選取常州金壇區(qū)朱林鎮(zhèn)紅旗圩村一片稻麥輪作田進(jìn)行采樣、監(jiān)測(cè)，并開(kāi)展田間入滲模擬。區(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)測(cè)剖面，在10 cm、20 cm、40 cm、60 cm和80 cm處埋設(shè)土壤水分探頭，各層利用環(huán)刀取土樣測(cè)定干容重，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)還建有氣象站與地下水位觀(guān)測(cè)井，詳見(jiàn)圖8。

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

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

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

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

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

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

表4 實(shí)驗(yàn)區(qū)測(cè)定的土壤水力參數(shù)

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

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

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

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

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

3.4 降雨入滲過(guò)程模擬

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

圖11 三場(chǎng)次降雨模擬的土壤含水量和實(shí)測(cè)值的對(duì)比Fig.11 Comparison of measured and simulated soil moistures of three rainfall events

由圖11可知，整體而言，土壤含水量模擬值與實(shí)測(cè)值擬合較好，各深度處土壤含水量的絕對(duì)誤差均小于0.03 cm3/cm3，可滿(mǎn)足野外模擬的精度需求。在雨強(qiáng)相對(duì)較大的20160402和20160406次降雨中，表層40 cm以上模擬值的響應(yīng)時(shí)間有明顯延遲，這很可能與表層土壤中較多的大孔隙路徑有關(guān)，在雨強(qiáng)較大時(shí)，部分水量通過(guò)大孔隙快速入滲。在初始地下水埋深相對(duì)較淺的20160420次降雨中，下層土壤水分的模擬也出現(xiàn)一定的延遲現(xiàn)象，這也暗示著大孔隙流可能會(huì)造成地下水位的快速響應(yīng)，致使毛管邊緣帶中的土壤水分響應(yīng)也較快，但變幅遠(yuǎn)小于表層土壤。此外，利用模型統(tǒng)計(jì)了3場(chǎng)次降雨末期(計(jì)算時(shí)間小于32 h)各埋深層中土壤水分增量，發(fā)現(xiàn)次降雨計(jì)算時(shí)段末期，超過(guò)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次降雨潛水補(bǔ)給比例最小，20160406次降雨補(bǔ)給比例最大。

4 結(jié)論

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

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

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

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