李彥龍, 王鐵行, 王娟娟

(西安建筑科技大學 土木工程學院, 陜西 西安 710055)

黃土高原旱區(qū)淺層黃土水分場的數(shù)值分析

李彥龍, 王鐵行, 王娟娟

(西安建筑科技大學 土木工程學院, 陜西 西安 710055)

摘要：[目的] 研究黃土高原淺層土體水分場在春季連續(xù)干旱條件下的動態(tài)變化。[方法] 通過非飽和黃土二維非穩(wěn)態(tài)流有限元控制方程對黃土高原淺層土體水分場在不同氣象條件下的動態(tài)變化進行數(shù)值計算。[結果] 在連續(xù)干旱條件下水分場的計算結果與實測結果較為一致，受蒸發(fā)影響的土層最大厚度為1.8 m，連續(xù)3個月的干旱使其平均含水率降低至7.9%，土壤蒸發(fā)強度隨著含水率的減小而減小。強度較小且分散的降雨對于緩解連續(xù)干旱期黃土高原土壤干旱基本無效，大強度集中降雨只能在短時間內(nèi)緩解土壤旱情。[結論] 數(shù)值計算結果可為實際工程中植被類型的選擇以及灌溉時間的把握所參考。

關鍵詞：連續(xù)干旱期; 黃土高原; 水分場; 數(shù)值分析

黃土高原路基、市政等工程領域的植物防護及綠化工程規(guī)模隨著西部大開發(fā)戰(zhàn)略的實施而日益增大[1]。以干旱、半干旱地區(qū)為主的黃土高原，作為植物生長和防護工程載體的淺層黃土，其土壤水資源是植被承載力的關鍵因素[2-4]。黃土高原水分循環(huán)以上行蒸發(fā)和降雨入滲為主，大氣環(huán)流特征使其降水量在年內(nèi)分布極不均勻，冬干春旱，夏秋季降雨集中[5]。對植被建設以及植物防護工程而言，春季連續(xù)干旱氣象條件下最容易導致植被萎蔫死亡。對黃土高原淺層土體水分場在春季連續(xù)干旱氣候條件下的動態(tài)變化進行研究，可以為實際工程中植被類型選擇和灌溉時間的把握上提供相關資料。

淺層土體水分場在氣象條件下的動態(tài)變化是SPAC系統(tǒng)土壤水動力學的一部分，同時也是巖土工程、環(huán)境科學、水利工程、地學和土壤學等學科的交叉研究課題，以土壤學為基礎的單一過程的裸地蒸發(fā)[6-8]和降雨入滲[9-11]已有較多的研究，包括經(jīng)驗法和機理法。經(jīng)驗法在使用過程中具有區(qū)域局限性，其雖然可以計算累計入滲或蒸發(fā)量，但是不能計算土體內(nèi)任一點的水頭；以土壤學為代表的機理法主要為IRE方法[12]和以SWMS-2D為代表的有限元方法。

但是通過有限元方法并針對黃土高原淺層土體水分場在氣象條件下的動態(tài)變化卻少有研究。因此，本文擬通過巖土工程中的二維非穩(wěn)態(tài)流有限元控制方程并結合實測資料對黃土高原淺層土體水分場在連續(xù)春旱條件下的水分場進行分析，并通過數(shù)值計算來探明不同條件下的降雨對緩解春旱的效用。

1研究區(qū)概況

實測場地為陜北米脂縣楊家溝鎮(zhèn)何家岔，地理坐標為東經(jīng)109°49′—110°29′，北緯37°39′—38°05′，土壤以黃綿土為主；年平均溫度為8.9 ℃，日照時數(shù)為2 716 h，年平均降雨量為420.2 m，無霜期160～170 d；海拔951～1 029 m，由于長期粗放式經(jīng)營導致該區(qū)自然植被稀少。實測時間為2013年3—5月(共12周)，期間該地區(qū)累計降雨為10 d，累計降雨量為45 mm。其中，有效降雨為5 d，有效日均降雨量9 mm/d。因此，實測期間的氣象狀況基本可以滿足連續(xù)干旱的條件。

土壤水的受力狀況決定了它并不能完全被植物吸收和利用，只有介于萎蔫濕度和田間持水量之間的土壤水分才能被植物根系吸收和利用。介于該區(qū)間的土壤水又可分為難效水(萎蔫濕度到生長阻滯含水率，小于田間持水量的60%)，中效水(田間持水量的60%～80%)和易效水(田間持水量的80%～100%)，低于凋萎濕度的土壤水則為無效水。陜北米脂地區(qū)土壤難效水(6.8%～9.1%)，中效水(9.1%～12.4%)和易效水(12.4%～15.2%)[13]。用土鉆在標準地內(nèi)取樣測定裸地的天然含水率，分別在3月1號。3月28號,4月28號，5月28號取一次樣，取樣深度為6 m，0—1 m土層以0.1 m為一個層次，1—6 m土層每0.2 m為一個層次，共計35個層次，通過烘干法測含水率。

2計算模型及參數(shù)選取

2.1　計算模型

連續(xù)干旱條件下，淺層黃土水分場中任一點的水頭和滲透系數(shù)均隨著隨時間而變化，因此淺層黃土中的水分場為非穩(wěn)態(tài)流，其控制方程為[14]：

(1)

式中：[D]——剛度矩陣； [E]——容量矩陣； [F]——流量列陣； {h}——水頭列陣。

采用向后差分格式，式(1) 即為：

(2)

基質(zhì)吸力模型采用van Genuchten模型，其表達式如下：

(3)

式中：θ——體積含水率；θs——飽和含水率；θr——殘余含水率；h——吸力水頭(m)；α，m，n——系數(shù)，m=1-1/n。

對黃土而言，式(3) 中相關參數(shù)可按參考文獻[15]確定：

θr=-0.38+0.36ρd

(4)

θs=1-0.38ρd

(5)

α=exp〔(8.98-0.08T)+(-9.36+0.07T)ρd〕

(6)

(7)

(8)

將公式(3) 帶到式(8) 中可得：

(9)

參考文獻[16]給出了非飽和黃土水分擴散率與含水率之間的關系：

D=exp(aθ2+bθ+c)

(10)

式中：D——非飽和黃土的水分擴散率；a，b，c——系數(shù)，其中，a=84.5ρd-57.1，b=-56ρd+61.8，c=5.7ρd-12。

非飽和黃土比水容C式(11)表示：

(11)

式中：ρw——水的密度；g——重力加速度。

非飽和土的滲透系數(shù)與擴散率和比水容有如下關系：

kw=D×C

(12)

將公式(10)，(11)帶入式(12)可得：

(13)

2.2　水頭的確定

土體中的水頭通?？梢员磉_為：

h=hg+hm+hp+hs+hT

(14)

式中：h——總水頭；hg——重力水頭；hm——基質(zhì)吸力水頭；hp——壓力水頭；hs——滲透水頭；hT——溫度水頭。

溫度水頭和溶質(zhì)水頭通?？梢院雎裕瑢Ψ秋柡屯林械乃诌w移而言，認為土體所有孔隙與大氣是聯(lián)通的，因此壓力水頭為零。因此只需考慮重力水頭和基質(zhì)吸力水頭，式(14)可以簡化為：

h=hg+hm

(15)

其中，重力水頭按位置確定，吸力水頭由VG模型確定。

2.3　有限元模型

計算深度為6.0 m，所采用的三節(jié)點三角形單元均為直角三角形，共計劃分60層，0—1.0 m深度內(nèi)，每層0.05 m；1.0—3.0 m深度內(nèi)，每層0.1 m，3—6 m深度內(nèi)，每層0.15 m；模型寬0.35 m，共計7列，每列0.05 m。該模型共計840個單元，488個節(jié)點。為了減小計算誤差，取每層各個節(jié)點的算術平均值作為最終計算結果。

2.4　計算方案及邊界條件

土體干密度統(tǒng)一取1.4 g/cm3，根據(jù)3月1號實測土層含水率，計算時初始含水率取12.4%(中效水上限)，計算深度內(nèi)土體溫度統(tǒng)一取15 ℃[17]。

3種計算方案如下：

(1) 1～12周連續(xù)干旱。

(2) 第5周連續(xù)降雨，日均降雨量25.0 mm/d，其他時間段內(nèi)連續(xù)干旱。

(3) 第5周連續(xù)降雨，日均降雨量50.0 mm/d，其他時間段內(nèi)連續(xù)干旱。

上述3種計算方案均采取流量邊界(Neumann type)。降雨入滲時，流量為正，其流量即為日均降雨量，認為所有雨水全被入滲土壤中，忽略地表徑流的影響；蒸發(fā)時，流量為負，其流量可按式(16)確定[18]：

(16)

3計算結果

圖1(a—c)為方案1的計算結果。(1) 連續(xù)干旱條件下，土層實測結果與數(shù)值計算結果較為一致。受蒸發(fā)影響的土層深度為0—1.8 m。這一計算結果與陜北米脂地區(qū)土壤蒸發(fā)規(guī)律相吻合，即蒸發(fā)影響深度為2.0 m，2.0 m以下的土層為相對穩(wěn)定層。

圖1　3-5月土層含水率計算方案1的計算值

同時植被根系也廣泛地分布在0—2.0 m土層中，因此0—2.0 m土層的含水率分布情況對植被的生長起到至關重要的作用。擬合得到了蒸發(fā)影響深度與隨時間的變化關系：

z=1.82-0.74·e-t/3.43

(17)

式中：z——影響深度(m)；t——蒸發(fā)時間(周)。

(2) 0—1.8 m土層的含水率隨著蒸發(fā)的不斷進行而減小，其中0—1.0 m土層的含水率變化最為顯著。連續(xù)3個月的干旱蒸發(fā)使土壤0—1.8 m平均含水率從蒸發(fā)初期的12.4%下降到蒸發(fā)末期的7.9%。當蒸發(fā)進行到第6周(4月中旬)時，0—0.85 m土層已經(jīng)進入難效水的區(qū)間。此時已形成輕度干層，對于根系發(fā)育較淺的植被而言，其生長開始受到阻礙。當蒸發(fā)進行到第8周時，進入難效水的土層厚度為0—1.2 m。當蒸發(fā)進行到第10周時，0—1.4 m的土層水分為難效水。擬合得到0～1.8 m平均含水率隨時間的變化關系：

(18)

(3) 蒸發(fā)強度隨著土壤含水率的下降而減小，同時其蒸發(fā)強度也隨著蒸發(fā)時間的增長而減小。主要原因是土壤中含水率隨著土壤中水分的不斷蒸發(fā)而下降，毛管水開始變的不連續(xù)，同時在表層0—0.1 m形成干層。當表層形成干層后，土壤表層導水率接近0，此時水分只能在干土層以下的表面氣化，并以氣態(tài)水的方式向周圍大氣中擴散。

圖2為方案2的計算結果。雨水入滲深度為0.55 m；降雨后1周內(nèi)(第6周)土壤入滲深度內(nèi)的水分大量蒸發(fā)，當蒸發(fā)進行到第7周(雨后2周)時，0—1.8 m土層平均含水率為9.2%已經(jīng)小于降雨前(第4周)的9.9%。

圖2　3-5月土層含水率計算方案2的計算值

由此表明，上述次降雨條件下裸地土壤水分只能得到短暫的補給，對緩解深層土壤的干旱基本無效，大部分入滲雨水將很快消耗于蒸發(fā)過程中，計算結果表明該蒸發(fā)消耗過程不超過1.5周。然而在以米脂縣為代表的黃土高原旱區(qū)，其春季實際累計平均降雨量只有45 mm(29 a統(tǒng)計結果)遠小于數(shù)值計算中所采用的175 mm降雨量，且其有效降雨時間比較分散，如此氣象條件下，本就稀少的降雨會很快地被蒸發(fā)消耗殆盡，對于緩解黃土高原旱區(qū)春季旱情效果甚微。

圖3為方案3的計算結果。從圖3中可以看出，雨水的入滲深度為1.1 m；雨后2周(第7周)土層0—1 m含水率仍處于中效水；雨后3周(第8周)土壤0—1.8 m平均含水率為10.1%仍高于雨前的9.8%；表明此種類型的降雨條件可以有效地緩解淺層土壤的干旱，尤其對于植被根系廣泛分布于0—1 m土層。

圖3　3-5月土層含水率計算方案3的計算值

4結 論

根據(jù)非飽和黃土二維非穩(wěn)態(tài)流有限元控制方程，對連續(xù)干旱條件下以及不同強度降雨條件下黃土高原旱區(qū)的水分場進行了數(shù)值計算。主要結論如下：對植被生長而言其最不利水分場為連續(xù)干旱條件下；對于黃土高原旱區(qū)而言，強度較小的降雨對于緩解土壤干旱基本無效；大強度集中降雨只能在一段時間內(nèi)緩解土壤旱情。上述結論是基于裸地蒸發(fā)和降雨入滲，對于有植被覆蓋的黃土高原淺層黃土而言，需要進一步研究；針對大強度集中降雨時有地表徑流和地表積水的情況也仍需進一步研究。計算時所采用的邊界條件是基于陜北米脂地區(qū)，對于黃土高原其他地區(qū)的裸地蒸發(fā)和降雨入滲，可以根據(jù)該地區(qū)的實際氣象條件來確定邊界條件。

[1]楊光,孫保平,趙廷寧,等.黃土丘陵溝壑區(qū)退耕還林工程植被恢復效益初步研究[J].干旱區(qū)資源與環(huán)境,2006,20(2):165-170.

[2]何永濤,李文華,李貴才,等.黃土高原地區(qū)森林植被生態(tài)需水研究[J].環(huán)境科學,2004,25(3):35-39.

[3]許炯心.黃土高原植被—降水關系的臨界現(xiàn)象及其在植被建設中的意義[J].生態(tài)學報,2005,25(6):1233-1239.

[4]張建興,馬孝義,趙文舉,等.黃土高原地區(qū)水資源承載力動態(tài)變化分析：以山西、陜西、寧夏、甘肅4省為例[J].干旱區(qū)研究,2009,26(1):115-119.

[5]邵明安,郭忠升,夏永秋,等.黃土高原土壤水分植被承載力研究[M].北京：科學出版社,2010.

[6]Campbell G S. Soil physics with BASIC: transport models for soil-plant systems[M]. Elsevier, 1985.

[7]Yanful E K, Mousavi S M. Estimating falling rate evaporation from finite soil columns[J]. Science of the Total Environment, 2003,313(1):141-152.

[8]Han H, Felker P. Estimation of daily soil water evaporation using an artificial neural network[J]. Journal of Arid Environments, 1997,37(2):251-260.

[9]Ma Ying, Feng Shaoyuan, Su Dongyuan, et al. Modeling water infiltration in a large layered soil column with a modified Green-Ampt model and HYDRUS-1D[J]. Computers and Electronics in Agriculture, 2010,71(S):40-47.

[10]Rao M D, Raghuwanshi N S, Singh R. Development of a physically based 1D-infiltration model for irrigated soils[J]. Agricultural water management, 2006,85(1):165-174.

[11]Gencoglan C, Gencoglan S, Merdun H, et al. Determination of ponding time and number of on-off cycles for sprinkler irrigation applications[J]. Agricultural Water Management, 2005,72(1):47-58.

[12]Lee D H, Abriola L M. Use of the Richards equation in land surface parameterizations[J]. Journal of Geophysical Research: Atmospheres (1984—2012), 1999,104(D22):27519-27526.

[13]楊文治,邵明安.黃土高原土壤水分硏究[M].北京：科學出版社,2000.

[14]Ng C W W, Menzies B. Advanced unsaturated soil mechanics and engineering[M]. CRC Press, 2007.

[15]王鐵行,盧靖,岳彩坤.考慮溫度和密度影響的非飽和黃土土—水特征曲線研究[J].巖土力學,2008,29(1):1-5.

[16]盧靖.非飽和黃土水分遷移問題的試驗研究[D].西安:西安建筑科技大學,2006.

[17]王鐵行,劉自成,岳彩坤.淺層黃土溫度場數(shù)值分析[J].西安建筑科技大學學報:自然科學版,2007,39(4):463-467.

[18]王鐵行,陳晶晶,李彥龍.非飽和黃土地表蒸發(fā)的試驗研究[J].干旱區(qū)研究,2014,31(6):985-990.

[19]Fredlund D G, Rahardjo H. Soil Mechanics for Unsaturated Soils[M]. John Wiley & Sons, 1993.

Numerical Simulation of Moisture Field in Shallow Loess of Loess Plateau During Continous Drought

LI Yanlong, WANG Tiehang, WANG Juanjuan

(SchoolofCivilEngineering,Xi’anUniversityofArchitecture&Technology,Xi’an,Shaanxi710055,China)

Abstract：[Objective] To analyze the dynamic change characteristics of moisture field in shallow loess of Loess Plateau during the continuous drought.[Methods] The dynamic changes of moisture field in shallow loess were calculated numerically by the unsaturated loess unsteady flow 2 dimension finite element control equation under the three different meteorological conditions. [Results] The calculation results were consistent with the measured results during the continuous drought, the thickness of loess affected by the evaporation was 0～1.8 m. Three consecutive months of drought makes its average moisture content dropped to 7.9%, the evaporation intensity of the loess decreased with the decrease of moisture content. The scattered and low intensity rainfall for alleviating continuous drought is invalid basically. The concentrated and intensive rainfall can alleviate soil drought only in a short period.[Conclusion] The results of numerical calculation can provide references for the choice of vegetation types and irrigation time.

Keywords:continuous drought; Loess Plateau; moisture field; numerical simulation

文獻標識碼：B

文章編號：1000-288X(2015)01-0148-05

中圖分類號：TU441

收稿日期：2014-01-10修回日期：2014-01-20

資助項目：國家自然科學基金項目“黃土高原旱區(qū)淺層土體水分場隨氣候變化機理及數(shù)值模擬”(51078309)

第一作者：李彥龍(1985—),男(漢族),河南省南陽市人,博士研究生,研究方向為黃土力學與工程。E-mail:liyanlong1229@163.com。