陆 峥，孙 景 博，何 源，李 宗 超，彭 书 艳，杨 晓 帆*

(1.北京师范大学地理科学学部，地表过程与资源生态国家重点实验室，北京 100875；2.南方海洋科学与工程广东省实验室(广州)，广州 广东 511458；3.中国环境监测总站，北京 100012)

0 引言

地表和大气间的动量、水热交换主导陆面过程及其与大气间的动力、热力相互作用[1-4]。区域尺度上地气间的相互作用能影响甚至改变气候变化响应机制，引发极端干旱和洪涝事件[5-7]。理解并准确表征地气间的动量、水热交换规律和传输过程对气候变化、农业、林业、水文和水资源管理等具有重要意义[8,9]。地下含水层系统物理参数(如渗透系数和孔隙度)的异质性对水分的交换和平衡起着重要控制作用，其中渗透系数是描述地下含水层系统空间异质性的重要参数，也是显著影响土壤和含水层水分、溶质运移模拟的参数之一[10]。研究表明，地表感热通量和潜热通量是关于浅层土壤水分的函数[11-13]；此外，地下水位对地表热通量和温度也有控制作用[14]，如农场灌溉和地下水开采可能会分别在短期和长期改变地下水位对地表热通量和温度的控制模式[15]。但在地下含水层系统渗透系数空间异质性对地表热通量的影响中，农业灌溉活动的作用尚不明确。

传统的分布式水文模型与地下水模型主要用于刻画水文循环现象及其物理过程，基于简化地下水模块的陆面过程模型则侧重表征地表土壤中及地气间的动量、水热交换过程，三者均难以精确描述地下水位和浅层土壤水分与地表热通量和温度之间的定量关系。因此，本研究采用地下水—陆面过程耦合模型(Parallel Flow-Common Land Model，ParFlow-CLM)[16]定量描述大气下边界层—地表径流—浅层土壤水分—深层地下水间的转换关系，以期揭示区域三维地下渗透系数场对地表热通量和温度的影响规律；同时考虑黑河流域中游密集农业灌区地表水灌溉和地下水回灌的作用，定量分析灌溉效应产生的浅层土壤水分对地表热通量的影响，以期为定量分析农业灌溉活动影响下区域渗透系数对地表热通量的影响提供参考，并为量化不同陆面过程对土壤和含水层系统异质性的响应提供借鉴。

1 研究区与数据

1.1 研究区概况

黑河流域地处西北干旱区，核心区面积约为13万km2，纵跨青海省、甘肃省和内蒙古自治区，是我国第二大内陆河流域，具有全球独特的、以水为桥梁和纽带的多层次自然景观，对我国中西部水文循环和生态环境影响巨大[17]。本文研究区为黑河流域中游地区(100.05°～100.69°E，38.73°～39.37°N)(图1)，该区水文地质结构较复杂，属中温带干旱型气候，下垫面类型主要为农田、草地和裸地，地带性植被多为温带小灌木、半灌木荒漠植被，地势东北、西南较高，中部为东南—西北走向、地形较平坦的狭长盆地[18-20]；区域年降雨量稀少且集中于6-10 月，平均潜在蒸(散)发量高于1 300 mm，主要水资源为南部祁连山区的冰雪融水和山区降水。特有的水土光热条件和平坦的地形使区内形成密集农业区，约900条主要灌渠构成错综复杂的人工灌溉渠系[21]。自西部大开发战略实施以来，农业灌溉和工业用水量大幅增长，地下水和地表水综合利用已成为该区农业和生活用水的主要形式[22]，河、渠、泉等地表水源灌溉及地下水回灌使地下渗透系数空间异质性对地表热通量和温度的影响作用更复杂[23]。因此，厘清该区域内地下含水层系统异质性对地表热通量和温度的影响，对我国战略水资源开发和生态环境保护意义重大[18,24]。

图1 黑河流域中游区域概况Fig.1 Sketch of the middle reaches of the Heihe River Basin

1.2 研究数据

研究数据包括：1)黑河流域大气驱动数据集[25]，可提供逐小时的气温、地表气压、水汽混合比、长短波辐射、二维风场和降水数据，原始分辨率为0.05°，采用文献[26]的方法降尺度至0.005°；2)地表覆盖类型和高程数据，分别源于黑河流域土地利用覆被数据集[27]和航天雷达地形测绘任务(SRTM)DEM数据集(http://srtm.csi.cgiar.org/srtmdata/)，原始分辨率分别为30 m和90 m，均重采样至0.005°；3)张掖灌溉渠系数据[28]，利用PriorityFlow[29]将DEM、河道、干渠数据处理为东西、南北向的坡度作为模型输入参数；4)土壤质地(粒径组成)空间分布数据集[30]；5)黑河流域30弧秒分辨率地表水及地下水灌溉量数据集[31]，利用黑河流域中游渠道流量测量数据集[32]中记录的灌溉时间，提取小时尺度的灌溉量；6)黑河流域地表过程综合观测网[33]中大满超级站(100.3722°E,38.8555°N)的涡动相关(EC)和大孔径闪烁仪(LAS)观测数据[34,35]。

2 模型与方法

2.1 ParFlow-CLM

ParFlow-CLM是一种开源、基于三维物理过程的地下水—陆面过程水文集成模型[36]，其将饱水带、包气带和地表水视为一个完整的水文连续体，能全面刻画不同尺度下基岩层到大气下边界层间的水文过程和水分、能量平衡[37]。对于地下部分，ParFlow基于三维地下水流动Richards方程求解饱水带—中间带—包气带的偏微分方程定解[38]，通过van Genuchten方程[39]描述各水动力变量的关系，并利用高效稳定的Newton-Krylov非线性求解器[40]进行矩阵求解；对于地表部分，则耦合CLM (Common Land Model,V2014修改版)[41]，利用二维地表水流动Saint-Venant方程的二阶动力波近似解模拟地表径流等过程，获得地表下渗、蒸(散)发、植被生理及雪动力等过程中的水热通量反馈到地下水与地表水的交界面，最终实现垂向的水分和能量平衡。此外，ParFlow-CLM采用的网格为结构化规则网格，可选正交网格或地形网格，模型详细的求解过程可参考文献[38,42,43]。

2.2 模型参数与情景设置

本文中模型参数设置(表1)为：水平分辨率为0.005°(约500 m)，垂直方向上网格结构分为多层，孔隙度、van Genuchten四参数(孔隙进气值参数、孔隙分布参数、残余含水量、饱和含水量)、曼宁粗糙度均设置为均质，其中孔隙度和van Genuchten四参数取自中国土壤数据集[44]，曼宁粗糙度修改自文献[45]，模拟区域的底部及四周均设置为零通量边界。使用2010年1月1日-2011年5月31日的驱动数据对模型进行起转率定(关闭地下水侧向流，计算步长为1 h)后，区域内的地下水储量变化范围稳定在5%之内，认为模型达到稳定[46,47]。最后，本文以起转率定获得的最终水位作为模型模拟的初始水位，并进行洪水测试和黑河核心干流径流量测试，结果表明本模型模拟地表水过程的精度较高。

表1 模型参数设置Table 1 Model parameter configurations

本文采用渗透系数表征土壤和含水层系统的异质性，并设计4种情景用于量化从土壤表层到基岩的渗透系数：1)Homogeneity情景，即模型的渗透系数恒定为0.01 m/h；2)GLHYMPS情景，仅在水平方向上设置异质性，渗透系数取自文献[48]的数据集，垂直方向上采用均质处理；3)CDHSP情景，在水平和垂直方向上均设置异质性，土壤渗透系数来源于中国土壤数据集[44,49]，包气带底部到基岩层渗透系数采用最底层土壤层的数值；4)MIX情景，表层土壤到包气带底部渗透系数的设置采用CDHSP情景，饱水带渗透系数的设置采用GLHYMPS情景。

为契合地表观测数据时段，模拟时段为2012年6月1日-10月31日，步长为1 h。ParFlow-CLM输出参数有地表感热通量(H)、潜热通量(LE)、上行长波辐射、下渗量、总蒸(散)发、土壤蒸发、植被蒸腾、不考虑升华的地表蒸发、雪水当量、地表温度(LST)、各土壤层温度(St)等。利用大满超级站涡动相关(EC)和大孔径闪烁仪(LAS)观测数据对感热通量、涡动相关观测数据对潜热通量、地表温度观测数据(采用文献[47]方法由大满超级站四分量仪器观测数据计算得到)对LST、地下土壤温(湿)度仪观测数据对St进行检验和情景对比分析。评价指标为均方根误差(RMSE)和Pearson相关系数(R)。

3 结果分析

3.1 地表热通量的验证与情景对比

如图2(彩图见封3)所示，地表感热通量和潜热通量模拟结果时序变化幅度均与观测数据相似，说明4种情景的模拟精度较高[15]，能在不同程度上捕捉到仪器观测数据的日内及季节变化趋势，较好地反映黑河流域中游水动力过程与能量传输过程。需要说明的是，相比潜热通量，感热通量与观测数据的一致性更高，R值达0.85以上(表2)，原因是潜热通量仅为点尺度观测，表征空间尺度的地表参量具有一定局限性。通过与观测数据对比可知，4种情景的模拟值均高估了感热通量、低估了潜热通量，这是由于模型中植被参数为固定值。从图2c 10月10-20日的数据可以看出，潜热通量模拟值在午后骤降，导致GLHYMPS情景的模拟结果与观测数据误差较小，其他3种情景误差较大，原因在于长期、持续的灌溉导致浅层土壤处于稳定的高饱和状态，打破了潜水埋深对潜热通量的控制，因而在饱和渗透系数情景(GLHYMPS)下模拟结果的误差较小。同时注意到，在4种情景的感热通量和潜热通量模拟中，Homogeneity情景的模拟精度最低，这充分说明在大尺度区域模拟中，考虑表层土壤—包气带—饱水带不同的水力性质并设置垂向和水平向非均质渗透系数是必要的，其能有效降低地表热通量的模拟误差。总体而言，Homogeneity情景的峰值较观测值有一定滞后性，GLHYMPS情景到达峰值后下降时间最早，而CDHSP和MIX情景的模拟值与观测值一致性很高，说明包气带的渗透系数对地表热通量具有决定性控制作用。

图2 ParFlow-CLM感热通量与潜热通量模拟值与大满超级站EC和LAS观测数据的验证与对比Fig.2 Validation and comparison of simulation values of ParFlow-CLM sensible heat flux and latent heat flux with EC and LAS observation data of Daman station

表2 感热通量(H)、潜热通量(LE)、地表温度(LST)和分层土壤温度模拟值与观测数据的验证及情景比较Table 2 Validation and scenario comparison of simulated and observed data for sensible heat flux,latent heat flux,land surface temperature and soil layer temperature

3.2 地表温度和分层土壤温度的验证与情景对比

由图3(彩图见封2)可知，4种情景模拟的地表温度总体上具有一定的契合度和相关性，R值均高于0.91，RMSE均小于3 K(表2)。其中，Homogeneity情景误差最大，MIX情景模拟值与观测数据误差最小、相关性最高，表明同时设置水平和垂直向异质性并考虑包气带—饱水带三维水力特征情景下地表温度的模拟精度最高。

图3 ParFlow-CLM地表温度模拟值与大满超级站观测数据的验证与对比Fig.3 Validation and comparison of simulation values of ParFlow-CLM land surface temperature with observation data of Daman station

利用大满超级站2、4、10、20、40、80、120、160(cm)多层土壤温度观测数据对4种情景模拟的不同土壤层温度进行日尺度比较(图4，彩图见封2)，可以看出，4种情景模拟的土壤温度廓线与观测值整体趋势一致，尤其是CDHSP和MIX情景下二者温度廓线形状最接近，但Homogeneity情景下土壤温度模拟值的垂向差异较大，尤其是在0.2～0.5 m土层深度范围内，GLHYMPS情景温度廓线的峰值出现更早且持续时间更久，与观测值差异最大。

图4 大满超级站观测数据与ParFlow-CLM土壤层温度模拟值对比Fig.4 Comparison of simulation values of ParFlow-CLM soil layer temperature with observation data of Daman station

进一步选取代表性较强的第一层、第四层和第七层土壤温度，与0.02 m、0.2 m和1.2 m土层深度观测值进行定量比较(图5)，可以看出，Homogeneity情景模拟值存在低估趋势，MIX情景模拟值与观测值最接近；结合表2中分层土壤温度模拟值与观测值的R和RMSE，可知MIX情景模拟的分层土壤温度精度最高，GLHYMPS情景精度最低，原因在于其使用的渗透系数难以反映真实的土壤水力传导情况。

3.3 空间分布特征

为探究土壤和含水层系统水力特征对地表热通量和温度空间分布的影响，以Homogeneity情景模拟结果为基准，进一步绘制GLHYMPS、CDHSP、MIX情景与Homogeneity情景模拟的地下水位埋深(WTD)、感热通量、潜热通量和地表温度空间差值图(图6，彩图见附录2)。1)从地下水位埋深的空间差值分布看，CDHSP情景模拟结果与Homogeneity情景较接近，MIX情景模拟结果在地形起伏较大的山区与Homogeneity情景差异较大，原因在于山区垂向渗透系数异质性较高，且存在地层分层现象。4种情景在平原区的整体差异较小。2)从地表热通量和温度的空间差值分布看，与Homogeneity情景相比，GLHYMPS情景的感热通量模拟值整体偏高，尤其是东北裸地和西南林草区，而潜热通量在上述两个区域基本为负值；CDHSP和GLHYMPS情景模拟结果与Homogeneity情景相反，在西部裸地及稀疏植被覆盖地区尤为明显，原因在于二者不同的土壤渗透系数参数化设置方案；相比Homogeneity情景，MIX情景稍显复杂，感热通量在东北裸地和西南林草区较高，在西北和西南裸地区偏低，地表温度和潜热通量值整体呈现偏高趋势，仅在西南黑河干、支流沿岸出现低值。综上，模型模拟的地下水埋深、地表热通量和地表温度的空间差值分布具有一定的敏感性，受下垫面、地形等环境因素影响而产生不同的陆面过程响应。整体上，4种情景在中部平坦的农田区差异较小，而在东北海拔较高的荒漠区差异明显。

图6 研究区平均地下水位埋深、感热通量、潜热通量和地表温度的空间差值对比Fig.6 Spatial difference comparison of averaged groundwater table depth,sensible heat flux,latent heat flux and land surface temperature

3.4 不同下垫面的影响

为深入探究下垫面因素通过土壤和含水层系统异质性对地表热通量的影响程度，本文绘制了模拟区域内裸地(荒漠及高山稀疏植被)和农田(农作物及农用地混合像元)两种下垫面对应的平均感热通量、潜热通量和地表温度箱型图(图7)。可以看出，整体上Homogeneity、GLHYMPS和MIX 3种情景结果较接近；CDHSP情景与其他情景模拟结果差异(标准差和中位数)较大，感热通量模拟值较高，潜热通量和地表温度模拟值偏低，这说明CDHSP情景下土壤和含水层系统的水力特性设置使得地表温度和热通量的模拟值域更广泛、分布更离散。横向比较，对于农田下垫面，4种情景下的平均热通量和地表温度模拟结果差异较小，这与前面分析结果一致；而对于裸地下垫面，Homogeneity和GLHYMPS情景的模拟结果中位数相近，MIX情景模拟的感热通量中位数比前两种情景高、潜热通量的中位数则较低。从模拟结果的值域范围看，MIX情景与Homogeneity情景最为接近。综上可知，裸地下垫面对地表热通量和温度的影响较大，MIX情景模拟结果的稳定性最好。

图7 裸地和农田下垫面的平均感热通量、潜热通量和地表温度箱型图Fig.7 Box plots of average sensible heat flux,latent heat flux and land surface temperature of bare land and cropland

4 结论

理解渗透系数的异质性对区域尺度水文集成模拟尤为重要。本文使用基于三维物理过程的地下水—陆面过程耦合模型ParFlow-CLM，依据不同土壤和含水层渗透系数情景设置，定量描述土壤和含水层系统异质性对地表热通量的影响。经过对地下水位埋深、地表感(潜)热通量、地表温度和分层土壤温度与观测数据的比较以及典型下垫面影响下时空分布特征的统计分析，主要结论如下：1)在黑河流域中游密集农业区，忽略水平向异质性会导致地表热通量模拟精度降低和不确定性增加；忽略垂直向异质性会导致地表温度模拟精度降低。同时，长期、持续的灌溉打破了潜水埋深对潜热通量的绝对控制，浅层水分成为主导土壤异质性影响地表能量分布及变化的因素之一。2)3种异质性情景与均值情景的地表热通量和温度差值空间分布对下垫面、地形等环境因素具有一定敏感性，下垫面因素对地表能量受土壤和含水层系统异质性的影响程度具有控制作用，尤其是对感热通量和地表温度的模拟影响较大。3)渗透系数的空间异质性在地形起伏较大的山区对地表热通量和温度的模拟影响远大于平原区。

本研究存在以下不足：1)模拟结果受输入数据(渗透系数数据集、驱动数据等)质量限制，需利用长期监测数据，进行不同尺度的对比和验证；同时，渗透系数数据集需利用野外勘察及辅助信息进行校核。2)模拟时长仍需继续扩展，最好能涵盖一个旱涝年周期，可深入探究渗透系数对地表热通量的年际及年内影响规律。3)黑河流域中游模拟区域渗透系数的异质性变化范围有限，仍需扩展到整个黑河流域甚至西北地区，以更全面地模拟土壤和含水层系统异质性对地下水—地表水交互的影响。

国家青藏高原科学数据中心提供黑河流域大气驱动数据集(2000-2018)、土地利用覆被数据集、30弧秒分辨率地表水和地下水灌溉量数据集、流域中游渠道流量测量数据集以及张掖灌溉渠系数据集、中国土壤数据集和祁连山综合观测网“黑河流域地表过程综合观测网数据集”，CGIAR-CSI提供SRTM DEM数据，此致谢忱！