王 策 張展羽 陳曉安 陳 于 朱成立 翟亞明

(1.河海大學(xué)農(nóng)業(yè)科學(xué)與工程學(xué)院， 南京 211100； 2.江西省水利科學(xué)院，南昌 330029；3.江蘇省農(nóng)村水利科技發(fā)展中心， 南京 210029)

0 引言

農(nóng)田土壤在動(dòng)物活動(dòng)、根系腐爛、風(fēng)化侵蝕和干濕凍融交替等條件下形成大孔隙結(jié)構(gòu)，導(dǎo)致灌溉水或施用養(yǎng)分未能被土壤充分吸收，在重力驅(qū)使下脫離基質(zhì)勢(shì)吸持作用向深層快速遷移，從而形成優(yōu)先流效應(yīng)[1-3]。優(yōu)先流現(xiàn)象可造成田間灌水效率降低、養(yǎng)分淋失[4-5]或地下水污染等風(fēng)險(xiǎn)。研究表明，95%的優(yōu)先流發(fā)生在孔徑大于250 μm的孔隙中，然而這些大孔隙僅占總孔隙的0.32%左右[6]。干縮裂隙隸屬膨縮性大孔隙結(jié)構(gòu)，發(fā)生于土壤脫濕收縮過程，其貫穿性和連通性結(jié)構(gòu)可為水分、懸浮顆粒等提供繞流路徑，使灌溉水大量損耗[7-8]，因此裂隙優(yōu)先流機(jī)理研究可為田間灌水效率提高和養(yǎng)分淋失控制提供理論依據(jù)，精準(zhǔn)模擬優(yōu)先流發(fā)育規(guī)律對(duì)農(nóng)田灌溉和水肥高效利用、地下水污染防治等具有重要意義[9]。

優(yōu)先流通常采用染色示蹤(碘化鉀、亮藍(lán)等)[10]、CT斷層掃描、穿透曲線等方法進(jìn)行觀測(cè)和機(jī)理研究。優(yōu)先流模擬過程中，可將土壤劃分為基質(zhì)域和裂隙域兩個(gè)孔隙域結(jié)構(gòu)，根據(jù)兩域形態(tài)或水流態(tài)同時(shí)模擬或分區(qū)模擬水流過程及其水分交換過程。目前，描述大孔隙優(yōu)先流模型主要有雙域模型或雙滲透模型等[11-15]。其中，雙滲透模型首次由GERKE等[16]提出，假定基質(zhì)域和大孔隙域(或裂隙域)為兩個(gè)相互獨(dú)立且統(tǒng)一的多孔介質(zhì)，分別采用Darcy-Richards方程以及兩套孔隙率、持水曲線和非飽和導(dǎo)水率進(jìn)行物理模擬，并采用兩域尺寸權(quán)重因子描述基質(zhì)域和裂隙區(qū)域比率，而兩域間水分交換則根據(jù)交界面平均導(dǎo)水率和水勢(shì)差推算。COPPLA等[17]根據(jù)裂隙愈合過程參數(shù)對(duì)雙滲透模型進(jìn)行了改進(jìn)，考慮了兩域水力參數(shù)及物理形態(tài)的動(dòng)態(tài)演變過程。雙域模型則是將土壤分割為兩個(gè)獨(dú)立區(qū)域，其中大孔隙流簡(jiǎn)化為重力驅(qū)使的重力流，其形式包括基于運(yùn)動(dòng)波的層流運(yùn)動(dòng)，哈根-泊肅葉的管狀流等[15]。例如，GRECO[18]構(gòu)建了考慮裂隙吸濕閉合的VIMAC模型，該模型將土壤研究區(qū)劃分為基質(zhì)域、膨脹裂隙域和永久大孔隙域，3個(gè)域分別采用Darcy-Richards流、考慮變形的裂隙邊壁層流方程和永久形態(tài)的大孔隙層流波動(dòng)方程模擬；LEPORE等[19]構(gòu)建的M&M模型將基質(zhì)域簡(jiǎn)化為方體多孔介質(zhì)，可采用Richards方程描述基質(zhì)中水分運(yùn)移，同時(shí)將裂隙域簡(jiǎn)化為基質(zhì)間的網(wǎng)狀/板狀通道，采用哈根-泊肅葉層流方程描述。兩域模型中，優(yōu)先流假定為由重力和粘滯性共同支配的粘滯層流運(yùn)動(dòng)[2]，采用薄壁型式的哈根-泊肅葉方程模擬，粘滯層流運(yùn)移參數(shù)可由時(shí)域反射儀測(cè)定，或由進(jìn)入大孔隙/裂隙內(nèi)的入滲通量推算。然而，針對(duì)不同類型土壤，其水力/物理特征、裂隙邊壁糙率、邊壁傾角、層流厚度等差異性顯著，導(dǎo)致雙域模型無(wú)法精確概化；此外，多數(shù)模型通常簡(jiǎn)化表層邊界或裂隙邊壁邊界條件，忽略了裂隙邊壁層流推進(jìn)過程[20]，存在一定的誤差。據(jù)此，本文基于表面入滲、層流運(yùn)移、裂隙邊壁水平吸滲與灌水強(qiáng)度相平衡的水量平衡原理，構(gòu)建優(yōu)先流雙域滲透模型，采用碘化鉀-淀粉染色示蹤試驗(yàn)進(jìn)行驗(yàn)證，并對(duì)該模型進(jìn)行應(yīng)用，以期進(jìn)一步深化對(duì)裂隙優(yōu)先流機(jī)理的認(rèn)識(shí)，為農(nóng)田水肥管理和深層滲漏抑制提供理論依據(jù)。

1 材料與方法

1.1 模型理論

1.1.1雙滲透模型

定義基質(zhì)域和裂隙域的下標(biāo)分別為m(matrix)和c(crack)。裂隙優(yōu)先流可采用經(jīng)典的雙滲透模型(Dual-permeability model)進(jìn)行描述，該模型假設(shè)基質(zhì)域和裂隙域?yàn)閮蓚€(gè)水分相互交換的多孔連續(xù)介質(zhì)域，兩域均可采用Darcy-Richards方程描述，則兩域中的水分運(yùn)移表示為

(1)

式中Cm、Cc——基質(zhì)域和裂隙域比水容量，cm-1

Km、Kc——基質(zhì)域和裂隙域的非飽和導(dǎo)水率，cm/min

t——入滲時(shí)間，min

z——土層深度位置坐標(biāo)，設(shè)定向下為正，cm

Γw——兩域間的一階水分交換項(xiàng)，min-1

hm、hc——基質(zhì)域和裂隙域的水勢(shì)，cm

ξm、ξc——基質(zhì)域和裂隙域體積比，cm3/cm3

然而，裂隙孔徑通常達(dá)到毫米甚至厘米級(jí)，因此在裂隙飽和情況下的水通量或?qū)士刹捎没诠?泊肅葉層流方程模擬，其飽和導(dǎo)水率Kc-s公式為[13]

(2)

式中wc——裂隙寬度，cm

g——重力加速度，m/s2

ν——運(yùn)動(dòng)粘度系數(shù)，m2/s

上述方法即為雙滲透模型理論，即基質(zhì)域水分運(yùn)移可采用Darcy-Richards方程模擬；針對(duì)大孔隙的裂隙域，需采用哈根-泊肅葉層流方程模擬。然而，不同裂隙流情境下的邊界條件具有差異，對(duì)于哈根-泊肅葉層流方程，其層流厚度隨灌溉(或模擬雨灌)強(qiáng)度而變化，因而難以獲取。本文擬針對(duì)差異性上邊界條件對(duì)模型進(jìn)行改進(jìn)。

1.1.2基于水量平衡原理的優(yōu)先流雙域滲透模型

(1)基質(zhì)域-裂隙域控制方程

擬采用二維Darcy-Richards方程對(duì)基質(zhì)域進(jìn)行描述，方程為

(3)

式中θ——土壤體積含水率，cm3/cm3

Km(θ)——關(guān)于含水率的非飽和導(dǎo)水率，cm/min

x——研究區(qū)域橫向位置坐標(biāo)，cm

裂隙域中的層流厚度由上邊界供水強(qiáng)度R0和表層入滲能力(入滲率)i(t)共同決定?；诠?泊肅葉在平行板間的層流推導(dǎo)，得出裂隙間柱狀層流單寬水流量與層流厚度的相關(guān)關(guān)系為[13]

(4)

式中w(R0,i(t))——關(guān)于上邊界恒定灌水強(qiáng)度R0和入滲能力i(t)的層流厚度，cm

ic——進(jìn)入裂隙的單寬水流量，cm/min

式(4)的前一項(xiàng)為平板狀層流通量，后一項(xiàng)為直接落入裂隙的降雨強(qiáng)度，而在實(shí)際運(yùn)算過程中，平板層流(優(yōu)先流量)遠(yuǎn)大于直接進(jìn)入裂隙內(nèi)的雨強(qiáng)，因此R0可忽略不計(jì)。

裂隙邊壁吸滲與層流下移同時(shí)發(fā)生，基于達(dá)西定律有

(5)

式中ih——水平吸滲率，cm/min

式(5)即為基質(zhì)域-裂隙域間水平吸滲控制方程，可根據(jù)Green-Ampt(GA)模型簡(jiǎn)化為

(6)

式中sf——濕潤(rùn)鋒處基質(zhì)吸力，cm

xf——濕潤(rùn)鋒與裂隙邊壁的距離，cm

根據(jù)GA模型及土壤導(dǎo)水率函數(shù)，可獲取裂隙邊壁水平吸滲率關(guān)于時(shí)間的函數(shù)。

式(3)～(5)分別為基質(zhì)域二維入滲控制方程、裂隙域邊壁層流方程、裂隙邊壁水平吸滲方程。進(jìn)入裂隙的水流通量由上邊界供水強(qiáng)度和表層入滲能力共同決定，因此需分析上邊界條件。

(2)上邊界條件

假設(shè)上邊界為恒定供水條件(模擬微噴)，R(t)=R0，且有R0>Ks。當(dāng)土壤表層干燥時(shí)，上表層入滲能力i(t)=-K(θ(z,t))?h/?z(h為土壤負(fù)壓水頭)隨時(shí)間逐漸由+∞降低至Ks。上邊界條件在不同階段的入滲率i(t)可采用Smith公式描述為

(7)

式中A、a——入滲參數(shù)，可采用薄層低水頭入滲條件測(cè)算

t0——漸近線橫坐標(biāo)，min

tp——徑流發(fā)生時(shí)刻，等效于積水時(shí)刻(Ponding point)或優(yōu)先流發(fā)生時(shí)刻，min

當(dāng)0≤t

(8)

式中iew——上邊界表面產(chǎn)生的多余水分流量，cm/min，其隨著時(shí)間t逐漸增大，后期R0與之比較可忽略不計(jì)

在已知土壤持水/導(dǎo)水率參數(shù)下，可采用積水時(shí)刻模型對(duì)積水發(fā)生時(shí)間tp(min)進(jìn)行模擬，公式為

(9)

(10)

式中Ip——積水時(shí)刻對(duì)應(yīng)的總?cè)霛B量，cm

θs、θi——飽和含水率和初始含水率，cm3/cm3

s——土壤基質(zhì)吸力，cm

Kr——相對(duì)導(dǎo)水率，(cm/min)/(cm/min)

土壤持水特征(VG模型)及非飽和導(dǎo)水率(Mualem模型)公式為

(11)

式中Se——有效飽和度

L、α、m、n——土壤持水特征/導(dǎo)水曲線參數(shù)，L常取0.5

其中m=1-1/n且0

(3)基于水量平衡的裂隙流模型

當(dāng)0≤t

(12)

式中 Δl1——裂隙邊壁上濕潤(rùn)長(zhǎng)度，cm

隨時(shí)間推移，邊壁Δl1上逐漸飽和，則層流進(jìn)入下個(gè)濕潤(rùn)段Δl2，則在tp+Δt≤t

(13)

同理，在tp+(n-1)Δt≤t

(14)

式(14)即為基于水量平衡原理的裂隙流水平吸滲方程，基于前期測(cè)得的土壤持水/導(dǎo)水率函數(shù)，根據(jù)經(jīng)驗(yàn)公式可得ih、iew和tp；其次設(shè)定步長(zhǎng)Δt，采用迭代法可獲取各時(shí)間段內(nèi)的濕潤(rùn)段Δl，以及裂隙邊壁水平吸滲濕潤(rùn)鋒距離xf，從而確定基質(zhì)及裂隙流濕潤(rùn)鋒分布。

上述模型即為基于水量平衡原理的雙域滲透物理模型及其相應(yīng)理念。然而在實(shí)際應(yīng)用過程中，其求解方法相對(duì)復(fù)雜，故擬采用數(shù)值模擬方法進(jìn)行簡(jiǎn)化，分段模擬優(yōu)先流過程。

1.2 試驗(yàn)方法及數(shù)值模擬

1.2.1研究區(qū)概況與供試土樣

采用碘化鉀-淀粉染色示蹤試驗(yàn)，對(duì)裂隙流模式進(jìn)行觀測(cè)，驗(yàn)證模型有效性。試驗(yàn)于河海大學(xué)節(jié)水園區(qū)節(jié)水與農(nóng)業(yè)生態(tài)試驗(yàn)場(chǎng)進(jìn)行，土壤試樣取自試驗(yàn)場(chǎng)內(nèi)稻麥輪作大田(31°86′N，118°60′E)，隸屬長(zhǎng)江中下游粘棕壤，總孔隙率53.5%，非毛管孔隙度約為6.8%。試驗(yàn)區(qū)年均降雨量和蒸發(fā)量為1 106 mm和900 mm，年均日照時(shí)長(zhǎng)為2 017.2 h。取表層0～40 cm旋耕土過10 mm篩后按原容重(1.25 g/cm3)逐層回填至尺寸為1.0 m×1.0 m×1.0 m的蒸滲儀內(nèi)，每層回填厚度10 cm且在層間打毛，保證層間水力均勻性和連續(xù)性。蒸滲儀前側(cè)可拆卸用于土壤縱剖，后側(cè)設(shè)置間隔10 cm、直徑1.0 cm的圓孔用以TDR監(jiān)測(cè)土壤水分水勢(shì)數(shù)據(jù)，見圖1a。裂隙培養(yǎng)前采用油灰刀平整表面。

1.2.2裂隙培養(yǎng)及染色示蹤試驗(yàn)

采用模擬噴灌預(yù)濕潤(rùn)處理試驗(yàn)，根據(jù)初始含水率和土層深40 cm內(nèi)田間持水量設(shè)計(jì)灌水定額，隨后隔絕降雨環(huán)境自然開裂，至表層含水率降至凋萎系數(shù)為止。裂隙發(fā)育過程中拍攝裂隙形態(tài)，裂隙深度采用帶有刻度尺的彈性塑料桿插入測(cè)定，土壤含水率采用TDR實(shí)時(shí)測(cè)定。裂隙培養(yǎng)完成后，采用碘化鉀-淀粉染色示蹤技術(shù)結(jié)合縱剖試驗(yàn)觀測(cè)水流運(yùn)移。在灌溉水中加入質(zhì)量濃度20 g/L的碘化鉀KI示蹤劑，對(duì)開裂土壤進(jìn)行灌溉處理，基于文獻(xiàn)[21-22]中優(yōu)先流試驗(yàn)灌水強(qiáng)度，設(shè)定60 mm/h(0.10 cm/min)作為灌水強(qiáng)度，歷時(shí)100 min，額定灌水定額為100 mm。由于灌溉水中的I-電荷與土顆粒外層電荷相同，因而不易被土壤吸附。灌水12 h后縱剖土壤試樣并對(duì)剖面噴灑硝酸鐵(Fe(NO3)3，質(zhì)量濃度20 g/L)和淀粉混合液(質(zhì)量濃度50 g/L)。三價(jià)鐵離子Fe3+的氧化性將灌溉水中的I-氧化為碘分子I2，I2與可溶性淀粉反應(yīng)生成靛藍(lán)色包合物，從而示蹤水流運(yùn)移模式。而后每4 cm逐層縱剖試樣并拍攝記錄優(yōu)先流染色模式。

1.2.3裂隙優(yōu)先流評(píng)價(jià)指標(biāo)

采用圖像處理技術(shù)對(duì)優(yōu)先流信息進(jìn)行提取(圖1b)，步驟包括圖像扭曲矯正、剪裁及修復(fù)、K-means聚類分割、二值分割等，可根據(jù)顏色深度將優(yōu)先流圖像分割為染色區(qū)域和非染色區(qū)，基于形態(tài)學(xué)算法提取相應(yīng)參數(shù)。采用的優(yōu)先流評(píng)價(jià)指標(biāo)如下：

(1)染色面積率Rdye(%)，定義為剖面染色面積Adye(cm2)占總面積At(cm2)的百分比，其中針對(duì)染色面積Adye，可統(tǒng)計(jì)二值圖像中染色像素?cái)?shù)量ndye，并將其轉(zhuǎn)換為尺寸量綱(cm)。

(2)剖面y內(nèi)任意點(diǎn)入滲深度，表示為

(15)

式中HI(x)——剖面y上表面坐標(biāo)(x，y)處所對(duì)應(yīng)的入滲深度，cm

M——剖面圖像像素行數(shù)

β(x,z)——剖面y上任意一點(diǎn)(x，z)處的染色值，其中染色時(shí)β(x,z)取1，未染色取0

λ——像素轉(zhuǎn)換為尺寸量綱的轉(zhuǎn)換系數(shù)，取為0.05 cm/像素

(16)

式中N——剖面圖像像素列數(shù)

(4)層深為z時(shí)的染色覆蓋率，表示為

(17)

式中，Rdye(z)為土壤剖面z-Δz至z+Δz層間的染色覆蓋率(%)，該指標(biāo)為大孔隙流或裂隙流常用指標(biāo)，量化了不同層間土壤染色深度，描述了垂直方向的染色均勻性[23]。

(5)基質(zhì)流深度Hm(cm)，定義為表層至z層方向連續(xù)染色覆蓋率不小于95%時(shí)對(duì)應(yīng)的染色深度，可根據(jù)Rdye(z)和深度z關(guān)系曲線確定。

(6)優(yōu)先流指數(shù)RPF，定義為優(yōu)先流染色面積APF占總?cè)旧娣eAdye的比率，估算式為

(18)

式中l(wèi)N——染色剖面寬度，cm

通過取樣可知，染色區(qū)試樣飽和度范圍為56.9%～66.3%，因此可忽略染色區(qū)域內(nèi)飽和度的差異，直接采用染色面積估算優(yōu)先流指數(shù)。

(7)灌水均勻度Eu，量化了灌溉區(qū)內(nèi)水量分布均勻度，是衡量灌水質(zhì)量的主要指標(biāo)，表達(dá)式為

(19)

其中

(20)

1.2.4雙滲透模型數(shù)值模擬分段求解

基于1.1.2節(jié)裂隙流雙域滲透物理模型，采用Hydrus 2D分段求解裂隙優(yōu)先流模式?？蓪⒘严毒W(wǎng)絡(luò)簡(jiǎn)化為方形網(wǎng)絡(luò)，將土壤分割為基質(zhì)域和裂隙域，取中心對(duì)稱區(qū)域作為研究區(qū)域進(jìn)行Hydrus 2D建模(圖2)。根據(jù)裂隙優(yōu)先流發(fā)生的時(shí)段，可根據(jù)時(shí)間節(jié)點(diǎn)0≤t≤tp和t>tp分別構(gòu)建數(shù)學(xué)模型。

當(dāng)0≤t≤tp時(shí)，未發(fā)生表面徑流和優(yōu)先流，OB為恒定強(qiáng)度邊界，OF、BC、DE和EF分別為零通量邊界(圖2)，則有定解條件

(21)

式中K(θ)——土壤變飽和導(dǎo)水率，cm/min

當(dāng)t>tp時(shí)，表面產(chǎn)生徑流并沿裂隙邊壁呈層流下移，形成優(yōu)先流；同時(shí)裂隙邊壁層流水平吸滲至基質(zhì)域內(nèi)，定解條件為

(22)

式中l(wèi)wf——裂隙邊壁上濕潤(rùn)段長(zhǎng)度，cm

需要注意的是，在t>tp階段，濕潤(rùn)段lwf隨時(shí)間逐漸向下推移，因此濕潤(rùn)邊界BG為動(dòng)態(tài)變化邊界，基于1.1.2節(jié)物理模型，在裂隙流發(fā)生后的任意時(shí)段(t,t+Δt)內(nèi)，存在水量平衡公式

(23)

式中i′h(t)——裂隙邊壁上單寬入滲率，cm/min

It——灌溉于表層邊界上的水分總量，cm

ΔVlar——裂隙內(nèi)增加的層流面積，cm2

通過Hydrus 2D軟件即可求得在時(shí)間步長(zhǎng)Δt(min)內(nèi)表層土壤入滲量Isur(cm)以及裂隙邊壁入滲量Ih(cm)，從而推算裂隙邊壁層流濕潤(rùn)鋒lwf的推移規(guī)律。

土壤持水特征曲線及非飽和導(dǎo)水率基于van Genuchten和Mualem模型，采用Hydrus神經(jīng)網(wǎng)絡(luò)模塊進(jìn)行預(yù)測(cè)，輸入?yún)?shù)包括機(jī)械組成(黏粒、粉粒和砂粒含量)和干容重。待試土壤機(jī)械組成和水力參數(shù)如表1所示。裂隙流Hydrus 2D模型的垂直剖面尺寸為70 cm，水平尺寸為40 cm，半裂隙寬度為1 cm(圖2)，水勢(shì)與試樣尺寸單位均為cm，初始時(shí)間與最小時(shí)間分別為0.144、0.014 4 min，最大迭代數(shù)設(shè)定10，目標(biāo)有限元網(wǎng)格尺寸1 cm，共計(jì)8 232個(gè)有限單元，入滲過程中暫不考慮裂隙愈合過程。

表1 裂隙優(yōu)先流試驗(yàn)各層土壤物理-水力參數(shù)Tab.1 Physical-hydraulic properties of soils in crack preferential flow experiments

1.2.5裂隙優(yōu)先流應(yīng)用

對(duì)裂隙流雙域滲透模型進(jìn)行應(yīng)用，分析了初始含水率(或裂隙形態(tài))和灌水強(qiáng)度對(duì)優(yōu)先流和灌水質(zhì)量的影響，選取初始體積含水率θi(裂隙寬度wc)和灌水強(qiáng)度R0進(jìn)行旋轉(zhuǎn)組合設(shè)計(jì)，設(shè)定6組灌水強(qiáng)度(0.10、0.08、0.06、0.04、0.03、0.02 cm/min)和4組

初始體積含水率(0.20、0.25、0.30、0.35 cm3/cm3)處理，各組處理灌水額定均為6.0 cm。由WANG等[24]研究可知，體積含水率在0.20～0.35 cm3/cm3之間時(shí)，裂隙面積率Rcr與含水率θi呈負(fù)相關(guān)關(guān)系，此時(shí)裂隙網(wǎng)絡(luò)骨架基本成型，含水率為0.20、0.25、0.30、0.35 cm3/cm3時(shí)所對(duì)應(yīng)裂隙面積率為9.5%、7.5%、5.5%和4.5%。數(shù)值模擬初始條件與邊界條件設(shè)定見表2。

表2 裂隙流數(shù)值模擬初始和邊界條件旋轉(zhuǎn)組合設(shè)計(jì)Tab.2 Rotation design of crack parameters (initial conditions) and boundary conditions for crack preferential flow

1.2.6模型評(píng)價(jià)指標(biāo)

采用決定系數(shù)(Coefficient of determination，R2)、均方根誤差(Root mean square error，RMSE)和平均絕對(duì)誤差(Mean absolute error，MAE)評(píng)價(jià)模型優(yōu)度[25]，其中R2越大，RMSE和MAE越小，模型優(yōu)度越高。采用皮爾遜相關(guān)系數(shù)r判別實(shí)測(cè)值和模擬值之間的相關(guān)性。其中，-1≤r≤1，當(dāng)r=1或r=-1時(shí)，實(shí)測(cè)值和模擬值之間呈極強(qiáng)性正相關(guān)或負(fù)相關(guān)，當(dāng)r趨向于0時(shí)，實(shí)測(cè)值和模擬值之間極弱相關(guān)或無(wú)相關(guān)。

2 結(jié)果與分析

2.1 裂隙優(yōu)先流試驗(yàn)分析

裂隙優(yōu)先流模式受裂隙率、彎曲度和連通性等因素影響，灌溉過程中超出表層入滲的水分通過裂隙優(yōu)先通道下滲，并同時(shí)水平吸滲至周圍基質(zhì)，形成以裂隙孔隙為中心、同時(shí)向周圍擴(kuò)散的滲流場(chǎng)[26]。染色示蹤劑中碘離子I-不易被土壤顆粒所吸附，因此在獲取染色劑濃度-吸附程度非線性特征曲線后，可采用碘化鉀淀粉示蹤法對(duì)優(yōu)先流進(jìn)行示蹤[27]。典型模式的裂隙優(yōu)先流如圖3所示。由圖可知，裂隙影響下的染色模式呈現(xiàn)出顯著的非均勻流和優(yōu)先流特征，由于微裂隙在吸濕過程中快速愈合，因此優(yōu)先流主要取決于貫穿裂隙[1]，其非均勻流特征表現(xiàn)為沿土壤深度方向染色覆蓋率逐漸減小，在裂隙影響下，染色模式可明顯區(qū)分為基質(zhì)流和裂隙流流態(tài)。剖面上層土壤染色模式均勻，10 cm內(nèi)的覆蓋率幾乎全部達(dá)到95%以上，呈現(xiàn)出基質(zhì)流流態(tài)，與裂隙形態(tài)無(wú)直接關(guān)系；土壤深度10 cm以下的染色覆蓋率隨深度急劇下降，呈現(xiàn)出優(yōu)先流流態(tài)，染色覆蓋率隨深度可分為3個(gè)階段：基質(zhì)流段(表層10 cm)、染色率驟降段(10～15 cm)、裂隙流段(15～60 cm)，如圖4所示。

不同深度染色覆蓋率直觀描述了優(yōu)先流在縱向剖面運(yùn)移軌跡的常用方法[28]，由圖4可知，各組試驗(yàn)的染色覆蓋率隨入滲深度逐漸降低。由圖3、4可知，表層邊界條件在入滲過程中的動(dòng)態(tài)變化促成了基質(zhì)流和染色流模式：①灌溉初期，開裂土壤基質(zhì)吸力較高，灌溉水全部轉(zhuǎn)化為入滲，此時(shí)裂隙對(duì)入滲未產(chǎn)生干擾。②隨表層含水率升高，其入滲能力逐漸降低，達(dá)到入滲臨界點(diǎn)，灌溉強(qiáng)度超過入滲能力且多余的水分進(jìn)入裂隙內(nèi)，以層流方式沿裂隙邊壁下移，并同時(shí)水平吸滲至基質(zhì)內(nèi)。優(yōu)先流模式與裂隙網(wǎng)絡(luò)形態(tài)密切相關(guān)，例如裂隙邊壁面積、裂隙彎曲度等均可改變?nèi)旧Ｊ剑疚乃芯康牧严稙閷?～2 cm、深20～40 cm的寬深裂隙。裂隙邊壁自上而下逐漸浸潤(rùn)至飽和狀態(tài)，同時(shí)伴隨著基質(zhì)域水平吸滲和吸濕膨脹過程。裂隙流模式呈現(xiàn)出表層基質(zhì)流流態(tài)和倒三角型的優(yōu)先流流態(tài)，如圖3所示。提取各剖面對(duì)應(yīng)表面裂隙率并與水分運(yùn)移參數(shù)進(jìn)行比對(duì)(表3)，對(duì)表面裂隙率和基質(zhì)流深度進(jìn)行相關(guān)性分析可知，皮爾遜相關(guān)系數(shù)r為0.299，決定系數(shù)R2為0.041，因此，兩變量呈現(xiàn)出弱相關(guān)性，基質(zhì)流深度與表層裂隙率不存在相關(guān)性，其主要取決于土壤物理/水力特征、表層初始含水率和灌水模式等。

表3 各剖層裂隙優(yōu)先流參數(shù)Tab.3 Parameters of crack preferential flow for nine profiles

2.2 裂隙流數(shù)值模擬與實(shí)測(cè)值

數(shù)值模擬基于圖3中剖面1裂隙參數(shù)進(jìn)行設(shè)定，其單裂隙與剖面呈垂直貫穿分布，且縱剖示蹤試驗(yàn)與灌水間隔較短，濕潤(rùn)鋒擴(kuò)散時(shí)間短，示蹤劑碘化鉀(KI)水動(dòng)力彌散過程相對(duì)微弱，故作為典型剖面進(jìn)行數(shù)值模擬。

基于剖面1表面裂隙寬度和深度參數(shù)作為邊界形態(tài)，采用Hydrus 2D對(duì)剖面1裂隙流進(jìn)行分階段數(shù)值模擬。土壤表層模擬灌溉強(qiáng)度為R0=0.10 cm/min，持續(xù)時(shí)間100 min，灌水定額100 mm。表層0～25 cm土壤的飽和導(dǎo)水率Ks為0.018 6 cm/min，小于0.10 cm/min，因此上邊界為變通量邊界，符合式(7)的Smith模型。模型分界點(diǎn)即積水時(shí)間點(diǎn)tp，該參數(shù)關(guān)系到裂隙流發(fā)生點(diǎn)、持續(xù)時(shí)間和進(jìn)入裂隙水通量等關(guān)鍵參數(shù)，具有顯著意義，因此擬采用試驗(yàn)記錄、物理推測(cè)和數(shù)值模擬3種方法對(duì)tp進(jìn)行測(cè)算。

試驗(yàn)監(jiān)測(cè)：因土壤表面糙率大、裂隙形態(tài)隨機(jī)性、收縮塊區(qū)不規(guī)則等原因，導(dǎo)致徑流時(shí)間點(diǎn)無(wú)法精確界定，基于試驗(yàn)觀測(cè)，估算徑流時(shí)間點(diǎn)tp為12.0～16.0 min。

數(shù)值模擬：采用Hydrus數(shù)值仿真表面入滲，監(jiān)測(cè)表面含水率變化，當(dāng)趨于飽和含水率時(shí)視為積水，數(shù)值模擬可得tp=14.4 min，并采用SPSS 25.0非線性擬合計(jì)算Smith公式參數(shù)A和a。據(jù)此可推算式(7)入滲參數(shù)t0、A和a分別為4.297 min、0.554和0.828，如圖5所示，其中Smith公式非線性回歸的決定系數(shù)R2為0.998，邊界條件模擬結(jié)果呈極顯著水平。

通過試驗(yàn)、物理推算及數(shù)值模擬計(jì)算可知，tp為15.0 min左右，在后續(xù)計(jì)算中按tp=14.4 min進(jìn)行數(shù)值模擬運(yùn)算。基于圖5可設(shè)定土壤上邊界變通量條件，以及進(jìn)入裂隙流量，據(jù)此根據(jù)水量平衡原理推算裂隙邊壁層流和水平吸滲參數(shù)。

獲取上邊界變通量條件后，基于優(yōu)先流雙域滲透模型理論，分階段對(duì)裂隙流過程進(jìn)行模擬。針對(duì)式(23)，在優(yōu)先流發(fā)生后選取Δt為10.0 min，裂隙流模擬過程及其相應(yīng)含水率分布如圖6所示。由圖可知，灌水初期濕潤(rùn)鋒呈基質(zhì)流流態(tài)(0≤t≤14.4 min)；當(dāng)表層入滲能力逐漸減小至R0時(shí)，到達(dá)積水時(shí)刻并初步發(fā)生優(yōu)先流，其中基質(zhì)流和優(yōu)先流可明顯區(qū)分，且裂隙對(duì)基質(zhì)流幾乎無(wú)影響(14.4 min85.0 min)，此時(shí)裂隙內(nèi)形成壓力水頭。裂隙流仿真模擬過程可分為3個(gè)階段：①供水控制階段。②表層入滲能力-裂隙邊壁層流控制階段。③表層/裂隙邊壁入滲能力控制階段。模型模擬結(jié)果與DICARLO[29]基于光傳輸?shù)膬?yōu)先流路徑圖像相一致，即優(yōu)先流/指流前端(Tip of the fingers)含水率趨于飽和(圖6h、6i)，而土壤上層優(yōu)先流路徑處含水率不高，可將該現(xiàn)象定義為優(yōu)先流指端超飽和度現(xiàn)象[29](Saturation overshoot at fingers of preferential flow)，該現(xiàn)象與基質(zhì)流含水率分布趨勢(shì)，即含水率沿深度由飽和逐漸減小的現(xiàn)象相反。

優(yōu)先流實(shí)測(cè)和模擬染色覆蓋率如圖7所示。結(jié)果表明，模擬流動(dòng)模式數(shù)據(jù)與實(shí)測(cè)值較為一致，實(shí)測(cè)染色覆蓋率與模擬值極顯著相關(guān)(P<0.01)，決定系數(shù)R2為0.981，說(shuō)明本文提出的裂隙流雙域滲透模型能較好地預(yù)測(cè)基質(zhì)流深度和裂隙流流態(tài)。由于裂隙深度測(cè)量技術(shù)限制，以及裂隙剖面形態(tài)假設(shè)偏差，導(dǎo)致裂隙內(nèi)入滲深度方向上存在一定差異。

2.3 優(yōu)先流模型應(yīng)用分析

裂隙骨架先于裂隙面積率提前到達(dá)穩(wěn)定[24,30]，因此假定各含水率下裂隙骨架形態(tài)相同，不同的是裂隙面積率或裂隙寬度?；谠撃Ｐ蛿?shù)值模擬仿真了4組初始含水率(裂隙形態(tài))和6組灌水強(qiáng)度下的優(yōu)先流模式，如圖8所示。橫向?qū)Ρ认嗤屎筒煌嗨畯?qiáng)度下的裂隙流可知，灌水強(qiáng)度增加，裂隙流愈加明顯，當(dāng)灌水強(qiáng)度R0=0.02 cm/min接近飽和導(dǎo)水率Ks=0.018 63 cm/min時(shí)，濕潤(rùn)體幾乎呈現(xiàn)基質(zhì)流狀態(tài)?？v向?qū)Ρ认嗤嗨畯?qiáng)度和不同初始含水率下的裂隙流模式可知，優(yōu)先流程度隨著初始含水率的增加而降低，且基質(zhì)流深度隨之增加，歸因于初始含水率升高使土壤導(dǎo)水能力呈數(shù)量級(jí)增加，而基質(zhì)吸力相比于導(dǎo)水率增幅偏小，因此基質(zhì)內(nèi)部水分傳導(dǎo)能力增強(qiáng)；其次，裂隙深度一定程度上決定了優(yōu)先流最大入滲深度，當(dāng)裂隙內(nèi)逐漸充滿水時(shí)，優(yōu)先流路徑達(dá)到深度下限，水流經(jīng)由裂隙底部緩慢下滲，因此初始含水率通過改變初始裂隙形態(tài)(寬度、深度和裂隙網(wǎng)絡(luò)骨架等)和基質(zhì)導(dǎo)水特性，而對(duì)優(yōu)先流產(chǎn)生顯著影響。此外，高初始含水率下的優(yōu)先流水平吸滲程度顯著增加，其原因?yàn)殡S著初始含水率增加，積水時(shí)刻或裂隙流發(fā)生點(diǎn)提前，裂隙內(nèi)充水并產(chǎn)生壓力勢(shì)，加速裂隙內(nèi)水分?jǐn)U散。

圖9分析了不同初始含水率(裂隙形態(tài))和灌水強(qiáng)度數(shù)值模擬下的灌水質(zhì)量指標(biāo)，包括裂隙優(yōu)先流指數(shù)、灌水均勻度、基質(zhì)流深度和最大入滲深度4個(gè)參數(shù)。結(jié)果表明，優(yōu)先流指數(shù)極值產(chǎn)生于灌水強(qiáng)度R0=0.10 cm/min和初始含水率θi=0.25 cm3/cm3時(shí)，當(dāng)灌水強(qiáng)度增加時(shí)，優(yōu)先流指數(shù)隨之增加，灌水均勻度降低。θi為0.20、0.25 cm3/cm3兩種情況下，初始含水率升高導(dǎo)致積水時(shí)刻提前，表面過剩水轉(zhuǎn)化為優(yōu)先流，然而此時(shí)裂隙寬度和深度較大，灌水結(jié)束后優(yōu)先流并未充滿裂隙，因此當(dāng)θi=0.25 cm3/cm3時(shí)優(yōu)先流指數(shù)偏高。由圖9可知，隨著初始體積含水率升高和灌水強(qiáng)度降低，基質(zhì)流深度顯著增大，優(yōu)先流大幅減弱?；诹严读鞣抡娼Y(jié)果可知，初始體積含水率、灌水強(qiáng)度和裂隙形態(tài)共同影響了裂隙優(yōu)先流程度，且很大程度上決定于灌水強(qiáng)度和裂隙深度。

3 討論

3.1 大孔隙/裂隙優(yōu)先流模型

優(yōu)先流在重力勢(shì)驅(qū)動(dòng)下使灌溉水沿優(yōu)先通道快速下滲，通常優(yōu)先流速相比于基質(zhì)流速高幾個(gè)數(shù)量級(jí)[17]。針對(duì)大孔隙或裂隙，可將土壤區(qū)域分割為基質(zhì)域和裂隙域，采用流態(tài)型式不同的Darcy-Richards基質(zhì)流模式，以及快速遷移的優(yōu)先流模式分別進(jìn)行物理模擬，同時(shí)考慮兩域間水分交換，即為雙域滲透模型[16]。然而，由于農(nóng)田大孔隙或裂隙復(fù)雜，難以精量化大孔隙率、基質(zhì)-裂隙界面特征、孔隙彎曲度和連通性等，且流態(tài)在基質(zhì)域和大孔隙間轉(zhuǎn)變復(fù)雜，因此加劇了優(yōu)先流模擬的難度[6]。在優(yōu)先流雙域滲透理論中，基質(zhì)域通常采用Darcy-Richards理論模擬是沒有爭(zhēng)議的，然而針對(duì)大孔隙域，BEVEN等[6]針對(duì)Richards方程在大孔隙流中的適用性提出質(zhì)疑：“多數(shù)物理模型將優(yōu)先流假定為雙連續(xù)介質(zhì)、雙孔隙或雙滲透媒介，這些概念或介質(zhì)構(gòu)想通常采用Richards方法進(jìn)行求解(不存在嚴(yán)格物理基礎(chǔ))，事實(shí)是Richards方程推導(dǎo)過程本身并不足以描述水分在異質(zhì)性土壤基質(zhì)中的遷移特征?！贝撕?，GERMANN提出了基于重力-粘滯力相平衡的優(yōu)先流Stokes層流運(yùn)移模式，且試驗(yàn)得出優(yōu)先流濕潤(rùn)鋒推移速率和灌水強(qiáng)度關(guān)系。該方法假定裂隙邊壁的水流動(dòng)量(或重力勢(shì))可因固體(土顆粒)和液體間的粘滯力而逐漸消散，然而邊壁層流厚度、水流通量、邊壁角、邊壁糙率等難以物理模擬，只能通過試驗(yàn)測(cè)定。其中，裂隙邊壁層流厚度主要取決于進(jìn)入裂隙的水流通量，通?？筛鶕?jù)平板型式的哈根-泊肅葉方程，將水流通量與層流厚度間建立相關(guān)關(guān)系[31]。

本文模型基于上述物理建模思路，對(duì)裂隙邊壁層流運(yùn)移過程中粘滯力-重力平衡過程進(jìn)行簡(jiǎn)化，可有效減少層流厚度等參數(shù)。該模型基本思路是水量平衡原理，即灌水強(qiáng)度可轉(zhuǎn)換為土壤表面入滲率、裂隙邊壁層流通量和裂隙邊壁水平吸滲率3部分，其中表面入滲率和水平吸滲率均可采用數(shù)值方法精確計(jì)算，亦可采用Smith入滲方程(表面入滲)和Green-Ampt物理模型估算。模型中基質(zhì)流入滲模擬較為精確，邊壁層流或優(yōu)先流的模擬主要取決于時(shí)間步長(zhǎng)的選擇，由于該模型是分段模擬和反復(fù)迭代過程，因此，時(shí)間步長(zhǎng)設(shè)定過小則會(huì)增加模擬工作量。借助于Hydrus 2D的裂隙流仿真模擬結(jié)果表明，本文提出的雙域滲透模型模擬結(jié)果與試驗(yàn)值吻合較好(R2>0.98)，基質(zhì)流數(shù)據(jù)模擬效果較好，裂隙流模擬在深度和水平吸滲范圍上存在微弱偏差，歸因于裂隙深度設(shè)定的誤差，以及裂隙內(nèi)積水狀態(tài)的差異，即由于試驗(yàn)中裂隙網(wǎng)絡(luò)存在高連通性，因此裂隙內(nèi)水流在重力驅(qū)使下在三維區(qū)域內(nèi)運(yùn)移，導(dǎo)致模型中85～100 min內(nèi)的模擬條件與試驗(yàn)存在一定差異。綜上，本文針對(duì)尺寸大且具有對(duì)稱特征的干縮裂隙結(jié)構(gòu)，基于水量平衡原理精準(zhǔn)模擬了優(yōu)先流過程，豐富和完善了雙域滲透模型理論。

3.2 初始體積含水率及灌水強(qiáng)度對(duì)優(yōu)先流的影響

土壤水分狀態(tài)對(duì)入滲深度及優(yōu)先流指數(shù)具有顯著影響。基于本文模型對(duì)裂隙流進(jìn)行應(yīng)用分析，模擬仿真了不同初始含水率及灌水強(qiáng)度下的優(yōu)先流模式。土壤初始體積含水率通過改變大孔隙/裂隙形態(tài)、非飽和導(dǎo)水特征和優(yōu)先流發(fā)生點(diǎn)，從而影響優(yōu)先流機(jī)理[32]。非穩(wěn)定流、指流或優(yōu)先流在干燥土壤中更加顯著，而濕潤(rùn)處理中優(yōu)先路徑并未完全發(fā)育為指流特征[32]。在原位或野外試驗(yàn)中，發(fā)現(xiàn)了未被染色的干縮裂隙，其表明了僅存在與指流水動(dòng)力連通的優(yōu)先路徑參與優(yōu)先流傳導(dǎo)過程，證實(shí)了大孔隙連通性在優(yōu)先流中的重要性。同樣地，GRANT等[33]通過田間原位染色示蹤試驗(yàn)得出，在初始含水率較低的粘性土中，優(yōu)先流經(jīng)由大孔隙結(jié)構(gòu)下滲深度更為顯著，相比于粉砂壤土具有較低的大孔隙域/基質(zhì)域水分交換作用。試驗(yàn)同時(shí)發(fā)現(xiàn)，干燥粘性土壤中優(yōu)先流入滲深度最大，初始水分條件通過影響干縮裂隙深度和連通性，從而導(dǎo)致次生的優(yōu)先流現(xiàn)象[33]。SALVE等[34]研究發(fā)現(xiàn)，在濕潤(rùn)季節(jié)的初期水分迅速下滲至深層基巖，此時(shí)非飽和帶處于最干燥狀態(tài)。NIMMO[3]將初始體積含水率機(jī)理歸結(jié)為以下原因：①土壤干燥時(shí)存在強(qiáng)疏水性，從而抑制大孔隙與基質(zhì)間水分交換作用，導(dǎo)致大孔隙流比例增加。②干燥土壤誘發(fā)裂隙，造成次生優(yōu)先流效應(yīng)，雖然濕潤(rùn)過程中裂隙逐漸愈合，然而愈合過程相比于優(yōu)先流發(fā)育速率較為緩慢。本文模擬結(jié)果與上述研究相一致，有效預(yù)測(cè)了初始水分狀況對(duì)優(yōu)先流影響機(jī)制；然而，本研究未考慮土壤親/疏水性和物理/水力特征在干燥程度間的變化，后續(xù)可進(jìn)一步試驗(yàn)驗(yàn)證。

有試驗(yàn)表明，優(yōu)先流隨著次降雨量或灌水強(qiáng)度的增加而增強(qiáng)，例如，UCHIDA等[35]根據(jù)多點(diǎn)原位試驗(yàn)比對(duì)發(fā)現(xiàn)，超過某閾值降雨定額后，優(yōu)先流效應(yīng)逐漸明顯，且土壤初始干燥程度越高則所需閾值越大；增加雨強(qiáng)則可有效激活優(yōu)先流路徑，在大孔隙土壤中灌水強(qiáng)度增加則優(yōu)先流效應(yīng)越顯著[3,36]，上述試驗(yàn)結(jié)論與本文仿真結(jié)果相一致。當(dāng)灌水強(qiáng)度低至飽和導(dǎo)水率時(shí)，灌溉水全部轉(zhuǎn)化為入滲；而隨著灌水強(qiáng)度逐漸增加，徑流時(shí)刻或優(yōu)先流發(fā)生點(diǎn)提前，表層入滲由Neumann 邊界(第2類)轉(zhuǎn)化為Dirichlet邊界條件(第1類)。高灌水強(qiáng)度下的優(yōu)先流效應(yīng)更為顯著，當(dāng)灌水定額相同時(shí)，灌水時(shí)間縮短，因此高灌水強(qiáng)度下基質(zhì)流穩(wěn)滲時(shí)間相比于低強(qiáng)度處理縮短，基質(zhì)流深度減小，而優(yōu)先流程度增加。本文從初始含水率(或裂隙形態(tài))及邊界條件(灌水強(qiáng)度)上有效模擬了優(yōu)先流效應(yīng)；然而，尚存在若干未知理論問題，例如土壤斥水性、濕潤(rùn)鋒不穩(wěn)定性、裂隙愈合過程、開裂程度及裂隙網(wǎng)絡(luò)形態(tài)等交互因素對(duì)優(yōu)先流的影響機(jī)理，仍需試驗(yàn)揭示或物理模擬。

4 結(jié)論

(1)采用染色示蹤試驗(yàn)揭示了9個(gè)剖面下的裂隙優(yōu)先流模式。結(jié)果表明，裂隙造成顯著優(yōu)先流效應(yīng)，裂隙流染色模式可明顯區(qū)分為基質(zhì)流態(tài)(土層深度0～10 cm)和優(yōu)先流態(tài)(土層深度大于10 cm)，染色覆蓋率在土層深度大于10 cm后急劇線性下降，基質(zhì)流深度與表面裂隙形態(tài)無(wú)顯著相關(guān)性，不受表層裂隙形態(tài)影響。

(2)基于表層入滲率、邊壁層流通量、裂隙邊壁水平吸滲3者與灌水強(qiáng)度相平衡原理，構(gòu)建了基于水量平衡理論的優(yōu)先流雙域滲透物理模型，并根據(jù)模型理念采用數(shù)值模擬方法進(jìn)行了驗(yàn)證。結(jié)果表明，該模型可有效模擬優(yōu)先流濕潤(rùn)鋒推移，模擬流態(tài)與實(shí)測(cè)流態(tài)相一致。該模型有效模擬了表層土壤變邊界條件特征、基質(zhì)流深度指標(biāo)和染色覆蓋率，實(shí)測(cè)和模擬染色覆蓋率極顯著相關(guān)(P<0.01，R2為0.981)。裂隙流過程可分為供水控制階段、表層入滲能力-裂隙邊壁層流控制階段和表層/裂隙邊壁入滲能力控制3個(gè)階段。

(3)采用初始含水率和灌水強(qiáng)度旋轉(zhuǎn)組合設(shè)計(jì)，基于本文模型進(jìn)行了應(yīng)用分析。結(jié)果表明，隨灌水強(qiáng)度增強(qiáng)，優(yōu)先流指數(shù)增加而灌水均勻度降低，積水時(shí)刻或徑流時(shí)刻提前；隨著初始含水率升高，基質(zhì)流深度增加，優(yōu)先流強(qiáng)度減弱。