基于水热耦合的青藏高原分布式水文模型——I. “积雪-土壤-砂砾石层”连续体水热耦合模拟

周祖昊,刘扬李,2,李玉庆,王鹏翔,,王 康,李 佳,朱熠明,,刘佳嘉,王富强

(1. 中国水利水电科学研究院流域水循环模拟与调控国家重点实验室,北京 100038；2. 中国电建集团贵阳勘测设计研究院有限公司,贵州贵阳 550081；3. 西藏农牧学院水利土木工程学院,西藏林芝 860000；4. 武汉大学水资源与水电工程科学国家重点实验室,湖北武汉 430072；5. 水利部南水北调规划设计管理局,北京 100038；6. 华北水利水电大学水利学院,河南郑州 450046)

目前,基于寒区水热耦合的特点,分布式水文模型已取得一定的进展[1-2]。国际上常用的分布式水文模型都已对寒区的积雪和冻土过程进行了一定程度的拓展表达,部分模型中已加入冻土、融雪模块来研究寒区水文循环,如CRHM(Cold Regions Hydrological model)、SWAT(Soil and Water Assessment Tool)、GEOtop (Geography and Topography model)、DWHC(Distributed Water-Heat Coupled model)、VIC (Variable Infiltration Capacity model)模型等[3-4]。当前较成熟的冻土水热耦合模型是基于非饱和土体水分迁移,用土壤未冻水含量随温度变化的函数来表征冻结机理的水动力学模型。其中以Harlan[5]提出的水热耦合模型应用最为广泛；Taylor和Luthin[6]在其基础上用未冻水含量结合土壤水分特征曲线计算土水势,确定水分迁移的驱动力；Horiguchi[7]又加入了温度梯度共同作用于水分迁移；美国农业部Flerchinger和Saxton[8]综合考虑气相、水流迁移引起的热量变化以及边界处的蒸发、辐射等,建立了垂直一维水热耦合模型,较全面地反映了影响冻土水热变化因素；雷志栋等[9]对各向同性土柱的一维冻结问题简化了模拟过程,忽略了水汽迁移、热量对流作用,提出了冻土非饱和土壤水分迁移方程,对地下水浅埋条件下的土壤冻融过程进行了模拟,模型模拟效果较好。水动力学模型参数简便易得,是较为成熟的水热耦合模型。

青藏高原是长江、黄河、澜沧江、雅鲁藏布江等多个大江大河的发源地。作为世界面积最大、海拔最高的高海拔寒区,平均海拔4 376 m[10],其广泛分布的积雪和冻土影响了该地区土壤介质的储水、导水、产水和热传导过程[11]。与中国其他地区相比,青藏高原土壤层较薄,土壤层下存在较厚的砂砾石层。与土壤相比,砂砾石具有不同的水、热性质[12],孔隙度、密度也有所不同,对土壤含水量和导水率有较大影响[13-14]。对青藏高原的水热运移模拟,部分陆面模式虽然考虑了砂砾石影响,通过砂砾石占比对土壤热传导率、土壤质地和地形等敏感参数进行调整,模拟效果会有一定提升[15-17],但由于模型将土壤概化为一维均质结构,仅通过参数调整难以体现该地区这种上下分层明显的地质结构对水热运移的影响,模拟结果均有一定偏差。因此,在研究青藏高原土壤冻融与水循环过程时,不仅要考虑土壤层内的水热迁移,对于土壤层下伏较厚的砂砾石层也应一并考虑。与此同时,积雪较好的隔热性和对太阳短波的反射性,也会影响土壤水热迁移。在水热耦合研究过程中,要把土壤和积雪作为一个复杂的、相互作用的联合体来统一考虑[18]。对积雪的研究,根据其模拟的复杂程度与否,大致可分为以下两类:一是利用强迫-恢复(force-restore)法[19-20],通过相对简单的单层积雪模型,区分积雪和土壤的热力学性质,模拟积雪-土壤复合层的温度变化,以度日因子法计算融雪量[21]；二是以物理机制为基础的复杂模型,此类模型考虑了积雪内部的质量守衡和能量平衡以及积雪表面与大气环境条件的相互作用,将积雪细分为多层,对积雪内部相变、水分运移等进行了精细描述[22-23]。但后一类模型太过复杂,影响了其推广与应用。

本文针对青藏高原的气候和地质特点,在土壤冻结融化期开展了野外土壤和砂砾石层水热监测试验,以WEP-COR模型[1]的土壤水热耦合模块为基础,构建基于“积雪-土壤-砂砾石层”连续体的水热耦合模型,并采用试验观测结果对模型进行验证。

1 研究区与资料

1.1 研究区选取

选取雅鲁藏布江的一级支流尼洋河流域作为本文研究的典型区域。尼洋河位于西藏东南林芝地区雅鲁藏布江中下游左岸,北纬29°28′—30°38′,东经92°10′—94°35′,源于西藏自治区米拉山西侧的错木梁拉,河源海拔约5 000 m,在林芝市巴宜区汇入雅鲁藏布江,河口海拔约2 920 m,落差2 080 m。尼洋河干流长309 km,流域平均坡降为0.73%;流域面积1.75万km2,在雅鲁藏布江众支流中排行第四。尼洋河流域地处中纬度地带,年均降水量1 295.1 mm,年均气温7.6 ℃左右[24]。

1.2 试验布设

1.2.1 土壤厚度和成分调查与取样分析

根据尼洋河流域实地考察(图1)和部分工程设计的地质勘测资料,该流域具有土壤层厚度较薄、土壤下方的砂砾石层厚度较大的特点。尼洋河流域的砂砾石层主要包括砾石层和卵石层,砾石层主要为圆状砾石,含卵石颗粒,砾石质量百分比约50%～65%,黏粒质量百分比5%～10%,孔隙间填充中细砂；卵石层孔隙间填充圆砾和中细砂,局部含大漂石、孤石。

针对尼洋河流域这一地质特点,为研究流域内土壤层从山脚到山顶的厚度变化规律,从尼洋河流域的源头到河口选取了海拔不同的24个取样点。其中点1—点16沿河道,点17—点24从山脚到山顶测量的土壤厚度(图2)。对24个取样点的土壤厚度及组分进行了测量和分析,尼洋河流域土壤主要为砂质壤土和壤质砂土,平均砂粒质量百分比为55.89%,粉砂粒质量百分比为31.2%,黏粒质量百分比为12.91%,厚度随着海拔升高而减小,平均厚度为53.7 cm。

图1 土壤-砂砾石层结构Fig.1Soil-sand gravel structure

图2 尼洋河流域土壤取样点、试验点及气象站位置Fig.2Soil sampling points and meteorological stations

1.2.2 冻结融化期“积雪-土壤-砂砾石层”监测试验

本次研究在尼洋河流域下游色季拉山的山腰上选择一处开展季节冻结融化期水热耦合过程监测试验。试验点(94°21′45″E,29°27′12″N)海拔为4 607 m,试验期为2016年11月—2017年4月,试验期间试验点降水8.3 mm,平均气温-3.4 ℃,最高6.3 ℃,最低-11.7 ℃。冻结融化期,试验点位置见图2。试验人员在试验点开挖了一个深度为1.6 m的试验坑,在地面以下垂直深度上每隔10 cm安装自动测定液态含水率的TDR传感器、测量温度的PT100传感器和测量基质势的TensionMark传感器,对下垫面冻结和融化过程中水、热、势能过程进行自动监测,并根据PT100传感器和TensionMark传感器测量结果推算土壤冻结深度,同时于开始降雪后间隔1～2周根据降雪情况通过米尺实地测量下垫面覆盖的雪层厚度。试验前,采用核磁共振方法对高原季节性冻融条件下水、热过程监测的仪器设备进行率定。仪器完成安装后进行原状土回填,保证回填土壤与实际土壤的参数一致。

1.3 气象资料

气象条件是水热耦合模拟的边界条件。由于流域内气象站点稀少,且高程差别较大,试验点的逐日气温采用考虑高程修正的RDS(Reversed Distance Squared)方法由气象站数据插值得到,降水数据采用平面插值、结合西藏自治区水资源公报[25]中降水等值线修正得到,逐日相对湿度、日照时速和风速采用RDS方法插值得到。本文用到的基本气象数据主要为尼洋河流域境内的林芝站和周边嘉黎站(图2)的降水、气温、相对湿度、日照时数和风速。气象数据来源于中国气象数据网(http:∥data.cma.cn),使用数据时间序列为2016—2017年,时间尺度为日。

2 “积雪-土壤-砂砾石层”连续体水热耦合模型构建

本文在WEP-COR模型中冻土水热耦合计算原理基础上,结合青藏高原特点加以改进,构建包含“积雪-土壤-砂砾石层”连续体的高寒山区水热耦合模型。

2.1 WEP-COR中的冻土水热耦合计算原理

WEP-COR模型在保留原WEP-L模型关于地表与大气间的能量交换计算的基础上,系统添加了土壤层间的水热通量和相关水热参数的计算过程,并且将活动土壤层细分为11层,见图3。其中,Rs为地表产流;Ri为第i层土壤的横向流或壤中流,Ri与坡度和土壤含水率有关;E1为土壤蒸发量;Er为植被蒸腾量;Qi为第i层土壤的重力排水;P为降水;Ta为大气温度;Ti为第i层土壤的温度;Gi为第i层土壤与相邻土壤层间由温度差异引起的热通量。水热耦合模拟采用显式差分进行数值迭代计算,模型模拟的时间步长为1 d。

图3 WEP-COR模型分层水热通量计算结构示意Fig.3WEP-COR model water and heat flux simulation structure

2.1.1 土壤水热运移方程

模型假设土壤冻融时只有液态水发生运移,土壤水分的运移主要受重力势、基质势和温度势的影响,因此添加的土壤一维垂直水分流方程可以写成[26]:

(1)

根据能量平衡原理,冻融系统中每一层土壤的能量变化都用于系统内的土壤温变和水分相变,温度势是水分相变的驱动力,而大气与表层土壤的温差则是热传导的源动力。假设土壤各向均质同性,并忽略土壤中的水汽迁移,添加的一维垂直流动的热流运动基本方程可变形为[27]

(2)

式中:θl、θi分别为土壤中液态水、冰的体积百分比,cm3/cm3；T为土壤温度,℃；t、z分别为时间、空间坐标(垂直向下为正)；D(θl)、K(θl)分别为非饱和土壤水分扩散率、导水率,cm2/s；Cv为土壤体积热容量,J/(m3·℃)；λ为土壤导热系数,W/(m·℃)；ρi、ρl分别为冰、水密度,kg/m3；Lf为融化潜热,Lf=3.35×106J/kg。

2.1.2 土壤水热联系方程

土壤冻融过程中,冻土水热运动间的联系主要表现在未冻水的含水率与土壤负温的动态平衡中[19-20]:

θl=θm(T)

(3)

式中:θm(T)为土壤负温对应的最大未冻水含水率。

2.2 “积雪-土壤-砂砾石层”连续体水热耦合模型原理

2.2.1 模型结构

积雪作为一种良好的隔热体,具有热容量大、导热系数较小等特点,对大气和土壤间的热传递具有明显的削弱作用。试验点所处位置积雪持续时间较长、厚度稳定,对土壤水热的影响不容忽视。根据尼洋河流域实地考察(图1)和基于26个点的测量结果建立的土壤厚度—高程关系图(图4),发现尼洋河流域土壤厚度整体较薄,在土壤层下方存在厚度较大、混合砂砾石和少量土壤的砂砾石层。而且土壤层厚度在河谷区基本稳定,在山坡上随着高程增加呈现变薄的趋势,在山顶附近厚度又基本保持不变。在试验点处土壤层为砂质壤土,砂粒质量百分比为61.92%,粉砂粒质量百分比为26.2%,黏粒质量百分比为11.88%,土壤层厚度为40 cm；土壤下砂砾石层砾石质量百分比为53.7%,孔隙间填充中细砂。土壤与砂砾石层水分特征参数见表1。

图4 尼洋河流域土壤厚度—海拔关系Fig.4Relationship between soil thickness and elevation in Niyang River basin

表1 尼洋河流域土壤水动力参数

根据尼洋河流域气候和地质特征,本文将水热耦合模拟的对象定义为“积雪-土壤-砂砾石层”连续体(图5)。模型计算结构上,在WEP-COR模型11层土壤层分层的基础上,于土壤层上方加入一层积雪层,同时将原有的11层土壤层细分为土壤和砂砾石层2类介质,计算结构如图6所示。第0层为积雪,根据等高带在计算单元中所处位置设置上部i层为土壤层,下部12-i层为砂砾石层。由于距离地表越近,水、热变化越快,故第1、2层厚度设为10 cm,第3～11层厚度设为20 cm。值得注意的是,土壤-砂砾石层各层的厚度需根据所在位置含水层的实际厚度进行修正,也就是说,当含水层总厚度小于11层土壤-砂砾石层厚度时,须根据含水层实际厚度确定土壤-砂砾石层的计算层数和最后一层的厚度。

图5 “积雪-土壤-砂砾石层”垂向分层概化Fig.5Vertical layering of “snow-soil-sand gravel layer”

图6 分层水热通量计算结构示意Fig.6Layered calculation structure of water and heat flux

2.2.2 计算原理

采用度日因子法[28]计算积雪融雪过程,具体如下:

M=k(Ta-T0)

(4)

式中:M为当日融雪当量,mm；k为积雪融化的度日因子,mm/(℃·d)；T0为积雪融化临界温度,℃,本文取-1 ℃。

对于积雪的度日因子,一般情况下在1～7 mm/(℃·d),本次研究根据下垫面情况取为4 mm/(℃·d)。由于青藏高原地区昼夜温差较其他地区偏大,本次研究的雨雪临界气温设为2 ℃。

积雪热量平衡公式[29]:

G=Le+Rn+Q1+Q2

(5)

式中:G为积雪层的能量变化,J/(m2·d)；Le为积雪的蒸发、融化和升华潜热通量,J/(m2·d)；Rn为积雪层表面的净辐射量,J/(m2·d),由到达积雪表面的净长波辐射与净短波辐射计算而得；Q1为积雪与大气间的热量交换量,J/(m2·d),采用强迫-恢复法计算；Q2为积雪与土壤间的热量交换量,J/(m2·d)。

积雪和土壤间的热量交换量采用以下公式[30-31]:

(6)

式中:ZS为积雪厚度,m；ZC为第一层土壤厚度,m；λS为冰雪层导热系数,W/(m·℃)；λC为土壤层导热系数,W/(m·℃)；RC为土壤层与积雪层的接触热阻,(m2·℃)/W；TC为第一层土壤层温度,℃；TS为积雪层温度,℃。

土壤层之间的水热运移方程和联系方程采用式(1)和式(2),砂砾石层之间以及砂砾石层和土壤层之间的水热通量计算采用与土壤层相同的计算原理,但砂砾石层的各项水热参数与土壤不同。

模型的数值模拟仍采用显式差分进行数值迭代计算,模拟步长为1 d。每个时间步长内模型先根据初始条件从积雪层到底层砂砾石层逐层计算热传导、温度和水分相变,然后以该层温度为判定条件进行迭代计算直至收敛(步骤1)；热量计算闭合收敛后,再根据水量平衡计算各层水分运移量,并修正各土壤和砂砾石层(积雪层不考虑)的含水率和含冰率(步骤2)；用液态含水率判定是否收敛,若不收敛则返回步骤1进行热量计算,直至两项迭代计算均闭合收敛后,“积雪-土壤-砂砾石层”连续体水热耦合计算完成。

2.2.3 模型边界确定

模型上、下边界的确定方式与WEP-COR模型类似。“积雪-土壤-砂砾石层”连续体的上边界为大气,但与上边界直接接触的不一定只是土壤。传入连续体的热量,在非积雪区由地表附近的大气和表层土壤的温度差以及表层土壤的水热参数决定,在积雪区由地表附近的大气和积雪的温度差以及积雪的水热参数决定,均采用强迫-恢复法计算。连续体的下边界为过渡层或者地下水层,将下边界温度假设为一恒定的温度作为底层边界温度。

2.2.4 积雪、土壤和砂砾石层水热参数确定

(7)

(8)

Cs=2.09×103ρs

(9)

式中:ρs为积雪密度,kg/m3；λs为积雪的导热系数,W/(m·℃)；Cs为积雪的体积热容量,J/(m3·℃)。

土壤、砂砾石层主要的水热参数也包括体积热容量、导热系数和土壤导水率等,各参数计算公式如下。

体积热容[35]:

Cv=(1-θs)Csg+θlCl+θiCi

(10)

导热系数计算参考了IBIS 模型,具体如下[36]:

λ=λst(56θl+224θi)
λst=1.500ωrock+0.300ωsand+0.265ωsilt+0.250ωclat

(11)

导水率[29]:

(12)

式中:Csg、Cl和Ci分别为土壤(砂砾石层)、水、冰的体积热容量,在0℃时土壤和砂砾石层分别为1.93×103J/(m3·℃)、3.1×103J/(m3·℃),水和冰分别为4.213×103J/(m3·℃)、1.94×103J/(m3·℃)；λ、λst为土壤或砂砾石层实际的导热系数、干燥状态下的导热系数,W/(m·℃)；ωrock、ωsand、ωsilt和ωclat分别为岩石(砾石和卵石)、砂粒、粉粒和黏粒的体积比；KS为土壤或砂砾石层经温度修正后的饱和含水率,cm/s。

参考陈仁升等的DWHC模型[37],不同温度条件下土壤或砂砾石层KS计算方法如下:

(13)

式中:K0为常温下的饱和导水率,cm/s；k0为冻结条件下最小的导水率,cm/s；Tf为土壤或砂砾石层最小导水率对应的临界温度,℃。本文考虑到土壤和砂砾石层水动力学特性的区别,对于土壤,k0按照原公式取值为0 cm/s；对于砂砾石层,由于孔隙较大,k0取值为大于0的值。

3 模型验证与结果分析

3.1 评价方法

本文采用决定系数(R2)及平均根方差(ERMS)评价模型模拟效果。R2和ERMS的计算公式如下:

(14)

(15)

3.2 土壤和砂砾石层水热耦合模拟

试验点冻结融化期间空气温度、模拟雪厚、实测雪厚如图7所示,试验点积雪从2016年12月3日开始,2017年4月1日完全融化,最大雪厚12.4 cm,模拟雪厚与实测值比较吻合,土壤冻结融化期间分别对试验点土壤温度、含水率进行了模拟对比。

图7 试验点冻融期间温度及积雪厚度Fig.7Temperature and snow thickness during freezing and thawing at the experimental site

3.2.1 温度模拟

图8是试验点2016—2017年冻结融化期垂向不同深度的温度模拟值和实测值的比较结果,总的来说,

图8 不同深度土壤和砂砾石层温度模拟与实测对比Fig.8Simulated and observed soil and loose rock strata temperatures at different depths

模拟温度与实测温度较为接近,但表层(10 cm、20 cm深)温度波动较大,模拟效果稍差,R2分别为0.83、0.89,ERMS较大分别为4.07 ℃、3.39 ℃；表层以下(40 cm深及以下)温度较为稳定,模拟效果较好,R2在0.96以上,ERMS均小于10 cm和20 cm层,平均ERMS为2.86 ℃。表层温度模拟效果较差,可能与气温波动较大有关,由于该流域气象站点较少,温度空间插值很难与当地实际气温相符。另外,80 cm深度以下的模拟温度略小于实际温度,可能与温度下边界条件的设置有关,需要在今后的研究中进一步改进。

3.2.2 含水率模拟

图9为试验点2016—2017年冻结融化期液态含水率模拟和实测值对比。由图中可以看出,在冻结阶段(11月—翌年1月),受温度降低的影响,首先是上层土壤-砂砾石层的液态含水量变小,随后下层砂砾石层的液态含水率变小,到1月下旬达到稳定。土壤层(40 cm以上)由于其持水能力大于砂砾石层,初始的含水率更大,且受空气温度影响也大,在冻结过程中含水率减小明显。2月下旬随着气温回暖,先是上层土壤-砂砾石层开始融化,液态含水量增加,随后土壤-砂砾石层底部也开始融化。在融化期,受积雪融化下渗影响,土壤-砂砾石层上部土壤层的含水量比下部砂砾石大,比冻结前也大。在快速冻结期(2016年12月18日前)R2均值为0.38,此后稳定冻结期与融化期R2均值为0.58。总的来说,模拟的分层土壤-砂砾石层的液态

图9 不同深度土壤液态含水率模拟与实测对比Fig.9Simulated and observed soil moisture contents at different depths

含水率与实测值基本接近,可以反映冻结融化期各阶段的变化趋势。由于土壤和砂砾石层中组成成分的不确定性较大,模型进行概化时,认为砂砾石层和土壤是均质结构,无法准确反映土壤-砂砾石层持水能力的不稳定,这也导致了模拟结果与实测值具有一定的误差。

3.2.3 冻结深度模拟

试验点2016—2017年冻结融化期冻结深度模拟结果见图10。试验点从11月下旬左右进入冻结期,3月上旬左右融通,冻结融化期总长在4个月左右,最大冻结深度约为67.2 cm。R2为0.76,ERMS为16.36 cm。根据实测冻结深度与模拟值的对比结果来看,试验点处模拟的冻结规律与实测结果基本一致。模拟的冻结融化期起始时间与实测结果相比略微偏早,表层土壤开始融化时间比实测略晚。

图10 土壤和砂砾石层冻结融化期冻结深度的模拟与实测对比Fig.10Simulated and observed soil and loose rock freezing depths

4 结 论

本文在WEP-COR模型冻土水热耦合计算原理基础上,结合青藏高原土壤层较薄、下伏砂砾石层较厚以及积雪期较长的特点加以改进,构建了“积雪-土壤-砂砾石层”连续体水热耦合模型。

(1) 模型主要具有以下特点:将均质的土壤结构细化分解成土壤和砂砾石层2类介质,水热耦合计算考虑了土壤和砂砾石层的冻融过程；在原模型土壤层的上方加入了一层积雪层,考虑了积雪对大气与土壤间热量传递的影响；将积雪、土壤和砂砾石层概化为一个紧密的一维连续结构,提出了“积雪-土壤-砂砾石层”连续体分层水热耦合模拟方法。

(2) 基于尼洋河流域下游色季拉山试验点的监测结果对本文构建的模型进行验证,各层温度模拟R2均值为0.91,冻结融化期内含水率模拟R2均值为0.52,土壤冻结深度模拟R2值为0.76,土壤和砂砾石层不同深度的温度、液态含水率和冻结深度的模拟结果与实测值接近。模型能较好地反映冻结融化期各阶段土壤和砂砾石层的分层水热变化趋势。

(3) 由于高寒山区试验条件恶劣,目前只设置了一个试验站点,且本次取得的监测资料时段较短。另外,因研究流域气象站点较少,资料较为缺乏,加上模型的下边界条件、“积雪-土壤-砂砾石层”连续体各项水热参数均存在一定程度的概化,导致模拟结果与实测值之间仍存在一定偏差,这些都需要在未来的研究中进一步完善。

(4) 青藏高原“积雪-土壤-砂砾石层”连续体的周期性冻融过程影响着流域产流的各个环节,对水文循环过程影响很大。为研究水热耦合运移对流域水循环的影响还需将其与水文模型结合进一步分析。