黎明揚,劉廷璽,2,羅艷云,2,段利民,2,張俊怡,周亞軍
(1.內蒙古農業(yè)大學 水利與土木建筑工程學院,內蒙古 呼和浩特 010018;2.內蒙古自治區(qū)水資源保護與利用重點實驗室,內蒙古 呼和浩特 010018)
土壤入滲是生態(tài)水文循環(huán)過程中水分運移的重要組成部分,是降水通過土壤向地下水轉化的必然過程[1]。我國對于土壤入滲的研究開展的較早,主要集中在東部和觀測設施較為完備的流域[2-3],例如Sun 等[4]收集了我國較為濕潤的東部和南部42 個地區(qū)不同植被和土地利用類型下的土壤入滲數據,分析了大尺度下土地利用變化對土壤入滲速率的影響,并指出建立混農業(yè)林生態(tài)系統(tǒng)有利于提高農田和人工林的土壤入滲能力;Wang 等[5]在黃土高原退耕還林后的土壤上進行了生物土壤結皮對土壤入滲特性時間的變化研究。而位于我國北部以錫林河、巴拉噶爾河等為首的草原型內陸河流域開展的研究不多且不深刻,尤其是近些年在氣候變化與過度放牧二者的雙重影響下,原本生態(tài)脆弱的草原更是面臨嚴峻挑戰(zhàn)[6-7],因此探索和揭示干旱半干旱地區(qū)草原型流域特殊的生態(tài)水文過程和水分轉化運移機制及其影響因素,對進一步開展生態(tài)環(huán)境的保護和退化草場的治理工作具有重要意義。
傳統(tǒng)的大規(guī)模試驗費時費力且破壞性大[8-9],而土壤傳遞函數(Pedo-Transfer Functions,簡稱PTFs)則可利用有機質含量、容重、粒徑組成等土壤基本特性指標,間接計算出常用的土壤水分特征參數和曲線[10],在較大尺度的土壤入滲特性研究中是十分必要的,該方法可以在提高試驗效率的同時減少取樣對環(huán)境的破壞。為此,本文采用野外原位土壤入滲試驗與原狀土采樣獲取試驗數據,通過參數擬合選擇最能精確反映地區(qū)入滲過程的模型并利用其模擬穩(wěn)定入滲率,利用主成分和相關性分析研究影響水分下滲的因素,建立了土壤入滲特性與入滲環(huán)境要素之間的土壤轉換函數,并在此基礎上聯(lián)合同期土壤生態(tài)調查數據在小流域尺度上進行了穩(wěn)定入滲率的大面積預測。
2.1 研究區(qū)域概況本研究選擇內蒙古高原半干旱的錫林河流域為研究區(qū),其位于內蒙古自治區(qū)中部錫林郭勒盟,發(fā)源于赤峰市克什克騰旗寶爾圖山,流經錫林郭勒盟阿巴嘎旗,在貝力克牧場轉向西北流經錫林浩特市,最終注入查干諾爾沼澤地自然消失[11]。開展野外試驗的區(qū)域位于錫林河支流浩勒圖郭勒與干流錫林高勒河交匯處以上流域(43°24″—44°4″N,116°2″—117°15″E),面積為1852 km2,區(qū)域地勢三面環(huán)山,超過90%的植被為天然牧草,以茅針和羊草最為常見。研究區(qū)屬中溫帶半干旱大陸性季風氣候,多年平均降水量為266.8 mm,其中,6—8月降水量占年總降水量的50%以上[12]。
2.2 試驗設計首先利用中國1∶100 萬土壤類型圖劃分出試驗區(qū)包含厚栗黃土、草甸沼澤土、草原風沙土、石灰性草甸土、淡黑土共五種土壤類型分區(qū),按照每種土壤類型占試驗區(qū)面積比例,設置數個具有地區(qū)代表性的試驗點,包括土壤入滲試驗點38 個,土壤生態(tài)調查點62 個。土壤入滲試驗點和生態(tài)調查點分為7 組,其中6 組布設成垂直于河流與等高線的斷面,另設1 組覆蓋前6 組未包含的土壤、植被及土地使用情況;所有點位包含了該區(qū)域9 種土地利用類型、4 種植被類型、90%以上的海拔、坡度坡向及生物生長狀況(圖1與表1)。
圖1 研究區(qū)位置及采樣點分布
表1 研究區(qū)土壤特性及所含采樣點個數
為避免晨露和降水對土壤初始含水率(Initial Moisture Content,簡稱IMC)的影響,所有試驗于天氣連續(xù)晴朗的2018年7月26 至30日的10 至16 時進行,使用手持GPS 記錄實際采樣點經緯度坐標,在去除地表浮土及植被后進行,試驗期間未發(fā)生降水。入滲試驗使用雙環(huán)入滲儀,入滲前,在入滲環(huán)四周利用土鉆取土,每10 cm 一層共取5 層,用于測定IMC;入滲過程中,記錄內環(huán)及兩個馬氏瓶水位變化,前20 min 每1 min 讀數一次,20 min 后每3 min 讀數一次;試驗結束后,在距入滲環(huán)30 cm 處挖深于入滲深度的剖面(對于入滲深度不足50 cm 的點位則按照50 cm 開挖剖面),從表層向下取樣,每10 cm 一層共取5 層,每層采用100 cm3的環(huán)刀和自封袋各重復取樣3 個,分別用于測定土壤容重(Bulk Density,簡稱ρ)、粒徑組成、有機質(Soil Organic Matters,簡稱SOM)、地下生物量(Underground Biomass,簡稱UB)及其他土壤理化特性。各土壤生態(tài)調查點按照入滲試驗中的標準開挖50 cm 深的剖面取土,用于測定IMC 及各土壤理化屬性。土壤粒徑利用HELOS & RODOS 激光粒度分析儀進行干法測定,分級標準采用美國制:砂粒(2~0.05 mm,Sand Content,簡稱SA),粉粒(0.05~0.002 mm,Silt Content,簡稱SI),黏粒(<0.002 mm,Clay Content,簡稱CL);平均粒徑大小(Average Particle Size,簡稱APS)通過SA、SI、CL 三種顆粒所含百分比乘以每種顆粒的平均粒徑計算獲得;ρ采用環(huán)刀法測定;SOM 采用濃硫酸-重鉻酸鉀外加熱法測定;IMC 采用恒溫箱烘干法測定;UB 采用烘干稱重法測定。
2.3 入滲模型及評價方法土壤入滲模型可以模擬時間與入滲量之間的關系,研究選擇Green-Ampt、Philip、Kostiacov、Horton、Mezencev 及USDA-NRCS 六種入滲模型對雙環(huán)入滲過程進行模擬[13-18],見表2,使用Matlab 的Curve Fitting tool 進行參數擬合,最佳模型參數是通過計算比較不同參數的最小平方誤差決定的,將模擬參數帶入選擇的6 個入滲模型中,分別計算與實測數據相同記錄時間的累計入滲量。為了更好地識別多元非線性擬合中自變量個數對模型R2的影響[19],本文使用Adj-R2和RMSE 對入滲模型的模擬精度及誤差進行評價。
表2 土壤入滲模型
3.1 試驗數據的統(tǒng)計特征對于整個研究區(qū)而言,土壤水分的入滲速率很快,單位面積平均穩(wěn)定入滲率達0.2 cm/(min·m2);入滲初期速度的衰弱比例較高,前5 min 和20 min 分別達到76.60%與92.55%;整個入滲過程穩(wěn)定較快,20 ~60 min 內即可達到穩(wěn)定,這也側面說明了研究區(qū)這種砂質草甸地土壤垂向的成分和結構比較穩(wěn)定。由于入滲中后期的入滲速度衰減明顯,研究使用入滲總量的對數形式表示其隨時間變化的關系(圖2),由圖可知荒漠風沙土的入滲總量和速率最大,草甸沼澤土最小,僅有荒漠風沙土的三分之一,而其余三種類型土壤有較為接近的入滲過程。
圖2 入滲時間與log 形式入滲總量的關系
研究區(qū)表層土壤基本理化特征如圖3所示,根據粒徑組成,研究區(qū)主要屬砂質土壤,APS 均在75 μm 以上,較高的SA 使得整體ρ偏高,荒漠風沙土在五種土壤類型中各層位的APS 最大為93.77 μm。土壤IMC 與采樣點位到河流的距離關系密切,分布在河流兩側的草甸沼澤土和石灰性草甸砂土的IMC明顯高于其他三種離河較遠的土壤,尤其是分布在干流的草甸沼澤土,其各層土壤的平均含水率最高達到26.6%。在隨深度變化方面,各類型土壤的APS、IMC 和ρ變化不大,該結果與Zhang 等[20]和Yu等[21]對半干旱區(qū)植被與土壤理化性質關系的研究結果近似。降雨次數和降水量稀少使得土壤長期處于干旱狀態(tài),主要分布在平坦地區(qū)的草甸沼澤土、石灰性草甸土以及淡黑土三種土壤的SOM 隨深度增加呈下降趨勢,而主要分布在山地的厚栗黃土,其SOM 在20 cm 土層有明顯增高,這是因為植物為了生長在畜水困難的坡地上,深層表土根系更加發(fā)達,擁有大量毛細根的20 cm 土層經過常年的積累SOM 有明顯提高[22]。
圖3 5 種土壤粒徑、含水率、容重及有機質含量隨深度變化規(guī)律
表3 五種土壤類型下6 個入滲模型模擬參數值
表4 入滲模型模擬精度及誤差
圖4 五種土壤類型下6 個入滲模型模擬累計入滲量與實測值對比
3.2 入滲過程模擬各入滲模型的模擬值與實測值對比、模型參數及模擬結果的精度和誤差評價如圖4、表3與表4所示。
結合圖表可知,經驗模型的擬合效果要明顯優(yōu)于理論模型,這是由于理論模型對于實際狀況進行了概化和一定程度的簡化[23-24],在模擬情況相對簡單的入滲過程時誤差較小,例如荒漠風沙土在研究區(qū)的5 種土壤類型中植被稀疏,SOM、IMC 和UB 都處于較低水平,Green-Ampt 模型對其入滲過程的模擬誤差是該模型模擬5 種土壤中最小的,而其他4 種類型的土壤情況則相對復雜,模擬誤差也隨之增大(圖4(a)),此外Green-Ampt 模型在研究區(qū)整體模擬效果較差,研究判斷導致該結果的原因可能有以下兩個:(1)濕潤鋒面在砂粒含量較高、入滲速率較快的砂質草甸地并不是垂直推進的,這與Clemmens 的看法相似[25];(2)研究用于模擬推導出的入滲時間t 與入滲總量I(t)間的隱函數不適合使用MATLAB 的Curve Fitting tool 進行參數模擬。相比于理論模型,經驗模型更注重時間與入滲量之間的數學關系,尤其是在入滲過程的中后期模擬效果都比較優(yōu)秀[19]。
Horton 和USDA-NRCS 模型的綜合表現(xiàn)最佳,二者擬合結果Adj-R2均高于0.9 且RMSE 均小于0.1,說明以上兩種模型在不同土壤類型的模擬中均具有較好的普適性,該結論同國內外部分學者的研究結果一致[19,26-27];Kostiakov 和Mezencev 模型在入滲前期的模擬有一定出入,尤其是在含砂量較高的地區(qū)表現(xiàn)較差,這是由于以上兩個模型缺少用來描述初始入滲速率的常數項,而Mezencev 模型作為Kostiakov 模型的改良,其模擬精度與誤差卻均略低于Kostiakov 模型,說明砂質土壤入滲前期的高滲透性和較快的入滲速率下降比例會影響Mezencev 模型對于K ′參數的估計,這與Furman 等[28]在改良Kostiakov 模型滲透性能的研究結果近似。
Horton 模型的fc參數雖然是經驗模擬值,但是比較接近真實的穩(wěn)定入滲率[29],同時由于Horton 模型考慮了初始入滲速率,對于刻畫干旱半干旱型草原土壤前期入滲速度快、穩(wěn)定入滲時間短的特殊入滲過程最為貼切,因此研究選用Horton 模型推求各點位的穩(wěn)定入滲率以供下文入滲特征值的變異性分析及土壤轉換函數的建立使用。
3.3 影響因素分析土壤入滲過程受粒徑、ρ、SOM、IMC、UB 等因素的影響[30-31]。研究對38 個入滲點位的土壤理化性質及入滲環(huán)境條件進行了主成分分析,并統(tǒng)計了觀測尺度上入滲特征值與影響因素的變異系數和相關性,由于瞬時入滲速率較小,難以精確測量,尤其是50 min 之后的入滲速率,其也是通過一次測量多分鐘的累計入滲量后取平均獲得的,因此研究選用30、50、80 min三個時刻的累計入滲量和穩(wěn)定入滲率作為入滲特征值從而代表入滲過程,主成分分析結果如圖5所示。
主成分分析可以獲得兩個包含77.2%原始信息的主成分,研究將其分別命名為與土壤密實程度和孔隙分布相關的土壤自身物理屬性及與入滲液體動能勢相關的外界環(huán)境成分(圖5)。結果表明,主要分布在排列圖左側的入滲點位水分運移迅速,較大的APS 和較高的SA 增加了土壤的ρ,同時由于試驗點地理位置的原因,分布在序列圖第四象限點位的試驗時間也基本在中午,較高的水溫也會降低水的黏滯性[32],土壤與溫度二者共同提高了入滲速率。分布在河道附近的點位主要分布在序列圖的第一象限及距其較近的二四象限,這些點位的IMC 均超過26%,這不僅會降低土壤的入滲速率,還會促進植物的生長,導致SOM 和UB 達到研究區(qū)其他干旱區(qū)域的2 倍以上。有機質可以促進土壤中的小顆粒形成水穩(wěn)性土壤團聚體,其中一部分直接填充在土壤毛管與孔隙中,阻礙了水分的滲透,這與Liu 等[32]和Wang 等[5]對有機質影響入滲的研究結果一致;另外干旱半干旱地區(qū)的植物根系發(fā)達,具有固水固土的作用,可以減緩土壤水分的運移,這與Wu 等[33]在干旱草原上植被對土壤水文過程的影響研究結果相同。
變異系數描述了變量組間變異性的大小,CV≤0.1 時為弱變異性,0.1<CV<1 時為中等變異性,CV≥1 時為強變異性。由表5可以看出,研究區(qū)土壤理化性質的整體變異性中等偏弱,不同試驗點的累計入滲量具有中等變異性,且該差異性有隨著時間推進逐漸增大的趨勢;除了由于含量較低所引起變異劇烈的CL 外,土壤自身物理屬性的變異性不大,而與入滲液體動能勢相關的外界環(huán)境成分,如IMC、UB 以及SOM 三者的變異系數相對較大,說明研究區(qū)的累計入滲量和穩(wěn)定入滲率可能主要是由土壤非自身物理屬性的空間變異造成的。
圖5 土壤理化性質及試驗條件的主成分序列圖
表5 入滲特征值與影響要素的統(tǒng)計特征
入滲特征值與影響要素的相關性分析結果顯示,累計入滲量和穩(wěn)定入滲率與SA、APS 及Temp 成顯著正相關,與SI、CL、ρ、SOM、IMC 和UB 成負相關關系,該結果與劉繼龍等[34]在觀測尺度上土壤入滲特性與影響因素的相關性分析結果趨勢一致。由于研究區(qū)SI 和CL 均處于較低水平,SA 和APS與所選入滲特征的相關性接近,研究認為在研究區(qū)APS 對于土壤的粒徑組成有較好的代表性;SOM和IMC 與累計入滲量的相關性較高,與穩(wěn)定入滲率的相關性相對較低;而Temp 的相關性趨勢則相反,其與穩(wěn)定入滲率的相關性最高。
3.4 土壤轉換函數基于大量有關入滲特征值與影響要素建立PTFs 的研究可知,入滲特征值與影響要素之間有較好的對數關系[35-36]。文章首先使用觀測尺度上所有的影響要素與累計入滲量和穩(wěn)定入滲率建立對數關系的多元非線性PTFs(式(1)),結果顯示穩(wěn)定入滲率的轉換函數模擬效果優(yōu)秀,而累計入滲量的模擬誤差略大。為了進一步簡化PTFs,研究基于表5中相關性分析得出的結果,篩選出ρ、APS、SOM、IMC、UB 及Temp 等6 個與入滲特征值顯著相關的影響要素重新建立簡化的PTF(式(2))。兩種不同形式PTFs 的具體形式及評價結果如表6所示。
草原型流域土壤顆??傮w較粗,表層風化層較薄土質相對均一,有機質含量低且土顆粒不易形成團粒結構,經過簡化的PTF 精度雖略有下降,但穩(wěn)定入滲率的模擬結果依然良好,說明在入滲環(huán)境相對簡單的干旱半干旱型草原,通過土壤理化性質及環(huán)境變量建立來模擬土壤入滲特征值所建立的PTFs 形式簡單,可以簡化且適應性較好。然而由于試驗中缺乏對相同點位不同土壤初始含水率的采樣,而土壤含水率相比于土壤理化屬性,其在短時間內的變化更明顯,討論其對PTFs 參數的影響程度在實際應用中尤為重要;初始入滲速率和累積入滲量會直接因其改變了入滲水流濕潤區(qū)的平均梯度而發(fā)生改變,而穩(wěn)定入滲率則更穩(wěn)定,在參數的敏感性分析中更有意義,因此研究分別選擇初始含水率增加或減少5、10、15、20%以及無變化9 種情況下,兩種PTFs 計算fc的參數進行分析,結果如圖6所示,分析中當土壤含水率低于殘余含水率時取殘余含水率,高于飽和含水率時取飽和含水率,殘余含水率與飽和含水率分別取研究區(qū)平均值3.5 與38.9%。
圖6 不同初始含水率對簡化PTF 參數的影響
表6 不同影響因子土壤轉換函數參數與評價
結果顯示,當含水率有較大幅度降低時,PTFs 參數大都隨之改變,其中常數項a 與土壤顆粒組成及粒徑相關的參數b、c、d、f 有明顯變化,表明土壤初始含水率對PTFs 參數影響劇烈,當在較低含水率的情況下使用PTFs 應當適量增加試驗,以減小模型誤差。整體來看,土壤初始含水率的變化對未簡化的PTF 計算fc的影響較大,對簡化的PTF 影響較小,說明簡化的PTF 可以在一定范圍內更好地適應土壤初始含水率的變化。
研究使用在入滲試驗同期對62 個土壤生態(tài)調查點的采樣數據(圖1),聯(lián)合38 個土壤入滲試驗點數據,繪制了研究區(qū)小流域尺度上較為精細的PTFs 輸入參量的面尺度插值拓展,并分別使用相同1 km×1 km 大小的柵格繪制了使用Horton 模型和PTFs 預測的研究區(qū)土壤穩(wěn)定入滲率分布圖(圖7(a)和圖7(b)),并對二者的預測結果進行了比較(圖7(c)),其中圖7(b)假設試驗水溫恒定為20 ℃。
結果顯示,PTFs 預測的研究區(qū)土壤穩(wěn)定入滲率分布圖更貼合實際,該分布趨勢與黎明揚等[22]在使用合成孔徑雷達和PTFs 聯(lián)合反演錫林河上游表層土壤飽和導水率的研究結果一致。由于Horton 模型只能模擬試驗點位的穩(wěn)定入滲率(圖7(a)),因此相較PTFs 預測的結果(圖7(b))其不僅損失了大量細節(jié),在區(qū)域模擬的極值表現(xiàn)上也存在較大差異。在兩河交匯的濕地地區(qū),PTFs 預測的穩(wěn)定入滲率低至不足0.001cm/(min·m2),而Horton 模型在此處的預測值過高;兩河中部雖屬荒漠風沙土,但該地區(qū)僅生長著少量貧瘠的牧草,其實際入滲率要高于Horton 模型的預測值(圖7(c)),尤其是在研究區(qū)中部入滲率極高的裸露沙地,Horton 模型有較大程度的低估;整體來看,Horton 模型表現(xiàn)過于片面,少量的點位不足以把控整個研究區(qū)的入滲特征。
圖7 兩種模型預測穩(wěn)定入滲率分布及對比圖
使用PTFs 大面積預測土壤入滲特征值也存在一些不足,例如對下墊面環(huán)境等影響要素考慮不全面、簡化PTFs 帶來的精度缺失以及原位采樣帶來的空間限制等。前文提到,影響入滲過程的要素多且相互聯(lián)系,氣壓、空氣溫濕度等氣候條件,土壤孔隙度與相同類型土壤下不同的壓實程度都會對入滲過程產生影響[37];粒徑組成、植物根系的量與分布狀態(tài)都與ρ息息相關,由此角度分析,單純通過相關性等要素分析篩選PTFs 的輸入參數略顯片面,通過構建一個或多個綜合性的影響要素指標或許可以從機理上描述入滲條件與入滲特征值之間的轉換關系;結合遙感技術可以在一定程度上彌補原位采樣的空間限制并精細刻畫土壤物理要素的分布,但其對于深層土壤參數的獲取仍表現(xiàn)出較大的局限性[22],尤其是在土壤屬性縱向變化較大的地區(qū)具有一定的不確定性。
(1)半干旱型草原土壤砂粒含量高,水分入滲迅速,單位面積平均穩(wěn)定入滲率達0.2 cm/(min·m2)。使用6 種入滲模型對5 種土壤類型下38 個點位的入滲過程進行模擬,結果表明,Horton 和USDA-NRCS 模型在研究區(qū)的綜合表現(xiàn)最佳。
(2)天然草原入滲過程受外界干擾相對較小且影響要素的空間變異性呈中低等水平;主成分分析可以將9 個影響要素壓縮成2 個包含77.2%原始信息的主成分,研究分別命名將其為與土壤密實程度和孔隙分布相關的土壤自身物理屬性及與入滲液體動能勢相關的外界環(huán)境成分;相關性分析結果顯示ρ、APS、SOM、IMC、UB 及Temp 等6 個影響要素與入滲特征值顯著相關。
(3)研究分別通過使用全部影響要素和6 個顯著相關的影響要素與入滲特征值建立PTFs,結果顯示入滲特征值與影響要素間有較好的對數關系,使用顯著相關影響要素建立的PTFs 精度雖略有下降,但穩(wěn)定入滲率的模擬結果依然良好;使用Horton 模型和簡化的PTF 分別對研究區(qū)土壤的穩(wěn)定入滲率進行模擬,結果表明,Horton 模型的預測結果不僅損失了大量細節(jié),在極值的表現(xiàn)方面也存在較大的局限性;使用PTFs 預測小流域尺度土壤入滲特征值有助于提高試驗效率,減少環(huán)境破壞,該方法在復雜入滲環(huán)境下的魯棒性尚不明晰,有待于進一步研究。