由溫度時(shí)序資料反演地下水流速的兩種解析解及其比較

吳志偉，宋漢周

1.河海大學(xué)地球科學(xué)與工程學(xué)院，南京 210098

2.東北電力設(shè)計(jì)院，長春 130033

0 引言

地表水與地下水交換是水文學(xué)、水文地質(zhì)學(xué)、環(huán)境科學(xué)、水資源管理等學(xué)科的研究重點(diǎn)之一，研究河床地表水與地下水交換的時(shí)空變化特征具有重要意義［1］。研究［2］發(fā)現(xiàn)，水作為熱量的優(yōu)良載體，在地表水與地下水交換過程中會(huì)攜帶能量變化的信息。利用熱量在地質(zhì)體中的傳遞規(guī)律，觀測河床地下水溫度場的時(shí)空分布，能夠用來確定地表水與地下水之間的交換特征［3］。

早在20世紀(jì)60年代，Stallman［4］就介紹了利用地下水溫度反演飽和多孔介質(zhì)中垂向流流向和流速的方法——假設(shè)上部溫度恒定，通過解析演算可以得到采用垂向溫度曲線形態(tài)計(jì)算地下水流速的典型曲線法。此后，根據(jù)鉆孔溫度梯度的變化，采用Stallman模型計(jì)算地下水流速得到了廣泛的應(yīng)用［5］。但是Stallman模型假設(shè)淺部溫度恒定，這顯然是過于簡化了，后續(xù)的研究開展了采用連續(xù)測溫資料的“突變點(diǎn)”確定水文地質(zhì)條件突變的方法［6］。而對(duì)于連續(xù)的溫度波動(dòng)資料，目前多數(shù)還是采用溫度曲線形態(tài)定性地解釋地下水流動(dòng)，確定地表水與地下水交換的時(shí)空分布［7－9］，但是相關(guān)的定量研究還很少。滲流場與溫度場耦合的數(shù)值模擬能有效地解決復(fù)雜的滲流傳熱問題［10］。數(shù)值模擬技術(shù)也在解釋測溫資料中得到應(yīng)用，通過“試錯(cuò)法”擬合實(shí)測的鉆孔溫度曲線，得到地下水垂向流速［11］。

研究［2］發(fā)現(xiàn)，淺部地質(zhì)體中的溫度受大氣溫度影響，具有周期性波動(dòng)特征，而不同深度處的溫度波動(dòng)是衰減的，且這種衰減特征與地下水活動(dòng)密切相關(guān)；由此建立了采用溫度時(shí)序資料計(jì)算垂向地下水流速的解析法。典型的解析模型有Hatch解［12］和Keery解［13］。 實(shí) 例 研 究［14－16］證 實(shí) 了 采 用 溫 度 時(shí) 序曲線解析法能方便快捷地計(jì)算地下水垂向流速、確定地表水與地下水交換的時(shí)空分布特征。

上述解析模型假設(shè)垂向上能觀測到溫度的衰減，但若上、下2個(gè)測點(diǎn)的溫度波動(dòng)曲線一致，那么解析模型就不再適用；此外，在建立基本模型時(shí)，Keery解忽略了熱彌散效應(yīng)的影響，這必然引入一定的誤差，已有研究尚未充分討論Hatch解與Keery解的差異。

鑒于此，筆者基于解析模型的基本假設(shè)，通過算例研究了解析模型的適用性和局限性，討論了Hatch解與Keery解的異同點(diǎn)；以期為有效地獲取溫度實(shí)測資料，并為由此解析地表水與地下水之間的交換而選用合適的求解方法提供依據(jù)。

1 溫度資料反演地下水流速解析模型

1.1 基本原理

地表水體及其底部河床中滲透性介質(zhì)的溫度受氣溫影響，溫度時(shí)序曲線具有日波動(dòng)、季節(jié)波動(dòng)等特征。河水與河床滲透性介質(zhì)之間的熱量交換受熱對(duì)流作用的控制，即與地表水與地下水的交換過程密切相關(guān)。當(dāng)?shù)乇硭c地下水發(fā)生交換時(shí)，水流大小及方向?qū)Χ咧g的能量傳遞過程具有顯著影響。

如圖1所示，在某測點(diǎn)垂向上可在不同位置（單井不同深度埋設(shè)探頭）連續(xù)測量，以獲得垂向上的溫度波動(dòng)差異。為了達(dá)到不同的研究目的，溫度記錄可以設(shè)置以時(shí)、天、月或年為時(shí)間間隔，得到觀測點(diǎn)連續(xù)變化的溫度波動(dòng)資料。即，淺部測點(diǎn)（埋深為zs）溫度時(shí)序曲線的振幅（As）必然大于深部測點(diǎn)（埋深為zd）的溫度振幅（Ad），與此同時(shí)，二者存在相位滯后（Δφ）。

圖1 溫度時(shí)序資料確定地下水流速（v）方法示意圖Fig.1 Sketch for determining seepage velocity（v）of streambed by temperature time－series records

如果沒有地下水活動(dòng)，河床淺部能觀測到溫度的日波動(dòng)，再深的部位能觀測到溫度的季節(jié)波動(dòng)或年波動(dòng)。而在垂直向下水流作用的部位，表層水溫的波動(dòng)能夠在巖土體中發(fā)生傳導(dǎo)和對(duì)流，且滲流速度越大，熱對(duì)流越強(qiáng)烈，在較深部位也能觀測到溫度波動(dòng)，而且溫度振幅衰減較慢，峰值延時(shí)較小；垂直向上的水流，對(duì)比于向下的水流作用，相同埋深點(diǎn)的溫度振幅較小，峰值延時(shí)較大［17］。采用1D對(duì)流－傳熱模型，假設(shè)河水溫度按正弦波動(dòng)，得到河床滲透性介質(zhì)的溫度波動(dòng)解析解；進(jìn)一步得到利用同一測點(diǎn)不同深度處的溫度時(shí)序曲線振幅比Ar（Ar＝As/Ad）和相位滯后Δφ來計(jì)算垂向地下水流速的解析法。

1.2 解析模型及2種解析解

1.2.1 1D瞬態(tài)對(duì)流－傳熱模型

在自然界中，河床地下水的垂向流動(dòng)是普遍存在的，尤其是在河流的河谷中部。因此可以采用1D對(duì)流－傳熱模型來描述河床溫度場。

首先對(duì)河床傳熱模型做如下簡化：假設(shè)河床頂面水平、下部空間無限大、河床為均質(zhì)各向同性的多孔介質(zhì)，那么該半無限空間的1D瞬態(tài)熱對(duì)流－傳熱方程［4］為

式中：T為深度z（m）處時(shí)刻t（s）的溫度（℃）；ke為飽和多孔介質(zhì)的等效熱擴(kuò)散系數(shù)（m2/s）；v為計(jì)算截面上的滲流流速（m/s，與z軸正方向相同為正，如圖1所示，即向下的流速為正），基于典型單元體假設(shè)，v與地下水實(shí)際流速u存在關(guān)系式v＝neu，ne為有效孔隙度；ρw和cw分別為水的密度（kg/m3）和質(zhì)量熱容（J/（kg·℃））；ρ和c分別為飽和多孔介質(zhì)的等效密度和等效質(zhì)量熱容。根據(jù)體積平均原理，ρ＝（1－ne）ρs＋neρw，c＝（1－ne）cs＋necw，ρs和cs分別為固體介質(zhì)的密度和質(zhì)量熱容。

飽和多孔介質(zhì)的等效熱擴(kuò)散系數(shù)ke定義為

式中：λe為多孔介質(zhì)的等效熱傳導(dǎo)系數(shù)（W/（m·℃））；λ0為基準(zhǔn)熱傳導(dǎo)系數(shù)（W/（m·℃）），即沒有地下水活動(dòng)的熱傳導(dǎo)系數(shù)；β為熱彌散系數(shù)（m）。對(duì)等效熱擴(kuò)散系數(shù)的定義考慮了水動(dòng)力條件對(duì)介質(zhì)傳熱的影響。

河水溫度與氣溫的相關(guān)性密切，可以簡化為按正弦波動(dòng)的簡單諧波信號(hào)：

式中：A為河水溫度波動(dòng)的振幅；P為溫度波動(dòng)周期；T0為上邊界溫度的非波動(dòng)因子，約等于平均溫度。在無窮遠(yuǎn)處的溫度邊界假定為地溫不受波動(dòng)因素及流體流動(dòng)的影響。Hatch和Keery分別給出了上述模型的解析解。

1.2.2 1D瞬態(tài)對(duì)流－傳熱模型的解析解

1）Hatch解

Hatch［12］指出式（1）的解為

其中，

上述各式中的vT（溫度前端運(yùn)移速度）與地下水實(shí)際流速u之間的轉(zhuǎn)換關(guān)系滿足下式：

因河床滲透性介質(zhì)滿足均質(zhì)各向同性假設(shè)，不同深度處的熱物理參數(shù)相同，那么深處溫度波動(dòng)曲線與淺部溫度波動(dòng)曲線的振幅比Ar和相位滯后Δφ可以由式（4）得到：

式中，Δz為深處測點(diǎn)與淺部測點(diǎn)間的距離（m），Δz＝zd－zs。

由式（7）和式（8）分別得到溫度前端運(yùn)移速度vT，進(jìn)而計(jì)算得到各計(jì)算時(shí)刻的滲流流速v。

2）Keery解

Keery［13］忽略地下水流動(dòng)的熱彌散效應(yīng)，認(rèn)為ke＝λ0/ρc，且直接采用滲流流速v計(jì)算，給出式（1）的解為

其中：

式（9）右側(cè)即為某深度處溫度曲線的波動(dòng)因子，它受飽和多孔介質(zhì)及流體的熱物理參數(shù)、流速、表層溫度波動(dòng)等因素的影響?？傮w而言，受地質(zhì)體的阻尼作用，當(dāng)?shù)叵滤魉俸愣〞r(shí)，不同深度處溫度曲線的振幅隨著深度增加而衰減，峰值相位存在一定的滯后。同樣由式（9）給出2個(gè)測點(diǎn)間的振幅比Ar和相位滯后Δφ，由此可得到由振幅比和相位滯后表示的地下水流速：

式中：H＝ρwcw/λ0。

1.3 工作方法

1.3.1 溫度測量

目前常用的溫度測量方法是將溫度傳感器布置在金屬或PVC材質(zhì)的護(hù)筒中，護(hù)筒管徑不宜過大。若護(hù)筒為金屬材質(zhì)，可直接錘入河床；若護(hù)筒為PVC材質(zhì)，可采用其他方法先成孔，然后埋設(shè)PVC管。護(hù)筒應(yīng)保持豎直狀態(tài)，且采取相應(yīng)措施防止水流沖刷破壞。護(hù)筒與被測試位置的接觸斷面應(yīng)設(shè)置濾孔和濾網(wǎng)，保證地下水能在護(hù)筒與巖土體中自由交換。測試護(hù)筒中地下水與測試點(diǎn)外圍滲透性介質(zhì)之間的熱量平衡過程受護(hù)筒導(dǎo)熱情況和濾網(wǎng)的滲透性影響。這種短暫的熱平衡過程對(duì)溫度時(shí)序資料測試成果有一定的影響，研究表明，護(hù)筒干擾對(duì)由振幅比計(jì)算地下水流速不會(huì)產(chǎn)生影響，而會(huì)給由相位滯后計(jì)算地下水流速帶來誤差［18］。

由于日周期波動(dòng)信息觀測起來較為經(jīng)濟(jì)，且波動(dòng)振幅較大，是理想的溫度時(shí)序資料，因此本文所述的溫度時(shí)序資料解析模型以日周期信號(hào)為研究對(duì)象。一般情況下，淺部河床能觀測到溫度的日波動(dòng)，傳感器埋深越大，溫度日波動(dòng)特征越弱。傳感器的埋置深度及埋設(shè)間距與地表水與地下水交換速率的大小有關(guān)：若地下水流速較大，則傳感器需埋在較深的部位，傳感器間距也應(yīng)較大；若地下水流速較小，則傳感器需埋在較淺的部位，傳感器間距也應(yīng)較小。

采用上述方法可以測得某一測點(diǎn)不同深度處的溫度時(shí)序資料，采用前述解析模型能夠計(jì)算該點(diǎn)的地下水流速在時(shí)間域上的變化。在平面上，可以布置多個(gè)測點(diǎn)，根據(jù)測點(diǎn)的勢能差，采用Darcy定律，能夠評(píng)價(jià)河床滲透性的空間差異。

1.3.2 溫度時(shí)序資料濾波處理

實(shí)測的溫度時(shí)序資料是多周期疊加的信號(hào)，為了應(yīng)用前述解析模型，必須通過對(duì)測試數(shù)據(jù)進(jìn)行濾波處理，得到按日周期波動(dòng)的正弦信號(hào)。最簡單的方法是通過快速傅里葉變換（FFT）的帶通濾波，提取日周期信號(hào)。此外，可以采用更加適合處理周期性信號(hào)的動(dòng)態(tài)諧波回歸分析法（DHR法）對(duì)實(shí)測資料進(jìn)行濾波處理。

1.3.3 地下水流速計(jì)算

深淺2個(gè)測點(diǎn)的溫度實(shí)測資料經(jīng)過濾波，得到按日波動(dòng)的溫度時(shí)序曲線，進(jìn)而讀取溫度時(shí)序曲線的振幅和峰值相位。理論上，深部測點(diǎn)的溫度振幅Ad總是小于等于淺部溫度振幅As，0.0≤Ar≤1.0；相位滯后Δφ也是大于0的數(shù)。當(dāng)巖土體熱物理參數(shù)、任意2點(diǎn)的垂向距離、溫度波動(dòng)振幅比或相位差等參數(shù)已知時(shí)，采用 Hatch解（式（7）和式（8））或Keery解（式（10）及式（11））計(jì)算地下水流速。上述過程均可采用Matlab軟件編寫相應(yīng)的程序來完成。式（7）及式（10）的根包含虛根，在實(shí)際應(yīng)用過程中，取其實(shí)數(shù)根為地下水流速。

2 Hatch解與Keery解對(duì)比分析

2.1 計(jì)算參數(shù)

Hatch解中的等效熱擴(kuò)散系數(shù)ke是流速的函數(shù)，需采用迭代法求解；而Keery解忽略了熱彌散效應(yīng)的影響，計(jì)算比較簡單。為了評(píng)價(jià)Keery解的這種處理是否有效，通過一個(gè)算例來分析該解析解由忽略熱彌散效應(yīng)所引入的誤差。

參照溶質(zhì)運(yùn)移研究中的彌散度經(jīng)驗(yàn)值，通常情況下，小尺度松散多孔介質(zhì)的熱彌散度值約為0.001～0.100m。因中砂最能代表河床滲透性介質(zhì)，取中砂作為典型介質(zhì)來評(píng)價(jià)上述兩種解析解之間可能存在的差異，相關(guān)的熱物理參數(shù)如表1所示。計(jì)算參數(shù)中的孔隙度0.37是中砂的典型值，介質(zhì)的密度、質(zhì)量熱容和熱導(dǎo)系數(shù)參照文獻(xiàn)［19］，水的熱物理參數(shù)參照相關(guān)手冊［20］給出常溫條件下的經(jīng)驗(yàn)值。

表1 模型對(duì)比分析計(jì)算參數(shù)Table 1 Parameters for calculation

取測點(diǎn)間距為0.3m，相位差為π/10～π/5，振幅比為0.1～0.6，分別采用Keery解與Hatch解法計(jì)算地下水流速vk和vh。

2.2 對(duì)比分析結(jié)果

對(duì)比相同條件下兩種解法的計(jì)算結(jié)果，得到Keery解不考慮熱彌散系數(shù)（β）引入的誤差如圖2所示。因Hatch解考慮了熱彌散效應(yīng)，認(rèn)為其計(jì)算值是相對(duì)精確的，圖中Keery解的誤差為（vk－vh）/vh。

當(dāng)采用相位滯后法（圖2a）計(jì)算時(shí)，對(duì)于1D瞬態(tài)熱－流模型（式（1））的上述兩種解法均不能確定地下水流向，得到的流速均為正值。圖2a中：誤差均為正值，說明Keery解析解計(jì)算得到的流速均大于Hatch解析解；當(dāng)β≤0.01m時(shí)，兩種解的計(jì)算結(jié)果很接近，β越大，Keery解帶來的誤差越大，當(dāng)β＝0.10m時(shí)，兩種解的計(jì)算結(jié)果最大偏差為12%；Δφ越小，兩種解的計(jì)算結(jié)果差異越小。在本算例的假設(shè)條件下，解析解不適用于Δφ大于π/5的情況，因此對(duì)比分析未給出對(duì)應(yīng)的情況。當(dāng)β取值小于0.10 m時(shí)，采用Keery解所帶來的誤差是可以接受的。

圖2 Keery解不考慮熱彌散系數(shù)β引入的誤差Fig.2 Error introduced by ignoring thermal diffusion coefficientβin Keery solution

當(dāng)采用振幅比法（圖2b）計(jì)算時(shí)，模型能夠確定地下水流向，向下的流速為正，向上的流速為負(fù)。在本算例的假設(shè)條件下，計(jì)算得到的流速均為正值。圖2b中：誤差為負(fù)表示Keery解計(jì)算的流速小于Hatch解，說明當(dāng)Ar＞0.3時(shí)，Keery解計(jì)算的流速偏?。划?dāng)Ar＝0.3時(shí)，Keery解的相對(duì)誤差最小。當(dāng)β≤0.01m時(shí)，兩種解析解的計(jì)算結(jié)果很接近，誤差小于±6%；β越大，Keery解帶來的誤差越大；當(dāng)β取值小于0.10m時(shí)，采用Keery解所帶來的誤差小于40%。當(dāng)Keery解計(jì)算值大于Hatch解時(shí)：Ar越小，Kerry解析解所產(chǎn)生的誤差越大；反之，Ar越大，Kerry解析解所產(chǎn)生的誤差越小。

由以上分析可知：β越大，Keery解所產(chǎn)生的誤差越大。當(dāng)β＜0.01m時(shí)，Keery解忽略彌散效應(yīng)的處理是可行的；當(dāng)β＞0.10m，Keery解所產(chǎn)生的誤差非常大。介質(zhì)的熱彌散度受介質(zhì)空隙特征影響，并具有一定的尺度效應(yīng)，可以結(jié)合流場介質(zhì)特征確定是否可以忽略熱彌散效應(yīng)。Keery解與Hatch解除了本節(jié)所論述的差異外，二者的理論基礎(chǔ)和計(jì)算方法是相同的。在后面的討論中，只采用精度更高的Hatch解來進(jìn)行分析。

3 模型解析解的適用性與局限性

3.1 解析解的適用性

采用溫度時(shí)序資料解析模型計(jì)算地下水流速時(shí)，2個(gè)測點(diǎn)間必須存在振幅衰減和相位滯后。為了評(píng)價(jià)模型解析解的適用性，采用表1所列的典型介質(zhì)的熱物理參數(shù)計(jì)算振幅比及相位滯后與計(jì)算流速的關(guān)系，結(jié)果如圖3所示。

圖3 模型解析解確定地下水流速的適用性Fig.3 Limitations and applicability of the analytic model to determine seepage velocity in streambed

對(duì)于振幅比法：1）采用常規(guī)的測點(diǎn)距離（1.0m左右），振幅比法能有效評(píng)價(jià)流速在（－1.0～8.0）×10－5m/s的地下水活動(dòng)。當(dāng)?shù)叵滤魉龠^小時(shí)，深處測點(diǎn)可能不存在溫度的日波動(dòng)；當(dāng)?shù)叵滤魉龠^大時(shí)，深淺2個(gè)測點(diǎn)的溫度波動(dòng)一致，同樣無法采用該解析模型。2）振幅比法可以得到地下水流向，負(fù)值表示地下水補(bǔ)給地表水，只有2個(gè)測點(diǎn)的距離很小且振幅比較小時(shí)，才可以通過溫度時(shí)序資料確定向上的地下水流動(dòng)。3）測點(diǎn)距離越大，相同的振幅比所反映的地下水流速越大，即地下水流速較大時(shí)，需采用較大的測點(diǎn)距離。4）相同的測點(diǎn)距離，振幅比越大，深部溫度波動(dòng)振幅衰減越小，則地下水流速越大。5）當(dāng)振幅比過大或過小時(shí)，振幅比的小幅變動(dòng)代表了較大范圍的流速，模型精度有所降低，因此需控制合理的測點(diǎn)間距，使得到的振幅比在合理范圍內(nèi)。

對(duì)于相位滯后法：1）采用常規(guī)的測點(diǎn)距離（1.0 m左右），相位滯后法能有效評(píng)價(jià)流速在（0.0～1.6）×10－4m/s的地下水活動(dòng)，能計(jì)算的地下水流速的范圍大于振幅比法。2）相位滯后法只能得到地下水流速，無法確定地下水流向。3）測點(diǎn)距離越大，相同的相位滯后所反映的地下水流速越大，對(duì)于大的地下水流速需采用大的測點(diǎn)距離。4）相同的測點(diǎn)距離，相位滯后越小，說明地下水流速越大。

3.2 解析解的優(yōu)勢

解析模型可以采用連續(xù)的溫度觀測資料評(píng)價(jià)研究時(shí)段內(nèi)地表水與地下水交換速率的變化，操作較為簡單。該模型特別適用于干旱地區(qū)或地下水集中開采地區(qū)的水資源研究，這些地區(qū)地表水與地下水的交換過程是動(dòng)態(tài)變化的，通過地下水流速可以很容易地得到地表水與地下水水量交換的動(dòng)態(tài)變化。

如果需要評(píng)價(jià)河床滲透性的空間分布，只需要測量測點(diǎn)水頭差，根據(jù)Darcy定律即可計(jì)算測點(diǎn)間的平均滲透系數(shù)。這種評(píng)價(jià)是基于原位測試資料得到的，避免了取樣擾動(dòng)，顯然評(píng)價(jià)結(jié)果更加可靠。

溫度時(shí)序資料解析解的計(jì)算只需要2個(gè)測點(diǎn)的溫度資料，溫度測量既經(jīng)濟(jì)又可靠，使在研究區(qū)平面上布置多個(gè)測點(diǎn)成為可能。通過多點(diǎn)的測試成果，可以評(píng)價(jià)研究區(qū)河床滲透性介質(zhì)滲透性的空間差異。

3.3 解析解的局限性

采用溫度時(shí)序資料確定地下水與地表水交換的解析模型是基于1D半無限均質(zhì)各向同性介質(zhì)中熱傳導(dǎo)－對(duì)流模型得到的，只有當(dāng)河床固相介質(zhì)可以簡化為均質(zhì)、各向同性介質(zhì)時(shí)才可以采用該模型。地質(zhì)體的非均質(zhì)性及各向異性對(duì)解析模型的精度會(huì)產(chǎn)生影響，如果淺部測點(diǎn)與深部測點(diǎn)的基本物理性質(zhì)不一致，就無法得到采用振幅比或相位滯后計(jì)算地下水流速的解析解。

解析解只能用來計(jì)算地下水垂向流速，而河床地下水的多維流動(dòng)是普遍存在的，多維流動(dòng)對(duì)河水與河床的熱量交換會(huì)產(chǎn)生一定的影響，也會(huì)給解析解帶來一定的誤差。因此在采用本文所述解析方法之前，需要初步查明河床滲透性介質(zhì)的物理性質(zhì)和地下水流場的基本特征，確定模型是否適用。

溫度時(shí)序資料日波動(dòng)信息的提取是決定模型精度的重要因素。采用合理的濾波方法提取有效的正弦信號(hào)，是處理溫度實(shí)測資料的關(guān)鍵。溫度時(shí)序資料曲線上常常含有突變點(diǎn)，這些突變點(diǎn)通常是由水動(dòng)力條件或河水溫度突變導(dǎo)致的。在濾波過程中，這些突變點(diǎn)可能會(huì)導(dǎo)致模型計(jì)算結(jié)果失真。

理論上采用豎向垂直剖面上任意2點(diǎn)的測溫資料都可以計(jì)算2點(diǎn)之間的平均地下水流速。但是測點(diǎn)間距過大而地下水流速較小時(shí)，振幅比會(huì)過小，相位滯后會(huì)過大，影響模型計(jì)算，因此需要預(yù)估區(qū)內(nèi)地下水流速以設(shè)置合理的測點(diǎn)間距。如果采用日波動(dòng)信息計(jì)算地下水流速，測點(diǎn)的埋深不宜過大，因深部地質(zhì)體中可能無法檢測到日波動(dòng)信息。如果需要評(píng)價(jià)河床的滲透性，可以同步測量測點(diǎn)間的水頭差，由地下水流速和水力梯度計(jì)算介質(zhì)滲透系數(shù)。而日波動(dòng)信息適用于評(píng)價(jià)淺部河床的滲透性，季節(jié)波動(dòng)或年波動(dòng)適用于評(píng)價(jià)深部河床的滲透性。

4 工程實(shí)例

4.1 實(shí)測資料及其處理

為評(píng)價(jià)美國Russian河地表水資源與地下水資源的交換，索諾馬水文局（Sonoma County Water Agency）以及美國地質(zhì)調(diào)查局西區(qū)分部（Branch of Regional Research，Western Region ，USGS）在Russian河河床埋設(shè)了溫度及水頭傳感器，連續(xù)觀測河水及地下水的壓力與溫度變化［21］。此處選取日周期波動(dòng)特征較為顯著的某一測點(diǎn)的資料來分析。

首先將帶濾網(wǎng)的PVC護(hù)筒打入河床滲透性介質(zhì)，然后將溫度探頭懸掛在護(hù)筒中，分別測量0.00 m和0.30m處的水溫。測試過程由數(shù)據(jù)記錄儀控制、記錄和存儲(chǔ)，每15min記錄一次。

采用動(dòng)態(tài)諧波回歸分析法（DHR法）［13］對(duì)實(shí)測資料進(jìn)行濾波處理，設(shè)置基頻為P＝12h，得到典型的按正弦簡諧波波動(dòng)的溫度時(shí)序資料，如圖4所示。DHR回歸分析的計(jì)算值與實(shí)測值的擬合效果非常好，非周期因素具有若干“突變點(diǎn)”，反映水動(dòng)力條件是不穩(wěn)定的。采用數(shù)值方法提取了濾波之后的溫度時(shí)序資料的振幅，即圖中“基頻振幅”所示，深部（0.30m）溫度振幅與淺部（0.00m）溫度振幅的比值為振幅比Ar。深淺2個(gè)測點(diǎn)濾波后的正弦曲線各極大值對(duì)應(yīng)的相位滯后同樣作為解析模型的輸入變量，用于相位滯后法計(jì)算地下水流速。

圖4 測點(diǎn)溫度時(shí)序資料DHR濾波分析成果Fig.4 Original and filtered temperature records

4.2 地表水與地下水交換速率

研究河段河床是未固結(jié)的砂礫石夾薄層淤泥或黏土。參照文獻(xiàn)［19］，取巖土體的熱物理參數(shù)如表1所示。采用2.2節(jié)所述之Hatch解析模型，得到觀測時(shí)段內(nèi)地表水與地下水交換速率的變化曲線，振幅比法的計(jì)算結(jié)果如圖5所示，相位滯后法的計(jì)算結(jié)果如圖6所示。研究時(shí)段內(nèi)測點(diǎn)附近的地下水流速是動(dòng)態(tài)波動(dòng)的，具有如下特征：

1）振幅法與相位滯后法計(jì)算的地下水垂向流速差異不大，波動(dòng)形態(tài)一致。當(dāng)?shù)叵滤a(bǔ)給河水時(shí)，相位滯后很小或超過π/2，無法采用解析模型計(jì)算，所以相位滯后法不能反映向上水流（圖6）。而且相位滯后法只能反映地下水流速的大小，不能確定地下水流動(dòng)方向，所以計(jì)算時(shí)間后段，振幅法確定的地下水流速存在負(fù)值，而相位滯后法無法確定正負(fù)。得到近似連續(xù)變化的河床地下水流速變化過程，對(duì)確定地表水與地下水的補(bǔ)排關(guān)系與補(bǔ)排量具有重要的應(yīng)用價(jià)值。

圖5 振幅比法計(jì)算的地下水垂向流速Fig.5 GW－SW exchange velocity calculated from amplitude damping method

圖6 相位滯后法計(jì)算的地下水垂向流速Fig.6 GW－SW exchange velocity calculated from phase lagging method

5 結(jié)語

1）對(duì)河床測溫資料進(jìn)行濾波處理，根據(jù)其振幅比和相位滯后可以計(jì)算地下水垂向流速，該方法可以方便快捷地評(píng)價(jià)地表水與地下水交換過程。

2）相對(duì)于Hatch解來說，Keery解忽略熱彌散效應(yīng)能引入較大的誤差。當(dāng)β＜0.01m時(shí)，Keery解忽略彌散效應(yīng)的處理是可行的；當(dāng)β＞0.10m時(shí)，Keery解所產(chǎn)生的誤差非常大。

3）只有研究區(qū)能簡化為均質(zhì)多孔介質(zhì)、地下水的1D垂向運(yùn)動(dòng)占主導(dǎo)地位時(shí)，才能采用解析法確定地下水流速；采用常規(guī)的測點(diǎn)距離（1.0m左右），振幅比法能有效評(píng)價(jià)流速在（－1.0～8.0）×10－5m/s的地下水活動(dòng)，而相位滯后法能有效評(píng)價(jià)流速在（0.0～1.6）×10－4m/s的地下水活動(dòng)。

4）后續(xù)的研究還需要進(jìn)一步討論多維流、巖土體非均質(zhì)性對(duì)解析模型的影響。

2）7月17日之前是地下水補(bǔ)給河水，地下水流速較小且小幅波動(dòng)，平均流速約為－1.0×10－6m/s。此時(shí)段河水水位較低，接受兩岸地下水的補(bǔ)給。

3）7月17日—8月8日地下水流速總體上緩慢增大，震蕩劇烈，大部分時(shí)段河水補(bǔ)給地下水，最大水流流速約為3.0×10－6m/s。地下水流速的動(dòng)態(tài)變化曲線上存在2個(gè)峰值，分別是7月22日和8月4日，此時(shí)段河水受降雨影響，水位上升，地表水與地下水之間的水頭差增大，地下水流速比較大。

上述分析表明，采用溫度時(shí)序曲線解析法能夠

（References）：

［1］王文科，楊澤元，程?hào)|會(huì)，等.面向生態(tài)的干旱半干旱地區(qū)區(qū)域地下水資源評(píng)價(jià)的方法體系［J］.吉林大學(xué)學(xué)報(bào)：地球科學(xué)版，2011，41（1）：137－144.Wang Wenke，Yang Zeyuan，Cheng Donghui，et al.Methods of Ecology－Oriented Groundwater Resource Assessment in Arid and Semi－Arid Area［J］.Journal of Jilin University：Earth Science Edition，2011，41（1）：137－144.

［2］方燕娜，廖資生，陳洪艷，等.吉林中部平原區(qū)地下水溫動(dòng)態(tài)變化影響因素的關(guān)聯(lián)度分析［J］.吉林大學(xué)學(xué)報(bào)：地球科學(xué)版，2006，36（1）：66－72.Fang Yanna，Liao Zisheng，Chen Hongyan，et al.Correlativity Analysis of the Affect Factors on Groundwater Thermal Regime in the Plain Area of Central Jilin［J］.Journal of Jilin University：Earth Science Edition，2006，36（1）：66－72.

［3］Constantz J.Heat as a Tracer to Determine Streambed Water Exchanges［J/OL］.Water Resources Research，2008，44，W00D10，doi：10.1029/2008WR006996.

［4］Stallman R W.Steady One－Dimensional Fluid Flow in a Semi－Infinite Porous Medium with Sinusoidal Surface Temperature［J］.Journal of Geophysical Research，1965，70（12）：2821－2827.

［5］Anderson M P.Heat as a Ground Water Tracer［J］.Ground Water，2005，43（6）：951－968.

［6］Taniguchi M.Evaluation of Vertical Groundwater Fluxes and Thermal Properties of Aquifers Based on Transient Temperature－Depth Profiles［J］.Water Resources Research，1993，29（7）：2021－2026.

［7］Christian A，Kerst B，Ronny V，et al.A Simple Thermal Mapping Method for Seasonal Spatial Patterns of Groundwater－Surface Water Interaction［J］.Journal of Hydrology，2011，397（1/2）：93－104.

［8］Niswonger R G，Prudic D E，F(xiàn)ogg G E，et al.Method for Estimating Spatially Variable Seepage Loss and Hydraulic Conductivity in Intermittent and Ephemeral Streams［J/OL］.Water Resources Research，2008，44： W05418，doi：10.1029/2007WR006626.

［9］Brewster C J.Delineating and Quantifying Groundwater Discharge Zones Using Streambed Temperatures［J］.Ground Water，2004，42（2）：243－257.

［10］吳志偉，宋漢周.壩址溫度場與變物性滲流場全耦合分析［J］.水利學(xué)報(bào)，2010，41（6）：703－710.Wu Zhiwei，Song Hanzhou.Fully－Coupled Analysis of Temperature Field and Variable Properties Seepage Field Around Dam－Site［J］.Journal of Hydraulic Engineering，2010，41（6）：703－710.

［11］Su G W，Jasperse J，Seymour D，et al.Estimation of Hydraulic Conductivity in an Alluvial System Using Temperatures［J］.Ground Water，2004，42（6）：890－901.

［12］Hatch C E，F(xiàn)isher A T，Revenaugh J S，et al.Quantifying Surface Water－Groundwater Interactions Using Time Series Analysis of Streambed Thermal Records：Method Development［J］.Water Resources Research，2006，42（10）：W10410.

［13］Keery J，Binley A，Crook N，et al.Temporal and Spatial Variability of Groundwater－Surface Water Fluxes：Development and Application of an Analytical Method Using Temperature Time Series［J］.Journal of Hydrology，2007，336（1/2）：1－16.

［14］Hatch C E，F(xiàn)isher A T，Ruehl C R，et al.Spatial and Temporal Variations in Streambed Hydraulic Conductivity Quantified with Time－Series Thermal Methods［J］.Journal of Hydrology，2010，389（3/4）：276－288.

［15］Cordon R P，Lautz L K，Briggs M A，et al.Automated Calculation of Vertical Pore－Water Flux from Field Temperature Time Series Using the VFLUX Method and Computer Program［J］.Journal of Hydrology，2012，420/421：142－158.

［16］Lautz L K.Observing Temporal Patterns of Vertical Flux Through Streambed Sediments Using Time－Series Analysis of Temperature Records［J］.Journal of Hydrology，2012，464/465：199－215.

［17］吳志偉，宋漢周.淺層地溫場中熱對(duì)流數(shù)值模擬［J］.巖土力學(xué)，2010，32（4）：1303－1308.Wu Zhiwei，Song Hanzhou.Numerical Simulation of Thermal Convection in Shallow Ground Temperature Field［J］.Rock and Soil Mechanics，2010，32（4）：1303－1308.

［18］Cardenas M B.Thermal Skin Effect of Pipes in Streambeds and Its Implications on Groundwater Flux Estimation Using Diurnal Temperature Signals［J/OL］.Water Resources Research，2010，12：W03536，doi：10.1029/2009WR008528.

［19］Healy R W，Ronan A D.Documentation of Computer Program VS2DH for Simulation of Energy Transport in Variably Saturated Porous Media［R］.Denver：US Geological Survey Information Center，1996.

［20］Weast R C.Handbook of Chemistry and Physics［M］.Boca Raton：Chemistry Rubra Pun Co Press，1982：261－263.

［21］Marisa H C，Hatch C E.Water Temperature，Stream Flow and Ground－Water Elevation in and Adjacent to the Russian River Between Hopland and Guerneville，CA，from 1998－2002［R］.Menlo Park：US Geological Survey Information Center，2003.