搜索
风蚀沙流层质量通量垂直分布的数值模拟及沙粒跃动特征参数的估算1 李万元2，吕世华，沈志宝（中科院寒区旱区环境与工程研究所西部气候环境与灾害重点实验室，甘肃兰州，730000）摘要风蚀地表上跃动沙流层的垂直质量通量结构与沙丘形成直接相关，也与沙流危害程度密不可分。Shao（2005）建立的沙流层相似模型能够很好地模拟出裸露、干燥、疏松、单一沙地地表上沙流层的质量通量垂直分布廓线，Liuetal.（2006）也针对该类地表作过一系列风洞试验，试验结果为相似模型的检验和修正提供了实践依据。本文将基于试验数据和模拟结果的对比分析，结合模型参数的物理意义和前人的相关定量表征方法，对沙粒跃离地表的平均初速度和地表粗糙度分别进行修正，从而改进模型对沙流层质量通量垂直分布廓线的模拟能力；然后利用改进模型对地表沙流通量、沙流层的特征高度和平均高度等重要参量进行模拟分析，以揭示其对摩擦速度、沙粒粒径和地表粗糙度的响应方式及其机制；昀后依据修正后的沙粒跃离地表的平均初速度合理估算出空气在垂直方向上对跃升沙粒的平均作用力，并详细剖析摩擦速度、沙粒粒径和地表粗糙度对初速度和作用力的影响方式，这些工作对研究沙粒跃动的特征轨迹、合理估算沙粒撞击地表的平均速度以及空气对沙粒跃动过程的平均作用皆具理论和实践意义。关键词：跃动沙流层；质量通量垂直结构；相似模型；沙流特征和平均高度；沙粒跃离地表的初速度；空气对跃升沙粒的垂直作用力1引言风蚀地表附近形成的“沙尘云团”中，绝大部分组分粒子（约占75％以上）是跃动着的沙粒[1-3]。跃动沙粒所致的危害并非仅仅通过总沙流通量的作用来实现[1,4-11]，同时也取决于沙流中质量通量的垂直分布状况。其实，后者正是无数单个沙粒跃动轨迹的总体表现，越来越多的科学家已将之视为研究风吹沙流的关键环节[3,12-15]。风蚀地表上沙流层的垂直 1本文的英文版已发表于Boundary-LayerMeteorology(2008)127:313-332。2E-mail:ywl@lzb.ac.cn.648 结构也与沙丘的形成直接相关。在野外观测或风洞试验中，总沙流通量一般不能直接获得，只能从不同高度上测得的沙粒质量和风速廓线出发，先给出表征质量通量垂直分布廓线的经验性函数，然后将该函数对高度积分，积分结果被视为总沙流通量[2,9,16]；可见，在这一实践意义上，沙流质量通量的垂直分布也理应受到足够的重视。多年来，大量的风洞试验和野外观测使人们对风蚀沙流质量通量的垂直分布状况已有了很深的认识。实践证明，风蚀过程中约90%的粒子是集中于近地面30-87厘米的范围内进行输送的，60-80%集中于5-10厘米的范围内，平均传输高度约为63厘米[2-3,14]。科学界已达成广泛共识的是，在地表的跃动沙流中，特别当风速处于中等或较大水平时，沙粒质量通量随高度增加而呈指数递减，且递减率随风速和沙粒的增大而减小，沙流的平均高度却随风速和沙粒的增大而增加[2-3,9,13,16,18]。也有一些科学家发现，沙流通量或沙粒浓度会随高度的增加以幂函数形式递减[2,14,18-19]。Chenetal.（1996）却发现塔克拉玛干沙漠地表上的沙流通量垂直分布廓线既不能用指数函数也不能用幂函数来描述[8,15,18]。Dongetal.[3,14]在研究风蚀地表纵向（顺风向）尺度对沙流通量廓线的影响时发现，在沙粒组成的地表上，沙流通量先是随高度增加而增加，然后很快变为随高度呈指数式递减；在砾石覆盖的地表上，沙流通量廓线却可用一个带有峰值的高斯函数来较恰当地表征。不同的科学家得出如此不同的结论，直接原因就是风蚀事件中影响沙粒起动过程及其运动轨迹的因素相当复杂[3,14]。作者以为，所有结论中昀具代表性的就是Liuetal.所给出的显式指数递减函数，根据这一函数可准确估算出相应的总沙流通量[16]。数值模拟是研究沙流通量垂直分布的另一重要而有效的途径。Owen针对单一的沙地表面，在假定风速遵守对数廓线律，且可完全忽略湍流的情况下，建立了一个风蚀跃动沙流模型，并给出沙粒跃动轨迹和顺风向沙流通量的解析表达式[4,12]。Anderson&Haff（1991）和McEwan&Willetts（1991）在充分考虑沙粒溅射效应的基础上，对Owen（1964）模型进行了有效改进[12]。Shaoetal.用自己发展的数值模型跟踪模拟了跃动沙流中风速和动量通量垂直廓线的发生和演变过程，模拟结果和风洞试验结果非常一致，并以此为基础给出有效粗糙度的估算方法[12]，这一方法与Raupach（1991）的模拟结果和Gilletteetal.（1998）649 的试验结果很接近[20-21]。Shao（2005）以沙粒运动方程为基点，将统计方法和相似理论结合起来，建立了专门描述单一沙地地表上跃动沙流层中质量通量垂直分布廓线的跃动沙流层相似模型，这一模型可利用沙粒粒径、地表粗糙度和摩擦速度定量给出跃动沙流层中沙粒浓度、质量通量和动量传输率随高度变化的垂直分布函数。该模型涵盖了风蚀过程的两种极限情况，一是弱湍流中较大沙粒形成的沙流层，一是强湍流中较小沙粒形成的沙流层；前者的沙流通量廓线是指数形式，后者的则是高斯函数形式，通常风蚀过程中跃动沙流层的通量廓线介于二者之间，并可表达为二者的某种合成[15]。Shao（2005）模型昀初是针对单一沙地而建立的，对于多种粒径沙粒组成的地表上产生的跃动沙流层，其通量廓线可看作各种单一粒径廓线的某种权重式合成[15]，正是基于这一理念，Shao&Mikami（2005）进一步建立了针对复杂地表的跃动沙流层相似模型，并应用该模型合理地模拟出2002年4月沙尘暴期间塔克拉玛干沙漠南缘的地表土壤粒度分布和总沙流通量[22-23]。在沙尘暴频发的荒漠区域，干燥、裸露、疏松的沙地比比皆是，这种地表上形成的跃动沙流层昀易在风洞中进行模拟，事实上，目前已积累了许多针对此类地表的可靠的相关风洞实验和野外观测数据[2-3,9,13,16-18]。另外，由Shao（2005）建立的跃动沙流层相似模型也能够很好地模拟出这类沙地地表上形成的沙流通量垂直分布结构，这在其稍后的工作中已得到证实[22-23]。然而，分析表明，这一模型还存在较大的误差和较多的不确定因素，利用实验数据对之改进是非常必要的。本文将重新并进一步深入地分析研究Liuetal.（2006）针对裸露、干燥、疏松、单一沙地地表上形成的跃动沙流层垂直结构所作的一系列风洞试验结果，然后基于试验数据和相似模型模拟结果的对比分析，同时结合模型参数（如沙粒跃离地表的平均初速度等）的物理意义，给出跃动沙流层相似模型的有效改进方案，并利用改进模型对跃动沙流层结构和沙粒特征跃动轨迹进行模拟研究，从而揭示其对各因子的响应方式及其机制[16]。下面即首先介绍跃动沙流层相似模型及其存在的缺陷，并有针对性地给出改善模型的思路和方法；同时，根据沙流通量垂直廓线模拟结果和实验结果间的对比情况证明这些思路和方法的合理性和有效性。然后，借助改善后的模型，详细探讨跃动沙流层的质量通量650 垂直结构及其随摩擦速度、地表粗糙度和沙粒粒径的变化规律。昀后，据修正后的沙粒跃离地表的平均初速度特别估算出空气在垂直方向上对沙流层中跃升沙粒的平均作用力，并据此判断摩擦速度、地表粗糙度和沙粒粒径对沙粒跃动状态的影响方式。2跃动沙流层相似模型简介在跃动沙流层模型的建立过程中，Shao（2005）首先假设相同粒径的沙粒跃离地表的初速度向量遵循一定的概率密度分布函数（PDF）),(00αvp（0v和0α分别表示沙粒初速度的大小和方向）；然后将表征湍流的拉格朗日统计学模型嵌入沙粒的运动方程，这样就可求出一定粒径的沙粒在一定风力条件下的所有可能跃动轨迹[15,24]。通过分析这些轨迹的综合效应，即可求出沙粒浓度c、沙流质量通量q以及沙粒动量传输pτ随高度z变化的函数关系。根据观测试验所得经验，合理假定各种风力、地表条件下沙粒浓度);(dzc、沙流质量通量);(dzq和沙粒动量传输);(dzpτ的垂直分布廓线均具有相当的动力相似性，并分别给出表征这些量垂直廓线的普适性相似函数);(dzJc、);(dzJq和);(dzJτ。毫无疑问，这些函数要受三个因子的左右：其一是对风速廓线和湍流强度均有重要影响的摩擦速度*u；其二是对风速廓线和湍流动能耗散率皆有影响的地表粗糙度0z；昀后即是沙粒粒径d，它既决定着粒子对湍流的响应时长，也决定着粒子降落的末速度。基于以上假定，通过对不同粒径d、不同摩擦速度*u、以及不同地表粗糙度0z情形下);(dzq的数百次模拟，应用统计和经验拟合技术给出相似性函数);(dzJq的昀终形式，其它两个相似性函数);(dzJc、);(dzJτ也可依照同样的方法求得。下面即是Shao（2005）给出的原始方程：0(;)exp()kqkmzJzdAL=−，（1）);,,(00*00pwdzuAA=，（2）);,,(00*pmmwdzuLL=，（3）);,,(00*pwdzukk=，（4）0(;)()(;)[()]exp()exp()kksqsmkkmmzzqzdFdJzdFdAALL==−=−，（5）651 10);()()(−∞⎥⎦⎤⎢⎣⎡=∫dzdzJdQdFqs，（6）);()();(dzJdFdzspττ=，（7）)()();0()(2*2*0ttppuuddd−=−==ρττττ，（8）102*2*);()()()();0();0(−∞⎥⎦⎤⎢⎣⎡−==∫dzdzJdQuudFddJqtspρττ，（9）这里，0pw是粒子跃离地表的平均垂直初速度，)(dFs是粒子剥蚀率，其它符号和Shao（2005）中的完全相同[15]。正如下面将证实的，其实0pw也决定于*u、d和0z，只是为了强调其重要性，就把0pw显式地置于方程（2）-（4）中。3跃动沙流层相似模型的改进思路和方法0v和0α决定沙粒的跃动轨迹，因而对沙流通量的垂直分布有重要影响，当然也直接影响着以上所引入的普适性相似函数。作为描述0v和0α平均状态的参数，0pw应该随着大气和地表条件的变化而变化，也就是说，0pw应该是*u、d和0z的函数。原始模型假定0pw是摩擦速度*u和地表临界摩擦速度tu*的线性函数，即*063.0uwp=，或tpuuw**0−=，不随0z或其它参数的改变而改变。显然，事实不会如此简单，也就是说，应该对0pw进行更深入的分析。下面这些方程摘自Shao（2005）（参数符号和原著完全相同）[15]：000000(;)(,)ctJzdpvdvdzπδαα∞∆=∆∑∫∫，（10）000000(;)(,)qptJzdupvdvdzπδαα∞∆=∆∑∫∫，（11a）000000(;)(,)pptJzduwpvdvdzπτδαα∞∆=−∆∑∫∫，（12a）方程11a和12a可变形为：652 000000000000(;)(,)(,)(;)qpppcttJzdupvdvdupvdvduJzdzzππδααδαα∞∞∆∆===∆∆∑∑∫∫∫∫，（11b）);();(dzJwudzJcpp−=τ，（12b）这里pu和ppwu−均是针对高度z处所有上升和下降粒子的总体平均。事实上，正如);(dzaτ通常表示空气对顺流方向动量的向下传输一样，);(dzpτ则表示处于高度z的沙粒对顺流方向动量的向下传输。ppwu实际是高度z处的某一粒子对顺流方向动量的垂直传输与其质量的比值，且pw决定传输的速率或快慢，而pu决定待传输的动量大小。当粒子向上运动时，0pw，动量向上传输，然而其对于);(dzpτ的贡献则应为ppppwuwu−=−；当0pw时，粒子对动量的向下传输为ppppwuwu−=，这种传输对);(dzpτ的贡献显然是积极或正面的。可见，不管粒子运动是向上还是向下，它对动量传输);(dzpτ的贡献始终可表示为ppwu−。因此，pu就能够代表高度z