原文信息

题目：基于冻土水文模拟的松花江流域水资源演变规律

作者：刘水清 周祖昊 刘佳嘉 李佳 谢新民 贾仰文 王浩   

期刊：《南水北调与水利科技（中英文）》23年1期

摘要 

为分析松花江流域水资源的演变规律，基于寒区水-热-氮素循环模型（the water and energy transfer processesand nitrogen cycle processes model in cold regions，WEP-N）和水资源评价方法，对径流发生突变的1998年前后（即1999—2018年和1956—1998年）进行比较，松花江流域年水资源总量减少217.0亿m3，减幅达到22.2%。其中，地表水资源量减少是水资源总量减少的主要组分，占水资源总量减少的比例为96.9%，不重复地下水资源减少量占3.1%。基于多因子归因分析方法分析可知，气候变化是松花江流域水资源减少的主要因素，对松花江流域全年水资源总量、地表水资源量、不重复地下水资源量减少的贡献率分别为81.6%、74.9%、286.6%，取用水的贡献率分别为18.4%、25.1%、－86.6%。从年内不同时期分析可知，非冻融期是全年水资源量减少的主要时期，占全年水资源总量减少的82.6%，冻融期占17.4%。和北方的海河流域、黄河流域相比，水资源减少幅度和主要影响因素各不相同，主要取决于气候变化和人类活动强度变化幅度的不同。与位于华北和西北的两大流域海河流域和黄河流域对比，气候变化对松花江流域水资源衰减的影响与海河流域相当，明显大于黄河流域，而人类活动对松花江流域水资源衰减的影响明显小于两大流域。

关键词  

松花江流域；水资源；演变规律；冻土水文模拟；归因分析

全球水资源衰减已受到广泛关注[1-3]。由于气候变化的影响和人类活动的日益加剧，寒区水循环正在发生巨大变化[4-8]。东北地区是我国主要的三大寒区之一，该区域广泛分布着季节性冻土和雪被，最北部的部分区域还存在永久性冻土[9-10]。松花江是我国的七大河流之一，是东北地区的主要河流。松花江流域作为我国重要的工农业基地，由于气候变化和人类活动影响，松花江流域的径流年际变化很大，经常出现连续丰水年和连续枯水年的情况[11-16]。松花江流域年内径流分配极不均匀，主要水文站点的全年径流主要集中在7−9月，可占到全年径流量的50%以上[11,17-18]。由于气候变化，松花江流域的冬季径流呈增加的趋势，夏季和秋季的径流呈减少的趋势[19]。

关于松花江流域径流变化规律的研究有很多。涂钢等[20]通过标准化的累积距平曲线的方法发现降水是影响松花江夏季径流的一个重要气象因子。李峰平[21]结合滑动t检验和降水-径流双累积曲线的方法发现1960−1974年属于人类活动影响很小的天然时期，1975−2009年属于大规模土地利用和水利工程建设等人类活动影响的时期，其中1975−1989年气候变化对松花江流域干流减少的贡献率为63%~65%，1990−2009年人类活动起主要作用[20]。王彦君等[22-23]基于累积量斜率变化率比较的方法分析了人类活动和降水变化对径流减少的贡献，发现在不考虑蒸散发影响时，人类活动对松花江干流径流减少的贡献率为57%~75%。史涵[24]基于Budyko假设的弹性系数法分析发现人类活动对松花江流域径流减少的贡献率高达96.1%，气候变化的影响仅为3.9%。现有的研究很少涉及流域地表水、地下水资源演变规律，更少关注到松花江作为一个寒冷地区的流域，其冻融期与非冻融期水资源究竟如何变化。在其他高强度人类活动寒区流域，现有的研究[25-26]对水资源的关注也很少。

因此，以东北松花江流域为例，应用寒区水-热氮素循环模型（the water and energy transfer processesand nitrogen cycle processes model in coldregions,WEP-N）[27]中的水文模块对流域水循环过程进行模拟，进而基于水资源评价方法对全年、冻融期和非冻融期的水资源进行评价，基于多因素归因分析方法，分析年内不同时期水资源总量、地表水资源量、不重复地下水资源量的演变规律及对气候变化和人类活动的响应，为松花江流域的水资源保护、利用以及高质量发展提供支撑。

1 研究区概况与数据来源

1.1 研究区概况

以中国东北部气候严寒的松花江流域作为典型，分析气候变化与人类取用水活动对水资源量的影响。松花江流域位于北纬41°42'~51°38'，东经76°52'~132°31'，海拔范围为50~2700m。松花江流域面积55.7万km2，跨黑龙江、吉林、辽宁和内蒙古4省（自治区）。松花江是黑龙江的主要支流，有南北两源，北源嫩江发源于大兴安岭支脉伊勒呼里山，南源发源于长白山天池。松花江流域存在季节性冻土，最大冻结深度大于2m，最长的冻结期为10月至次年7月，无冰川[28]。

1.2 数据来源

研究需要的数据主要包括：数字高程数据DEM，采用30″分辨率的SRTM30DEM数据，由美国太空总署（NASA）和国防部国家测绘局（NIMA）联合测定；气象数据，采用51个气象站1956−2018年逐日降水、气温、湿度、风速、日照时数作为基础输入，并采用反距离平方法将站点数据展布到子流域形心上作为模型驱动；土地利用数据，以中科院地理所提供的30m分辨率的土地利用数据为基础，共1990、2000、2005年3期；土壤及其特征数据，采用全国第二次土壤普查资料；1956−1999年松花江流域取用水数据与水资源量数据来源于《中国水资源及其开发利用调查评价》[29]，2000−2018年数据来源于松辽水资源公报（http://www.slwr.gov.cn/）；松花江流域的水库参数与调度资料来源于松辽水利网（http://www.slwr.gov.cn/）；WEP-N模型水文模块验证数据，江桥（1956−2018年）、扶余（1956−2000年、2006−2018年）和佳木斯站（1956−2018年）的实测流量数据来源于松辽流域水资源年鉴。

2 研究方法

2.1 WEP-N模型中的冻土水文模拟原理

WEP-N是寒区水-热-氮素循环模型，该模型已经应用于寒区流域并取得了很好的模拟效果[27]。WEP-N模型可以综合考虑气象、下垫面、人类取用水、水利水保工程等因素对水循环过程的影响，实现“自然-社会”二元水循环过程耦合模拟和分析，给出水循环要素时间和空间变化过程以及流域水循环通量。

WEP-N模型将流域划分为子流域套等高带单元，同时根据土地利用不均匀性又进行了下垫面细分，根据各下垫面面积所占比例分别进行产流计算[30-31]。每个等高带单元在垂直方向上分为7层，即植被冠层或建筑物截留层、地表洼地储留层、3层根系层、过渡层和地下水层。为了进行土壤水热耦合计算，将地面以下土壤分为11层。水循环模拟中，截留层蒸发、土壤层蒸发、水面蒸发和植被蒸腾等采用Penman-Monteith公式进行计算；积雪融化过程采用温度指数法；地表产流根据水量平衡计算；入渗根据降雨强度采用Green-Ampt模型或Richards方程计算；土壤热传递采用一维垂直流动的热流运动基本方程计算；土壤水运移采用一维垂直水分流动方程计算，土壤水运移由重力势、基质势、温度势、溶质势4种势能驱动；地下水的流动分山丘区与平原区，采用BOUSINESSQ方程进行浅层地下水二维数值计算，同时地下水与河流的水量交换受温度影响；坡地汇流根据各等高带的高程、坡度与Manning糙率系数采用一维运动波方程计算，由流域的最上游等高带依次计算至最下游等高带；河道的汇流计算采用一维运动波，河流运动过程中考虑河流冰的形成与融化过程，根据上下游汇流关系由上游子流域依次计算至下游子流域[30-31]。

一维垂直流动的热流运动基本方程[32]为

式中：z表示每一层土壤的深度，m；λs表示土壤的热传导率，W/(m•℃)；T表示每一层土壤的温度，℃；Cv表示土壤体积热容，J/(m3•℃)；t表示时间，s；Li表示融化潜热，3.35×105J/kg；ρi表示土壤中冰的密度，920kg/m3；θi表示体积含冰率，下标i表示与冰有关的参数。

土壤水分在土水势作用下发生运动。Richards方程计算土壤水的运移[33-34]为

式中：θ表示体积含水率；ρl表示水的密度，1000kg/m3；K(θ)表示非饱和水力传导度，m/s；H表示土水势，m。

总土水势H的计算公式为

式中：z表示重力势，m；hm表示基质势，m；hs表示溶质势，m；ht表示温度势，m。

2.2 水资源评价方法

水资源评价计算方法[35]为

式中：W表示水资源总量，亿m3；Rs表示地表水资源量，亿m3；Rg表示不重复的地下水资源量，亿m3。其中：Rs为来自降水的地表产流量、壤中流产流量、基流产流量之和；Rg为来自降水的地下水补给量减去基流产流量；地表产流量、壤中流产流量、基流产流量与地下水补给量可由WEP-N模型输出。

2.3 多因子归因分析方法

采用多因子归因分析方法进行归因分析，该方法采用某个因素所有情景的分解贡献量的平均值作为该因素在多因素综合影响下的影响贡献量，所有情景的权重是一样的。以平均权重为基础，可以使得计算所得各因素影响贡献量之和恰好等于所有因素综合影响下的变化量，而且计算结果唯一。根据影响水资源的气候（气温与降水）与取用水分别处于基准期和变化期，共设置4个模拟情景；使用寒区水-热-氮素循环模型WEP-N对每种情景分别进行模拟，求得各情景下某个水循环变量的多年平均值，采用以下公式进行贡献率分解[36]公式为

式中：∆Xj表示第j个因素的影响贡献量；αi, j表示第j个因素对应情景i的权重系数，由情景设置方案表推导得出，即如果使用变化期输入数据则等于1，使用基准期输入数据则等于−1；Si表示对应情景i的模拟结果（这里为方便就直接使用情景序号表示对应情景模拟值）；n表示所考虑的因素个数；A表示所有因素影响贡献量之和；βin表示在考虑n个影响因素的前提下第i个情景所有因素权重系数之和，即情景状态序列字符串中“1”的个数减去“0”的个数的差值。

各因素对水循环要素改变的贡献率计算公式为

3 结果与讨论

3.1 WEP-N模型验证

WEP-N模型在黑顶子河小流域结合冻融期（2017−2018年、2018−2019年）土壤水热耦合实验和河道流量监测实验对分层土壤含水量和温度以及河道流量进行了验证[37]。基于WEP-N模型，对松花江流域干流上主要水文站点1956−2018年的流量过程进行了验证。其中，1956−1990年为率定期，1991−2018年为验证期。各站流量过程模拟结果和效率系数见图1和表1。率定期与验证期松花江主要干流水文站月流量过程Nash-Sutcliffe效率系数(ENS)均在0.7以上，相对误差(ER)均在5%以内。结果表明，WEP-N模型总体上能够较好地描述松花江流域的径流过程。

3.2 基于WEP-N模型的松花江流域水资源评价

松花江流域水资源评价结果对比见图2。与《中国水资源及其开发利用调查评价》和松辽水资源公报中的水资源数据相比，松花江流域1956−2018年地表水资源量、不重复地下水资源量、水资源总量年过程评价的相对误差分别为−5.0%、7.2%、−3.3%，相关系数分别为0.88、0.62、0.88，模型评价效果较好。

3.3 松花江流域水资源演变规律分析

采用Mann-Kendall检验方法分析可知，松花江流域1956−2018年的水资源总量整体呈减少的趋势，变化率为−0.33亿m3/a（图3）。采用Pettitt检验方法可知，松花江流域1956−2018年水资源总量在0.05显著性水平下具有突变，突变年份为1998年。因此，以1998年为节点，对1998年前后水资源量变化及影响因素进行分析，即将基准期和变化期分别设置为1956−1998年和1999−2018年。

气候变化与人类活动是影响流域水资源的两类主要因素。在松花江流域，下垫面类型发生了一定变化，但是各下垫面类型的变化量占松花江流域总面积的比例均小于6%[37]，因此在分析过程中忽略下垫面变化对水资源的影响，人类活动主要考虑取用水的影响。突变期（1999−2018年）相对基准期（1956−1998年）水资源量减少，采用多因子归因分析方法对气候变化与人类活动的贡献率进行分析。情景设置见表2，由于考虑两类因素，因此设置4个情景，主要影响因素多年平均值见表3。采用WEP-N模型对各情景分别进行模拟，并基于水资源动态评价方法得到松花江流域的水资源总量、地表水资源量、不重复地下水资源量，然后采用公式（5）~（8）进行贡献率计算，松花江流域突变期1999−2018年相对基准期1956−1998年水资源变化的归因分析结果见表4。由表4可知，基准期松花江流域全年水资源总量、地表水资源量、不重复地下水资源量分别为977.4亿、837.3亿、140.1亿m3，变化期分别为760.4亿、627.0亿、133.4亿m3。由于气候和取用水变化，松花江流域全年水资源总量的减少量为217.0亿m3，变化率为−22.2%；其中：地表水资源的减少量为210.3亿m3，变化率为−25.1%；不重复地下水资源的减少量为6.7亿m3，变化率为−4.8%。松花江流域地表水资源是水资源总量减少的主要组成部分，占全年水资源减少量的96.9%，不重复地下水资源量的减少量占3.1%。气候变化是松花江流域水资源减少的主要因素。基于多因子归因分析方法计算可知，气候变化对松花江流域全年水资源总量、地表水资源量、不重复地下水资源量减少的贡献率分别为81.6%、74.9%、286.6%，取用水的贡献率分别为18.4%、25.1%、−186.6%。

情景S1中气温、降水、取用水均处于基准期（1956−1998年），通过改变单一因素，使其处于突变期（1999−2018年），然后基于WEP-N模型与水资源动态评价方法得到单一因素影响下的水资源分项（表5）。由表5可知，气温升高1.3℃、降水减少140.2亿m3、取用水增多85.6亿m3分别导致水资源总量减少56.8亿、103.5亿、23.2亿m3。单位气温升高、单位降水减少、单位取用水增多导致水资源总量分别减少43.70亿、0.74亿、0.27亿m3。在气候与取用水影响下，水资源演变结果大体是合理的。从水资源分项变化可知，各因素对地表水资源的影响大，对不重复地下水资源的影响小。取用水增多造成不重复地下水资源量增多，主要是因为取用水特别是地下取用水造成地下水位下降进而导致地下水补给量增加，同时由于地下水位下降导致基流减少，最终造成不重复地下水资源量增加。因此，结果是合理的。

由表6可知：在冻融期，由于气候和取用水变化，松花江流域水资源总量减少37.7m3，其中地表水资源量的减少量为46.0亿m3，不重复地下水资源量增加量为8.3亿m3。冻融期水资源总量、地表水资源量、不重复地下水资源量的减少量分别占全年减少量的17.4%、21.9%、−121.2%。变化期1999−2018年与基准期1956−1998年相比，虽然冻融期降水增多，但是由于气温升高导致蒸散发增大，同时取用水增加，所以冻融期松花江流域的水资源总量与地表水资源量整体呈减少的趋势。对于不重复地下水资源而言，由于土壤冻融时间变短、地下取水量增加，同时降水增多，共同导致地下水补给量增大。

在非冻融期，气温升高、降水减少，同时取用水量增加，导致松花江流域水资源总量减少179.3亿m3，其中，地表水资源量减少164.2亿m3，不重复地下水资源量减少15.1亿m3。非冻融期是全年水资源变化的主要时期，其水资源总量、地表水资源量、不重复地下水资源量的减少量分别占全年减少量的82.6%、78.1%、221.2%。

3.4 松花江流域与北方其他流域水资源演变对比分析

受气候变化与人类活动影响，北方大部分流域的水资源总量呈减少的趋势。海河是中国华北地区的最大水系。据曹建廷等[38]关于海河水资源总量演变的分析[39]，1980−2010年与1956−1979年相比，海河流域的水资源总量减少了37.2mm，减幅为28.2%[39]，而与1956−1998年相比，1999−2018年松花江流域的水资源总量减少了40.2mm，减幅为22.2%，小于海河流域的减少幅度。由分析可知，受气候变化影响，松花江流域（1999−2018年与1956−1998年相比）与海河流域（1980−2010年与1956−1979年相比）水资源总量分别减少32.8、26.7mm，松花江流域与海河流域水资源总量的减幅分别为18.1%和20.6%，松花江流域水资源总量的减幅略微小于海河流域。在人类活动影响下松花江流域与海河流域水资源总量分别减少7.4、10.5mm，水资源总量减幅分别为4.1%和7.6%，松花江流域的减幅明显小于海河流域。这主要是由于海河流域取用水量占水资源总量的128%（2000−2010年海河流域水资源公报http://www.hwcc.gov.cn/hwcc/static/szygb/hhlyszygb.htm），远远大于松花江流域的26.1%[39]。另外，在对比两个流域水资源变化的过程中，两个流域采用的时间系列不一致。2010−2018年海河流域依然存在气温升高、降水减少与取用水增多的趋势[40-42]，因此，若对比海河流域1999−2018年相对1956−1998年，在气候变化与取用水影响下海河流域水资源总量的减幅会更大。

黄河是中国北方重要的淡水资源，是中国水资源短缺非常严重的流域。黄河流域地表水资源量占水资源总量的84.2%[43]，据周祖昊等[44]关于黄河流域地表水资源演变的分析，与1956−1979年相比，1956−2016年黄河流域地表水资源量减少了15.6mm，减幅为20.2%。松花江流域地表水资源量占水资源总量的比例略微大于黄河流域，为85.4%，与1956−1998年相比，1999−2018年松花江流域的地表水资源量减少了39mm，减幅为25.1%，减少幅度明显大于黄河流域。受气候变化影响，松花江流域（1999−2018与1956−1998相比）与黄河流域的地表水资源量（1956−2016年与1956−1979年相比）分别减少29.2、3.8mm，减幅分别为18.8%和4.9%。气候变化对松花江流域水资源的影响远大于黄河流域，是因为变化期相对基准期，松花江流域气温增加1.3℃，远大于黄河流域的0.45℃，同时松花江流域降水的减幅4.8%又明显大于黄河流域1.7%。在取用水影响下，松花江流域与黄河流域地表水资源量分别减少9.8、7.9mm，减幅分别为6.3%和10.1%。取用水对松花江流域地表水资源的影响明显小于黄河流域，主要是因为松花江流域的取用水量只占水资源总量的26.1%，远小于黄河流域取用水量在水资源总量中的占比50.0%[43-44]。黄河还有一个特殊问题，水土保持规模巨大，对水循环影响大。结果显示，在下垫面变化影响下，黄河流域的地表水资源量减少了3.9mm[44]。另外，黄河流域的基准期和变化期两个系列有24年的重叠，如果对比黄河流域1980−2016年相对于1956−1979年或1999−2018年相对于1956−1998年，在气候与取用水影响下地表水资源的减幅会更大。

由以上对比可知：气候变化对东北寒区典型流域松花江流域水资源衰减的影响与华北最大的水系海河流域相当，明显大于黄河流域；人类活动对松花江流域水资源衰减的影响明显小于海河和黄河流域。

4 结论

基于WEP-N模型对松花江流域水循环过程进行模拟，松花江主要干流水文站江桥、扶余、佳木斯月流量过程Nash-Sutcliffe效率系数均在0.7以上，相对误差均在5%以内。

根据突变性分析，1998年为径流衰减的突变点。1999−2018年和1956−1998年比较，由于气候变化和人类活动影响，松花江流域水资源总量减少217.0亿m3，减幅达到22.2%。其中，地表水资源量减少是水资源总量减少的主要组分，不重复地下水资源的减少量远远小于地表水资源的减少量。从年内不同时期分析可知，非冻融期是松花江流域水资源减少的主要时期。基于多因子归因分析方法计算可知，气候变化是导致松花江流域水资源减少的主要因素。

与位于华北和西北的两大流域海河流域和黄河流域对比，气候变化对松花江流域水资源衰减的影响与海河流域相当，明显大于黄河流域，而人类活动对松花江流域水资源衰减的影响明显小于两大流域。