-
目前, 基于寒区水热耦合的特点, 分布式水文模型已取得一定的进展[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]的土壤水热耦合模块为基础, 构建基于“积雪-土壤-砂砾石层”连续体的水热耦合模型, 并采用试验观测结果对模型进行验证。
-
选取雅鲁藏布江的一级支流尼洋河流域作为本文研究的典型区域。尼洋河位于西藏东南林芝地区雅鲁藏布江中下游左岸, 北纬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)和部分工程设计的地质勘测资料, 该流域具有土壤层厚度较薄、土壤下方的砂砾石层厚度较大的特点。尼洋河流域的砂砾石层主要包括砾石层和卵石层, 砾石层主要为圆状砾石, 含卵石颗粒, 砾石质量百分比约50%~65%, 黏粒质量百分比5%~10%, 孔隙间填充中细砂；卵石层孔隙间填充圆砾和中细砂, 局部含大漂石、孤石。
图 1 土壤-砂砾石层结构
Figure 1. Soil-sand gravel structure
针对尼洋河流域这一地质特点, 为研究流域内土壤层从山脚到山顶的厚度变化规律, 从尼洋河流域的源头到河口选取了海拔不同的24个取样点。其中点1-点16沿河道, 点17-点24从山脚到山顶测量的土壤厚度(图 2)。对24个取样点的土壤厚度及组分进行了测量和分析, 尼洋河流域土壤主要为砂质壤土和壤质砂土, 平均砂粒质量百分比为55.89%, 粉砂粒质量百分比为31.2%, 黏粒质量百分比为12.91%, 厚度随着海拔升高而减小, 平均厚度为53.7 cm。
图 2 尼洋河流域土壤取样点、试验点及气象站位置
Figure 2. Soil sampling points and meteorological stations
-
本次研究在尼洋河流域下游色季拉山的山腰上选择一处开展季节冻结融化期水热耦合过程监测试验。试验点(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周根据降雪情况通过米尺实地测量下垫面覆盖的雪层厚度。试验前, 采用核磁共振方法对高原季节性冻融条件下水、热过程监测的仪器设备进行率定。仪器完成安装后进行原状土回填, 保证回填土壤与实际土壤的参数一致。
-
气象条件是水热耦合模拟的边界条件。由于流域内气象站点稀少, 且高程差别较大, 试验点的逐日气温采用考虑高程修正的RDS(Reversed Distance Squared)方法由气象站数据插值得到, 降水数据采用平面插值、结合西藏自治区水资源公报[25]中降水等值线修正得到, 逐日相对湿度、日照时速和风速采用RDS方法插值得到。本文用到的基本气象数据主要为尼洋河流域境内的林芝站和周边嘉黎站(图 2)的降水、气温、相对湿度、日照时数和风速。气象数据来源于中国气象数据网(http://data.cma.cn), 使用数据时间序列为2016-2017年, 时间尺度为日。
-
本文在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模型分层水热通量计算结构示意
Figure 3. WEP-COR model water and heat flux simulation structure
-
模型假设土壤冻融时只有液态水发生运移, 土壤水分的运移主要受重力势、基质势和温度势的影响, 因此添加的土壤一维垂直水分流方程可以写成[26]:
$$ \frac{{\partial {\theta _1}}}{{\partial t}} = \frac{\partial }{{\partial z}}\left[ {D\left( {{\theta _1}} \right)\frac{{\partial {\theta _1}}}{{\partial z}} - K\left( {{\theta _1}} \right)} \right] - \frac{{{\rho _{\rm{i}}}}}{{{\rho _1}}}\frac{{\partial {\theta _{\rm{i}}}}}{{\partial t}} $$ (1) 根据能量平衡原理, 冻融系统中每一层土壤的能量变化都用于系统内的土壤温变和水分相变, 温度势是水分相变的驱动力, 而大气与表层土壤的温差则是热传导的源动力。假设土壤各向均质同性, 并忽略土壤中的水汽迁移, 添加的一维垂直流动的热流运动基本方程可变形为[27]
$$ {C_{\rm{v}}}\frac{{\partial T}}{{\partial t}} = \frac{\partial }{{\partial z}}\left[ {\lambda {\kern 1pt} {\kern 1pt} \frac{{\partial T}}{{\partial z}}} \right] + {L_{\rm{f}}}{\rho _{\rm{i}}}\frac{{\partial {\theta _{\rm{i}}}}}{{\partial t}} $$ (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×106 J/kg。
-
土壤冻融过程中, 冻土水热运动间的联系主要表现在未冻水的含水率与土壤负温的动态平衡中[19-20]:
$$ {\theta _1} = {\theta _{\rm{m}}}(T) $$ (3) 式中:θm(T)为土壤负温对应的最大未冻水含水率。
-
积雪作为一种良好的隔热体, 具有热容量大、导热系数较小等特点, 对大气和土壤间的热传递具有明显的削弱作用。试验点所处位置积雪持续时间较长、厚度稳定, 对土壤水热的影响不容忽视。根据尼洋河流域实地考察(图 1)和基于26个点的测量结果建立的土壤厚度-高程关系图(图 4), 发现尼洋河流域土壤厚度整体较薄, 在土壤层下方存在厚度较大、混合砂砾石和少量土壤的砂砾石层。而且土壤层厚度在河谷区基本稳定, 在山坡上随着高程增加呈现变薄的趋势, 在山顶附近厚度又基本保持不变。在试验点处土壤层为砂质壤土, 砂粒质量百分比为61.92%, 粉砂粒质量百分比为26.2%, 黏粒质量百分比为11.88%, 土壤层厚度为40 cm；土壤下砂砾石层砾石质量百分比为53.7%, 孔隙间填充中细砂。土壤与砂砾石层水分特征参数见表 1。
图 4 尼洋河流域土壤厚度-海拔关系
Figure 4. Relationship between soil thickness and elevation in Niyang River basin
表 1 尼洋河流域土壤水动力参数
Table 1. Soil moisture characteristics
类型 θs/(cm3·cm-3) θf/(cm3·cm-3) θr/(cm3·cm-3) Sf/cm n 土壤层 0.43 0.25 0.06 6.1 3.37 砂砾石层 0.28 0.13 0.05 4.8 3.12 注:θs、θf和θr分别为土壤或砂砾石层的饱和含水率、田间持水率、残留含水率；Sf为湿润峰土壤吸力；n为Mualem常数。 根据尼洋河流域气候和地质特征, 本文将水热耦合模拟的对象定义为“积雪-土壤-砂砾石层”连续体(图 5)。模型计算结构上, 在WEP-COR模型11层土壤层分层的基础上, 于土壤层上方加入一层积雪层, 同时将原有的11层土壤层细分为土壤和砂砾石层2类介质, 计算结构如图 6所示。第0层为积雪, 根据等高带在计算单元中所处位置设置上部i层为土壤层, 下部12-i层为砂砾石层。由于距离地表越近, 水、热变化越快, 故第1、2层厚度设为10 cm, 第3~11层厚度设为20 cm。值得注意的是, 土壤-砂砾石层各层的厚度需根据所在位置含水层的实际厚度进行修正, 也就是说, 当含水层总厚度小于11层土壤-砂砾石层厚度时, 须根据含水层实际厚度确定土壤-砂砾石层的计算层数和最后一层的厚度。
图 5 “积雪-土壤-砂砾石层”垂向分层概化
Figure 5. Vertical layering of "snow-soil-sand gravel layer"
图 6 分层水热通量计算结构示意
Figure 6. Layered calculation structure of water and heat flux
-
采用度日因子法[28]计算积雪融雪过程, 具体如下:
$$ M = k\left( {{T_{\rm{a}}} - {T_0}} \right) $$ (4) 式中:M为当日融雪当量, mm；k为积雪融化的度日因子, mm/(℃·d)；T0为积雪融化临界温度, ℃, 本文取-1 ℃。
对于积雪的度日因子, 一般情况下在1~7 mm/(℃·d), 本次研究根据下垫面情况取为4 mm/(℃·d)。由于青藏高原地区昼夜温差较其他地区偏大, 本次研究的雨雪临界气温设为2 ℃。
积雪热量平衡公式[29]:
$$ G = {L_{\rm{e}}} + {R_{\rm{n}}} + {Q_1} + {Q_2} $$ (5) 式中:G为积雪层的能量变化, J/(m2·d)；Le为积雪的蒸发、融化和升华潜热通量, J/(m2·d)；Rn为积雪层表面的净辐射量, J/(m2·d), 由到达积雪表面的净长波辐射与净短波辐射计算而得；Q1为积雪与大气间的热量交换量, J/(m2·d), 采用强迫-恢复法计算；Q2为积雪与土壤间的热量交换量, J/(m2·d)。
$$ {Q_2} = \frac{1}{{\frac{{{Z_{\rm{S}}}}}{{2{\lambda _{\rm{S}}}}} + \frac{{{Z_{\rm{C}}}}}{{2{\lambda _{\rm{C}}}}} + {R_{\rm{C}}}}}\left( {{T_{\rm{C}}} - {T_{\rm{S}}}} \right) $$ (6) 式中:ZS为积雪厚度, m；ZC为第一层土壤厚度, m；λS为冰雪层导热系数, W/(m·℃)；λC为土壤层导热系数, W/(m·℃)；RC为土壤层与积雪层的接触热阻, (m2·℃)/W；TC为第一层土壤层温度, ℃；TS为积雪层温度, ℃。
土壤层之间的水热运移方程和联系方程采用式(1)和式(2), 砂砾石层之间以及砂砾石层和土壤层之间的水热通量计算采用与土壤层相同的计算原理, 但砂砾石层的各项水热参数与土壤不同。
模型的数值模拟仍采用显式差分进行数值迭代计算, 模拟步长为1 d。每个时间步长内模型先根据初始条件从积雪层到底层砂砾石层逐层计算热传导、温度和水分相变, 然后以该层温度为判定条件进行迭代计算直至收敛(步骤1)；热量计算闭合收敛后, 再根据水量平衡计算各层水分运移量, 并修正各土壤和砂砾石层(积雪层不考虑)的含水率和含冰率(步骤2)；用液态含水率判定是否收敛, 若不收敛则返回步骤1进行热量计算, 直至两项迭代计算均闭合收敛后, “积雪-土壤-砂砾石层”连续体水热耦合计算完成。
-
模型上、下边界的确定方式与WEP-COR模型类似。“积雪-土壤-砂砾石层”连续体的上边界为大气, 但与上边界直接接触的不一定只是土壤。传入连续体的热量, 在非积雪区由地表附近的大气和表层土壤的温度差以及表层土壤的水热参数决定, 在积雪区由地表附近的大气和积雪的温度差以及积雪的水热参数决定, 均采用强迫-恢复法计算。连续体的下边界为过渡层或者地下水层, 将下边界温度假设为一恒定的温度作为底层边界温度。
-
积雪主要的水热参数包括导热系数、体积热容量和降雪密度, 各参数计算公式如下[32-34]:
$$ {\rho _{\rm{s}}} = \left\{ {\begin{array}{*{20}{l}} {67.9 + 51.3\exp \left( {{T_{\rm{a}}}/2.6} \right)}&{{T_{\rm{a}}} \le 0}\\ {119.2 + 20{T_{\rm{a}}}}&{{T_{\rm{a}}} > 0} \end{array}} \right. $$ (7) $$ {\lambda _{\rm{s}}} = \left\{ {\begin{array}{*{20}{l}} {0.138 - 1.01 \times \frac{{{\rho _{\rm{s}}}}}{{1{\kern 1pt} {\kern 1pt} {\kern 1pt} 000}} + 3.233 \times {{\left( {\frac{{{\rho _{\rm{s}}}}}{{1{\kern 1pt} {\kern 1pt} {\kern 1pt} 000}}} \right)}^2}}&{156 < {\rho _{\rm{s}}} \le 600}\\ {0.023 + \frac{{0.234{\rho _{\rm{s}}}}}{{1{\kern 1pt} {\kern 1pt} {\kern 1pt} 000}}}&{{\rho _{\rm{s}}} \le 156} \end{array}} \right. $$ (8) $$ {C_{\rm{s}}} = 2.09 \times {10^3}{\rho _{\rm{s}}} $$ (9) 式中:ρs为积雪密度, kg/m3；λs为积雪的导热系数, W/(m·℃)；Cs为积雪的体积热容量, J/(m3·℃)。土壤、砂砾石层主要的水热参数也包括体积热容量、导热系数和土壤导水率等, 各参数计算公式如下。体积热容[35]:
$$ {C_{\rm{v}}} = \left( {1 - {\theta _{\rm{s}}}} \right){C_{{\rm{sg}}}} + {\theta _{\rm{l}}}{C_{\rm{l}}} + {\theta _{\rm{i}}}{C_{\rm{i}}} $$ (10) 导热系数计算参考了IBIS模型, 具体如下[36]:
$$ \begin{array}{*{20}{c}} {\lambda = {\lambda _{{\rm{st}}}}\left( {{{56}^{{\theta _l}}} + {{224}^{{\theta _{\rm{i}}}}}} \right)}\\ {{\lambda _{{\rm{st}}}} = 1.500{\omega _{{\rm{rock}}}} + 0.300{\omega _{{\rm{sand}}}} + 0.265{\omega _{{\rm{silt}}}} + 0.250{\omega _{{\rm{clat}}}}} \end{array} $$ (11) 导水率[29]:
$$ K\left( {{\theta _{\rm{l}}}} \right) = \left\{ {\begin{array}{*{20}{l}} {{K_{\rm{S}}}}&{{\theta _{\rm{l}}} = {\theta _{\rm{s}}}}\\ {{K_{\rm{S}}}{{\left( {\frac{{{\theta _{\rm{l}}} - {\theta _{\rm{r}}}}}{{{\theta _{\rm{s}}} - {\theta _{\rm{r}}}}}} \right)}^n}}&{{\theta _{\rm{l}}} \ne {\theta _{\rm{s}}}} \end{array}} \right. $$ (12) 式中:Csg、Cl和Ci分别为土壤(砂砾石层)、水、冰的体积热容量, 在0℃时土壤和砂砾石层分别为1.93×103 J/(m3·℃)、3.1×103 J/(m3·℃), 水和冰分别为4.213×103 J/(m3·℃)、1.94×103 J/(m3·℃)；λ、λst为土壤或砂砾石层实际的导热系数、干燥状态下的导热系数, W/(m·℃)；ωrock、ωsand、ωsilt和ωclat分别为岩石(砾石和卵石)、砂粒、粉粒和黏粒的体积比；KS为土壤或砂砾石层经温度修正后的饱和含水率, cm/s。
参考陈仁升等的DWHC模型[37], 不同温度条件下土壤或砂砾石层KS计算方法如下:
$$ {K_{\rm{S}}} = \left\{ {\begin{array}{*{20}{l}} {{K_0}}&{T > 0}\\ {{K_0}(0.54 + 0.023T)}&{{T_{\rm{f}}} \le T \le 0}\\ {{k_0}}&{T < {T_{\rm{f}}}} \end{array}} \right. $$ (13) 式中:K0为常温下的饱和导水率, cm/s；k0为冻结条件下最小的导水率, cm/s；Tf为土壤或砂砾石层最小导水率对应的临界温度, ℃。本文考虑到土壤和砂砾石层水动力学特性的区别, 对于土壤, k0按照原公式取值为0 cm/s；对于砂砾石层, 由于孔隙较大, k0取值为大于0的值。
-
本文采用决定系数(R2)及平均根方差(ERMS)评价模型模拟效果。R2和ERMS的计算公式如下:
$$ {R^2} = 1 - \frac{{\sum\limits_{i = 1}^N {{{\left( {{O_i} - {S_i}} \right)}^2}} }}{{\sum\limits_{i = 1}^N {{{\left( {{O_i} - \overline {{O_i}} } \right)}^2}} }} $$ (14) $$ {E_{{\rm{RMS}}}} = \sqrt {\frac{1}{n}\sum\limits_{i = 1}^N {{{\left( {{O_i} - {S_i}} \right)}^2}} } $$ (15) 式中:N代表样本数；Oi为观测值；Oi为观测值的平均数；Si为模拟值。对于ERMS, 越小说明模拟与实测的差值越小, 模拟效果越好。对于R2, 数值越接近1模拟效果越好。
-
试验点冻结融化期间空气温度、模拟雪厚、实测雪厚如图 7所示, 试验点积雪从2016年12月3日开始, 2017年4月1日完全融化, 最大雪厚12.4 cm, 模拟雪厚与实测值比较吻合, 土壤冻结融化期间分别对试验点土壤温度、含水率进行了模拟对比。
图 7 试验点冻融期间温度及积雪厚度
Figure 7. Temperature and snow thickness during freezing and thawing at the experimental site
-
图 8是试验点2016-2017年冻结融化期垂向不同深度的温度模拟值和实测值的比较结果, 总的来说, 模拟温度与实测温度较为接近, 但表层(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深度以下的模拟温度略小于实际温度, 可能与温度下边界条件的设置有关, 需要在今后的研究中进一步改进。
图 8 不同深度土壤和砂砾石层温度模拟与实测对比
Figure 8. Simulated and observed soil and loose rock strata temperatures at different depths
-
图 9为试验点2016-2017年冻结融化期液态含水率模拟和实测值对比。由图中可以看出, 在冻结阶段(11月-翌年1月), 受温度降低的影响, 首先是上层土壤-砂砾石层的液态含水量变小, 随后下层砂砾石层的液态含水率变小, 到1月下旬达到稳定。土壤层(40 cm以上)由于其持水能力大于砂砾石层, 初始的含水率更大, 且受空气温度影响也大, 在冻结过程中含水率减小明显。2月下旬随着气温回暖, 先是上层土壤-砂砾石层开始融化, 液态含水量增加, 随后土壤-砂砾石层底部也开始融化。在融化期, 受积雪融化下渗影响, 土壤-砂砾石层上部土壤层的含水量比下部砂砾石大, 比冻结前也大。在快速冻结期(2016年12月18日前)R2均值为0.38, 此后稳定冻结期与融化期R2均值为0.58。总的来说, 模拟的分层土壤-砂砾石层的液态含水率与实测值基本接近, 可以反映冻结融化期各阶段的变化趋势。由于土壤和砂砾石层中组成成分的不确定性较大, 模型进行概化时, 认为砂砾石层和土壤是均质结构, 无法准确反映土壤-砂砾石层持水能力的不稳定, 这也导致了模拟结果与实测值具有一定的误差。
图 9 不同深度土壤液态含水率模拟与实测对比
Figure 9. Simulated and observed soil moisture contents at different depths
-
试验点2016-2017年冻结融化期冻结深度模拟结果见图 10。试验点从11月下旬左右进入冻结期, 3月上旬左右融通, 冻结融化期总长在4个月左右, 最大冻结深度约为67.2 cm。R2为0.76, ERMS为16.36 cm。根据实测冻结深度与模拟值的对比结果来看, 试验点处模拟的冻结规律与实测结果基本一致。模拟的冻结融化期起始时间与实测结果相比略微偏早, 表层土壤开始融化时间比实测略晚。
图 10 土壤和砂砾石层冻结融化期冻结深度的模拟与实测对比
Figure 10. Simulated and observed soil and loose rock freezing depths
-
本文在WEP-COR模型冻土水热耦合计算原理基础上, 结合青藏高原土壤层较薄、下伏砂砾石层较厚以及积雪期较长的特点加以改进, 构建了“积雪-土壤-砂砾石层”连续体水热耦合模型。
(1) 模型主要具有以下特点:将均质的土壤结构细化分解成土壤和砂砾石层2类介质, 水热耦合计算考虑了土壤和砂砾石层的冻融过程；在原模型土壤层的上方加入了一层积雪层, 考虑了积雪对大气与土壤间热量传递的影响；将积雪、土壤和砂砾石层概化为一个紧密的一维连续结构, 提出了“积雪-土壤-砂砾石层”连续体分层水热耦合模拟方法。
(2) 基于尼洋河流域下游色季拉山试验点的监测结果对本文构建的模型进行验证, 各层温度模拟R2均值为0.91, 冻结融化期内含水率模拟R2均值为0.52, 土壤冻结深度模拟R2值为0.76, 土壤和砂砾石层不同深度的温度、液态含水率和冻结深度的模拟结果与实测值接近。模型能较好地反映冻结融化期各阶段土壤和砂砾石层的分层水热变化趋势。
(3) 由于高寒山区试验条件恶劣, 目前只设置了一个试验站点, 且本次取得的监测资料时段较短。另外, 因研究流域气象站点较少, 资料较为缺乏, 加上模型的下边界条件、“积雪-土壤-砂砾石层”连续体各项水热参数均存在一定程度的概化, 导致模拟结果与实测值之间仍存在一定偏差, 这些都需要在未来的研究中进一步完善。
(4) 青藏高原“积雪-土壤-砂砾石层”连续体的周期性冻融过程影响着流域产流的各个环节, 对水文循环过程影响很大。为研究水热耦合运移对流域水循环的影响还需将其与水文模型结合进一步分析。
Distributed hydrological model of the Qinghai Tibet Plateau based on the hydrothermal coupling: Ⅰ: hydrothermal coupling simulation of "snow-soil-sand gravel layer" continuum
-
摘要: 青藏高原号称"亚洲水塔"，是典型的高寒山区，其广泛存在的积雪、多年和季节冻土，影响了整个区域的水循环过程；青藏高原具有土层较薄、下伏砂砾石层较厚的特点，形成了特殊的"积雪-土壤-砂砾石层"水热介质结构。为深入研究青藏高原的水循环机理，本文选取尼洋河流域作为典型区，基于野外冻土水热耦合试验，结合青藏高原地质及气候特点，构建了包含12层"积雪-土壤-砂砾石层"连续体的青藏高原水热耦合模型，描述了完整的水热耦合模拟方程和参数计算方法。采用2016-2017年冻结融化期0~160 cm深度内的土壤和砂砾石层的温度、液态含水率和冻结深度的实测结果对模型进行了验证，各层温度模拟R2均值为0.91，冻结融化期内液态含水率模拟R2均值为0.52，土壤冻结深度模拟R2值为0.76。结果表明该模型在青藏高原地区有较好的适用性，可反映该地区冻结融化过程中土壤和砂砾石层水分与温度的特殊变化规律。
-
关键词:
- "积雪-土壤-砂砾石层"连续体 /
- 水热耦合模拟 /
- 青藏高原 /
- 高寒山区 /
- 积雪
Abstract: The Qinghai-Tibet Plateau, known as the "Asian Water Tower", is a typical alpine mountain area. Snow-covered permafrost and seasonally frozen soil widely exist on this plateau, profoundly affecting the water cycle process of the entire region. In addition, the Qinghai-Tibet Plateau has the characteristics of a thinner soil layer and a thicker underlying sand and gravel layer, forming a special "snow-soil-sand gravel layer" hydrothermal medium structure. To further study the water cycle mechanism of the Qinghai-Tibet Plateau, this study selected the Niyang River basin as a typical area to perform field-heat-coupling experiments in the field. Combined with the geological and climatic characteristics of the Qinghai-Tibet Plateau, a 12-layer "snow-soil-sand gravel layer" continuum was constructed to describe the hydro-thermal coupling model of the Qinghai-Tibet Plateau. A complete hydro-thermal coupling simulation equation and parameter calculation method are described. The model was validated using the measured results of the temperature, liquid water content, and freezing depth of the soil and gravel layer in the freezing and thawing period of 0-160 cm from 2016 to 2017. The mean values of the simulated R2 of the temperature of each layer and the moisture content during freezing and thawing were 0.91 and 0.52, respectively. The R2 of the soil freezing depth is 0.76. The results show that the model can better reflect the hydrothermal variation of stratified soil and sand gravel during freezing and thawing processes. The results show that the model has adequate applicability in the Qinghai-Tibet Plateau and can reflect the special change process of moisture and temperature of the soil and gravel layer during the freezing and melting processes in this area. -
图 3 WEP-COR模型分层水热通量计算结构示意
Figure 3. WEP-COR model water and heat flux simulation structure
图 4 尼洋河流域土壤厚度-海拔关系
Figure 4. Relationship between soil thickness and elevation in Niyang River basin
图 7 试验点冻融期间温度及积雪厚度
Figure 7. Temperature and snow thickness during freezing and thawing at the experimental site
图 8 不同深度土壤和砂砾石层温度模拟与实测对比
Figure 8. Simulated and observed soil and loose rock strata temperatures at different depths
图 9 不同深度土壤液态含水率模拟与实测对比
Figure 9. Simulated and observed soil moisture contents at different depths
图 10 土壤和砂砾石层冻结融化期冻结深度的模拟与实测对比
Figure 10. Simulated and observed soil and loose rock freezing depths
表 1 尼洋河流域土壤水动力参数
Table 1. Soil moisture characteristics
类型 θs/(cm3·cm-3) θf/(cm3·cm-3) θr/(cm3·cm-3) Sf/cm n 土壤层 0.43 0.25 0.06 6.1 3.37 砂砾石层 0.28 0.13 0.05 4.8 3.12 注:θs、θf和θr分别为土壤或砂砾石层的饱和含水率、田间持水率、残留含水率；Sf为湿润峰土壤吸力；n为Mualem常数。 
下载: 导出CSV
-
[1] LI J, ZHOU Z H, WANG H, et al. Development of WEP-COR model to simulate land surface water and energy budgets in a cold region[J]. Hydrology Research, 2019, 50(1):99-116. doi: 10.2166/nh.2017.032 [2] ZHANG Y L, CHENG G D, LI X, et al. Coupling of a simultaneous heat and water model with a distributed hydrological model and evaluation of the combined model in a cold region watershed[J]. Hydrological Processes, 2013, 27(25):3762-3776. doi: 10.1002/hyp.9514 [3] POMEROY J W, GRAY D M, BROWN T, et al. The cold regions hydrological model:a platform for basing process representation and model structure on physical evidence[J]. Hydrological Processes, 2007, 21(19):2650-2667. doi: 10.1002/hyp.6787 [4] RIGON R, BERTOLDI G, OVER T M. GEOtop:a distributed hydrological model with coupled water and energy budgets[J]. Journal of Hydrometeorology, 2006, 7(3):371-388. doi: 10.1175/JHM497.1 [5] HARLAN R L. Analysis of coupled heat-fluid transport in partially frozen soil[J]. Water Resources Research, 1973, 9(5):1314-1323. doi: 10.1029/WR009i005p01314 [6] TAYLOR G S, LUTHIN J N. A model for coupled heat and moisture transfer during soil freezing[J]. Canadian Geotechnical Journal, 1978, 15(4):548-555. doi: 10.1139/t78-058 [7] HORIGUCHI K. An osmotic model for soil freezing[J]. Cold Regions Science and Technology, 1987, 14(1):13-22. doi: 10.1016/0165-232X(87)90040-1 [8] FLERCHINGER G N, SAXTON K E. Simultaneous heat and water model of a freezing snow-residue-soil system:Ⅱ:Field verification[J]. Transactions of the ASAE, 1989, 32(2):573-576. doi: 10.13031/2013.31041 [9] 雷志栋, 尚松浩, 杨诗秀, 等.地下水浅埋条件下越冬期土壤水热迁移的数值模拟[J].冰川冻土, 1998, 20(1):51-54. https://www.cnki.com.cn/Article/CJFDTOTAL-BCDT801.010.htm LEI Z D, SHANG S H, YANG S X, et al. Numerical simulation on simultaneous soil moisture and heat transfer under shallow ground water table in winter[J]. Journal of Glaciology and Geocryology, 1998, 20(1):51-54.(in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-BCDT801.010.htm [10] 李炳元, 潘保田, 程维明, 等.中国地貌区划新论[J].地理学报, 2013, 68(3):291-306. https://www.cnki.com.cn/Article/CJFDTOTAL-DLXB201303003.htm LI B Y, PAN B T, CHENG W M, et al. Research on geomorphological regionalization of China[J]. Acta Geographica Sinica, 2013, 68(3):291-306.(in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-DLXB201303003.htm [11] 程根伟, 范继辉, 彭立.高原山地土壤冻融对径流形成的影响研究进展[J].地球科学进展, 2017, 32(10):1020-1029. doi: 10.11867/j.issn.1001-8166.2017.10.1020 CHENG G W, FAN J H, PENG L. Progresses in soil freezing-thawing effects on the runoff generation in plateau-mountain regions[J]. Advances in Earth Science, 2017, 32(10):1020-1029.(in Chinese) doi: 10.11867/j.issn.1001-8166.2017.10.1020 [12] ZHOU B B, SHAO M A, SHAO H B. Effects of rock fragments on water movement and solute transport in a Loess Plateau soil[J]. Comptes Rendus Geoscience, 2009, 341(6):462-472. doi: 10.1016/j.crte.2009.03.009 [13] 何玉洁, 宜树华, 郭新磊.青藏高原含砂砾石土壤导热率实验研究[J].冰川冻土, 2017, 39(2):343-350. https://www.cnki.com.cn/Article/CJFDTOTAL-BCDT201702014.htm HE Y J, YI S H, GUO X L. Experimental study on thermal conductivity of soil with gravel on the Qinghai-Tibet Plateau[J]. Journal of Glaciology and Geocryology, 2017, 39(2):343-350.(in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-BCDT201702014.htm [14] 曹伟, 盛煜, 吴吉春, 等.青藏高原坡面冻土土壤水分空间变异特性[J].水科学进展, 2017, 28(1):32-40. doi: 10.14042/j.cnki.32.1309.2017.01.004 CAO W, SHENG Y, WU J C, et al. Spatial variability of permafrost soil-moisture on the slope of the Qinghai-Tibet Plateau[J]. Advances in Water Science, 2017, 28(1):32-40.(in Chinese) doi: 10.14042/j.cnki.32.1309.2017.01.004 [15] 朱元骏, 邵明安.含砾石土壤降雨入渗过程模拟[J].水科学进展, 2010, 21(6):779-787. http://skxjz.nhri.cn/article/id/200 ZHU Y J, SHAO M A. Simulation of rainfall infiltration in stony soil[J]. Advances in Water Science, 2010, 21(6):779-787.(in Chinese) http://skxjz.nhri.cn/article/id/200 [16] PAN Y J, LYU S H, LI S S, et al. Simulating the role of gravel in freeze-thaw process on the Qinghai-Tibet Plateau[J]. Theoretical and Applied Climatology, 2017, 127(3/4):1011-1022. doi: 10.1007/s00704-015-1684-7?Code=363c8d46-3ea1-473e-8b37-05b80940a84e [17] YANG K, CHEN Y Y, QIN J. Some practical notes on the land surface modeling in the Tibetan Plateau[J]. Hydrology and Earth System Sciences, 2009, 13(5):687-701. doi: 10.5194/hess-13-687-2009 [18] 王子龙, 何馨, 姜秋香, 等.气候变化下东北中等流域冬季径流模拟和预测[J].水科学进展, 2020, 31(4):575-582. doi: 10.14042/j.cnki.32.1309.2020.04.011 WANG Z L, HE X, JIANG Q X, et al. Simulation and forecast of winter runoff in medium basin of Northeast China under climate change[J]. Advances in Water Science, 2020, 31(4):575-582. (in Chinese) doi: 10.14042/j.cnki.32.1309.2020.04.011 [19] PITMAN A J, YANG Z L, HENDERSON-SELLERS A. Description of bare essentials of surface transfer for the Bureau of Meteorology Research Centre AGCM[R]. Melbourne:Bureau of Meteorology Research Centre, 1991:56-57. [20] DOUVILLE H, ROYER J F, MAHFOUF J F. A new snow parameterization for the Météo-France climate model[J]. Climate Dynamics, 1995, 12(1):21-35. [21] FLERCHINGER G N, COOLEY K R, DENG Y. Impacts of spatially and temporally varying snowmelt on subsurface flow in a mountainous watershed:1:snowmelt simulation[J]. Hydrological Sciences Journal, 1994, 39(5):507-520. doi: 10.1080/02626669409492771 [22] LOTH B, GRAF H F, OBERHUBER J M. Snow cover model for global climate simulations[J]. Journal of Geophysical Research:Atmospheres, 1993, 98(D6):10451-10464. doi: 10.1029/93JD00324 [23] LYNCH-STIEGLITZ M. The development and validation of a simple snow model for the GISS GCM[J]. Journal of Climate, 1994, 7(12):1842-1855. doi: 10.1175/1520-0442(1994)007<1842:TDAVOA>2.0.CO;2 [24] 达瓦次仁, 巴桑赤烈, 白玛, 等.尼洋河流域水文特性分析[J].水文, 2008, 28(4):92-94. doi: 10.3969/j.issn.1000-0852.2008.04.026 DAWA C R, BASANG C L, BAI M, et al. Analysis of hydrological characteristics of Niyang River basin[J]. Journal of China Hydrology, 2008(4):92-94.(in Chinese) doi: 10.3969/j.issn.1000-0852.2008.04.026 [25] 西藏自治区水利厅.西藏自治区水资源公报(2017)[R].拉萨: 西藏自治区水利厅, 2017. Tibet Autonomous Region Water Conservancy Department. Tibet water resources bulletin(2017)[R]. Lhasa: Tibet Autonomous Region Water Conservancy Department, 2017.(in Chinese)) [26] 尚松浩, 雷志栋.冻结条件下土壤水热耦合迁移数值模拟的改进[J].清华大学学报(自然科学版), 1997, 37(8):62-64. https://www.cnki.com.cn/Article/CJFDTOTAL-QHXB708.017.htm SHANG S H, LEI Z D. Numerical simulation improvement of coupled moisture and heat transfer during soil freezing[J]. Journal of Tsinghua University(Science and Technology), 1997, 37(8):62-64. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-QHXB708.017.htm [27] 雷志栋, 杨诗秀, 谢森传.土壤水动力学[M].北京:清华大学出版社, 1988:74-75. LEI Z D, YANG S X, XIE S C. Soil hydrodynamics[M]. Beijing:Tsinghua University Press, 1988:74-75.(in Chinese) [28] 张勇, 刘时银.度日模型在冰川与积雪研究中的应用进展[J].冰川冻土, 2006, 28(1):101-107. doi: 10.3969/j.issn.1000-0240.2006.01.014 ZHANG Y, LIU S Y. Progress of the application of degree-day model to study glaciers and snow cover[J]. Journal of Glaciology and Geocryology, 2006, 28(1):101-107. (in Chinese) doi: 10.3969/j.issn.1000-0240.2006.01.014 [29] JIA Y W, NI G H, KAWAHARA Y, et al. Development of WEP model and its application to an urban watershed[J]. Hydrological Processes, 2001, 15(11):2175-2194. doi: 10.1002/hyp.275 [30] 藤本明宏, 渡邊洋, 福原輝幸.多層路面雪氷状態モデル[J].日本雪工学会誌, 2008, 24(1):3-11. FUJIMOTO A, WATANABE H, FUKUHARA T. Multilayered snow/ice model on road[J]. Journal of Snow Engineering of Japan, 2008, 24(1):3-11.(in Japanese) [31] 藤本明宏, 渡邊洋, 福原輝幸.路面雪氷層-舗装間の接触熱抵抗[J].土木学会論文集, 2007(63):156-165. FUJIMOTO A, WATANABE H., FUKUHARA T. Thermal contact resisitance between snow/ice layer and pavement surface[J]. Journal of Japan Society of Civil Engineers, 2007(63):156-165.(in Japanese) [32] HEDSTROM N R, POMEROY J W. Measurements and modelling of snow interception in the boreal forest[J]. Hydrological Processes, 1998, 12(10/11):1611-1625. [33] GOODRICH L E. The influence of snow cover on the ground thermal regime[J]. Canadian Geotechnical Journal, 1982, 19(4):421-432. doi: 10.1139/t82-047 [34] LING F, ZHANG T J. Sensitivity of ground thermal regime and surface energy fluxes to tundra snow density in northern Alaska[J]. Cold Regions Science and Technology, 2006, 44(2):121-130. doi: 10.1016/j.coldregions.2005.09.002 [35] CHEN R S, LU S H, KANG E S, et al. A distributed water-heat coupled model for mountainous watershed of an inland river basin of Northwest China:Ⅰ:model structure and equations[J]. Environmental Geology, 2008, 53(6):1299-1309. doi: 10.1007/s00254-007-0738-2 [36] FOLEY J A, PRENTICE I C, RAMANKUTTY N, et al. An integrated biosphere model of land surface processes, terrestrial carbon balance, and vegetation dynamics[J]. Global Biogeochemical Cycles, 1996, 10(4):603-628. doi: 10.1029/96GB02692 [37] CHEN R S, KANG E S, LU S H, et al. A distributed water-heat coupled model for mountainous watershed of an inland river basin in Northwest China:Ⅱ:using meteorological and hydrological data[J]. Environmental Geology, 2008, 55(1):17-28. doi: 10.1007/s00254-007-0960-y -
点击查看大图
图(10) / 表 (1)
计量
- 文章访问数: 579
- HTML全文浏览量: 99
- PDF下载量: 108
- 被引次数: 0
