考慮土體非均質(zhì)傳熱的河岸帶水熱耦合模型

張文兵，沈振中，任 杰，徐力群，周聰聰，楊金孟

(1. 河海大學(xué)水文水資源與水利工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210098；2. 上海海事大學(xué)海洋科學(xué)與工程學(xué)院，上海 201306；3. 河海大學(xué)水利水電學(xué)院，江蘇 南京 210098；4. 西安理工大學(xué)省部共建西北旱區(qū)生態(tài)水利國家重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710048)

河岸帶是指介于水生生態(tài)系統(tǒng)與陸地生態(tài)系統(tǒng)的過渡帶[1-2]，是河流系統(tǒng)的重要屏障，能有效調(diào)蓄洪水、削減污染和保護(hù)水土環(huán)境，對維護(hù)河流生態(tài)系統(tǒng)健康具有重要作用[3-6]。河流地表水通過河岸帶沉積層與地下水發(fā)生水熱交換的區(qū)域被稱為河岸帶潛流層[7]，它是河岸帶邊緣效應(yīng)和功能的重要體現(xiàn)區(qū)[8]。河岸帶潛流層地下水與地表水交換過程時(shí)刻伴隨著能量的傳遞和交換，溫度作為能量的直觀載體，是反映地表水與地下水交換過程時(shí)空變化的主要表征因子，具有無污染、易觀測和受擾動小等特點(diǎn)，是理想的天然示蹤劑[9]。近年來，隨著溫度自動化觀測設(shè)備的革新及水熱耦合機(jī)理研究的深入，對河岸帶地表水與地下水交換速率及過程模式的分析已從傳統(tǒng)的水文學(xué)及水動力學(xué)方法逐漸發(fā)展到溫度示蹤法[10]?，F(xiàn)有的河岸帶水熱交換過程研究多側(cè)重于試驗(yàn)手段[11-15]，對河岸帶水熱耦合模型的研究鮮見報(bào)道。因此，有必要構(gòu)建合適的河岸帶水熱耦合模型以評估其內(nèi)部水熱動態(tài)變化過程，有助于揭示河岸帶內(nèi)水分運(yùn)移和熱量交換規(guī)律，對反演河岸帶污染物的遷移過程以及揭示河岸帶的環(huán)境效應(yīng)具有重要意義。

土體熱性質(zhì)參數(shù)作為水熱耦合模型的重要驅(qū)動因子，是決定模型是否成功和準(zhǔn)確的關(guān)鍵[16]，其中，有效導(dǎo)熱系數(shù)是影響和決定土體在傳熱過程中溫度分布的重要參數(shù)[17]。在現(xiàn)有的水熱耦合模擬中，土體有效導(dǎo)熱系數(shù)均被視作固定值，例如：Su等[18]基于GA-VS2DH開發(fā)的水熱耦合模擬方法，研究了大汶河和秦淮河2個(gè)不同介質(zhì)類型的河岸帶水熱交換過程；Ren等[19]基于COMSOL構(gòu)建了二維河岸帶飽和-非飽和水熱耦合模型，分析了河岸帶溫度場在不同季節(jié)的變化規(guī)律；Ren等[20]基于洞庭湖河岸帶某典型斷面的實(shí)測溫度和水位資料，驗(yàn)證了所提水熱耦合模型的合理性，這些模型均未考慮水熱耦合模擬中土體的非均質(zhì)傳熱問題。已有研究表明[16-17]，土體有效導(dǎo)熱系數(shù)與土體類型、質(zhì)地、粒徑分布、孔隙度以及飽和度等因素均有關(guān)，并且孔隙度和飽和度的影響最大，這對涉及非飽和問題的河岸帶水熱交換研究至關(guān)重要。因此，對于河岸帶水熱耦合模擬研究，還需考慮因土體各部分含水率不同而導(dǎo)致的非均質(zhì)傳熱問題。土體有效導(dǎo)熱系數(shù)模型是通過建立導(dǎo)熱系數(shù)與含水率之間的關(guān)系，進(jìn)而實(shí)現(xiàn)對不同飽和狀態(tài)土體有效導(dǎo)熱系數(shù)的預(yù)測，若將其引入河岸帶水熱耦合模型中，則能實(shí)現(xiàn)對河岸帶水熱交換模擬過程中土體非均質(zhì)傳熱問題的考慮。

本文在飽和-非飽和滲流及多孔介質(zhì)傳熱理論基礎(chǔ)上，引入土體有效導(dǎo)熱系數(shù)模型，構(gòu)建考慮土體非均質(zhì)傳熱的河岸帶水熱耦合模型，給出其在COMSOL商業(yè)有限元軟件中的實(shí)現(xiàn)方法，并通過野外原型觀測資料驗(yàn)證和對比分析不同有效導(dǎo)熱系數(shù)模型下河岸帶水熱耦合模型的模擬效果，以期為河岸帶水熱耦合模型的構(gòu)建及有效導(dǎo)熱系數(shù)模型的選取提供借鑒。

1 河岸帶水熱耦合數(shù)學(xué)模型

1.1 流動方程

河岸帶飽和-非飽和瞬態(tài)滲流場方程采用Richards方程描述：

(1)

式中：ρw為水的密度，kg/m3；Cm為容水度，m-1；g為重力加速度，m/s2；Se為土體的相對飽和度；Ss為彈性貯水率，Pa-1；p為壓力，Pa；t為時(shí)間，s；?為拉普拉斯算子；θ為體積含水率，m3/m3；ks為飽和滲透系數(shù)，m/s；kr(θ)為非飽和帶相對滲透系數(shù)，m/s，是體積含水率的函數(shù)；μ(T)為水的動力黏度，Pa·s，μ(T)=0.000 024 24×10247.8/(T+133.16)，T為溫度， ℃；z為計(jì)算點(diǎn)高程，m；Qm為滲流場源匯項(xiàng)，kg/(m3·s)。

土壤水分特征采用van Genuchten模型描述：

θ=θr+Se(θs-θr)

(2)

(3)

(4)

(5)

式中：θr為殘余體積含水率，m3/m3；θs為飽和體積含水率，m3/m3；hc為壓力水頭，m；α為水分特征曲線進(jìn)氣值的倒數(shù)，m-1；β為水分特征曲線坡度的指示參數(shù)，通過擬合土壤水分特征曲線得到，m=1-1/β。

1.2 熱量交換方程

河岸帶中水體與土體之間的熱量交換可采用下式表示：

(6)

式中：Cw為水的體積熱容，J/(m3·℃)；n為孔隙度；Cs為土體的體積熱容，J/(m3·℃)；Keff為土體有效導(dǎo)熱系數(shù)，W/(m·℃)；DHi,j為水動力彌散系數(shù)，m2/s；u為平均流速，m/s，u=v/θ，v為Darcy滲流速度；Qs為溫度場源匯項(xiàng)，W/m3。

式(6)左邊表示溫度在變飽和條件下隨時(shí)間的變化，即為非穩(wěn)態(tài)項(xiàng)；等式右邊第1項(xiàng)表示熱傳導(dǎo)項(xiàng)，第2項(xiàng)表示熱彌散項(xiàng)，第3項(xiàng)表示熱對流項(xiàng)。

近兩年外部環(huán)境不好，但對于黑龍江銷售來說，最大的影響來自當(dāng)?shù)胤欠ǖ褂蛯医恢埂！皷|北地區(qū)高峰時(shí)一天能進(jìn)來50臺油槽車、合計(jì)1500噸的量，算下來對黑龍江銷售整體銷量影響能達(dá)到10%?！焙樗蓾f，雖然今年年初各個(gè)片區(qū)聯(lián)合當(dāng)?shù)貓?zhí)法部門對其實(shí)施打壓，但風(fēng)聲一過又出來了。

其中，水動力彌散系數(shù)可表示為

(7)

式中：αT為橫向彌散度，m；|v|為流速矢量的大小，m/s；δij為克里格常量，當(dāng)i=j時(shí)為1，否則為0；αL為縱向彌散度，m；vi、vj分別為速度矢量的第i個(gè)和第j個(gè)分量，m/s。

1.3 有效導(dǎo)熱系數(shù)模型

目前，有關(guān)預(yù)測土體有效導(dǎo)熱系數(shù)的模型眾多，這些模型可分為經(jīng)驗(yàn)?zāi)Ｐ秃屠碚撃Ｐ?。相較于經(jīng)驗(yàn)?zāi)Ｐ?，理論模型往往預(yù)測精度不高、公式復(fù)雜且部分參數(shù)不易確定，難以直接應(yīng)用[16-17]。在經(jīng)驗(yàn)?zāi)Ｐ头矫妫K李君等[21]總結(jié)了部分具有代表性的土體有效導(dǎo)熱系數(shù)模型，并提出了優(yōu)化改進(jìn)模型。表1總結(jié)和歸納了文獻(xiàn)[21]中所涉及的土體有效導(dǎo)熱系數(shù)模型，以期應(yīng)用于河岸帶水熱耦合模型中，并比較在河岸帶水熱交換模擬中的表現(xiàn)。

表1 土體有效導(dǎo)熱系數(shù)經(jīng)驗(yàn)?zāi)Ｐ涂偨Y(jié)Table 1 Summary of thermal conductivity empirical model of soils

2 河岸帶水熱耦合模型實(shí)現(xiàn)方法

針對構(gòu)建的河岸帶水熱耦合模型，采用適合多場耦合計(jì)算的COMSOL軟件實(shí)現(xiàn)模型的開發(fā)及應(yīng)用[26-27]。對于河岸帶飽和-非飽和滲流問題，選用COMSOL內(nèi)置地下水流模塊中的Richards方程求解計(jì)算，而對于河岸帶傳熱問題，COMSOL內(nèi)置的多孔介質(zhì)傳熱模塊未能考慮土體有效導(dǎo)熱系數(shù)模型，因此需借助相關(guān)二次開發(fā)接口來實(shí)現(xiàn)對土體非均質(zhì)傳熱的考慮。這里采用COMSOL提供的自定義偏微分方程(PDE)模塊用于等效代替多孔介質(zhì)傳熱模塊求解河岸帶溫度變化，COMSOL提供的系數(shù)型PDE如下：

(8)

式中：M為因變量；ea、da、c、α、?、β、a和f均為自定義系數(shù)。

上述自定義系數(shù)可根據(jù)用戶的需求定義為常數(shù)或者不同類型的函數(shù)，并且函數(shù)可以是不同階次的導(dǎo)數(shù)，既可以時(shí)間相關(guān)，也可以空間相關(guān)，其自變量可以是獨(dú)立變量，也可以是求解變量本身。為利用系數(shù)型PDE等效熱量交換方程，需將式(6)按照式(8)的形式進(jìn)行整理，然后在相應(yīng)的編輯框中定義各系數(shù)(即ea=0,M=T,da=θCw+(1-n)Cs,c=Keff+θCwDHi,j,α=0,γ=0,β=θCwu,a=0,f=0)，則修改后的式(8)可等效熱量交換方程。需要注意的是，Keff主要通過θ(θ=nSr)建立與孔隙度和飽和度(Sr)的關(guān)系，反映非均質(zhì)參數(shù)的空間分布。因此，需增加預(yù)置儲存模塊用于計(jì)算和儲存每個(gè)時(shí)間步內(nèi)不同空間位置的有效導(dǎo)熱系數(shù)，該過程可以通過在模型樹的“組件”選項(xiàng)下定義局部變量實(shí)現(xiàn)。

河岸帶水熱耦合模型的有限元求解方法按照滲流場和溫度場的直接耦合和間接耦合，可分為直接求解算法和分步求解算法。直接求解算法是通過Newton迭代直接求解式(1)和式(8)，即在每個(gè)時(shí)間步內(nèi)同時(shí)求解單元節(jié)點(diǎn)的壓力和溫度，但由于水力特性對飽和度和壓力的依賴性，導(dǎo)致滲流場方程高度非線性，并且有效導(dǎo)熱系數(shù)模型的引入使得這種非線性程度增強(qiáng)，因此在利用直接求解算法時(shí)，結(jié)果往往不易收斂。為避免模型計(jì)算結(jié)果收斂性差的問題，可采用分步式求解算法對壓力(p)和溫度(T)進(jìn)行解耦計(jì)算，即將滲流場與溫度場劃分為2個(gè)獨(dú)立的場進(jìn)行求解，并通過變量交換實(shí)現(xiàn)耦合作用。分步求解過程具體可描述為：首先，初始化變量，以獲得初始壓力(p0)、溫度(T0)和有效導(dǎo)熱系數(shù)(Keff0)，然后在tn+1

3 模型驗(yàn)證及對比分析

3.1 試驗(yàn)數(shù)據(jù)及模型設(shè)置

為研究美國沃克湖流域相關(guān)河流的滲漏損失，美國墾務(wù)局以及弗吉尼亞州雷斯頓地質(zhì)調(diào)查局的研究人員在沃克湖流域的部分河床及河岸帶埋設(shè)了溫度和水位監(jiān)測裝置，對該區(qū)域地表水與地下水溫度和水位進(jìn)行長期動態(tài)監(jiān)測[28](圖1)。首先將帶有濾網(wǎng)包裹的PVC管打入河岸帶沉積層，然后將3～4個(gè)獨(dú)立的防水溫度傳感器分別懸掛在距離河道中心約2.8 m和5.8 m處的PVC管中，分別測量距離地表0.95 m、1.40 m、2.00 m和1.50 m、1.95 m、2.30 m處沉積層的溫度變化。此外，在河道中還布設(shè)有水溫、水位監(jiān)測儀用于實(shí)時(shí)監(jiān)測河道水位及水溫變化。地表溫度則是通過埋設(shè)在河岸帶淺層土體中的溫度傳感器進(jìn)行記錄。上述溫度與水位數(shù)據(jù)在收集過程中由數(shù)據(jù)記錄儀控制、記錄和存儲，均為每1 h記錄一次。

圖1 水位、溫度監(jiān)測儀器布設(shè)示意Fig.1 Layout diagram of water level and temperature monitoring instruments

圖2 2012年福克斯1號河道水位、水溫和地表溫度時(shí)序資料Fig.2 Time-series data of water level,water temperature and land surface temperature of the Fox 1 River

基于構(gòu)建的河岸帶水熱耦合模型及其在COMSOL軟件中的實(shí)現(xiàn)方法，結(jié)合河岸帶水位和溫度原型觀測資料，可以實(shí)現(xiàn)對考慮土體非均質(zhì)傳熱的河岸帶水熱耦合模型有限元求解。圖3為河岸帶二維模型計(jì)算區(qū)域示意圖。根據(jù)滲透性的不同，計(jì)算區(qū)域大致可分為圖示的3個(gè)區(qū)域。對于河岸帶飽和-非飽和滲流場，模型左邊界(AF)、右邊界(BC)和空氣與土壤接觸界面(CD、EF)為無流動邊界；河水與土壤接觸界面(DE)為變水位邊界條件，是通過將實(shí)測水位以插值函數(shù)的形式定義在模型邊界上；底部邊界(AB)為透水層邊界；滲流場的初始條件是通過將模擬前一時(shí)段的實(shí)測壓力水頭進(jìn)行線性插值施加于計(jì)算域。對于河岸帶溫度場，模型左邊界(AF)、右邊界(BC)和底部邊界(AB)均假定為絕熱邊界；空氣與土壤接觸界面(CD、EF)及河水與土壤接觸界面(DE)分別為地表溫度邊界和變水溫邊界，同樣是通過實(shí)測地表溫度和河道水溫以插值函數(shù)的形式施加于模型邊界上；溫度場的初始條件為模擬前一時(shí)段各區(qū)域平均溫度。在網(wǎng)格劃分方面，采用COMSOL預(yù)定義的三角形網(wǎng)格單元，并將網(wǎng)格單元大小校準(zhǔn)為較細(xì)化的流體動力學(xué)尺寸，共產(chǎn)生7 877個(gè)網(wǎng)格節(jié)點(diǎn)和15 371個(gè)網(wǎng)格單元。

圖3 模型計(jì)算區(qū)域示意Fig.3 Schematic diagram of model calculation area

參考相關(guān)文獻(xiàn)[28-30]，給出河岸帶水熱耦合計(jì)模型算參數(shù)，如表2所示。其中，參數(shù)ks、θr、α、β、n、Com、Cs和Cw參照文獻(xiàn)[28]給出；Cclay、Csand、Csilt和θs根據(jù)文獻(xiàn)[29]中的砂壤土取值；ρb取自文獻(xiàn)[30]中土壤質(zhì)地、孔隙率和粒徑分布相近的土體；橫向、縱向熱彌散度均取0.01 m；水的導(dǎo)熱系數(shù)取0.598 W/(m·℃)；水的密度取1 000 kg/m3。

表2 河岸帶水熱耦合模型計(jì)算參數(shù)Table 2 Calculation parameters for hydrothermal coupling model of riparian zone

3.2 非均質(zhì)傳熱建模有效性驗(yàn)證

為驗(yàn)證利用PDE模塊代替多孔介質(zhì)傳熱模塊計(jì)算河岸帶土體非均質(zhì)傳熱過程的有效性，將前述6種土體有效導(dǎo)熱系數(shù)模型分別考慮至河岸帶水熱耦合模型中，計(jì)算得出不同模型土體有效導(dǎo)熱系數(shù)與體積含水率變化關(guān)系的數(shù)值解，并與相應(yīng)的模型解析解對比，如圖4所示。由圖4可見，不同模型下土體有效導(dǎo)熱系數(shù)與含水率變化關(guān)系的數(shù)值解與解析解吻合度高，并且具有較好的一致性。

此外，從圖4顯示的土體有效導(dǎo)熱系數(shù)與體積含水率的關(guān)系還可以看出，不同飽和狀態(tài)土體所對應(yīng)的有效導(dǎo)熱系數(shù)明顯不同，并且當(dāng)土體處于非飽和狀態(tài)時(shí)(即θ<0.41)，不同模型預(yù)測的有效導(dǎo)熱系數(shù)亦存在顯著差異，表明所構(gòu)建的河岸帶水熱耦合模型能夠考慮因土體處于不同飽和狀態(tài)而導(dǎo)致的非均質(zhì)傳熱問題，并且能夠?qū)崿F(xiàn)對不同土體有效導(dǎo)熱系數(shù)模型的準(zhǔn)確計(jì)算。

圖4 土體有效導(dǎo)熱系數(shù)模型的解析解與數(shù)值解對比Fig.4 Comparison between analytical solution and numerical solution of soil effective thermal conductivity models

3.3 河岸帶水熱耦合模型驗(yàn)證及對比分析

為驗(yàn)證河岸帶水熱耦合模型模擬效果，擬將模擬的滲流場和溫度場結(jié)果與現(xiàn)場監(jiān)測結(jié)果進(jìn)行對比。由于試驗(yàn)過程中河岸帶壓力傳感器發(fā)生故障，未能有效獲得計(jì)算時(shí)段內(nèi)河岸帶監(jiān)測井的水頭數(shù)據(jù)，給直接驗(yàn)證河岸帶水頭變化帶來困難。這里借鑒Ren等[19]及Naranjo和Smith[28]對滲流場的驗(yàn)證方法，即在缺失實(shí)測壓水頭變化數(shù)據(jù)的情況下，通過比較模型計(jì)算的河岸帶地下水滲流速度與實(shí)測河流水位變化規(guī)律的一致性來定性反映模擬的滲流場合理與否。此外，還可以將模擬的河岸帶地下水滲流速度與水動力學(xué)法計(jì)算結(jié)果進(jìn)行比較，進(jìn)而分析滲流場計(jì)算結(jié)果的可靠性。圖5給出了河流水位及河岸帶地下水滲流速度時(shí)序變化曲線。由圖5可見，基于水熱耦合模型模擬得到的河岸帶地下水滲流速度與河流水位變化規(guī)律基本一致，并且模擬的地下水滲流速度與水動力學(xué)法計(jì)算結(jié)果較為吻合。因而，所構(gòu)建的水熱耦合模型能合理反映河岸帶內(nèi)滲流場變化。

圖5 2012年福克斯1號河流水位及河岸帶地下水滲流速度時(shí)序變化曲線Fig.5 Time-series variation curves of river level and riparian groundwater velocity of the Fox 1 River

為對河岸帶溫度場作進(jìn)一步驗(yàn)證，并比較不同有效導(dǎo)熱系數(shù)模型下河岸帶水熱耦合模型的模擬效果，圖6給出了不同模型下河岸帶各測點(diǎn)溫度模擬值與實(shí)測值對比。由圖6可見，不同有效導(dǎo)熱系數(shù)模型下河岸帶水熱耦合模型模擬的各測點(diǎn)溫度變化效果存在差異，其中，Campbell有效導(dǎo)熱系數(shù)模型下的河岸帶水熱耦合模型模擬的各測點(diǎn)溫度值明顯低于實(shí)測結(jié)果，表明該模型對河岸帶土體傳熱效果預(yù)測偏低。由圖4顯示的土體有效導(dǎo)熱系數(shù)與體積含水率的關(guān)系可知，在相同體積含水率下，Campbell模型預(yù)測的有效導(dǎo)熱系數(shù)明顯低于其他幾種模型，這是導(dǎo)致該模型對河岸帶溫度變化預(yù)測不準(zhǔn)的最主要原因。與其他土體有效導(dǎo)熱系數(shù)模型相比，Campbell模型相對簡單，所涉及的模型參數(shù)僅包括土體堆積密度、含水率和黏土質(zhì)量分?jǐn)?shù)，而土體有效導(dǎo)熱系數(shù)主要與孔隙度和飽和度有關(guān)，該模型并未考慮孔隙度和飽和度的影響，使得土體有效導(dǎo)熱系數(shù)預(yù)測精度不高。此外，當(dāng)土體有效導(dǎo)熱系數(shù)按飽和狀態(tài)取為固定值時(shí)(Keff=2.01 W/(m·℃))，模擬得到的各測點(diǎn)溫度值與實(shí)際偏差較大。相比之下，Johansen、Cté、Lu、改進(jìn)Cté和改進(jìn)Lu模型均表現(xiàn)出了較好的模擬效果，這些模型均是在考慮孔隙度和飽和度基礎(chǔ)上做出的改進(jìn)，并且考慮的影響因素較Campell模型更加豐富，因此在預(yù)測土體有效導(dǎo)熱系數(shù)時(shí)具有較高的精度，進(jìn)而使得模擬的河岸帶溫度值與實(shí)測值較為接近。由此可以看出，在利用水熱耦合模型模擬河岸帶水熱交換過程時(shí)，選擇合適的土體有效導(dǎo)熱系數(shù)模型至關(guān)重要。

圖6 不同模型下各測點(diǎn)溫度模擬值與實(shí)測值對比Fig.6 Comparison of simulated and measured temperature values at each point under different models

為定量評價(jià)不同有效導(dǎo)熱系數(shù)模型下河岸帶水熱耦合模型的模擬效果，引入均方根誤差(ERMS)、決定系數(shù)(R2)和平均相對誤差(EMR)作為模型性能評價(jià)指標(biāo)[31]。表3給出了不同模型下河岸帶各測點(diǎn)溫度模擬值與實(shí)測值的對比統(tǒng)計(jì)結(jié)果。由表3可見，Campbell模型溫度模擬結(jié)果的ERMS值為0.57～2.72 ℃，R2值為0.55～0.96，EMR值為3.91%～14.39%，各項(xiàng)性能評價(jià)指標(biāo)均處于相對較低水平，該模型下的水熱耦合模型低估了河岸帶水熱交換作用。當(dāng)不考慮土壤熱性質(zhì)發(fā)生變化時(shí)(即土體有效導(dǎo)熱系數(shù)取固定值)，其溫度模擬結(jié)果的ERMS值為0.74～2.05 ℃，R2值為0.43～0.96，EMR值為3.54%～19.44%，各評價(jià)指標(biāo)結(jié)果表現(xiàn)較差，同樣不能較好地反映河岸帶水熱交換過程。相比之下，其他幾種土體有效導(dǎo)熱系數(shù)模型模擬結(jié)果的各項(xiàng)性能指標(biāo)表現(xiàn)較好，其中，Johansen模型溫度模擬結(jié)果的ERMS值為0.59～1.16 ℃，R2值為0.91～0.95，EMR值為3.35%～7.72%，各性能評價(jià)指標(biāo)表現(xiàn)最好。

表3 不同模型下各測點(diǎn)溫度模擬結(jié)果的評價(jià)指標(biāo)統(tǒng)計(jì)Table 3 Statistics of evaluation indexes of temperature simulation results of each point under different models

由于不同模型模擬的河岸帶溫度變化在不同測點(diǎn)效果表現(xiàn)不一，難以綜合評價(jià)各有效導(dǎo)熱系數(shù)模型的性能表現(xiàn)。因此，將河岸帶6個(gè)監(jiān)測點(diǎn)處的溫度實(shí)測值與模擬值進(jìn)行整體比較(圖7)，得到了反映模型整體性能表現(xiàn)的評價(jià)指標(biāo)統(tǒng)計(jì)結(jié)果(表4)。由圖7和表4可見，在整體分析情況下，Johansen、Cté、Lu、改進(jìn)Cté和改進(jìn)Lu模型的模擬結(jié)果偏差總體較小，這些模型均能較好地反映河岸帶水熱動態(tài)變化過程。相比之下，當(dāng)不考慮土體非均質(zhì)傳熱或采用Campbell有效導(dǎo)熱系數(shù)模型時(shí)，河岸帶水熱耦合模型模擬效果較差，進(jìn)一步表明河岸帶水熱耦合模型的模擬效果與有效導(dǎo)熱系數(shù)的預(yù)測準(zhǔn)確與否密切相關(guān)?？偟膩砜矗贘ohansen模型的河岸帶水熱耦合模型在模擬河岸帶水熱交換過程中效果表現(xiàn)最佳，該模型能較為真實(shí)地反映河岸帶內(nèi)部水熱動態(tài)變化過程。

表4 整體分析下不同模型溫度模擬結(jié)果的評價(jià)指標(biāo)統(tǒng)計(jì)Table 4 Statistics of evaluation indexes of temperature simulation results of different models under global analysis

圖7 不同有效導(dǎo)熱系數(shù)模型下河岸帶溫度模擬值 與實(shí)測值的整體比較Fig.7 Global comparison of measured and simulated temperature values of riparian zone under different models

圖8給出了基于Johansen模型模擬得到的河岸帶體積含水率及溫度分布圖。由圖8可見，河岸帶地下水與河流地表水進(jìn)行側(cè)向潛流交換作用并發(fā)生熱量交換，其溫度場在太陽輻射和地表水入侵雙重影響下重新分布。河岸帶溫度在水平方向上大致可以分為高溫區(qū)、中溫區(qū)和低溫區(qū)，并且在垂直方向上呈現(xiàn)“上暖下涼”的現(xiàn)象，這與李英玉等[7]和劉東升等[13]此前野外試驗(yàn)所揭示規(guī)律一致。

圖8 河岸帶不同時(shí)刻的體積含水率及溫度分布Fig.8 Distribution of water content and temperature in the riparian zone at different time

4 結(jié) 論

有效導(dǎo)熱系數(shù)是影響和決定土體在傳熱過程中溫度分布的重要參數(shù)，對河岸帶水熱交換的模擬有直接影響。在飽和-非飽和滲流及多孔介質(zhì)傳熱理論基礎(chǔ)上，通過引入土體有效導(dǎo)熱系數(shù)模型，構(gòu)建了考慮土體非均質(zhì)傳熱的河岸帶水熱耦合模型，主要結(jié)論如下：

(1) COMSOL軟件在材料屬性、邊界條件及求解器設(shè)置上極具靈活性，在利用地下水流模塊模擬河岸帶水流運(yùn)動的基礎(chǔ)上，通過自定義偏微分方程來實(shí)現(xiàn)考慮土體非均質(zhì)傳熱的河岸帶水熱耦合模型開發(fā)及應(yīng)用，可避免模型開發(fā)的編程求解。

(2) 考慮土體非均質(zhì)傳熱的河岸帶水熱耦合模型能夠?qū)崿F(xiàn)對流場的合理反映，準(zhǔn)確、可靠地模擬河岸帶溫度時(shí)空變化，對刻畫河岸帶水熱動態(tài)變化過程具有優(yōu)越性。

(3) Johansen有效導(dǎo)熱系數(shù)模型是反映河岸帶水熱交換過程中土體非均質(zhì)傳熱的最佳模型，而對于更廣泛的土體有效導(dǎo)熱系數(shù)模型性能表現(xiàn)，需進(jìn)一步分析。

(4) 由于缺失河岸帶監(jiān)測井水頭變化數(shù)據(jù)，目前僅對河岸帶滲流場進(jìn)行了間接驗(yàn)證。另外，因研究流域氣象站點(diǎn)較少，缺乏降雨資料，加上模型的土-水接觸邊界存在一定程度的概化，導(dǎo)致模擬結(jié)果與實(shí)測值仍存在一定偏差，這些都需要在未來研究中進(jìn)一步完善。