低覆蓋度下單個(gè)粗糙元保護(hù)作用的變化特征

摘要：圓柱體和長(zhǎng)方體是地表風(fēng)蝕相關(guān)研究中最常見(jiàn)的2種粗糙元.基于已有關(guān)于這2種粗糙元流向風(fēng)速變化模型，推導(dǎo)了低覆蓋度下均勻分布的單個(gè)粗糙元的保護(hù)區(qū)域最大流向長(zhǎng)度、保護(hù)區(qū)域面積和保護(hù)效率的表達(dá)式，并考察和比較了采用這2種粗糙元時(shí)，這些物理量隨面率和粗糙密度的變化規(guī)律.結(jié)果表明，隨著粗糙密度的增加，單個(gè)粗糙元保護(hù)區(qū)域最大流向長(zhǎng)度先增加后減小，最后逐漸趨于穩(wěn)定；粗糙元保護(hù)效率隨寬度呈單調(diào)減小后逐漸趨于穩(wěn)定，而隨著高度的增加呈現(xiàn)先增大后減小而后趨于穩(wěn)定；相同面率和粗糙密度情況下，長(zhǎng)方體粗糙元的保護(hù)效率比圓柱體粗糙元的更高.這些結(jié)果可為防治地表侵蝕提供理論指導(dǎo).

關(guān)鍵詞：面率；粗糙密度；最大保護(hù)長(zhǎng)度；保護(hù)效率；保護(hù)面積

中圖分類(lèi)號(hào)：X513文獻(xiàn)標(biāo)志碼：A

0引言

大氣邊界層中的植被（樹(shù)與灌木）和建筑物2種主要地表粗糙元，不僅可以通過(guò)減少地表裸露程度，而且還可以通過(guò)對(duì)其周?chē)鲌?chǎng)的影響降低其下風(fēng)向一定區(qū)域內(nèi)的風(fēng)速，來(lái)降低粉塵釋放.基于形狀特征，植被通常被簡(jiǎn)化為圓柱體，而建筑物通常被簡(jiǎn)化為長(zhǎng)方體粗糙元.目前已有許多理論、實(shí)驗(yàn)和數(shù)值研究對(duì)這2種粗糙元的影響隨覆蓋密度和分布規(guī)律等的變化進(jìn)行了討論［1-17］.不僅對(duì)如何利用和設(shè)計(jì)人工形成的草垛來(lái)降低風(fēng)沙輸運(yùn)［18］進(jìn)行了研究，而且還提高了對(duì)城市顆粒物質(zhì)擴(kuò)散的認(rèn)識(shí)［19-20］，更重要的是，這些研究結(jié)果可以為風(fēng)蝕預(yù)報(bào)模式提供重要參數(shù)，比如剪切力分配模式參數(shù).

盡管方法在不斷改進(jìn)，基于剪切力分配的風(fēng)蝕預(yù)報(bào)模式仍然會(huì)產(chǎn)生不確定性.在數(shù)值預(yù)報(bào)中，沙塵濃度估計(jì)不準(zhǔn)確仍然是一個(gè)常見(jiàn)的問(wèn)題［21］.研究表明，采用剪切力分配模式使得估計(jì)得到的侵蝕量與實(shí)驗(yàn)測(cè)量之間存在量級(jí)上的誤差［10］.這也表明剪切力分配模式的不完善可能是沙塵大尺度預(yù)報(bào)中出現(xiàn)有些地方并沒(méi)有發(fā)生沙塵事件但預(yù)報(bào)中有，有的地方發(fā)生了沙塵事件預(yù)報(bào)卻沒(méi)有的一個(gè)重要原因.

即使經(jīng)過(guò)幾十年的研究，粗糙元對(duì)流場(chǎng)的相互作用仍沒(méi)有解釋清楚［22-24］，也不能通過(guò)某一物理統(tǒng)計(jì)量準(zhǔn)確描述粗糙元對(duì)地表的保護(hù)作用.因此，本文將基于已有研究［1，15-16］，比較圓柱體粗糙元和長(zhǎng)方體粗糙元在不同情況下對(duì)地表的保護(hù)作用.一方面可以考察已有理論在陰影保護(hù)區(qū)域上的差異，從而為降低剪切力分配帶來(lái)的誤差有所幫助；另一方面也可以提供一種新的思路和方法，為數(shù)值模式中的地表參數(shù)化方案提供更多參考，從而使得數(shù)值模擬更準(zhǔn)確.

1研究方法

粗糙元的形狀和排列密度的不同將使得周?chē)鲌?chǎng)分布發(fā)生改變，從而在粗糙元下流方向形成不同范圍的速度降低區(qū)域.本研究介紹了在粗糙元均勻分布時(shí)不同密度和來(lái)流風(fēng)速下投影形狀為圓形的圓柱體粗糙元和投影形狀為正方形的長(zhǎng)方體粗糙元形成的保護(hù)區(qū)域的估算方法.

1.1基本概念和定義

圖1給出了圓柱體和長(zhǎng)方體粗糙元形成的陰影保護(hù)區(qū)域投影形狀示意圖.Raupach［1］假設(shè)圓柱體下風(fēng)向的陰影保護(hù)區(qū)域形狀為三角形，如圖1（A）所示.對(duì)于長(zhǎng)方體粗糙元，其陰影保護(hù)區(qū)域形狀與面率（AR，高度與寬度之比）相關(guān).當(dāng)面率lt;1時(shí)，長(zhǎng)方體粗糙元形成的陰影保護(hù)區(qū)域形狀可以假設(shè)為梯形［15］，如圖1（B）所示；當(dāng)面率≥1時(shí)，陰影保護(hù)區(qū)域形狀可由梯形和三角形構(gòu)成［16］，如圖1（C）所示.下標(biāo)R，Y和S分別對(duì)應(yīng)Raupach［1］、Yang［15］和Sadique等［16］的工作.As，R、 As，Y 和 As，S為對(duì)應(yīng)研究工作給出的保護(hù)區(qū)域面積， Lx，R、 Lx，Y 和 Lx，S 為對(duì)應(yīng)研究工作給出的最大流向保護(hù)長(zhǎng)度.Lx，S由le和lc 2部分組成，hy （hy =h）和 lw 為展向的擴(kuò)展寬度.

在進(jìn)行進(jìn)一步分析和討論之前，需要給出一些基本定義.對(duì)于一個(gè)粗糙元，其對(duì)地表總的保護(hù)區(qū)域面積（A0）由粗糙元投影面積（Ab）和陰影保護(hù)區(qū)域面積（As）2部分組成.對(duì)于圓柱體粗糙元，其投影面積為πw2/4，以符號(hào)Ab，cir表示；而對(duì)于長(zhǎng)方體粗糙元，其投影面積為w2，以符號(hào)Ab，sq表示.對(duì)于地表粗糙程度的描述，常用地表粗糙覆蓋度λc和地表粗糙密度λd來(lái)表示.根據(jù)圖2，圓柱體粗糙元地表粗糙覆蓋度定義為πw2/（4l2），長(zhǎng)方體粗糙元地表粗糙覆蓋度定義為w2/l2.而地表粗糙密度則定義為單位地表面積上的粗糙元迎風(fēng)面積.高度為h的圓柱體和長(zhǎng)方體粗糙元的地表粗糙密度均為wh/l2.

當(dāng)?shù)乇泶植谠母采w度逐漸增加，單個(gè)粗糙元形成的保護(hù)區(qū)域形狀和面積可能發(fā)生變化，如圖3所示.由于實(shí)際地表粗糙元的排列隨機(jī)且復(fù)雜，本研究主要定量地研究均勻分布的粗糙元的影響.

1.2圓柱體粗糙元陰影保護(hù)區(qū)域投影面積

基于Raupach［1］的假設(shè)，如果相鄰粗糙元之間沒(méi)有相互影響，即粗糙元密度很低時(shí)，單個(gè)粗糙元形成的總的保護(hù)區(qū)域投影面積A0，R為，

A0，R=c*whUh/u*（1）

式中，c*=c′*（w/h）f，c′*取值為1，而f是w/h的函數(shù)［1］，f=exp（-w/h）.

由公式（1）可得到單個(gè)粗糙元的陰影保護(hù)區(qū)域下風(fēng)向長(zhǎng)度Lx，R為，

Lx，R=2（AO，R-0.5Ab，cir）/w（2）

此時(shí)粗糙元密度足夠低，Lx，R≤l-w/2成立，粗糙元陰影保護(hù)區(qū)域投影面積與底面覆蓋面積的比值（即為保護(hù)效率P0）為，

P0，R=AsAb=A0，R-Ab，cirAb，cir=4c*hUhπwu*-1（3）

P0，R=AsAb=A0，R-Ab，cirAb，cir=4c*hUhπwu*-4θR（Lx，R-l+w/2）2πw2-1（4）

隨著粗糙元密度的增加，相鄰粗糙元之間的距離變小，即有Lx，Rgt;l-w/2，此時(shí)單個(gè)粗糙元的實(shí)際陰影保護(hù)區(qū)域下風(fēng)向長(zhǎng)度為L(zhǎng)1x，R=l-w/2.此時(shí)減少的保護(hù)區(qū)域面積實(shí)際上為扇形，其對(duì)應(yīng)半徑長(zhǎng)為L(zhǎng)*=Lx，R-L1x，R=Lx，R=（l-w/2）.因而，總的實(shí)際保護(hù)面積變?yōu)锳10，R=A0，R-θRL2*=A0，R-（Lx，R-l+w/2）2，相應(yīng)的保護(hù)效率則可以由公式（4）給出，其中θR=tan-1（u*/Uh）.考慮到由經(jīng)驗(yàn)公式得到的準(zhǔn)確度（尤其是在粗糙密度較大時(shí)），本研究實(shí)際操作中取近似角度θR=tan-1（0.5w/Lx，R）.

γ=Uh/u*=exp（cγλd/2）/Cs+λdCr（5）

這里，定義γ=Uh/u*，按照Raupach［1］的推導(dǎo)，當(dāng)?shù)乇泶植诿芏鹊陀讦薲，m時(shí)，γ與地表粗糙密度之間的關(guān)系可以表達(dá)為公式（5）；而大于這個(gè)值時(shí)，γ的值保持不變.其中，λd，m=（cr+c2r+4cs（ec/2）2/［2（ec/2）2］，cr=0.3，cs=0.003，c=0.37，e≈2.718.

1.3長(zhǎng)方體粗糙元陰影保護(hù)區(qū)域投影面積

相比于圓柱體粗糙元陰影保護(hù)區(qū)域三角形態(tài)對(duì)所有面率都適用，長(zhǎng)方體粗糙元陰影保護(hù)區(qū)域形態(tài)隨著面率的改變會(huì)發(fā)生變化.這里根據(jù)已有研究主要分高和低面率2種情況討論.

1.3.1低面率情況

面率低于1的時(shí)候?yàn)榈兔媛是闆r.在這種情況下，陰影保護(hù)區(qū)域投影為如圖1（B）所示的梯形.如果相鄰粗糙元之間沒(méi)有相互影響（粗糙元密度很低），單個(gè)粗糙元形成的陰影保護(hù)區(qū)域下風(fēng)向長(zhǎng)度Lx，Y=h/tanθY［15］，其中，θY=tan-1（cθu*/Uh），而cθ=1/3+2h/3w.因而，最終的陰影保護(hù)區(qū)域下風(fēng)向長(zhǎng)度為，

Lx，y=hUh（1/3+2h/3w）u*（6）

單個(gè)長(zhǎng)方體粗糙元形成的總的保護(hù)區(qū)域最大面積As，Y為，

As，Y=（w+h）Lx，Y=（w+h）hUh（1/3+2h/3w）u*（7）

當(dāng)Lx，Y≤l-w成立時(shí)，單個(gè)粗糙元的保護(hù)效率為，

P0，Y=AsAb=As，YAb，sq=（w+h）hUh（1/3+2h/3w）w2u*（8）

當(dāng)l-wlt;Lx，Y≤l時(shí)，實(shí)際陰影保護(hù)區(qū)域下風(fēng)向長(zhǎng)度為L(zhǎng)1x，Y=l-w.減少的保護(hù)區(qū)域面積實(shí)際上為矩形，其寬度為w，而沿流向的尺度為L(zhǎng)*=Lx，Y-L1x，Y=Lx，Y-（l-w）.總的實(shí)際保護(hù)面積變?yōu)锳1s，Y=As，Y-w（Lx，Y-l+w），那么保護(hù)效率則為，

P0，Y=AsAb=A1s，YAb，sq=h2Uh（1/3+2h/3w）w2u*+lw-1（9）

而當(dāng)Lx，Ygt;l時(shí)，下風(fēng)向的粗糙元會(huì)整體浸入上風(fēng)向粗糙元的陰影保護(hù)區(qū)域中，這樣2個(gè)粗糙元形成的陰影保護(hù)區(qū)域就會(huì)發(fā)生重疊.那么，經(jīng)過(guò)簡(jiǎn)單幾何分析，得到此時(shí)單個(gè)粗糙元總的實(shí)際保護(hù)面積為A2s，Y=（l-w）w+h［2l-（1/3+2h/3w）l2u*/（hUh）］，這樣保護(hù)效率則為，

P0，Y=AsAb=A2s，YAb，sq=

h［2l-（1/3+2h/3w）l2u*/（hUh）］w2+l-ww（10）

1.3.2高面率情況

當(dāng)面率大于1時(shí)， Sadique 等［16］認(rèn)為影響陰影保護(hù)區(qū)域形態(tài)的主要因素不再是低面率下的高度，而是粗糙元的寬度.此時(shí)，長(zhǎng)方體粗糙元形成的保護(hù)區(qū)域由梯形和三角形共同構(gòu)成（見(jiàn)圖1（C）），而單個(gè)粗糙元形成的陰影保護(hù)區(qū)域下風(fēng)向長(zhǎng)度應(yīng)該是矩形和三角形沿風(fēng)向長(zhǎng)度尺寸之和.

Lx，S=le+lc=2.5 le=2.5 wUh/u*（11）

式中，le=wUh/u*，lc=1.5 le

單個(gè)粗糙元形成的最大保護(hù)區(qū)域面積為，

As，S=（w+lw）le+（w+2lw）lc/2=4.25wle=4.25w2Uh/u*（12）

式中，lw=w.

當(dāng)Lx，S≤l-w成立時(shí)，保護(hù)效率可表述為，

P0，S=AsAb=As，SAb，sq=4.25Uhu*（13）

相比于低面率，在高面率情況下，當(dāng)Lx，Sgt;l-w成立時(shí)，陰影保護(hù)面積減少可能的理論情況存在多樣性.如果l-wlt;Lx，S≤lc/3+l-w，實(shí)際的陰影面積保護(hù)面積減小為A1s，S=5w（l-w）-2w2Uh/u*-（l-w）2u*/Uh，相應(yīng)的保護(hù)效率為，

P0，S=AsAb=A1s，SAb，sq=5（l-w）/w-2Uh/u*-

（l-w）2u*/w2Uh（14）

而當(dāng)lc/3+l-wlt;Lx，S≤lc/3+l時(shí)，單個(gè)粗糙元在下風(fēng)向形成的實(shí)際的陰影保護(hù)面積減小為A2s，S=w（l-w）+2w2Uh/u*.此時(shí)，保護(hù)效率為，

P0，S=AsAb=A2s，SAb，sq=（l-w）/w+2Uh/u*（15）

而當(dāng)Lx，Sgt;lc/3+l時(shí)，單個(gè)粗糙元在下風(fēng)向形成的實(shí)際的陰影保護(hù)面積減小為A3s，S=w（l-w）+2w2Uh/u*-lw［2-lu*/（2wUh）］.那么，保護(hù)效率為，

P0，S=AsAb=A3s，SAb，sq=（l-w）/w+2Uh/U*-l［2-lu*/（2wUh）］/w（16）

1.4比值Uh/u*的確定

由前文可知，要預(yù)測(cè)和比較粗糙元的陰影保護(hù)區(qū)域投影面積及其保護(hù)效率，除了知道粗糙元本身的幾何特征以外，確定Uh/u*的比值是關(guān)鍵.對(duì)于極低粗糙元密度情況，比如只有1個(gè)粗糙元存在于較大面積的地表上時(shí)，可以認(rèn)為粗糙元對(duì)來(lái)流風(fēng)速的影響可以忽略，此時(shí)地表時(shí)空平均的水平風(fēng)速垂向廓線(xiàn)近似滿(mǎn)足對(duì)數(shù)規(guī)律，即公式（17）.這樣，只要知道h，u*和z0，即可求得Uh/u*.

U（z）=u*kln（zz0）（17）

隨著粗糙密度的增大，風(fēng)速對(duì)數(shù)廓線(xiàn)規(guī)律不再適用，且圓柱體粗糙元和長(zhǎng)方體粗糙元情況下的Uh/u*求解方法有所不同.對(duì)于圓柱體粗糙元，結(jié)合Raupach［1］提出的λd，m和公式（5）即可求得Uh/u*.對(duì)于長(zhǎng)方體粗糙元，通過(guò)Yang等［15］提出的循環(huán)迭代方法可以求得Uh/u*.

1.5其他設(shè)置

本研究涉及的邊界層厚度為4倍粗糙元高度［15］.為了較為詳細(xì)地比較2種形狀粗糙元對(duì)地表保護(hù)作用的差異，本研究主要考慮在地表顆粒臨界啟動(dòng)風(fēng)速下的情況.將覆蓋粗糙元的地表假設(shè)為均勻粒徑組成，其粒徑D大小為0.1～0.6 mm.這些粒徑流體啟動(dòng)的臨界啟動(dòng)風(fēng)速u(mài)*ft由公式（18）給出［25］.其中，重力加速度g=9.81 m/s2，空氣密度ρa(bǔ)=1.225 kg/m3，沙粒密度ρs=2 650 kg/m3，常數(shù)β0取值為0.000 29.本研究涉及到的面率AR值為0.05～20.地表粗糙覆蓋度λc值的范圍為0.000 2～0.111 1，粗糙元密度λd值的范圍為0.002 7～0.666 7.

u*ft=0.111（ρs-ρa(bǔ)）gD/ρa(bǔ)+β0（ρa(bǔ)D）（18）

2結(jié)果與討論

2.1單個(gè)粗糙元情況下的比較

本研究首先進(jìn)行了單個(gè)粗糙元保護(hù)效率的比較.由公式（17）可知，對(duì)于高度相同的單個(gè)粗糙元，計(jì)算保護(hù)區(qū)域面積的關(guān)鍵參數(shù)Uh/u*前需確定地表粗糙度z0.對(duì)于由不同粒徑組成的地表床面，其z0在地球環(huán)境下可近似估計(jì)為D/30.因此，本研究首先討論了單個(gè)粗糙元保護(hù)效率隨著地表粗糙水平的變化，結(jié)果如圖4所示.圖中給出了高和低2種面率比，即h/w=4.0和h/w=0.25時(shí)保護(hù)效率的變化.可以看到，隨著粒徑的增加，2種形狀的粗糙元的保護(hù)效率都在逐漸降低.但是，長(zhǎng)方體粗糙元的保護(hù)效率在2種面率情況下都大約是圓柱體粗糙元的保護(hù)效率的2倍，即使此時(shí)2種粗糙元所形成的粗糙密度λd都是一樣的.同時(shí)，考慮到兩者寬度w是相同的，長(zhǎng)方形粗糙元的底面積比圓柱體要大.因此實(shí)際上估算得到的總的保護(hù)區(qū)域面積（投影底面積和陰影保護(hù)區(qū)域面積之和，見(jiàn)圖1）也是長(zhǎng)方體粗糙元比圓柱體粗糙元的面積大.

隨著床面粒徑的變化隨著面率的變化，長(zhǎng)方體粗糙元的保護(hù)效率依然比圓柱體粗糙元的保護(hù)效率高，如圖5所示.圖5中，所有面率下的圓柱體粗糙元的保護(hù)效率都由公式（3）求得；而對(duì)于長(zhǎng)方體粗糙元，在面率lt;1時(shí)保護(hù)效率由公式（8）求得，在面率≥1時(shí)保護(hù)效率由（13）求得.圖5表明，對(duì)于圓柱體粗糙元，隨著面率的增大，保護(hù)效率都呈先增大后減小并逐漸趨于穩(wěn)定的變化規(guī)律.然而，雖然在低面率時(shí)圖5（A）和圖5（B）都表明長(zhǎng)方體粗糙元保護(hù)效率在逐漸增大；但是，當(dāng)面率≥1時(shí)，圖5（A）和圖5（B）顯示的長(zhǎng)方體粗糙元保護(hù)效率變化規(guī)律有所不同.即圖5（A）顯示保護(hù)效率隨著面率近似線(xiàn)性增加規(guī)律，而圖5（B）的保護(hù)效率不再隨面率變化.由于圖5（A）中w是固定的，而h是變化的，而圖5（B）中h不變而w是變化的.表明長(zhǎng)方形粗糙元的保護(hù)效率對(duì)高度變化比對(duì)寬度顯得更加敏感.

2.2低蓋度情況下的比較

隨著粗糙元密度的增加，流體與粗糙元的相互作用變得顯著，粗糙層內(nèi)及附近高度區(qū)域內(nèi)的風(fēng)速將發(fā)生改變.同時(shí)，由于幾何形狀的差異，圓柱體粗糙元和長(zhǎng)方體粗糙元對(duì)流場(chǎng)的擾動(dòng)有一定差異，因而他們對(duì)風(fēng)速的影響也會(huì)有一定區(qū)別.因此，這里首先比較了2種形狀粗糙元在相同粗糙密度λd時(shí)估算得到的Uh/u*比值的差異，結(jié)果如圖6所示.圖6表明，無(wú)論是圓柱體還是長(zhǎng)方體粗糙元，Uh/u*隨著粗糙密度的增加呈逐漸減少并趨于穩(wěn)定的變化規(guī)律.但2種粗糙元造成的Uh/u*隨寬度和高度的變化在定量上存在一定差異.圓柱體粗糙元的Uh/u*隨粗糙密度增加的減緩速率比長(zhǎng)方體粗糙元小.在粗糙密度小于0.01時(shí)，圓柱體粗糙元情況下的Uh/u*比長(zhǎng)方體粗糙元時(shí)的值??；而粗糙密度大于0.01時(shí)，圓柱體情況下的值比長(zhǎng)方體情況下的值大.

本研究還比較了2種粗糙元情況下流向最大保護(hù)長(zhǎng)度Lx隨粗糙密度的變化如圖7所示.總的來(lái)說(shuō)，相比于Uh/u*，Lx的變化顯得更加復(fù)雜.與Uh/u*類(lèi)似，存在一個(gè)臨界粗糙密度，當(dāng)密度小于0.01時(shí)，圓柱體情況下的Lx比長(zhǎng)方體情況下?。欢芏却笥?.01時(shí)，圓柱體情況下的值比長(zhǎng)方體時(shí)的值大.在圓柱體情況下，Lx隨粗糙元密度的變化（無(wú)論隨高度還是隨寬度）都呈現(xiàn)為先增大后減小然后逐漸趨于穩(wěn)定的規(guī)律.然而，在長(zhǎng)方體情況下，估計(jì)得到的Lx隨粗糙元密度的變化出現(xiàn)顯著的不連續(xù)性，并且Lx隨高度和寬度的變化規(guī)律并不一致.圖7（A）顯示，隨著粗糙密度的增加（源自寬度的增加），Lx先近似線(xiàn)性增加而后降低并保持為恒定值.圖7（B）顯示，隨著粗糙密度的增加（源自高度的增加），Lx先增加后略微減小，然后再突然跳躍增大到某一值最后呈逐漸減小并趨于穩(wěn)定趨勢(shì).

隨著粗糙密度的增加，2種粗糙元情況下保護(hù)效率的變化是一致的如圖8所示.隨著寬度的增加，粗糙元保護(hù)效率呈單調(diào)減少并逐漸趨于穩(wěn)定的變化規(guī)律（見(jiàn)圖8（A））；隨著高度的增加，粗糙元保護(hù)效率呈先增大后減小最后趨于穩(wěn)定的變化規(guī)律（見(jiàn)圖8（B））.

就粗糙元形狀而言，2種形狀的影響是不同的.在多粗糙元情況下，長(zhǎng)方體粗糙元的保護(hù)效率比圓柱體粗糙元的效率高，這與單個(gè)粗糙元情況時(shí)的結(jié)論是一致的.隨著寬度的增加，長(zhǎng)方體粗糙元的保護(hù)效率衰減速率比圓柱體大；隨著高度的增加，長(zhǎng)方體粗糙元的保護(hù)效率在增加和減小階段的變化都比圓柱體劇烈.

2.3討論

本研究基于已有粗糙元陰影保護(hù)區(qū)域模型［1，15-16］，分析和推導(dǎo)了均勻分布粗糙元情況下保護(hù)區(qū)域尺寸以及相對(duì)保護(hù)效率的表征.基于這些理論表征，討論了單個(gè)粗糙元及低覆蓋度下多個(gè)粗糙元情況下圓柱體和長(zhǎng)方體粗糙元對(duì)保護(hù)區(qū)域的影響.

通常，一些植物類(lèi)的地表粗糙元被簡(jiǎn)化為圓柱形［4］.這些粗糙元的主要作用是防止或降低地表侵蝕.在沙漠邊緣區(qū)域，地表植被覆蓋度非常低，這時(shí)就可以將地表植被考慮為單個(gè)圓柱粗糙元來(lái)計(jì)算他們的保護(hù)作用.雖然沙粒在傳輸過(guò)程中會(huì)出現(xiàn)分選現(xiàn)象［26］，從而使得地表粒徑分布產(chǎn)生一些差異.但由圖4可看出，雖然粒徑有一定影響，但其影響是可以忽略的.另外，如果要在地表人工種植或鋪設(shè)單個(gè)類(lèi)圓柱體植被或粗糙元，最佳的面率比（高寬比）在5左右，這樣的效率最好.雖然以上結(jié)論是基于指數(shù)形式的f函數(shù)，但不同的f函數(shù)對(duì)結(jié)果的影響相對(duì)較小，不同的f函數(shù)得到的最佳面率比也在5左右（見(jiàn)圖9）.

此外，相比于圓柱體粗糙元，長(zhǎng)方體粗糙元的保護(hù)效率P0和保護(hù)長(zhǎng)度Lx隨變量的變化呈現(xiàn)出不連續(xù)的問(wèn)題.這些不連續(xù)并不表示物理事實(shí)是這樣的（即物理上本身可能是連續(xù)的），而主要是由在考慮面率對(duì)保護(hù)區(qū)域影響時(shí)的理論簡(jiǎn)化造成.比如，當(dāng)面率大于1時(shí)，將所有面率情況下的保護(hù)區(qū)域lw近似為w［16］，而實(shí)際上lw可能仍然隨著面率在緩慢變化.盡管圖8中保護(hù)效率變化的不連續(xù)性似乎不存在了，但圖5和圖7的結(jié)果暗示著目前對(duì)于多粗糙元情況下的保護(hù)效率在面率等于1附近極有可能存在著定量上的高估.許多研究將城市建筑簡(jiǎn)化為長(zhǎng)方體粗糙元［14-16］，因此，提升不同面率下保護(hù)區(qū)域的定量描述，可以提高對(duì)于城市地表粉塵釋放及這些粉塵在城市大氣邊界層中擴(kuò)散的估計(jì)和預(yù)報(bào).

對(duì)于圓柱體粗糙元，雖然理論上考慮的是剪切應(yīng)力小于臨界值的區(qū)域都是保護(hù)區(qū)域，但Raupach［1］確定保護(hù)區(qū)域面積的判斷依據(jù)是區(qū)域內(nèi)切應(yīng)力均為0.因此實(shí)際上基于Raupach［1］的模型很有可能是低估了圓柱體的保護(hù)效率.植被的多樣性使得一些植被并不能被視為不透風(fēng)的圓柱體，而應(yīng)該考慮其為孔隙結(jié)構(gòu)，其保護(hù)區(qū)域形狀和面積也與固體粗糙元存在一定差異［27］.因而，由固體粗糙元轉(zhuǎn)變到不同孔隙率的粗糙元時(shí)保護(hù)區(qū)域的變化可能是后續(xù)的研究方向.對(duì)于長(zhǎng)方體粗糙元，Yang［15］和Sadique等［16］確定保護(hù)區(qū)域的形狀及面積的標(biāo)準(zhǔn)是來(lái)流風(fēng)速的0.2～0.3倍所包括的范圍.如果來(lái)流風(fēng)速恰如本研究選取的臨界風(fēng)速，那么他們的模型給出的結(jié)果應(yīng)該是低估了長(zhǎng)方體粗糙元的保護(hù)效率.這樣，他們的模型其實(shí)是對(duì)來(lái)流風(fēng)速依賴(lài)的.因此，在選用他們的模型來(lái)計(jì)算保護(hù)區(qū)域時(shí)，應(yīng)該注意地表?xiàng)l件和來(lái)流風(fēng)速情況.

3結(jié)論

本研究根據(jù)已提出的粗糙元保護(hù)區(qū)域形態(tài)及其尺寸模型，推導(dǎo)出了低覆蓋度（3wlt;l）情況下，均勻分布粗糙元的保護(hù)區(qū)域最大流向長(zhǎng)度、保護(hù)區(qū)域面積和保護(hù)效率的公式，然后考察和比較了2種粗糙元情況下這些物理量隨幾何特征參數(shù)（面率）和空間分布特征參數(shù)（粗糙密度）的變化規(guī)律，得到以下結(jié)論：

1）對(duì)于單個(gè)或粗糙元非常稀少的地方，沙質(zhì)地表平均粒徑的增加可降低粗糙元的保護(hù)效率但影響相對(duì)較小.

2）相同面率和粗糙密度情況下，長(zhǎng)方體粗糙元的保護(hù)效率比圓柱體粗糙元的要大.

3）隨著粗糙密度的增加，特征風(fēng)速Uh/u*呈現(xiàn)為單調(diào)減少并趨于穩(wěn)定的變化規(guī)律，而保護(hù)區(qū)域最大流向長(zhǎng)度Lx則先增加后減小最后趨于穩(wěn)定值.

4）在粗糙密度小于0.01時(shí)，基于長(zhǎng)方體粗糙元估計(jì)得到的Uh/u*和Lx比圓柱體粗糙元大，而粗糙密度大于0.01時(shí)，圓柱體的值比長(zhǎng)方體的值大.

5）粗糙元均勻分布情況下，粗糙元保護(hù)效率隨寬度呈單調(diào)減少并逐漸趨于穩(wěn)定的變化規(guī)律，而隨著高度的增加，粗糙元保護(hù)效率呈先增大后減小最后趨于穩(wěn)定的變化規(guī)律.考慮到真實(shí)植被多為孔隙結(jié)構(gòu)，后續(xù)研究應(yīng)該關(guān)注粗糙元保護(hù)區(qū)域形狀和尺寸隨植被或結(jié)構(gòu)孔隙度的變化.

參考文獻(xiàn)：

［1］Raupach M.Drag and drag partition on rough surfaces［J］.Bound Lay Meteorol，1992，60（4）：375-395.

［2］Raupach M R，Gillette D A，Leys J F.The effect of roughness elements on wind erosion threshold［J］.J Geophys Res Atmos，1993，98（D2）：3023-3029.

［3］Musick H，Trujillo S M，Truman C R.Wind-tunnel modelling of the influence of vegetation structure on saltation threshold［J］.Earth Surf Process Landf，1996，21（7）：589-605.

［4］Crawley D M，Nickling W G.Drag partition for regularly-arrayed rough surfaces［J］.Bound Lay Meteorol，2003，107（2）：445-468.

［5］Gillies J A，Nickling W G，King J.Shear stress partitioning in large patches of roughness in the atmospheric inertial sublayer［J］.Bound Lay Meteorol，2007，122（2）：367-396.

［6］Hagishima A，Tanimoto J，Nagayama K，et al.Aerodynamic parameters of regular arrays of rectangular blocks with various geometries［J］.Bound Lay Meteorol，2009，132（2）：315-337.

［7］Walter B，Gromke C，Leonard K C，et al.Spatio-temporal surface shear-stress variability in live plant canopies and cube arrays［J］.Bound Lay Meteorol，2012，143（2）：337-356.

［8］Walter B，Gromke C，Lehning M.Shear-stress partitioning in live plant canopies and modifications to Raupach’s model［J］.Bound Lay Meteorol，2012，144（2）：217-241.

［9］Walter B，Voegeli C，Horender S.Estimating sediment mass fluxes on surfaces sheltered by live vegetation［J］.Bound Lay Meteorol，2017，163（2）：273-286.

［10］Webb N P，Okin G S，Brown S.The effect of roughness elements on wind erosion：The importance of surface shear stress distribution［J］.J Geophys Res Atmos，2014，119（10）：6066-6084.

［11］McKenna N C，Bédard O.A wind tunnel study of flow structure adjustment on deformable sand beds containing a surface-mounted obstacle［J］.J Geophys Res Earth Surf，2015，120（9）：1824-1840.

［12］Kastner K P，Rotach M W.Mean flow and turbulence characteristics in an urban roughness sublayer［J］.Bound Lay Meteorol，2004，111（1）：55-84.

［13］Roth M，Inagaki A，Sugawara H，et al.Small-scale spatial variability of turbulence statistics，（co） spectra and turbulent kinetic energy measured over a regular array of cube roughness［J］.Environ Fluid Mech，2015，15（2）：329-348.

［14］Huang G，Simons S，Vinkovic I，et al.Law-of-the-wall in a boundary-layer over regularly distributed roughness elements［J］.J Turbul，2016，17（5）：518-541.

［15］Yang X I，Sadique J，Mittal R，et al.Exponential roughness layer and analytical model for turbulent boundary layer flow over rectangular-prism roughness elements［J］.J Fluid Mech，2016，789：127-165.

［16］Sadique J，Yang X I，Meneveau C，et al.Aerodynamic properties of rough surfaces with high aspect-ratio roughness elements：effect of aspect ratio and arrangements［J］.Bound Lay Meteorol，2017，163（2）：203-224.

［17］Kang L，Zhang J，Zou X，et al.Experimental investigation of the aerodynamic roughness length for flexible plants［J］.Bound Lay Meteorol，2019，172（3）：397-416.

［18］Gillies J A，Green H，McCarley H G，et al.Using solid element roughness to control sand movement：Keeler Dunes，Keeler，California［J］.Aeolian Res，2015，18：35-46.

［19］Shen L，Han Y，Xu G，et al.A novel canopy drag coefficient model for analyzing urban wind environments based on the large eddy simulation［J］.Energies，2021，14（4）：796-1-796-10.

［20］Yao L，Liu C H，Mo Z，et al.Statistical analysis of the organized turbulence structure in the inertial and roughness sublayers over real urban area by building-resolved large-eddy simulation［J］.Build Environ，2022，207：108464-1-108464-20.

［21］Kang J Y，Yoon S C，Shao Y，et al.Comparison of vertical dust flux by implementing three dust emission schemes in WRF/Chem［J］.J Geophys Res Atmos，2011，116（D9）：D09202-1-D09202-18.

［22］Hamed A M，Peterlein A M，Randle L V.Turbulent boundary layer perturbation by two wall-mounted cylindrical roughness elements arranged in tandem：Effects of spacing and height ratio［J］.Phys" Fluids，2019，31（6）：065110-1-065110-12.

［23］Hamed A M，Peterlein A M.Turbulence structure of boundary layers perturbed by isolated and tandem roughness elements［J］.J Turbul，2020，21（1）：17-33.

［24］Yang X，Ge M.Revisiting Raupach’s Flow-Sheltering Paradigm［J］.Bound Lay Meteorol，2021，179（2）：313-323.

［25］Shao Y，Lu H.A simple expression for wind erosion threshold friction velocity［J］.J Geophys Res Atmos，2000，105（D17）：22437-22443.

［26］Livingstone I.Grain-size variation on a ‘complex’linear dune in the Namib Desert［J］.Geol Soc Spec Publ，1987，35（1）：281-291.

［27］Okin G S.A new model of wind erosion in the presence of vegetation［J］.J Geophys Res Earth Surf，2008，113（F2）：F02S10-1-F02S10-11.

（實(shí)習(xí)編輯：姚運(yùn)秀）

On Variation Character of Protecting Role of Single

Roughness Element in Low Coverage

FU Lintao（School of Architecture and Civil Engineering，Chengdu University，Chengdu 610106，China）Abstract：The circular cylinder and rectangular solid are two commonly used roughness elements in aeolian research.Based on the existing models of wind distribution around the two roughness elements，the maximum protecting length in stream-wise direction，protecting area，and protecting efficiency were theoretically derived for uniformly distributed roughness cases under low coverage conditions.The variations of these three physical quantities with area ratio and roughness density were evaluated and compared between the two roughness elements of different shapes.The research result reveals that the maximum protecting length increases firstly and then decreases gradually with the increase in roughness density.The protecting efficiency decreases gradually with the increase in roughness width，while it increases firstly and then decreases gradually with the increase in roughness height.Under identical area ratio and roughness density，rectangular solid elements have a higher protecting efficiency than circular cylinder elements.These results may provide theoretical guide for the control of erosion on the surface of earth.

Key words：area ratio;roughness density;maximum protecting length;protecting efficiency;protecting area