土壤侵蝕臨界切應(yīng)力的計(jì)算研究評述

(汕頭大學(xué)環(huán)境與土木工程系 廣東 汕頭 515063)

一、前言

土壤侵蝕一直都是世界范圍內(nèi)影響人類生存和發(fā)展環(huán)境的重要問題之一，由于土壤侵蝕導(dǎo)致的水土流失已成為限制許多地方經(jīng)濟(jì)發(fā)展的主要障礙。從20世界30年代開始，國內(nèi)外的學(xué)者們對土壤可蝕性的研究經(jīng)歷了從研究某種單一因素[1]～[3]到確立可蝕性因子(K)評價(jià)指標(biāo)[4]，最終建立圖解法(諾謨圖)[5]、經(jīng)驗(yàn)土壤理化性質(zhì)測定法[6]或各類預(yù)測模型公式[7]～[9]并在實(shí)際中加以運(yùn)用的過程，經(jīng)典的USLE模型[7]、RUSE模型[8]、EPIC模型[9]等都是對于該問題的重要研究成果。

隨著眾多學(xué)者不斷深入地研究侵蝕機(jī)理及計(jì)算機(jī)技術(shù)的快速進(jìn)步，土壤侵蝕預(yù)報(bào)研究也由一般的因子相關(guān)分析向以土壤侵蝕發(fā)生過程為基礎(chǔ)的過程模型發(fā)展。人們以泥沙動(dòng)力學(xué)及河床演變理論為依據(jù)，分析侵蝕發(fā)生過程中各因子間的力學(xué)關(guān)系，用較為嚴(yán)密的控制方程計(jì)算土壤侵蝕量，出現(xiàn)了一些具有一定物理意義的土壤侵蝕預(yù)報(bào)模型和實(shí)驗(yàn)裝置[10]～[12]。

本文從土壤侵蝕模型的重要參數(shù)之一——臨界剪切應(yīng)力τc著手，介紹了目前侵蝕測算和泥沙動(dòng)力學(xué)中對該參數(shù)的不同理論與計(jì)算方法，綜合分析現(xiàn)有方法的不足，并指出可能的研究方向及重點(diǎn)領(lǐng)域，以期能為今后的土壤水動(dòng)力學(xué)研究起到一定的借鑒和幫助。

二、臨界切應(yīng)力的概念及其計(jì)算

對于土壤，侵蝕的過程就是土壤顆粒在外營力(水力、風(fēng)力、融凍等)的作用下，分離、搬運(yùn)、沉積的過程[13]。當(dāng)土顆粒與土體開始分離的臨界狀態(tài)時(shí)，由水體施加給土體的切應(yīng)力即是臨界切應(yīng)力，以τc表示。

Shields[14]在分析了床面上泥沙顆粒的受力平衡條件后，推導(dǎo)出無粘性均勻沙的臨界切應(yīng)力公式，并在此基礎(chǔ)上得出無因次切應(yīng)力Θc與沙粒雷諾數(shù)Re*的關(guān)系曲線，其值范圍為0.040～0.060[14]。后來學(xué)者在Shields的研究基礎(chǔ)上陸續(xù)提出了一些改進(jìn)型經(jīng)驗(yàn)公式，比如Briaud[10]在這基礎(chǔ)上通過力學(xué)分析和實(shí)驗(yàn)，認(rèn)為τc與平均粒徑D50成線性關(guān)系，為簡化運(yùn)算在應(yīng)用中可將τc近似取值為平均粒徑D50；美國80年代開始著手研發(fā)的被稱為WEPP(Water Erosion Prediction Project)[11]的新一代土壤水蝕預(yù)報(bào)模型最早將臨界剪切力(τc)作為直接的參數(shù)之一，在模型中使用時(shí)通過簡單的公式推算取值；Hanson等[12]的JETs實(shí)驗(yàn)也提出了基于其實(shí)驗(yàn)裝置對τc的經(jīng)驗(yàn)計(jì)算方式。

(一)EFA實(shí)驗(yàn)

Briaud在文章中指出，只有在水流切應(yīng)力大于高于臨界剪切應(yīng)力時(shí)土顆粒才會(huì)發(fā)生分離并能觀測到一定的沖刷速率[15]，在此分析基礎(chǔ)上，Briaud提出將與標(biāo)準(zhǔn)化的小侵蝕速率相對應(yīng)的剪應(yīng)力作為臨界剪切應(yīng)力τc的一個(gè)可接受的新定義。他定義使用erosion function apparatus(EFA)裝置進(jìn)行為期一周的水波沖刷而僅造成168mm侵蝕時(shí)，1mm/hr的速率即是該侵蝕速率標(biāo)準(zhǔn)閾值[15]。

Briaud[16]以理論分析砂和礫石的情況下，沖刷過程是由顆粒的重量控制的，計(jì)算τc可看做直接計(jì)算克服兩個(gè)顆粒間摩擦力所需的水流切應(yīng)力。水平方向受力平衡方程為

τcAe=Wtanφ

其中Ae為水流作用在顆粒上的有效面積；W為土顆粒在水中的重量；φ為兩顆粒間的界面摩擦角。若將土壤顆粒視作球形(圖3(a))，公式可寫作

或

式中α=有效摩擦面積與球形顆粒最大橫截面之比，約6.14；D50為代表土壤粒徑分布的平均粒徑，mm；ρs、ρw為土顆粒的密度，水的密度，g/cm3。

圖3 沖刷過程中作用于土顆粒上的力

對于情形如圖3(b)時(shí)，水對土壤顆粒所施加的合力是平行于侵蝕表面的剪切力，相鄰的顆粒不會(huì)阻礙這一過程，旋轉(zhuǎn)是圍繞著與底層顆粒的接觸點(diǎn)進(jìn)行的。由于分析是針對砂顆粒進(jìn)行的，因此忽略了顆粒間的電磁力和靜電力。在開始運(yùn)動(dòng)時(shí)，接觸點(diǎn)O附近的力矩平衡

τcAea=Wb

或化簡可得

對滑動(dòng)和滾動(dòng)機(jī)制的簡單分析有助于闡明影響粗粒土初期運(yùn)動(dòng)的重要因素，可以看出臨界剪應(yīng)力τc與平均粒徑成線性關(guān)系，且其比例系數(shù)可能取決于相對密度。這與希爾茲[14]對平坦砂層做一系列水槽實(shí)驗(yàn)后所畫出的結(jié)果圖形(the Shields diagram)相一致。通過對于實(shí)驗(yàn)數(shù)據(jù)的對比分析，在粗粒土中臨界剪應(yīng)力τc可以經(jīng)驗(yàn)性地視作平均粒徑D50。

EFA實(shí)驗(yàn)是一種簡單高效的測試方法，可用標(biāo)準(zhǔn)的取樣進(jìn)行原位測試，但應(yīng)用范圍僅限于侵蝕率在1mm/d(軟巖)至數(shù)m/h(干凈細(xì)沙)范圍且流速0.1m/s—0.6m/s的情況，實(shí)驗(yàn)給出的臨界切應(yīng)力對應(yīng)于運(yùn)動(dòng)的開始和超過這一點(diǎn)的侵蝕率，侵蝕率切(mm/h)切應(yīng)力(N/m2)的相對誤差在10%左右，并不適用于精度要求較高的侵蝕測量或預(yù)測使用。

(二)WEPP模型

WEPP模型中，將土壤可蝕性進(jìn)一步劃分為細(xì)溝間(interrill)可蝕性(ki)、細(xì)溝(rill)可蝕性(kr)和土壤臨界切應(yīng)力τc。其中的細(xì)溝侵蝕，就是指由于水流的沖刷力大于臨界剪切力而造成的土體分離[17]。Gilley[18]等人通過實(shí)驗(yàn)統(tǒng)計(jì)了細(xì)溝侵蝕與剪切應(yīng)力之間的經(jīng)驗(yàn)數(shù)據(jù)，在對臨界剪切應(yīng)力和土壤特性進(jìn)行逐步多元回歸分析后，發(fā)現(xiàn)二者具有顯著的線性相關(guān)。

對于粘土含量小于7.5%的土壤(soil with water dispersible clay content of less than 7.5%)，

τc=0.216Cla-183(coefficientoflinearextensibility)+0.412ω1.5MPa+0.780

式中Cla為粘粒含量,%；ω1.5MPa為1.5MPa條件下土壤的含水率，%。

對于粘土含量大于等于7.5%的土壤(soils with water dispersible clay content of 7.5% or greater)，

τc=0.296CECCa+1.53Fe+7.75C-11.4CECK-0.535Sand-0.208

式中CECCa為鈣含量，cmol/kg；Fe為鐵含量，%；C為有機(jī)碳含量，%；CECK為鉀含量，cmol/kg；Sand 為細(xì)沙含量，%。

這種以兩個(gè)經(jīng)驗(yàn)公式推算不同土壤臨界切應(yīng)力τc的方法并不是一種基于機(jī)理分析而是實(shí)驗(yàn)數(shù)據(jù)分析。事實(shí)上為了得到WEPP模型的可蝕性參數(shù)，1987年和1988年[19]在全美24個(gè)州做了2年人工降雨和放水沖刷試驗(yàn)，但結(jié)果并不理想。因?yàn)樵囼?yàn)得到的土壤流失量難于區(qū)分是分離的結(jié)果還是搬運(yùn)的結(jié)果。

此后，雷廷武[20]等曾在理論分析的基礎(chǔ)上設(shè)計(jì)了一種測量臨界剪切應(yīng)力的實(shí)驗(yàn)裝置，用于測量沉積疏松土壤中細(xì)溝再生的力學(xué)參數(shù)并得到與溝坡相關(guān)的臨界剪切應(yīng)力測量值。張晴雯[21]利用室內(nèi)細(xì)溝侵蝕模擬沖刷實(shí)驗(yàn)，由實(shí)驗(yàn)來直接計(jì)算細(xì)溝侵蝕可蝕性參數(shù)及土壤臨界抗剪切應(yīng)力。其認(rèn)為，對于細(xì)溝侵蝕產(chǎn)沙數(shù)學(xué)模型[22]：

式中Dr為細(xì)溝剝蝕分散率，kg·m-2·s-1；Kr為單位寬度上的細(xì)溝可蝕性參數(shù)，kg·N-2·s-1；τ、τc為水流剪切應(yīng)力和土壤臨界剪切應(yīng)力，N·m-2；q為流量，m2·s-1；c為泥沙含量，kg·m-3；Tc為水流的輸沙能力，kg·m-1·s-1。

當(dāng)水流中的含沙量c=0時(shí)，細(xì)溝侵蝕過程具有最大可能剝蝕分散率Dr，max

Dr,max=Kr(τ-τc)

可以看出概念上WEPP基于物理過程，且在計(jì)算侵蝕時(shí)幾乎涵蓋了影響土壤侵蝕的各種因素，但一些過程的表達(dá)仍是經(jīng)驗(yàn)公式，對侵蝕機(jī)理仍需要進(jìn)一步研究和驗(yàn)證。WEPP中所需輸入的參數(shù)τc的確定尚必須采用非測量的“估算”方法，多數(shù)情況下需要采用經(jīng)驗(yàn)數(shù)據(jù)對所得參數(shù)進(jìn)行率定和修正。因此，在加強(qiáng)典型區(qū)域的侵蝕研究的同時(shí)，還需要強(qiáng)化土壤侵蝕模型中參數(shù)區(qū)域規(guī)律的研究，建立多重尺度上的模型參數(shù)與其相關(guān)環(huán)境因子的關(guān)系。

(三)JET實(shí)驗(yàn)

JET(Jet Erosion Test)是一種可用于現(xiàn)場和實(shí)驗(yàn)室評估土壤可蝕性的通用實(shí)驗(yàn)[23]。Hanson等認(rèn)為土壤侵蝕速率ε與有效剪切應(yīng)力τe和臨界剪切應(yīng)力τc的差值成正相關(guān)[24]。

ε=Kd(τe-τc)

式中ε為可蝕性系數(shù)，m3/N-s；τe為有效剪切應(yīng)力，Pa；τc為臨界剪切應(yīng)力，Pa。

圖2 JET實(shí)驗(yàn)示意圖與參數(shù)定義

實(shí)驗(yàn)中，臨界剪切應(yīng)力τc主要是根據(jù)平衡沖刷深度Je確定的，其表達(dá)式為:

式中τ0為噴嘴處射流速度引起的最大應(yīng)力，Pa；Jp為等速核長度，m；Jp為平衡沖刷深度，m；U0為噴嘴處的速度，m/s；g為重力加速度常數(shù)，m/s2；h為水頭差，m。通過實(shí)驗(yàn)測量所需數(shù)據(jù)，代入式中即得出侵蝕的臨界剪切應(yīng)力[24]。

綜上所述，不論是哪種預(yù)測模型，在計(jì)算土壤侵蝕率或可蝕性時(shí)臨界切應(yīng)力τc作為一種計(jì)算可蝕性的必要參數(shù)，在對其進(jìn)行處理時(shí)，都是以實(shí)驗(yàn)回歸分析和實(shí)驗(yàn)經(jīng)驗(yàn)分析得出，并簡化了坡度、顆粒均勻性等因素的影響。最終計(jì)算出的土壤可蝕性與真實(shí)結(jié)果有著一定的誤差，需要在每次不同的實(shí)際應(yīng)用時(shí)進(jìn)行相應(yīng)的模型修正。

三、臨界切應(yīng)力的理論分析計(jì)算方法

在泥沙動(dòng)力學(xué)中，泥沙顆粒由靜止進(jìn)入運(yùn)動(dòng)狀態(tài)的過程稱為泥沙的起動(dòng)[25]，是泥沙動(dòng)力學(xué)研究的基本問題之一，對于臨界切應(yīng)力τc又被稱為起動(dòng)切應(yīng)力[26]，在國內(nèi)外學(xué)者中都有廣泛的研究。

Engelund[27]在研究水流阻力時(shí)得到的規(guī)律

Θ′=0.06+0.4Θ2

其中Θ′與Θ分別是沙粒阻力與總阻力有關(guān)的無因次切力，隨著Θ的減小，Θ′趨向于定值0.06，因此Θc=0.06可認(rèn)為相當(dāng)于泥沙的起動(dòng)條件，此時(shí)的水流切應(yīng)力即為臨界切應(yīng)力τc。

葉基阿扎洛夫[28]對Shields的結(jié)果作了論證，建立顆粒的平衡條件，并假設(shè)泥沙起動(dòng)時(shí)，顆粒的作用流速等于顆粒的自由沉降速度，公式變換可得

對于粗泥沙，粗糙紊流，χ=1，顆粒繞流阻力系數(shù)CD=0.4。取Ks=d，y=0.63d，代入上式可得Θc=0.0616?0.06。

愛因斯坦在研究推移質(zhì)輸沙率時(shí)[29]，將上舉力看做隨機(jī)變量，泥沙起跳概率即是上舉力FL大于顆粒重力W的概率P(W/FL<1)，顆粒的起跳條件公式可改寫為

τc=0.715(γs-γ)d

當(dāng)gb=0時(shí)或φb=0時(shí)，得θ′=Θc=0.047。這一結(jié)果已為許多研究者所引用，如格斯勒在研究非均勻沙的起動(dòng)條件時(shí)[31]，取個(gè)別顆粒的起動(dòng)條件是Θc=0.047；萬兆惠在研究泥沙起動(dòng)條件時(shí)也采用了這一結(jié)果[32]。

相對于國外的研究，國內(nèi)學(xué)者起步較晚。褚君達(dá)[33]對國內(nèi)外數(shù)十種不同形式的泥沙起動(dòng)公式進(jìn)行了分析，統(tǒng)一轉(zhuǎn)化為無量綱起動(dòng)切應(yīng)力Θc=τc/(γs-γ)D(τc為水流施加的臨界切應(yīng)力，γs和γ分別為泥沙、水的容重，D為泥沙粒徑)的形式，并計(jì)算得出無粘性均勻沙的Θc為0.0231～0.0716；何文社[34]從推移質(zhì)泥沙起動(dòng)的特點(diǎn)出發(fā)，引入附加質(zhì)量力，采用滾動(dòng)起動(dòng)模型推導(dǎo)出了泥沙起動(dòng)無因次切應(yīng)力起動(dòng)公式，結(jié)果表明弱動(dòng)(個(gè)別運(yùn)動(dòng))情況下無因次切應(yīng)力為0.0291～0.0582，中動(dòng)(少量運(yùn)動(dòng))的情況下為0.353～0.706；

解剛等[35]考慮了水流條件及床沙組成對泥沙起動(dòng)的影響，采用滾動(dòng)起動(dòng)模型推導(dǎo)出了非均勻沙分級起動(dòng)的切應(yīng)力公式:

式中di為分級粒徑，mm；α為系數(shù)，可取α=2/3；η為床面起動(dòng)流速與切應(yīng)力間的關(guān)系系數(shù)，可取η=7.055；dm為平均粒徑。

黃偉[36]等以長江口取樣的細(xì)顆粒泥沙為樣本進(jìn)行了環(huán)形水槽起動(dòng)試驗(yàn)，給出3組不同中值粒徑泥沙的起動(dòng)流速和臨界起動(dòng)切應(yīng)力公式，同時(shí)用湍動(dòng)能法估算了不同含沙量水體的床面切應(yīng)力，給出了泥沙起動(dòng)過程中床面切應(yīng)力與含沙量之間的關(guān)系式，得出中值粒徑為0.082mm、0.035mmhe 0.008mm的泥沙起動(dòng)流速分別為0.29m/s、0.55m/s和0.78m/s，臨界切應(yīng)力分別為0.19Pa、0.34Pa、0.46Pa。可見中值粒徑在0.1mm以下的細(xì)顆粒泥沙，臨界切應(yīng)力和起動(dòng)流速隨著粒徑的增大反而減小(即粒徑越大，越容易起動(dòng))。

四、結(jié)論與展望

(1)目前在計(jì)算土壤侵蝕時(shí)，工程應(yīng)用中對于土壤侵蝕的臨界切應(yīng)力計(jì)算主要是利用經(jīng)驗(yàn)公式或?qū)嶒?yàn)?zāi)Ｐ瓦M(jìn)行，因此方法本身存在一定的誤差，對于如何找到更準(zhǔn)確表征臨界切應(yīng)力的計(jì)算方法將是一項(xiàng)長期任務(wù)。

(2)相對于各類侵蝕預(yù)測模型的經(jīng)驗(yàn)計(jì)算方法，泥沙動(dòng)力學(xué)中對于臨界切應(yīng)力的研究更注重侵蝕過程的機(jī)理分析，但主要是借鑒明渠水流參數(shù)進(jìn)行。由于坡面侵蝕發(fā)生條件及過程復(fù)雜，徑流特征與明渠水流存在較大差異，研究結(jié)果帶有一定的不確定性。未來研究中影響因素設(shè)定有復(fù)雜化、全面化的趨勢，從而以更準(zhǔn)確的徑流水動(dòng)力機(jī)制的角度解釋土壤侵蝕的發(fā)生和發(fā)展。

(3)目前大多研究所獲得的結(jié)論都是整個(gè)坡面的平均狀況，對于復(fù)雜坡面侵蝕演化過程中的分布規(guī)律及其特性研究十分稀少，是今后研究中的重要內(nèi)容。

(4)加強(qiáng)解析計(jì)算方法與實(shí)際工程的結(jié)合，特別是應(yīng)加強(qiáng)土壤侵蝕水動(dòng)力學(xué)特性的原位研究，以更好地揭示土壤侵蝕水動(dòng)力特性的變化規(guī)律。