基于多種土體導(dǎo)熱系數(shù)模型的潛流帶溫度場模擬

任 杰，門利利, 2，陳 波, 3，王大博，陳建琪，王 帆

(1.西安理工大學 省部共建西北旱區(qū)生態(tài)水利國家重點實驗室，陜西 西安710048； 2.中國電建集團 西北勘測設(shè)計研究院有限公司，陜西 西安 710065; 3.四川省水利水電勘測設(shè)計研究院有限公司，四川 成都 610072)

引言

潛流帶“Hyporheic zone”是指河流河床以及河岸帶飽和沉積層，是地表水與地下水相互作用、物質(zhì)和能量交換、生物群落生長和繁殖的區(qū)域[1-3]。類似于氮、磷等污染物直接或間接排放到水體中，會降低水體有效利用且使得水生系統(tǒng)退化[4-6]，而潛流帶中地表水與地下水的交換過程對污染物的凈化、河流生態(tài)健康的調(diào)控起著關(guān)鍵作用。生態(tài)系統(tǒng)的健康程度依賴于完善的生態(tài)系統(tǒng)功能，健康的河流生態(tài)系統(tǒng)不僅是國家生態(tài)安全體系的重要組成部分，也是經(jīng)濟社會可持續(xù)發(fā)展的重要基礎(chǔ)[7-8]。溫度通過影響水的密度、動力粘度以及土壤孔隙度，進而對滲流場造成影響[9-10]。使用水溫作為天然的示蹤劑，研究河床或河岸帶的物理性質(zhì)和變化過程，對河流水生態(tài)系統(tǒng)具有重要的意義[11-12]。

土體導(dǎo)熱系數(shù)是土體熱特性的關(guān)鍵指標之一，是影響土體傳熱過程中溫度分布以及能量變化的主要因素，也是研究流-熱耦合作用、物質(zhì)傳輸、氣體擴散等土壤物理過程的基礎(chǔ)[13-14]。Kersten等[15]于1949年提出土體導(dǎo)熱系數(shù)與容重有關(guān)的經(jīng)驗?zāi)Ｐ?，是最早用來預(yù)測凍土導(dǎo)熱系數(shù)的方法之一。Johansen等[16]于1975年提出了歸一化導(dǎo)熱系數(shù)的概念，它在預(yù)測不同質(zhì)地的土壤中普遍具有較好的性能，已被廣泛應(yīng)用。Chung和Horton[17]于1987年提出一個較為簡單的經(jīng)驗?zāi)Ｐ?，模型只考慮了含水率對導(dǎo)熱系數(shù)的影響。Lu等[18]于2007年在Johansen研究的基礎(chǔ)上將土體大致分為粗質(zhì)土壤和細質(zhì)土壤，引進一個與砂粒含量相關(guān)的經(jīng)驗常數(shù)。Lu和Dong[19]于2015年提出了閉合形式的土體導(dǎo)熱系數(shù)經(jīng)驗方程，將土壤水分特征曲線與熱傳導(dǎo)系數(shù)相聯(lián)系，能夠較為準確的預(yù)測不同類型土壤的熱傳導(dǎo)系。Ren等[20]于2019年考慮有機質(zhì)對導(dǎo)熱系數(shù)的影響，同時考慮土壤粒組成分對形狀因子的影響，在已有模型的基礎(chǔ)上提出一個新的經(jīng)驗?zāi)Ｐ?，相比于以往模型，該模型適用范圍較更廣，擬合效果更好，精度更高。

Witherspoon[21]在20世紀70年代首先提出耦合理論，經(jīng)過數(shù)十年的不斷改進，耦合理論被廣泛應(yīng)用于各個領(lǐng)域。Molina-Giraldo等[22]通過建立地表水補給地下水的流-熱耦合二維模型，指出非飽和帶傳熱對河岸帶的溫度分布具有明顯影響。Ren等[23]基于室內(nèi)水槽試驗，建立飽和-非飽和流-熱耦合模型，研究不同入滲水頭、不同水溫條件下，砂槽溫度的變化規(guī)律，結(jié)果發(fā)現(xiàn)，水槽低溫區(qū)域隨滲入水頭的增大而增大，且水溫對水槽溫度場影響較為顯著，模型模擬結(jié)果與試驗一致。吳志偉等[24]基于Lu模型構(gòu)建淺部地區(qū)溫度場和飽和-非飽和滲流場的耦合模型，發(fā)現(xiàn)溫度場總體上增大時，滲流速度減小，且溫度場較滲流場的變化幅度偏小，表明了滲流場對溫度場的影響程度較大[25]。于友斌[26]通過Comsol軟件構(gòu)建流-熱耦合模型，通過求解得到滲流對凍結(jié)溫度場與滲流場的影響，并將結(jié)果與試驗進行比較，結(jié)果發(fā)現(xiàn)數(shù)值模擬溫度場與實測溫度場規(guī)律相同。土壤導(dǎo)熱系數(shù)作為流-熱耦合理論的重要驅(qū)動因子，是直接決定模型是否成功和準確的關(guān)鍵，雖然不少學者開展了關(guān)于流-熱耦合理論模型的研究，但對于基于土壤導(dǎo)熱系數(shù)與流-熱耦合理論之間相互作用并應(yīng)用于河流地表水-地下水相互作用機理方面鮮有報道。

基于此，本文建立基于HYDRUS軟件的Chung & Horton模型以及基于COMSOL Multiphysics軟件的Lu(2007)模型、Lu(2014)模型和Ren(2019)模型的河床、河岸帶飽和-非飽和流-熱耦合模型，通過比較野外試驗測得的溫度場和4種方法建立的模型計算得到的溫度場進行對比分析，以此驗證構(gòu)建的潛流層飽和-非飽和流熱全耦合模型的精度和效果，并進行了耦合作用下河床、河岸帶的溫度、水位、垂向流速、環(huán)境溫度、水溫之間的相關(guān)性分析，研究影響溫度及其滲流的關(guān)鍵要素，為通過流-熱耦合模型研究地表水和地下水相互作用提供理論依據(jù)。

1 數(shù)學模型

1.1 滲流場模型

飽和-非飽和滲流場采用Richards方程表示[24]：

(1)

式中：ρw為水的密度；Cm為土體容水度；g為重力加速度；Se為土壤的相對飽和度；Ss為彈性貯水率；P為壓強； ▽為拉普拉斯算子；θ為土壤含水量；ks為介質(zhì)飽和滲透系數(shù)；kr(θ)為非飽和帶相對滲透率(0≤kr(θ)≤1)，是含水率θ的函數(shù)；μ(T)為水的動力粘滯系數(shù)，為溫度T的函數(shù)[27]；z為計算點位置高程；Qs為水流源匯項。

采用van Genuchten模型(簡稱VG模型)來描述土壤水分特征曲線[28]：

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

(2)

式中：Se為土壤相對飽和度；θr為殘余含水率；θs為飽和含水率。其中Se為：

(3)

式中：hp為壓力水頭，hp=pw/ρw，在非飽和帶等于負壓水頭Hc；α和nv為VG模型參數(shù)，其中，α為水分特征曲線進氣值的倒數(shù)，nv為水分特征曲線斜率的指示參數(shù)，也是反映土壤孔隙分布的參數(shù)；通過擬合土壤水分特征曲線得到，m= 1-1/nv。

土壤容水度Cm和相對滲透率kr分別采用下述經(jīng)驗公式表示：

(4)

(5)

1.2 溫度場模型

溫度場采用熱對流傳熱方程描述[25]：

(6)

式中：λeq，ρeq，Ceq分別為整個含水層系統(tǒng)的等效導(dǎo)熱系數(shù)，等效密度以及等效比熱容，T為地下水溫度，Cw、ρw分別為流體的密度和比熱容，u為滲流速度，QG為熱量源匯項。

飽和-非飽和土體等效密度以及比熱容遵循體積平均定律，其表達式為：

ρeq=(1-n)ρs+θρw+(n-θ)ρg.

(7)

Ceq=(1-n)Cs+θCw+(n-θ)Cg.

(8)

式中：n為多孔介質(zhì)的孔隙度，當土體處于飽和狀態(tài)時，即θ=n；ρs、ρw、ρg分別為土、水和空氣的密度，Cs、Cw、Cg分別為土、水和空氣的比熱容。

1.3 土體導(dǎo)熱系數(shù)模型

1.3.1 Chung & Horton(1987)模型

Chung和 Horton在1987年提出一個較為簡單的經(jīng)驗?zāi)Ｐ?，模型只考慮了含水率對導(dǎo)熱系數(shù)的影響，因此在土壤其他物理參數(shù)發(fā)生變化的情況下，模擬誤差就會較大，模型公式如下[17]：

λ=b1+b2·θ+b3·θ0.5.

(9)

式中：b1，b2，b3為經(jīng)驗常數(shù)，其值與土壤的類型有關(guān)。

1.3.2 Lu(2007)模型

Lu和Ren在Johansen研究的基礎(chǔ)上，采用Johansen的歸一化模型，提出一個比較適合低含水率條件下土壤導(dǎo)熱系數(shù)的計算模型，但因其將石英含量等量為砂的含量，因此不適合用來計算砂土含量較高的土壤導(dǎo)熱系數(shù)，計算值較實測值偏低[20]，模型公式如下[18]：

λ=(λsat-λdry)·Ke+λdry.

(10)

(11)

式中：λdry和λsat分別為干土和飽和土體導(dǎo)熱系數(shù)，單位W·m-1·K-1；Ke為Kersten系數(shù)；1.33指的是形狀參數(shù)；m是由砂粒含量確定的經(jīng)驗參數(shù)，對于砂粒含量大于40%的粗質(zhì)土壤和砂粒含量小于40%的細質(zhì)土壤，m分別取值為：0.96、0.27。λdry通過與孔隙度n線性關(guān)系擬合得到：

λdry=-a·n+b.

(12)

式中：a，b為經(jīng)驗系數(shù)；n為孔隙度，%。在0.2

飽和土體導(dǎo)熱系數(shù)λsat的計算公式如下：

(13)

1.3.3 Lu(2014)模型

Lu(2014)模型在Lu(2007)模型的基礎(chǔ)上，進一步證實了土壤顆粒成分對λ的影響，模型公式簡單，參數(shù)易于獲取，但是對于干燥土體的導(dǎo)熱系數(shù)只是粗略地通過線性擬合得到，沒有全面考慮不同土壤質(zhì)地，從而產(chǎn)生誤差，模型精度有待進一步提高。模型公式如下[19]：

λ=λdry+exp(φ-θ-τ),θ>0.

(14)

式中：τ和φ是與砂粒含量、黏土含量、體積密度有關(guān)的λ(θ)曲線的形狀因子。λdry采用Lu(2007)的關(guān)于n的線性方程(1-12)計算。

τ與黏粒質(zhì)量分數(shù)Cclay滿足如下關(guān)系：

τ=0.67Cclay+0.24.

(15)

參數(shù)φ具有以下多元回歸方程，其值由砂土質(zhì)量分數(shù)Csand和體積密度ρb確定：

φ=1.97Csand+1.87ρb-1.36Csandρb-0.95.

(16)

1.3.4 Ren(2019)模型

為了改變早期模型對土質(zhì)的粗略劃分，以及參數(shù)的不確定性，Ren等將早期模糊的參數(shù)取值改進為輸入土體本身所擁有的精確物理參數(shù)。其次，考慮有機質(zhì)含量對土體導(dǎo)熱系數(shù)的影響，并且考慮土壤質(zhì)地和顆粒組成對形狀因子(表征導(dǎo)熱系數(shù)增長速率和曲線斜率)的影響。在以往模型的基礎(chǔ)上，采用Johansen(1975)的歸一化模型，提出了在整個含水率范圍內(nèi)的新導(dǎo)熱系數(shù)經(jīng)驗?zāi)Ｐ蚚20]：

λ=(λsat-λdry)·Ke+λdry.

(17)

λsat=0.53Csand+0.1γd.

(18)

λdry=0.087Csand+0.019γd.

(19)

λsat和λdry與砂土含量和干土容重(γd)有關(guān)，當砂含量0

γd=g·ρb

(20)

式中：g= 9.8 m2·s-1，是重力加速度。

Ke與θ的新關(guān)系如下：

Ke=exp(α-θ-β).

(21)

α和β是λ(θ)曲線的形狀因子，其值會影響到導(dǎo)熱系數(shù)曲線的斜率和增長速度，表達式如下：

α=0.493Csand+0.86Csilt+0.014Com+0.778.

(22)

β=0.736Cclay+0.006Com+0.222.

(23)

式中：Cclay、Csilt、Csand分別為黏粒、粉粒、砂粒的質(zhì)量分數(shù)(%)；Com為有機質(zhì)質(zhì)量比(g·kg-1)。

1.4 模型評估

采用均方根誤差(root-mean-square error，RMSE)、決定系數(shù)(coefficient of determination，R2)和相對誤差(relative error，Re)對模型的模擬精度進行評估。

(24)

(25)

(26)

2 實例應(yīng)用

2.1 研究區(qū)域概況

引用美國地質(zhì)勘探局于2012年3月至2013年10月在沃克河流域的實地研究項目，選取梅森谷的福克斯灌渠2號監(jiān)測點的資料進行分析[31]，該監(jiān)測點位于38°54′06″N，119°09′30″E，海拔約為 1 354.45 m。渠道水位和水溫分別采用Micro-Diver型壓力傳感器和溫度自動監(jiān)測儀進行測量，測量精度分別為±1.0 cm和±0.1℃，分辨率分別為0.2 cm和0.01℃。各監(jiān)測井分別位于河道中心位置0 m、2.8 m和5.8 m處。各監(jiān)測井傳感器距離地表位置分別為T6(0.1 m、0.2 m、0.5 m、0.75 m、1 m)，BP1(0.50 m、0.70 m、1.50 m、2.15 m)，BP2(1.20 m、1.70 m、2.30 m、2.75 m)，各數(shù)據(jù)每間隔1 h記錄一次。

2.2 模型邊界條件與參數(shù)設(shè)置

2.2.1 邊界條件

本研究選取?？怂构嗲?012年4月1日至2012年10月26日的溫度和水位實測數(shù)據(jù)進行分析，依據(jù)土壤物理特性將研究區(qū)域劃分為區(qū)域1和區(qū)域2。區(qū)域1和區(qū)域2的初始溫度取監(jiān)測數(shù)據(jù)第一個星期的平均溫度，分別為6.6 ℃和6.25 ℃。

在COMSOL Multiphysics建模中，模型邊界條件及尺寸如圖1(a)所示。對于溫度場，將兩邊AB、CD以及底部BC均設(shè)置為絕熱邊界，河道中與水流接觸的EF為水溫波動邊界，最高水位之上與大氣接觸的AF和ED設(shè)置為氣溫波動邊界，由于大氣溫度及水溫是隨時間變化的溫度序列，故采用三次樣條插值函數(shù)施加在邊界上，圖2是2012年4月1日至2012年10月26日的氣溫與水溫實測數(shù)據(jù)動態(tài)變化圖。對于滲流場，水位之上AF、ED以及斷面左右兩邊AB邊界和CD邊界均為無流動邊界，底部BC為透水邊界，EF設(shè)置為變水位邊界，同樣采用三次樣條插值函數(shù)進行施加，河道水位數(shù)據(jù)如圖3所示。

圖1 模型計算區(qū)域及邊界條件(單位：m)Fig.1 Model calculation area and boundary conditions

圖2 氣溫與水溫動態(tài)變化圖Fig.2 Dynamic changes intemperature and water temperature

圖3 河道水位時程曲線圖Fig.3 Time-history chart of river water level

相比于COMSOL Multiphysics，在HYDRUS建模中其滲流場邊界條件設(shè)置與COMSOL Multiphysics軟件相同，這里不再贅述。熱傳導(dǎo)邊界對最高水位至河床表面EF設(shè)置為第三類邊界(Cauchy邊界)，即就是隨時間變化的水溫邊界，這與COMSOL Multiphysics軟件中的水溫邊界一致，都是插入實際水溫值；AF和ED設(shè)置為第一類邊界條件(Dirichlet邊界)，即就是指定溫度邊界，可插入實際大氣溫度；左右兩邊AB、CD以及斷面底部BC設(shè)置為第二類邊界條件(Neumann邊界)，即就是COMSOL Multiphysics軟件中的絕熱邊界，具體設(shè)置如圖1(b)所示。

2.2.2 模型參數(shù)

結(jié)合現(xiàn)場勘查資料并參考[31]野外試驗相關(guān)數(shù)據(jù)，給出了飽和-非飽和流-熱耦合模型滲流場及溫度場的計算參數(shù)，其中滲透系數(shù)采用室內(nèi)試驗中的常水頭法進行測量，土壤孔隙度通過環(huán)刀法測量，殘余含水率通過烘干法測量，土壤粒徑級配(PSD)采用移液管法測定[32]，有機質(zhì)含量采用Walkley-Black滴定法測量[33]，水的比熱容為默認值，各參數(shù)的具體數(shù)值如表1所示。

表1 河床、河岸帶飽和-非飽和流-熱耦合模型計算參數(shù)Table1 Calculation parameters of saturated-unsaturated hydro-thermal coupling model for riverbed and riparian zone

續(xù)表

3 結(jié)果與討論

3.1 模型對比分析

圖4和圖5分別為各模型對3個監(jiān)測井BP1、BP2和T6的13個監(jiān)測點在2012年4月1日至2012年10月26日的5 000個小時的溫度擬合效果圖，表2為各模型對各監(jiān)測點的擬合誤差結(jié)果。

圖4 各模型對BP1和BP2監(jiān)測井的擬合溫度曲線圖Fig.4 The fitting temperature curve of BP1 and BP2 monitoring tube by each model

圖5 各模型對T6監(jiān)測井的擬合溫度曲線圖Fig.5 The fitting temperature curve of T6 monitoring tube by each model

表2 各模型對各監(jiān)測點的擬合誤差表Table2 The fitting error table of each model for each monitoring point

由圖4、圖5和表2可以看出，對于河道中心的T6監(jiān)測井，越靠近河床表面，溫度的正弦日變化越明顯，溫度波動越大，隨著深度的增加，溫度日變化減小，呈季節(jié)性周期變化。這可能是因為越靠近河床表面，其受水位變化以及大氣溫度變化的影響越明顯。相比于T6監(jiān)測井，距離河道中心2.8 m處的BP1監(jiān)測井和距離河道中心5.8 m處的BP2監(jiān)測井各點溫度變化較為穩(wěn)定，溫度日變化不明顯。

從擬合效果可以看出，相比于其他模型，基于COMSOL Multiphysics軟件的Ren(2019)模型更適合應(yīng)用于河床以及河岸帶飽和-非飽和流-熱耦合研究中，整體擬合精度較高，適用范圍更廣，能夠較為精確地刻畫飽和-非飽和溫度場的動態(tài)變化過程。其均方根誤差RMSE的變化范圍在0.824～2.691°C之間，均值為1.539℃，明顯小于其他模型，決定系數(shù)R2的變化范圍在0.589～0.974之間，均值為0.851，相對誤差Re變化范圍在4.5%～15.5%之間，均值為8.6%，均優(yōu)于其他模型；其次是Chung & Horton(1987)模型，均方根誤差RMSE均值為2.010°C，決定系數(shù)R2均值為0.764，相對誤差Re均值為11.1%，Lu(2007)和Lu(2014)模型擬合效果一般。

這是因為有研究表明有機質(zhì)含量也會影響土體的導(dǎo)熱系數(shù)[34]，有機質(zhì)含量增加了土體的疏水性[35]，即就是改變了土壤中的水分分布，從而影響土的導(dǎo)熱系數(shù)。而Ren(2019)模型相比于其他模型考慮因素更加全面，其不僅考慮了有機質(zhì)含量對導(dǎo)熱系數(shù)的影響，同時也考慮了土壤質(zhì)地對用來表征導(dǎo)熱系數(shù)增長速度和曲線斜率的形狀因子的影響；相比于Ren(2019)模型，Chung & Horton(1987)模型只考慮了含水率對導(dǎo)熱系數(shù)的影響，因此在土壤其他物理參數(shù)發(fā)生變化的情況下，模擬誤差就會較大；Lu(2007)將土體大致分為粗質(zhì)土壤和細質(zhì)土壤，引進一個與砂粒含量相關(guān)的經(jīng)驗常數(shù)，但由于沒有定量或定性的區(qū)分土質(zhì)，導(dǎo)致計算時對該常數(shù)取值不易確定，從而產(chǎn)生偏差，Lu(2014)模型[36]對于干燥土體的導(dǎo)熱系數(shù)只是粗略地通過線性擬合得到，沒有全面考慮不同土壤質(zhì)地，因此誤差仍然很大。

3.2 河床、河岸帶溫度場的時空分布規(guī)律

為研究溫度在河床以及河岸帶剖面的空間分布規(guī)律及動態(tài)變化特征，根據(jù)上述已驗證的河床及河岸帶飽和-非飽和流-熱耦合模型，通過比較2012年 6 月 22 日和2012年10月5日不同時刻的溫度云圖，發(fā)現(xiàn)在該時間點研究區(qū)域分別符合季“上暖下涼”，冬季“上涼下暖”的溫度時空分布規(guī)律，為了進一步表明變化的趨勢和規(guī)律，選取6月22日和10月5日為典型時間點進行分析。在t= 2:00、t= 8:00、t= 14:00和t= 20:00四個不同時刻模擬數(shù)據(jù)，繪制不同季節(jié)河流橫斷面的溫度日變化云圖，如圖6所示。所選區(qū)域從地表面開始，同時包含了3個監(jiān)測井的所有監(jiān)測點，橫軸表示距離BP2監(jiān)測井的距離，縱軸表示高程。

圖6 不同季節(jié)河岸帶溫度變化特征Fig.6 Characteristics of temperature change in riparian zone in different seasons

由圖6(a)可知：在夏季總體上呈現(xiàn)表層 > 中層 > 底層的趨勢，河床、河岸帶溫度大致可分為低溫區(qū)、中溫區(qū)和高溫區(qū)，即在遠離河道的淺層區(qū)域溫度較高，河岸帶深層區(qū)域溫度較低。通過比較6月22日不同時刻的溫度云圖發(fā)現(xiàn)：1) 在t= 2:00至t= 08:00時段內(nèi)，氣溫逐漸上升，使得河岸帶地表溫度隨之升高，而河床溫度相對較低，這種現(xiàn)象同樣出現(xiàn)在Ren等[37-38]的研究中，這是由于此階段河流水位上升，達到一年內(nèi)最高水位，水溫相對較低導(dǎo)致；2)在t= 8:00至t= 14:00時段內(nèi)，太陽輻射作用逐漸增強，河岸帶沉積物因比熱容相對較低，表層溫度迅速上升，與河道水溫形成高溫度差，河岸帶內(nèi)的垂向溫度分層明顯，表層溫度明顯高于中層和底層；3)當t= 20:00時，整個研究區(qū)域低溫區(qū)面積逐漸增大，高溫區(qū)面積開始減小，由于空氣和水二者熱容相差較大，因此河岸帶沉積層溫度降低較快，而河床部分溫度相對穩(wěn)定。

由圖6(b)可以看出：在冬季總體上呈現(xiàn)表層 < 中層 < 底層的趨勢，河床、河岸帶溫度大致也可以分為低溫區(qū)、中溫區(qū)和高溫區(qū)，然而與夏季相反，低溫和中溫區(qū)域主要保持在河岸帶淺層區(qū)域，而河床內(nèi)部與河岸帶深層區(qū)域溫度分布相對均勻。通過比較10月5日的四個不同時刻發(fā)現(xiàn)：1) 在t= 2:00至t= 8:00，河床、河岸帶溫度分布相對均勻，受夜間溫度影響，河岸帶表層溫度較低，而河床區(qū)域溫度相對較高且比較穩(wěn)定；2)在t= 8:00至t= 14:00，整個研究區(qū)域溫度較之前升高，變化最為明顯的是河岸帶表層區(qū)域；3)當t= 20:00時，由于環(huán)境溫度的逐漸減小，河岸帶淺層溫度逐漸降低，但河道水溫因滯后以及比熱容較大的原因，并未立即降低，整體溫度場分布逐漸恢復(fù)至t= 2:00時的狀態(tài)。

通過上述分析可以看出，季節(jié)性溫度變化以及水位對河床、河岸帶的溫度場影響較大，在同一季節(jié)，河岸帶表層沉積物的溫度變化較為顯著，深層區(qū)域相對穩(wěn)定，這是因為超過地表以下1 m深度的土壤對大氣溫度以及太陽輻射不太敏感，溫度的年波動最多可延伸至地表以下10 m[39]，此外Engelhardt等[40]研究表明，晝夜溫度變化作用對較深層的潛流帶溫度作用不明顯，這一點也被Lee等[41]的研究所證實。因此更加說明利用COMSOL Multiphysics軟件建立的基于Ren(2019)模型的飽和-非飽和流-熱耦合模型是可靠的。

3.3 流速空間變化規(guī)律分析

基于上述驗證后的飽和-非飽和流-熱耦合模型，反演推算出流速場信息。圖7為河床以及河岸帶內(nèi)部流速場隨水位波動發(fā)生規(guī)律性變化圖。文中給出T6，BP1，BP2監(jiān)測井中兩兩監(jiān)測點中心位置的垂向流速，即T6-(a，b，c，d)、BP1-(a，b，c)和BP2-(a，b，c)。當水位上升時，負值表示地表水補給地下水，正值表示地下水補給地表水。由圖7可以看出，不同監(jiān)測井不同深度處的流速在時間和空間上都有很大的差異。對位于河岸帶上的BP1和BP2監(jiān)測井，監(jiān)測點越靠近河床表面，流速波動范圍越小，流速值隨水位上升增長較慢，且流速較小；水位下降，流速朝相反方向變化，地下水補給地表水，但補給速度仍然較小，整個監(jiān)測過程中流速有正有負(BP1-a和BP2-a)。而隨深度的增加，流速波動范圍增大，隨水位上升，流速迅速朝負方向增大，流速值較大，且流速主要以負值為主，呈現(xiàn)出地表水補給地下水的現(xiàn)象。

圖7 各監(jiān)測井不同深度流速與水位關(guān)系圖Fig.7 Diagram of relationship between flow velocity and water level at different depths of each monitoring tube

而位于河床中心的T6監(jiān)測井則呈現(xiàn)出與河岸帶監(jiān)測點相反的趨勢。具體表現(xiàn)為，離河床越近，流速隨水位波動就越大，流速值越大，隨著監(jiān)測點深度的增加，流速波動逐漸減緩。但與BP1和BP2監(jiān)測井同樣，在淺層區(qū)域，流速有正有負；隨深度增加，流速主要以負方向為主，即主要為地表水補給地下水。

3.4 相關(guān)性分析

表3是3個監(jiān)測井各自不同影響因素之間的相關(guān)矩陣與顯著性分析，由表3可以看出：

(1)BP1監(jiān)測井淺層區(qū)域溫差(ΔTBP1-a，ΔTBP1-b)與BP1各監(jiān)測點溫度呈弱負相關(guān)，而深層區(qū)域溫差(ΔTBP1-c)與各監(jiān)測點溫度近似呈中度程度正相關(guān)，對于BP2監(jiān)測井，其淺層區(qū)域溫差(ΔTBP2-a)與BP2各監(jiān)測點溫度同樣呈弱負相關(guān)，隨深度增加，轉(zhuǎn)化為正相關(guān)，達到弱正相關(guān)和中度程度正相關(guān)的范圍。而對于T6監(jiān)測井，可明顯看出，其各區(qū)域溫差與各監(jiān)測點溫度均呈正相關(guān)，因此可以推斷，各區(qū)域溫差與各監(jiān)測點的相關(guān)性可能受地下水位影響，與該區(qū)域距離河床或者地下水位的距離有關(guān)。

(2)環(huán)境溫度與BP1、BP2監(jiān)測井的各監(jiān)測點溫度均為弱負相關(guān)，而與T6監(jiān)測井的各監(jiān)測點溫度隨深度呈消減式弱正相關(guān)，這是因為，監(jiān)測點埋置在地表以下，即使氣溫波動較大，但因為“衰減”和“滯后”效應(yīng)，位于河岸帶內(nèi)部的監(jiān)測點溫度波動幅度仍相對較小，所以呈現(xiàn)出弱相關(guān)，這也體現(xiàn)了河床、河岸帶夏季“上暖夏涼”，冬季“上涼下暖”的溫度分層現(xiàn)象。

(3)水溫與T6監(jiān)測井各監(jiān)測點溫度呈強正相關(guān)，與BP1各監(jiān)測點溫度呈中度程度正相關(guān)，與BP2各監(jiān)測點呈弱正相關(guān)，由此可以看出，距離河床越近，其溫度場受河水影響越大。

(4)水位對距離地表較近的淺層區(qū)域流速影響較小，而對深層區(qū)域流速影響較大，呈強負相關(guān)；水位升高，地表水開始補給地下水，流速逐漸減小并朝負方向增大。

(5)水位與河岸帶區(qū)域的BP1各監(jiān)測點的溫度均呈中度程度負相關(guān)，與BP2各監(jiān)測點溫度同樣呈中度程度負相關(guān)，與T6各監(jiān)測點的溫度呈弱負相關(guān)，相比于河床附近的監(jiān)測點，水位對距離河床較遠的監(jiān)測點的溫度影響更大，這與水溫對各監(jiān) 測井的溫度影響恰恰相反。

(6)對于河岸帶淺層區(qū)域，垂向流速與流速點所對應(yīng)的區(qū)域溫差呈弱負相關(guān)，隨深度的增加，垂向流速與溫差逐漸呈正相關(guān)，且相關(guān)性逐漸增大。而對于河床中心的T6監(jiān)測井，其各點垂向流速與流速點所對應(yīng)的區(qū)間溫差呈弱負相關(guān)，由此可以看出，溫度梯度會引起水體流動，水體流動又會導(dǎo)致溫度變化，且該現(xiàn)象對于河岸帶深層區(qū)域較為明顯。

(7)對于河岸帶區(qū)域，除了離地表較近的監(jiān)測點以外，隨著深度的增加，流速與監(jiān)測點溫度的相關(guān)性由弱正相關(guān)增強到中度程度正相關(guān)，隨著地下水流速的增大，各監(jiān)測點溫度也隨之增大。而對于T6監(jiān)測井，其各中心點流速與其他各監(jiān)測點溫度幾乎都呈弱正相關(guān)，相比于河岸帶，流速的變化對各監(jiān)測點溫度的影響較小。另外，從表3還可以看出，環(huán)境溫度和水溫與 BP1、BP2和T6各監(jiān)測點流速均呈弱負相關(guān)，有一定的抑制作用。

表3 各監(jiān)測井中各影響因素之間的相關(guān)矩陣與顯著性分析Table3 Correlation matrix and significance analysis of each influencing factor in each monitoring tube

(a)

(b)

4 結(jié)論

本文通過運用兩款商業(yè)軟件嵌入四種土體導(dǎo)熱系數(shù)模型，基于內(nèi)華達州梅森谷的?？怂构嗲?號監(jiān)測點的數(shù)據(jù)資料，分析了不同模型對河流潛流帶溫度場的模擬效果：

(1)相比于基于HYDRUS軟件的Chung & Horton(1987)模型以及基于COMSOL Multiphysics軟件的Lu(2007)、Lu(2014)模型來說，基于COMSOL Multiphysics軟件的Ren(2019)模型更適合應(yīng)用于河流潛流帶的飽和-非飽和流-熱耦合研究中，模型考慮因素更加全面，使得整體擬合精度較高，適用范圍更廣，能夠較為精確地刻畫飽和-非飽和溫度場的動態(tài)變化過程。

(2)水溫以及地表以下各監(jiān)測井的溫度相比于環(huán)境溫度具有明顯的“衰減”和“滯后”現(xiàn)象，環(huán)境溫度的日變化對河床區(qū)域的溫度場分布影響較大，而對河岸帶區(qū)域影響較小，距離河床越遠，埋置深度越深，受大氣溫度與水溫的影響就越小，溫度場變化越穩(wěn)定。

(3)對于河岸帶區(qū)域，隨著深度增加，流速與其上下兩個監(jiān)測點以及兩點溫差逐漸由極顯著弱相關(guān)變?yōu)闃O顯著中度程度正相關(guān)，與水位呈極顯著負強相關(guān)，且越靠近河床，其相關(guān)系數(shù)越大，即隨深度增加河岸帶受側(cè)向潛流交換以及地下水影響增大；且各監(jiān)測點溫度與水位呈極顯著中度程度負相關(guān)。對于河床區(qū)域，各測點溫度受水溫影響較大，呈極顯著強相關(guān)；流速與水位呈極顯著極強負相關(guān)。