岸滩侧蚀崩塌是冲积河流河床横向变形的重要表现形式，同时也是一种危害较大的自然灾害^[1-2]。对于我国北方地区的季节性冰冻河流而言，河岸侵蚀除了受常规水沙条件影响外，还受冻融作用影响, 水动力、冻融耦合作用特征明显^[3]。以黑龙江流域为例，冻融作用使得土体结构发生变化，抗剪强度降低；在水流淘刷下岸滩土体更易于被冲刷，河岸也更容易变陡，甚至发生失稳崩塌；据不完全统计，黑龙江岸滩侵蚀造成我国农田土地面积平均每年坍失约7 km^2[4]，影响河势稳定和国家土地安全。因此，开展季节性冰冻河流岸滩崩塌研究具有重要的理论意义和实际应用价值。

岸滩崩塌形成机理和发生过程复杂，影响因素众多。针对此问题，国内外很多专家学者开展了崩岸机理研究工作，揭示了影响河岸崩塌的不同因素(包括水流冲刷、不同土体组成、潜水位影响及植被条件等)；与此同时，很多学者从土力学边坡稳定性理论出发，构建了不同崩塌模式下的河岸稳定性计算模型。例如：在揭示河岸崩塌的影响因素方面，Papanicolaou等^[5]以及舒安平等^[6]分析了弯道环流对岸滩侵蚀的影响及其变化过程；文献[7-8]在总结影响长江中下游河道崩岸的基础上，发现长江中下游河段80 % 的崩岸均发生于受水流剧烈冲刷的凹岸弯顶位置或迎流顶冲点，指出近岸水流冲刷是河岸发生崩塌的最主要因素。同时，不同土体组成、潜水位及河道水位等同样对河岸崩退过程产生重要影响。余明辉等^[9]开展了不同近岸河床组成情况下岸坡崩塌试验，并分析了岸坡崩塌与河床冲淤的交互影响过程。假冬冬等^[10]通过引入侧向凝聚力，分析了相邻土体对黏性岸滩稳定性的影响。现有的河岸稳定性分析模型通常是基于土力学中的边坡稳定理论，如：Osman和Thorne^[11]基于平面滑动提出了黏性岸滩失稳崩塌的计算模式。Darby和Thorne^[12]考虑孔隙水和静水压力的影响，对黏性岸滩崩塌模式进行了补充和完善。王延贵等^[13]通过河岸崩塌体的稳定分析，导出了折线岸滩初次崩塌及二次崩塌的临界崩塌高度计算公式。在具体季节性冰冻河流应用方面，李宗华^[14]通过松花江干流治理工程的案例，分析了松花江干流治理时的生态护坡实施过程，并通过对比几种不同形式的护坡，认为生态护坡具有较好的固岸作用。曲洋等^[15]根据丰满水库坝下至河口段河流地形特征选取典型断面，通过河床稳定性指标计算，结合平面套绘结果，分析了松花江该河段河流演变特征及演变趋势，但此类研究并未深入考虑季节性冰冻河流冻融作用的影响，尤其是冻融作用对岸坡崩退的定量贡献，还有待深入研究。

本文以松花江大顶子山航电枢纽下游近坝段典型断面为研究对象，同时考虑土体冻融影响、侧向水压力以及孔隙水压力的作用，结合坡脚冲刷计算模块，采用BSTEM断面尺度模型计算其崩退过程，并定量分析冻融作用对崩岸过程的影响，研究结果可为中国北方季节性冰冻河流岸滩崩塌及河道演变研究提供参考。

1 研究区域概况

松花江是我国七大江河之一，也是黑龙江在我国境内的最大支流。松花江干流属于典型的平原冲积河流，全长约940 km，区间集水面积18.64万km²，河道平均坡降0.08 % 左右。根据地形、地势及河性，松花江干流可分为上、中、下游3段，其中上游段为三岔河哈尔滨，中游段为哈尔滨佳木斯，下游段为佳木斯同江。因流域内植被腐殖质含量丰富，河道两岸表层为呈黑色的腐殖土，表层多为黏土或沙质黏土，下层多为粉沙或细沙，质地疏松抗冲能力较差。松花江河流横向变形主要是通过洲滩或岸滩的侧蚀冲刷来体现。通过现场查勘发现，松花江干流沿程均存在不同程度的洲滩或岸滩侧蚀现象，主槽两侧土质疏松，抗冲性较差，尤其是弯道顶点附近较为明显，多处汊道弯顶出现明显的侧蚀现象，而中游大顶子山枢纽下游近坝段尤为明显，年际最大冲蚀幅度在50 m左右。

对于松花江而言，由于冬季低气温处于主导地位，河道内虽有较大流量、流速，但每年均会封冻，其封冻期一般长4~5个月，因此这类河流也称为季节性冰冻河流，简称“季冻区河流”。与常规的冲积河流相比，季冻区河流岸滩崩塌及其演变过程具有显著差异，水动力、冻融耦合作用特征明显^[3]。

本次所选河段计算断面位于大顶子山航电枢纽下游近坝段(见图 1中的DM1~DM3)。2008年大顶子山航电枢纽工程运行后，上游来沙量减小，使得枢纽下游近坝段出现较为明显的冲刷，局部河段岸滩崩塌现象较为突出，深槽摆动，局部河段洲滩崩退现象时有发生。依据实测断面资料结果，选取该河段3个典型的冲刷崩塌断面作为此次计算断面，分析冻融作用对河岸崩退的影响，计算水文年为2009年。

2 研究方法

美国国家泥沙实验室开发的BSTEM (Bank Stability and Toe Erosion Model)模型，是在Osman等^[11]提出的河岸稳定性分析方法基础上发展起来的，可同时考虑侧向水压力、孔隙水压力及河岸土体分层等影响，结合坡脚冲刷及河岸稳定性模型模拟河岸崩退过程，是运用较广泛的模型之一^[16-17]。

该模型主要由坡脚冲刷模块(TEM)和河岸稳定性模块(BSM)组成。通过录入典型河岸断面轮廓地形、河道水位、比降、糙率、分层土体物理力学参数等基础资料，运行坡脚冲刷模块并根据冲刷后的新地形导入到河岸稳定性模块中计算河岸稳定安全系数F_s。

2.1 坡脚横向冲刷宽度计算

坡脚横向冲刷基于剩余切应力法计算^[18]，用来预测坡脚及岸坡表面的冲刷和下切。作用在河岸土体上的横向冲刷宽度由河岸横向冲刷速率和冲刷时间决定。其中，河岸横向冲刷速率主要由水流冲刷强度及土体抗冲能力决定，只有当作用于河岸土体的近岸水流切应力大于河岸土体的起动切应力时，河岸土体才能起动。河岸横向冲刷宽度E表示为：

$ E = k\cdot({\tau _0} - {\tau _c})\cdot\Delta t $

(1)

式中，E为河岸横向冲刷宽度(m)；k为冲刷系数(m³/(N ·s))，与土体自身特性有关；τ₀为水流切应力(N/m²)；τ_c为河岸土体的起动切应力(N/m²)，一般通过土体冲刷试验获得或根据粒径大小按照经验公式计算得到；Δt为时间(s)。

对于水流切应力的计算，采用水力半径分割法^[19]计算作用在河岸土体上每个节点的平均边界剪切应力：

$ {\tau _0} = {\gamma _\omega }RS $

(2)

式中，γ_ω为水体容重(9.81 kN/m³)；R为水力半径(m)，表示子区域面积与湿周之比；S为水力坡度。

河岸是否发生冲刷由水流剪切力与河岸土体起动切应力的相对大小决定，土体起动切应力被广泛用于代表河岸土体抗冲刷能力，Hanson等^[20]通过大量实测数据得到一个经验公式用于估算非黏性土体临界起动切应力：

$ {\tau _c} = 0.044 \times 16.2 \times {d_{50}} $

(3)

式中，d₅₀为泥沙颗粒中值粒径，大顶子山航电枢纽下游近坝段断面处河岸细沙中值粒径为0.09 mm，通过式(3)计算得到下部细沙临界起动切应力为0.06 N/m²。冲刷系数k与土体自身特性及起动切应力有关，此处非黏性土冲刷系数采用经验关系式^[21]$k=1 \times 10^{-7} \tau_{\mathrm{c}}^{-0.5}$。

黏性土体的起动切应力，根据宗全利等^[22]的分析结果，可用如下经验关系式计算：

$ \tau_{\mathrm{c}}=0.265 \times \rho_{\mathrm{d}}^{3.51} $

(4)

式中，ρ_d为土体干密度(t/m³)，此外，此处黏性土冲刷系数采用水冲刷试验计算公式^[22]计算。$k = 7.677 \times 10^{-6} \tau_{\mathrm{c}}^{-1.949}$计算。

2.2 河岸安全系数计算

河岸稳定安全系数F_s主要采用基于极限平衡法的水平层法计算。首先通过运行边坡形态参数生成河岸剖面，然后结合河岸土体特征参数、河道水位以及孔隙水压力等计算得到岸坡安全系数F_s，如图 2所示。

水平层法由Simon等^[23]研发的楔式崩塌模型发展而来，进一步完善了Osman-Thorne^[11]模型。水平层法可将河岸土体最多分为5层，每层均可定义土体的力学性质指标，安全系数F_s的计算公式为：

$ \begin{array}{c} {F_{\rm{s}}} = \sum\limits_{i = 1}^I {\left( {c_i^\prime {L_i} + \left( {{\mu _{{\rm{a}}i}} - {\mu _{{\rm{wi}}}}} \right){L_i}\tan \varphi _i^{\rm{b}} + \left[ {{W_i}\cos \beta - {\mu _{{\rm{a}}i}}{L_i} + {P_i}\cos (\alpha - \beta )} \right]\tan \varphi _i^\prime } \right)} /\\ \sum\limits_{i = 1}^I {\left( {{W_i}\sin \beta - {P_i}\sin [\alpha - \beta ]} \right)} \end{array} $

(5)

式中，L_i为第i层土体中崩塌破坏面长度(m)；μ_ai为第i层土体孔隙气压力(kN/m²)；μ_wi为第i层土体孔隙水压力(kN/m²)；P_i为由外界水流施加给第i层土体的静水压力(kN/m²)；W_i为第i层土体重量(kN)；c′ _i为第i层土体的有效黏聚力(kN/m²)；φ′ _i为第i层土体有效内摩擦角(°)；φ_i^b为第i层土体表观凝聚力随基质吸力增加而增加的快慢程度(°)；α为河岸坡度(°)；β为崩塌面角度(°)；I为河岸崩塌体总层数。

一般情况下，对于F_s值，若F_s＞1.0，说明河岸稳定；若F_s＜1.0，说明河岸不稳定。本次模拟中，河岸稳定安全系数F_s临界值取值为1.0，作为崩塌与否的判断条件。

2.3 计算参数的选取 2.3.1 土体力学参数

为了获得模型所需要的河岸土体组成及力学性质计算参数，在收集典型河岸断面实测岸坡形态、相应水沙条件等资料基础上，对崩岸典型断面进行了实地查勘，并根据断面土体组成、结构及性质不同，进行了分层取样，由室内土工试验得到了每层土体的物理及力学性质指标。河道两岸表层为呈黑色的腐殖土，表层多为黏土或沙质黏土，下层为粉沙或细沙，质地疏松抗冲能力较差。针对岸滩泥沙样品，通过现场环刀取样测定原状土体的密度；依据《土工试验方法标准》(GB/T 50123-2019)对土样进行剪切试验，得出相应土体的剪切强度(包括内摩擦角、黏聚力)，岸滩密度一般为1.5~1.7 g/cm³，内摩擦角一般为20°~28°，黏聚力一般为7~18 kN/m²。

DM1、DM2、DM3断面位于大顶子山航电枢纽下游近坝段，沿程土体组成基本相同，以DM1断面为例，河岸土体的物理力学特性参数见表 1。

对于季节性冰冻区河流而言，河岸因冻融作用使得土体结构发生变化，土体黏聚力会逐渐下降。土体因冻融作用的影响不仅改变了土体的抗剪强度，而且使得岸滩在水流淘刷作用下更容易发生失稳坍塌，降低岸滩稳定性。黏聚力和内摩擦角是决定河岸土体抗剪强度的重要指标，根据文献[24]的研究，对岸滩崩塌过程进行模拟，即单次冻融后4月份黏聚力减幅为30 % (取均值)，内摩擦角保持不变。同时根据文献[25]的研究，考虑4月份冻融作用对土体结构的影响，在坡脚冲刷模块中，土体干密度降低10 %，以适应天然状态下土体因冻融引起的衰减值。

2.3.2 水位变化情况

松花江流域多年平均降水量一般在500 mm左右，汛期(6-9月)的降水量占全年的60 % ~80 %，冬季(12-2月)的降水量仅为全年的5 % 左右。根据水位变化，参考以往研究划分方法，将2009水文年划分为4个不同时期，分别为枯水期(12月16日-3月15日)、涨水期(3月16日-5月31日)、洪水期(6月1日-9月30日)和退水期(10月1日-12月15日)。同时，根据不同时期崩岸发生频率、水位涨落情况，再将每个时期划分为若干时段计算。模拟时段内水位变化过程如图 3所示。本次模拟中考虑了潜水位的影响，松花江典型河段上层黏性土较薄，下层沙性土较厚，且由于沙性土渗透系数大，潜水位滞后河道水位变化时间很短，故假定河道水位与潜水位同步变化。

3 结果与讨论 3.1 不同时期河岸稳定性

为分析河岸不同时期的稳定性，考虑冻融作用的影响，选取松花江大顶子山航电枢纽下游近坝段3个典型断面进行河岸稳定性计算。根据初始的实测河岸形态、断面处水动力条件及土体力学性质，采用BSTEM模型计算各时期的河岸稳定性，并得到相应安全系数F_s的变化过程(经稳定性计算，DM1断面2009水文年全年稳定无崩塌，故此处未给出)，如图 4所示，不同时期的计算结果分述如下：

1) 枯水期(12月16日-3月15日)，大顶子山水库下泄流量相对较小且变化不大，通过该断面的平均流量较小，河道内水位较低，近岸水流切应力较小，但仍大于下部沙土层的起动切应力，此时河岸下部沙土层存在一定的冲刷。经计算统计，3个断面枯水期坡脚冲刷量分别为22.04、7.42和90.37 m³/m；但由于该时期河岸坡度较缓、孔隙水压力较小，因此河岸稳定程度较高，稳定安全系数F_s值相应较大，断面的平均F_s值基本都在1.4以上，发生崩岸的概率较小。

2) 涨水期(3月16日-5月31日)，冰冻融化，流量相应增加，河道内水位上涨，河岸下部沙土层进一步冲刷。经统计，3个断面涨水期冲刷量分别为17.77、6.77和48.13 m³/m，此时稳定安全系数F_s值在4月下降较快，这是因为季冻区河流在汛前涨水期土体内潜水位升高且存在冻融作用影响土体的抗剪强度，导致河岸稳定性降低。经计算，此时DM2、DM3断面的F_s值分别下降到1.29和1.10左右。实际上，每年4月初左右冰冻融化期，会在短时间内产生独特的春汛现象且随着开河的进行，产生的冰凌亦会对岸滩稳定性造成不利影响，本次模拟尚未考虑冰凌对河岸冲击的影响，因此在实际情况中，汛前涨水阶段仍有崩岸发生的可能性。

3) 洪水期(6月1日-9月30日)，河道内流量增大，水位达到全年最高，河岸下部沙土层大量冲刷、河岸坡度变陡，且洪峰期河道内水位时涨时落，潜水位稳定上升，使得孔隙水压力较大，河岸总体稳定性较低。经计算统计，3个断面洪水期冲刷量分别为49.89、15.37和134.77 m³/m，其中DM2、DM3断面分别发生了1和6次崩塌，崩塌量分别为26.06和251.68 m³/m。在每次崩塌后，由于岸坡坡度变缓以及崩塌土体坡脚堆积的影响，河岸稳定性都会相应增大，并随着坡脚冲刷和水位持续下降，F_s值减小直至再次崩塌，整个洪水期河岸基本表现为此周期性的崩塌。总体来说，洪水期坡脚冲刷最剧烈，崩塌量最多，河岸稳定性较差，为崩岸多发时期。

4) 退水期(10月1日-12月15日)，河道内水位迅速回落，随着水位的下降，作用于河岸侧向水压力逐渐减小，同时坡脚处的冲刷累积作用，使得河岸坡度和坡脚角度增加，安全系数下降，增加了河岸崩塌的可能性。计算过程中，DM2和DM3断面退水期冲刷量为8.85和64.86 m³/m，分别发生了1和4次崩塌，崩塌量分别为16.65和131.74 m³/m。从DM2断面来看，落水期冲刷量较少，但崩塌量相对较大，说明通过前期水流冲刷和冻融影响的累积作用，退水期发生崩岸的可能性相对较大，亦为崩岸多发时期。

3.2 河岸岸坡形态变化

依据前述的河岸稳定性计算结果，得到如图 5所示的河岸岸坡形态变化过程。需要指出的是，BSTEM模型计算时认为河岸发生崩塌后所有崩塌土体均被水流立刻冲走而没有在坡脚堆积，这与实际情况不符，所以在应用BSTEM模型对河岸崩塌进行模拟时，还需考虑崩塌后土体在坡脚的堆积。本文认为崩塌后土体堆积在坡脚，并且根据每次崩塌量的多少假设崩塌土体堆积在坡脚的角度不同，同时不考虑上游崩塌土体对该计算断面的影响，崩塌后土体堆积在坡脚，并对坡脚起到了一定保护作用，从而对后续河岸形态产生影响。

因此，本文首先用模型计算了未考虑崩塌土体在坡脚堆积的岸坡形态，然后根据上述处理方法，最终得到了考虑崩塌土体在坡脚堆积的河岸岸坡形态(概化为坡脚处岸坡角度因堆积变缓)。图 5给出了3个典型断面的河岸岸坡形态变化过程。从各典型断面的冲淤变化来看，DM1断面2009水文年河岸稳定系数均大于1，经计算，DM1断面边界水流切应力较小，坡脚冲刷量和横向冲刷宽度较少，河岸稳定性相对较高，无崩塌，坡脚冲刷总量为99.03 m³/m，其中洪水期为49.89 m³/m，约占总冲刷量的50 % 左右，说明河岸形态主要由洪水期冲刷决定。DM2断面坡脚冲刷总量为38.41 m³/m，从图 5可看出，2009水文年在洪水期和退水期分别发生过1次崩塌，崩塌宽度分别为0.88和1.96 m，崩塌总量为42.71 m³/m，约占冲刷与崩塌总量的53 % 左右；DM2在2009水文年计算的崩塌总宽度为2.84 m，与实测的平均崩塌宽度吻合较好(实测平均崩塌宽度约为3 m)。DM3坡脚冲刷量为338.13 m³/m，总体来看冲刷量较大；DM3在2009水文年洪水期和退水期分别发生6和4次崩塌，崩塌宽度分别为18.21和16.58 m，崩塌总量为383.42 m³/m，约占冲刷和崩塌总量的53 % 左右；DM3崩塌总宽度为34.79 m，与实测崩塌宽度吻合较好(实测崩塌宽度约为36 m)。总体而言，考虑冻融作用后，整个水文年计算的河岸岸坡形态变化和河岸崩塌宽度与实测结果基本吻合。

3.3 冻融作用对河岸崩退影响的定量分析

为进一步分析冻融作用对河岸崩塌的定量影响，针对4月考虑冻融和4月未考虑冻融2种情况，分别进行了河岸崩塌的模拟分析，前者简称工况1，后者简称工况2。考虑DM1未发生崩塌，故选取DM2、DM3计算结果，说明如下：

1) 在其它条件相同的情况下，工况1的冲刷量大于工况2，且冲刷崩塌总量大于工况2。从图 6、图 7可看出，2个断面4月以后不同时期工况1的河岸冲刷量均大于工况2，且工况1的崩塌总量大于工况2。从冲刷崩塌总量方面(图 8)来看，DM2断面工况1和工况2分别为81.12和57.48 m³/m，即考虑冻融影响后，冲刷崩塌总量增幅约为41 %；DM3断面分别为721.52和671.84 m³/m，即考虑冻融影响后，冲刷崩塌总量增幅约为7 %。这也从侧面反映出考虑冻融作用导致岸滩土体抗剪强度降低，坡脚冲刷量变大，河岸安全系数F_s减小，增加了崩塌发生的可能性，这与实际情况相符。考虑冻融影响后，DM2冲刷崩塌总量增幅相对较大，一方面是因为2009水文年DM2断面总体冲刷量相对较小；另一方面是考虑冻融影响后河岸崩塌次数(2次)比不考虑冻融影响时(1次)多1次，崩塌量增加了15.89 m³/m；因而，考虑冻融影响后，DM2断面冲刷崩塌总量增幅相对较大。

2) 在其它条件相同的情况下，冻融作用使河岸提前发生崩塌，且累计崩塌宽度大于未考虑冻融时的宽度。如图 9所示，在4月份以后，每次工况1发生崩塌以后，工况2随后才发生崩塌；DM2断面在涨水期，考虑冻融作用后计算的河岸最小安全系数(F_s =1.55)比不考虑的计算值(F_s=1.90)降低18.5 % 左右。从累计崩塌宽度(图 10)来看，未考虑冻融作用时，DM2、DM3两个断面崩塌宽度分别为1.04和30.24 m，考虑冻融作用影响后，累计崩塌宽度分别增加到2.84和34.79 m，分别增加了1.80和4.55 m，考虑冻融作用后的崩塌宽度大于未考虑时的宽度，且更接近实测河岸崩塌宽度(实测分别约为3和36 m)。值得注意的是，4月因冻融作用，土体抗剪强度发生变化后，河岸并未立刻发生崩塌，原因是河岸崩塌是累积性过程，此时来水量相对较小、河道内水位较低、冲刷幅度有限，但随着坡脚淘刷不断累积、岸坡坡度变陡，河岸安全系数减小，更易达到临界值发生崩塌，所以对后续崩塌埋下安全隐患。上述分析表明，冻融作用降低了河岸土体的力学强度，尤其是降低了土体抗剪强度，增大了河岸发生崩塌的可能，说明对于季冻区河流岸滩崩塌的模拟而言，需考虑冻融影响。

同时也应指出，我国国土幅员辽阔，北方不同纬度河流冻融过程存在一定差异，包括封冻时期、开河日期、岸滩土体类型等方面。对于河岸崩塌需重点关注的冻融作用土体抗剪强度变化而言，虽然具体数值可能存在差异，但冻融作用使得土体黏聚力呈降低趋势这一规律是基本一致的，因此，本研究以我国典型季节性冰冻河流松花江为例，定量分析了冻融作用对崩岸过程的影响，研究结果可为中国北方季节性冰冻河流岸滩崩塌及河道演变研究提供参考。另一方面，季节性冰冻河流水、冰、沙、土相互作用，崩岸机理复杂，本文仅对冻融的影响开展了初步概化研究，关于岸滩崩塌的问题还有待进一步深入研究。

4 结论

以松花江大顶子山航电枢纽下游近坝段典型断面为例，考虑土体冻融作用影响，结合坡脚冲刷，采用BSTEM模型计算了2009水文年断面的崩退过程，定量分析了冻融作用对河岸崩退的影响。主要研究结论为：

1) 涨水期河岸稳定性相对较高，洪水期和退水期河岸稳定性相对较低；DM1断面全年稳定无崩塌，DM2和DM3断面全年分别发生了2次和10次崩塌，其中洪水期和退水期为崩岸多发时期。

2) 在其它条件相同的情况下，冻融作用会使崩岸提前发生；考虑冻融作用累计崩塌宽度大于未考虑冻融影响时的宽度，且考虑冻融作用后的崩塌宽度更接近实测的河岸崩塌宽度，同时坡脚冲刷量、冲刷崩塌总量均大于未考虑冻融时的量；考虑冻融作用后，DM2和DM3全年累计冲刷崩塌总量增幅分别约为41 % 和7 %，与实际情况较接近

3) 综合考虑坡脚冲刷、冻融影响、土体分层和崩塌土体在坡脚的堆积等因素，断面岸坡形态和河岸崩塌总宽度模拟结果与实测情况基本吻合，该方法可为深入研究季冻区河流岸滩崩塌及河道治理提供一种有效技术手段。