潘維艷，楊姍姍，劉俊鋒，錢秀紅，徐征和

（1.濟南大學(xué)水利與環(huán)境學(xué)院，山東 濟南 250022；2.山東省海河淮河小清河流域水利管理服務(wù)中心，山東 濟南 250014）

0 引 言

潛流帶是流域系統(tǒng)動態(tài)水文過程與生物地球化學(xué)反應(yīng)的關(guān)鍵帶，潛流帶內(nèi)地表水與地下水的相互作用直接影響水在河床沉積帶內(nèi)的滯留時間、沉積帶營養(yǎng)物質(zhì)的通量、河床孔隙水的酸堿性及氧化還原條件［1，2］。在河岸潛流帶內(nèi)，地表水通過河岸沉積層與地下水發(fā)生物質(zhì)和能量交換，生物活動強烈。地表水入滲過程中，水分運移過程能夠攜帶污染物進入地下水中，而潛流帶對地表水中的污染物具有重要的截污作用，因此，潛流帶對地下水水量和水質(zhì)變化過程發(fā)揮著重要作用［1，3］。潛流帶也已成為近年來國內(nèi)外學(xué)者研究的熱點［4-6］。因此，定量分析潛流帶水量與流速的時空變化特征是充分認(rèn)識地表水-地下水交互作用規(guī)律的關(guān)鍵問題。

潛流帶中時刻伴隨著水熱交換，溫度作為載體，可以直觀反映地表水和地下水流動過程及其交換關(guān)系，地表水與地下水交互作用的研究已經(jīng)從水文學(xué)和水動力學(xué)法發(fā)展到溫度示蹤法［7，8］，溫度示蹤法檢測成本低、靈敏度高、數(shù)據(jù)獲得直接、穩(wěn)定且無污染，在示蹤地下水流速方面具有明顯優(yōu)勢［9-11］。近年來，溫度監(jiān)測手段的進步和溫度場研究的深入，利用溫度示蹤法進行流速定量計算的研究越來越多，如Caissie 等［12］運用溫度示蹤法評估并分析了季節(jié)性冰蓋河流的地下水的流速和方向，張文兵等［13］通過溫度示蹤的潛流交換通量解析模型對比研究發(fā)現(xiàn)，淺層潛流帶中潛流交換頻繁，且深度越大，交換作用減弱。喬曉英等［14］將溫度示蹤應(yīng)用于湖水與地下水交互作用，發(fā)現(xiàn)基于溫度示蹤的數(shù)值模擬法結(jié)果更理想。董林垚等［15］利用溫度示蹤法測算地下水流速發(fā)現(xiàn)地質(zhì)體擴散率和地表溫度是影響結(jié)果準(zhǔn)確性的主要因子。溫度是研究河流等地表水與地下水交換的較為有效的方法之一，但不同研究帶的水文地質(zhì)特征和水動力條件的不同對地下水流速有不同影響，目前，基于溫度示蹤法探討河岸潛流帶的研究成果還相對較少，仍需加強。

因此，本文選取再生水補給河湖的典型河道為研究區(qū)，定期監(jiān)測地表水、地下水的水位和水溫，利用水溫和水位數(shù)據(jù)資料分析潛流帶溫度場的變化規(guī)律及其影響因素，并利用水動力法和溫度示蹤法計算流速，進行對比分析，探討采用溫度示蹤法用于潛流帶的適用性和合理性。再生水回灌區(qū)潛流帶滲漏研究是再生水安全、高效回灌的關(guān)鍵，可以為再生水資源配置及地下水環(huán)境保護提供理論依據(jù)和決策支持。

1 研究區(qū)概況

本研究的實驗場地選擇在北京市朝陽區(qū)奧林匹克森林公園內(nèi)的清洋河河岸（40°00′～40°02′N，116°22～116°24′E）（圖1）。研究區(qū)多年平均氣溫11～12 °C，平均水面蒸發(fā)量在1 200 mm 左右，平均降雨量約600 mm，降雨年內(nèi)分配不均，85%以上主要集中在汛期。清洋河河道補水全部來自再生水、降雨和徑流，再生水補水期為每年的3-11月，年供給水量約1.5×106m3。清洋河是人工河道，采樣分析發(fā)現(xiàn)研究區(qū)河岸帶的土壤屬于壤土，河床底部有防滲處理。

圖1 試驗區(qū)位置圖Fig.1 Location of study area

2 材料與方法

2.1 試驗布置與數(shù)據(jù)采集

在河岸帶選取一處典型斷面，沿垂直河水方向，在斷面上布置1 處河水溫度觀測點，在距離河岸位置2.5 和5 m 處分別布置 2 組觀測井（W1 和W2），每眼井分別在距離地面0.5、1、1.5、2、2.5、3、4、5、6、7、8、9 m 處安裝熱電偶熱線，以觀測包氣帶∕飽和帶各深度處的溫度，并在2 眼井附近布設(shè)水位探頭（HOBO U20L-01，美國）。利用數(shù)據(jù)采集器（CR3000，美國）采集并記錄各觀測點溫度，采集頻率設(shè)置為4 次∕h，自2013年9月25日至2014年1月20日進行了連續(xù)觀測。試驗布置和現(xiàn)場試驗圖如圖2所示。采用土鉆法每月采集一次包氣帶土壤含水量，采樣深度為0.5、1、1.5和2 m，采集時間為2013年9月-2014年1月。

圖2 試驗監(jiān)測井布置和試驗現(xiàn)場圖Fig.2 Test monitoring well layout and test site diagram

2.2 入滲速率計算方法

采用的熱運移擴散方程為［16-18］：

式中：T為溫度，℃；t為時間，s；q為流體的流速，m∕s；Cw為水的體積熱熔，J∕（m3·℃）；C是飽和沉積物體積熱容，J（m3·℃）；z是深度，m；ke為飽和介質(zhì)有效熱擴散系數(shù)，m2∕s。

Hatch 等［16］提出了基于上述方程的流速解。Hatch 通過振幅或相位差計算的流速解析解分別為：

式中：q為垂向向下流速，m∕s；Ar為下側(cè)傳感器數(shù)據(jù)振幅與上側(cè)之比；Δz為傳感器間距離，m；v為溫度鋒速度，m∕s；Δt為信號傳播的時間差，s；P為溫度信號波動周期，s；α定義為：α=ke為等效熱擴散系數(shù)；ke定義為：β|Vf|；λ0為水流的熱傳導(dǎo)系數(shù)；β為熱彌散度；Vf為溫度前端運移速度，m∕s；H定義為H=C∕λ0，在上述解析解中通過振幅法計算的解析解包含流速的大小和方向，通過相位計算得到的解析解只有流速大小。

本文利用Hatch 相位和振幅兩種方法對溫度時間序列進行分析并計算不同深度處的水流流速。計算過程中所需的參數(shù)主要通過試驗和相關(guān)文獻料獲得，參數(shù)見表1。采用VFLUX2.0程序針對上述方法進行流速計算。VFLUX2.0 是一種基于matlab平臺，利用熱傳輸方程計算飽和多孔介質(zhì)中一維垂直流體流動（滲流通量）的程序［19］。文中選擇溫度波動相對強烈的時段為例，進行流速計算，計算時段為2013年11月7日至11月17日，共計240 h，Hatch位相法計算所用的溫度采用頻率為0.5 次∕h。水動力學(xué)法基于達西定律進行計算［20］，根據(jù)多點原位取樣發(fā)現(xiàn)研究區(qū)河岸潛流帶質(zhì)地比較均一，通過變水頭入滲試驗測得平均滲透系數(shù)，然后按照達西定律計算監(jiān)測井附近的流速［21］。

表1 相關(guān)計算參數(shù)Tab.1 Related calculation parameters

3 結(jié)果與討論

3.1 河水深度與地下水埋深變化特征

試驗期間河水深及地下水井W1 和W2 的地下水位埋深隨時間的變化如圖3所示。河水水深波動范圍為1.0～1.8 m，W1的地下水埋深波動范圍為在1.5～3.0m，W2 處的地下水埋深波動范圍為2.5～5.5 m。在觀測期內(nèi)，W1和W2處地下水埋深整體上呈增大趨勢。其中，W1和W2處地下水埋深在2016年10月5日至15日呈先減小后增大的趨勢，在2017年1月3日至1月23日地下水埋深呈現(xiàn)逐漸增大趨勢，W1 處地下水埋深從2.2 m 增大到2.8 m，W2 處地下水埋深從3.5 m 增大到5.1 m。從圖3中可看出，W1和W2的地下水埋深隨時間變化呈現(xiàn)相似的變化規(guī)律，地下水埋深的變化主要取決于河水位的波動，即河水水位減小導(dǎo)致地下水水位抬升，并呈現(xiàn)一定的滯后性。結(jié)果表明，在觀測期間，河水位持續(xù)高于地下水位，河道再生水持續(xù)補給地下水。

圖3 地下水埋深與河水水深隨時間變化曲線Fig.3 Temporal variation of groundwater table and river table

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

本次試驗對W1 和W2 兩處監(jiān)測井進行了持續(xù)近4 個月的溫度觀測，由于設(shè)備故障等原因造成部分時段數(shù)據(jù)出現(xiàn)缺失（如圖4、5 中兩處缺口部分）。根據(jù)采集的溫度數(shù)據(jù)顯示，2 眼井中的溫度數(shù)據(jù)變化規(guī)律和趨勢一致，僅在數(shù)據(jù)上有差異，選擇數(shù)據(jù)更加完善的W2 數(shù)據(jù)進行分析。圖4是W2 中不同各深度處溫度隨時間變化的曲線，其中，0.5、1、1.5和2 m觀測點長期處于地下水位以上，屬于包氣帶層的地溫，6、7、8 和9 m 四個觀測點長期處于地下水位以下，為地下水溫度。

由圖4可知，河岸包氣帶地溫在時間上，2016年9月-2017年1月，不同深度的包氣帶地溫隨時間逐漸減小，溫度波動范圍分別為23.4～3.4 ℃（0.5 m），24.5～6.0 ℃（1 m），24.8～7.8 ℃（1.5 m）和25.0～10.2 ℃（2 m），且整體上呈現(xiàn)線性下降趨勢（各R2均大于0.96）。

河岸包氣帶溫度在空間上，包氣帶平均地溫隨著深度的增加呈現(xiàn)逐漸增加的趨勢，平均溫度從0.5 m 處的10.8 ℃增加到2 m 處的16.9 ℃。觀測期內(nèi)，氣溫波動劇烈（圖4），包氣帶地溫（2 m 以內(nèi)）隨氣溫波動變化明顯，兩者呈顯著正相關(guān)關(guān)系（R2=0.87～0.89，圖5）。結(jié)果表明，氣溫是影響包氣帶土壤溫度的主要因素。在4 個月的觀測期內(nèi)（9月-次年1月），除個別時間點外，地溫均高于氣溫，其中，表層地溫最接近氣溫，隨著深度增大，地溫與氣溫的差異越大，地溫和氣溫的平均差異在5～11 ℃范圍內(nèi)。

圖4 W2中不同深度處溫度隨時間變化規(guī)律Fig.4 Variation of temperature at different depths with time in W2

圖5 W2處包氣帶地溫與氣溫之間的相關(guān)性Fig.5 Correlation between ground and air temperature in the vadose zone in W2

包氣帶不同深度處的土壤含水量與包氣帶地溫之間的相關(guān)性分析如圖6所示，結(jié)果發(fā)現(xiàn)，在0.5～2.0 m 范圍包氣帶土壤質(zhì)量含水率與包氣帶地溫之間呈現(xiàn)顯著的正相關(guān)關(guān)系（R2=0.680 4），包氣帶地溫隨著含水量的增加而增大，在觀測期內(nèi)，當(dāng)含水率在17%～26%內(nèi)變化時，相應(yīng)深度處的包氣帶溫度在1.9～25.0 ℃范圍內(nèi)變化，可見，土壤含水量是影響包氣帶溫度的因素之一。有研究指出降水也是影響包氣帶溫度的重要因素之一，由于試驗期間研究區(qū)降水量較少，故未考慮其影響。

圖6 土壤質(zhì)量含水率與包氣帶地溫之間的關(guān)系Fig.6 Relationship between soil mass moisture content and ground temperature in vadose zone

由圖7可知，河岸飽和帶地下水溫度呈現(xiàn)出與包氣帶地溫不同的時空變化特征。從時間上看，與包氣帶地溫變化規(guī)律不同，飽和帶地下水溫度整體上呈先現(xiàn)逐漸增加后逐漸下降的變化趨勢，但不同深度呈溫度變化的拐點時間略有差異，其中，6、7 和8 m 三處的地下水溫度自2016年9月25 至10月23 呈緩慢增加趨勢，隨后至次年1月中旬，溫度緩慢下降；9 m 處，地下水溫度自2016年9月25日至12月14日，呈緩慢增加趨勢，隨后至1月中旬又緩慢下降，不同時間節(jié)點溫度波動范圍分別為24.0～24.2～16.9 ℃（6 m），21.8～22.3～17.6 ℃（7 m），19.0～20.1～18.0 ℃（8 m），16.2～18.8～18.1 ℃（9 m）。4 處觀測點處的地下水平均溫度波動范圍為17.4～21.4 ℃，非飽和帶平均地溫波動范圍為6.8～24.4 ℃，可見，地下水溫度隨時間的波動性比非飽和有所減弱。

圖7 W2處地下水溫度隨時間的變化過程Fig.7 Variation of groundwater temperature with time in W2

空間上，從包氣帶到過渡帶到飽和帶，不同深度處平均溫度分別為10.8 ℃（0.5 m）、13.7 ℃（1 m）、15.4 ℃（1.5 m）、16.9 ℃（2 m）、18.2 ℃（3 m）、20.0 ℃（4 m）、21.1 ℃（5 m）、20.9 ℃（6 m），20.3 ℃（7 m），19.2 ℃（8 m），18.1 ℃（9 m）。在不同季節(jié)呈現(xiàn)出不同的變化趨勢，在秋季，溫度隨著深度增加呈現(xiàn)先增后減的趨勢，表現(xiàn)為上冷中暖下冷的分層現(xiàn)象；在冬季，溫度隨著深度增加呈現(xiàn)逐漸增加的趨勢，表現(xiàn)為上冷下暖的分層現(xiàn)象。在觀測期內(nèi)，地下水平均溫度隨深度變化規(guī)律與包氣帶平均地溫相反，即地下水平均溫度呈隨深度的增加而下降的趨勢。

3.3 水動力學(xué)法與溫度示蹤法計算結(jié)果對比

水動力學(xué)法和Hatch 相位法、振幅法計算得到的流速結(jié)果如圖8所示。水動力學(xué)法計算得到的流速值變化范圍為6.74～7.74×10-5m∕s，最大值出現(xiàn)在第102 h，流速整體變化幅度不大，這主要是因為該河段在試驗期間河水水位波動不大。Hatch 相位法計算的流速值變化范圍為6.34～8.36×10-5m∕s，最大值出現(xiàn)在114 h。

圖8 水動力學(xué)法和溫度示蹤法結(jié)果Fig.8 Results of hydrodynamic method and temperature tracer method

從數(shù)值上來看，水動力學(xué)法的計算結(jié)果與Hatch 相位法的計算結(jié)果相似。但是，與水動力學(xué)法相比，Hatch 相位法對計算的流速值變化波動范圍和變化幅度更大，且波動趨勢呈現(xiàn)一定的滯后性，這是因為溫度在多孔介質(zhì)中的傳到和擴散速度要滯后于地下流速［21］。Hatch 振幅法計算的流速值變化范圍為0.61～1.01×10-5m∕s，相位法計算的流速值要明顯大于振幅法，相差接近一個數(shù)量級，兩者的變化趨勢相似。振幅法計算結(jié)果為正值，表明河岸帶的滲透流速均為垂向向下，即河水補給地下水［22］，這與實際情況相符。對比河水位、地下水位與流速的關(guān)系還發(fā)現(xiàn)，當(dāng)河水與地下水水位差較大時，對應(yīng)的流速值相對較大，結(jié)果表明，當(dāng)河水補地下水時，河水與地下水的水位差對研究河段河水入滲率的變化有一定影響。根據(jù)水動力學(xué)法的計算結(jié)果可以發(fā)現(xiàn)，Hatch 相位法計算流速的準(zhǔn)確性要高于Hatch 振幅法，這與李英玉等［9］的研究結(jié)果相似，這主要取決于兩種計算方法解析解表達式的差異。因此，后文的比較分析均采用Hatch相位法計算的流速值結(jié)果。

3.4 河岸帶不同深度處流速對比

用Hatch 相位法計算得到的河岸帶不同深度處的q1流速（距地面2.5 m）、q2流速（距地面3 m）和q3流速（距地面4 m）如圖9所示。由圖9可知，3 個深度處的流速變化范圍為5.5～8.3×10-5m∕s，最大值出現(xiàn)在淺層q1流速處，最小值出現(xiàn)在深層q3流速處。不同深度流速的時程變化趨勢在大部分時段相似，但不同深度處流速值存在一定差異，總體來看，流速值隨著深度的增加呈減小趨勢，這是因為隨著深度增大，水壓消散和水頭損失增大，從而導(dǎo)致流速減小。

圖9 不同深度處流速隨時間的變化曲線Fig.9 Variation curves of flow velocity with time at different depths

對比q1流速處和q3流速處的變化趨勢發(fā)現(xiàn)，雖然不同深度處溫度的監(jiān)測頻率一樣，但溫度示蹤法對淺層流速變化的刻畫描述要比深層好，這是由熱傳導(dǎo)和熱擴散過程中溫度波動信號的衰減導(dǎo)致的［21］。因此，溫度示蹤法用于潛流層地下水流速的計算能夠較好的反映有水位波動引起的淺層地下水流速時程變化趨勢，并均有較好的準(zhǔn)確性。

4 結(jié) 論

（1）從包氣帶到過渡帶到飽和帶，溫度在不同秋季和冬季呈現(xiàn)不同的變化特征，在秋季，溫度隨著深度增加呈現(xiàn)先增后減的趨勢，表現(xiàn)為上冷中暖下冷的分層現(xiàn)象；在冬季，溫度隨著深度增加呈現(xiàn)逐漸增加的趨勢，表現(xiàn)為上冷下暖的分層現(xiàn)象。在觀測期內(nèi)，地下水溫度隨深度增加呈下降趨勢，這與河岸包氣帶地溫變化規(guī)律相反，且地下水溫度隨時間的波動性比包氣帶有所減弱。氣溫和土壤含水量與包氣帶地溫成正相關(guān)關(guān)系。

（2）分別采用水動力學(xué)法和Hatch 相位法、振幅法進行地下水流速計算，結(jié)果發(fā)現(xiàn)，水動力學(xué)法計算得到的流速值變化范圍為6.74～7.74×10-5m∕s，Hatch 相位法計算得到的流速值變化范圍為6.34～8.36×10-5m∕s，Hatch 振幅法計算得到的流速值變化范圍為0.61～1.01×10-5m∕s，在計算時段內(nèi)，3 種方法計算得到的流速的變化趨勢在大致相似，溫度示蹤法計算結(jié)果呈現(xiàn)一定的滯后性。對比發(fā)現(xiàn)，對于研究區(qū)來說，Hatch 相位法用于河岸帶潛流層地下水流速的計算更合理更準(zhǔn)確。

（3）比較不同深度處的地下水流速發(fā)現(xiàn)，深度越大，地下水流速越小，且溫度示蹤法對流速的刻畫描述效果越差。