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

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

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

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

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

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

1.1 流動(dòng)方程

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

(1)

式中：ρw為水的密度，kg/m3；Cm為容水度，m-1；g為重力加速度，m/s2；Se為土體的相對(duì)飽和度；Ss為彈性貯水率，Pa-1；p為壓力，Pa；t為時(shí)間，s；?為拉普拉斯算子；θ為體積含水率，m3/m3；ks為飽和滲透系數(shù)，m/s；kr(θ)為非飽和帶相對(duì)滲透系數(shù)，m/s，是體積含水率的函數(shù)；μ(T)為水的動(dòng)力黏度，Pa·s，μ(T)=0.000 024 24×10247.8/(T+133.16)，T為溫度， ℃；z為計(jì)算點(diǎn)高程，m；Qm為滲流場(chǎng)源匯項(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ù)，通過(guò)擬合土壤水分特征曲線得到，m=1-1/β。

1.2 熱量交換方程

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

(6)

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

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

其中，水動(dò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ù)測(cè)土體有效導(dǎo)熱系數(shù)的模型眾多，這些模型可分為經(jīng)驗(yàn)?zāi)Ｐ秃屠碚撃Ｐ汀Ｏ噍^于經(jīng)驗(yàn)?zāi)Ｐ?，理論模型往往預(yù)測(cè)精度不高、公式復(fù)雜且部分參數(shù)不易確定，難以直接應(yīng)用[16-17]。在經(jīng)驗(yàn)?zāi)Ｐ头矫?，蘇李君等[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)方法

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

(8)

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

上述自定義系數(shù)可根據(jù)用戶(hù)的需求定義為常數(shù)或者不同類(lèi)型的函數(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主要通過(guò)θ(θ=nSr)建立與孔隙度和飽和度(Sr)的關(guān)系，反映非均質(zhì)參數(shù)的空間分布。因此，需增加預(yù)置儲(chǔ)存模塊用于計(jì)算和儲(chǔ)存每個(gè)時(shí)間步內(nèi)不同空間位置的有效導(dǎo)熱系數(shù)，該過(guò)程可以通過(guò)在模型樹(shù)的“組件”選項(xiàng)下定義局部變量實(shí)現(xiàn)。

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

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

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

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

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

圖2 2012年?？怂?號(hào)河道水位、水溫和地表溫度時(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é)合河岸帶水位和溫度原型觀測(cè)資料，可以實(shí)現(xiàn)對(duì)考慮土體非均質(zhì)傳熱的河岸帶水熱耦合模型有限元求解。圖3為河岸帶二維模型計(jì)算區(qū)域示意圖。根據(jù)滲透性的不同，計(jì)算區(qū)域大致可分為圖示的3個(gè)區(qū)域。對(duì)于河岸帶飽和-非飽和滲流場(chǎng)，模型左邊界(AF)、右邊界(BC)和空氣與土壤接觸界面(CD、EF)為無(wú)流動(dòng)邊界；河水與土壤接觸界面(DE)為變水位邊界條件，是通過(guò)將實(shí)測(cè)水位以插值函數(shù)的形式定義在模型邊界上；底部邊界(AB)為透水層邊界；滲流場(chǎng)的初始條件是通過(guò)將模擬前一時(shí)段的實(shí)測(cè)壓力水頭進(jìn)行線性插值施加于計(jì)算域。對(duì)于河岸帶溫度場(chǎng)，模型左邊界(AF)、右邊界(BC)和底部邊界(AB)均假定為絕熱邊界；空氣與土壤接觸界面(CD、EF)及河水與土壤接觸界面(DE)分別為地表溫度邊界和變水溫邊界，同樣是通過(guò)實(shí)測(cè)地表溫度和河道水溫以插值函數(shù)的形式施加于模型邊界上；溫度場(chǎng)的初始條件為模擬前一時(shí)段各區(qū)域平均溫度。在網(wǎng)格劃分方面，采用COMSOL預(yù)定義的三角形網(wǎng)格單元，并將網(wǎng)格單元大小校準(zhǔn)為較細(xì)化的流體動(dòng)力學(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ì)傳熱過(guò)程的有效性，將前述6種土體有效導(dǎo)熱系數(shù)模型分別考慮至河岸帶水熱耦合模型中，計(jì)算得出不同模型土體有效導(dǎo)熱系數(shù)與體積含水率變化關(guān)系的數(shù)值解，并與相應(yīng)的模型解析解對(duì)比，如圖4所示。由圖4可見(jiàn)，不同模型下土體有效導(dǎo)熱系數(shù)與含水率變化關(guān)系的數(shù)值解與解析解吻合度高，并且具有較好的一致性。

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

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

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

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

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

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

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

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

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

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

表4 整體分析下不同模型溫度模擬結(jié)果的評(píng)價(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í)測(cè)值的整體比較Fig.7 Global comparison of measured and simulated temperature values of riparian zone under different models

圖8給出了基于Johansen模型模擬得到的河岸帶體積含水率及溫度分布圖。由圖8可見(jiàn)，河岸帶地下水與河流地表水進(jìn)行側(cè)向潛流交換作用并發(fā)生熱量交換，其溫度場(chǎng)在太陽(yáng)輻射和地表水入侵雙重影響下重新分布。河岸帶溫度在水平方向上大致可以分為高溫區(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ù)是影響和決定土體在傳熱過(guò)程中溫度分布的重要參數(shù)，對(duì)河岸帶水熱交換的模擬有直接影響。在飽和-非飽和滲流及多孔介質(zhì)傳熱理論基礎(chǔ)上，通過(guò)引入土體有效導(dǎo)熱系數(shù)模型，構(gòu)建了考慮土體非均質(zhì)傳熱的河岸帶水熱耦合模型，主要結(jié)論如下：

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

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

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

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