风沙跃移是风沙地貌演化、沙尘释放的主要驱动力。沙尘天气的预报依赖于跃移通量模式的精度。由于粗糙床面（覆盖有稀疏植被、砾石等粗糙元）具有复杂的阻力分解机制，预报这类床面的跃移通量和沙尘通量还存在不确定性。例如，当以粉尘释放模型（Dust Production Model， DPM）预报日常沙尘天气时^［1-4］，发现对蒙古高原戈壁及戈壁草原的沙尘天气出现了漏报的情形^［5］，此类漏报可能与跃移起动模式高估了粗糙床面的跃移起动风速有关^［6］。

对粗糙床面跃移通量的计算而言，DPM模型主要考虑了空气动力学粗糙度等因素对湍流切应力分解的影响^［1-4］，并通过空气动力学粗糙度影响起动摩阻风速和跃移通量。近年来研究表明，粗糙床面湍流切应力分解不仅受粗糙元几何特征和空气动力学性质的影响^［7］，还受到粗糙元分布格局^［8-18］及边界层自由风速^［17-18］的影响（即床面气流剪切力与整个粗糙床面总湍流切应力的比值受边界层自由风速变化的影响，表1、2），这意味着粗糙床面的跃移通量不仅与起动摩阻风速和摩阻风速有关，还可能受边界层自由风速的影响，而DPM模型还没有考虑自由风速的影响。为此，本文设计了6组密实粗糙元和3组孔隙粗糙元覆盖的床面，在边界层风洞中同步观测了不同自由风速下的风速廓线和跃移通量，以理解边界层自由风速等因素对跃移通量的影响。通过建立基于观测结果的跃移通量模型，以期为改进沙尘天气的预报提供参考。

风洞试验在中国科学院西北生态环境资源研究院沙漠与沙漠化重点实验室的野外风洞中完成，其几何尺寸和空气动力学性质见文献［19］。

试验沙样来自于沙漠与沙漠化重点实验室的筛分沙样。试验前用筛分法对其进行了测试。测试结果表明其为中细沙（图1），其由2个中值粒径分别为170 μm和320 μm对数正态分布群体构成，且以后者主导。由于是筛分样品，其几乎不存在100 μm以下细颗粒，这对提高实验结果的可重复性有作用。

本研究包括1个裸沙床面和9个粗糙床面用于风洞试验。粗糙元几何参数设计主要考虑了粗糙元形状（用风向比率来度量）及侧影盖度等因素。表1为2类密实粗糙元覆盖的床面，11~14号床面为细高粗糙元，其方向比率AR（单个粗糙元高度与直径比）为12；26和28号床面的粗糙元属于粗矮粗糙元，其AR分别为0.5和1，显著地低于细高粗糙元。

孔隙粗糙元由直径为2 mm、高度为50 mm 铁钉环绕而成（图2、3，表2）。环绕铁钉数量不同而对应着不同的孔隙率。孔隙率P定义如公式（1）：

式中：H、D_k为虚拟的单个密实粗糙元的高度和直径；n_s、h_s、d_s分别是围成孔隙粗糙元个体的数目、高度及直径。具体地，孔隙粗糙元孔隙率分别为0.15、0.55和0.75，对应着36~38号床面（表2）。

裸沙床面和粗糙床面的覆沙厚度均为10 mm，试验开始前用木板刮平。跃移通量采用阶梯式积沙仪测定。积沙仪的开口高度分别为2、4、6、8、10、12、14、16、18、20 cm，开口面积均为2 cm×2 cm。

风洞试验分别测定6个自由风速（10、12、14、16、18、20 m·s^-1）下的跃移通量。利用固定在10个高度上（3、6、10、15、30、60、120、200、350、500 mm）的10个比托管构成的风速廓线仪同步测定输沙过程的风速廓线及对应的摩阻风速。考虑到积沙仪的容量，所用的积沙时间随自由风速增大而调整。一般地，10 m·s^-1对应的采集时间为10 min，12 m·s^-1对应5 min，14 m·s^-1对应4 min，16 m·s^-1对应2 min，18 m·s^-1和20 m·s^-1对应1 min。考虑到输沙过程风洞边界层底部风速偏离对数廓线的问题，这里对10 mm以上的风廓线用对数方程拟合（若拟合的相关系数的平方值≥0.9000，可用来计算跃移边界层的摩阻风速，若拟合的相关系数的平方值低于此数值则考虑15 mm以上高度风速廓线拟合）。

图3展示了不同粗糙床面风洞试验的场景。

光滑床面起动摩阻风速计算采用Marticorena模式^［3-4］。粗糙床面起动摩阻风速的计算分别采用Marticorena模式^［3-4］和Alfaro模式^［2］，其分别对应公式（2a）和（2b）：

式中：R₁和R₂分别为前述模型对应的起动摩阻风速比；U_*ts为光滑床面的起动摩阻风速，其主要取决于床面颗粒的粒径分布特征；U_*t 是粗糙元床面的起动摩阻风速，取决于U_*ts及光滑与粗糙床面的空气动力学粗糙度Z_0s和Z₀。这里的Z_0s是根据沙粒粒径分布群体中最小的中值粒径的1/30计算^［3-4］，而Z₀是根据测定粗糙床面的风廓线并利用对数律拟合而来^［20］。

跃移通量F_hc的计算采用公式（3）^［3-4］：

式中：C为经验常数，这里取2.61；U_* 为摩阻风速；ρ_a为气流密度；g为重力加速度。若采用Marticorena模式，则用R₁代替R；若用Alfaro模式，用R₂代替R。

图4为Marticorena模式预测和风洞试验观测跃移通量的比较。在高风速下该模式明显地低估了裸沙床面的跃移通量，特别是当自由风速16 m·s^-1情况下模拟值约比观测值低50%。

对Z₀1 mm的粗糙床面而言，NO26床面的观测与模拟跃移通量有较好一致性，模拟比观测值约大10%，而NO11床面模拟比观测值大80%，此差异可能反映了粗糙元的方向比率对跃移强度的影响。虽然这里观测与模拟的一致性有差异，但观测与模拟数值仍保留在同一个数量级。

对Z₀1 mm的粗糙床面而言，除个别数据外NO37床面的观测值比模拟值高20%~50%，而NO38床面模拟比观测值大80%。对NO13、NO14、NO28和NO36等床面而言，计算的跃移通量为零而观测数值为10^-3~10¹ g·m^-1·s^-1，此差异与Marticorena模式过高地估计起动摩阻风速有关。

一般地，Marticorena 模式在低空气动力学粗糙度情况下预测的起动摩阻风速偏小，而在较高的空气动力学粗糙度情境下预测的起动摩阻风速偏大，这可能和模式本身只考虑了空气动力学粗糙度而没有考虑粗糙元结构参数如方向比率等对跃移起动的影响有关^［6］。同时，高估起动摩阻风速也是导致模拟与观测的跃移通量明显差异的重要原因^［6］。

图5展示了Alfaro模式的预测结果^［2］。对裸沙床面而言，Alfaro模式和Marticorena模式预测的跃移与观测的通量之间的线性关系相似，即总体上预测的跃移通量明显低于观测值。对Z₀1.5 mm的粗糙床面而言，其中NO26、NO38粗糙床面模拟数值比观测值高40%，而NO11床面模拟值比实测值高80%。对Z₀1.5 mm的粗糙床面而言，它们（NO12， NO13，NO14，NO28， NO36和NO37）的观测值仅是模拟的2%~30%。和Marticorena模式不同，Alfaro模式能敏感地捕捉到较大空气动力学粗糙度床面的风蚀强度。

总体来说，Marticrena模式预测的跃移通量或高于或低于观测值；而Alfaro模式预测数值明显比风洞观测值高，对Z₀1 mm床面而言，该模式预报比风洞观测数据高40%~80%；对Z₀1 mm粗糙床面而言，该模式预报比风洞观测数据高一个数量级（这可看作是异常情况）。相比较而言，Marticrena 模式的可靠性更高一些。Alfaro模式高估跃移通量的主要原因是低估了起动摩阻风速。鉴于此，在计算起动摩阻风速时，建议采用改进的模式^［6］。

为了增强模式的普适性，这里对裸沙床面和粗糙床面的跃移通量数据进行归一化处理。

式中：NFh为无量纲跃移通量；F_hr与F_h分别为相同自由风速下粗糙床面的跃移通量和裸沙床面的跃移通量。在此基础上来寻求粗糙床面的侧影盖度λ和空气动力学粗糙度Z₀与NFh的定量关系。

图6展示了不同自由风速下NFh随λ和Z₀的增加而以负指数函数衰减的趋势，可以分别用公式（5）和公式（6）表示：

式中：m、n、m'、n'为经验系数，且这些经验系数随着自由风速的变化而变化。

从表3可看到m和n分别可以表示为自由风速的二次函数和一次函数。根据相关系数数值可看出，公式（6）的拟合效果好于公式（5）。与侧影盖度相比较，空气动力学粗糙度能更好地反映湍流动量在粗糙元和跃移颗粒之间的分配，因为后者为空气动力学参数，其数值与风廓线特征有关，属于空气动力学参数，而前者仅反映了床面粗糙元几何分布的特征。另外，此类负指数的衰减趋势随着自由风速的增大而有所减缓（图6，表3），其对应的n和n' 随自由风速增大而明显减小。

本文观测的NFh随λ的变化趋势与前人观测数据具有相似的趋势，但在具体的数值上有差异。例如，与前人的风洞试验数据相比较^［20］，当0.04λ0.1时前人测定的NFh为0.25~0.55，其数值明显地高于本文数据。另外，本文观测的数据大于Jornada Experimental Range的野外试验数据^［10］，而与Owens Lake试验数据^［21］接近。NFh随λ或者随Z₀变化的本质是湍流切应力分解机理，其与前人关于湍流切应力分解系数的变化趋势相似^{［9，11，14，17-18］}，也与跃移通量随λ或者随Z₀的变化趋势相似^［12-13］。

本文观测的NFh随自由风速的变化趋势与前人观测的湍流切应力随自由风速的变化趋势相似，即随着自由风速的增加，裸露床面获得的湍流切应力略有增加^［17-18］。这也意味着在模拟跃移通量时不仅要考虑摩阻风速的作用还需要考虑自由风速的影响。造成此趋势的可能原因是风沙边界层属于非常应力层，即粗糙床面跃移边界层的湍流切应力随着高度的增大而增大。为了便于刻画湍流切应力与跃移通量的关系，一般都把风沙边界层如某个临界高度以上范围看作是常应力层，而实际上并非如此。当采用对数廓线方程近似拟合的摩阻风速很难完全反映边界层切应力的状况，因而当同时考虑自由风速的因素时，有可能更接近风沙边界层动量传递的实际情形。这样无论是理解粗糙床面的湍流切应力分解问题^［17-18］，还是理解本文观测的跃移通量和摩阻风速的关系时，都发现在刻画边界层动量输送时，除了考虑摩阻风速外，自由风速也是十分必要的。

综合地看，粗糙元与边界层气流相互作用过程涉及到自由风速、起动摩阻风速和摩阻风速等3个指标。为了综合地反映此过程，这里引入调整摩阻风速U_*m：

公式（7）实际上是公式（3）的变形，旨在利用该参数反映粗糙床面与边界层气流的相互作用效果。进一步，引进基于调整摩阻风速的跃移通量的定义γ：

若γ是不随自由风速变化的常数，则意味着公式（3）不需要同时考虑自由风速对跃移通量的影响，即只需考虑起动摩阻风速和摩阻风速对跃移通量的影响，否则就意味着必须考虑自由风速的作用，而图7所展示的结果意味着γ随着摩阻风速显著地变化。

图7表明，在相同的λ（λ0.15）或Z₀（Z₀3 mm）情况下，γ具有随着自由风速的增大而增大的趋势，这表明跃移通量受到了调整摩阻风速的明显影响。另一方面，γ还随着λ和Z₀增大而以负指数的形式下降，这表明γ具有与跃移通量相似的功能，反映了不同边界层流状况下的跃移强度。同时，对比负指数函数拟合的相关系数可知，γ与Z₀的相关程度更高，原因与NFh与Z₀相关性高的原因相似。

据图7，γ可表示为λ和Z₀ 的负指数函数公式（9）和（10）。相应地其对应的经验系数见表4。

ω和ω'随着自由风速的增大而减小而χ和χ '随自由风速增加则分别呈现先增后减及线性增加的趋势（表4）。上述趋势表明粗糙床面的跃移通量不仅受到起动摩阻风速、摩阻风速的影响，还受到起动摩阻风速、摩阻风速和自由风速等因素协同作用的影响，因而与公式（3）、（5）、（6）相比较，公式（9）和公式（10）更能反映粗糙床面的跃移强度与边界层风动量输送的关系，因而可以认为公式（10）是目前值得推荐的跃移通量模式。另外，由于跃移通量不仅受起动摩阻风速和摩阻风速的影响，还受到调整摩阻风速的影响，则意味着即使采用修正的起动摩阻风速模式，若仍然采用公式（3）来预测跃移通量则必然带来较大的不确定性。因而这里建议把改进的起动摩阻风速模式和公式（10）搭配，可以取得更好的预测结果。

前人的研究表明^［17-18］，粗糙床面的风动量分解不仅受粗糙元密度、单个粗糙元的阻力系数及床面的空气动力学粗糙度等因素的控制，还与粗糙元的异质格局及导致的湍流切应力的空间分布不均有很大的关系，这是在未来的大尺度跃移通量建模中需要进一步考虑的。

风洞试验测量的粗糙床面的跃移通量随自由风速和调整摩阻风速的变化趋势表明，Marticorena模式和Alfaro模式的预测结果与试验观测数据存在明显的差距，此差距不仅与上述模式预报起动摩阻风速具有不确定性，也与它们没有考虑自由风速对跃移通量的影响有关。新模式综合地考虑了自由风速、调整摩阻风速及空气动力学粗糙度等对跃移通量的影响，其计算结果与风洞试验结果更为接近，这为跃移通量预报提供了新的参考。自由风速对跃移通量的影响可能与粗糙床面的风沙边界层属于非常应力层有关，这意味着刻画跃移边界层风动量时，不仅需要考虑摩阻风速的大小，还需要考虑自由风速的大小。