全球泥沙循环是地球系统的基本过程之一，体现了气候变化、造山运动、物理化学侵蚀和人类活动的相互作用^[1]。陆地系统中，河流泥沙对河流形态、水质与生态系统、河道冲淤、水利设施的建设与管理、冲积平原与河口三角洲的发展等方面有很大的影响^[2-3]。在气候变化和人类活动影响下，全球大多数河流输沙量呈现无明显变化或减少的趋势，而寒区河流的输沙量变化整体以增加趋势为主^[4]。例如，喜马拉雅高山区河流输沙量在近60 a增加速率达到1.299%/a，对河道疏浚以及下游水利工程的正常运行等造成了不利影响^[5-6]；在北极地区、欧洲山脉、亚洲高山区和安第斯山脉的研究^[7]同样表明，实测点位中侵蚀输沙通量呈增加趋势的比例高达95%。寒区输沙量变化的特殊性主要在于气温升高导致寒区出现冰川退缩、积雪减少以及多年冻土退化等现象^[8]，使冰冻环境(如冰川区、多年冻土区和冰冻河床)封存的泥沙在冰雪融水携带下逐渐释放^[1]。然而，由于冰川径流拐点的存在以及冻融作用对水沙过程的影响尚不明确，寒区河流输沙量的未来变化存在极大不确定性^[9]。因此，量化分析气候变化背景下冰川消融、融雪以及土壤冻融这3种冰冻圈过程对寒区流域输沙量的影响，对进一步理解气候变化下寒区侵蚀产沙过程具有重要意义^[10-11]。

目前寒区侵蚀产沙过程研究，主要包括模型模拟、统计类方法分析、同位素示踪以及室内外实验观测，流域及区域尺度的研究多采用前2类方法。其中，统计类方法包括多重双累积曲线法^[12]、泥沙恒等式^[13]以及经验回归模型^[14]等。虽然部分研究计算出了冰冻圈要素对输沙量变化的贡献，但这些方法对于过程和机理的认识存在一定的局限性。同时在降雨和冰雪融水的影响下，寒区流域输沙量季节变化与年际变化很大，只用观测数据进行归因得出的结论代表性有限^[15-16]。模型模拟方面的研究多以数据需求量较少的坡面降雨侵蚀模型应用为主，极少数研究采用耦合了流域泥沙输送过程模拟的水沙综合模型。例如，Aygün等^[17]采用Modified Universal Soil Loss Equation(MUSLE)模型对位于加拿大东部的Acadie流域进行模拟得出，随着融雪期的提前，该流域土壤侵蚀最严重的时期从春季提前到了冬季；Wang等^[18]将分布式水文模型Jena Adaptable Modelling System 2000(J2000)与坡面降雨侵蚀模型Revised Universal Soil Loss Equation(RUSLE)耦合，计算得到黑河流域降雨期和融雪期对全年土壤侵蚀率的贡献率分别为90.67%和9.33%；Tian等^[19]采用不包含冰川模块的SWAT模型模拟并预测了拉萨河流域水沙过程，发现在气温升高驱动下蒸散发增加导致径流减少，使20世纪末的河流输沙量显著减少。然而，这些水沙模型仅部分考虑了冰冻圈要素的影响，难以进行具体的贡献划分以解析侵蚀产沙过程的物理机制。因此，在现有水沙模型的基础上完善模型结构，是定量研究气候变化下冰冻圈要素如何影响侵蚀产沙过程的重要方向。

为此，笔者分别梳理冰川消融、融雪和土壤冻融影响下的侵蚀产沙过程及相关模型算法，对比现有流域尺度分布式水沙模型在寒区的适用性，并基于现有模型的不足及考虑单个冰冻圈要素影响下水沙过程的算法提出模型结构改进的建议，以期未来能够采用改进后的模型进行寒区水沙过程的模拟，以降低其不确定性，并更具体地量化各冰冻圈要素对流域输沙量的影响，进一步明确寒区流域侵蚀产沙过程的变化机制及未来变化。

在冰川运动作用下，冰川区积累了大量的岩石碎屑物质。当冰川退缩时，这些碎屑物逐渐暴露并在冰川融水的冲刷下输送至河道^[7]。因此，冰川对侵蚀产沙过程的影响主要包括2方面：一是冰川融化形成的冰川径流是冰川流域径流侵蚀力的重要来源；二是长时间尺度下的冰川基岩侵蚀(主要包括磨蚀、拔蚀)产生了大量碎屑物。这些碎屑物或通过冰下水文过程输送至冰川末端，或由于冰川退缩而出露，使得冰川区成为重要的泥沙源区^[20-21]。此外，冰川表面和内部也会有少量碎屑物随着冰川融水进入河道。

河流泥沙的输送以径流为载体。冰川流域的径流来源包括冰川融水、积雪融水和降雨等，其中冰川径流是十分重要的径流组分。现有水文模型对于冰川径流量的计算较为成熟，主要分为冰川消融量的计算、冰川动态变化考虑和冰川融水的汇流方式3部分。

冰川消融量的计算包括温度指数算法、能量平衡算法以及修正的温度指数算法^[22]。其中，温度指数算法假设冰川消融量与气温线性正相关，因其计算简单得到了广泛应用。模型中的系数以度日因子表示，可通过在冰川表面实际观测消融量和气温获得。有些模型将度日因子设置为固定值，例如，Spatial Processes in Hydrology(SPHY)模型采用固定度日因子的方案^[23]。此外，Luo等^[24]对Soil & Water Assessment Tool(SWAT)模型与动态冰川HRU方案进行松耦合得到可以计算冰川动态变化的SWAT-glacier，其中的度日因子设置为年积日的正弦函数来反映季节性变化。基于物理过程的能量平衡算法是指通过计算冰川表面能量收支情况，考虑冰川表面的辐射平衡、热传导等过程，最后得到冰川消融量。其计算公式较为复杂，对数据要求较高，基于Variable Infiltration Capacity(VIC)模型改进的VIC-CAS模型以及Cold Regions Hydrological Model(CRHM)模型等采用这种方案^[25-26]。修正的温度指数算法介于以上2种方案之间，即在温度指数算法的基础上加入能量平衡算法中的一些因素。比如引入辐射因子、风速因子或水汽压因子等，其中辐射因子是最常考虑的因素，分布式水文模型J2000采用这种方案^[27-28]。总的来说，温度指数算法及其修正版本适合模拟年尺度或月尺度的冰川消融量，能量平衡算法则更适用于解析冰川消融过程，在日尺度上的模拟效果更好^[29]。此外还有研究表明，虽然使用能量平衡算法还需要辐射等多种数据，但即便输入同样的观测数据，能量平衡算法的模拟效果仍然优于温度指数模型，而后者在不同地区需要采用不同的度日因子也表明其局限性，因此一般认为在分布式模型中采用能量平衡算法更为合适^[30]。

近些年为了考虑气候变化下冰川退缩的影响，一些研究在冰川消融量计算的基础上耦合了冰川动态算法以定期更新冰川冰量，从而使冰川消融量的计算结果更准确。冰川动态变化算法包括最常用的体积-面积关系算法，此外还有变平衡线高度法、参数化法和Weertman滑动定律等^[31]。体积-面积法更新冰川的过程通常包括基于现有的观测数据或模型预测确定冰川的初始面积和体积，并建立冰川体积-面积的相关关系，然后根据冰川的动态变化(如冰体的融化或者积累)计算冰川变化量(消融量或者积累量)更新冰川体积，再应用体积-面积关系得到更新后的面积值，模拟过程中不断重复上述步骤，直至模拟结束^{[24, 32]}。变平衡线高度法首先通过遥感影像目视解译获得10 a尺度上冰川物质平衡线的变化，并按比例计算冰川物质积累区面积和冰川体积，该方法适用于较长时间尺度^[33]。参数化法的核心是根据大量资料得出冰川冰水当量的标准化变化量与不同海拔的相关关系，利用该关系将冰川质量变化量分配到不同海拔的冰川上，再计算冰川变化量^[34]。Weertman滑动定律假设冰川滑动速度与冰川底部的剪切力正相关，同时剪切力还与冰川的厚度和坡度正相关，从而计算出冰川滑动速度，得到冰川运动过程中向下游的输送量，再将流域内气温>0℃的区域结合度日因子模型计算消融量^[35]。

通过以上算法计算冰川消融量后，水文模型中主要假设3种分配方式来计算冰川径流。第一种是假定所有冰川消融量直接汇入地表径流，如J2000、VIC-CAS、VIC-glacier^[32]和SWAT-glacier；第二种是假定冰川消融量部分参与土壤入渗过程，例如SPHY模型中一部分进入土壤中的地下水层，其他部分形成冰川径流进入河道；第三种是假定冰川消融量全部参与土壤入渗过程，与其他径流组分一起进行产汇流计算，比如CRHM模型。不同分配方式体现了模型开发人员对冰川融水汇流路径的不同观点，与开发模型所处的地理环境息息相关。

一般认为，地质时间尺度上的冰川基岩侵蚀形成的碎屑物质是冰川流域主要泥沙来源^[20]。在侧重于地质学研究的冰川基岩侵蚀相关论文中，基岩侵蚀速率可以通过分别计算磨蚀和拔蚀，或者直接采用冰川侵蚀定律计算来获得^{[21, 36]}。其中，磨蚀速率主要与冰川移动速度、冰川底部的碎屑物含量以及冰川对冰川床的压力正相关^[37]；拔蚀速率主要与冰川-冰川床的接触面积、冰川床中突出岩石的高度和承受的差异应力正相关，与在冰川床上发生侵蚀的长度负相关^[38]。这类算法对冰川基岩侵蚀的过程考虑得较为详细，物理意义较强，但由于数据难以获取，目前只有极少数研究在小流域进行了试验模拟。冰川基岩侵蚀定律则认为冰川磨蚀和拔蚀过程主要通过冰川底部的滑动实现，即冰川的底部滑动速度直接决定了冰川对地表的塑造能力。因此基于统计数据建立了冰川基岩侵蚀速率e与冰川滑动速度u_s间的指数关系式e=K_g|u_s|^l(式中K_g和l均为系数)^{[21, 36]}。这些算法中侵蚀速率的常用单位是mm/a，因此理论上可以通过冰川基岩侵蚀速率e和冰川面积计算出冰川基岩侵蚀产生的泥沙总量。

还有一些研究认为，冰川基岩侵蚀形成的泥沙先在冰川底部沉积形成冰碛层，再由冰川底部径流携带至冰川末端。Delaney等^[39-40]基于此假设提出了冰下泥沙输送算法，对发源于阿尔卑斯山脉的3条冰川底部的输沙过程进行模拟。冰下泥沙输送模型的计算过程包括3种情况：当输沙能力小于侵蚀量时，泥沙输送受限制，实际输沙量取决于输沙能力；当冰川下的冰碛层高度达最大值或输沙能力为0时，实际输沙量取0；其他情况下，实际输沙量为前2种情况计算值的平滑过渡。该算法有很多前提假设，且计算公式比较复杂，对观测数据要求较高，应用案例较少。

总的来说，由于冰川基岩侵蚀过程发生缓慢，其算法更多用于大陆演化模型和冰川地貌模型中，而冰下泥沙输送模型对观测资料要求很高，目前只有极少数研究在小流域进行了尝试，再加上冰川区泥沙来源除了基岩侵蚀外还受到冰川周围山坡及冰碛物的影响，因此冰川基岩侵蚀输沙量的计算方法在水沙模型中并未得到广泛应用。

在现有研究中，将冰川区作为一个整体建立降雨量、径流量等可以通过观测或水文模型获得的变量与含沙量或输沙量的关系来计算输沙量是更为常用的方案，相关算法见表 1。

表 1(Tab. 1)

表 1 冰川区输沙量计算方案总结 Tab. 1 Summary of methods for calculating sediment flux in glacial areas 方法类型 Method type 具体算法 Detailed algorithm 参考文献 Reference 建立Q-SSC回归模型Developing regression model of Q-SSC 根据观测数据建立含沙量SSC和径流量Q的线性模型。Developing linear model between suspended sediment concentration SSC and runoff Q based on observation data: SSC=aQ+b [41] 根据观测数据建立水沙关系曲线。Developing sediment rating curve based on observation data：SSC=aQb [42] 分时期(包括按主导侵蚀力划分、按涨退水过程划分、按洪水频率划分或结合多种划分方式)建立水沙关系曲线。Developing sediment rating curve by period (includes by dominant erosive force, by rising and falling process of flood, by flood frequency, or by a combination of the above) [43-44] 建立不同侵蚀力和SSC的幂函数模型再求和。Developing a power function model between different erosive forces and SSC and having a sum： SSCt=a1ERt-l1b1+a2IMt-l2b2+a3SMt-l3b3 式中：ER为日平均侵蚀性降雨；IM为冰川消融量；SM为积雪融水量；t为计算时刻；l1、l2和l3为各侵蚀力对应的滞后时间。In the formula：ER is the average daily erosive rainfall. IM is the glacial meltwater. SM is the snowmelt. t is a certain moment of calculation. l1, l2, l3 are the corresponding time lags. [45] 建立不同泥沙来源和SSC的幂函数模型再求和。Developing a power function model between different sediment sources and SSC and having a sum：                 SSCi=Aia1Qia2Ti+b2+Ai(a3Ii+b3)Qib1+a4Qi 这3个分项分别表示热融过程可供给的泥沙、冰雪融水和降雨可供给的泥沙及基流可供给的泥沙。式中：Qi为日流量；Ai为地貎耗竭指数；Ti为第i天与前7天的气温平均值；Ii为第i天相较与前1天的流量增加量。These three equations represent the sediment that can be supplied by thermal thaw processes, ice melt and rainfall, baseflow, respectively. Where Qi is the daily discharge. Ai is a dimensionless geomorphic exhaustion index. Ti is the 8-day average temperature. Ii is the 2-day discharge increase. [14] 建立Q-SF回归模型Developing regression model of Q-SF 建立不同径流组分和SF的幂函数模型再求和。Developing a power function model between different runoff components and SF：SF=a1Q1b1+a2Q2b2+a3Q3b3…… Q1、Q2、Q3……可以是直接径流、壤中流和基流，或降雨径流、融雪径流、冰川径流和基流等组合。Q1, Q2, Q3…… can be a combination of direct runoff, lateral flow and baseflow, or rainfall runoff, snowmelt runoff, glacier runoff, and baseflow, or others [46-47] 采用机器学习算法，以Qt、Qt-1、Rainfallt和Rainfallt-1为输入项，SF为输出项，其中Qt和Qt-1分别为第t和t-1时刻的径流量，Rainfallt和Rainfallt-1分别为第t和t-1时刻的降雨量。Based on the machine learning algorithm, Qt, Qt-1, Rainfallt and Rainfallt-1 are the input, and SF is the output, where Qt and Qt-1 are runoff at time t and t-1 respectively, Rainfallt and Rainfallt-1 are rainfall at time t and t-1 respectively [48] 注：a、b、a1、b1……为待定系数。Notes: a, b, a1, b1…… are coefficients need to be determined.

表 1 冰川区输沙量计算方案总结 Tab. 1 Summary of methods for calculating sediment flux in glacial areas

其中，分侵蚀力拟合、分径流组分拟合和机器学习这3类算法恰好能与水文模型的输出耦合，具有很好的发展潜力。前2类方法中的参数可以采用最小二乘法、梯度下降法、高斯牛顿法、牛顿法、拟牛顿法以及Levenberg-Marquardt算法等进行估计，再采用交叉验证、正则化等方法进行模型选择和优化^{[11, 44, 49]}。

积融雪过程是寒区秋末到翌年春末十分重要的地表过程，对侵蚀产沙有十分重要的影响。与降雨期(非寒区流域侵蚀产沙主要时期)相比，春末融雪期的解冻土层更薄，加上积雪融化增加了土壤含水量，更容易产生地表径流，从而增加坡面侵蚀量^[50]。同时，融雪侵蚀在土壤冻融作用下会进一步加剧，因此在我国东北部、青藏高原地区以及欧洲和北美洲的高寒山区都有研究表明部分年份融雪期的侵蚀速率可以达到甚至超过降雨期^[51-53]。

现有模型中计算融雪量的方案与计算冰川消融量的基本类似，主要包括温度指数算法和能量平衡算法。在温度指数算法中，由于雪的反照率较大，因此融雪度日因子一般小于融冰度日因子^[54-55]。此外，融雪径流侵蚀力的计算方法多样，在各个涉及该过程的模型中也有所不同，总结见表 2^[56]。

表 2(Tab. 2)

表 2 涉及融雪侵蚀的计算方案总结 Tab. 2 Summary of methods involving snowmelt erosion 模型 Model 具体算法 Detailed algorithm 参考文献 Reference USLE1 直接修正降雨侵蚀力因子：R′=R+1.01P12-3。式中：R和R′分别为修正前后的降雨侵蚀力；P12-3为12月—翌年3月的降水量。Correcting rainfall erosivity factors directly: R′=R+1.01P12-3, R and R′ are the erosive force of rainfall before and after correction, respectively, and P12-3 is the amount of precipitation from December to March of the next year [57-58] RUSLE1 直接修正降雨侵蚀力因子，冬季细沟降雨侵蚀力(Req)wr=-1 877.3+7.22P，式中P为降水量；冬季总降雨侵蚀力(即细沟和细沟间的降雨侵蚀力之和) (Req)wt=(Req)wr/k，式中k为根据经验得到的冬季细沟侵蚀与细沟间侵蚀的比例。Correcting rainfall erosivity factors directly, the empirical formula for the rainfall erosive force of rill in winter, (Req)wr, is: (Req)wr=-1 877.3+7.22P, where P is annual precipitation; the total rainfall erosive force in winter, (Reg)wt (i.e., the sum of the rainfall erosive forces of rill and interrill): (Req)wt=(Req)wr/k, where k is the empirical ratio of the rainfall erosive force of rill to interrill [58-59] WEPP2 用能量平衡算法计算融雪量，然后将融雪视作降雨进行侵蚀计算。包括参与入渗及产汇流过程计算融雪径流量，汇入总径流后基于一系列物理过程公式计算细沟侵蚀量和细沟间侵蚀量。Snowmelt is calculated by energy balance algorithm, and then treated as rainfall for erosion calculations: Including participation in infiltration, runoff generation and routing to calculate snowmelt runoff, then followed by a series of physical process equations to calculate rill and interrill erosion [60-61] SPHY 用温度指数算法计算融雪量，并认为融雪量即融雪径流量，直接汇入总径流后用MMF、MUSLE、INCA、SHETRAN、DHSVM或HSPF模块3计算侵蚀量。此外，在计算降雨侵蚀力时SPHY直接采用总降水量，不区分降雨和降雪。Snowmelt is calculated by temperature index algorithm, it is regarded to be snowmelt runoff and included in the total runoff directly, then erosion will be calculated by MMF, MUSLE, INCA, SHETRAN, DHSVM, or HSPF modules. Besides, when calculating the rainfall erosive force, SPHY uses precipitation directly without dividing rainfall and snowfall [23, 62] SWAT 用引入积雪覆盖比例和积雪温度改进的温度指数算法计算融雪量。然后将融雪等同于动能为0的降雨4进行侵蚀计算。包括参与入渗及产汇流过程计算融雪径流量，同时模型假设24h内积雪均匀消融来估算洪峰流量，最后将总径流量和洪峰流量带入MUSLE5或USLE模型计算侵蚀量。Snowmelt is calculated by an improved temperature index algorithm which introduces snow cover percentage and snow temperature, then usually treated as rainfall with zero kinetic energy to calculate erosion: Including participation in infiltration, runoff generation and routing to calculate snowmelt runoff, meanwhile the model assumes that the snow melts uniformly within 24 h to estimate the peak discharge, and finally the total runoff and peak flood flow are taken into the MUSLE or USLE model to calculate erosion [63] EPIC 用温度指数算法计算逐日融雪量，然后将融雪等同于动能为0的降雨进行侵蚀计算。包括参与入渗及产汇流过程来计算融雪径流量和洪峰流量，最后将总径流量、洪峰流量带入Onstad-Foster版本的USLE模型计算侵蚀量。Snowmelt is calculated by temperature index algorithm, then usually treated as rainfall with zero kinetic energy to calculate erosion: Including participation in infiltration, runoff generation and routing to calculate snowmelt runoff and peak discharge, and finally the total runoff and peak flood flow are taken into the USLE model which modified by Onstad-Foster to calculate erosion [64-65] CRHM 用能量平衡算法计算融雪量，然后将融雪等同于动能为0的降雨进行侵蚀计算。包括参与入渗及产汇流过程后得到融雪径流量，汇入总径流后再输入MMF模块计算总侵蚀量。Snowmelt is calculated by energy balance algorithm, then usually treated as rainfall with zero kinetic energy to calculate erosion: Including participation in infiltration, runoff generation and routing to calculate snowmelt runoff, finally the total runoff will be taken into the MMF model to calculate erosion [26] 注：1USLE(Universal Soil Loss Equation)和RUSLE的计算形式均为A=RKLSCP；2WEPP: Water Erosion Prediction Project model；3MMF(Morgan-Morgan-Finney model)、DHSVM(Distributed Hydrology-Soil-Vegetation Model)、HSPF(Hydrological Simulation Program-FORTRAN)、INCA(Integrated Catchments model)和SHETRAN(System Hydrological European (SHE)-Transport)侵蚀模块的计算公式较多，主要原理是分别计算植被影响下的雨滴击溅侵蚀和径流冲刷侵蚀，再与水流挟沙力比较计算侵蚀量，此外HSPF还简单考虑了土壤板结对侵蚀量的减小；4将融雪等同于动能为0的降雨即认为融雪不发生雨滴溅蚀；5MUSLE的计算形式为A=a(Qqp)KLSCP。Notes: 1Both USLE and RUSLE are calculated in the form of A=RKLSCP. 2 WEPP: Water Erosion Prediction Project model. 3 The formulas of MMF, DHSVM, HSPF, INCA, SHETRAN are too numerous to list. The main principle of these erosion modules is calculating the raindrop splash erosion (under the influence of vegetation) and runoff erosion, and then calculate the erosion by comparing with the sediment transport capacity, besides, HSPF also briefly considers the decrease of erosion caused by soil crusting. 4Snowmelt is treated as rainfall with zero kinetic energy means snowmelt cannot cause raindrop splash erosion. 5MUSLE is calculated in the form of A=a(Qqp)KLSCP

表 2 涉及融雪侵蚀的计算方案总结 Tab. 2 Summary of methods involving snowmelt erosion

土壤冻融对寒区侵蚀产沙过程的影响包括3方面：一是在冻融过程中土壤解冻层和冻结层之间会形成不透水层，当表层未冻结土壤厚度较小时，融水或降雨入渗后则更容易形成地表径流或壤中流，因此径流侵蚀力更大；二是从秋季到次年春季频繁的冻融循环会降低土壤团聚体稳定性，从而增加土壤可蚀性和土壤分离能力；三是当未冻结土壤厚度较大时，坡面可供侵蚀的土壤更充足^[66-67]。

土壤冻融对产流过程的影响在已有模型中得到了较好地刻画，从最初的在传统水文模型中引入冻融参数，逐渐发展到构建基于冻融物理过程的分布式水文模型，在一定程度上能够反映土壤冻融作用对产汇流过程的影响。

在传统水文模型中引入土壤冻融相关算法的研究有很多。例如，关志成^[68]假设土壤冻结程度形成是累积负气温的函数改进了新安江模型；李明亮等^[69]引入冻结土壤导水率随气温呈指数变化的函数改进了分布式水文模型Geomorphology-Based Hydrological Model(GBHM)。引入较为复杂的概念模型可以对土壤冻融对产流过程的影响进行更细致地刻画，比如，Wang等^[70]在Precipitation Runoff Modeling System(PRMS)模型基础上加入了温度变源产流模型，间接考虑了气温降低土壤冻结对产流量的影响；Hydrological Predictions for the Environment(HYPE)模型首先计算了土壤温度和土壤湿度，再通过概念性公式计算冻融深度，从而分类计算入渗产流过程^[71]；SWAT模型中通过地表温度区分土壤冻结和未冻结时期，当土壤冻结时SCS径流方程中的滞留系数变小，产流量和径流响应速度随之增大。Qi等^[72]还开发了一个概念性土壤温度模块来改进SWAT中土壤温度和冻融循环的计算，从而进一步改善其对土壤冻融状态的判断。

直接在模型中计算了冻融过程的有VIC、Cryospheric Basin Hydrological Model(CBHM)、CRHM及Water and Energy Budget-based Distributed Hydrological Model(WEB-DHM)模型等。VIC模型通过求解指定土层之间的热通量来计算土壤含冰量，然后确定土壤含水量，从而修正产流量^[73]。CBHM模型根据地表温度将土壤分成未冻结、部分冻结或完全冻结3种情况，然后再根据经验公式求解水热耦合传输方程得到土壤含水量，从而修正土壤冻融时的产流量^[74]。CRHM模型采用XG算法近似求解Stefan方程的近似解来计算冻结锋的位置，并将水力传导度和孔隙度等参数设置为分层变化，温度同一土壤层内表现为线性变化，然后分冻结、未冻结和中间状态3类情况考虑不同土壤冻结程度的入渗过程^[26]。WEB-DHM模型则是通过土壤与大气、冻结层之间的热量交换计算土壤中可用于融化或冻结的水当量，从而计算出活动层厚度，然后修正不同冻结状态下土壤的导热系数和水力传导度，并根据冻融、产流情况更新土壤热量、土壤含水量等状态量，最后得到产流量。其计算过程中较为特殊的是使用焓而不是土壤温度作为预测变量，以避免在土壤相态变化计算过程中产生的不稳定因素^[75-76]。

土壤冻融作用通过影响土壤结构、土壤水分和活动层厚度等使坡面侵蚀产沙过程发生变化^[77]。实验^[78-79]表明，冻融循环次数过多时会破坏土壤团聚体稳定性，增加土壤孔隙度，降低土壤抗剪强度，使土壤颗粒更容易被水流冲刷携带，从而增加土壤可蚀性。冻融循环与土壤可蚀性相关指标的关系得到了室内实验研究的量化。Liu等^[80]发现，土壤分离能力随冻融循环次数增加而增加，在冻融循环次数达10次左右趋于稳定，并基于实验数据给出了土壤分离能力和冻融循环次数的指数回归关系；Xia等^[81]采用不同含水量下的土壤进行实验发现，类似的回归关系在冻融循环次数与土壤孔隙度之间也存在，而孔隙度也是判断土壤分离能力的重要指标之一。

模型中冻融循环对土壤侵蚀的影响常采用分时期设置土壤可蚀性，或引入年积日或冻融循环次数等因素来修正土壤可蚀性相关参数，以考虑冻融循环对土壤可蚀性的影响。如RUSLE可以选择将逐月的土壤可蚀性因子K值在多年平均K值的基础上设为随时间的余弦函数以综合考虑冻融作用等季节性影响因素，不过计算公式仅限于有原始实验数据支撑的美国部分地区^{[58, 82]}。在物理机制较强的WEPP模型中，土壤可蚀性因子的计算形式为K=K₀f，式中：K₀为修正前的土壤可蚀性因子；f为包含冻融循环次数的函数关系式。公式对参数的需求量较大，计算较为复杂^[60]。

现有模型通常默认坡面土壤的厚度是无限的，只考虑了土壤侵蚀量受到水流挟沙力的限制，未考虑泥沙供应的限制，然而在实际情况中流域泥沙供应受限的情况发生普遍。通常情况下，寒区土壤发育不充分，土层较薄，因此冻结期可被径流冲刷的土壤颗粒十分有限，能够发生侵蚀的土壤厚度和流域面积(本文中将这个三维的空间命名为“土壤可侵蚀范围”)明显小于非冻结期^[83-84]。部分模型对于冻融时期土壤厚度对侵蚀的限制通过调整土壤可蚀性表现，而未直接考虑冻融深度影响下土壤可侵蚀范围的变化，使模型在寒区，尤其是多年冻土区的应用存在一定局限性。

已有很多实验研究表明，由于土壤侵蚀量受到地表径流所占比例和土壤可蚀性的影响，未冻结土壤厚度和土壤侵蚀相关变量之间存在复杂的非线性关系。Wang等^[85]基于室内实验发现，土壤侵蚀量随解冻深度先增加后减小；郑颖^[86]在黑河上游冻土小流域进行野外原位冲刷试验，进一步发现浅解冻条件下(解冻深度Z < 30 mm)的侵蚀量明显大于深解冻状态的(Z≥30 mm)侵蚀量。因此，通过在寒区水沙模型中纳入土壤可侵蚀范围的概念来探究土壤冻融深度与土壤侵蚀量的关系是十分必要的。

此外，虽然水沙模型没有考虑坡面土壤可侵蚀范围的影响，但WEPP模型在计算细沟侵蚀和河道侵蚀时对土壤厚度有一定的限制，其细沟侵蚀计算中以耕作层深度为可侵蚀层厚度，河道侵蚀的计算中也在一定深度设置了“不可侵蚀层”^[60]。这表明在模型中通过未冻结土壤厚度来考虑坡面土壤可侵蚀范围，从而限制水流可携带的土壤颗粒是可行的，可以用于完善寒区水沙模型结构以更好的模拟土壤冻融影响下的侵蚀过程。

目前寒区水沙相关模型模拟研究以坡面降雨侵蚀模型如USLE或RUSLE的应用为主，只有极少数研究采用综合了流域产汇流和侵蚀输沙过程计算的分布式水沙模型^{[19, 74, 87]}。因此，笔者调研了多个分布式水沙模型发表论文中的模型结构、计算公式以及公布的模型代码，选择其中考虑冰冻圈要素影响的模型进行对比，总结见表 3。

表 3(Tab. 3)

表 3 寒区分布式水沙模型对比总结1 Tab. 3 Comparative summary of distributed runoff-sediment models in cold regions1 模型 Model 产流部分 Runoff generation 侵蚀产沙部分 Erosion and sediment generation 冰川融水 Glacial meltwater 积雪融水 Snowmelt 土壤冻融 Soil freeze-thaw 冰川区 Glacial area 土壤冻融 Soil freeze-thaw WEPP — 能量平衡算法Energy balance algorithm 通过能量平衡算法计算土壤冻融深度，再分类计算入渗及产流过程Calculating soil freezing or thawing depth by energy balance algorithm, then calculating infiltration and runoff production by different types — 修正了土壤可蚀性因子Modified soil erodibility HYPE — 温度指数算法Temperature index algorithm 根据土壤温度和湿度通过概念性公式计算冻融深度，再分类计算入渗及产流过程Calculating the depth of freezing or thawing by conceptual formulas based on soil temperature and moisture, then calculating infiltration and runoff production by different types — — SPHY 温度指数算法Temperature index algorithm 温度指数算法Temperature index algorithm — — — SWAT2 温度指数算法Temperature index algorithm 温度指数算法Temperature index algorithm 根据地表温度划分土壤冻结和未冻结时期，再分类计算入渗及产流过程Classifying frozen and unfrozen soil periods according to surface temperature, then calculating infiltration and runoff production by different types — — CRHM 能量平衡算法Energy balance algorithm 能量平衡算法Energy balance algorithm 采用XG算法求解Stefan方程计算冻结锋，再分类计算入渗及产流过程Computing freezing fronts by solving Stefan's equation using the XG algorithm, then calculating infiltration and runoff production by different types — — 注：1由于融雪过程对侵蚀产沙过程的影响主要在于促进地表径流，因此其影响归类在产流部分进行对比；2这里的SWAT指的是Luo等[24]提出的耦合了动态冰川HRU方案得到的SWAT-glacier模型。—表明对应的模型未考虑该过程的计算。Notes: 1 Since snowmelt process mainly affects erosion by increasing surface runoff, it is categorized in the runoff production section for comparison. 2 SWAT here refers to the SWAT-glacier model proposed by Luo et al[24]. which couples the dynamic glacier HRU scheme. — indicates that the corresponding model does not calculate this process.

表 3 寒区分布式水沙模型对比总结1 Tab. 3 Comparative summary of distributed runoff-sediment models in cold regions1

从表 3可以发现，目前分布式水沙模型对寒区侵蚀产沙过程的计算都不够完善，尤其是对冰川、土壤冻融影响下的侵蚀产沙过程考虑不足。具体而言，WEPP模型没有考虑冰川融水对产流产沙过程的影响，且模型目前多用于面积较小的坡面侵蚀过程模拟；HYPE模型没有考虑冰川融水对产流产沙过程的影响，且模型更侧重径流和水质模拟；SPHY模型近些年在寒区径流模拟和非寒区侵蚀输沙模拟得到了较为广泛地应用，但没有考虑冰川对侵蚀产沙的影响以及土壤冻融对产流产沙的影响。相对而言，SWAT和CRHM模型在产流部分的模型结构更完善，都在寒区径流模拟研究中得到了广泛应用，但仍然没有充分考虑冰川区和土壤冻融作用在寒区流域侵蚀过程中的特殊性。

综上所述，发展适用于寒区的分布式水沙模型从而量化气候变化下冰川、融雪和土壤冻融对径流量及输沙量的影响，是研究气候变化下寒区侵蚀产沙过程变化的重要方向。已有很多研究分别从冰川、融雪以及土壤冻融这3个角度探讨了其对侵蚀产沙过程的影响及对应的算法。笔者对此进行了总结，并比较了现有分布式水沙模型在寒区的适用性。具体结论如下：

1) 冰川对侵蚀产沙过程的影响主要包括通过冰川融水改变径流侵蚀力和水流挟沙力，以及通过冰川基岩侵蚀作用提供冰川区侵蚀输沙的泥沙物源。其中冰川融水量的计算包括能量平衡算法、温度指数算法及其改进版本，部分水文模型还会辅以冰川动态算法及对冰川融水量采用不同分配方式得到冰川径流；基岩侵蚀的影响则可以通过经验方法将冰川区作为整体得到总输沙量后再与水文模型耦合。

2) 融雪对侵蚀产沙过程的影响主要是通过促进地表产流增加融雪期的径流侵蚀力。现有模型同样采用能量平衡算法、温度指数算法或其改进版本计算融雪量，然后被视作动能为0的降雨参与下一步的侵蚀计算。

3) 土壤冻融作用对侵蚀产沙过程的影响包括冻融深度影响产流量，频繁的冻融循环增加土壤可蚀性以及冻融深度限制坡面土壤可侵蚀范围，通常采用概念性或经验性的方法对这些过程进行计算。

4) 现有分布式水沙模型对寒区侵蚀产沙过程的考虑各有不足，其中较为完善的SWAT和CRHM模型考虑了冰川消融、融雪径流以及土壤冻融影响下的径流过程，但未全面考虑冰川消融和土壤冻融对侵蚀产沙过程的影响。

因此，为了提高现有流域水沙模型在寒区的适用性，结合前文总结的冰川消融、融雪和土壤冻融影响下的侵蚀产沙相关算法，同时考虑到模型改进的实用性与可行性，建议未来可以对SWAT或CRHM的模型结构做出以下改进：

1) 在水沙模型对冰川区径流模拟的基础上，利用分侵蚀力、分径流组分或者机器学习这3种方案建立相关变量与输沙量的关系，结合冰川区末端实测输沙数据进行单独校验后，以此为上边界条件带入模型中进一步模拟下游区域的输沙量，从而反映冰川区输沙模数较大且对气候变化敏感的特点。

2) 在水沙模型对土壤冻融过程计算的基础上，采用表层土壤的冻融循环次数修正模型侵蚀模块中的土壤可蚀性因子，同时基于冻融深度引入“土壤可侵蚀范围”的概念对发生侵蚀的流域面积和土壤厚度加以限制，从而捕捉叠加土壤冻融的融雪期侵蚀输沙量较大的特点。具体而言，对于单位面积土壤，未冻结土壤厚度决定了可被侵蚀的土壤有多少，未冻结土壤厚度>0的流域范围才能发生侵蚀。因此，冻融深度决定了土壤可侵蚀范围，土壤可侵蚀范围再与水流挟沙力一起限制土壤侵蚀量。

未来研究应当采用改进后的水沙模型进行寒区水沙过程的模拟，以减少模型结构带来的不确定性，同时能够更具体地量化冰川消融、融雪及土壤冻融这3类冰冻圈过程对侵蚀产沙过程的影响，从而明确寒区侵蚀产沙过程的物理机制以及流域输沙量的未来变化趋势，为寒区流域的生态保护和可持续发展提供科学依据。