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

（1.北京師范大學地理科學學部地表過程與資源生態(tài)國家重點實驗室，北京100875；2.北京師范大學地理科學學部自然資源學院，北京100875）

0 引言

冰凍圈是世界上許多大江大河的發(fā)源地，是干旱區(qū)內(nèi)陸河流域的“水塔”［1-3］。凍土作為陸地冰凍圈的重要組成部分［4］，不僅對寒區(qū)水文過程有顯著影響［5-7］，而且對全球氣候變化極其敏感［8-10］。凍土本身具有相對隔水層的特點，使地下水文過程更加復(fù)雜［11-14］，進而影響地表水文過程、地-氣之間的能量交換和碳循環(huán)［15-18］。而凍融循環(huán)對寒區(qū)水文過程（特別是地下部分）的影響，包括土壤水分的變化［19］、徑流路徑的變化［20-21］、地下水的分布［22-23］以及陸地水儲量［24-26］等。此外，全球氣候變暖導(dǎo)致土壤溫度升高，多年凍土的退化加劇，季節(jié)凍結(jié)深度減小，融化深度增大；而地下存儲的碳以溫室氣體的形式被釋放到大氣中［27-30］，又對全球變暖產(chǎn)生正反饋［21，31］。因此，研究寒區(qū)土壤水文過程，即土壤水熱傳輸過程，是寒區(qū)水文學的關(guān)鍵問題之一。

由于寒區(qū)環(huán)境惡劣、人跡罕至，開展野外觀測較為困難，因此，野外水文氣象、水文地質(zhì)和地球物理觀測的范圍、頻率和數(shù)據(jù)量、數(shù)據(jù)精度等均受限［32-34］。隨著數(shù)值模擬和計算機科學的快速發(fā)展，基于近幾十年積累的觀測數(shù)據(jù)，構(gòu)建數(shù)值模型研究氣候和季節(jié)變化下的寒區(qū)土壤水熱傳輸過程，是理解寒區(qū)土壤水文物理過程，揭示其動力學機制的重要工具［35-37］。寒區(qū)土壤水熱傳輸?shù)奈锢磉^程較為復(fù)雜，涉及多相流動、冰水相態(tài)變化、凍脹等，使得描述水熱傳輸過程的本構(gòu)方程也具有高度的復(fù)雜性和非線性，這些都為建模和數(shù)值模擬帶來了巨大的科學挑戰(zhàn)［38］。目前，寒區(qū)水文模型主要包括基于陸地水文學的寒區(qū)分布式水文模型（distributed hydrologic model）［39］和基于水文地質(zhì)學的寒區(qū)水文地質(zhì)模型（groundwater/hydrogeologic model）［40］。寒區(qū)分布式水文模型，如CRHM［41］、GBEHM［42］、VIC［43］、HydroSiB2-SF［44］、WEB-DHM-pF［45］，是針對寒區(qū)水文循環(huán)中子過程/子單元分別建立子模型/模塊，并集成至統(tǒng)一的模型框架中。這類模型已盡可能集成了寒區(qū)水文循環(huán)中的所有環(huán)節(jié)，但存在空間分辨率較低（千米級）、參數(shù)眾多難以率定、需要大量觀測數(shù)據(jù)加以驗證的問題［46］，且側(cè)重于陸面地表過程，對于地下水文過程仍停留在參數(shù)化方案或一維地下水模型層面［47］。因此，此類模型無法精細刻畫寒區(qū)土壤水分運移和熱傳遞的物理過程，對凍融等寒區(qū)特殊現(xiàn)象的本質(zhì)成因和動力學機理研究不夠成熟［48］。源自地下水科學、水文地質(zhì)學和凍土水文學發(fā)展的寒區(qū)水文地質(zhì)模型（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］等，利用多孔介質(zhì)滲流和傳熱理論建立地下水多相流動和熱傳遞的本構(gòu)方程，并通過土壤水分凍融曲線或基于物理過程的本構(gòu)關(guān)系（例如Clausius-Clapyeron方程）來表征凍融循環(huán)對孔隙水相態(tài)、土壤孔隙率和滲透性的影響，進而模擬計算土壤水分和溫度、地下徑流量及其對流域地表徑流的補給量等［59］。但是，此類模型較為復(fù)雜，大多為封裝化或商業(yè)化軟件，較難進行二次開發(fā)，且計算精度和效率有待提高。

綜上所述，氣候變化影響下的寒區(qū)土壤水文過程較為復(fù)雜，包括多相流動、傳熱和冰水相態(tài)變化等。雖然寒區(qū)分布式水文模型近年來發(fā)展較快，但多將土壤水文過程進行了不同程度的簡化，且在模型建立、計算精度、并行計算效率等方面有待改進。此外，幾類主流的寒區(qū)水文地質(zhì)模型多為封裝化或商業(yè)化軟件，大大限制了模型的應(yīng)用和推廣。因此，亟需自主研發(fā)更穩(wěn)定高效、高精度、支持并行計算、便于二次開發(fā)的開源寒區(qū)土壤水熱耦合模型，精細刻畫土壤溫度、含水量與含冰量的時空動態(tài)變化，深入理解寒區(qū)土壤水文過程。此外，不同類別的寒區(qū)水文模型的研究對象和用途不盡相同，同一類別內(nèi)的模型所選用的水文物理過程表征方式、參數(shù)化方案、數(shù)值離散方法和網(wǎng)格類型等也千差萬別，其準確性和物理意義上的真實性均需校驗［58］。本研究旨在基于寒區(qū)土壤水文物理過程和計算流體力學方法（computational fluid dynamics，CFD），利用大型并行開源軟件OpenFOAM?（www.openfoam.com）自主研發(fā)一套高精度、高效率、適用于飽和狀態(tài)下的寒區(qū)土壤水熱耦合模型，用于描述飽和狀態(tài)下的寒區(qū)土壤水分與溫度的動態(tài)變化過程，并通過一維傳熱方程解析解、二維簡單基準測試算例和室內(nèi)凍結(jié)實驗對模型進行系統(tǒng)的驗證，綜合評估其物理意義上的準確性與完整性、計算精度及計算效率。

1 寒區(qū)土壤水熱耦合模型的構(gòu)建

1.1 計算流體力學軟件OpenFOAM簡介

OpenFOAM?（https：//openfoam.org）是一個完全由C++編寫的面向?qū)ο蟮挠嬎懔黧w力學（CFD）類庫［60］（約100個），用于創(chuàng)建可執(zhí)行文件（如應(yīng)用程序），其數(shù)值模擬理念是將偏微分方程進行有限體積離散化后獲得數(shù)值解。OpenFOAM自帶了約250個應(yīng)用程序，用戶也可根據(jù)自己需求自行編寫。應(yīng)用程序主要分為兩類：求解器與工具。其中，求解器用于解決特定的流體（或連續(xù)體）力學中的特定問題；工具用于執(zhí)行涉及數(shù)據(jù)操作等任務(wù)。此外，OpenFOAM本身的工具包括前處理、后處理接口，以確保不同環(huán)境之間數(shù)據(jù)傳輸?shù)膮f(xié)調(diào)性。Open-FOAM的整體結(jié)構(gòu)如圖1所示，包括核心的解算部分，以及前處理和后處理環(huán)境（如ParaView（www.paraview.org），Tecplot（www.tecplot.com）等）。

圖1 OpenFOAM總體結(jié)構(gòu)示意圖Fig.1 Overall structure of the OpenFOAM

作為目前最強大的計算流體力學類庫，Open-FOAM采用的數(shù)值離散方法為有限體積法（finite volume method，F(xiàn)VM），其自帶的網(wǎng)格生成工具snappyHexMesh可以快速高效的劃分六面體及多面體網(wǎng)格，因而可以處理復(fù)雜的幾何外形，且所生成的網(wǎng)格質(zhì)量較高；支持大型并行計算，計算效率高；研發(fā)環(huán)境良好，有利于對模型進行開發(fā)或二次開發(fā)；接口方式友好，便于耦合其他模型或求解器。因此，OpenFOAM已經(jīng)被越來越多地用于地球科學中與環(huán)境流體力學相關(guān)的模擬計算［61-64］。

1.2 寒區(qū)土壤水熱耦合模型

本研究針對寒區(qū)土壤水分運移和傳熱過程，建立了基于多相流動與傳熱理論的寒區(qū)土壤水熱耦合模型，其控制方程包括質(zhì)量守恒方程、動量方程和能量守恒定律。模型中包含的變量及參量定義匯總在表1中。

表1 模型參數(shù)Table 1 List of model parameters

飽和條件下（即不考慮氣相）的質(zhì)量守恒方程［65］為

考 慮 水 的 可 壓 縮 性［49］，引 入 壓 縮 系 數(shù)β=，其中P為壓力，則式（1）變?yōu)?/p>

若暫不考慮土壤的凍脹和融沉［66］，且在完全飽和情況下，θi=ε-θw，同時將地下水頭定義為h=，式（2）則變?yōu)?/p>

假設(shè)地下水流動滿足達西定律［67］，

式中：Kw為導(dǎo)水率，定義為

式中：K（rT）為與土壤凍結(jié)相關(guān)的相對導(dǎo)水率，可表示［51，68］為

式（6）描述了相對導(dǎo)水率Kr（T）與固態(tài)含冰量θi（T）之間的關(guān)系，原理是在凍結(jié)期由于土壤水成冰導(dǎo)致土壤孔隙變小，進而影響土壤水遷移［69］。此外，固態(tài)含冰量θi(T)=ε-θw(T)，而液態(tài)含水量θw是與溫度T相關(guān)的變量θw（T），通常用經(jīng)驗函數(shù)來表示兩者之間的關(guān)系，即土壤凍結(jié)特征曲線［70］。本模型中采用最為常用的指數(shù)形式的經(jīng)驗函數(shù)［49，51，64］來表示。

將式（4）代入式（3）中，并除以液態(tài)水的密度ρw，則最終的流動方程為

描述傳熱過程的能量守恒方程［49］為

式中：λt為土壤總導(dǎo)熱率，可表示為

由此可見，描述寒區(qū)地下水熱傳輸過程的流動和傳熱方程，由于多相流和凍融引起的水冰相態(tài)變化，具有高度的非線性。

目前，該模型已被成功編譯至開源CFD軟件OpenFOAM?（v1712版本）中，成為一類新的求解器，命名為darcyTHFOAM。模型的具體計算流程圖如圖2所示，在每個時間步長中，首先利用Open-FOAM中不完全的對角Cholesky預(yù)條件共軛求解器（DIC-PCG）［60］來求解流動方程（8），更新速度場（4）及土壤的水熱性質(zhì)，進而利用對角不完全LU穩(wěn)定化預(yù)條件雙共軛求解器（DILC-PBiCGStab）［60］來求解傳熱方程（9），用Picard循環(huán)來處理方程中的非線性問題。其中，Δt為時間步長，NPicP為求解壓力的Picard迭代次數(shù)，NIterP為求解壓力的最大Picard迭代次數(shù)，P為壓力，ErrP為求解壓力的誤差，PicP為求解壓力的最大誤差，U為速度，T為溫度，NPicT為求解溫度的Picard迭代次數(shù)，NIterT為求解溫度的最大Picard迭代次數(shù)，ErrT為求解溫度的誤差，PicT為求解溫度的最大誤差，tfactor為自動調(diào)整時間步長的參數(shù)。每個時間步長Δt中，在求解壓力的流動方程（8）時，如果誤差ErrP超過最大誤差PicP，則在Picard最大迭代次數(shù)NIterP內(nèi)繼續(xù)求解，直至誤差ErrP低于最大誤差PicP，進而求解速度；同樣的，在求解傳熱方程（9）時，如果誤差ErrT超過最大誤差PicT，則在Picard最大迭代次數(shù)NPicT內(nèi)繼續(xù)求解，直至誤差ErrT低于最大誤差PicT；如果求解溫度與壓力的誤差均滿足最大誤差的限制，則繼續(xù)下一時間步長的計算，否則，自動調(diào)整時間步長，重新展開計算，直至滿足誤差需求。另外，darcyTHFOAM的主要性能特點如表2所示，適用于Linux系統(tǒng)，編程語言為C++，時間離散方法為歐拉格式，空間離散方法為線性插值格式，采用自適應(yīng)的時間步長策略來平衡計算精度和計算效率［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模型的基本假設(shè)為：①不考慮氣相，含水層處于完全飽和狀態(tài)下；②不考慮土壤的凍脹和融沉；③土壤水流動滿足達西定律；④采用指數(shù)形式的土壤凍結(jié)特征曲線，來表示液態(tài)含水量與溫度的關(guān)系。而darcyTHFOAM模型是在開源的計算流體力學軟件OpenFOAM的基礎(chǔ)上開發(fā)的，故支持高效的并行計算（適用于mm至km級的模擬計算）；時間步長采用自適應(yīng)的時間步長策略（可從μs到h）；OpenFOAM采用六面體來劃分網(wǎng)格，空間分辨率可精細至mm級；網(wǎng)格本身是三維，但通過調(diào)整各個維度的長度，進而可使模型適用于一維、二維、三維的完全飽和算例。該模型適用于含水層完全飽和狀態(tài)下的凍結(jié)與融化過程，如飽和狀態(tài)下的基準測試算例研究（冰包裹體融化算例、融區(qū)貫通/閉合算例）、飽和狀態(tài)下的室內(nèi)土柱凍結(jié)實驗等等，但是，darcyTHFOAM模型的應(yīng)用范圍不僅限于此，其適用于一維到三維的問題，時空分辨率較廣（從mm到km，從μs到h），邊界條件靈活可控，并且，除了土壤含水層外，還可用于其他多孔介質(zhì)。此外，darcyTHFOAM模型接口非常靈活，用戶可根據(jù)需求耦合新的方程（例如溶質(zhì)運移）、引入復(fù)雜的邊界條件與新的源項等。最后，darcyTHFOAM模型也可根據(jù)需求在現(xiàn)有的本構(gòu)方程基礎(chǔ)上進行修改，使其適用于非飽和條件。

1.3 軟件

為方便用戶使用，利用Python 3.5（www.python.org）和Qt設(shè)計器（www.qt.io）設(shè)計了對應(yīng)的可視化界面軟件。該軟件支持Linux系統(tǒng)，包括8個主界面：登錄、打開算例文件、劃分網(wǎng)格、設(shè)置初始值、設(shè)置水熱參數(shù)、設(shè)置特殊區(qū)域、求解、可視化，如圖3所示。其中，用戶輸入正確的用戶密碼方可登陸，在“劃分網(wǎng)格”的界面，定義算例的邊界條件類型和計算區(qū)域；壓力場、速度場和溫度場的初始值、在特殊區(qū)域的值，分別在“設(shè)置初始值”和“設(shè)置特殊區(qū)域”中定義；土壤中各相的水熱性質(zhì)、土壤凍結(jié)特征曲線、相對導(dǎo)水率函數(shù)均在“設(shè)置水熱參數(shù)”中定義；最后，用戶可以在“求解”中自定義時間步長、模擬時間、輸出時間間隔，所有的仿真模擬計算將在后臺完成，而結(jié)果展示與分析將通過調(diào)用Para-View實現(xiàn)。

圖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 寒區(qū)土壤水熱耦合模型的驗證

在模型研究中，選擇合適穩(wěn)定且具有物理意義的基準測試算例來驗證，是非常必要的。2014年底，為全面驗證和優(yōu)化寒區(qū)地下水文模型，由美國、加拿大、瑞典、德國、英國、荷蘭和法國等國的科學家組成了科學小組，并啟動了“寒區(qū)地下水熱耦合模型比較計劃（InterFrost）”（wiki.lsce.ipsl.fr/interfrost）。InterFrost計劃的主旨是提出一系列從簡單到復(fù)雜，且具有代表性和實用性的基準測試算例，用于模型比較和基準驗證，進而優(yōu)化和發(fā)展模型。因此，在本研究中采用一維解析解、InterFrost的兩個基準測試算例（冰包裹體融化和融區(qū)貫穿/閉合）以及室內(nèi)凍結(jié)實驗對darcyTHFOAM模型逐步進行驗證。

2.1 一維Stefan傳熱方程

一維Stefan方程存在解析解［34，75-76］，可根據(jù)溫度來估算凍結(jié)深度隨時間的變化。

式中：Zf為凍結(jié)深度；t為時間；λu為熱導(dǎo)率。

為更好地與Stefan方程的解析解對比，Schilling等［34］采用HGS模擬了一維垂直土柱的融化過程。根據(jù)他們研究中的算例設(shè)置，假定土柱是均質(zhì)的，高10 m，起初處于-5℃的完全凍結(jié)狀態(tài)，土柱的頂部和底部為恒定為10℃和-5℃。土壤固體顆粒的導(dǎo)熱率λs=0.017 J·s-1·cm-3·K-1，比熱容Cs=35.07 J·cm-3·K-1，孔隙度ε=0.4。為將模擬結(jié)果與解析解對比，將土壤凍結(jié)特征曲線中的擬合參數(shù)W設(shè)置為4，其他參數(shù)與以下基準測試算例相同。

整個計算區(qū)域在垂直方向上劃分了400個網(wǎng)格，起始步長和最大步長分別為0.1 s和30 min，利用自主開發(fā)的darcyTHFOAM模擬計算了90天的土柱融化過程，在本機（Window系統(tǒng)/4G RAM/1T內(nèi)存）下的虛擬機（Ubuntu18.04系統(tǒng)/2G RAM/40G內(nèi)存）（下同）上運行時間約為48 s。土柱中凍結(jié)鋒面與土柱頂部的距離變化，最能反映土柱的融化過程，且能夠與解析解比較。如圖4所示，darcyTHFOAM模擬的凍結(jié)深度與Stefan方程的解析解和HGS的模擬結(jié)果相比，結(jié)果均非常準確。

圖4 darcyTHFOAM模擬的凍結(jié)深度與Stefan方程的解析解和HGS模擬結(jié)果的相互比較Fig.4 Inter-comparison of the frost depth calculated using Stefan’s equation and simulated by the HGS and darcyTHFOAM

2.2 冰包裹體融化算例

如圖5所示，整個計算區(qū)域長為3 m，寬為1 m，在內(nèi)部有一邊長為0.333 m的正方形冰塊。左側(cè)入口處設(shè)置為恒定的5℃，其他邊界均為絕熱邊界，除了藍色冰塊為-5℃，其他區(qū)域溫度均為5℃。而針對壓力水頭邊界條件，在左側(cè)入口與右側(cè)出口邊界分別設(shè)置為固定值0.09 m、0 m，上下邊界均為無通量邊界，內(nèi)部初始值與右側(cè)出口邊界的值相同，其他相關(guān)參數(shù)見表3。在壓力和溫度的驅(qū)動下，初始處于凍結(jié)狀態(tài)的冰塊將慢慢融化。

表3 兩個基準測試算例中的參數(shù)設(shè)置Table 3 Parameters used in two benchmark cases

圖5 冰包裹體融化算例的計算區(qū)域及邊界條件［49］（灰色字為流動邊界條件，黑色字為傳熱邊界條件）Fig.5 Melting of ice inclusion case：computational domain，initial and boundary conditions［49］（gray characters：flow，black characters：heat transfer）

在通過網(wǎng)格分辨率測試后，整個計算區(qū)域最終被劃分為300×100個網(wǎng)格，起始時間步長為0.1 s，隨后在計算過程中可自動調(diào)整。每個時間步長，壓力和溫度的收斂標準均設(shè)置為1×10-10。在求解壓力的Picard循環(huán)中，最大迭代次數(shù)和最大誤差分別設(shè)置為20和1×10-3，而在求解溫度的Picard循環(huán)中，最大迭代次數(shù)和最大誤差分別設(shè)置為50和1×10-5。為提高計算效率，將計算區(qū)域劃分為120個部分，在天河二號超級計算器上調(diào)用了5個核（共120個節(jié)點）進行仿真模擬計算，最終用時3 min完成了5天的冰塊融化過程。

圖6 和圖7分別展示了溫度和液態(tài)水飽和度（液態(tài)水含量/孔隙度）隨時間的變化，進而反映了冰塊在不同時刻的融化情況。在冰塊融化過程中，沿計算區(qū)域中心線的液態(tài)水飽和度如圖8（a）所示，從上到下是從初始到最終時刻的液態(tài)水飽和度，最后全部融化，與頂部坐標軸完全重合；沿計算區(qū)域中心線的溫度如圖8（b）所示，從下到上表示初始和最終時刻的溫度，黑色箭頭表示溫度剖面的變化。在壓力水頭梯度的驅(qū)動下，融化加速，較冷的流體通過平流和熱擴散向下游輸送，最終達到與初始條件相同的溫度（5℃）。為了更好地將溫度的時空變化與Interfrost中給出的其他模型的結(jié)果進行量化對比，參照Grenier等［49］中的評價指標，選取了計算區(qū)域的最低溫度進行對比。此外，其他模型或軟件的數(shù)值方案及網(wǎng)格劃分情況已經(jīng)總結(jié)在表4中。結(jié)果表明，darcyTHFOAM模擬的最低溫度與其他模型預(yù)測的最低溫度一致［圖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 冰包裹體融化算例在不同時刻的液態(tài)水飽和度演變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℃的階段呈現(xiàn)出一定差異，特別是在-1～0℃（對應(yīng)相變階段），溫度呈現(xiàn)出急速上升趨勢的開始時刻，darcyTHFOAM，分別比COMSOL和darcyTools晚2 400 s和1 300 s，而又比Cast3M和permaFOAM提早3 500 s和5 400 s，這可能是由于不同的數(shù)值算法和網(wǎng)格離散方法造成的。與同類模型相比，darcyTHFOAM采用更少的網(wǎng)格數(shù)量便可達到相同的模擬結(jié)果（如表4所示），且支持并行計算，計算效率更高，接口較為靈活，便于二次開發(fā)。但目前僅適用于完全飽和狀態(tài)下的水熱耦合過程，未來可拓展至模擬變飽和狀態(tài)下的土壤水熱傳輸過程。

表4 darcyTHFOAM與其他模型或軟件數(shù)值算法和網(wǎng)格離散方法的比較Table 4 Mesh and numerical schemes used by darcyTHFOAM，Cast3M，permaFOAM，COMSOL and DarcyTools

2.3 融區(qū)貫穿/閉合算例

貫穿融區(qū)（“天窗”）存在于湖泊或河流下面的多年凍土中，是多年凍土內(nèi)的局部融化區(qū)［77-80］。該基準測試算例取材于自然界中凍土的典型特征，兼具代表性與適用性［49］，能夠反映出多年凍土中貫穿融區(qū)的演變過程。兩個起始凍結(jié)的半圓形區(qū)域（半徑為0.5099 m）溫度都是-5℃，代表了兩部分凍土，均位于一個正方形計算域內(nèi)（1 m×1 m）。最初，計算區(qū)域的背景溫度為5℃，左側(cè)入口邊界賦予固定值5℃，上下邊界賦予固定值-5℃（圖9），右側(cè)出口邊界為絕熱邊界，壓力邊界條件與上一個冰包裹體融化算例的相同，其他相關(guān)參數(shù)見表3。

圖9 融區(qū)貫穿/閉合算例的計算區(qū)域及邊界條件［49］（灰色字為流動邊界條件，黑色字為傳熱邊界條件）Fig.9 Talik opening/closure case：computational domain，initial and boundary conditions［49］（gray characters：flow，black characters：heat transfer）

經(jīng)過網(wǎng)格分辨率測試，整個計算區(qū)域最終被劃分為100×100個網(wǎng)格，起始時間步長為0.1 s，在計算的過程可自動調(diào)整，最大時間步長設(shè)置為1 000 s。對于每個時間步長，壓力和溫度的收斂標準均設(shè)置為1×10-10。在求解壓力的Picard循環(huán)中，最大迭代次數(shù)和最大誤差分別設(shè)置為20和1×10-3，而在求解溫度的Picard循環(huán)中，最大迭代次數(shù)和最大誤差分別設(shè)置為50和1×10-5。為提高計算效率，將計算區(qū)域分為了4部分，在天河二號超級計算器上采用1個核（共4個節(jié)點）進行仿真模擬計算，最終耗費2 min完成了40 h的融區(qū)演變過程模擬。

不同時刻的溫度和液態(tài)水飽和度分布情況如圖10和圖11所示，計算區(qū)域中心垂直線上的液態(tài)水飽和度和溫度變化如圖12（a）和圖12（b）所示，從上到下依次代表，從初始到終止時刻，其中，黑色箭頭表示溫度剖面的旋轉(zhuǎn)。在溫度與壓力的驅(qū)動下，兩個最初凍結(jié)的半圓區(qū)域逐漸融化。為了將darcyTHFOAM的模擬結(jié)果與Grenier等［49］中其他求解器的模擬數(shù)據(jù)進行定量比較，在圖12（c）中繪制了兩個選定位置（Pt1點和Pt2點）的溫度曲線。這兩個點都接近凍結(jié)區(qū)和非凍結(jié)區(qū)之間的初始邊界，對融區(qū)的貫通/閉合過程非常敏感。結(jié)果表明，這兩點處的溫度變化與其他模型模擬的結(jié)果非常一致，再次證明了模型的可靠性。

圖10 融區(qū)貫穿/閉合算例在不同時刻的溫度演變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 融區(qū)貫穿/閉合算例在不同時刻的液態(tài)水飽和度演變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 融區(qū)貫通/閉合算例沿中軸線的演變過程及對比驗證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 室內(nèi)土柱凍結(jié)實驗

Mizoguchi［82］早在1991年就開展了室內(nèi)土柱凍結(jié)實驗，并在實驗過程中進行了溫度觀測。隨著數(shù)值模型的發(fā)展，已有越來越多的寒區(qū)水文地質(zhì)模型和室內(nèi)實驗進行比對［83］，因此將darcyTHFOAM模擬此室內(nèi)凍結(jié)實驗的結(jié)果與實測的溫度數(shù)據(jù)比對，是模型驗證的有效途徑。

室內(nèi)土柱凍結(jié)實驗是在石英砂填充的飽和環(huán)境中進行的，為方便比較，在數(shù)值模擬時采用軸對稱區(qū)域，初始和邊界條件如圖13所示，整個土柱的初始溫度為5℃，兩側(cè)均為絕熱邊界，而頂部和底部分別暴露于溫度為-10℃和5℃的循環(huán)流體。針對壓力邊界條件，在土柱的底部、頂部和兩側(cè)均采用零通量邊界條件。因此，在溫度梯度的驅(qū)動下，土柱將從上到下開始凍結(jié)。

圖13 室內(nèi)凍結(jié)實驗計算域設(shè)置（灰色字為流動邊界條件，黑色字為傳熱邊界條件）Fig.13 Laboratory freezing experiment：computational domain，initial and boundary conditions（gray characters：flow，black characters：heat transfer）

表5 中總結(jié)了實驗參數(shù)，其中，土壤固體顆粒Cs的導(dǎo)熱系數(shù)是從Mochizuki等［84］中獲得，相對導(dǎo)水率函數(shù)中的阻抗系數(shù)Ω，參考前面基準測試算例設(shè)置為50。通過對比模擬結(jié)果和實驗數(shù)據(jù)，將土壤凍結(jié)函數(shù)中的擬合參數(shù)W設(shè)為5.25。考慮到計算效率和精度，經(jīng)過網(wǎng)格分辨率測試，將計算域離散為80×500個網(wǎng)格。該算例的收斂準則和Picard循環(huán)的設(shè)置，也與基準測試算例相同。初始時間步長設(shè)置為0.1 s，最大時間步長設(shè)置為30 min。采用darcyTHFOAM在本機上模擬計算50 h的凍結(jié)過程，耗時約3 min（CPU時間）。

表5 室內(nèi)凍結(jié)實驗中的參數(shù)Table 5 Parameters used in simulating laboratory freezing experiment

利用darcyTHFOAM展開仿真模擬計算后，1、3、6、12、24 h的溫度和液態(tài)水飽和度分布情況分別如圖14（a）和圖14（b）所示。從頂部到底部，溫度逐漸下降，同時液態(tài)水飽和度逐漸減少，直觀反映出了土柱的凍結(jié)過程。為了進一步驗證模型，從文獻［82］的圖3.1-5中提取出實驗觀測的溫度數(shù)據(jù)（Engauge Digitize，http：//digitizer.sourceforge.net）。圖15中展示了在不同時刻沿對稱軸的溫度分布，darcyTHFOAM模擬的數(shù)值結(jié)果與實驗數(shù)據(jù)非常吻合，進一步證明了模型的準確性。Mizoguchi［82］在設(shè)計該室內(nèi)凍結(jié)實驗時，整個土柱的初始溫度為5℃，而頂部和底部分別暴露于溫度為-10℃和5℃的循環(huán)流體。但從初始0 h的觀測數(shù)據(jù)來看，整個土柱內(nèi)部溫度平均在7℃，只有頂部溫度為5℃，說明在土柱凍結(jié)實驗開始時內(nèi)部并沒有完全達到5℃，進而導(dǎo)致觀測數(shù)據(jù)在早期（1 h）呈現(xiàn)出高估的趨勢（高達6℃）。

圖14 在凍結(jié)開始后1、3、6、12、24小時的溫度（a）和液態(tài)水飽和度（b）的演變Fig.14 Evolution of temperature（a）and liquid water saturation（b）at 1，3，6，12 and 24 h after freezing started

圖15 室內(nèi)凍結(jié)實驗數(shù)據(jù)與數(shù)值模擬結(jié)果對比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 結(jié)論與展望

本研究基于多孔介質(zhì)滲流和傳熱理論、寒區(qū)土壤水文物理過程和計算流體力學方法，建立了高精度、高效率的寒區(qū)土壤水熱耦合模型（darcyTHFOAM），綜合使用Stefan方程解析解、基準測試算例（冰包裹體融化、融區(qū)貫穿/閉合）和室內(nèi)實驗對模型進行了系統(tǒng)的校驗。該模型能夠精細刻畫寒區(qū)土壤中水分遷移和熱量傳輸過程，尤其是冰水相態(tài)的動態(tài)變化，為探究寒區(qū)土壤水熱傳輸?shù)奈锢磉^程、動力學機理和影響機制提供強大的數(shù)值模擬工具，有望為氣候變化引起的凍土退化、季節(jié)凍土變化和水資源短缺等地球系統(tǒng)科學問題提供理論依據(jù)和可靠預(yù)測。此外，模型的研發(fā)環(huán)境較為友好，具有一定的魯棒性，方便二次開發(fā)，未來將繼續(xù)耦合地表過程（如積雪消融）和生物地球化學反應(yīng)動力學過程，以研究寒區(qū)地下水-地表水相互作用和生物地球化學循環(huán)。