基于計算流體力學的寒區土壤水熱耦合模型研究

胡錦華， 陸 崢， 仝金輝， 李小雁， 劉紹民， 楊曉帆

（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方程解析解、基準測試算例（冰包裹體融化、融區貫穿/閉合）和室內實驗對模型進行了系統的校驗。該模型能夠精細刻畫寒區土壤中水分遷移和熱量傳輸過程，尤其是冰水相態的動態變化，為探究寒區土壤水熱傳輸的物理過程、動力學機理和影響機制提供強大的數值模擬工具，有望為氣候變化引起的凍土退化、季節凍土變化和水資源短缺等地球系統科學問題提供理論依據和可靠預測。此外，模型的研發環境較為友好，具有一定的魯棒性，方便二次開發，未來將繼續耦合地表過程（如積雪消融）和生物地球化學反應動力學過程，以研究寒區地下水-地表水相互作用和生物地球化學循環。