土壤冻融过程中动态蓄水容量产流计算方法与流程

本发明涉及一种新的流域土壤冻融过程中动态蓄水容量产流计算方法，属于地球物理下水文分支技术领域。

背景技术：

流域产流作为水循环中最重要的一个环节，产流计算方法是流域水文模拟、水资源计算和洪水预报的重要理论基础。上世纪三十年代horton提出了超渗产流的概念，即当雨强大于下渗能力时产生地表径流，二十世纪六十年代河海大学赵人俊提出了蓄满产流的概念，即当降雨量满足土壤蓄水容量后产生径流(包括地表径流、壤中流和基流)，自此以后发展的水文模型都以超渗或蓄满作为产流计算的理论方法，近年来也有人认为流域内在不同时间会有两种产流机制共存的产流方式，提出混合产流计算方法。

多年和季节冻土广布于高寒地区，其水热传输贯穿于寒区流域的产流、入渗和蒸散发过程中，是寒区水文过程的核心环节。而常用的单一蓄满产流、单一超渗产流以及蓄超混合产流方式都不能涵盖高寒地区土壤冻融过程中的产流机理，新的流域土壤冻融过程中动态蓄水容量产流计算方法不仅可以提高高寒地区水文模拟和水资源计算精度，同时也大大丰富了寒区水文学的内涵。已有研究表明冻土不会阻止融雪/降雨水分下渗，融雪水/降雨在土壤冻结和不冻结情况下都会下渗到土壤中。关于冻土下渗产流问题，已有许多试验和理论研究，但是目前还没有能够用于实际且满足寒区流域径流预报的下渗强度和冻土透水程度的试验资料。在国内外广泛应用的水文模型(如vic模型、swat模型、新安江模型等)和陆面模式(如clm)中，至今没有考虑土壤冻融动态蓄水容量与产流的关系和相应的冻土产流计算方法。

技术实现要素：

本发明所要解决的技术问题是针对上述现有技术存在的问题，而提供一种提高高寒地区水文模拟和水资源计算精度的土壤冻融过程中动态蓄水容量产流计算方法。

为解决上述技术问题，本发明采用的技术方案是：

土壤冻融过程中动态蓄水容量产流计算方法，主要步骤如下：

(a)根据土壤温度分布计算流域网格内不同时间的冻土深度；

(b)分析冻土深度的逐日空间分布状况，得到流域不同网格内包气带达到的随时间分布的田间持水量w’m；

(c)找出流域不同网格内最大的田间持水量w’mm，并计算流域最大持水量wm，流域最大持水量wm为：

(d)根据得到的不同网格内最大的田间持水量w’mm和流域最大持水量wm，计算流域土壤冻融产流：

当p-e＞0，则产流，否则不产流，

产流量计算方法为：

若p-e+a＜w’mm则局部产流，有

若p-e+a≥w’mm，则全流域产流，有

r＝p-e-(wm-w0)(13)

式中，r为产流量；p为降水量；e是蒸发量；w0为流域初始土壤蓄水量；a、b为参数。

采用土壤水热耦合迁移模型计算步骤a中的冻土深度，其中土壤水热耦合迁移模型为：

其中，θu、θi分别为土壤中未冻水、冰的体积含量，t、z分别为时间及空间坐标，d(θu)、k(θu)分别为非饱和冻土水分扩散率与导水率，ρi、ρw分别为冰和水的密度，t为土壤温度，cvs、λ分别为土壤体积热容量、热导率，l为融化潜热，θmax(t)为相应土壤负温度(t)条件下可能的最大未冻水含量；

冻土深度由温度小于0位置处的空间坐标确定。

基于土壤水热耦合迁移模型推求的网格点冻土深度，通过方程(2)转换为土壤蓄水量：

w0＝h*θu(2)

其中，w0为土壤蓄水量，h为网格消融层深度。

所述步骤(b)包括流域单点动态蓄水容量计算和流域空间动态蓄水容量曲线获取：

流域单点动态蓄水容量计算：依据土壤冻融深度计算结果和土壤活动层上层土壤冻结和消融状态，获取流域单点逐日蓄水容量过程；

流域空间动态蓄水容量曲线：利用克里金插值法，对土壤特征参数比水容量cw(θu)、导水率k(θu)及扩散率d(θu)和土壤热特征参数体积比热容cvs和热导率λ进行空间插值分析，计算流域内每个网格内的土壤冻融深度和土壤活动层上层土壤冻结和消融状态，统计分析土壤冻融深度的逐日空间分布状况，绘制流域空间动态逐日蓄水容量曲线。

步骤(c)中，基于动态蓄水容量曲线的冻土产流计算方法，利用蓄满产流原理，计算得到流域入渗到土壤中的水量和产流两部分。

本发明基于多点冻土观测资料以及冻土下渗实验，分析土壤冻融过程中蓄水容量变化和下渗产流机理，进而推求流域土壤冻融动态蓄水容量曲线，提出一种新的流域土壤冻融过程中动态蓄水容量产流计算方法。

高寒地区是气候变化的敏感区，主要是因为随着气温升高，冰川、融雪和冻土对水文过程影响机理更加复杂，考虑到目前用于流域水文过程模拟的模型都没有考虑到土壤冻融过程动态蓄水容量对产流的影响，不能详尽刻画冻土产流过程、春汛洪水的模拟以及水资源演变计算，本发明把土壤水热耦合迁移运动与下渗产流有机结合，提出一种新的流域土壤冻融过程中动态蓄水容量产流计算方法，既可以计算土壤非冻结期产流，也可以计算土壤冻融过程中的产流。

本发明可以依据观测气温模拟逐日土壤冻结/消融、冻土深度和土壤温度，依据降雨观测计算融雪/降雨径流的逐日过程，提高了春季土壤消融期的径流模拟精度，为春汛防洪决策提供科学依据，同时也填补了现行国内外水文模型中冻土区产流计算的空白。

本发明采用以上技术方案与现有技术相比，具有以下技术效果：本发明利用土壤水热耦合迁移数值模拟土壤冻结和消融过程、土壤不同深度温度变化，研究冻土变化与包气带蓄水容量关系，提出土壤冻融动态土壤蓄水容量产流计算方法，发展动态土壤蓄水容量的产流模块，为高寒地区土壤冻融过程中的产流计算提供了一种新的方法，该流域土壤冻融过程中动态蓄水容量产流计算方法可以提高高寒地区水文模拟和水资源计算精度，推进寒区水文学的发展。

附图说明

图1为本发明流域土壤冻融过程中动态蓄水容量曲线计算方法流程图；

图2为本发明产流计算技术路线图；

图3为流域蓄水容量曲线概化图；

图4为基于流域蓄水容量曲线的产流量计算示意图；

图5为多组变动蓄水容量曲线示意图；

图6为本发明实施例中的冻土观测点冻土深度逐日变化过程图；

图7为本发明实施例中的冻土观测点逐日气温变化过程图；

图8为本发明实施例中土壤冻融深度模拟结果与实测结果对比图；

图9为本发明实施例中不同深度下(5cm)的土壤温度模拟结果与实测结果对比图；

图10为本发明实施例中不同深度下(20cm)的土壤温度模拟结果与实测结果对比图；

图11为本发明实施例中流域逐日流量模拟结果与实测结果对比图。

具体实施方式

下面结合附图和具体实施例，进一步阐明本发明。

本文以黄河源区某地区为例，采用本发明的方法对该地区土壤冻融过程中的产流进行计算。

具体包括如下步骤：

第一步：从中国气象数据网(http：//data.cma.cn)上下载研究区域气象站点逐日降水(雪)、日平均气温、日最高温度、日最低温度和0cm地温资料，利用克里金插值法，对下载资料进行空间插值分析，生成流域内每个网格的逐日资料系列。高寒地区多以融雪与降水为主要补给源的山区性流域，其中融雪量计算采用度日因子模型：：

m＝cm×(ti-tb)+ceer(3)

式中，m为日平均积雪消融量，cm为融雪的度日因子，ti为第i个网格融雪的日平均温度(℃)；tb为雪消融的临界温度，ce为雪的辐射系数，er为太阳短波辐射或净辐射。

第二步：利用土壤水热耦合迁移方程构建数值模拟模型，基于土壤特征参数、土壤热特征参数和黄河源区冻土深度和温度实验观测值，计算流域网格内不同深度的土壤温度和冻融深度随时间变化过程。根据计算土壤温度分布变化，识别不同时间低于0℃的土壤剖面分布，可以得到不同时间的土壤冻结厚度、冻结位置和冻结锋面，从而得到土壤活动层上层土壤冻结和消融状态。

1)在冻融过程中非饱和土水热耦合迁移过程的研究中，认为冻土中水分迁移规律与非饱和土壤水运动规律类似，可用含相变项的richards方程表示，以自变量为θ的richards方程：

式中，θu、θi分别为土壤中未冻水、冰的体积含量，t、z分别为时间及空间坐标(垂直向下为正)，d(θu)、k(θu)非饱和冻土水分扩散率与导水率，ρi、ρw为冰和水的密度。

该方程其特点是便于用数值模拟方法进行求解，适用于均质非饱和水分运动。

将相变潜热作为内热源的传导方程为：

式中，t为土壤温度，cvs、λ为土壤体积热容量、热导率，l为融化潜热。

上述(4)和(5)为土壤冻融过程中水热耦合迁移的两组基本方程，但是需要求解的是三个未知函数，即θu(z，t)、θi(z，t)和t(z，t)。因此还必须补充一个联系方程，即土壤中未冻水含水率θu与温度t的关系方程。在一定的负温度下，冻土中总含有部分未冻结水θu，并与负温度、压力等条件下处于动力平衡状态，在冻土研究中，当外界压力一定时，未冻水含量是温度的函数，可表示为冻土中水与热运动之间的联系：

θu≤θmax(t)(6)

式中，θmax(t)为相应土壤负温度(t)条件下可能的最大未冻水含量。

2)设定初始和边界条件。初始条件中的含水量分布θ0(z)和温度分布t(z)是已知的，土壤水分迁移的上边界条件是融雪/降雨入渗或者土壤蒸散发，下边界条件可为定水位和无限深度。土壤热流的边界条件在第一边界条件时，已知地表(z＝0)处温度随时间的变化过程t(t)及下边界处温度维持不变，为t(l)＝c。

3)土壤水分特征参数计算。与土壤水分迁移相关特征参数包括土壤水分特征曲线(土水势ψ或吸力s与土壤含水量关系)、比水容量cw(θu)、导水率k(θu)、扩散率d(θu)，各参数存在以下关系：

利用实验或者理论方法得到其中两个参数，其它参数即可计算得到。土壤水分特征曲线可以在野外或实验室内测定，也可以通过vg模型计算得到。

4)土壤热特征参数计算。土壤热特性参数包括体积比热容cvs和热导率λ，可以由试验测定，也可用半经验半理论公式计算。

5)冻土水热耦合迁移方程离散化要求。利用有限差分法求解，将计算区域离散化。由于土壤冻结锋面处水的相变引起释放大量潜热和冻土消融过程中吸收潜热，离散化时应取合适的距离步长和时间步长，在冻结锋面处距离步长取小一些。

6)不同深度土壤温度、未冻水含量和冰含量的计算。利用中心差分格式对冻土水热耦合迁移方程进行数值求解，可以计算出不同深度的土壤温度、未冻水含量和冰含量随时间变化过程。

7)土壤冻融深度的计算；根据计算土壤温度分布变化，识别不同时间低于0℃的土壤剖面分布，可以得到不同时间的土壤冻结厚度、冻结位置和冻结锋面，从而计算出流域网格内不同深度的土壤温度随时间变化过程，为动态蓄水容量计算提供数据。

第三步：统计分析土壤冻融深度的逐日空间分布状况，可得到不同网格下包气带达到田间持水量随时间的分布，绘制若干组不同时间下的流域动态蓄水容量曲线(图5)。

基于土壤水热耦合迁移模型推求的单个网格点冻土深度，可通过方程(2)转换为土壤蓄水容量。

流域上各个网格点包气带厚薄及土壤特性一般不相同，当全流域处于最干旱状态时，各处的包气带的缺水量不一定，即各处的包气带达到田间持水量不一样，其中最大的田间持水量为w’mm。将全流域面积看做1，以包气带田间持水量为纵坐标，小于等于某一田间持水量所占的流域面积比重为横坐标α，所得到的曲线(如图3)称为流域蓄水容量曲线：

曲线所包围的全部面积等于流域平均蓄水容量或最大持水量wm。

式中w’m为流域某处包气带达到的田间持水量，α值表示流域中≤w’m的流域面积所占的比重，b为流域蓄水容量曲线的方次，一般取值0.2～0.4，表征蓄水容量分布不均匀性的参数，b越大代表流域蓄水容量分布越不均匀。

第五步：利用蓄满产流原理，基于动态蓄水容量曲线的流域产流计算方法，计算得到流域入渗到土壤中的水量δw和产流两部分。如图4，若初始土壤含水

量为w。，则

当p-e＞0，则产流，否则不产流，产流量计算方法为：

若p-e+a＜w’mm则局部产流，有

δw＝p-e-r(12)

若p-e+a≥w’mm，则全流域产流，有

r＝p-e-(wm-w0)(13)

式中，w0～为流域初始土壤蓄水量(mm)；r～为产流量(mm)。

在本实施例中，选择黄河源区某区域作为研究区域，黄河源区一般是指河源至唐乃亥之间的区域，海拔高度在3000m以上，地处青藏高原的东北部，地理位置在95°50′～103°30′，32°20′～35°50′n之间。流域内属高原大陆性气候，主要为湿润半湿润气候区，多年平均气温为-4-5.2℃，年日照时数为2250-3131小时，平均风速3-4.5m/s。

为了验证本发明方法的实施，选择7月1日到翌年的6月30日为一周期，这个时间段能完全把黄河源区气象站的冻融期包括在内，实测数据包括该地区1997-2007年的冻土资料、地表温度和径流资料，图6和图7分别是该地区冻土观测点的冻土深度和气温的部分逐日过程，基于该实测资料对构建的土壤水热耦合迁移模型进行参数率定和验证，从而模拟区域单个网格点的土壤冻融深度和不同深度的土壤温度变化(图8-10)，从模拟结果来看，模型模拟效果较好，可为绘制流域空间动态逐日蓄水容量曲线提供数据。

利用本发明提出的土壤冻融动态土壤蓄水容量产流计算方法，发展动态土壤蓄水容量的产流模块，将模拟的地表径流量和实际观测值作对比，如图11所示，模拟的地表径流量与实测的流量比较接近，相对误差为4％，确定性系数为0.89，模拟精度较高，说明该发明提出的研究方法在高寒区具有较好的适用性。