潘夢(mèng)綺，黃權(quán)中，馮 榕，黃冠華

·農(nóng)業(yè)水土工程·

采用不同監(jiān)測(cè)數(shù)據(jù)組合反演飽和均質(zhì)石英砂水熱運(yùn)移參數(shù)

潘夢(mèng)綺，黃權(quán)中※，馮 榕，黃冠華

（1. 中國(guó)農(nóng)業(yè)大學(xué)中國(guó)農(nóng)業(yè)水問題研究中心，北京 100083；2. 中國(guó)農(nóng)業(yè)大學(xué)中國(guó)-以色列國(guó)際農(nóng)業(yè)研究培訓(xùn)中心，北京 100083）

土壤及含水層的水力參數(shù)與熱參數(shù)對(duì)于定量描述土壤水、地下水遷移規(guī)律及其伴生的熱運(yùn)移過程十分重要。為探討不同監(jiān)測(cè)數(shù)據(jù)類型組合對(duì)多孔介質(zhì)水熱參數(shù)估計(jì)的影響，該研究基于熱示蹤方法，開展了3種不同粒徑條件下的飽和均質(zhì)石英砂的熱示蹤試驗(yàn)，并結(jié)合HYDRU-2D模型對(duì)介質(zhì)的飽和導(dǎo)水率、導(dǎo)熱系數(shù)和縱向、橫向熱彌散度進(jìn)行反演。參數(shù)估計(jì)時(shí)分別設(shè)置3種情景對(duì)介質(zhì)水熱參數(shù)進(jìn)行估計(jì)：僅采用觀測(cè)點(diǎn)溫度（R1）、觀測(cè)點(diǎn)溫度+水流通量（R2）、觀測(cè)點(diǎn)溫度+水流通量+熱量損失（R3）。并對(duì)R1情景設(shè)置3種不同參數(shù)反演組合，即同時(shí)對(duì)2組參數(shù)（飽和導(dǎo)水率和導(dǎo)熱系數(shù)）、3組參數(shù)（飽和導(dǎo)水率、縱向和橫向熱彌散度）和4組參數(shù)（飽和導(dǎo)水率、導(dǎo)熱系數(shù)、縱向和橫向熱彌散度）進(jìn)行估計(jì)。研究結(jié)果表明：同時(shí)對(duì)介質(zhì)飽和導(dǎo)水率、導(dǎo)熱系數(shù)與熱彌散度進(jìn)行估計(jì)，有利于提高介質(zhì)水熱參數(shù)的估計(jì)精度；對(duì)導(dǎo)熱系數(shù)的合理估計(jì)可減小R1情景中介質(zhì)飽和導(dǎo)水率的估計(jì)誤差。4組參數(shù)中飽和導(dǎo)水率是敏感性最高的參數(shù)，增加用于參數(shù)反演的水流運(yùn)動(dòng)和熱量傳遞信息時(shí)，粗砂、中砂、細(xì)砂的累積流量相對(duì)誤差分別減少了9.74、6.65和12.53個(gè)百分點(diǎn)，顯著提高了介質(zhì)飽和導(dǎo)水率的反演精度。飽和導(dǎo)水率的估計(jì)值隨介質(zhì)粒徑增大而增大，而縱向熱彌散度隨粒徑的變化則呈相反的變化規(guī)律，橫向熱彌散度估值基本不變。增加水流和熱量傳遞信息還能顯著提高中砂的導(dǎo)熱系數(shù)反演精度，導(dǎo)熱系數(shù)的估計(jì)值隨著介質(zhì)孔隙度增大而逐漸降低。研究可為基于不同數(shù)據(jù)類型的均質(zhì)介質(zhì)參數(shù)反演提供。

數(shù)據(jù)；參數(shù)；熱示蹤；飽和石英砂；水熱運(yùn)移模型；HYDRUS軟件

0 引 言

土壤和地下含水層中的水分與熱量遷移在農(nóng)業(yè)生產(chǎn)、地下水污染評(píng)價(jià)等方面具有重要意義[1]。其中，多孔介質(zhì)水力參數(shù)與熱參數(shù)的確定，對(duì)于定量描述土壤水、地下水運(yùn)移規(guī)律及其伴生的熱遷移規(guī)律十分重要。在工程實(shí)際中，滲流介質(zhì)和邊界條件往往都比較復(fù)雜，利用解析公式求解滲透系數(shù)或滲透張量則更加困難，因此數(shù)值模擬方法逐漸成為了研究飽和/非飽和多孔介質(zhì)中水熱遷移規(guī)律的重要手段。

傳統(tǒng)測(cè)量滲透系數(shù)的方法有經(jīng)驗(yàn)公式法、室內(nèi)試驗(yàn)法、野外抽水實(shí)驗(yàn)、壓水實(shí)驗(yàn)法、數(shù)學(xué)模型反演計(jì)算法等[2-6]。溫度采集設(shè)備的智能化發(fā)展使得利用溫度數(shù)據(jù)來補(bǔ)充和替代測(cè)壓管水頭并分析多孔介質(zhì)流場(chǎng)及其水力特性成為可能。大量的研究表明，熱示蹤可在一定程度上較準(zhǔn)確對(duì)介質(zhì)的水、熱參數(shù)進(jìn)行反演[7-11]。并且在利用溫度數(shù)據(jù)反演參數(shù)的基礎(chǔ)上同時(shí)考慮其他不同類型的數(shù)據(jù)信息，有利于提高參數(shù)估計(jì)的精度。Nakhaei等[12]采用HYDRUS模型對(duì)土壤溫水入滲過程進(jìn)行模擬，結(jié)果表明利用多種反演數(shù)據(jù)可以獲得具有較高精度的水熱估計(jì)參數(shù)。Giambastiani等[13]則將相同水力條件下的溶質(zhì)遷移和熱示蹤試驗(yàn)相結(jié)合，達(dá)到了準(zhǔn)確估計(jì)縱向熱彌散的目的。近年來，為獲得更準(zhǔn)確的非均質(zhì)介質(zhì)的滲透系數(shù)分布，斷層掃描也逐步應(yīng)用于熱示蹤試驗(yàn)中[14-15]。

本文采用3種不同顆粒尺寸的石英砂分別開展了飽和穩(wěn)定流場(chǎng)條件下的均質(zhì)介質(zhì)的熱示蹤實(shí)驗(yàn)，并結(jié)合HYDRU-2D模型對(duì)介質(zhì)的水熱參數(shù)進(jìn)行反演。參數(shù)估計(jì)時(shí)設(shè)置了3種不同類型數(shù)據(jù)的反演情景，利用不同類型的數(shù)據(jù)信息對(duì)介質(zhì)的水熱參數(shù)進(jìn)行估計(jì)，并分析數(shù)據(jù)類型對(duì)參數(shù)估計(jì)的影響，從而獲得基于熱示蹤的飽和均質(zhì)介質(zhì)水、熱參數(shù)反演方法。

1 材料與方法

1.1 試驗(yàn)裝置

本文構(gòu)建了基于飽和多孔介質(zhì)水熱運(yùn)移過程的熱示蹤試驗(yàn)平臺(tái)，分別由進(jìn)水裝置、有機(jī)玻璃砂箱、數(shù)據(jù)采集裝置及出水裝置4部分構(gòu)成，如圖1所示。

圖1 熱示蹤試驗(yàn)裝置示意圖

有機(jī)玻璃砂箱的長(zhǎng)、寬、高分別為50、40、2.5 cm。砂箱左側(cè)中部為進(jìn)水室，長(zhǎng)為10 cm；右側(cè)全部邊界為出水室，長(zhǎng)為40 cm。分別在進(jìn)、出水室和石英砂介質(zhì)間設(shè)置厚度為5 cm的多孔隔板和反濾層，從而防止細(xì)小的石英砂堵塞多孔隔板。在砂箱背面均勻設(shè)置20個(gè)熱電偶觀測(cè)點(diǎn)，分別在進(jìn)、出水室底部設(shè)置2個(gè)熱電偶觀測(cè)點(diǎn)。T型熱電偶線（美國(guó)Omega工程公司，0.01 ℃）一端埋設(shè)于介質(zhì)內(nèi)部，一端與CR3000型數(shù)據(jù)采集器（美國(guó)Campbell科技公司，3次/min）相連。進(jìn)水裝置由溢流水槽、加熱水箱及微型水泵組成。溢流水槽連接砂箱進(jìn)水室為試驗(yàn)提供恒定水頭；采用控溫加熱器為試驗(yàn)提供恒溫?zé)嵩?；出水室連接排水水箱，并且進(jìn)、出水室均安裝壓力計(jì)。

試驗(yàn)采用3種不同粒徑的白色精制工業(yè)石英砂，砂箱從進(jìn)水口到出水口分層填裝，介質(zhì)的飽和含水率與導(dǎo)熱系數(shù)均通過一維均質(zhì)砂柱試驗(yàn)獲得，如表1所示。

表1 石英砂的基本物理性質(zhì)

1.2 試驗(yàn)方法

試驗(yàn)開始前在恒定水頭條件下先注入常溫水（23 ℃），當(dāng)出水流量穩(wěn)定和各觀測(cè)點(diǎn)溫度不再變化后，保持水頭不變，將恒溫水流（45 ℃）通過溢流水槽立即注入進(jìn)水室。試驗(yàn)采用熱電偶對(duì)石英砂內(nèi)部溫度進(jìn)行監(jiān)測(cè)，測(cè)量頻率為1次/min。試驗(yàn)過程中分別對(duì)進(jìn)、出水室的壓力計(jì)讀數(shù)和排水箱累積出流量進(jìn)行測(cè)定。當(dāng)砂箱內(nèi)各觀測(cè)點(diǎn)溫度均達(dá)到穩(wěn)定不再變化后，試驗(yàn)結(jié)束。當(dāng)砂箱冷卻至常溫，保持該試驗(yàn)水頭恒定，重復(fù)試驗(yàn)。

1.3 介質(zhì)水熱參數(shù)反演

1.3.1 二維飽和介質(zhì)中水、熱運(yùn)移控制方程

本文采用HYDRUS-2D軟件對(duì)飽和穩(wěn)定流條件下的介質(zhì)水、熱遷移參數(shù)進(jìn)行反演。二維飽和介質(zhì)中水流運(yùn)動(dòng)的控制方程為[16]

式中為壓力水頭，cm；x為方向笛卡爾坐標(biāo)，cm；K為各向異性張量的無量綱分量；K為石英砂的飽和導(dǎo)水率，cm/min。

二維飽和介質(zhì)中熱運(yùn)移控制方程為[17]

式中為石英砂的溫度，℃；λ()為石英砂的表觀熱導(dǎo)率，W/(m·℃)；()和C分別表示多孔介質(zhì)和水的體積比熱容，J/(m·℃)。該方程忽略水汽彌散對(duì)熱運(yùn)移的影響；q為方向達(dá)西流速，cm/min；為時(shí)間，min。

多孔介質(zhì)的體積比熱容()具有廣義可加性，可以表示為

式中為體積百分比，cm3/cm3；下標(biāo)、分別表示多孔介質(zhì)中固相和液相。由于本試驗(yàn)采用飽和狀態(tài)下的石英砂，故認(rèn)為有機(jī)質(zhì)及氣相兩部分可忽略。

根據(jù)de Marsily（1986）[18]，表觀熱導(dǎo)率λ()可以表示為

其中，多孔介質(zhì)導(dǎo)熱系數(shù)0()為

式中為達(dá)西流通量密度，cm/min；δ為Delta函數(shù)；λ、λ分別表示縱、橫向熱彌散度，cm；1、2、3均為計(jì)算導(dǎo)熱系數(shù)的經(jīng)驗(yàn)參數(shù)。

1.3.2 初始條件與邊界條件

模擬區(qū)域?yàn)?0 cm×40 cm，采用有限元方法將區(qū)域劃分為邊長(zhǎng)為1 cm的均勻三角形網(wǎng)格，共計(jì)4 000個(gè)。模擬區(qū)bc為供水邊，ef為排水邊，均定義為定水頭邊界，上下邊界的水頭差為6.5 cm。與水頭邊界相對(duì)應(yīng)的溫度邊界設(shè)置為第1類溫度邊界，采用進(jìn)、出水室內(nèi)部實(shí)測(cè)的熱電偶溫度。其他邊界設(shè)置為零通量邊界。模擬區(qū)域初始含水率為各粒徑石英砂的飽和含水率，初始溫度為常溫。

1.3.3 參數(shù)的設(shè)置及反演

不同粒徑石英砂的體積比熱容和殘余含水率的取值參照文獻(xiàn)[19]，詳細(xì)取值如表2所示。待估參數(shù)分別為介質(zhì)的飽和導(dǎo)水率（K）、導(dǎo)熱系數(shù)0的經(jīng)驗(yàn)參數(shù)（1、2、3）、縱向熱彌散度（）和橫向熱彌散度（）。分別設(shè)置3種基于不同數(shù)據(jù)類型的反演情景對(duì)水熱參數(shù)進(jìn)行估計(jì)：僅采用觀測(cè)點(diǎn)溫度（R1）、觀測(cè)點(diǎn)溫度+水流通量（R2）和觀測(cè)點(diǎn)溫度+水流通量+熱量損失（R3）。并對(duì)R1情景設(shè)置3種不同參數(shù)反演組合，即同時(shí)對(duì)2組參數(shù)（K和（1、2、3））、3組參數(shù)（K和λ、λ）和4組參數(shù)（K、（1、2、3）、λ和λ）進(jìn)行估計(jì)。利用HYDRUS-2D中的反問題模型對(duì)參數(shù)進(jìn)行估計(jì)，該模型采用Levenberg-Marquardt非線性優(yōu)化算法，參數(shù)優(yōu)化的目標(biāo)函數(shù)為[12]

式中為、某方向坐標(biāo)，cm；為實(shí)測(cè)值；v為待估參數(shù)類型；為待估參數(shù)。y為實(shí)測(cè)溫度，y為相同測(cè)點(diǎn)、相同時(shí)刻（t）下的模擬溫度，℃；m為每類待估參數(shù)的個(gè)數(shù)，n、w分別為實(shí)測(cè)溫度的數(shù)量及其權(quán)重。本文采用石英砂內(nèi)部溫度對(duì)其待估的水熱參數(shù)進(jìn)行反演，初始值如表3所示。參數(shù)優(yōu)化的目標(biāo)函數(shù)采用的求解方法為最速下降法[20]。

表2 模型中已知的水熱參數(shù)設(shè)置

表3 待估參數(shù)的初值設(shè)置

進(jìn)一步對(duì)介質(zhì)的待估參數(shù)（飽和導(dǎo)水率、導(dǎo)熱系數(shù)和縱向熱彌散度）在模型中的參數(shù)敏感性進(jìn)行分析。假設(shè)每個(gè)參數(shù)都是獨(dú)立變化的，并加以±5%的參數(shù)擾動(dòng)，并以式（7）中的作為評(píng)價(jià)指標(biāo)[21]，對(duì)HYDRUS-2D模型中的3個(gè)參數(shù)進(jìn)行敏感性分析。

=(O?O)/O（7）

式中O為參數(shù)擾動(dòng)后的值；O為敏感性分析時(shí)參數(shù)的初值。

2 結(jié)果與分析

2.1 基于觀測(cè)點(diǎn)溫度的參數(shù)反演

在相同水力梯度和熱源溫度條件下，得到不同時(shí)刻的模擬溫度空間分布隨時(shí)間的變化如圖2所示。結(jié)果表明，由于熱源邊界為部分線源邊界，當(dāng)熱量隨著水流運(yùn)動(dòng)進(jìn)入砂箱后，砂箱中部升溫較快，兩側(cè)溫度變化較慢，并且在熱運(yùn)移鋒前緣運(yùn)動(dòng)至出水邊界之前，熱運(yùn)移鋒都將保持這一溫度分布趨勢(shì)逐漸向前運(yùn)動(dòng)。當(dāng)熱運(yùn)移鋒前緣運(yùn)動(dòng)至出水邊界后，砂箱中部與兩側(cè)的溫差減小，整個(gè)砂箱溫度也逐漸升高至接近進(jìn)水室溫度。試驗(yàn)期間由于熱量損失始終存在且無法避免，距熱源越近則熱量損失越大。由于熱量損失逐漸累積，因此當(dāng)試驗(yàn)結(jié)束時(shí)砂箱溫度仍然呈現(xiàn)出沿水流運(yùn)動(dòng)方向逐漸下降的變化趨勢(shì)。3種不同粒徑的均質(zhì)石英砂具有相似的溫度變化過程。由同一時(shí)刻不同介質(zhì)中的溫度分布可知，熱運(yùn)移鋒的運(yùn)動(dòng)速率隨石英砂粒徑的增大而逐漸加快，表明輸入相同水頭條件下水流運(yùn)動(dòng)的平均速度是飽和均質(zhì)粗砂最高，均質(zhì)中砂次之，均質(zhì)細(xì)砂最低。

圖2 不同均質(zhì)石英砂的模擬溫度隨時(shí)間（t）的變化

利用砂箱中20個(gè)熱電偶的實(shí)測(cè)溫度數(shù)據(jù)對(duì)飽和均質(zhì)介質(zhì)的水熱參數(shù)進(jìn)行反演。由于同時(shí)反演介質(zhì)的水流和熱運(yùn)動(dòng)參數(shù)具有較高的精度[12]，因此在僅采用觀測(cè)點(diǎn)溫度對(duì)介質(zhì)水熱參數(shù)進(jìn)行反演（R1）的條件下設(shè)置3種反演情景：S1（對(duì)K和1、2、3進(jìn)行估計(jì)）、S2（對(duì)K、和進(jìn)行估計(jì)）和S3（對(duì)K、1、2、3、和同時(shí)估計(jì)），其中未參與估計(jì)的參數(shù)取初始值。飽和導(dǎo)水率、導(dǎo)熱系數(shù)和熱彌散度3個(gè)參數(shù)之間的相關(guān)性不強(qiáng)，Nakhaei 和?im?nek[12]的研究結(jié)果表明，當(dāng)采用HYDRUS的反問題模型對(duì)參數(shù)進(jìn)行反演時(shí)，參與反演的參數(shù)個(gè)數(shù)小于6個(gè)時(shí)，能得到可靠的結(jié)果。不同反演情景下得到的均質(zhì)中砂的水熱參數(shù)如表4所示。

表4 基于溫度數(shù)據(jù)的參數(shù)反演結(jié)果

注：S1表示介質(zhì)熱彌散度已知，估計(jì)介質(zhì)飽和導(dǎo)水率和導(dǎo)熱系數(shù)經(jīng)驗(yàn)參數(shù)；S2表示導(dǎo)熱系數(shù)經(jīng)驗(yàn)參數(shù)已知，估計(jì)介質(zhì)飽和導(dǎo)水率和熱彌散度；S3表示同時(shí)估計(jì)飽和導(dǎo)水率、熱彌散度與導(dǎo)熱系數(shù)經(jīng)驗(yàn)參數(shù)；#表示該參數(shù)未參與反演。

Note: S1 represents the thermal dispersivities are known and the empirical parameters of thermal conductivity and the saturated hydraulic conductivity are estimated; S2 represents the thermal conductivities are known and the thermal dispersivity and the saturated hydraulic conductivity are estimated; S3 represents the empirical parameters of thermal conductivity, thermal dispersivities and saturated hydraulic conductivity are estimated; # represents the parameter is not involved in estimation.

結(jié)果表明（表4），縱向彌散度的取值變化不大，基本穩(wěn)定在1.00～3.44 cm之間，橫向彌散度取值均為0.10 cm。將導(dǎo)熱系數(shù)經(jīng)參數(shù)1、2、3代入式（5）可計(jì)算得到飽和均質(zhì)石英砂的導(dǎo)熱系數(shù)（0）。S1情景下的中砂的0為1.10 W/(m·℃)，與表1中參考值（1.950 W/(m·℃)）差異較大，而S3情景下中砂0為1.93 W/(m·℃)，與參考值基本一致，表明同時(shí)對(duì)介質(zhì)飽和導(dǎo)水率、導(dǎo)熱系數(shù)和縱向、橫向彌散度進(jìn)行估計(jì)可以提高介質(zhì)導(dǎo)熱系數(shù)的估計(jì)精度。此外，在導(dǎo)熱系數(shù)估計(jì)較為合理的條件下（即S2和S3情景），介質(zhì)飽和導(dǎo)水率的估計(jì)精度較S1明顯提高。3種反演情景下的溫度時(shí)序曲線均具有較高的決定系數(shù)（2>0.85），表明模型模擬溫度與實(shí)測(cè)溫度擬合較好，且S3情景下的擬合溫度精度最高，模擬累積流量的相對(duì)誤差最小，因此選擇S3反演情景對(duì)粗砂和細(xì)砂的水熱參數(shù)進(jìn)行反演，結(jié)果見表4。石英砂飽和導(dǎo)水率的估計(jì)值隨介質(zhì)粒徑增大而增大，而介質(zhì)縱向熱彌散度則隨粒徑減小而增大，而橫向熱彌散度取值不變，均為0.10 cm。

以均質(zhì)中砂S3情景下的參數(shù)值作為敏感性分析的參數(shù)初值，以整個(gè)石英砂層吸收的累積熱量為目標(biāo)，對(duì)以上3個(gè)水熱參數(shù)進(jìn)行敏感性分析，結(jié)果如圖3所示。飽和導(dǎo)水率最為敏感，因此在模型模擬過程中要尤其關(guān)注飽和導(dǎo)水率的取值，以及其對(duì)溫度空間分布的影響。導(dǎo)熱系數(shù)和縱向彌散度的參數(shù)敏感性依次減小。

圖3 待估參數(shù)敏感性分析

在利用觀測(cè)點(diǎn)溫度對(duì)介質(zhì)參數(shù)反演的過程中，S3反演情景下參數(shù)具有比S1和S2更高的估計(jì)精度。盡管如此，模擬累積出流量與實(shí)測(cè)的累積流量間仍存在較大誤差（如表4所示），說明僅利用溫度數(shù)據(jù)對(duì)介質(zhì)水熱特性參數(shù)進(jìn)行反演會(huì)造成較大的估計(jì)誤差，尤其對(duì)參數(shù)敏感性高的參數(shù)——飽和導(dǎo)水率的估計(jì)影響顯著。因此，需要結(jié)合試驗(yàn)過程中水流運(yùn)動(dòng)和熱量傳遞的信息進(jìn)一步分析以達(dá)到對(duì)介質(zhì)水熱參數(shù)準(zhǔn)確估計(jì)的目的。

2.2 考慮水流通量和熱量損失的參數(shù)反演

在僅利用觀測(cè)溫度對(duì)介質(zhì)水熱參數(shù)反演的基礎(chǔ)上，分別增加實(shí)測(cè)水流通量和試驗(yàn)過程中的熱量損失2種不同類型的數(shù)據(jù)對(duì)介質(zhì)的水熱參數(shù)進(jìn)行反演，分別設(shè)置R2（觀測(cè)點(diǎn)溫度+水流通量）和R3（觀測(cè)點(diǎn)溫度+水流通量+熱量損失）2個(gè)情景。其中熱量損失通過計(jì)算觀測(cè)點(diǎn)虛擬溫度來實(shí)現(xiàn)，具體方法為：1）假設(shè)試驗(yàn)過程中所有損失熱量（包括有機(jī)玻璃砂箱吸收的熱量和散失到環(huán)境中的熱量）均被石英砂層吸收；2）計(jì)算無熱量損失條件下砂層應(yīng)達(dá)到的溫度，即虛擬溫度。該虛擬溫度可以認(rèn)為是觀測(cè)點(diǎn)實(shí)測(cè)溫度與熱量損失共同作用的結(jié)果，具體計(jì)算原理及方法可參考文獻(xiàn)[22]。

不同介質(zhì)中，基于水流通量和熱量損失的參數(shù)反演結(jié)果如表5所示。結(jié)果表明，不同介質(zhì)中，R2和R3情景下的飽和導(dǎo)水率估計(jì)值均大于S3情景。R2和R3情景下的介質(zhì)縱向彌散度估值變化不大，橫向彌散度的估計(jì)值則有所降低，表明垂直于水流方向的熱量遷移有所減少，更多的熱量交換發(fā)生在沿水流運(yùn)動(dòng)方向?？v橫向的熱彌散比變化范圍在16.67～81.20之間，該結(jié)果與普遍認(rèn)同的縱向彌散度與橫向彌散度的倍數(shù)比例（10~100）基本一致[21,23]。

圖4為基于不同類型數(shù)據(jù)的典型觀測(cè)點(diǎn)（2-3）的實(shí)測(cè)與模擬溫度隨時(shí)間的變化圖。砂箱內(nèi)其余各觀測(cè)點(diǎn)與典型觀測(cè)點(diǎn)呈相似的溫度變化過程。均質(zhì)砂介質(zhì)中，不同反演情景的模擬溫度在試驗(yàn)前期具有較大差異，飽和導(dǎo)水率估值越高，前期溫度升高越快；在試驗(yàn)開始20 min之后，不同反演情景的模擬溫度間無明顯差異。均質(zhì)粗砂介質(zhì)中，不同反演情景的模擬溫度始終高于實(shí)測(cè)觀測(cè)點(diǎn)溫度，也具有飽和導(dǎo)水率估值越高，觀測(cè)點(diǎn)溫度升高越快的變化趨勢(shì)。均質(zhì)細(xì)砂中，由于飽和導(dǎo)水率估計(jì)值較小，且不同反演情景間差異不大，因此模擬溫度時(shí)序曲線無明顯差異。不同反演情景的模擬溫度與實(shí)測(cè)溫度均具有較高的擬合精度，2均大于0.93。

表5 基于不同類型數(shù)據(jù)的參數(shù)反演結(jié)果

注：R2表示用于參數(shù)反演的數(shù)據(jù)類型包括觀測(cè)點(diǎn)溫度和實(shí)測(cè)水流通量；R3表示用于參數(shù)反演的數(shù)據(jù)類型包括觀測(cè)點(diǎn)溫度、實(shí)測(cè)水流通量和熱量損失。

Note: R2 represents measured temperature and water flux are used for estimation, R3 represents measured temperature, water flux and heat loss are used for estimation.

圖4 基于不同類型數(shù)據(jù)的典型觀測(cè)點(diǎn)模擬及實(shí)測(cè)溫度穿透曲線

2.3 介質(zhì)水熱參數(shù)估計(jì)誤差分析

對(duì)不同反演情景下的累積出流量的模擬結(jié)果進(jìn)行分析，與實(shí)測(cè)值間的相對(duì)誤差如圖5所示。

圖5 不同反演情景下累積出流量的估計(jì)誤差

如圖5可知，在任何情景設(shè)置中累積出流量的估計(jì)值均低于實(shí)測(cè)累積出流量，且隨著用于估計(jì)參數(shù)的數(shù)據(jù)類型不斷增多，估計(jì)誤差逐漸減小。與S3情景相比，R3情景下粗砂、中砂、細(xì)砂的累積流量相對(duì)誤差分別減少了9.74、6.65和12.53個(gè)百分點(diǎn)，對(duì)細(xì)砂飽和導(dǎo)水率的估計(jì)精度提升最顯著。盡管如此，在R3情景下3種均質(zhì)介質(zhì)的模擬累積流量仍存在10%～15%左右的估計(jì)誤差。相關(guān)的試驗(yàn)和模擬研究中也存在水流通量估值偏低的問題：Hopmans等[24]認(rèn)為熱脈沖幾何形狀對(duì)熱傳導(dǎo)的影響以及模擬過程未考慮熱彌散作用是導(dǎo)致估算水流通量偏小的主要原因；Ochsner等[25]提出估算誤差產(chǎn)生是某種原因增大了對(duì)流項(xiàng)，但未解釋產(chǎn)生這種現(xiàn)象的物理成因；而Gao等[26]則通過室內(nèi)試驗(yàn)證明，在土柱試驗(yàn)中，壁面流會(huì)導(dǎo)致水流通量增大10%～20%，因此通過設(shè)置擴(kuò)大系數(shù)（1.1～1.2）可以校正模型估算的水流通量。本試驗(yàn)采用的砂箱為有機(jī)玻璃材料，表面光滑，與石英砂顆粒接觸不完全，易出現(xiàn)較大的孔隙，從而在砂箱內(nèi)壁和石英砂層間形成水流的優(yōu)先流動(dòng)，導(dǎo)致實(shí)測(cè)累積出流量偏高。因此壁面流可能是本試驗(yàn)中造成模擬與實(shí)測(cè)水流通量存在差異的主要因素。

進(jìn)一步對(duì)飽和均質(zhì)介質(zhì)的導(dǎo)熱系數(shù)估計(jì)值進(jìn)行分析（見圖6）。

圖6 不同反演情景下導(dǎo)熱系數(shù)的估計(jì)誤差

圖6a表明隨用于估計(jì)參數(shù)的數(shù)據(jù)類型不斷增多，飽和中砂的導(dǎo)熱系數(shù)估計(jì)精度逐漸提高（S2情景下導(dǎo)熱系數(shù)未估計(jì)）；且除S1情景外，其余情景下的估計(jì)誤差均低于5%，在合理估計(jì)的范圍內(nèi)。如圖6b所示，粗砂和細(xì)砂的導(dǎo)熱系數(shù)估計(jì)值始終與參考值存在較大差異，其中S3情景下的導(dǎo)熱系數(shù)與參考值差異最小，R2情景下的導(dǎo)熱系數(shù)與參考值差異最大。S3情景下的飽和粗砂導(dǎo)熱系數(shù)估計(jì)值為3.76 W/(m·℃)，較參考值高估了54.28%；飽和細(xì)砂的導(dǎo)熱系數(shù)估計(jì)值為1.43 W/(m·℃)，較參考值低估了22.78%。介質(zhì)導(dǎo)熱系數(shù)主要受介質(zhì)自身結(jié)構(gòu)、固體顆粒的導(dǎo)熱系數(shù)、飽和度、壓實(shí)程度等多種因素影響[27]。本試驗(yàn)中，細(xì)砂的粒徑均勻統(tǒng)一，壓實(shí)程度較高，而粗砂則由于存在較大顆粒，因此壓實(shí)程度較低。不同的壓實(shí)程度可能造成飽和多孔介質(zhì)導(dǎo)熱系數(shù)的估計(jì)誤差。雖然粗砂與細(xì)砂的導(dǎo)熱系數(shù)估計(jì)值與參考值存在較大差異，但飽和石英砂導(dǎo)熱系數(shù)隨孔隙度（飽和含水率）的變化趨勢(shì)與參考值保持一致，即隨著介質(zhì)孔隙度增大，飽和石英砂的導(dǎo)熱系數(shù)逐漸降低。

3 結(jié) 論

本文在開展3種不同粒徑的飽和均質(zhì)石英砂的熱示蹤試驗(yàn)的基礎(chǔ)上，分別將不同類型的數(shù)據(jù)與HYDRUS模型相結(jié)合，對(duì)二維飽和均質(zhì)介質(zhì)的水、熱運(yùn)移參數(shù)進(jìn)行了反演與分析。得到主要結(jié)論如下：

1）溫度數(shù)據(jù)能反映均質(zhì)介質(zhì)中水熱遷移過程，將其與HYDRUS-2D模型相結(jié)合可用來反演介質(zhì)的水力與熱力學(xué)參數(shù)。同時(shí)對(duì)介質(zhì)飽和導(dǎo)水率、導(dǎo)熱系數(shù)以及熱彌散度進(jìn)行估計(jì)，有利于提高水熱參數(shù)估計(jì)精度；并且在導(dǎo)熱系數(shù)估計(jì)較為合理的條件下可以減小溫度數(shù)據(jù)對(duì)介質(zhì)飽和導(dǎo)水率的估計(jì)誤差。在恒定熱源條件下，當(dāng)砂箱吸收的熱量不再變化，溫度數(shù)據(jù)也不再體現(xiàn)不同介質(zhì)間的滲透性差異。

2）飽和導(dǎo)水率是反演參數(shù)中敏感性最高的參數(shù)，增加用于參數(shù)反演的水流運(yùn)動(dòng)和熱量傳遞信息時(shí)，粗砂、中砂、細(xì)砂的累積流量相對(duì)誤差分別減少了9.74、6.65和12.53個(gè)百分點(diǎn)，顯著提高了介質(zhì)飽和導(dǎo)水率的反演精度。反演所得的飽和導(dǎo)水率估值隨粒徑的增大而顯著增大。在觀測(cè)點(diǎn)溫度擬合較好的情況下仍存在10%～15%左右的水流通量估計(jì)誤差，在考慮熱量損失的條件下，壁面流可能是造成累積流量誤差的主要原因。

3）導(dǎo)熱系數(shù)和熱彌散度在水熱遷移模型中的敏感性均很低，縱向熱彌散度隨粒徑的減小而增大，橫向熱彌散度則基本無變化。增加用于參數(shù)反演的水流運(yùn)動(dòng)和熱量傳遞信息能顯著提高中砂的導(dǎo)熱系數(shù)反演精度。

[1]吳志偉，宋漢周. 地下水溫度示蹤理論與方法研究進(jìn)展[J]. 水科學(xué)進(jìn)展，2011，22(5)：733-739.

Wu Zhiwei, Song Hanzhou. Temperature as a groundwater tracer: Advances in theory and methodology[J]. Advances in Water Science, 2011, 22(5): 733-739. (in Chinese with English abstract)

[2]秦耀軍，周曉勇，楊亞賓，等. 基于數(shù)據(jù)同化技術(shù)的地質(zhì)參數(shù)反演分析研究[J]. 水科學(xué)與工程技術(shù)，2017(6)：78-82.

Qin Yaojun, Zhou Xiaoyong, Yang Yabin, et al. Back analysis of geological parameters based on data assimilation[J]. Water Sciences and Engineering Technology, 2017(6): 78-82. (in Chinese with English abstract)

[3]閔濤，周孝德，馮民權(quán). 非線性布西尼斯克方程的直線解法及滲透系數(shù)反演計(jì)算[J]. 水利學(xué)報(bào)，2004，35(7)：21-25.

Min Tao, Zhou Xiaode, Feng Minquan. Numerical solution of Boussinesq equation based on straight-line method and inverse calculation of permeability coefficient[J]. Journal of Hydraulic Engineering, 2004, 35(7): 21-25. (in Chinese with English abstract)

[4]李守巨，上官子昌，劉迎曦，等. 地下水滲流模型參數(shù)識(shí)別的模擬退火算法[J]. 巖石力學(xué)與工程學(xué)報(bào)，2005，24(A01)：5031-5036.

Li Shouju, Shangguan Zichang, Liu Yingxi, et al. Parameter identification procedure for groundwater flow model with simulated annealing[J]. Chinese Journal of Rock Mechanics and Engineering, 2005, 24(A01): 5031-5036. (in Chinese with English abstract)

[5]劉曉楠，張樸，孔力，等. 基于SVD和DE反演方法的滲透率估計(jì)[J]. 華中科技大學(xué)學(xué)報(bào)：自然科學(xué)版，2014，42(9)：5-9，18.

Liu Xiaonan, Zhang Pu, Kong Li, et al. Permeability estimation based on singular value decomposition and differential evolution inversion method[J]. Journal of Huazhong University of Science and Technology: Natural Science Edition, 2014, 42(9): 5-9,18. (in Chinese with English abstract)

[6]王兵賢，王澤文，徐定華，等. 二維流Boussinesq方程滲透系數(shù)反演的變分伴隨方法[J]. 水利水電科技進(jìn)展，2010，30(6)：11-14.

Wang Binxian, Wang Zewen, Xu Dinghua, et al. Variational adjoint method for inversion of permeability coefficient of two-dimensional flow Boussinesq equation[J]. Advances in Science and Technology of Water Resources, 2010, 30(6): 11-14. (in Chinese with English abstract)

[7]任圖生，邵明安，巨兆強(qiáng). 利用熱脈沖時(shí)域反射技術(shù)測(cè)定土壤水熱動(dòng)態(tài)和物理參數(shù)Ⅱ. 應(yīng)用[J]. 土壤學(xué)報(bào)，2004，41(4)：523-529.

Ren Tusheng, Shao Ming’an, Ju Zhaoqiang. Measurement of soil physical properties with thermo-time domain reflectometry Ⅱ. application[J]. Acta Pedologica Sinica, 2004, 41(4): 523-529. (in Chinese with English abstract)

[8]任圖生，邵明安，巨兆強(qiáng)，等. 利用熱脈沖-時(shí)域反射技術(shù)測(cè)定土壤水熱動(dòng)態(tài)和物理參數(shù)Ⅰ. 原理[J]. 土壤學(xué)報(bào)，2004，41(2)：225-229.

Ren Tusheng, Shao Ming’an, Ju Zhaoqiang, et al. Measurement of soil physical properties with thermo-time domain reflectometry Ⅰ.theory[J]. Acta Pedologica Sinica, 2004, 41(2): 225-229. (in Chinese with English abstract)

[9]Anderson M P. Heat as a ground water tracer[J]. Ground Water, 2005, 43(6): 951-968.

[10]Constantz J. Heat as a tracer to determine streambed water exchanges[J]. Water Resource Research, 2008, 44: W00D10.

[11]Steenpass C, Vanderborght J, Herbst M, et al. Estimating soil hydraulic properties from infrared measurement of soil surface temperature and TDR data[J]. Vadose Zone Journal, 2010, 9: 910-924.

[12]Nakhaei M, ?im?nek J. Parameter estimation of soil hydraulic and thermal property functions for unsaturated porous media using the HYDRUS-2D code[J]. Journal of Hydrology, 2014, 62(1): 7-15.

[13]Giambastiani B M S N C, Mastrocicco M. Limitation of using heat as a groundwater tracer to define aquifer properties: experiment in a large tank model[J]. Environmental Earth Science, 2013, 70: 719-728.

[14]Saley A D, Jardani A, Ahmed A S, et al. Hamiltonian Monte Carlo algorithm for the characterization of hydraulic conductivity from the heat tracing data[J]. Advance in Water Resource, 2016, 97: 120-129.

[15]Somogyvari M, Bayer P. Field validation of thermal tracer tomography for reconstruction of aquifer heterogeneity[J]. Water Resources Research 2017, 53(6): 5070-5084.

[16]Bear J. Dynamics of Fluids in Porous Media[M]. New York: Academic Press, 1972.

[17]Sophocleous M. Analysis of water and heat flow in unsaturated-saturated porous media[J]. Water Resources Research, 1979, 15(5): 1195-1206.

[18]de Marisly G. Quantitative Hydrogeology[M]. London: Academic Press, 1986.

[19]潘夢(mèng)綺，黃權(quán)中，馮榕，等. 基于熱示蹤的飽和層狀介質(zhì)水熱參數(shù)反演[J]. 水利學(xué)報(bào)，2017，48(11)：1363-1372.

Pan Mengqi, Huang Quanzhong, Feng Rong, et al. Estimation of hydraulic and thermal parameters in saturated layered porous media based on heat tracing method[J]. Journal of Hydraulic Engineering, 2017, 48(11): 1363-1372. (in Chinese with English abstract)

[20]Donald W M. An algorithm for least-squares estimation of nonlinear parameters[J]. Journal of Society for Industrial and Applied Mathematics, 1963, 11(2): 431-441.

[21]Rau G C, Anderson M S, McCallum A M, et al. Heat as a tracer to quantify water flow in near-surface sediments[J]. Earth Science Reviews, 2014, 129: 40-58.

[22]Pan M Q, Huang Q Z, Feng R, et al. Improving the estimation of hydraulic and thermal properties of heterogeneous media via the addition of heat loss[J]. Vadose Zone Journal, 2019, 18(1): 180149.

[23]Constanz J, Cox M H, Su G W. Comparison of heat and bromide as ground water tracers near stream[J]. Ground Water, 2003, 41(5): 647-656.

[24]Hopmans J W, Simunek J, Bristow K L. Indirect estimation of soil thermal properties and water flux using heat pulse probe measurements: Geometry and dispersion effects[J]. Water Resources Research, 2002, 38(1): 1006-1019.

[25]Ochsner T E，Horton R，Kluitenberg G J, et al. Evaluation of the heat pulse ratio method for measuring soil water flux[J]. Soil Science Society of America Journal, 2005, 69(3): 757-765.

[26]Gao J Y, Ren T S, Gong Y S. Correcting wall flow effect improves the heat-pulse technique for determining water flux in saturated soils[J]. Soil Science Society of America Journal, 2006, 70(3): 711-717.

[27]Ochsner T E, Horton R, Ren T S. A new perspective on soil thermal properties[J]. Soil Science Society of America Journal, 2001, 65(6): 1641-1647.

Estimation of water and heat transfer parameters of saturated silica sand by using different types of data

Pan Mengqi, Huang Quanzhong※, Feng Rong, Huang Guanhua

(1.,,100083,; 2.,,100083,)

The hydraulic and thermal parameters of soil and aquifer are very important for quantitative description of soil water, groundwater migration and its accompanying heat and salt transport. In order to investigate the influence of different types of data on the hydraulic and thermal parameter estimation, three heat tracing experiments packed with saturated homogeneous silicon sand of three different particle size (coarse, medium and fine sands) were conducted under steady state. Twenty T-thermocouple probes were installed uniformly in the sandbox to record the temperature of the silicon sand. Two additional T-thermocouple probes were installed in the inflow and outflow chambers to record the water temperature. And then, the measured temperatures of silicon sand and water were applied for the inverses model of HYDRUS-2D software to estimate the saturated hydraulic conductivity, the thermal conductivity and the longitudinal and transverse thermal dispersivity of three saturated silicon sands. In this study, three scenarios based on different types of data were designed to estimate these parameters, i.e. R1 (including the measured temperature at the observation point alone), R2 (including the measured temperature at the observation point and the cumulated outflow) and R3 (including the measured temperature at the observation point, the cumulated outflow and the heat loss). In the addition, three more scenarios under scenario R1 consisted of different numbers of parameters were set, i.e. S1 (the thermal dispersivities were known and the empirical parameters of thermal conductivity and the saturated hydraulic conductivity were estimated); S2 (the thermal conductivities were known and the thermal dispersivity and the saturated hydraulic conductivity were estimated); S3 (the empirical parameters of thermal conductivity, thermal dispersivities and saturated hydraulic conductivity were estimated). The results showed that the thermal conductivity, the longitudinal and transverse thermal dispersivities and the saturated hydraulic conductivity of silicon sand estimated at the same time could significantly improve the accuracy of parameter estimation, and the accuracy of the saturated hydraulic conductivity was improved when the thermal conductivity was reasonable under scenario R1. The saturated hydraulic conductivity was the parameter with the highest sensitivity in parameter estimation, followed by the thermal conductivity and the longitudinal thermal dispersivity. When the simulated temperatures were consistent with the measured temperatures at observation points, there was still 10%-15% estimated error of cumulated outflow. When considering heat loss, wall flow may be the main reason for the estimated error of cumulated outflow. The additional information of water flow and heat loss was helpful to reduce the estimated error of saturated hydraulic conductivity, and then the relative error of cumulated outflow was significantly decreased for coarse, medium and fine sands. The estimated value of saturated hydraulic conductivity increased with the increasing of particle size of silicon sands while the longitudinal thermal dispersivity showed the opposite trend. And the value of transverse thermal dispersivity was same for all the sands. The additional information of water flow and heat loss improved the estimation of thermal conductivity of medium sand as well. The estimated value of thermal conductivitydecreased with increasing particle size of silicon sands. This study can help for parameter estimation of homogeneous porous media based on different type of data.

data; parameters; heat tracing; saturated silica sand; water flow and heat transport models; HYDRUS-2D software

潘夢(mèng)綺，黃權(quán)中，馮榕，等. 采用不同監(jiān)測(cè)數(shù)據(jù)組合反演飽和均質(zhì)石英砂水熱運(yùn)移參數(shù)[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2020，36(10)：75-82.doi：10.11975/j.issn.1002-6819.2020.10.009 http://www.tcsae.org

Pan Mengqi, Huang Quanzhong, Feng Rong, et al. Estimation of water and heat transfer parameters of saturated silica sand by using different types of data[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2020, 36(10): 75-82. (in Chinese with English abstract) doi：10.11975/j.issn.1002-6819.2020.10.009 http://www.tcsae.org

2020-01-10

2020-05-10

國(guó)家自然科學(xué)基金項(xiàng)目（51779256、51639009）；國(guó)家重點(diǎn)研發(fā)計(jì)劃項(xiàng)目（2016YFC0501304）

潘夢(mèng)綺，博士生，主要從事多孔介質(zhì)中水熱遷移規(guī)律的研究。Email：panmengqi@cau.edu.cn

黃權(quán)中，博士，教授，主要從事多孔介質(zhì)水-熱-鹽運(yùn)移機(jī)制及定量表征的研究。Email：huangqzh@cau.edu.cn

10.11975/j.issn.1002-6819.2020.10.009

P641

A

1002-6819(2020)-10-0075-08