敦煌莫高窟顶灌木林带和草方剂沙障固沙效应研究

敦煌莫高窟顶灌木林带和草方剂沙障固沙效应研究

敦煌是中国乃至世界古代文化艺术的瑰宝。但由于地处我国西北干旱荒漠地带,气候干旱少雨,这一方面为石窟长期完整保存创造了得天独厚的自然环境条件,另一方面,莫高窟窟顶以西约1km处即是高大鸣沙山的东段,沙源丰富,在西北、东北以及西南风的作用下,洞窟长期受到风沙的危害。据统计,上世纪90年代以前,每年要从窟区清除3000～4000m3积沙。1987年,莫高窟被列入世界文化遗产清单,随后,对莫高窟的保护工作得到国际社会的广泛关注,从多学科的角度对其保存状况进行了较为深入系统的研究,特别是在遗存环境方面取得了历史性突破。凌裕泉等、屈建军等和汪万福等在分析莫高窟风沙危害的来源、类型及其风沙流运动特点的基础上,对崖顶的风沙防治问题进行了一系列研究。凌裕泉等对尼龙网栅栏防沙工程的效益问题进行了初步分析。李最雄等开展了莫高窟化学固沙实验,并提出采用高摩数硅酸钾(简称PS)可进行崖面防风蚀加固。薛娴等根据风洞试验和野外观测的结果,研究了砾石覆盖对风蚀的抑制效应。汪万福等对崖顶生物固沙进行了初步研究,并对灌木带的防护效应进行了风洞模拟和野外监测对比实验。朱震达、王涛等和张伟民等对莫高窟风沙危害综合防护体系建立的基本原理进行了探讨。随着敦煌莫高窟风沙危害防治研究工作的不断深入,一个以“固”为主,“固”、“阻”、“输”、“导”相结合,由工程、生物和化学措施组成的多层次、多功能的风沙危害综合防护体系方案已通过专家论证。其中核心内容是在流动沙丘区实施较大规模的草方格沙障固沙和在草方格沙障后方建立一定规模的灌木固沙带,固定流沙,切断沙源,达到根本治理的目的。为了对项目实施后的固沙效应做出科学预测,依据试验阶段已建立的灌木带和草方格沙障的实测参数,结合中尺度大气数值模式MM5和Shao的风蚀起沙方案,计算模拟了2002年4月12日0时到14日0时(UTC)莫高窟顶种植灌木林带和设置草方格沙障前后地表临界摩擦速度、顺风向沙粒通量和垂直尘粒通量等风蚀起沙参量的变化。通过两组结果比较,说明灌木林带和草方格沙障的固沙效应。研究结果为莫高窟风沙防治综合体系的建立提供必要的理论依据,并为防沙工程综合防护效应的评价体系提供借鉴。1研究区植被分布研究区位于莫高窟顶,根据窟顶地貌形态及地表组成物质等因素,从东到西划分为:固定戈壁、砂砾质戈壁、平沙地和沙山4种类型,地表物质粒级含量分别为:砾石依次为40.13%,23.92%,2.70%,0%,细砂为19.23%,19.82%,19.97%,90.60%,微砂为24.60%,29.09%,42.33%,5.27%。研究区植被稀少,戈壁区主要有零星分布的白刺等,盖度不到1%;平沙地和流动沙丘区植被种类相对较多,主要有梭梭、沙拐枣、羽毛三芒草、泡泡刺及一些短命植物等,盖度在1%左右,随着不同季节及年降水量的变化有所波动。在戈壁平坦开阔的地方设有为中日合作“风送沙尘的形成、输送机制及其对气候与环境影响(ADEC)的研究”项目服务的自动气象站(40.04°N,94.80°E),其主要观测项目见表1。自动气象站西南约800m的平沙地上是已种植的灌木固沙试验林带,密度为1株(丛)/4m2,平均高度约1.5m,平均冠幅(迎风面宽度)约1m,平均覆盖度约40%。灌木带以西的中小沙丘上是已设置的草方格沙障,方格规格为1m×1m,露头高度10～15cm,地表麦草覆盖度约40%。实测结果表明,4月份戈壁(气象站附近)、灌木固沙试验林带及草方格沙障(坡底)地表(0～5cm)平均含水量分别为0.063,0.070,0.010m3/m3。2地表土壤水分含量及临界沙粒通量的确定大气数值模式是由PSU/NCAR发展的最新一代非静力中尺度模式MM5v3.6,可对不同区域不同尺度的大气状况进行模拟和预报。模式对于各物理过程提供许多不同的方案,包括降水的物理过程、地表能量的收支、地表边界层的处理和大气辐射的物理过程等,使用者可根据不同要求去选择。其中,增加的陆面过程模式(OSU/EtaLSM)可以为风蚀起沙方案的计算提供地表(0～10cm)土壤水分含量。风蚀起沙方案采用Shao的模型,该模型在充分考虑地表土壤风蚀起沙物理机制的基础上提出,并得到风洞实验和野外观测的验证,是一种计算地面风蚀起沙量的有效方法。在土壤类型(粒度分布)、土壤表面覆盖因子、地表土壤水分含量和摩擦速度等4个基本参量已知的条件下,通过计算可得到地表土壤风蚀起沙的临界摩擦速度、顺风向沙粒通量以及垂直沙尘通量等。本文将MM5和Shao的方案相结合,把实测的模拟区域(照片1)下垫面特征参数(土壤类型、植被覆盖和类型等)在MM5的地形模块中进行替换,从MM5的后处理部分输出地表土壤水分含量和摩擦速度。将获得的参数代入Shao的方案中求解,得到模拟区域内的临界摩擦速度、顺风向沙粒通量和垂直尘粒通量。MM5的大尺度资料采用NCEP提供的再分析资料,取2002年4月12日00:00(UTC)作为初值,采用张弛逼近边界条件,积分时间2d,积分时间步长3s,每1h输出一次结果,区域中心经纬度(94.77°E,40.08°N),东西向和南北向的格点数分别为29和21,水平格距1km,垂直方向分23层,模式大气顶层气压100hpa。3结果与分析3.1风速及温度的计算图1是莫高窟顶ADEC自动气象站附近MM5模拟的摩擦速度u*和利用该站的观测资料(1m和2m的风速及温度)计算得到的u*的比较。从图1可以看出,观测结果和模拟结果在6:00时均达到了峰值,两者的变化趋势基本一致,模拟值与观测的每小时的平均值和最大值呈较好的线性相关,相关系数(R)分别为0.88和0.80,说明MM5对该区域u*的模拟结果是比较可靠的。3.2不同地表土壤水分含量和级联反应表及强度因子地面风蚀起沙的u*t根据Shao等给出的方法计算:u*t=Η(w)R(λ)√AΝ(σρgd+ε/ρd)(1)u*t=H(w)R(λ)AN(σρgd+ε/ρd)−−−−−−−−−−−−−−√(1)式中:d——土壤微粒的直径;g——重力加速度,取9.81m/s2;AN,ε——经验常数,分别近似取0.0123,0.0003kg/s2;σρ——土壤微粒密度(2650kg/m3)和空气密度的比值;ρ——空气密度,对敦煌取1.1kg/m3;H(w)——表征地表土壤水分含量w对风蚀起沙的阻碍作用,其表达式分别根据Fecan的研究得到,对w的取值由于MM5的模拟结果与实际值出入较大,这里的计算以观测值为准,对草方格以西(主要是沙山)近似取0.010m3/m3,对灌木林带以东的戈壁地表近似取0.063m3/m3,对设置草方格前后的地表分别取0.010,0.011m3/m3,对种植林带前后的地表分别取0.010,0.070m3/m3;R(λ)——表征粗糙元密集度λ对风蚀起沙的阻碍作用,其表达式根据Raupach的研究得到,灌木林带和草方格因其高度、冠幅及密度的不同而对风蚀起沙的影响也不同,其R(λ)的大小也有较大差别。根据(1)式可得,灌木林带种植前后该区平均的u*t分别为0.28m/s和1.15m/s,后者是前者的4.1倍。通常情况下,u*不会大于1m/s。因此,灌木林带的种植使地表流动沙粒几乎不会因风吹而脱离地表,从而起到完全固沙的作用;草方格设置前后该区平均的u*t分别为0.28m/s和0.53m/s,后者是前者的1.9倍,说明草方格沙障的设置可明显地减少起沙量,但其阻滞效果不如灌木林带明显。3.3地表细度的计算风沙物理学中,顺风向沙粒通量Q和垂直尘粒通量F用来表征地表的风蚀起沙量。Q反映的是进行跳跃运动的沙粒量,其物理意义是指单位时间内单位宽度从地表到积分高度处所形成的垂直平面(与水平风向垂直)内通过的沙粒质量,量纲为g/(m·s)。但具体计算时常用Owen方程:Q={csρu3*/g[1-(u*t/u*)2]u*>u*t0u*<u*t(2)Q={csρu3∗/g[1−(u*t/u∗)2]0u∗>u*tu∗<u*t(2)式中:cs——Owen系数,在较典型的情况下(d≈150μm,u*≈0.8m/s),cs≈0.8,具体应用时,可通过Owen的经验关系式确定。F表征的是风蚀过程中地表尘粒的排放量,其物理意义是指单位时间内某一高度处与地表面平行的单位面积内通过的尘粒的质量,量纲为g/(m2·s)。本文采用Shao的方案进行计算:F(di‚ds)=cY[(1-γ)+γpm(di)pf(di)]Qgu*2m(ρbηfiΩ+ηcim)(3)式中:cY=17cs;di是第i个尘粒的直径;pm(di)是地表土壤遭受最小程度破坏时粒子尺度分布的概率密度函数;pf(di)是地表土壤遭受最大程度破坏时粒子尺度分布的概率密度函数;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等的计算方法确定。图2是结合MM5和Shao的方案得到的2002年4月13日06时(UTC)模拟区域内F的空间分布情况。从图2中可以看出,种植灌木林带和设置草方格沙障前后地表垂直尘粒通量有明显的变化。分别对灌木林带区和草方格区的Q和F求平均值(表2)后可以看出,设置草方格前后的Q分别为1.25×10-2kg/(m·s)和4.95×10-3kg/(m·s),后者约为前者的2/5,F分别为6.36×10-7kg/(m2·s)和8.54×10-8kg/(m2·s),后者约为前者的1/8,说明草方格的设置使其地表的风蚀起沙量有了明显减少;灌木林带建立前后Q与F的变化更为显著,林带建立前Q与F分别为1.12×10-2kg/(m·s)和5.91×10-7kg/(m2·s),林带建立后则Q与F均减小为0,充分说明灌木林带的建立完全抑制了其地表的风蚀起沙,达到了阻固流沙的良好效果。4地表风蚀起沙的临界速度mm5本文结合中尺度大气模式MM5和Shao的风蚀起沙方案,对敦煌莫高窟顶风沙危害综合防护试验研究中的灌木林带和草方格沙障的固沙效应进行了数值模拟研究,得到以下结论:(1)MM5模拟的摩擦速度和观测的摩擦速度在变化趋势上呈较好的线性相关,相关系数(R)达0.8以上,说明MM5对该区域摩擦速度的模拟结果是可靠的。(2)研究结果表明,灌木林带的建立,使其地表风蚀起沙的临界摩擦速度(u*t)达到1.15m/s,是林带建立前的4.1倍,地表的顺风向沙粒通量