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

李彥龍, 王鐵行, 王娟娟

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

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

李彥龍, 王鐵行, 王娟娟

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

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

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

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

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

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

1研究區(qū)概況

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

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

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

2.1　計(jì)算模型

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

(1)

式中：[D]——?jiǎng)偠染仃嚕?[E]——容量矩陣； [F]——流量列陣； {h}——水頭列陣。

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

(2)

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

(3)

式中：θ——體積含水率；θs——飽和含水率；θr——?dú)堄嗪剩籬——吸力水頭(m)；α，m，n——系數(shù)，m=1-1/n。

對(duì)黃土而言，式(3) 中相關(guān)參數(shù)可按參考文獻(xiàn)[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)

參考文獻(xiàn)[16]給出了非飽和黃土水分?jǐn)U散率與含水率之間的關(guān)系：

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

(10)

式中：D——非飽和黃土的水分?jǐn)U散率；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ù)與擴(kuò)散率和比水容有如下關(guān)系：

kw=D×C

(12)

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

(13)

2.2　水頭的確定

土體中的水頭通常可以表達(dá)為：

h=hg+hm+hp+hs+hT

(14)

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

溫度水頭和溶質(zhì)水頭通?？梢院雎?，對(duì)非飽和土中的水分遷移而言，認(rèn)為土體所有孔隙與大氣是聯(lián)通的，因此壓力水頭為零。因此只需考慮重力水頭和基質(zhì)吸力水頭，式(14)可以簡(jiǎn)化為：

h=hg+hm

(15)

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

2.3　有限元模型

計(jì)算深度為6.0 m，所采用的三節(jié)點(diǎn)三角形單元均為直角三角形，共計(jì)劃分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，共計(jì)7列，每列0.05 m。該模型共計(jì)840個(gè)單元，488個(gè)節(jié)點(diǎn)。為了減小計(jì)算誤差，取每層各個(gè)節(jié)點(diǎn)的算術(shù)平均值作為最終計(jì)算結(jié)果。

2.4　計(jì)算方案及邊界條件

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

3種計(jì)算方案如下：

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

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

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

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

(16)

3計(jì)算結(jié)果

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

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

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

z=1.82-0.74·e-t/3.43

(17)

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

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

(18)

(3) 蒸發(fā)強(qiáng)度隨著土壤含水率的下降而減小，同時(shí)其蒸發(fā)強(qiáng)度也隨著蒸發(fā)時(shí)間的增長(zhǎng)而減小。主要原因是土壤中含水率隨著土壤中水分的不斷蒸發(fā)而下降，毛管水開(kāi)始變的不連續(xù)，同時(shí)在表層0—0.1 m形成干層。當(dāng)表層形成干層后，土壤表層導(dǎo)水率接近0，此時(shí)水分只能在干土層以下的表面氣化，并以氣態(tài)水的方式向周?chē)髿庵袛U(kuò)散。

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

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

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

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

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

4結(jié) 論

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

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

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

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

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

[5]邵明安,郭忠升,夏永秋,等.黃土高原土壤水分植被承載力研究[M].北京：科學(xué)出版社,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].北京：科學(xué)出版社,2000.

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

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

[16]盧靖.非飽和黃土水分遷移問(wèn)題的試驗(yàn)研究[D].西安:西安建筑科技大學(xué),2006.

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

[18]王鐵行,陳晶晶,李彥龍.非飽和黃土地表蒸發(fā)的試驗(yàn)研究[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

文獻(xiàn)標(biāo)識(shí)碼：B

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

中圖分類號(hào)：TU441

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

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

第一作者：李彥龍(1985—),男(漢族),河南省南陽(yáng)市人,博士研究生,研究方向?yàn)辄S土力學(xué)與工程。E-mail:liyanlong1229@163.com。