一种基于土壤冻融过程的流域涵蓄水能力测算方法

摘要

本发明公开了一种基于土壤冻融过程的流域涵蓄水能力测算方法，首先进行基本数据处理；其次计算流域范围的土壤冻融水热状态，获得土壤冰含量比重；之后定量表达流域内各栅格点涵蓄水能力的时间分布，计算栅格型土壤涵蓄水能力指标；最后计算每一阶段子流域内栅格点的涵蓄水能力,绘制各子流域离散化的涵蓄水能力累积分布曲线，从而得到冻融区流域土壤涵蓄水能力的动态测算模型，对土壤冻融过程的流域涵蓄水能力进行测算。本发明将点尺度的融雪-冻土水、热耦合系统建立数值模型，并与流域栅格土壤涵蓄水能力建立联系，实现在空间尺度上流域涵蓄水能力指标动态设定并与实际情况吻合，减少寒区流域均一化所导致的误差。

说明书

技术领域

本发明属于水利工程的水循环技术领域，特别是涉及一种基于土壤冻融过程的流域涵蓄水能力测算方法。

背景技术

影响自然系统、生物系统、人类社会等方面的全球气候变化已经发生，并且还将加剧，它的任何走向都会对自然生态系统以及社会经济产生不可忽视的影响。在受到气候影响的各种环境系统中，高寒区生态系统是最敏感的区域之一，其活动过程正在深刻地影响着整个地球系统。高寒区水循环过程对人类经济活动及社会活动都有极其重要的影响，水循环活动造成的洪灾、旱灾以及水资源的分配不均给全世界造成了巨大损失。为了能够了解水循环的过程，进而通过科学手段尽量避免水旱灾害，水文学对水循环过程做出了很多深入研究。现代水文学通过实验、物理过程分析总结水循环规律，并在此基础上结合数学物理方法建立能描述水循环过程的水文学模型，期待通过模型能准确地了解水循环的各相关过程。从二十世纪50年代到现代，水文学模型研究越来越趋近于理论化，模型从最初的统计模型，逐步过渡到半物理的概念性模型以及全动力过程描述的物理模型，通过这些模型的研究有力地推动了人类对水文循环过程的了解，对指定趋利避害的规则提供了大量的科学依据。目前水文模型主要的研究思路是分别研究水文产流和汇流，其中产流过程是核心，它决定了降雨后形成的洪水总量。而在产流过程中流域土壤的涵蓄水能力是决定性因素，降雨、积雪融化，冻土升温冰融化扣除土壤涵养及持蓄水的能力后即可得到产流量。基于此，当前的水文模型均引入了土壤涵蓄水能力相关的参数，蓄水容量是反映土壤涵蓄水能力的有效表达。

流域下垫面情况是确定土壤涵蓄水能力的重要因子，我国著名的新安江模型提出了反映土壤涵蓄水能力空间分布的蓄水容量曲线，隐含地表达了地形变化和土壤下垫面变化对流域蓄水容量的影响；寒区流域冻土层土壤的季节性冻融变化是下垫面变化的主要影响因素之一，其冰含量的季节变化过程是土壤涵蓄水能力的重要因子；SHAW、SiB、BATS、VIC等模型针对寒区冻土冻融过程模型中的土壤结构参数、水热动力学参数等，描述多年冻土区活动层水热输运情况。

虽然目前有关于寒区冻融机理的研究，但多着重于水热耦合过程描述，较少关注冻融转化对流域水文过程特别是产流过程的影响。有关于土壤冻融活动对流域土壤涵蓄水能力影响的研究，但研究的成果多采用表现为经验性曲线来描述，比如新安江模型提出了一个抛物型的蓄水容量曲线，其中有经验性的参数需要率定，模型结构理论化程度不高。寒区涵蓄水能力本质上与流域气候和土壤质地密切相关，但目前还未见结合融雪―冻土水、热耦合与土壤类型来共同确定土壤涵蓄水能力的报道。

发明内容

本发明所要解决的技术问题是为克服现有技术的不足，而提供一种基于土壤冻融过程的流域涵蓄水能力测算方法，可通过建立地理地形数据、土壤类型信息及冻土冻融季节性活动与流域栅格土壤涵蓄水能力之间的关系，实现在时间、空间尺度上流域涵蓄水能力指标动态设定并与实际情况的吻合，从而得到流域涵蓄水能力测算模型，将有效地推动流域水文模型的科学发展。

本发明为解决上述技术问题，提出以下技术方案：

一种基于土壤冻融过程的流域涵蓄水能力测算方法，包括以下步骤：

步骤一，进行冻融区流域气象及下垫面基本数据处理，包括：气象数据采集、地理数据采集、计算栅格指标、核算栅格土壤信息；

步骤二，计算流域范围的土壤冻融水热状态，具体为：根据水分传递、热量传导理论，分别计算流域栅格点土壤温度、土壤液态水含量，获得土壤冰含量比重；

步骤三，根据步骤二获得的土壤冰含量比重，定量表达流域内各栅格点涵蓄水能力的时间分布，计算栅格型土壤涵蓄水能力指标；

步骤四，根据步骤三的结果，计算每一阶段子流域内栅格点的涵蓄水能力,按照从小到大的顺序排列，再根据栅格面积比绘制各子流域离散化的涵蓄水能力累积分布曲线，从而得到冻融区流域土壤涵蓄水能力的动态测算模型，对土壤冻融过程的流域涵蓄水能力进行测算。

进一步的，本发明的一种基于土壤冻融过程的流域涵蓄水能力测算方法，在步骤一中：

所述的气象数据采集，包括收集流域内气温、降水、湿度、气压数据；

所述的地理数据采集，包括收集流域地形数据、栅格高程数据、栅格土壤特性数据；

所述的计算栅格指标，包括计算栅格内的坡度、坡向以及河长，划分流域边界，生成子流域拓扑关系与流域水系；

所述的核算栅格土壤信息，包括计算栅格内各类土壤的比例，叠加栅格内各类土壤的组成，量定土壤水热参数。

进一步的，本发明的一种基于土壤冻融过程的流域涵蓄水能力测算方法，步骤二中所述的计算流域范围的土壤冻融水热状态具体包括以下步骤：

首先，基于能量平衡、水量平衡原理，得到雪盖温度和析出水量：

$>T_{\mathrm{snow}}^{t+1}=T_{\mathrm{snow}}^t+\left(\frac{Q_{\mathrm{net}}-Q_{\mathrm{phase}}}{{\rho }_w{C}_i{W}_{\mathrm{ice}}}\right)---\left(1\right)$>

$>W_{\mathrm{liq}}^{t+1}=W_{\mathrm{liq}}^t+[P-E+|\frac{Q_{\mathrm{net}}-Q_{\mathrm{change}}}{{\rho }_w{\lambda }_f}\left|\right]---\left(2\right)$>

式中：W_liq为雪层内液态水含量；T_snow为雪层内温度；Q_net为雪层内的净输入能量；Q_phrase为水分相变能量；Q_change为积雪融化能量；P为降雨量；E为蒸发损失量；ρ_w为水密度；C_i为冰比热容；W_ice为雪水当量；λ_f为融化热，t为计算时刻；

其次，利用试验方式获得栅格内各类土壤水热参数C_S、L_f、λ、K、D，结合式(1)～(2)作为输入条件，参与土壤水热联立式(3)～(4)计算，得到土壤冰含量比重θ_Bi；

$>-\frac{\mathrm{\Delta t}{\lambda }_{i-1/2}}{{\left({\mathrm{\Delta Z}}_i\right)}^2}{T}_{i-1}^{t+1}+(C_{\mathrm{Si}}+\frac{\mathrm{\Delta t}{\lambda }_{i-1/2}}{{\left(\Delta Z_i\right)}^2}+\frac{\Delta {\mathrm{t\lambda }}_{i+1/2}}{{\left({\mathrm{\Delta Z}}_i\right)}^2}){T_i}^{t+1}-\frac{\mathrm{\Delta t}{\lambda }_{i+1/2}}{{\left({\mathrm{\Delta Z}}_i\right)}^2}{T}_{i+1}^{t+1}=C_S{T_i}^t+L_f{\rho }_B({\theta }_{\mathrm{Bi}}^{t+1}-{\theta }_{\mathrm{Bi}}^t)---\left(3\right)$>

$>\left(\begin{array}{c}-\frac{\mathrm{\Delta t}{D}_{i-1/2}}{{\left(\Delta Z_i\right)}^2}{\theta }_{l_{i-1}}^{t+1}+(1+\frac{\mathrm{\Delta t}{D}_{i-1/2}}{{\left(\Delta Z_i\right)}^2}+\frac{\mathrm{\Delta t}{D}_{i+1/2}}{{\left({\mathrm{\Delta Z}}_i\right)}^2}){\theta }_{l_i}^{t+1}-\frac{\mathrm{\Delta t}{D}_{i+1/2}}{{\left(\Delta Z_i\right)}^2}{\theta }_{l_{i+1}}^{t+1}\\ ={\theta }_{l_i}^t+\frac{\mathrm{\Delta t}}{\Delta Z_i}(K_{i-1/2}-K_{i+1/2})-\frac{{\rho }_B}{{\rho }_l}({\theta }_{\mathrm{Bi}}^{t+1}-{\theta }_{\mathrm{Bi}}^t)\end{array}\right)---\left(4\right)$>

式中：T为土壤温度；ρ_B为冰密度；ρ_l为水密度；θ_l为土壤液态水含量；Z为栅格空间坐标，L_f为相变潜热，C_S为土壤体积热容量，D为土壤水扩散系数，K为土壤水导水系数，λ为土壤导热系数，i=1,2,…M，M为栅格型土壤上至下划分的层数。

进一步的，本发明的一种基于土壤冻融过程的流域涵蓄水能力测算方法，步骤三中所述的计算栅格型土壤涵蓄水能力指标，具体过程包括：

1)利用时段栅格内各层土壤冰含量比重，得到土壤冰含量比重的年内变化：

$>f_B(t,x,y,h)=\underset{i=1}{\overset{i=M}{\Sigma }}{\theta }_{\mathrm{Bi}}^t---\left(5\right)$>

式中：t为计算时刻；x、y、h分别为栅格行、列及高程数据，i=1,2,…M，M为栅格型土壤上至下划分的层数，θ_Bi为土壤冰含量比重；

2)利用步骤1）的土壤冰含量比重的年内变化得到土壤冻融状态冰含量比重分布曲线，获得计算深度H_soil内土壤蓄水容量的变化量dW：

dW=[f_B(t+1,x,y,h)-f_B(t,x,y,h)]H_soil          (6)

3)利用式(6)计算栅格型土壤涵蓄水能力指标W_M：

W_M(t+1,x,y,h)=W_M(t,x,y,h)+dW          (7)

流域中各栅格点土壤水热条件不同，土壤特征各异，通过式(7)得到每个栅格点上相应的涵蓄水能力指标。

本发明与现有技术相比的显著优点为：一是本发明能够实际计算出流域中各点、各季节的涵蓄水能力，反应流域内涵蓄水能力的时间、空间分布特征；二是本发明能够考虑流域中土壤的类型，水、热条件，冻融情况将导致不同的涵蓄水能力曲线；三是在基于现代测量技术所得的前期资料支持下，本发明所得模型应用于季节冻土区实时洪水预报将可大幅度提高寒区流域洪水预报的精度，为流域防洪预警等提供可靠的科学支撑。本发明既适用于集总式水文模型，又适用于分布式水文模型，将有效地推动流域水文科学研究的深入发展。

附图说明

图1是本发明提出的流域涵蓄水能力测算方法的流程示意图。

图2是格网点水流方向编码示意图。

图3是寻找洼地区域或平坦区域出口的方法示意图。

图4是最短流程算法示意图。

图5是一个简单的DEM及其计算结果。

图6是DEM水文分析流程图。

图7是三江源流域水系图。

图8是三江源子流域分割及流域边界图。

图9是冻土活动层内土壤冰含量比重示意图。

图10是分层涵蓄水能力计算流程图。

图11是三江源流域涵蓄水能力曲线图。

图12是三江源流域径流精度分析图。

具体实施方式

下面结合附图和具体应用实施例对本发明作进一步的详细说明。

结合图1，本发明提出的一种基于土壤冻融过程的流域涵蓄水能力测算方法，总体思路是：步骤1，流域气象及下垫面基本数据处理；步骤2，计算流域范围的土壤冻融水热状态；步骤3，定量表达流域内各栅格点涵蓄水能力的时间分布；步骤4，绘制流域涵蓄水能力动态累计分布曲线，最终构建寒区流域的涵蓄水能力测算模型。

实施本发明的具体步骤包括：

步骤1，流域气象及下垫面基本数据处理；包括气象数据采集、地理数据采集、计算栅格指标、核算栅格土壤信息。

1、气象数据采集，包括收集流域内气温、降水、湿度、气压数据；在此不再赘述。

2、地理数据采集，包括收集流域地形数据、栅格高程数据、栅格土壤特性数据；其中土壤特性数据收集包括：搜集流域范围内土壤的相关资料，包括土壤质地、土壤类型、土壤饱和含水量，土壤孔隙度、土壤导水性质等方面数据。

3、计算栅格指标，包括计算栅格内的坡度、坡向以及河长，划分流域边界，生成子流域拓扑关系与流域水系；以DEM(Digital Elevation Model，数值高程模型)数据为基础，生成子流域拓扑关系、流域水系、划分流域边界，按照以下方式进行：

1)水流方向矩阵的计算：

结合图2，采用D8(Deterministic Eight-neighbours)算法进行计算，将被处理的格网点K同其最邻近的8个格网点之间的坡降进行比较，其中与落差最大的一个格网点中心之间连线的方向便定义为被处理格网点K的水流方向，并且规定一个格网点的水流方向用一个特征码表示。有效的水流方向定义为东、东南、南、西南、西、西北、北和东北，并分别用1、2、3、4、5、6、7和8这8个有效特征码表示。

被处理格网点K同相邻8个格网点之间单位距离的高差为：

MD=Z/D               (8)

式中：MD为两个格网点之间的单位距离的高差，表示地形坡度；Z为两个格网点之间的高程差；D为两个格网点中心之间的距离。

确定水流方向的具体步骤如下：

①对所有DEM边缘的格网，赋以指向边缘的方向值0；

②对所有在第一步中未赋方向值的格网，计算其对8个邻域格网的单位距离的高差值，确定具有最大落差值的格网，执行以下步骤：

A.如果该格网与相邻8个邻域格网的最大落差小于0，则赋以负值以表明此格网方向未定（这种情况在经洼地处理的DEM中不会出现）。

B.如果该格网的高程与相邻8个邻域格网最大落差大于或等于0，且最大落差只有一个，则该格网的水流方向赋以指向最大落差的方向。

C.如果该格网的高程与相邻8个邻域格点最大落差大于0，且最大落差有多个，则该格网的水流方向在逻辑上以查表方式确定，也就是说，如果中心格网在一条边上的三个邻域点有相同的落差，则中间的格网方向被作为中心格网的水流方向，如果中心格网的相对边上有两个邻域格网落差相同，则任选一格网方向作为水流方向。

D.如果该格网的高程与相邻8个邻域格点最大落差等于0，且最大落差有多个，则以这些0值所对应的方向值相加。

③对①、②步没有赋以负值，0，1，2，3，…，8的每一格网，检查对中心格网有最大落差值的邻域格网。如果邻域格网的水流方向值为1，2，3，…，8，且此方向没有指向中心格网，则以此格网的方向值作为中心格网的方向值；

④重复③，直到所有格网都被赋值。

2)洼地处理：

被高程较高的区域围绕的洼地是应用数字高程流域水系生成模型生成流域水系的一大障碍，因为在决定水流方向以前，必须先将洼地填充。有些洼地是DEM生成过程中带来的数据错误，但另外一些却表示了真实的地形。一些研究试图通过平滑处理消除洼地，但平滑处理只能处理浅的洼地，较深的洼地仍然无法处理。处理洼地的另一种方法是通过将洼地的每一格网点赋以洼地边缘最小值，从而达到消除洼地的目的。这些方法都需要对地形数据进行修改，可能会产生一些不合理的方向阵，下面将具体的介绍此种情况的处理。

根据水流特点，通过对洼地区域和平坦区域标志，利用最短流程算法，修改洼地区域和平坦区域水流方向矩阵，使研究区域中的水流能够通过洼地区域和平坦区域。算法的步骤如下：

①结合图3，先找出洼地区域或平坦区域的边缘格网，然后在边缘格网中找出高程最低的格网点。如果不存在最低格网，则将洼地区域或平坦区域向外围扩大一个格网继续寻找，直到找到这个最低格网；

②结合图4，如果边缘格网中高程最低的格网点不为洼地区域内的格网点，则将边缘格网中高程最低格网点作为洼地区域或平坦区域的水流出口点；

③结合图4，以洼地区域或平坦区域的水流出口点为起点，修改洼地区域或平坦区域内高程低于水流出口格网点的水流方向阵。

3)水流累计矩阵的计算：

结合图5，区域水流累计矩阵表示区域地形每一点的流水累积量。其基本思想是，假定以规则格网表示的数字地面高程模型的每一点有一个单位的水量，按照水流从高往低流的规律，根据区域地形的水流方向矩阵，计算每一格网点流过的水量数值，便可得到该区域的水流累积矩阵。下面给出一个从原始DEM矩阵计算出相应水流方向矩阵及水流累积矩阵的范例。

从(a)到(b)的变换方法前面已经介绍，下面简要解释从(b)到(c)的变换方法：

以(b)左上单元格(2)为例，因为周围没有那个单元格的水流流入该单元格（即该单元格为流域边界），所以在(c)中相应填入0，同理可知(c)中第一行所填的数均为0；以(b)中第二行第二单元格为例，由于其左上方单元格的水流将流入此单元格，因此，在(c)中相应为指填入1；以(b)中第二行第四列单元格为例，其左上方和正上方单元格的水流都将流入其中，即有两个单位的水流汇入此单元格，因此，在(c)中相应的单元格填入2；再以(b)中第三行第四列单元格为例，其左上方和正上方的水流将汇入其中；另外，左上方已经有一个单位的入流，正上方单元格则有两个单位的入流，再加上两个单元格自身的水量，此时，(c)中相应单元格所填入的数应该是5。以此类推，即可由(b)得到(c)，从而生成水流累积矩阵。

4)基本水文分析，其中包括汇水面积计算，流域分水线识别，河网生成，建立流域水系拓扑关系等四部分。

①汇水面积计算：

按最大落差规则确定的水流路径可以非常方便地计算出指定格网点以上的流域汇水面积，若以格网数表示面积的话，则该格网点的汇水面积数值是该格网点以上汇入该格网点的格网点数目。在算法上可利用递归算法实现，由指定点开始按逆水流方向向上搜索迭代，即可得到该集水流域内任意格网点的汇水面积，其结果如图5(c)表示的水流累积矩阵。

②流域分水线的识别：

给定流域的主要入口和出口的断面位置，即它们所在栅格单元的行列坐标。一旦两者的位置确定，根据它们的汇水面积大小，程序可自动搜索从而勾划出流域边界并计算出流域面积。

③河网生成：

河网生成分为三个步骤：确定流域边界内的水道；裁减小于某一临界长度的河段；生成河道编码。

首先，给定最小河道给养面积阀值，小于该值的汇水面积不可能产生足够的径流而形成水道。流域范围内汇水面积超过该阀值的那些格网点即定义为水道。

其次，给定最小河道长度，若一级河道的累积长度小于该长度，则该水道被裁减。由第一步生成的一些河道可能会很短，那些很短的一级水道很可能是伪水道、位于河谷两边的凹痕或沟壑的出口，需要将它们裁减除去。

最后，确定河道级数和河段长度。根据流域出口断面确定干流河道，流入干流的河道被定为一级支流，流入一级支流的河道定为二级支流，依此类推，确定所有河道的编码。同时可以确定各级支流汇入上一级河道的节点，对河网所有节点进行编码，这样定义的节点编码可用于水文汇流计算或河网数据库的构建。

④流域水系拓扑关系：

一旦生成联结完好的河网，根据各河网节点，可以确定相应各支流的流域边界线，从而建立河网节点、河段和子流域的拓扑关系，包括河段坡度、高程值、上游集水面积与侧向集水面积及相互连接的拓扑信息。一方面，河网与子流域边界等空间信息是以栅格形式存储，这样易于用GIS软件作可视化显示；另一方面，河段或子流域的拓扑关系还以表格形式存储，有利于数字水文模型的调用。图6为利用DEM进行水文分析的流程。

4、核算栅格土壤信息，包括计算栅格内各类土壤的比例，叠加栅格内各类土壤的组成，量定土壤水热参数。即针对研究区同一栅格范围内存在多种土壤类型，且各成分所占比例不同，将多种土壤类型性状参数数据与其对应的成分比例在同一栅格内叠加，得到栅格内综合的土壤水热参数。

步骤2，流域范围内核算土壤冻融水热状态，是指以DEM数据分析过程中栅格内土壤类型及其水热参数为基础，核算土壤中的冰含量比重，其计算遵循以下步骤进行：

土壤冻融水热状态的计算：

①计算雪盖温度和析出水量：

根据能量平衡原理，计算雪盖内获得的净能量，再利用水量平衡原理，估算积雪层内的雪水当量，算法的步骤如下：

(a)运用能量平衡法，逐日计算雪盖内的净能量

Q_net=Q_r+Q_s+Q_l+Q_p+Q_g        (9)

式中：Q_net为雪层内的净输入能量；Q_r为净辐射通量；Q_s为感热通量；Q_l为潜热通量；Q_p为降水带来的热通量；Q_g为地热通量。

当Q_net>0时，雪盖内无水分相变发生，Q_phrase为零，Q_net提供雪盖内的温度升高；

当Q_net<0时，积雪层内持续降温（Q_phrase=ρ_wλ_fW_liq-old），将积雪层内液态水冻结。

(b)雪盖内温度为

$>T_{\mathrm{snow}}^{t+1}=T_{\mathrm{snow}}^t+\left(\frac{Q_{\mathrm{net}}-Q_{\mathrm{phase}}}{{\rho }_w{C}_i{W}_{\mathrm{ice}}}\right)---\left(1\right)$>

式中：T_snow为雪层内温度；ρ_w为水密度；C_i为冰比热容；W_ice为雪水当量；Q_phrase为水分相变能量；W_liq-old为雪层内上一时段末的液态含水量；λ_f为融化热；t为计算时刻。

(c)运用水量平衡法，得到计算时段内，积雪层析出的水量

$>W_{\mathrm{liq}}^{t+1}=W_{\mathrm{liq}}^t+[P-E+|\frac{Q_{\mathrm{net}}-Q_{\mathrm{change}}}{{\rho }_w{\lambda }_f}\left|\right]---\left(2\right)$>

式中：W_liq为雪层内液态水含量；Q_change为积雪融化能量；P为降雨量；E为蒸发损失量。

②计算土壤冰含量比重

根据非饱和土壤水流运动理论以及热传导理论，结合①中(b)，(c)步骤得到建立水热耦合模型模拟冻土内的水分迁移过程。

(a)冻土水热耦合模型的定解条件为初始条件和边界条件，对应模型的土壤含水量上边界条件，考虑①中式(2)积雪层析出的水量以及蒸发量；对应地表温度上边界条件，以①中式(1)雪层内温度代替；下边界条件为通量为零的第二类边界条件给定；初始条件中土层的含水量和温度做均一化处理。

(b)利用试验方式获得栅格内各类土壤水热参数C_S、L_f、λ、K、D等，结合定解条件，参与土壤水、热联立计算，得到土壤冰含量比重θ_Bi

$>-\frac{\mathrm{\Delta t}{\lambda }_{i-1/2}}{{\left({\mathrm{\Delta Z}}_i\right)}^2}{T}_{i-1}^{t+1}+(C_{\mathrm{Si}}+\frac{\mathrm{\Delta t}{\lambda }_{i-1/2}}{{\left(\Delta Z_i\right)}^2}+\frac{\Delta {\mathrm{t\lambda }}_{i+1/2}}{{\left({\mathrm{\Delta Z}}_i\right)}^2}){T_i}^{t+1}-\frac{\mathrm{\Delta t}{\lambda }_{i+1/2}}{{\left({\mathrm{\Delta Z}}_i\right)}^2}{T}_{i+1}^{t+1}=C_S{T_i}^t+L_f{\rho }_B({\theta }_{\mathrm{Bi}}^{t+1}-{\theta }_{\mathrm{Bi}}^t)---\left(3\right)$>

$>\left(\begin{array}{c}-\frac{\mathrm{\Delta t}{D}_{i-1/2}}{{\left(\Delta Z_i\right)}^2}{\theta }_{l_{i-1}}^{t+1}+(1+\frac{\mathrm{\Delta t}{D}_{i-1/2}}{{\left(\Delta Z_i\right)}^2}+\frac{\mathrm{\Delta t}{D}_{i+1/2}}{{\left({\mathrm{\Delta Z}}_i\right)}^2}){\theta }_{l_i}^{t+1}-\frac{\mathrm{\Delta t}{D}_{i+1/2}}{{\left(\Delta Z_i\right)}^2}{\theta }_{l_{i+1}}^{t+1}\\ ={\theta }_{l_i}^t+\frac{\mathrm{\Delta t}}{\Delta Z_i}(K_{i-1/2}-K_{i+1/2})-\frac{{\rho }_B}{{\rho }_l}({\theta }_{\mathrm{Bi}}^{t+1}-{\theta }_{\mathrm{Bi}}^t)\end{array}\right)---\left(4\right)$>

式中：T为土壤温度；ρ_B为冰密度；ρ_l为水密度；θ_l为土壤液态水含量；Z为栅格空间坐标。

步骤3，以栅格为单位，结合土层内冰含量比重的季节性变化，建立土壤涵蓄水能力数学模型，计算栅格型土壤涵蓄水能力指标，其计算遵循以下步骤进行：

①利用时段栅格内各层（1～M层）土壤冰含量比重，得到土壤冰含量比重的年内变化：

$>f_B(t,x,y,h)=\underset{i=1}{\overset{i=M}{\Sigma }}{\theta }_{\mathrm{Bi}}^t---\left(5\right)$>

式中：x、y、h分别为栅格行、列及高程数据；f_B为栅格内逐层土壤冰含量比重。

②利用土壤冻融状态冰含量比重分布曲线，获得计算深度内土壤蓄水容量的变化量：

dW=[f_B(t+1,x,y,h)-f_B(t,x,y,h)]H_soil       (6)

式中：dW为相邻时刻土壤蓄水容量的变化量；S_zm为参数；H_soil为计算深度。

③利用式(6)计算栅格型土壤涵蓄水能力指标：

W_M(t+1,x,y,h)=W_M(t,x,y,h)+dW        (7)

式中：W_M为栅格型土壤涵蓄水能力指标。

流域中各栅格点土壤水热条件不同，土壤特征各异，通过式(7)得到每个栅格点上相应的涵蓄水能力指标。

步骤4，结合流域分割，得到各子流域内栅格涵蓄水能力序列，根据栅格面积比建立各子流域离散化的涵蓄水能力累积分布曲线，其计算遵循以下步骤进行：

①结合公式(3)～(4)，土壤水热耦合计算模型，计算流域范围内各栅格的涵蓄水能力；

②结合流域水系构成、子流域划分范围，以及各子流域的拓扑关系，划分各子流域内所组成的栅格；

③统计各子流域内栅格涵蓄水能力，按照从小到大顺序排列，根据栅格面积比建立各子流域离散化的涵蓄水能力累积分布曲线，得到流域内的涵蓄水能力空间分布。

下面以青藏高原三江源流域为应用实施区域实例，进一步对本发明的具体实施方式作如下详细说明。

三江源区为典型的高原大陆性气候，表现为冷热两季交替、干湿两季分明，年温差小、日温差大、日照时间长、辐射强烈等气候特征。三江源区全年平均气温通常在-5.6℃到-3.8℃之间。降水量高度集中，水热同季，降水地域差异很大，年平均降水量262.2-772.8mm，6-9月的降水量约占全年总降水量的75%，年蒸发量达到730-1700mm。黄河干流唐乃亥以上的黄河源区、长江干流直门达以上长江源区。

选取三江源流域为应用实施区域，利用本发明的全流域涵蓄水能力模型，结合DEM数据驱动流域水系生成、子流域范围分割以及地形指数算法，利用三江源流域的水文气象数据，土壤类型栅格数据和数字高程数据等，驱动全流域涵蓄水能力模型，模拟三江源流域水文过程。

步骤1：流域下垫面基本数据处理

①采集三江源流域DEM资料

三江源流域，经度89.45°～102.23°，纬度31.39°～36.12°，海拔2800～6564m。由公开的30米分辨率的ASTGTM数据，采集研究区范围内地形栅格数据；

②水流方向矩阵计算

采用D8算法按照如下步骤进行计算：

首先，对所有DEM边缘的格网，赋以指向边缘的方向值0；

其次，对所有在前一步中未赋方向值的格网，计算其对8个邻域格网的单位距离的高差值，确定具有最大落差值的格网；

再次，对前两步没有赋以负值，0，1，2，3，…，8的每一格网，检查对中心格网有最大落差值的邻域格网。如果邻域格网的水流方向值为1，2，3，…，8，且此方向没有指向中心格网，则以此格网的方向值作为中心格网的方向值；

最后，重复第三步，直到所有格网都被赋值。

③洼地处理

在处理水流方向之前，根据水流特点，通过对洼地区域和平坦区域标志，利用最短流程算法，修改洼地区域和平坦区域水流方向矩阵，使研究区域中的水流能够通过洼地区域和平坦区域，按如下步骤执行洼地处理，计算水流方向：

首先，结合图3，先找出洼地区域或平坦区域的边缘格网，然后在边缘格网中找出高程最低的格网点。如果不存在最低格网，则将洼地区域或平坦区域向外围扩大一个格网继续寻找，直到找到这个最低格网；

其次，结合图4，如果边缘格网中高程最低的格网点不为洼地区域内的格网点，则将边缘格网中高程最低格网点作为洼地区域或平坦区域的水流出口点；

最后，结合图4，以洼地区域或平坦区域的水流出口点为起点，修改洼地区域或平坦区域内高程低于水流出口格网点的水流方向阵，最终得到三江源流域水流方向矩阵形式。

④水流累计矩阵计算

根据三江源流域栅格水流方向阵计算结果，计算相应的水流累计矩阵，其计算步骤如下：

以图5(b)左上单元格(2)为例，因为周围没有那个单元格的水流流入该单元格（即该单元格为流域边界），所以在5(c)中相应填入0，同理可知5(c)中第一行所填的数均为0；以5(b)中第二行第二单元格为例，由于其左上方单元格的水流将流入此单元格，因此，在5(c)中相应为指填入1；以5(b)中第二行第四列单元格为例，其左上方和正上方单元格的水流都将流入其中，即有两个单位的水流汇入此单元格，因此，在5(c)中相应的单元格填入2；再以5(b)中第三行第四列单元格为例，其左上方和正上方的水流将汇入其中；另外，左上方已经有一个单位的入流，正上方单元格则有两个单位的入流，再加上两个单元格自身的水量，此时，5(c)中相应单元格所填入的数应该是5。以此类推，即可由5(b)得到5(c)，从而生成水流累积矩阵，最终得到三江源流域累计矩阵。

⑤流域水系生成

首先，给定三江源流域的主要入口和出口的断面位置（长江源流域出口—直门达：经度：97°13′，纬度：33°02′；黄河源流域出口—唐乃亥：经度：100°09′，纬度：35°30′）。根据流域的汇水面积大小，按照流域累计矩阵搜索从而勾划出流域水系；

其次，给定最小河道给养面积阀值，小于该值的汇水面积不可能产生足够的径流而形成水道。流域范围内汇水面积超过该阀值的那些格网点即定义为水道。

再次，给定最小河道长度，若一级河道的累积长度小于该长度，则该水道被裁减。由前一步生成的一些河道可能会很短，那些很短的一级水道很可能是伪水道、位于河谷两边的凹痕或沟壑的出口，需要将它们裁减除去。

最后，确定河道级数：四级，河段长度：1497～42275m不等。根据流域出口断面：直门达水文站和唐乃亥水文站，确定干流河道，流入干流的河道被定为一级支流，流入一级支流的河道定为二级支流，依此类推，确定所有河道的拓扑关系，三江源流域水系，如图7所示。

⑥流域水系拓扑关系

生成联结完好的河网水系，根据河网节点，确定相应各支流的流域边界线，从而建立河网节点、河段和子流域的拓扑关系，包括河段坡度、高程值、上游集水面积与侧向集水面积及相互连接的拓扑信息，从而勾划出流域边界并计算出流域面积，三江源流域共划分出83个子流域，如图8所示。

步骤2：流域范围内核算土壤冰含量比重

①计算雪盖温度和析出水量

首先，运用能量平衡法，逐日计算雪盖内的净能量

Q_net=Q_r+Q_s+Q_l+Q_p+Q_g         (9)

式中：Q_net为雪层内的净输入能量；Q_r为净辐射通量；Q_s为感热通量；Q_l为潜热通量；Q_p为降水带来的热通量；Q_g为地热通量。

当Q_net>0时，雪盖内无水分相变发生，Q_phrase为零，Q_net提供雪盖内的温度升高；

当Q_net<0时，积雪层内持续降温（Q_phrase=ρ_wλ_fW_liq-old），将积雪层内液态水冻结。

此时，雪盖内的温度达到

$>T_{\mathrm{snow}}^{t+1}=T_{\mathrm{snow}}^t+\left(\frac{Q_{\mathrm{net}}-Q_{\mathrm{phase}}}{{\rho }_w{C}_i{W}_{\mathrm{ice}}}\right)---\left(1\right)$>

式中：T_snow为雪层内温度；ρ_w为水密度；C_i为冰比热容；W_ice为雪水当量；Q_phrase为水分相变能量；W_liq-old为雪层内上一时段末的液态含水量；λ_f为融化热；t为计算时刻。

其次，运用水量平衡法，得到计算时段内积雪层析出的水量

$>W_{\mathrm{liq}}^{t+1}=W_{\mathrm{liq}}^t+[P-E+|\frac{Q_{\mathrm{net}}-Q_{\mathrm{change}}}{{\rho }_w{\lambda }_f}\left|\right]---\left(2\right)$>

式中：W_liq为雪层内液态水含量；Q_change为积雪融化能量；P为降雨量；E为蒸发损失量。

然后，将上述雪层内温度和积雪层析出的水量扣除蒸发量后，作为土壤表层温度及含水量的输入条件；下边界条件为通量为零的第二类边界条件给定；初始条件中土层的含水量和温度做均一化处理。基于非饱和土壤水流运动理论以及热传导理论，建立水热耦合模型模拟冻土内的水分迁移过程。利用试验方式获得栅格内土壤以粘土为主，其中相变潜热L_f为334.56J/kg；土壤体积热容量C_S、土壤水扩散系数D和土壤水导水系数K与土壤的冻融状态有关：

未冻土：

K(θ_l)⁺=2.9657×10^-3θ_l^11.2593          (10)

D(θ_l)⁺=3.03×10^-2θ_l^4.4413           (11)

$>C_u=\frac{C_{\mathrm{su}}+{\mathrm{\theta C}}_w}{1+\theta }{\rho }_u---\left(12\right)$>

冻土：

$>K{\left({\theta }_l\right)}^{-}=K{\left({\theta }_l\right)}^{+}/{10}^{{10\theta }_i}---\left(13\right)$>

$>D{\left({\theta }_l\right)}^{-}=D{\left({\theta }_l\right)}^{+}/{10}^{{10\theta }_i}---\left(14\right)$>

$>C_f=\frac{C_{\mathrm{df}}+{\theta }_i{C}_i+{\theta }_l{C}_w}{1+\theta }{\rho }_f---\left(15\right)$>

式中：θ_l为土壤液态水含量；ρ_f和ρ_u分别为冻、融状态下土壤密度；θ为土壤含水率（液态水和冰）；C_su为未冻时土壤体积热容量(0.84kJ/(kg·℃))；C_df为冻时土壤体积热容量(0.77kJ/(kg·℃))；C_i为冰比热(2.135kJ/(kg·℃))。

当土壤含水量在15～41%之间时，土壤导热系数λ取值与季节有关，其中12月至次年5月，基本维持在1.61W/(m.℃)左右不变（表1所示）。

表1土壤导热系数季节变化列表(单位:W/(m.℃))

土壤类型5月6月7月8月9月10月11月12月粘土1.611.530.651.141.040.861.441.16

最后，构建土壤水、热耦合模型，计算土壤冰含量比重θ_Bi

$>-\frac{\mathrm{\Delta t}{\lambda }_{i-1/2}}{{\left({\mathrm{\Delta Z}}_i\right)}^2}{T}_{i-1}^{t+1}+(C_{\mathrm{Si}}+\frac{\mathrm{\Delta t}{\lambda }_{i-1/2}}{{\left(\Delta Z_i\right)}^2}+\frac{\Delta {\mathrm{t\lambda }}_{i+1/2}}{{\left({\mathrm{\Delta Z}}_i\right)}^2}){T_i}^{t+1}-\frac{\mathrm{\Delta t}{\lambda }_{i+1/2}}{{\left({\mathrm{\Delta Z}}_i\right)}^2}{T}_{i+1}^{t+1}=C_S{T_i}^t+L_f{\rho }_B({\theta }_{\mathrm{Bi}}^{t+1}-{\theta }_{\mathrm{Bi}}^t)---\left(3\right)$>

$>\left(\begin{array}{c}-\frac{\mathrm{\Delta t}{D}_{i-1/2}}{{\left(\Delta Z_i\right)}^2}{\theta }_{l_{i-1}}^{t+1}+(1+\frac{\mathrm{\Delta t}{D}_{i-1/2}}{{\left(\Delta Z_i\right)}^2}+\frac{\mathrm{\Delta t}{D}_{i+1/2}}{{\left({\mathrm{\Delta Z}}_i\right)}^2}){\theta }_{l_i}^{t+1}-\frac{\mathrm{\Delta t}{D}_{i+1/2}}{{\left(\Delta Z_i\right)}^2}{\theta }_{l_{i+1}}^{t+1}\\ ={\theta }_{l_i}^t+\frac{\mathrm{\Delta t}}{\Delta Z_i}(K_{i-1/2}-K_{i+1/2})-\frac{{\rho }_B}{{\rho }_l}({\theta }_{\mathrm{Bi}}^{t+1}-{\theta }_{\mathrm{Bi}}^t)\end{array}\right)---\left(4\right)$>

式中：T为土壤温度；ρ_B为冰密度；ρ_l为水密度；Z为栅格空间坐标。

将雪盖-土壤中水、热数据转化为相应土壤冰含量比重。

步骤3：流域涵蓄水能力分布的定量表达

①利用时段栅格内各层（1～M层）土壤冰含量比重，得到土壤冰含量比重的年内变化：

$>f_B(t,x,y,h)=\underset{i=1}{\overset{i=M}{\Sigma }}{\theta }_{\mathrm{Bi}}^t---\left(5\right)$>

式中：x、y、h分别为栅格行、列及高程数据；f_B为栅格内逐层土壤冰含量比重。

根据步骤2中，收集三江源内各点的土壤冰含量比重随季节的变化过程（此处选取长江源及黄河源区六个典型位置为例），如图9所示。

②利用土壤冻融状态冰含量比重分布曲线，分别核算相邻时刻计算深度内的冰含量比重，将其记做该深度内土壤蓄水容量的变化量：

dW=[f_B(t+1,x,y,h)-f_B(t,x,y,h)]H_soil        (6)

式中：dW为相邻时刻土壤蓄水容量的变化量；S_zm为参数；H_soil为计算深度。

③利用式(22)计算栅格型土壤涵蓄水能力指标：

W_M(t+1,x,y,h)=W_M(t,x,y,h)+dW       (7)

式中：W_M为栅格型土壤涵蓄水能力指标。

其具体累加分上、中、下三层，如图10所示。

对于正在发生解冻过程的栅格土壤：

首先，衡量土壤蓄水容量的变化量增加到上层土壤中，是否会超过上层土壤的最大涵蓄水能力，如果未超出，则仅上层土壤的涵蓄水能力得以更新，否则，将超出部分折算到中层土壤中，同时考察是否超出中层土壤的最大涵蓄水能力，如果未超出，则仅中层土壤的涵蓄水能力得以更新，否则，将超出部分折算到下层土壤中。

对于正在发生冻结过程的栅格土壤：

首先，衡量上层土壤的涵蓄水能力是否能够满足冻结量的需要，如果满足冻结需要，则仅上层土壤的涵蓄水能力得到更新，否则，上层土壤冻结，含蓄水能力为零，且冻结状态转移至中层，并考察中层的涵蓄水能力是否能够满足余下冻结量的需要，如果能满足冻结需要，则仅中层土壤的涵蓄水能力得到更新，否则，中层土壤冻结，含蓄水能力为零，冻结状态转移至下层。

步骤4：流域涵蓄水能力累积分布曲线

同步得出各栅格内涵蓄水能力后，结合流域分割，得到各子流域内栅格涵蓄水能力序列，根据栅格面积比建立各子流域离散化的涵蓄水能力累积分布曲线，其计算遵循以下步骤进行：

①结合土壤水热耦合计算模型，得出流域范围内各栅格的涵蓄水能力；

②结合流域水系构成、子流域划分范围，以及各子流域的拓扑关系，划分各子流域内所属栅格；

③统计各子流域内栅格涵蓄水能力，按照从小到大顺序排列，根据栅格面积比建立各子流域离散化的涵蓄水能力累积分布曲线，得到流域内的涵蓄水能力空间分布。

三江源流域共划分83个子流域，其中长江源流域有54个子流域，黄河源有29个子流域。

采用所提方法计算出两个个子流域的涵蓄水能力曲线见图11中的（a）、（b）所示。

结合涵蓄水能力指标，针对三江源流域内长江源及黄河源分别进行径流模拟，选取2005年资料率定，2006年汛期水文资料用于模型验证，以直门达、唐乃亥水文站流量过程为目标，分别对洪水的洪峰、确定性系数以及误差均方根等进行统计，得到的结果见表2。

表22005年模型率定结果

水文站洪峰相对误差/%确定性系数误差均方根/（m³·s^-1）唐乃亥-1.30.93144.3直门达0.980.9662.2

在率定过程中，通过根据冻融过程而动态调整流域涵蓄水能力后，流域出口流量模拟精度高。通过表2可知，流量计算过程确定性系数均大于0.9，洪峰误差在正负2%以内，误差均方根为60-150m³/s区间。

以相似的模型参数计算2006年长江源及黄河源流域径流过程，各选择两场洪水过程为例，对比其计算值精度，其计算结果如图12所示。由图中可见，2006年确定性系数均大于0.85，与实测过程拟合效果最好；由误差均方根显示，长江源直门达站的径流误差最大——101.4m³/s，其主要原因在于参与模拟土壤冻融活动规律分析的水文气象测点太少，且海拔均在4000m左右，无法完全反映流域冻融活动规律。

验证结果表明基于土壤季节性冻融规律建立起来的流域涵蓄水能力测算方法适用于寒区季节冻土区的产流计算。该方法不仅能计算流域各栅格点的涵蓄水能力，而且可以用于确定流域的涵蓄水能力曲线。实现将单点的涵蓄水能力确定方法成功应用于大尺度寒区水文模型中，提高产流模型的精度。