RWEQ中一个容易被忽略的反直觉现象侵蚀因子越大实际风蚀量为何反而下降在使用修正风蚀方程Revised Wind Erosion EquationRWEQ计算防风固沙服务时我们通常认为风力越强、土壤越易蚀、植被覆盖越低实际风蚀量就应越大。按照这种直觉综合侵蚀因子增大后模型结果理应持续升高。然而在固定地块长度为50 m进行计算时我发现了一个反直觉现象RWEQ计算的实际风蚀量并不是一直随综合侵蚀因子增大而是先增加、达到峰值然后下降。具体来说当综合侵蚀因子约为61.73时50 m位置的实际风蚀量达到最大此后综合侵蚀因子继续增大实际风蚀量反而逐渐降低。这究竟是模型错误、程序错误还是RWEQ公式本身的特征一、从RWEQ公式开始RWEQ首先用以下公式计算最大输沙能力Qmax⁡109.8P Q_{\max}109.8PQmax​109.8P其中综合侵蚀因子为PWF×EF×SCF×K′×C PWF\times EF\times SCF\times K\times CPWF×EF×SCF×K′×C这里WFWFWF为气象因子EFEFEF为土壤可蚀性因子SCFSCFSCF为土壤结皮因子K′KK′为地表粗糙度因子CCC为植被或地表覆盖因子。随后计算临界地块长度S150.71P−0.3711 S150.71P^{-0.3711}S150.71P−0.3711距离迎风边界xxx处的局地土壤流失量为SL(x)2xS2Qmax⁡exp⁡[−(xS)2] SL(x)\frac{2x}{S^2}Q_{\max} \exp\left[-\left(\frac{x}{S}\right)^2\right]SL(x)S22x​Qmax​exp[−(Sx​)2]固定x50 mx50\ \mathrm{m}x50m将Qmax⁡Q_{\max}Qmax​和SSS代入可得SLa(50,P)0.48341P1.7422exp⁡(−0.11007P0.7422) SLa(50,P)0.48341P^{1.7422} \exp\left(-0.11007P^{0.7422}\right)SLa(50,P)0.48341P1.7422exp(−0.11007P0.7422)这个函数同时包含两个方向相反的部分P1.7422P^{1.7422}P1.7422促使实际风蚀量增加而指数项exp⁡(−0.11007P0.7422)\exp(-0.11007P^{0.7422})exp(−0.11007P0.7422)促使实际风蚀量衰减。因此最终结果并不是简单的线性上升而是典型的单峰曲线。二、转折点在哪里对SLa(50,P)SLa(50,P)SLa(50,P)求导并令导数等于0可以得到峰值P∗61.73 P^*61.73P∗61.73此时SLamax⁡≈60.86 SLa_{\max}\approx60.86SLamax​≈60.86因此曲线可分为三个阶段阶段综合侵蚀因子实际风蚀量变化上升阶段0P61.730P61.730P61.73随PPP增大而增加峰值P61.73P61.73P61.73达到最大值约60.86下降阶段P61.73P61.73P61.73随PPP增大反而下降部分计算结果如下PPPSLa(50,P)SLa(50,P)SLa(50,P)10.43355.5491014.5382032.3104054.53761.7360.8588058.10110051.33820017.9425000.373从结果可以看到PPP从1增加到61.73时实际风蚀量快速上升但PPP达到100以后实际风蚀量开始明显下降。当P500P500P500时计算值甚至已接近0。三、RWEQ为什么会采用这种函数理解这一问题首先要区分“水平输沙通量”和“局地土壤流失量”。RWEQ的基础函数并不是一个先上升后下降的输沙函数而是Q(x)Qmax⁡{1−exp⁡[−(xS)2]} Q(x)Q_{\max}\left\{1-\exp\left[-\left(\frac{x}{S}\right)^2\right]\right\}Q(x)Qmax​{1−exp[−(Sx​)2]}Q(x)Q(x)Q(x)表示距迎风边界xxx处的水平输沙通量。它始终随距离增加只是增加速度逐渐减慢最终趋近于最大输沙能力Qmax⁡Q_{\max}Qmax​。当相对“干净”的风刚进入可蚀地表时空气中运动的土壤颗粒较少输沙量较低。随着风沿地表向下游运动跃移颗粒不断撞击地表、破坏土壤团聚体并溅起更多颗粒。这种“颗粒越多—撞击越强—释放颗粒越多”的磨蚀放大作用使输沙通量在前期迅速增加。但是风的携沙能力并不是无限的。当输沙量逐渐接近Qmax⁡Q_{\max}Qmax​时风继续携带新增颗粒的能力越来越弱。大量颗粒仍然在运动但每向下游前进1 m由地表新加入输沙流的颗粒越来越少。RWEQ中的局地土壤流失量实际上是水平输沙通量沿下风向距离的变化率SL(x)dQ(x)dx SL(x)\frac{\mathrm dQ(x)}{\mathrm dx}SL(x)dxdQ(x)​所以Q(x)Q(x)Q(x)虽然始终增加并逐渐饱和但它的导数SL(x)SL(x)SL(x)会先增加、达到最大值然后逐渐趋近于0。这正是RWEQ描述的地块长度效应。原作者也指出发生风蚀时最大的土壤损失往往出现在迎风边界下游几十米范围内而不是随着地块长度无限增加(Fryrear, Chen Lester)。四、为什么固定50 m后SLaSLaSLa随PPP也会下降需要特别注意SLSLSL随距离xxx先升后降具有明确的风蚀过程含义但固定x50 mx50\ \mathrm{m}x50m后SLSLSL随PPP先升后降是Qmax⁡(P)Q_{\max}(P)Qmax​(P)、S(P)S(P)S(P)和固定观察位置共同作用的数学结果。随着PPP增大Qmax⁡109.8P Q_{\max}109.8PQmax​109.8P不断增大说明最大输沙能力提高与此同时S150.71P−0.3711 S150.71P^{-0.3711}S150.71P−0.3711不断减小说明输沙流在更短距离内就会接近饱和。换句话说PPP越大局地风蚀峰值越靠近迎风边界。因此当我们始终观察50 m位置时较大的PPP可能意味着风蚀峰值早已出现在50 m以前。此时50 m位置已经处于局地风蚀率的下降段于是就会出现“PPP继续增大但SLa(50)SLa(50)SLa(50)反而减小”的结果。这里下降的是50 m处新增进入输沙流的土壤量而不是整个地块的累计输沙量。不能据此得出“侵蚀动力越强地块总风蚀量越小”的结论。五、为什么这一现象会导致防风固沙量出现负值防风固沙服务通常表示为SRSLp−SLa SRSL_p-SL_aSRSLp​−SLa​Here,SLpSL_pSLp​is the potential wind erosion amount under conditions without vegetation protection, andSLaSL_aSLa​is the actual wind erosion amount under actual vegetation conditions. Since the potential scenario typically setsC1C1C1, it generally holds that:PpPa P_pP_aPp​Pa​如果两种情景都位于SLa(P)SLa(P)SLa(P)的上升段那么PpPaP_pP_aPp​Pa​自然意味着SLpSLaSL_pSL_aSLp​SLa​防风固沙量为正。但是当潜在情景的Pp61.73P_p61.73Pp​61.73时它可能已经位于曲线下降段。此时虽然潜在情景的综合侵蚀因子更大但50 m处的局地风蚀量却可能小于实际情景最终出现SLpSLa SL_pSL_aSLp​SLa​进而得到SR0 SR0SR0这种负值并不一定表示植被促进了风蚀而更可能是潜在与实际两种情景在局地风蚀曲线上的位置不同所产生的数学交叉。六、全国尺度应用需要特别谨慎RWEQ最初针对具有明确迎风边界的农业地块建立。地块长度xxx具有清楚的物理含义即从不可蚀边界或迎风边界到计算位置的沿主风向距离。在全国尺度栅格计算中如果简单地把所有像元统一设置为x50 mx50\ \mathrm{m}x50m实际上隐含了以下假设每个像元都是一个独立、均质的50 m地块每个像元都从一个无输沙输入的迎风边界开始不考虑相邻像元之间的输沙输入和输出不考虑不同地区真实地块长度和主风向差异。这些假设在区域或全国尺度上通常并不完全成立。因此PPP特别高的地区可能因为SSS很小在50 m处已经进入输沙饱和段从而得到异常偏低的实际风蚀量甚至进一步导致防风固沙服务出现负值。更值得注意的是RWEQ的事件尺度验证表明模型预测的Qmax⁡Q_{\max}Qmax​和SLSLSL与观测值存在显著关系但临界地块长度SSS的预测值与观测值没有达到显著相关且模型存在低估Qmax⁡Q_{\max}Qmax​和SLSLSL、高估SSS的倾向。这意味着在超出原始农业地块条件后不能把SSS及其引起的曲线转折视为毫无不确定性的物理事实。七、结语固定x50 mx50\ \mathrm{m}x50m时RWEQ实际风蚀量在P61.73P61.73P61.73达到峰值此后随PPP增大而下降。这不是因为更强的侵蚀动力真的减少了整个地块的风蚀而是因为较大的PPP一方面提高了最大输沙能力另一方面缩短了达到输沙饱和所需的距离使局地风蚀峰值向迎风边界移动。固定观察50 m位置时该位置可能已经越过峰值因而呈现下降。这一发现提醒我们将地块尺度RWEQ直接用于全国栅格计算时不能只关注参数是否齐全还必须理解每个公式的空间过程含义。尤其是地块长度xxx和临界长度SSS它们不只是普通参数而是决定局地风蚀曲线形态、潜在与实际情景关系以及防风固沙量是否出现负值的关键变量。另外这也告诉我们计算RWEQ因子的时候可能还不能直接用每年的WF因子进行计算可能还是得基于每15天的数据进行计算。WF总的来说是个累计因子如果WF用每个月的因素来算的话P值会很大而它能不能用于这个公式还有待探究。参考文献Fryrear, D. W., Saleh, A. (1996). Wind erosion: Field length.Soil Science, 161(6), 398–404. https://doi.org/10.1097/00010694-199606000-00007Fryrear, D. W., Bilbro, J. D., Saleh, A., Schomberg, H. M., Stout, J. E., Zobeck, T. M. (2000). RWEQ: Improved wind erosion technology.Journal of Soil and Water Conservation, 55, 183–189.Fryrear, D. W., Saleh, A., Bilbro, J. D., Schomberg, H. M., Stout, J. E., Zobeck, T. M. (1998).Revised Wind Erosion Equation (RWEQ). USDA-ARS Technical Bulletin No. 1.Zobeck, T. M., Van Pelt, S., Stout, J. E., Popham, T. W. (2001). Validation of the Revised Wind Erosion Equation (RWEQ) for single events and discrete periods. InSoil Erosion Research for the 21st Century, 471–474. https://doi.org/10.13031/2013.4579Merrill, S. D., et al. (2022). Field measurement of wind erosion flux and soil erodibility factors as affected by tillage and seasonal drought.Soil Science Society of America Journal. https://doi.org/10.1002/saj2.20436