敦煌地区风蚀起沙量的计算方案

敦煌地区风蚀起沙量的计算方案

1气调与净化方法研究沙颗粒是表面活性剂的主要成分之一。在干旱、半干旱地区,一次风蚀起沙过程中,富含营养和有机物质的较小的土壤微粒悬浮于空中,并在大气湍流的作用下向四周扩散,最后输送到几千公里以外的地方。悬浮于大气中尘粒的辐射强迫作用会影响到区域和全球气候,因此,估算地表土壤的风蚀起沙量以便估计出大气中的沙尘浓度对气候变化和预测而言就显得极为重要。对此,国内外已有很多理论和试验方面的研究。Dong等通过风洞试验研究了风沙流的通量廓线;沈志宝等根据大气中粒子数浓度的观测资料,分别估算了敦煌地区戈壁地表在沙尘暴和扬沙阶段的沙尘排放率;Gillette等根据地表土壤风蚀起沙的物理机制,认为沙尘排放量是风速(或摩擦速度)和地表状况(包括土壤水分和地表植被等)的函数;Lu等在考虑土壤的可侵蚀性和塑性压力的基础上,从理论上计算了沙粒碰撞土壤表面时所产生的凹坑的体积,并建立了一个沙尘排放模式;Shao等通过风洞试验研究,证实了地表沙尘排放的主要因素是沙粒的跳跃碰撞,然后又在充分考虑风蚀起沙的物理机制和风蚀起沙过程中地表土壤粒子尺度分布的基础上,进一步完善了已有的沙尘排放模式。本文利用Gillette等、Lu等和Shao等三种不同的计算方案,根据中日合作“风送沙尘的形成、输送机制及其对气候与环境影响(ADEC)的研究”敦煌站的观测资料,分别计算了2002年4月8日戈壁地表的风蚀起沙量(垂直尘粒通量),并与文献的观测估算结果进行比较,最后简要讨论了三种计算方案的优缺点。2地面风蚀起沙的物理机制地表土壤的风蚀起沙过程是运动的空气流与地表上的粒子在界面上相互作用的一种动力过程。风沙物理学中根据地表土壤微粒的尺度(直径d)和所受的合力大小,将其运动划分为三种形式:尺度较大(d>1000μm)的微粒的蠕移运动、沙粒(60μm<d<1000μm)的跳跃碰撞运动和尘粒(d<60μm)的悬浮运动。沙粒在一定的空气动力作用下,可以被抬升到离地表几十厘米的地方,但当空气动力不足以超过重力时,它们又落回到地表面并与地表发生碰撞,这种碰撞很容易使尘粒克服内部粘性力的束缚,脱离地表进入到大气中。尘粒一旦进入大气,内部粘性力的束缚作用完全消失,由于其重力远小于空气动力,从而导致尘粒悬浮于大气中并在湍流的作用下扩散到很远的地方。由此,对于一定类型的地表土壤,其能否因风蚀而引起尘粒的排放,决定于地表的沙粒能否进行跳跃碰撞运动,从宏观上而言,也就是决定于实际的摩擦速度u*(即空气动力)和沙粒脱离地表的临界摩擦速度u*t(即临界空气动力)。当u*>u*t时,空气动力大于沙粒的重力和内部粘性力的合力,跳跃运动发生,地面会因此而起尘;当u*<u*t时,空气动力小于或等于微粒的重力和内部粘性力的合力,跳跃运动不能发生,地面也不会因此而起尘。通常情况下,尘粒不仅直接存在于地表土壤中,而且还覆盖于较大粒子(沙粒等)的表面,形成一个集合体。研究表明,沙粒的跳跃运动,一方面可以使直接存在于地表土壤中的尘粒因碰撞而悬浮于空中;另一方面,在碰撞的过程中,集合体受撞击而破裂,覆盖在其表面的尘粒也会被剥落而悬浮于空中。因此,沙粒的跳跃运动就成为地面风蚀起沙最关键的运动形式,尘粒的排放量取决于能够进行跳跃运动的沙粒的量。风沙物理学中用顺风向沙粒通量Q表征进行跳跃运动的沙粒的量。Q是沙粒漂移强度q(z)沿高度的积分,量纲为g·m-1·s-1,其物理意义是指单位时间内单位宽度从地表到积分高度处所形成的平面(与水平风向垂直)内的沙粒质量。具体计算时,对于只包含单一粒子尺度的土壤,最常用的是Owen方程:Q={csρu3∗g[1−(u*tu∗)2]u∗>u*t0u∗≤u*t,(1)Q={csρu*3g[1-(u*tu*)2]u*>u*t0u*≤u*t,(1)其中cs是Owen系数,在较典型的情况下(d≈150μm,u*≈0.8m·s-1),cs≈0.8,具体应用时,可通过Owen的经验关系式确定;ρ是空气密度;g是重力加速度,取9.81m·s-2。由于沙尘气溶胶巨大的辐射效应,关于地表土壤尘粒排放量的计算和预报一直是风蚀起沙研究的一个重点,对此,风沙物理学中用垂直尘粒通量F表征。F的量纲为g·m-2·s-1,是指单位时间内某一高度处与地表面平行的单位面积内的尘粒的质量。根据上述地面风蚀起沙的物理机制可知,F的大小决定于顺风向沙粒通量Q,即:F=f(Q).(2)F=f(Q).(2)研究表明,F主要包含三个部分:(1)在微弱的风蚀条件下,没有跳跃运动,地表土壤中的尘粒直接在空气动力的作用下被抬升到大气中的部分Fa;地表土壤中的尘粒由于沙粒的跳跃碰撞而悬浮于空中的部分Fb;(3)覆盖在集合体表面的尘粒由于沙粒的跳跃碰撞而悬浮于空中的部分Fc。在较强的风蚀起沙过程中,Fa相对很小,可以忽略不计,Fb和Fc则量级相当,是地表尘粒排放的主要部分。Gillette等在考虑地表土壤风蚀起沙物理机制的基础上,试验研究了在各种条件下下垫面风蚀起沙的临界摩擦速度,并按照Bagnold指出的指数关系,通过拟合试验资料,给出了相应的起沙模型,即:FG={C2u4∗(1−u∗t/u∗)u∗>u∗t0u∗≤u∗t,(3)FG={C2u*4(1-u*t/u*)u*>u*t0u*≤u*t,(3)其中C2是经验常数,近似取1.4×10-15g·cm-6·s3。Lu等在考虑土壤的可侵蚀性和塑性压力的基础上,从理论上计算了沙粒碰撞土壤表面时所产生的凹坑的体积,并建立了一个沙尘排放模式:FL=cagfρb2p(0.24+cβu∗ρpp−−√)Q,(4)FL=cagfρb2p(0.24+cβu*ρpp)Q,(4)其中ca和cβ是两个经验常数,近似取5.0和1.37;f是地表土壤中尘粒的质量分数,由粒子尺度分布情况决定;ρb和ρp分别是土壤的体积密度和微粒密度,ρp=2650kg·m-3;p是土壤的塑性压力,是表征地面阻力大小的一个物理量,由土壤的致密程度决定,松软的土壤阻力较小,p也较小,致密的土壤阻力较大,p也较大;ρb和p对于不同的土壤类型有不同的取值。Shao在充分考虑地表土壤风蚀起沙的物理机制和风蚀起沙过程中地表土壤粒子尺度分布的基础上,提出了一个更完善的计算方案:FS(di,ds)=cY[(1−γ)+γpm(di)pf(di)]Qgu2∗m(ρbηfiΩ+ηcim),(5)FS(di,ds)=cY[(1-γ)+γpm(di)pf(di)]Qgu*2m(ρbηfiΩ+ηcim),(5)其中cY=17cscY=17cs;di是第i个尘粒的直径;pm(d)是地表土壤遭受最小程度破坏时粒子尺度分布的概率密度函数,pf(d)是地表土壤遭受最大程度破坏时粒子尺度分布的概率密度函数,pm(d)和pf(d)表征的是地表土壤两种理想的粒子尺度分布状态,对不同类型的地表土壤其分布特征也不同,具体计算时可通过3到4种高斯对数分布拟合,其形式在文献中有较详细叙述;γ和(1-γ)分别是pm(d)和pf(d)在风蚀起沙时所占的权重,γ=e-t(u*-u*t)n,t和n是两个经验系数,近似取27.3和3;m是直径为ds的沙粒的质量,具体计算时将沙粒近似作为球形;ηf指地表土壤中能够排放到大气中的尘粒的质量分数,ηc指覆盖于集合体表面的尘粒的质量分数,二者均决定于地表土壤的粒子尺度分布;Ω是沙粒碰撞土壤表面时所产生的凹坑的体积,根据Lu等的计算方法确定。上述(4)和(5)式的计算方案中均包含了顺风向沙粒通量Q,对于只包含单一粒子尺度的土壤,可用Owen方程确定。事实上,自然界几乎所有的土壤都包含着多种尺度的粒子。试验研究表明,Owen方程也可用于计算多粒子尺度分布的土壤,在这种情况下,Q的计算可利用(5)式根据土壤中沙粒的尺度分布进行积分来确定。相应地,FL也需要根据沙粒的尺度分布进行积分,FS则需根据沙粒和尘粒的尺度分布进行二重积分,具体的计算方法详见文献。3对气候条件的支撑西北地区的戈壁和沙漠是我国沙尘天气的主要源区之一。为研究戈壁地表的风蚀起沙情况,下面利用中日合作“风送沙尘的形成、输送机制及其对气候与环境影响(ADEC)的研究”项目的观测资料(对该资料的详细描述参见文献),根据以上关于垂直尘粒通量F的三种计算方案,以敦煌地区2002年4月8日这一沙尘天气为例,计算戈壁地表的风蚀起沙量(垂直尘粒通量)。综合以上三种计算方案可知,垂直尘粒通量的确定需要以下4个参数:土壤类型、风蚀起沙过程中地表土壤的粒子尺度分布p(d)、摩擦速度u*和临界摩擦速度u*t。3.1地表上的土壤指土壤的机械组成分类。根据采样分析得到的土壤的平均粒子尺度分布可知,敦煌戈壁的地表土壤中,砂粒(50μm<d≤2000μm)、粉粒(2μm<d≤50μm)和粘粒(d≤2μm)所占的质量百分含量分别为89.8%,10.1%和0.1%。对应于(4)和(5)式两种方案中所采用的美国农业部(USDA)的划分标准,这样的地表土壤为砂土(sand)。由此,(4)式中土壤的体积密度和塑性压力可以确定:ρb=1000kg·m-3、p=2×106N·m-2。3.2土壤粒子尺度分布随时间的变化这是一种随时间变化的瞬时分布,不同于上述采样分析得到的土壤的平均粒子尺度分布,其原因在于,在风蚀起沙过程中,地表的一部分土壤微粒被抬升脱离地表,导致地表土壤的粒子尺度分布随时间发生变化。对此,Shao采用加权的方法进行确定:p(d)=γpm(d)+(1−γ)pf(d),(6)p(d)=γpm(d)+(1-γ)pf(d),(6)其中γ随u*和u*t也即随时间变化。(4)式中的f,(5)式中的ηf和ηc以及最后对Q,FL和FS进行积分计算时所需的沙粒和尘粒的尺度分布都可由此式得以确定。3.3空气动力学相似理论在风蚀起沙过程中,u*表征的是风对地表土壤微粒抬升的能力。根据敦煌站的观测资料和Monin-Obukhov相似理论,u*可利用两层(2m和1m)的风速和温度资料通过空气动力学方法确定:kzu∗∂u∂z=φm,(7)kzu*∂u∂z=φm,(7)其中k是Karman常数,一般取为0.4。∂u/∂z为风速的垂直梯度,取为差分格式Δu/Δz。φm是风的Monin-Obukhov相似性函数,其具体形式参见文献。3.4土壤水分及植被覆盖对风蚀起沙的影响临界摩擦速度表征的是地表对风蚀起沙的阻碍能力。u*t的大小随粒子尺度的变化而变化,对于一定尺度的沙粒,其脱离地表的临界摩擦速度还决定于土壤水分含量w以及地表植被覆盖度λ等。对此,Shao等曾给出计算方法:u*t=H(w)R(λ)AN(σρgd+ερd)−−−−−−−−−−−−−√,(8)u*t=Η(w)R(λ)AΝ(σρgd+ερd),(8)其中H(w)和R(λ)分别表征土壤水分和植被覆盖对风蚀起沙的阻碍作用,其表达式根据Fecan和Raupach的研究得到,w和λ均由观测得到;AN和ε是经验常数,分别近似取0.0123和3×10-4kg·s-2;σρ是土壤微粒密度(2650kg·m-3)和空气密度(1.23kg·m-3)的比值,ρ是空气密度。4敦煌壁下地表细度和风蚀起沙过程中生长的垂直尺度f图1是利用上述三种方法计算的2002年4月8日敦煌戈壁地表的顺风向沙粒通量Q、垂直尘粒通量F以及摩擦速度u*的时间变化。从图中可以看出:Q和F的变化与u*的变化非常一致,均较好地反映了该日地表土壤风蚀起沙强度的强弱变化,约在15:00(北京时,下同)风蚀起沙达到最强,这与申彦波等依据其它沙尘天气特征参量分析得到的结果一致。三种方案对垂直尘粒通量的计算结果与Q的变化基本吻合,这反映了Q对于F的决定性作用。在垂直尘粒通量的计算中,Gillette等的结果变化幅度最小,Shao的结果变化幅度最大,其最大值比最小值(不包括0)高出7个量级以上,其主要原因在于Gillette等的方案中只考虑了u*和u*t两个因素对F的影响,而Shao的方案中则详细地考虑了地表土壤风蚀起沙的物理机制,除在计算Q时直接考虑了u*和u*t的影响外,还考虑了风蚀起沙过程中地表土壤粒子尺度分布随时间(或随u*和u*t)的变化对F的影响。从(5)式可以看出,FS与u*之间不再只是Lu等或Gillette等的3次方或4次方的关系,而是更高次方的关系,因而使得F的计算结果对u*的变化更敏感。表1给出了利用上述三种方法计算的敦煌戈壁地表2002年4月8日沙尘天气过程中顺风向沙粒通量Q及垂直尘粒通量F的平均值,作为比较,也给出了沈志宝等根据观测得到的2m高度上的沙尘粒子数浓度估算的沙尘排放率(即垂直尘粒通量,该结果有两种,本文只引用其d>5.0μm的尘粒取平均直径d¯=7.5μmd¯=7.5μm时的结果)的平均值。由于沈志宝等只计算了11:31～16:28(沙尘天气较严重)的值,为了便于比较,三种计算结果也只对这段时间求平均。另外,文献的观测估算结果是假定2m高度上的粒子数垂直通量等于地面的起沙量,忽略了尘粒的水平输送,因此与真实的垂直尘粒通量会有一定差异,这里的引用只是作为一种量级上的参考。从表中可以看出,Q的量级为10-4kg·m-1·s-1;F的三个计算结果在数量级上是一致的,为10-8kg·m-2·s-1。三种计算结果中,Lu等的平均值最大,主要原因在于该方案所需的土壤的塑性压力p实质上是随时间变化的,风蚀起沙过程中,随着地表较松软的微粒被抬升到空中,地表土壤的致密性增强,p也随之增大,但由于p的这种变化不易确定,在实际计算中只能取常数,因而导致了FL偏大。另外,参数ca的不易确定性可能也是一个原因。图2是2002年4月8日沙尘天气过程中敦煌戈壁地表的垂直尘粒通量F随u*的变化。从图中可以看出:敦煌戈壁地表土壤风蚀起沙的临界摩擦速度u*t≈0.43m·s-1,这个结果与申彦波等根据沙尘浓度指数的临界值得到的结果比较一致;F的三种计算结果随u*的变化近似以u*=0.7m·s-1为界。当u*t<u*≤0.7m·s-1时,F随u*增大较快,且三者之间有量级上的差异,其中Gillette等的结果最大,但随u*的变化幅度最小,Shao的结果最小,但随u*的变化幅度最大,这与图1中关于FS对u*的敏感性的分析也是一致的。而当u*>0.7m·