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

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

（1.北京师范大学地理科学学部地表过程与资源生态国家重点实验室，北京100875；2.北京师范大学地理科学学部自然资源学院，北京100875）

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本身的工具包括前处理、后处理接口，以确保不同环境之间数据传输的协调性。Open-FOAM的整体结构如图1所示，包括核心的解算部分，以及前处理和后处理环境（如ParaView（www.paraview.org），Tecplot（www.tecplot.com）等）。

图1 OpenFOAM总体结构示意图Fig.1 Overall structure of the OpenFOAM

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

1.2 寒区土壤水热耦合模型

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

表1 模型参数Table 1 List of model parameters

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

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

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

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

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

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

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

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

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

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

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

目前，该模型已被成功编译至开源CFD软件OpenFOAM®（v1712版本）中，成为一类新的求解器，命名为darcyTHFOAM。模型的具体计算流程图如图2所示，在每个时间步长中，首先利用Open-FOAM中不完全的对角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 darcyTHFOAM计算流程图Fig.2 Flow chart of the darcyTHFOAM

表2 darcyTHFOAM的性能特点Table 2 Specific features of the darcyTHFOAM

总体来说，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所示。其中，用户输入正确的用户密码方可登陆，在“划分网格”的界面，定义算例的边界条件类型和计算区域；压力场、速度场和温度场的初始值、在特殊区域的值，分别在“设置初始值”和“设置特殊区域”中定义；土壤中各相的水热性质、土壤冻结特征曲线、相对导水率函数均在“设置水热参数”中定义；最后，用户可以在“求解”中自定义时间步长、模拟时间、输出时间间隔，所有的仿真模拟计算将在后台完成，而结果展示与分析将通过调用Para-View实现。

图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为冻结深度；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 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。在压力和温度的驱动下，初始处于冻结状态的冰块将慢慢融化。

表3 两个基准测试算例中的参数设置Table 3 Parameters used in two benchmark cases

图5 冰包裹体融化算例的计算区域及边界条件［49］（灰色字为流动边界条件，黑色字为传热边界条件）Fig.5 Melting of ice inclusion case：computational domain，initial and boundary conditions［49］（gray characters：flow，black characters：heat transfer）

在通过网格分辨率测试后，整个计算区域最终被划分为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 冰包裹体融化算例在不同时刻的温度演变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 冰包裹体融化算例在不同时刻的液态水饱和度演变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 冰包裹体融化算例沿中轴线的演变过程及对比验证Fig.8 Melting of ice inclusion case：simulated instantaneous liquid water saturation profiles along the central line of the computationaldomain（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）

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

表4 darcyTHFOAM与其他模型或软件数值算法和网格离散方法的比较Table 4 Mesh and numerical schemes used by darcyTHFOAM，Cast3M，permaFOAM，COMSOL and DarcyTools

2.3 融区贯穿/闭合算例

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

图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 融区贯穿/闭合算例在不同时刻的温度演变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 融区贯穿/闭合算例在不同时刻的液态水饱和度演变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 融区贯通/闭合算例沿中轴线的演变过程及对比验证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 室内冻结实验计算域设置（灰色字为流动边界条件，黑色字为传热边界条件）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

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