冰川冻土, 2021, 43(4): 948-963 doi: 10.7522/j.issn.1000-0240.2021.0054

第二次青藏高原综合科学考察研究

基于计算流体力学的寒区土壤水热耦合模型研究

胡锦华,1,2, 陆峥1,2, 仝金辉1,2, 李小雁1,2, 刘绍民1,2, 杨晓帆,1,2

1.北京师范大学 地理科学学部 地表过程与资源生态国家重点实验室, 北京 100875

2.北京师范大学 地理科学学部 自然资源学院, 北京 100875

Simulating thermo-hydrologic processes in cold region soil system： a computational fluid dynamics study

HU Jinhua,1,2, LU Zheng1,2, TONG Jinhui1,2, LI Xiaoyan1,2, LIU Shaomin1,2, YANG Xiaofan,1,2

1.State Key Laboratory of Earth Surface Processes and Resource Ecology，Faculty of Geographical Science，Beijing Normal University，Beijing 100875，China

2.School of Natural Resources，Faculty of Geographical Science，Beijing Normal University，Beijing 100875，China

通讯作者: 杨晓帆，教授，主要从事地下水科学、环境流体力学、多尺度建模与数值模拟研究. E-mail： xfyang@bnu.edu.cn

编委: 武俊杰

收稿日期: 2021-02-08   修回日期: 2021-06-20  

基金资助: 国家自然科学基金项目.  42077172.  41730854
第二次青藏高原综合科学考察研究项目.  2019QZKK0306

Received: 2021-02-08   Revised: 2021-06-20  

作者简介 About authors

胡锦华，博士研究生，主要从事寒区水文模拟研究.E-mail：hujinhua@mail.bnu.edu.cn , E-mail：hujinhua@mail.bnu.edu.cn

摘要

土壤水文过程（水分运移和传热）及其对气候变化的响应是寒区水文学的前沿问题。然而，冻土的存在使得寒区土壤水文过程变得极其复杂。此外，寒区自然环境恶劣，较难获取长时间序列和高分辨率的野外观测资料。近年来，充分利用已有的观测数据，构建寒区土壤水热耦合模型，并开展相应的数值模拟研究，已成为理解寒区土壤水文物理过程，揭示其动力学机制的重要途径。基于寒区土壤水文物理过程和计算流体力学方法，构建了高分辨率、适用于完全饱和状态下的寒区土壤水热耦合模型，且自主研发了相应的数值求解器和软件包。随后，通过一系列完全饱和状态下的验证算例，如经典的一维传热方程解析解、被广泛应用的二维基准测试算例和室内土柱冻结实验等，对已构建的模型进行了系统的检验。模型模拟结果与解析解、基准算例的结果以及实验数据相比，均有较好的一致性，表明该模型较为准确且高效地模拟了寒区土壤在完全饱和状态下的水分运移和传热过程，尤其能够精细刻画冻土水-冰相态变化等关键过程，有望成为研究寒区土壤水文过程的有力工具。

关键词： 寒区水文学 ; 地下水模型 ; 水热耦合过程 ; 计算流体力学 ; 基准测试算例

Abstract

Cold regions， which supply freshwater to downstream communities， have been significantly disturbed by anthropogenic climate change that is pronounced to continue in future. As unique and critical elements of the cold regions， frozen soil plays a vital role in cold region hydrology. Although frozen soil is considered relatively impermeable to groundwater flow， the freezing and thawing of the frozen soils under seasonal and climate change may influence the hydrologic components such as surface water infiltration， groundwater recharge， and hydrogeologic connectivity that are important for water resources. Under the pronounced climate warming， increased soil temperature expedites the thawing of frozen soil in cold regions. Therefore， understanding the complex thermo-hydrologic processes in frozen soils is one of key issues in studying cold region hydrology. However， due to the harsh environment， it is difficult to conduct hydrometeorological and geophysical field observations yet collect gap-free， high-resolution datasets. With the recent development of computational methods and resources and based on the observational data in recent decades， it is essential to enable modeling and simulation of the thermo-hydrologic processes in frozen soils to study cold region hydrology. The thermo-hydrologic processes in frozen soils is fundamentally interpreted into multi-phase flow and heat transfer in porous media with phase change. The current challenges in simulating such processes， especially in cold regions include but not limited to the interactions among liquid water， vapor， ice and soil grains and the phase change between ice and liquid water， which introduce high non-linearities and uncertainties in the processes. However， traditional distributed or hydrogeological models were limited in terms of resolution， accessibility and capability of solving the coupled thermo-hydrologic processes with high fidelity. To account for these feedbacks as far as possible， a physically-based， massively-parallel， fine-resolution and fully-saturated cryo-hydrogeological model was in-house developed in OpenFOAM® （http：//www.openfoam.com）， an open-source computational fluid dynamics （CFD） solver， which has been widely used in groundwater modeling and computational hydrology. The cryo-hydrogeological model， named darcyTHFOAM， was formulated by mass conservation， Darcy’s law， energy conservation and soil freezing function. To ensure the convenience and joy for users， an interface-based software was developed correspondingly using Python 3.5 （www.python.org） and designed with QT Designer （www.qt.io） in Linux environment. In general， rigorous validation and benchmarking of the numerical models are prerequisites and critical for their validation and applications. Hence， a series of benchmarking simulations under fully-saturated condition were performed to validate the current model， including the classical analytical solution of Stefan’s equation， two community-recognized benchmark cases （under thawing） using synthetic porous media samples in the Interfrost Project （wiki.lsce.ipsl.fr/interfrost） and laboratory freezing experiment with a 25 cm-long soil column. Simulated results include the evolution of temperature， liquid water content and ice content. Frozen depth， temperature profiles at fixed points or lines， the minimum of the temperature field were selected for quantitative comparisons. Simulated results were in excellent agreements with those obtained from the analytical solution， other cryo-hydrogeological models and experimental data， respectively， which demonstrated the reliability and accuracy of the current model. The study revealed that the proposed model is capable of simulating coupled thermo-hydrologic processes under fully-saturated condition in cold regions with high spatial resolutions and efficiency. Overall， the in-house developed cryo-hydrogeological model is expected to serve as a powerful tool for studying subsurface hydrology in cold regions. Further code developments involve coupling the surface processes， especially snow， and transport processes （e.g. contaminants） for studying groundwater-surface water interactions and hydro-biogeochemical processes in cold regions.

Keywords： cold region hydrology ; groundwater model ; coupled thermo-hydrologic processes ; computational fluid dynamics ; benchmark case

PDF (6222KB) 元数据 多维度评价 相关文章 导出 EndNote| Ris| Bibtex  收藏本文

本文引用格式

胡锦华, 陆峥, 仝金辉, 李小雁, 刘绍民, 杨晓帆. 基于计算流体力学的寒区土壤水热耦合模型研究[J]. 冰川冻土, 2021, 43(4): 948-963 doi:10.7522/j.issn.1000-0240.2021.0054

HU Jinhua, LU Zheng, TONG Jinhui, LI Xiaoyan, LIU Shaomin, YANG Xiaofan. Simulating thermo-hydrologic processes in cold region soil system： a computational fluid dynamics study[J]. Journal of Glaciology and Geocryology, 2021, 43(4): 948-963 doi:10.7522/j.issn.1000-0240.2021.0054

0 引言

冰冻圈是世界上许多大江大河的发源地，是干旱区内陆河流域的“水塔”［1-3］。冻土作为陆地冰冻圈的重要组成部分［4］，不仅对寒区水文过程有显著影响［5-7］，而且对全球气候变化极其敏感［8-10］。冻土本身具有相对隔水层的特点，使地下水文过程更加复杂［11-14］，进而影响地表水文过程、地-气之间的能量交换和碳循环［15-18］。而冻融循环对寒区水文过程（特别是地下部分）的影响，包括土壤水分的变化［19］、径流路径的变化［20-21］、地下水的分布［22-23］以及陆地水储量［24-26］等。此外，全球气候变暖导致土壤温度升高，多年冻土的退化加剧，季节冻结深度减小，融化深度增大；而地下存储的碳以温室气体的形式被释放到大气中［27-30］，又对全球变暖产生正反馈［21，31］。因此，研究寒区土壤水文过程，即土壤水热传输过程，是寒区水文学的关键问题之一。

由于寒区环境恶劣、人迹罕至，开展野外观测较为困难，因此，野外水文气象、水文地质和地球物理观测的范围、频率和数据量、数据精度等均受限［32-34］。随着数值模拟和计算机科学的快速发展，基于近几十年积累的观测数据，构建数值模型研究气候和季节变化下的寒区土壤水热传输过程，是理解寒区土壤水文物理过程，揭示其动力学机制的重要工具［35-37］。寒区土壤水热传输的物理过程较为复杂，涉及多相流动、冰水相态变化、冻胀等，使得描述水热传输过程的本构方程也具有高度的复杂性和非线性，这些都为建模和数值模拟带来了巨大的科学挑战［38］。目前，寒区水文模型主要包括基于陆地水文学的寒区分布式水文模型（distributed hydrologic model）［39］和基于水文地质学的寒区水文地质模型（groundwater/hydrogeologic model）［40］。寒区分布式水文模型，如CRHM［41］、GBEHM［42］、VIC［43］、HydroSiB2-SF［44］、WEB-DHM-pF［45］，是针对寒区水文循环中子过程/子单元分别建立子模型/模块，并集成至统一的模型框架中。这类模型已尽可能集成了寒区水文循环中的所有环节，但存在空间分辨率较低（千米级）、参数众多难以率定、需要大量观测数据加以验证的问题［46］，且侧重于陆面地表过程，对于地下水文过程仍停留在参数化方案或一维地下水模型层面［47］。因此，此类模型无法精细刻画寒区土壤水分运移和热传递的物理过程，对冻融等寒区特殊现象的本质成因和动力学机理研究不够成熟［48］。源自地下水科学、水文地质学和冻土水文学发展的寒区水文地质模型（cryo-hydrogeological model）［49］，如SUTRA-ICE［50-52］、FEFLOW［53-54］、ATS［55-56］、Cast3M（www-cast3m.cea.fr）、PFLOTRAN-ICE（www.pflotran.org）、HydroGeoSphere（HGS）［57-58］等，利用多孔介质渗流和传热理论建立地下水多相流动和热传递的本构方程，并通过土壤水分冻融曲线或基于物理过程的本构关系（例如Clausius-Clapyeron方程）来表征冻融循环对孔隙水相态、土壤孔隙率和渗透性的影响，进而模拟计算土壤水分和温度、地下径流量及其对流域地表径流的补给量等［59］。但是，此类模型较为复杂，大多为封装化或商业化软件，较难进行二次开发，且计算精度和效率有待提高。

综上所述，气候变化影响下的寒区土壤水文过程较为复杂，包括多相流动、传热和冰水相态变化等。虽然寒区分布式水文模型近年来发展较快，但多将土壤水文过程进行了不同程度的简化，且在模型建立、计算精度、并行计算效率等方面有待改进。此外，几类主流的寒区水文地质模型多为封装化或商业化软件，大大限制了模型的应用和推广。因此，亟需自主研发更稳定高效、高精度、支持并行计算、便于二次开发的开源寒区土壤水热耦合模型，精细刻画土壤温度、含水量与含冰量的时空动态变化，深入理解寒区土壤水文过程。此外，不同类别的寒区水文模型的研究对象和用途不尽相同，同一类别内的模型所选用的水文物理过程表征方式、参数化方案、数值离散方法和网格类型等也千差万别，其准确性和物理意义上的真实性均需校验［58］。本研究旨在基于寒区土壤水文物理过程和计算流体力学方法（computational fluid dynamics，CFD），利用大型并行开源软件OpenFOAM®（www.openfoam.com）自主研发一套高精度、高效率、适用于饱和状态下的寒区土壤水热耦合模型，用于描述饱和状态下的寒区土壤水分与温度的动态变化过程，并通过一维传热方程解析解、二维简单基准测试算例和室内冻结实验对模型进行系统的验证，综合评估其物理意义上的准确性与完整性、计算精度及计算效率。

1 寒区土壤水热耦合模型的构建

1.1　计算流体力学软件OpenFOAM简介

OpenFOAM®（https：//openfoam.org）是一个完全由C++编写的面向对象的计算流体力学（CFD）类库［60］（约100个），用于创建可执行文件（如应用程序），其数值模拟理念是将偏微分方程进行有限体积离散化后获得数值解。OpenFOAM自带了约250个应用程序，用户也可根据自己需求自行编写。应用程序主要分为两类：求解器与工具。其中，求解器用于解决特定的流体（或连续体）力学中的特定问题；工具用于执行涉及数据操作等任务。此外，OpenFOAM本身的工具包括前处理、后处理接口，以确保不同环境之间数据传输的协调性。OpenFOAM的整体结构如图1所示，包括核心的解算部分，以及前处理和后处理环境（如ParaView（www.paraview.org），Tecplot（www.tecplot.com）等）。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1   OpenFOAM总体结构示意图

Fig.1   Overall structure of the OpenFOAM


作为目前最强大的计算流体力学类库，OpenFOAM采用的数值离散方法为有限体积法（finite volume method，FVM），其自带的网格生成工具snappyHexMesh可以快速高效的划分六面体及多面体网格，因而可以处理复杂的几何外形，且所生成的网格质量较高；支持大型并行计算，计算效率高；研发环境良好，有利于对模型进行开发或二次开发；接口方式友好，便于耦合其他模型或求解器。因此，OpenFOAM已经被越来越多地用于地球科学中与环境流体力学相关的模拟计算［61-64］。

1.2　寒区土壤水热耦合模型

本研究针对寒区土壤水分运移和传热过程，建立了基于多相流动与传热理论的寒区土壤水热耦合模型，其控制方程包括质量守恒方程、动量方程和能量守恒定律。模型中包含的变量及参量定义汇总在表1中。

表1   模型参数

Table 1  List of model parameters

参数单位定义参数单位定义
λwW·m-1·K-1水的导热率Kwm·s-1导水率
λiW·m-1·K-1冰的导热率ρwkg·m-3水的密度
λsW·m-1·K-1土壤颗粒的导热率ρikg·m-3冰的密度
λtW·m-1·K-1总导热率ρskg·m-3土壤颗粒的密度
CwJ·kg-1·K-1水的比热容W冻结特征曲线的拟合参数
CiJ·kg-1·K-1冰的比热容Kintm2固有渗透系数
CsJ·kg-1·K-1土壤颗粒的比热容Ω阻抗系数
βm·s-2·kg-1多孔介质压缩系数θi固态含冰量
μkg·m-1·s-1水的动力黏度θw液态含水量
θr残余含水量TK温度
zm向上的垂直坐标hm压力水头
ε孔隙度PPa压力
LJ·kg-1融化潜热u⃗m·s-1速度
gm·s-2重力加速度ts时间

新窗口打开| 下载CSV


饱和条件下（即不考虑气相）的质量守恒方程［65］为

∂θwρw∂t+∇∙ρwu⃗=- ∂θiρi∂t

考虑水的可压缩性［49］，引入压缩系数β=1ρw∂ρw∂P，其中P为压力，则式（1）变为

θwβρw∂P∂t+ρw∂θw∂t+ρi∂θi∂t+∇∙ρwu⃗=0

若暂不考虑土壤的冻胀和融沉［66］，且在完全饱和情况下，θi=ε-θw，同时将地下水头定义为h=Pρwg，式（2）则变为

θwβρw2g∂h∂t+(ρw-ρi)∂θw∂t+∇∙ρwu⃗=0

假设地下水流动满足达西定律［67］，

u⃗=- Kw(∇h+∇z)

式中：Kw为导水率，定义为

Kw=KrTKintρwgμ

式中：Kr（T）为与土壤冻结相关的相对导水率，可表示［51，68］为

KrT=max(10- 6,10- ΩθiT)

式（6）描述了相对导水率Kr（T）与固态含冰量θi（T）之间的关系，原理是在冻结期由于土壤水成冰导致土壤孔隙变小，进而影响土壤水迁移［69］。此外，固态含冰量θiT=ε-θwT，而液态含水量θw是与温度T相关的变量θw（T），通常用经验函数来表示两者之间的关系，即土壤冻结特征曲线［70］。本模型中采用最为常用的指数形式的经验函数［49，51，64］来表示。

θwT=ε                                                     T≥273.15 Kθr+ε-θrexp- T-273.15W2    T<273.15 K

将式（4）代入式（3）中，并除以液态水的密度ρw，则最终的流动方程为

θwβg∂h∂t+(ρw-ρi)ρw∂θw∂t=∇∙(Kw(∇h+∇z))

描述传热过程的能量守恒方程［49］为

           ρwθwCw+ρiθiCi+1-ερsCs-ρiL∂θi∂T∂T∂t =∇∙λt∇T+ρwCw∇∙u⃗T

式中：λt为土壤总导热率，可表示为

λt=θwλw+θiλi+(1-ε)λs

由此可见，描述寒区地下水热传输过程的流动和传热方程，由于多相流和冻融引起的水冰相态变化，具有高度的非线性。

目前，该模型已被成功编译至开源CFD软件OpenFOAM®（v1712版本）中，成为一类新的求解器，命名为darcyTHFOAM。模型的具体计算流程图如图2所示，在每个时间步长中，首先利用OpenFOAM中不完全的对角Cholesky预条件共轭求解器（DIC-PCG）［60］来求解流动方程（8），更新速度场（4）及土壤的水热性质，进而利用对角不完全LU稳定化预条件双共轭求解器（DILC-PBiCGStab）［60］来求解传热方程（9），用Picard循环来处理方程中的非线性问题。其中，Δt为时间步长，NPicP 为求解压力的Picard迭代次数，NIterP 为求解压力的最大Picard迭代次数，P为压力，ErrP为求解压力的误差，PicP为求解压力的最大误差， U 为速度，T为温度，NPicT 为求解温度的Picard迭代次数，NIterT 为求解温度的最大Picard迭代次数，ErrT为求解温度的误差，PicT为求解温度的最大误差，tfactor 为自动调整时间步长的参数。每个时间步长∆t中，在求解压力的流动方程（8）时，如果误差ErrP超过最大误差PicP，则在Picard最大迭代次数NIterP 内继续求解，直至误差ErrP低于最大误差PicP，进而求解速度；同样的，在求解传热方程（9）时，如果误差ErrT超过最大误差PicT，则在Picard最大迭代次数NPicT 内继续求解，直至误差ErrT低于最大误差PicT；如果求解温度与压力的误差均满足最大误差的限制，则继续下一时间步长的计算，否则，自动调整时间步长，重新展开计算，直至满足误差需求。另外，darcyTHFOAM的主要性能特点如表2所示，适用于Linux系统，编程语言为C++，时间离散方法为欧拉格式，空间离散方法为线性插值格式，采用自适应的时间步长策略来平衡计算精度和计算效率［71-74］，可在本机上或高性能计算集群（例如，中国广州的国家超级计算中心Tianhe-2A-TH-IVB-FEP集群，www.nscc-gz.cn）完成仿真模拟计算。

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2   darcyTHFOAM计算流程图

Fig.2   Flow chart of the darcyTHFOAM


表2   darcyTHFOAM的性能特点

Table 2  Specific features of the darcyTHFOAM

求解器darcyTHFOAM
编程语言C++
操作系统Linux
计算机x86
开发平台OpenFOAM®（version v1712）
数值方法有限体积法
时间步长策略自适应
非线性Picard循环
线性求解器PCG/PBiCG
预条件DIC/DILU
自动划分网格否
并行计算支持

新窗口打开| 下载CSV


总体来说，darcyTHFOAM模型的基本假设为：①不考虑气相，含水层处于完全饱和状态下；②不考虑土壤的冻胀和融沉；③土壤水流动满足达西定律；④采用指数形式的土壤冻结特征曲线，来表示液态含水量与温度的关系。而darcyTHFOAM模型是在开源的计算流体力学软件OpenFOAM的基础上开发的，故支持高效的并行计算（适用于mm至km级的模拟计算）；时间步长采用自适应的时间步长策略（可从μs到h）；OpenFOAM采用六面体来划分网格，空间分辨率可精细至mm级；网格本身是三维，但通过调整各个维度的长度，进而可使模型适用于一维、二维、三维的完全饱和算例。该模型适用于含水层完全饱和状态下的冻结与融化过程，如饱和状态下的基准测试算例研究（冰包裹体融化算例、融区贯通/闭合算例）、饱和状态下的室内土柱冻结实验等等，但是，darcyTHFOAM模型的应用范围不仅限于此，其适用于一维到三维的问题，时空分辨率较广（从mm到km，从μs到h），边界条件灵活可控，并且，除了土壤含水层外，还可用于其他多孔介质。此外，darcyTHFOAM模型接口非常灵活，用户可根据需求耦合新的方程（例如溶质运移）、引入复杂的边界条件与新的源项等。最后，darcyTHFOAM模型也可根据需求在现有的本构方程基础上进行修改，使其适用于非饱和条件。

1.3　软件

为方便用户使用，利用Python 3.5（www.python.org）和Qt设计器（www.qt.io）设计了对应的可视化界面软件。该软件支持Linux系统，包括8个主界面：登录、打开算例文件、划分网格、设置初始值、设置水热参数、设置特殊区域、求解、可视化，如图3所示。其中，用户输入正确的用户密码方可登陆，在“划分网格”的界面，定义算例的边界条件类型和计算区域；压力场、速度场和温度场的初始值、在特殊区域的值，分别在“设置初始值”和“设置特殊区域”中定义；土壤中各相的水热性质、土壤冻结特征曲线、相对导水率函数均在“设置水热参数”中定义；最后，用户可以在“求解”中自定义时间步长、模拟时间、输出时间间隔，所有的仿真模拟计算将在后台完成，而结果展示与分析将通过调用ParaView实现。

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3   darcyTHFOAM软件界面

Fig.3   Interfaces of the darcyTHFOAM software： login （a）， open folder （b）， create mesh （c）， initial value （d）， basic properties （e）， set fields （f）， solve （g）， and visualization （h）


2 寒区土壤水热耦合模型的验证

在模型研究中，选择合适稳定且具有物理意义的基准测试算例来验证，是非常必要的。2014年底，为全面验证和优化寒区地下水文模型，由美国、加拿大、瑞典、德国、英国、荷兰和法国等国的科学家组成了科学小组，并启动了“寒区地下水热耦合模型比较计划（InterFrost）”（wiki.lsce.ipsl.fr/interfrost）。InterFrost计划的主旨是提出一系列从简单到复杂，且具有代表性和实用性的基准测试算例，用于模型比较和基准验证，进而优化和发展模型。因此，在本研究中采用一维解析解、InterFrost的两个基准测试算例（冰包裹体融化和融区贯穿/闭合）以及室内冻结实验对darcyTHFOAM模型逐步进行验证。

2.1　一维Stefan传热方程

一维Stefan方程存在解析解［34，75-76］，可根据温度来估算冻结深度随时间的变化。

Zf=(2ρiθiL)12·[t∑(λuT)]12

式中：Zf为冻结深度；t为时间；λu为热导率。

为更好地与Stefan方程的解析解对比，Schilling等［34］采用HGS模拟了一维垂直土柱的融化过程。根据他们研究中的算例设置，假定土柱是均质的，高10 m，起初处于-5 ℃的完全冻结状态，土柱的顶部和底部为恒定为10 ℃和-5 ℃。土壤固体颗粒的导热率λs=0.017 J·s-1·cm-3·K-1，比热容Cs=35.07 J·cm-3·K-1，孔隙度ε=0.4。为将模拟结果与解析解对比，将土壤冻结特征曲线中的拟合参数W设置为4，其他参数与以下基准测试算例相同。

整个计算区域在垂直方向上划分了400个网格，起始步长和最大步长分别为0.1 s和30 min，利用自主开发的darcyTHFOAM模拟计算了90天的土柱融化过程，在本机（Window系统/4G RAM/1T内存）下的虚拟机（Ubuntu18.04系统/2G RAM/40G内存）（下同）上运行时间约为48 s。土柱中冻结锋面与土柱顶部的距离变化，最能反映土柱的融化过程，且能够与解析解比较。如图4所示，darcyTHFOAM模拟的冻结深度与Stefan方程的解析解和HGS的模拟结果相比，结果均非常准确。

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4   darcyTHFOAM模拟的冻结深度与Stefan方程的解析解和HGS模拟结果的相互比较

Fig.4   Inter-comparison of the frost depth calculated using Stefan’s equation and simulated by the HGS and darcyTHFOAM


2.2　冰包裹体融化算例

如图5所示，整个计算区域长为3 m，宽为1 m，在内部有一边长为0.333 m的正方形冰块。左侧入口处设置为恒定的5 ℃，其他边界均为绝热边界，除了蓝色冰块为-5 ℃，其他区域温度均为5 ℃。而针对压力水头边界条件，在左侧入口与右侧出口边界分别设置为固定值0.09 m、0 m，上下边界均为无通量边界，内部初始值与右侧出口边界的值相同，其他相关参数见表3。在压力和温度的驱动下，初始处于冻结状态的冰块将慢慢融化。

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5   冰包裹体融化算例的计算区域及边界条件［49］（灰色字为流动边界条件，黑色字为传热边界条件）

Fig.5   Melting of ice inclusion case： computational domain， initial and boundary conditions［49］ （gray characters： flow， black characters： heat transfer）


表3   两个基准测试算例中的参数设置

Table 3  Parameters used in two benchmark cases

参数数值单位参数数值单位
λw0.60W·m-1·K-1Kint1.3×10-10m2
λi2.14W·m-1·K-1L334 000J·kg-1
λs9.00W·m-1·K-1g9.81m·s-2
Cw4 182J·kg-1·K-1ρw1 000kg·m-3
Ci2 060J·kg-1·K-1ρi920kg·m-3
Cs835J·kg-1·K-1ρs2 650kg·m-3
β1.0×10-8m·s-2·kg-1ε0.37
μ1.793×10-3kg·m-1·s-1W0.50
Ω50θr0.0185

新窗口打开| 下载CSV


在通过网格分辨率测试后，整个计算区域最终被划分为300×100个网格，起始时间步长为0.1 s，随后在计算过程中可自动调整。每个时间步长，压力和温度的收敛标准均设置为1×10-10。在求解压力的Picard循环中，最大迭代次数和最大误差分别设置为20和1×10-3，而在求解温度的Picard循环中，最大迭代次数和最大误差分别设置为50和1×10-5。为提高计算效率，将计算区域划分为120个部分，在天河二号超级计算器上调用了5个核（共120个节点）进行仿真模拟计算，最终用时3 min完成了5天的冰块融化过程。

图6和图7分别展示了温度和液态水饱和度（液态水含量/孔隙度）随时间的变化，进而反映了冰块在不同时刻的融化情况。在冰块融化过程中，沿计算区域中心线的液态水饱和度如图8（a）所示，从上到下是从初始到最终时刻的液态水饱和度，最后全部融化，与顶部坐标轴完全重合；沿计算区域中心线的温度如图8（b）所示，从下到上表示初始和最终时刻的温度，黑色箭头表示温度剖面的变化。在压力水头梯度的驱动下，融化加速，较冷的流体通过平流和热扩散向下游输送，最终达到与初始条件相同的温度（5 ℃）。为了更好地将温度的时空变化与Interfrost中给出的其他模型的结果进行量化对比，参照Grenier等［49］中的评价指标，选取了计算区域的最低温度进行对比。此外，其他模型或软件的数值方案及网格划分情况已经总结在表4中。结果表明，darcyTHFOAM模拟的最低温度与其他模型预测的最低温度一致［图8（c）］。

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6   冰包裹体融化算例在不同时刻的温度演变

Fig.6   Melting of ice inclusion case： evolution of temperature at t=7 200 s （a）， t=21 600 s （b）， t=43 200 s （c）， and t=64 800 s （d）


图7

新窗口打开| 下载原图ZIP| 生成PPT

图7   冰包裹体融化算例在不同时刻的液态水饱和度演变

Fig.7   Melting of ice inclusion case： evolution of liquid water saturation at t=7 200 s （a）， t=21 600 s （b）， t=43 200 s （c）， and t=64 800 s （d）


图8

新窗口打开| 下载原图ZIP| 生成PPT

图8   冰包裹体融化算例沿中轴线的演变过程及对比验证

Fig.8   Melting of ice inclusion case： simulated instantaneous liquid water saturation profiles along the central line of the computational domain （profiles from bottom to top： initial and final moments， the line coincide with the top axis in the end） （a）， simulated instantaneous temperature profiles along the central line of the computational domain （profiles from bottom to top： temperature at initial and final moments， the black arrow indicates the revolution of the profiles） （b）， and inter-comparison of the minimum of the temperature field （c）


表4   darcyTHFOAM与其他模型或软件数值算法和网格离散方法的比较

Table 4  Mesh and numerical schemes used by darcyTHFOAM， Cast3M， permaFOAM， COMSOL and DarcyTools

模型darcyTHFOAMCast3MpermaFOAMCOMSOLDarcyTools
数值算法有限体积法有限体积法有限体积法有限元法有限体积法
网格数量30 00031 609480 00035 00014 786
网格类型六面体四边形六面体四边形笛卡尔网格
计算区域全部一半全部全部全部
并行计算

支持，

1~1 000核

不支持

支持，

100~1 000核

共享内存，

8核

支持，

64核

时间步长

策略

自适应调整

时间步长

指定时间

步长

自适应调整

时间步长

自适应调整

时间步长

指定时间

步长

文献来源本研究www-cast3m.cea.fr［64］www.comsol.com/comsol-multiphysics［81］

新窗口打开| 下载CSV


然而，不同模型在温度小于0 ℃的阶段呈现出一定差异，特别是在-1~0 ℃（对应相变阶段），温度呈现出急速上升趋势的开始时刻，darcyTHFOAM，分别比COMSOL和darcyTools晚2 400 s和1 300 s，而又比Cast3M和permaFOAM提早3 500 s和5 400 s，这可能是由于不同的数值算法和网格离散方法造成的。与同类模型相比，darcyTHFOAM采用更少的网格数量便可达到相同的模拟结果（如表4所示），且支持并行计算，计算效率更高，接口较为灵活，便于二次开发。但目前仅适用于完全饱和状态下的水热耦合过程，未来可拓展至模拟变饱和状态下的土壤水热传输过程。

2.3　融区贯穿/闭合算例

贯穿融区（“天窗”）存在于湖泊或河流下面的多年冻土中，是多年冻土内的局部融化区［77-80］。该基准测试算例取材于自然界中冻土的典型特征，兼具代表性与适用性［49］，能够反映出多年冻土中贯穿融区的演变过程。两个起始冻结的半圆形区域（半径为0.5099 m）温度都是-5 ℃，代表了两部分冻土，均位于一个正方形计算域内（1 m×1 m）。最初，计算区域的背景温度为5 ℃，左侧入口边界赋予固定值5 ℃，上下边界赋予固定值-5 ℃（图9），右侧出口边界为绝热边界，压力边界条件与上一个冰包裹体融化算例的相同，其他相关参数见表3。

图9

新窗口打开| 下载原图ZIP| 生成PPT

图9   融区贯穿/闭合算例的计算区域及边界条件［49］（灰色字为流动边界条件，黑色字为传热边界条件）

Fig.9   Talik opening/closure case： computational domain， initial and boundary conditions［49］ （gray characters： flow， black characters： heat transfer）


经过网格分辨率测试，整个计算区域最终被划分为100×100个网格，起始时间步长为0.1 s，在计算的过程可自动调整，最大时间步长设置为1 000 s。对于每个时间步长，压力和温度的收敛标准均设置为1×10-10。在求解压力的Picard循环中，最大迭代次数和最大误差分别设置为20和1×10-3，而在求解温度的Picard循环中，最大迭代次数和最大误差分别设置为50和1×10-5。为提高计算效率，将计算区域分为了4部分，在天河二号超级计算器上采用1个核（共4个节点）进行仿真模拟计算，最终耗费2 min完成了40 h的融区演变过程模拟。

不同时刻的温度和液态水饱和度分布情况如图10和图11所示，计算区域中心垂直线上的液态水饱和度和温度变化如图12（a）和图12（b）所示，从上到下依次代表，从初始到终止时刻，其中，黑色箭头表示温度剖面的旋转。在温度与压力的驱动下，两个最初冻结的半圆区域逐渐融化。为了将darcyTHFOAM的模拟结果与Grenier等［49］中其他求解器的模拟数据进行定量比较，在图12（c）中绘制了两个选定位置（Pt1点和Pt2点）的温度曲线。这两个点都接近冻结区和非冻结区之间的初始边界，对融区的贯通/闭合过程非常敏感。结果表明，这两点处的温度变化与其他模型模拟的结果非常一致，再次证明了模型的可靠性。

图10

新窗口打开| 下载原图ZIP| 生成PPT

图10   融区贯穿/闭合算例在不同时刻的温度演变

Fig.10   Talik opening/closure case： evolution of temperature at t=7 200 s （a）， t=21 600 s （b）， and t=86 400 s （c）


图11

新窗口打开| 下载原图ZIP| 生成PPT

图11   融区贯穿/闭合算例在不同时刻的液态水饱和度演变

Fig.11   Talik opening/closure case： evolution of liquid water saturation and with an imposed pressure head gradient of 3% at t=7 200 s （a）， t=21 600 s （b）， and t=86 400 s （c）


图12

新窗口打开| 下载原图ZIP| 生成PPT

图12   融区贯通/闭合算例沿中轴线的演变过程及对比验证

Fig.12   Talik opening/closure case： simulated instantaneous liquid water profiles along the central line of the computational domain （profiles from top to bottom： initial and final moments） （a）， simulated instantaneous temperature profiles along the central line of the computational domain （profiles from top to bottom： temperature at initial and final moments， the black arrow indicates the revolution of the profiles） （b）， and temperature profiles at point Pt1 and point Pt2 （c）


2.4　室内土柱冻结实验

Mizoguchi［82］早在1991年就开展了室内土柱冻结实验，并在实验过程中进行了温度观测。随着数值模型的发展，已有越来越多的寒区水文地质模型和室内实验进行比对［83］，因此将darcyTHFOAM模拟此室内冻结实验的结果与实测的温度数据比对，是模型验证的有效途径。

室内土柱冻结实验是在石英砂填充的饱和环境中进行的，为方便比较，在数值模拟时采用轴对称区域，初始和边界条件如图13所示，整个土柱的初始温度为5 ℃，两侧均为绝热边界，而顶部和底部分别暴露于温度为-10 ℃和5 ℃的循环流体。针对压力边界条件，在土柱的底部、顶部和两侧均采用零通量边界条件。因此，在温度梯度的驱动下，土柱将从上到下开始冻结。

图13

新窗口打开| 下载原图ZIP| 生成PPT

图13   室内冻结实验计算域设置（灰色字为流动边界条件，黑色字为传热边界条件）

Fig.13   Laboratory freezing experiment： computational domain， initial and boundary conditions （gray characters： flow， black characters： heat transfer）


表5中总结了实验参数，其中，土壤固体颗粒Cs的导热系数是从Mochizuki等［84］中获得，相对导水率函数中的阻抗系数Ω，参考前面基准测试算例设置为50。通过对比模拟结果和实验数据，将土壤冻结函数中的拟合参数W设为5.25。考虑到计算效率和精度，经过网格分辨率测试，将计算域离散为80×500个网格。该算例的收敛准则和Picard循环的设置，也与基准测试算例相同。初始时间步长设置为0.1 s，最大时间步长设置为30 min。采用darcyTHFOAM在本机上模拟计算50 h的冻结过程，耗时约3 min（CPU时间）。

表5   室内冻结实验中的参数

Table 5  Parameters used in simulating laboratory freezing experiment

参数数值单位参数数值单位
λw0.56W·m-1·K-1Kint1.3×10-10m2
λi2.24W·m-1·K-1L334 000J·kg-1
λs1.90W·m-1·K-1g9.81m·s-2
Cw4 180J·kg-1·K-1ρw1 000kg·m-3
Ci2 050J·kg-1·K-1ρi920kg·m-3
Cs835J·kg-1·K-1ρs2 650kg·m-3
β1.0×10-8m·s-2·kg-1ε0.4025
μ1.793×10-3kg·m-1·s-1θr0.02
Ω50W5.25

新窗口打开| 下载CSV


利用darcyTHFOAM展开仿真模拟计算后，1、3、6、12、24 h的温度和液态水饱和度分布情况分别如图14（a）和图14（b）所示。从顶部到底部，温度逐渐下降，同时液态水饱和度逐渐减少，直观反映出了土柱的冻结过程。为了进一步验证模型，从文献［82］的图3.1-5中提取出实验观测的温度数据（Engauge Digitize，http：//digitizer.sourceforge.net）。图15中展示了在不同时刻沿对称轴的温度分布，darcyTHFOAM模拟的数值结果与实验数据非常吻合，进一步证明了模型的准确性。Mizoguchi［82］在设计该室内冻结实验时，整个土柱的初始温度为5 ℃，而顶部和底部分别暴露于温度为-10 ℃和5 ℃的循环流体。但从初始0 h的观测数据来看，整个土柱内部温度平均在7 ℃，只有顶部温度为5 ℃，说明在土柱冻结实验开始时内部并没有完全达到5 ℃，进而导致观测数据在早期（1 h）呈现出高估的趋势（高达6 ℃）。

图14

新窗口打开| 下载原图ZIP| 生成PPT

图14   在冻结开始后1、3、6、12、24小时的温度（a）和液态水饱和度（b）的演变

Fig.14   Evolution of temperature （a） and liquid water saturation （b） at 1， 3， 6， 12 and 24 h after freezing started


图15

新窗口打开| 下载原图ZIP| 生成PPT

图15   室内冻结实验数据与数值模拟结果对比

Fig.15   Comparison between experimental and numerical results of temperature distribution along the symmetry axis at 0， 1， 3， 6， 12 and 24 h after freezing started


3 结论与展望

本研究基于多孔介质渗流和传热理论、寒区土壤水文物理过程和计算流体力学方法，建立了高精度、高效率的寒区土壤水热耦合模型（darcyTHFOAM），综合使用Stefan方程解析解、基准测试算例（冰包裹体融化、融区贯穿/闭合）和室内实验对模型进行了系统的校验。该模型能够精细刻画寒区土壤中水分迁移和热量传输过程，尤其是冰水相态的动态变化，为探究寒区土壤水热传输的物理过程、动力学机理和影响机制提供强大的数值模拟工具，有望为气候变化引起的冻土退化、季节冻土变化和水资源短缺等地球系统科学问题提供理论依据和可靠预测。此外，模型的研发环境较为友好，具有一定的鲁棒性，方便二次开发，未来将继续耦合地表过程（如积雪消融）和生物地球化学反应动力学过程，以研究寒区地下水-地表水相互作用和生物地球化学循环。

参考文献

Qin Dahe， Yao Tandong， Ding Yongjian， et al. Introduction to cryospheric science［M］. Beijing： Science Press， 2017.

[本文引用: 1]

秦大河， 姚檀栋， 丁永建， 等. 冰冻圈科学概论［M］. 北京： 科学出版社， 2017.

[本文引用: 1]

Yao Tandong， Thompson L， Yang Wei， et al.

Different glacier status with atmospheric circulations in Tibetan Plateau and surroundings

［J］. Nature Climate Change， 2012， 2（9）： 663-667.

Immerzeel W W， Lutz A F， Andrade M， et al.

Importance and vulnerability of the world’s water towers

［J］. Nature， 2020， 577（7790）： 364-369.

[本文引用: 1]

Ming-ko Woo. Permafrost hydrology［M］. Berlin： Springer， 2012.

[本文引用: 1]

Gao Hongkai， Wang Jingjing， Yang Yuzhong， et al.

Permafrost hydrology of the Qinghai-Tibet Plateau： a review of processes and modeling

［J/OL］. Frontiers in Earth Science， 2021， 8 ［2021-06-06］. .

[本文引用: 1]

Ding Yongjian， Zhang Shiqiang， Chen Rensheng， et al.

Hydrological basis and discipline system of cryohydrology： from a perspective of cryospheric science

［J/OL］. Frontiers in Earth Science， 2020， 8 ［2021-06-06］. .

Wang Genxu， Lin Shan， Hu Zhaoyong， et al.

Improving actual evapotranspiration estimation integrating energy consumption for ice phase change across the Tibetan Plateau

［J/OL］. Journal of Geophysical Research： Atmospheres， 2020， 125 ［2021-06-06］. .

[本文引用: 1]

Kang Shichang， Guo Wanqin， Zhong Xinyue， et al.

Changes in the mountain cryosphere and their impacts and adaptation measures

［J］. Climate Change Research， 2020， 16（2）： 143-152.

[本文引用: 1]

康世昌， 郭万钦， 钟歆玥， 等.

全球山地冰冻圈变化、影响与适应

［J］. 气候变化研究进展， 2020， 16（2）： 143-152.

[本文引用: 1]

IPCC.

Summary for policymakers

［R/OL］ // IPCC special report on the ocean and cryosphere in a changing climate. Geneva， Switzerland： IPCC， 2019 ［2021-06-06］. .

Ding Yongjian， Mu Cuicui， Wu Tonghua， et al.

Increasing cryospheric hazards in a warming climate

［J/OL］. Earth-Science Reviews， 2020， 213 ［2021-06-06］. .

[本文引用: 1]

McKenzie J M， Kurylyk B L， Walvoord M A， et al.

Invited perspective： what lies beneath a changing arctic？

［J］. The Cryosphere， 2021， 15（1）： 479-484.

[本文引用: 1]

Lemieux J M， Fortier R， Murray R， et al.

Groundwater dynamics within a watershed in the discontinuous permafrost zone near Umiujaq （Nunavik， Canada）

［J］. Hydrogeology Journal， 2020， 28（3）： 833-851.

Chen Rensheng， Wang Guanxing， Yang Yong， et al.

Effects of cryospheric change on alpine hydrology： combining a model with observations in the upper reaches of the Hei River， China

［J］. Journal of Geophysical Research： Atmospheres， 2018， 123： 3414-3442.

Wang Genxu， Mao Tianxu， Chang Juan， et al.

Processes of runoff generation operating during the spring and autumn seasons in a permafrost catchment on semi-arid plateaus

［J］. Journal of Hydrology， 2017， 550： 307-317.

[本文引用: 1]

Nitze I， Grosse G， Jones B M， et al.

Remote sensing quantifies widespread abundance of permafrost region disturbances across the Arctic and Subarctic

［J/OL］. Nature Communications， 2018， 9（1） ［2021-06-06］. .

[本文引用: 1]

Chang Juan， Wang Genxu， Guo Linmao.

Simulation of soil thermal dynamics using an artificial neural network model for a permafrost alpine meadow on the Qinghai-Tibetan Plateau

［J］. Permafrost and Periglacial Processes， 2019， 30（3）： 195-207.

Wang Qingfeng， Zhang Tingjun， Wu Jichun， et al.

Investigation on permafrost distribution over the upper reaches of the Heihe River in the Qilian Mountains

［J］. Journal of Glaciology and Geocryology， 2013， 35（1）： 19-29.

王庆峰， 张廷军， 吴吉春， 等.

祁连山区黑河上游多年冻土分布考察

［J］. 冰川冻土， 2013， 35（1）： 19-29.

Zhao Lin， Zou Defu， Hu Guojie， et al.

Changing climate and the permafrost environment on the Qinghai-Tibet （Xizang） Plateau

［J］. Permafrost and Periglacial Processes， 2020， 31（3）： 395-405.

[本文引用: 1]

Zhao Lin， Sheng Yu. Permafrost changes on Qinghai-Tibet Plateau［M］. Beijing： Science Press， 2019.

[本文引用: 1]

赵林， 盛煜. 青藏高原多年冻土及变化［M］. 北京： 科学出版社， 2019.

[本文引用: 1]

Cheng Guodong， Zhao Lin， Li Ren， et al.

Characteristic， changes and impacts of permafrost on Qinghai-Tibet Plateau

［J］. Chinese Science Bulletin， 2019， 64（27）： 2783-2795.

[本文引用: 1]

程国栋， 赵林， 李韧， 等.

青藏高原多年冻土特征、变化及影响

［J］. 科学通报， 2019， 64（27）： 2783-2795.

[本文引用: 1]

Kurylyk B L， Macquarrie K T B， McKenzie J M.

Climate change impacts on groundwater and soil temperatures in cold and temperate regions： implications， mathematical theory， and emerging simulation tools

［J］. Earth-Science Reviews， 2014， 138： 313-334.

[本文引用: 2]

Bense V F， Ferguson G， Kooi H.

Evolution of shallow groundwater flow systems in areas of degrading permafrost

［J/OL］. Geophysical Research Letters， 2009， 36（22） ［2021-06-06］. .

[本文引用: 1]

Walvoord M A， Kurylyk B L.

Hydrologic impacts of thawing permafrost： a review

［J/OL］. Vadose Zone Journal， 2016， 15（6） ［2021-06-06］. .

[本文引用: 1]

St.

Jacques J-M

SauchynD J. Increasing winter baseflow and mean annual streamflow from possible permafrost thawing in the Northwest Territories， Canada［J/OL］. Geophysical Research Letters， 2009， 36（1） ［2021-06-06］. .

[本文引用: 1]

Rennermalm A K， Wood E F， Troy T J.

Observed changes in pan-arctic cold-season minimum monthly river discharge

［J］. Climate Dynamics， 2010， 35（6）： 923-939.

Somers L D， McKenzie J M.

A review of groundwater in high mountain environments

［J/OL］. WIREs Water， 2020， 7（6） ［2021-06-06］. .

[本文引用: 1]

Turetsky M R， Abbott B W， Jones M C， et al.

Carbon release through abrupt permafrost thaw

［J］. Nature Geoscience， 2020， 13（2）： 138-143.

[本文引用: 1]

Wang Taihua， Yang Dawen， Yang Yuting， et al.

Permafrost thawing puts the frozen carbon at risk over the Tibetan Plateau

［J/OL］. Science Advances， 2020， 6（19） ［2021-06-06］. .

Song Xiaoyan， Wang Genxu， Fan Fei， et al.

Soil moisture as a key factor in carbon release from thawing permafrost in a boreal forest

［J/OL］. Geoderma， 2020， 357 ［2021-06-06］. .

Walvoord M A， Voss C I， Ebel B A， et al.

Development of perennial thaw zones in boreal hillslopes enhances potential mobilization of permafrost carbon

［J/OL］. Environmental Research Letters， 2019， 14（1） ［2021-06-06］. .

[本文引用: 1]

Karra S， Painter S L， Lichtner P C.

Three-phase numerical model for subsurface hydrology in permafrost-affected regions （PFLOTRAN-ICE v1.0）

［J］. The Cryosphere， 2014， 8（5）： 1935-1950.

[本文引用: 1]

Burt T P， Mcdonnell J J.

Whither field hydrology？ The need for discovery science and outrageous hydrological hypotheses

［J］. Water Resources Research， 2015， 51（8）： 5919-5928.

[本文引用: 1]

Che Tao， Li Xin， Liu Shaomin， et al.

Integrated hydrometeorological， snow and frozen-ground observations in the alpine region of the Heihe River Basin， China

［J］. Earth System Science Data， 2019， 11（3）： 1483-1499.

Schilling O S， Park Y， Therrien R， et al.

Integrated surface and subsurface hydrological modeling with snowmelt and pore water freeze-thaw

［J］. Groundwater， 2019， 57（1）： 63-74.

[本文引用: 3]

Ding Yongjian， Zhang Shiqiang， Chen Rensheng. Introduction to hydrology in cold regions［M］. Beijing： Science Press， 2017.

[本文引用: 1]

丁永建， 张世强， 陈仁升. 寒区水文导论［M］. 北京： 科学出版社， 2017.

[本文引用: 1]

Yang Meixue， Wang Xuejia， Pang Guojin， et al.

The Tibetan Plateau cryosphere： observations and model simulations for current status and recent changes

［J］. Earth-Science Reviews， 2019， 190： 353-369.

Bui M T， Lu Jinmei， Nie Linmei.

A review of hydrological models applied in the permafrost-dominated arctic region

［J/OL］. Geosciences， 2020， 10（10） ［2021-06-06］. .

[本文引用: 1]

Lamontagne-Hallé P， McKenzie J M， Kurylyk B L， et al.

Guidelines for cold-regions groundwater numerical modeling

［J/OL］. WIREs Water， 2020， 7（6） ［2021-06-06］. .

[本文引用: 1]

Yang Dawen， Gao Bing， Jiao Yang， et al.

A distributed scheme developed for eco-hydrological modeling in the upper Heihe River

［J］. Science China： Earth Sciences， 2015， 58（1）： 36-45.

[本文引用: 1]

Tian Yong， Zheng Yi， Han Feng， et al.

A comprehensive graphical modeling platform designed for integrated hydrological simulation

［J］. Environmental Modelling & Software， 2018， 108： 154-173.

[本文引用: 1]

Marsh C B， Pomeroy J W， Wheater H S.

The Canadian Hydrological Model （CHM） v1.0： a multi-scale， multi-extent， variable-complexity hydrological model-design and overview

［J］. Geoscientific Model Development， 2020， 13（1）： 225-247.

[本文引用: 1]

Li Hongyi， Li Xin， Yang Dawen， et al.

Tracing snowmelt paths in an integrated hydrological model for understanding seasonal snowmelt contribution at basin scale

［J］. Journal of Geophysical Research： Atmospheres， 2019， 124： 8874-8895.

[本文引用: 1]

Hamman J J， Nijssen B， Bohn T J， et al.

The Variable Infiltration Capacity model version 5 （VIC-5）： infrastructure improvements for new applications and reproducibility

［J］. Geoscientific Model Development， 2018， 11： 3481-3496.

[本文引用: 1]

Wang Lei， Zhou Jing， Qi Jia， et al.

Development of a land surface model with coupled snow and frozen soil physics

［J］. Water Resources Research， 2017， 53（6）： 5085-5103.

[本文引用: 1]

Song Lei， Wang Lei， Li Xiuping， et al.

Improving permafrost physics in a distributed cryosphere-hydrology model and its evaluations at the upper Yellow River basin

［J/OL］. Journal of Geophysical Research： Atmospheres， 2020， 125（18） ［2021-06-06］. .

[本文引用: 1]

Wang Lei， Li Xiuping， Zhou Jing， et al.

Hydrological modelling over the Tibetan Plateau： current status and perspective

［J］. Advances in Earth Science， 2014， 29（6）： 674-682.

[本文引用: 1]

王磊， 李秀萍， 周璟， 等.

青藏高原水文模拟的现状及未来

［J］. 地球科学进展， 2014， 29（6）： 674-682.

[本文引用: 1]

Yu Lianyu， Zeng Yijian， Su Zhongbo.

Understanding the mass， momentum， and energy transfer in the frozen soil with three levels of model complexities

［J］. Hydrology and Earth System Sciences， 2020， 24（10）： 4813-4830.

[本文引用: 1]

Ye Renzheng， Chang Juan.

Study of groundwater in permafrost regions of China： status and process

［J］. Journal of Glaciology and Geocryology， 2019， 41（1）： 183-196.

[本文引用: 1]

叶仁政， 常娟.

中国冻土地下水研究现状与进展综述

［J］. 冰川冻土， 2019， 41（1）： 183-196.

[本文引用: 1]

Grenier C， Anbergen H， Bense V， et al.

Groundwater flow and heat transport for systems undergoing freeze-thaw： intercomparison of numerical simulators for 2D test cases

［J］. Advances in Water Resources， 2018， 114： 196-218.

[本文引用: 11]

Evans S G， Godsey S E， Rushlow C R， et al.

Water tracks enhance water flow above permafrost in upland Arctic Alaska hillslopes

［J/OL］. Journal of Geophysical Research： Earth Surface， 2020， 125（2） ［2021-06-06］. .

[本文引用: 1]

Evans S G， Ge Shemin.

Contrasting hydrogeologic responses to warming in permafrost and seasonally frozen ground hillslopes

［J］. Geophysical Research Letters， 2016， 44： 1803-1813.

[本文引用: 2]

McKenzie J M， Voss C I.

Permafrost thaw in a nested groundwater-flow system

［J］. Hydrogeology Journal， 2013， 21（1）： 299-316.

[本文引用: 1]

Huang Kewei， Dai Junchen， Wang Genxu， et al.

The impact of land surface temperatures on suprapermafrost groundwater on the central Qinghai-Tibet Plateau

［J］. Hydrological Processes， 2020， 34（6）： 1475-1488.

[本文引用: 1]

Diersch H O G. FEFLOW： finite element modeling of flow， mass and heat transport in porous and fractured media［M］. Berlin： Springer， 2014.

[本文引用: 1]

Jan A， Coon E T， Painter S L.

Evaluating integrated surface/subsurface permafrost thermal hydrology models in ATS （v0.88） against observations from a polygonal tundra site

［J］. Geoscientific Model Development， 2020， 13（5）： 2259-2276.

[本文引用: 1]

Atchley A L， Painter S L， Harp D R， et al.

Using field observations to inform thermal hydrology models of permafrost dynamics with ATS （v0.83）

［J］. Geoscientific Model Development， 2015， 8（9）： 2701-2722.

[本文引用: 1]

Cochand F， Therrien R， Lemieux J M.

Integrated hydrological modeling of climate change impacts in a snow-influenced catchment

［J］. Groundwater， 2019， 57（1）： 3-20.

[本文引用: 1]

Therrien R. HydroGeoSphere： a three-dimensional numerical model describing fully-integrated sub-surface and surface flow and solute transport［M］. Waterloo， Ontario， Canada： University of Waterloo， 2017.

[本文引用: 2]

Kurylyk B L， Watanabe K.

The mathematical representation of freezing and thawing processes in variably-saturated， non-deformable soils

［J］. Advances in Water Resources， 2013， 60： 160-177.

[本文引用: 1]

The OpenFOAM Foundation.

OpenFOAM user guide

［EB/OL］. ［2021-06-06］. .

[本文引用: 3]

Mortensen M， Langtangen H P， Wells G N.

A FEniCS-based programming framework for modeling turbulent flow by the Reynolds-averaged Navier-Stokes equations

［J］. Advances in Water Resources， 2011， 34（9）： 1082-1101.

[本文引用: 1]

Miller C T， Dawson C N， Farthing M W， et al.

Numerical simulation of water resources problems： models， methods， and trends

［J］. Advances in Water Resources， 2013， 51： 405-437.

Yang Xiaofan， Mehmani Y， Perkins W A， et al.

Intercomparison of 3D pore-scale flow and solute transport simulation methods

［J］. Advances in Water Resources， 2016， 95： 176-189.

Orgogozo L， Prokushkin A S， Pokrovsky O S， et al.

Water and energy transfer modeling in a permafrost-dominated， forested catchment of Central Siberia： the key role of rooting depth

［J］. Permafrost and Periglacial Processes， 2019， 30（2）： 75-89.

[本文引用: 3]

Horgue P， Soulaine C， Franc J， et al.

An open-source toolbox for multiphase flow in porous media

［J］. Computer Physics Communications， 2015， 187： 217-226.

[本文引用: 1]

Bear J， Bachmat Y. Introduction to modeling of transport phenomena in porous media［M］. Berlin： Springer， 1990.

[本文引用: 1]

Sharp J J， Simmons C T.

The compleat Darcy： new lessons learned from the first English translation of Les Fontaines Publiques de la Ville de Dijon

［J］. Ground Water， 2005， 43（3）： 457-460.

[本文引用: 1]

McKenzie J M， Voss C I， Siegel D I.

Groundwater flow with energy transport and water-ice phase change： numerical simulations， benchmarks， and application to freezing in peat bogs

［J］. Advances in Water Resources， 2007， 30（4）： 966-983.

[本文引用: 1]

Jame Y W， Norum D I.

Heat and mass transfer in a freezing unsaturated porous medium

［J］. Water Resources Research， 1980， 16（4）： 811-819.

[本文引用: 1]

Hu Guojie， Zhao Lin， Zhu Xiaofan， et al.

Review of algorithms and parameterizations to determine unfrozen water content in frozen soil

［J/OL］. Geoderma， 2020， 368 ［2021-06-06］. .

[本文引用: 1]

Orgogozo L， Renon N， Soulaine C， et al.

An open source massively parallel solver for Richards equation： mechanistic modelling of water fluxes at the watershed scale

［J］. Computer Physics Communications， 2014， 185（12）： 3358-3371.

[本文引用: 1]

Weill S， Mouche E， Patin J.

A generalized Richards equation for surface/subsurface flow modelling

［J］. Journal of Hydrology， 2009， 366（1/2/3/4）： 9-20.

Šimůnek J， van Genuchten M T， Šejna M.

Development and applications of the HYDRUS and STANMOD software packages and related codes

［J］. Vadose Zone Journal， 2008， 7（2）： 587-600.

Williams G A， Miller C T.

An evaluation of temporally adaptive transformation approaches for solving Richards’ equation

［J］. Advances in Water Resources， 1999， 22（8）： 831-840.

[本文引用: 1]

Lunardini V J. Heat transfer with freezing and thawing［M］. New York： Elsevier， 1991.

[本文引用: 1]

M-k Woo， Arain M A， Mollinga M， et al.

A two-directional freeze and thaw algorithm for hydrologic and land surface modelling

［J/OL］. Geophysical Research Letters， 2004， 31（12） ［2021-06-06］. .

[本文引用: 1]

Rowland J C， Travis B J， Wilson C J.

The role of advective heat transport in talik development beneath lakes and ponds in discontinuous permafrost

［J/OL］. Geophysical Research Letters， 2011， 38（17） ［2021-06-06］. .

[本文引用: 1]

Bense V F， Kooi H， Ferguson G， et al.

Permafrost degradation as a control on hydrogeological regime shifts in a warming climate

［J/OL］. Journal of Geophysical Research： Earth Surface， 2012， 117（F3） ［2021-06-06］. .

Wellman T P， Voss C I， Walvoord M A.

Impacts of climate， lake size， and supra- and sub-permafrost groundwater flow on lake-talik evolution， Yukon Flats， Alaska （USA）

［J］. Hydrogeology Journal， 2013， 21（1）： 281-298.

Jafarov E E， Coon E T， Harp D R， et al.

Modeling the role of preferential snow accumulation in through talik development and hillslope groundwater flow in a transitional permafrost landscape

［J/OL］. Environmental Research Letters， 2018， 13（10） ［2021-06-06］. .

[本文引用: 1]

Svensson U， Ferry M.

DarcyTools： a computer code for hydrogeological analysis of nuclear waste repositories in fractured rock

［J］. Journal of Applied Mathematics and Physics， 2014， 2（6）： 365-383.

[本文引用: 1]

Mizoguchi M. Water， heat and salt transport in freezing soil

［D］. Tokyo， Japan： University of Tokyo， 1990.

[本文引用: 3]

Hansson K， Šimůnek J， Mizoguchi M， et al.

Water flow and heat transport in frozen soil： numerical solution and freeze-thaw applications

［J］. Vadose Zone Journal， 2004， 3（2）： 693-704.

[本文引用: 1]

Mochizuki H， Mizoguchi M， Miyazaki T.

Effects of NaCl concentration on the thermal conductivity of sand and glass beads with moisture contents at levels below field capacity

［J］. Soil Science and Plant Nutrition， 2008， 54（6）： 829-838.

[本文引用: 1]

/

〈 〉

陇ICP备2021001824号-10
　 甘公网安备 62010202000676号
版权所有 © 《冰川冻土》编辑部
地址：兰州市天水中路8号　
邮编：730000　电话：0931-8260767　
E-mail：edjgg@lzb.ac.cn
总访问量 　 今日访问 　 在线人数
技术支持：北京玛格泰克科技发展有限公司