沼液穴灌入滲影響因素分析及數值模擬研究

鄭 健，周鵬杰，張彥寧，殷李高，張恩繼

(1.蘭州理工大學西部能源與環境研究中心，甘肅 蘭州 730050；2.甘肅省生物質能與太陽能互補供能系統重點實驗室，甘肅 蘭州 730050；3.蘭州理工大學能源與動力工程學院，甘肅 蘭州 730050；4.陜西省土地工程建設集團，陜西 西安 710075)

沼液是禽畜糞便等有機物經過厭氧發酵后的殘留液體，其總固體含量小于1%[1]，含有水溶性及多種養分的速效肥料，保留了90%以上的發酵原料中的氮、磷和鉀，且氮素結構較原污水更優，施用可以顯著改善土壤的理化性質，提高土壤的保水保肥能力，緩解土壤的板結酸化[2]。與清水和常規施肥相比，施用沼液會降低土壤細菌值，增加微生物數量并提高土壤酶活性[3]，同時還能促進作物生長，提高農產品產量和品質[4-6]，降低病蟲害發生率[7]，沼液浸種還可提高種子發芽率和生長水平[8-10]。目前關于沼液施用的研究主要集中在作物產量、品質及土壤肥力等方面。

沼液穴灌技術是從設施農業生產中總結出的一種節肥技術，目前該技術的實施主要根據經驗，與其相關的試驗研究和理論分析尚少。土壤水分運動的分析是研究農田灌溉技術的基礎，準確掌握穴灌的土壤水分運動變化規律對于促進穴灌技術的發展具有極其重要的意義。國內外眾多學者對多種灌溉方式下的土壤水分分布問題做了研究，但從文獻內容來看，多集中在滴灌[11-12]、溝灌[13]、膜孔灌[14-15]等方面。而將沼液與穴灌技術結合，關于沼液灌溉入滲規律和灌溉技術參數等方面的研究鮮有報道。

基于此，本文通過設置不同沼液配比、穴孔直徑和土壤容重，研究這些因素對沼液穴灌入滲的影響，并運用COMSOL Multiphysics軟件模擬分析沼液穴灌入滲濕潤鋒運移及土壤含水率分布規律，以期為水—沼液一體化灌溉技術提供依據，使沼液穴灌技術建立在更為科學合理的基礎上。

1 材料與方法

1.1 試驗材料及裝置

試驗土壤選取甘肅省蘭州市周邊正常使用的溫室大棚內表層0～40 cm的土壤，土壤經過自然風干、粉碎后，過2 mm孔徑土壤篩備用，土壤顆粒級配見表1。按照中國土壤分類標準，其土壤類型為粉質壤土，田間持水量為22.23%(質量含水率)。沼液取自發酵物總固體質量分數(TS)為8%，以牛糞和麥秸為主要原料，正常使用的戶用沼氣池，靜置2個月，待其理化性質穩定后，用紗布過濾掉沼液中較大的顆粒備用。

試驗裝置主要由土箱和馬氏瓶組成，考慮到穴孔入滲斷面的對稱性，采用實際穴灌土壤濕潤體的1/12作為研究對象[16]，土箱為30°扇形有機玻璃裝置，高和寬分別為50 cm和40 cm，馬氏瓶用來提供恒定水頭，容量520 mL。試驗裝置如圖1所示。

1.2 試驗設計及方法

試驗設置土壤容重為1.35 g·cm-3和1.45 g·cm-3；依據前期預試驗結果設置3個沼液與水的配比(體積比)，即1∶4、1∶6、1∶8和1個純水處理(對照)，3個穴孔直徑(3 cm、5 cm和7 cm)；設定孔深為5 cm。每組試驗做3次，共計72組試驗，取平均值作為試驗結果，試驗方案見表2。

采用室內土箱試驗，土體采用分層填裝(每層5 cm)，以設定土壤容重計算每層裝土質量，采用精度0.01 g電子秤稱量；最后5 cm填土時，在土箱30°角處預留試驗設計尺寸的穴孔；將按比例配比好的沼液裝入馬氏瓶，用乳膠軟管分別將馬氏瓶出水孔和排氣孔與裝完土的試驗土箱連接，為保證恒定水頭，將排氣管固定在入滲孔與土體表面平齊的位置。為方便記錄試驗數據，在土箱外壁繪制間距5 cm的網格線。然后同時打開馬氏瓶的出水孔和排氣孔閥門，當沼液進入穴孔時開始計時，隨著試驗的進行，當濕潤鋒運移到與外壁水平向(H)和垂向(V)的網格線相切時，記錄一次入滲時間和馬氏瓶的讀數，并用記號筆在土箱外壁繪制濕潤鋒運移曲線。當入滲量達到520 ml時或入滲時間約12 h時同時關閉馬氏瓶的進水孔和排氣孔閥門，待穴孔內沼液滲完，試驗結束。為避免土壤水分再分布造成土壤含水率的變化，在試驗結束后立即沿著土槽中軸線水平方向和垂直方向每3 cm的交叉位置用土鉆取土一次，土樣編號并放入烘箱內烘干(105℃，8 h)后測定土壤質量含水率。

表1 土壤顆粒級配

1.馬氏瓶； 2.進水管； 3.排氣管； 4.穴孔；5.土槽； 6.濕潤鋒； 7.取樣點； 8.土體1. Markov bottle; 2. Inlet pipe; 3. Venting pipe; 4. The hole;5. Trough; 6. Wetting front; 7. Sampling point; 8. Soil.圖1 試驗裝置及取樣示意圖Fig.1 Schematic diagram of the test equipment and sampling

處理Treatment土壤容重/( g·cm-3)Soil bulk density沼液配比Biogas slurry:water ratio穴孔直徑/cmHole diameter處理Treatment土壤容重/( g·cm-3)Soil bulk density沼液配比Biogas slurry:water ratio穴孔直徑/cmHole diameterT11.351∶43T131.451∶43T21.351∶45T141.451∶45T31.351∶47T151.451∶47T41.351∶63T161.451∶63T51.351∶65T171.451∶65T61.351∶67T181.451∶67T71.351∶83T191.451∶83T81.351∶85T201.451∶85T91.351∶87T211.451∶87T101.35純水 Water only3T221.45純水 Water only3T111.35純水 Water only5T231.45純水 Water only5T121.35純水 Water only7T241.45純水 Water only7

1.3 基于COMSOL軟件模型建立

本文采用多物理場耦合軟件COMSOL Multiphysics進行水、沼液滲流模擬分析。COMSOL以有限元為基礎，通過求解偏微分方程(單場)或者偏微分方程組(多場)來實現真實物理現象的仿真[17]，具有大規模的計算能力。求解器是基于C++程序采用最新的計算技術編寫，它的求解器類型有直接、迭代、分離、多重格點；求解類型有穩態、瞬態、參數等；處理器支持多核與集群。根據需要，用戶可以通過MATLAB GUI界面選擇單場求解或者多物理場求解。內嵌豐富的CAD建模工具的有限元軟件COMSOL，可以直接進行二維以及三維建模。COMSOL在網格剖分方面具有強大的功能，它的網格生成器可以將腔體劃分為三角形或者四面體的網格單元，通過仿真系統的交互界面用戶可以根據自己的需要選擇剖分的精細程度。

圖2 直徑3 cm穴孔模型及網格劃分圖Fig.2 Cavity model and mesh map of 3 cm hole

穴孔灌的土壤水分運動屬于三維點源入滲條件下的土壤水分運動，具有軸對稱性。在此次模擬分析中，假定土壤均勻且各向同性，同時將穴孔灌的水分運動簡化為軸對稱條件下的二維土壤入滲問題進行處理，本次模擬將沼液與水配比后的液體視為混雜了顆粒物的渾水，模擬對象統一定為邊長40 cm的正方形軸對稱模型，在對稱軸頂點處留有不同尺寸的穴孔，見圖2。

土壤水分運動控制方程可表達為[18]:

(1)

式中，θ為土壤體積含水率；h為土壤負壓水頭，cm；Kh為土壤非飽和導水率，cm·min-1；t為時間，min；x為徑向坐標，cm，x向右為正；z為縱向坐標，cm。

1.4 模型參數

根據COMSOL模型模擬需要，對土壤水分特征曲線參數進行擬合，土壤水分特性曲線θh可用van Genuchten模型來表示，即

(2)

式中，θr為殘余土壤體積含水率；θs為飽和土壤體積含水率；α、n、m為土壤水分特征曲線的形狀參數，m=1-1/n，n>1。

根據表1土壤級配表，采用ROSETTA軟件擬合了兩種容重土壤的水力參數，結果見表3，其中Ks為土壤飽和導水率。

沼液中含有有機懸浮顆粒，且粘性較水高，在入滲過程中會影響沼液入滲[18]。因此，本文在模擬中沼液入滲的飽和導水率采用變量設置。為保證變量設置的合理性，本文采用課題組實測數據擬合后獲得的函數[19]。函數為冪函數的形式：

Kst=Ks1·t-α

(3)

式中，Kst為t時刻的入滲速率(cm·min-1)；t為入滲歷時(min)；Ks1為第一時間單位末的入滲速率，根據試驗測得；α為經驗入滲指數，通過Origin軟件對試驗數據進行擬合獲得。

1.5 定解條件的確定

1.5.1 初始條件 假定土壤初始含水量在所研究的區域內分布均勻，初始條件可表示為:

hx,z,t=h00≥z≥Z,0≤x≤X，t=0

(4)

式中，h0為土壤初始負壓水頭(cm)；Z、X分別指研究區域垂直方向、水平方向最大距離(cm)。

1.5.2 邊界條件 參照圖3可以看出：上邊界AB，試驗過程中不考慮蒸發，水流通量為零，則

(5)

表3 土壤水力參數

圖3 穴孔入滲模型Fig.3 Hole infiltration model

右邊界BC和下邊界DC，由于所選研究區域較大，可認為水分無法到達該邊界，邊界條件描述為：

右邊界BC，

hX,z,t=h0t>0

(6)

下邊界DC，

hx,Z,t=h0t>0

(7)

邊界OD為對稱軸,水流通量為零，

(8)

邊界ab、Ab為穴孔中入滲面，與入滲液體接觸，邊界ab，

hx,z,t=5cm (z=-5 cm)

(9)

邊界Ab，

hx0,z,t=hi

(10)

式中，x0為穴孔半徑，hi為作用于孔壁各點的水頭，cm，其數值等于孔壁各點到水面的距離。

2 結果及分析

2.1 孔徑對入滲的影響

圖4中黑色實線表示濕潤鋒的實際運移情況，不同顏色色塊代表濕潤鋒不同時段運移的模擬結果。從T7、T8、T9處理的濕潤鋒運移情況可以看出，在垂向入滲距離達到15 cm時，所需入滲時間分別為214.8、130.4、83.2 min，T8相對于T7、T9相對于T8的入滲耗時衰減率分別為39.3%和36.2%，即孔徑越大，水分入滲速率越快，且在相同時間內，孔徑越大，濕潤體的體積也就越大。這主要是由于水分在土壤入滲過程中初始階段主要受到土壤基質勢作用，在入滲的后期，土壤基質勢作用逐漸減弱，而重力勢作用逐漸加強，隨穴孔直徑的增大，穴孔內的沼液體積增加，導致重力場的變化，因而對入滲過程的作用愈加顯著，使入滲速率隨著穴孔直徑的增大而增大。為進一步分析孔徑對入滲的影響，本文以實際垂向入滲距離達到15 cm時的濕潤鋒為研究對象，分析實際水平入滲距離與垂直入滲距離的比值變化規律，結果表明T7、T8和T9處理的比值分別為1.007、0.952和0.945，T10、T11和T12處理的比值分別為1.050、0.950和0.902，T19、T20和T21處理的比值分別為1.061、0.938和0.948，T22、T23和T24處理的比值分別為1.127、0.959和0.886，由比值變化可以看出孔徑越大，比值越小，即濕潤體縱向長度隨孔徑增大而逐漸大于橫向長度(穴孔尺寸不計)。初步分析產生的原因是：隨著孔徑的增大，穴孔中入滲面面積也隨之變大，在孔深不變的情況下，即孔壁入滲面面積不變，孔底入滲面面積變大，相應的過水通道變多，使垂向入滲距離逐漸大于水平入滲距離。

圖4 不同處理濕潤鋒模擬值與實測值對比圖Fig.4 Comparison of simulated and measured values of wetting front under different treatments

2.2 土壤容重對土壤水分入滲的影響分析

為說明土壤容重對入滲的影響，選取T7與T19處理、T8與T20處理、T9與T21處理進行分析，從圖4中可以看出，T7與T19處理、T8與T20處理、T9與T21處理實際垂向入滲距離達到15 cm處所用時間分別為214.8 min與335.3 min、130.4 min與154.8 min、83.2 min與94.7 min，可以看出土壤容重1.45 g·cm-3的入滲速率明顯低于1.35 g·cm-3，即土壤容重大的入滲速率慢，這一試驗現象與很多學者的研究成果相符合[20-23]。張振華等[24]的研究指出，在入滲過程中，土壤容重對入滲率的影響具有穩定性和一致性，表現為對應于同一時刻，容重越大入滲率越低。分析形成的主要原因是隨著土壤容重的增大，土壤的團粒結構喪失，土壤變得緊密，從而土體內的大孔隙比例降低，導致土壤的基質吸附作用和毛管作用降低，使得入滲速率也隨之而降低。圖5中(a)、(b)、(c)和(d)分別表示1∶4(沼液∶水)、1∶6、1∶8和純水CK對照條件下相同孔徑在不同土壤容重條件下的累積入滲量對比，從圖5中可以看出，累積入滲量在相同時間內土壤容重呈1.45 g·cm-3<1.35 g·cm-3。例如圖5(c)中累積入滲量呈T7>T19處理，T8>T20處理，T9>T21處理。

2.3 沼液配比對入滲的影響

入滲液體濃度的差異會對土壤水分入滲產生影響[25]。圖4中T14、T17和T20處理表示孔徑5 cm，土壤容重1.45 g·cm-3，沼液水配比分別是1∶4、1∶6和1∶8的處理，從圖4濕潤鋒運移規律可以看出，當垂向入滲距離達到15 cm時，分別耗時244.1、210.0、154.8 min，即入滲速率隨著沼液配比度的升高而減慢。從圖6累積入滲量的變化規律可以看出，隨著水沼液混合液中水體積的升高，累積入滲量逐漸增大，并且混合液的入滲量衰減速度明顯快于純水，初步分析產生上述規律的主要原因為沼液溶液中微小的懸浮顆粒含量隨著沼液配比的增加而增多，在入滲過程中，懸浮顆粒堵塞了土壤中的孔隙[26]，降低了土壤的基質吸附作用和毛管作用，從而降低了沼液的入滲速率。

2.4 入滲模擬模型有效性驗證

2.4.1 土壤含水率模擬結果驗證 對24個處理進行模擬后共得到466個觀測點質量含水率數據，與實測值對比分析后，相對誤差在0～5%之間的點有290個，占62.2%，5%～10%之間的點有162個，占34.8%，10%～15%之間的有9個，占1.9%，大于15%的有5個，占1.1%。圖7和圖8分別顯示了相同土壤容重在距離OD軸6 cm和9 cm處取樣點含水率實測值與模擬值對比，從圖中可以看出模擬與實測值吻合程度較好，純水處理的模擬結果精確性要高于沼液；整個模擬結果中，誤差較大處大多出現在第一個和最后一個取樣點。實測值與模擬值的對比分析結果表明所建立的模型能較好反映穴孔自由入滲的土壤含水率分布情況。

圖5 不同處理累積入滲量曲線Fig.5 Cumulative infiltration curve of different treatments

2.4.2 濕潤鋒模擬結果驗證 圖4中，實線表示不同時刻實際濕潤鋒位置，不同顏色區域表示不同時刻模擬濕潤鋒位置，從圖中可以看出：(1)各時刻濕潤體形狀的模擬結果與實測結果一致性較好，濕潤鋒的模擬值與實測值可以較好的吻合；(2)由于沼液中顆粒物的存在，入滲速度衰減較快，模擬精度比純水要低；(3)在入滲后期實際濕潤鋒運移速度快于模擬濕潤鋒運移速度。

2.4.3 影響模擬的因素分析 模擬結果表明，土壤含水率和濕潤鋒的模擬與實測結果存在一定誤差，初步分析是受到土壤條件、沼液配比、試驗裝置和邊界條件設置的影響。主要體現在：(1)實際試驗中土壤采用分層裝填，可能存在局部的不均勻性；(2)模擬中對上邊界條件的簡化以及采用的土壤水分運動控制方程時忽略了溫度與土壤氣體對水分運動的影響，同時試驗土箱的側壁為光滑的有機玻璃，與土體接觸性較差，試驗過程中可能存在一定程度的優先流也造成了實測數值與模擬值之間的差異；(3)濕潤鋒邊界附近含水率變化較大，用土鉆取樣較容易產生誤差，并且濕潤鋒附近土體含水率值較小，測量容易產生偏差；(4)水和沼液的混合液在實際入滲過程中，土壤水分運動參數是動態變化的，本文采用實測數據的擬合方程代入，對模擬精度也造成了一定的影響。

圖6 不同沼液配比累積入滲量Fig.6 Relationship of accumulated infiltration volume withinfiltration time under different biogas slurry ratios

注：M表示實測值，S表示模擬值，下同。Note: M represents the measured value, S represents the simulated value, the same below.圖7 1.35 g·cm-3含水率模擬值與實測值對比圖Fig.7 Comparison of the simulated values of soil water content with the measured valuesunder soil bulk density of 1.35 g·cm-3

圖8 1.45 g·cm-3含水率模擬值與實測值對比圖Fig.8 Comparison of the simulated values of soil water content with the measured valuesunder soil bulk density of 1.45 g·cm-3

3 結 論

1)水-沼液一體穴孔自由入滲能力受穴孔孔徑、土壤容重和沼液配比的影響較為顯著。累積入滲量和土壤濕潤體隨穴孔直徑的增加而增大，隨沼液配比的增加而減小；水平入滲距離與垂向入滲距離的比值隨穴孔孔徑增大而減小，濕潤體形狀隨穴孔直徑的增加逐漸由水平軸大于垂直軸的橢球體趨向于水平軸小于垂直軸；1.45 g·cm-3土壤容重條件下的入滲速率和累積入滲量明顯低于1.35 g·cm-3土壤容重，表明隨土壤容重增大，相同時間內入滲速率、累積入滲量和濕潤體面積呈減小趨勢；隨沼液配比的增大，入滲速率和累積入滲量均隨之減小。在實際操作過程中，當土壤容重較大時可通過增加穴孔直徑和降低沼液水配比來增加土壤入滲量。

2)采用COMSOL Multiphysics建立穴孔入滲模型，并對室內試驗條件下的穴孔入滲試驗進行了數值模擬研究。各處理條件下濕潤鋒運移距離的實測值與模擬值吻合較好，均表現出初期運移速度較快，隨試驗時間的延長漸漸減緩，并趨于穩定。土壤含水率的實測值與模擬值對比來看，相對誤差<10%的占97%，說明所建立的模型合理，利用COMSOL Multiphysics可以準確描述水-沼液一體穴孔自由入滲過程中的水分運動。