運用水熱運移模擬定量評價地表水與地下水的轉化

盧昱含，鮮 陽，靳孟貴，白宏偉，劉亞磊，劉延鋒

(1.中國地質大學(武漢)盆地水文過程與濕地生態恢復學術創新基地，湖北 武漢 430074；2.中國地質大學(武漢)環境學院，湖北 武漢 430074；3.中國地質調查局武漢地質調查中心,湖北 武漢 430205)

地表水與地下水聯系緊密，兩者之間的相互作用與轉化關系復雜，密切影響著潛流帶水分、物質與能量的交換過程[1-3]。因此，研究地表水與地下水的相互作用對于水資源評價、水資源開發管理、水污染防治和生態系統保護均具有重要的意義[4-5]。

干旱地區人類活動引起的土地利用和渠系滲漏的變化對地表水與地下水相互作用的影響強烈，而引水灌溉的渠道入滲是干旱地區地表水與地下水相互作用的主要形式[6]。近年來，干旱地區地表水與地下水的相互作用研究受到國內外研究機構的廣泛重視[7-8]。干旱地區水資源開發利用與調度配置的關鍵和基礎是判定河流和渠道與地下水的轉化關系，量化兩者間的轉化量及其時空變化規律[9]。

干旱地區淺層地下水水溫相對恒定，而河流等地表水水溫受氣溫的影響，一般具有較為明顯的季節性變化和日變化規律，地表水與地下水發生相互作用的過程中，潛流帶溫度場會發生明顯的變化[10-12]。近年來，隨著溫度場與滲流場耦合研究的深入、溫度數據監測成本的降低[13]以及多孔介質水熱耦合運移模擬軟件的大量開發，運用水熱運移數值模擬定量評價地表水與地下水相互轉化關系的研究有了較大的發展。如Su等[14]和Constantz[15]利用VS2DHI軟件模擬了河岸剖面二維溫度場，并結合觀測井中監測的季節性地下水溫度估算了含水層的水力傳導系數和地表水與地下水的相互轉化速率；任杰等[16]和潘維艷等[9]利用水熱運移模擬分別研究了大壩與農渠的滲漏問題；Ju等[17]結合沙箱試驗與數值模擬，分別采用實測溫度、水頭數據分析評價了數值模擬結果，結果表明水頭與溫度數據結合可提高滲流模擬結果的可信度和精度。但目前針對干旱渠道水文動態劇烈變化條件下渠道地表水與地下水相互轉化關系和轉化水量的時空變化規律方面的研究成果較少，且轉化量評價存在很大的不確定性。

本文以新疆瑪納斯河流域炮臺鎮附近頭道溝渠道的排水干渠作為研究對象，通過野外原位試驗實時監測渠道橫斷面地下水溫度、地下水水位、渠道地表水水位和不同埋深處河床沉積物的溫度變化，分析地表水與地下水相互轉化關系變化前后溫度的規律特征，并建立二維均質多孔介質變飽和土壤水熱運移耦合數值模型，模擬了渠道附近地下水的水熱運移，反演土壤飽和水力參數，定量計算渠道地表水與地下水的轉換量，為渠道滲漏定量估算和地表水與地下水聯合調度提供參考。

1 研究區域概況

瑪納斯河流域位于新疆天山北麓中段、準噶爾盆地南緣，地處北緯43°27′～45°21′、東經85°01′～86°32′，地勢自東南向西北緩緩傾斜。流域由東向西分別由塔西河、瑪納斯河、寧家河、金溝河、巴音溝河以及相連的五片沖積平原組成。從南向北地貌具有明顯的分帶性，可以劃分為5個地貌帶：山地帶、山前褶皺丘陵帶、串珠狀沖積扇帶、平原曲流帶、瑪納斯河尾閭湖泊和沙漠。

流域屬典型的干旱內陸盆地氣候，冬冷夏熱，年(日)溫差大，光照充足，蒸發強烈[18]。該流域年平均氣溫為9.1℃；年平均降水量為110～200 mm，降水分布不均,垂直分帶特征明顯，主要集中在5～8月份，約占年降雨量的67%；年均蒸發量為1 500～2 000 mm[19]。流域地表水補給主要來自于降水及冰川融化，地下水補給主要來自于通過河流的垂直滲透、渠道滲漏、灌溉回灌、山前側向地下水補給[20-21]。該地區平原區地表水與地下水的相互作用同時受到自然條件和人類活動的影響，兩者之間的關系極為復雜。

頭道溝渠道的排水干渠位于炮臺鎮附近(N44°45′～44°50′，E85°30′～85°35′)，地處平原灌溉區，見圖1。

圖1 研究區地理位置Fig.1 Location of the research area

根據研究區左右岸2m深鉆孔的土壤巖性數據，對比白宏偉等[22]對瑪納斯河流域二維巖性結構剖面的研究成果，認為研究區土壤質地分布較為均一，大部分為粉質壤土，受渠道及氣候的影響，入滲水水溫較低，低溫水入滲形成的溫度變化信號較為明顯。頭道溝渠道兼顧排水與灌溉功能，每年4月初上游開始向渠道放水，9月底停止供水，渠道地表水水位呈現規律性的升降變化。

2 研究方法

2.1 溫度和水位監測

頭道溝渠道屬未襯砌的典型農渠，本研究選定垂直頭道溝渠道的斷面，測定渠道斷面幾何形態，并布置監測井、水位和溫度傳感器。將監測井G1底部以下3.94 m處水平線定為基準線，頭道溝渠道監測斷面和監測井布置見圖2。

圖2 頭道溝渠道監測斷面和監測井布置Fig.2 Layout of the monitoring cross sections and wells of Toudaogou canal

監測井管為4 m長的PVC管，濾水管長為0.5 m。河岸監測井(G1、G2、G3、G4)的濾水管布置在監測井底部以上1 m處，每個監測井內均布置一個Solinst Solinst Levelogger 2參數探頭(壓力、溫度)自計水位和溫度。水位測量精度為±0.05 m，分辨率為0.001 m；溫度測量精度為±0.2℃，分辨率為0.01℃。

河床監測井(G5、G6)的濾水管位于河床下1 m處，其中監測井G5內布置一個Solinst Levelogger 2參數探頭自計水位和溫度；監測井G6內布置一個HOBO溫度記錄儀記錄30 cm、50 cm、70 cm、100 cm埋深處河床沉積物的溫度。水位測量精度為±0.05 m，分辨率為0.001；溫度測量精度為±0.2℃，分辨率為0.02℃。

除圖示監測井外，另布置一個Solinst Levelogger 2參數探頭監測渠道地表水水位和溫度，并使用數據采集儀記錄各測點的溫度和水位，采集時間間隔為15 min，自2015年7月27日22∶00至9月2日22∶00連續監測。

2.2 水文地質概念模型

由于模型水力參數受河床沉積物巖性結構和飽和度變化的影響可能會呈現幾個數量級的變化[10]，單純依靠水位校正數值模型會使得模型水力參數率定存在較大的不確定性，而熱傳導系數的變化范圍較小。因此，本文通過建立水熱運移耦合模型，利用水頭和溫度數據約束和校正模型水力參數，以降低模型的不確定性。

地表水與地下水相互作用時，土壤中水分運移與熱量傳輸會產生相互影響，土壤中水分運移會引起土壤的體積熱容量和熱傳導系數的變化，改變熱量傳輸與交換，進而改變土壤中溫度的分布。由于兩者之間的耦合作用過程非常復雜，為了研究渠道地表水與地下水的水熱運移規律，可忽略溫度變化引起水的性質變化，將監測斷面概化為二維均質多孔介質變飽和水熱運移耦合模型。

以頭道溝渠道斷面形態為上邊界，以監測井G1埋深10.06 m處水平線為底邊界，模型深度約為10 m；以監測井G1、G4為模型左、右邊界，模型水平長度為31.56 m。建立的頭道溝渠道研究斷面水文地質概念模型，見圖3。

圖3 頭道溝渠道研究斷面水文地質概念模型Fig.3 Conceptual hydrogeological model of the research cross section of Toudaogou canal

滲流場邊界條件：地表邊界B1、B3為通量邊界 (Neumann邊界)，監測期間幾乎無降雨，邊界賦值為結合經驗公式和實測水面蒸發E0數據計算出的邊界蒸發通量，由付秋萍等[23]研究可知，清華大學雷氏公式在新疆地區具有廣泛的適用性，參考該研究中粉壤土參數擬合成果，以函數形式賦值輸入；渠道邊界B2為已知變水頭邊界，輸入實測渠中水位值(Dirichlet邊界)；左右邊界B4、B5為已知變水頭邊界，分別輸入井孔G1、G4的實測水頭值(Dirichlet邊界)；底邊界B6為隔水邊界(Neumann邊界)。初始條件為實測壓力水頭數據的線性插值。

溫度場的邊界條件：地表邊界B1、B3為變溫度邊界，實測日平均氣溫數據(Dirichlet邊界)；渠道邊界B2為變溫度邊界，輸入實測渠道水溫數據(Dirichlet邊界)；底邊界B6為定溫度邊界，T=10℃(Dirichlet邊界)。初始條件為實測溫度數據的線性插值。

2.3 數學模型

2.3.1 土壤水分運動方程

假定研究斷面土壤為均質各向同性，設z軸垂直向上，x軸水平向右，二維均質各向同性土壤水分運動的偏微分方程及其定解條件可表示如下：

(1)

式中：t為時間變量(s)；x、z均為空間坐標(m)；h(x,z,t)為壓力水頭(m)；c(h)為比水容量,c(h)=?/?h；K(h)為非飽和土壤水力傳導系數(m/d)，是h的函數；h0(x,z)為滲流場初始地下水水位(m)；h1(x,z)為第一類邊界Γ1水位值(m)；ε(x,z,t)為第二類邊界Γ2通量(地表蒸發量)(m/d)；Ω為滲流區域。

2.3.2 土壤熱運移方程

采用熱傳導和對流方程來表征土壤熱運移過程，忽略氣態水擴散，僅考慮液態水的運動對土壤熱量傳輸的影響，二維均質多孔介質變飽和土壤的熱運移偏微分方程及其定解條件可表示如下：

(2)

2.3.3 土壤水力參數

土壤含水率θ、土壤有效飽和度Se、非飽和土壤水力傳導系數K這些參數均隨著壓力水頭h的變化而發生改變，需要通過數學模型來表征它們的本構關系。本文選用Mualem-van Genuchten公式[24]進行描述，具體公式可表示如下：

(3)

(4)

(5)

式中：θr為土壤殘余體積含水率(m3/m3)；n、m和α為經驗擬合參數，其中m=1-1/n，α為受土壤物理性質影響的參數(1/m)；l為經驗擬合參數，通常取平均值0.5；Se為土壤有效飽和度；Ks為土壤飽和水力傳導系數(m/d)。

2.4 數值模型建立與參數識別

采用地下水數值模擬軟件FEFLOW建立二維均質多孔介質變飽和土壤水熱運移耦合數值模型。使用advancing front剖分網格，得到總網格數為3 178、節點數為1 683。設定邊界條件，并設置土壤水力參數和土壤介質熱特性參數初值，采用自動時間步長劃分離散時間，使用線性插值結果作為變飽和土壤非穩定滲流場和溫度場的初始值。

數值模擬采用的土壤水力學參數：在距離監測井G2和G3正北1 m處，鉆孔采集0 cm、-30 cm、-50 cm、-70 cm、-100 cm、-130 cm、-150 cm、-170 cm、-200 cm深度處的土壤樣品，采用中國地質大學(武漢)材料與化學分析測試中心的Mastersizer 3000E激光衍射粒度分析儀測定土壤的粒度組成，結果表明土壤的粒度組成主要為粉粒(51.8%)和砂粒(43.84%)，并參考美國農業部土壤質地分類方法，將其定名為粉質壤土；環刀采集左右岸0 cm、-50 cm、-100 cm深度處土壤樣品，測定土壤的干容重為1 650 kg/m3，利用Hydrus1D中內嵌的Rosetta Lite V1.1模塊，使用土壤粒度組成和干容重數據確定初始土壤水力參數。

本次選用為期12 d(2015年7月27日22∶00至8月8日22∶00)的實測數據率定模型參數，使用水頭和溫度數據反演土壤水力參數，校正的模型參數見表1。

表1 反演的土壤水力參數Table 1 Inversed soil hydraulic parameters

數值模擬采用的土壤介質熱特性參數有：干土的比熱容Cs、干土的熱導率λs、水的比熱容Cw、水的熱導率λw、土壤縱向熱擴散度DL、土壤橫向熱擴散度DT。利用Lu等[25]研究模型中的常溫土壤熱導率公式計算干土熱導率；參考Bridger等[26]和龍燾等[27]的研究成果設定干土比熱容的數值；參考張盼盼[28]的研究數據設定土壤熱運移的縱向擴散度和橫向擴散度，土壤介質的熱特性參數見表2。

表2 土壤介質的熱特性參數Table 2 Soil thermal parameters

3 結果與討論

保持率定后的模型參數不變，使用為期37 d(2015年7月27日22∶00至9月2日22∶00)的野外監測數據驗證模型的可靠性。本文利用均方根誤差RMSE和決定系數R2來評價模擬值與實測值的吻合程度，其計算公式如下：

(6)

(7)

不同埋深處河床沉積物溫度模擬值與實測值的對比結果，見圖4。

由圖4可見，30 cm、50 cm、70 cm、100 cm埋深處河床沉積物溫度模擬值與實測值的差異較小，總體趨勢一致，且隨著埋深的增加，溫度波動的振幅減小；溫度模擬值與實測值的總均方根誤差RMSE為0.220，決定系數R2均大于0.95，認為河床沉積物溫度梯度的模擬結果合理可靠。

圖4 不同埋深處河床沉積物溫度模擬值與實測值的對比Fig.4 Comparisons between simulated and observed bed sediment temperature at different depths

河岸監測井G2、G3水頭模擬值與實測值的對比結果，見圖5。

圖5 河岸監測井G2、G3水頭模擬值與實測值的對比Fig.5 Comparisons between simulated and observed hydraulic head for monitoring wells G2 and G3

由圖5可見，盡管河岸監測井G2、G3處水頭的模擬值與實測值在某些觀測點上存在一定的誤差，但總體趨勢基本一致；水頭模擬值與實測值的總均方根誤差RMSE為0.033，決定系數R2均大于0.98，認為滲流場的模擬結果合理可靠，部分誤差可能是由實際渠道斷面土壤的空間變異性造成的。

河岸監測井G2、G3溫度模擬值與實測值的對比結果，見圖6。

由圖6可見，河岸監測井G2、G3溫度模擬值與實測值的差異較小，總體趨勢一致；總均方根誤差RMSE為0.115，決定系數R2均大于0.98，認為溫度場的模擬結果合理可靠。

綜上分析可見，研究區地溫和水溫的波動大，有明顯的晝夜變化。渠道地表水與地下水相互作用時，渠道地表水的溫度信號影響著河床沉積物的溫度分布。本文針對不同埋深處河床沉積物的溫度數據，提取其日周期波動的正弦函數擬合振幅值見圖7。

圖6 河岸監測井G2、G3溫度模擬值與實測值的對比Fig.6 Comparison of simulated and monitored temperatures for monitoring wells G2 and G3

圖7 不同埋深處河床沉積物溫度曲線正弦函數擬合振幅Fig.7 Sinusoidal fitting amplitudes of bed sediment temperature curves at different depths

通過分析圖4和圖7可見，河床沉積物的平均溫度從上至下逐漸降低，埋深小于30 cm處河床沉積物的溫度受到渠道水溫波動的影響較大，在埋深為50 cm及以下處河床沉積物的溫度受渠道水溫波動的影響不明顯；當渠道地表水補給地下水時，渠道水溫日周期波動信號向下傳遞，河床沉積物溫度日周期波動的振幅(正弦函數擬合)隨埋深增加而減小；當地下水補給渠道地表水時，低溫地下水向上進行熱量交換，緩沖了由渠道地表水引起的河床沉積物溫度波動，50 cm埋深處河床沉積物的溫度曲線信號隨之變化，其日周期正弦波動的振幅減小。可見，通過分析50 cm埋深以上處河床沉積物的溫度信號，即使在無水位數據的情況下，也能分析得出渠道地表水與地下水的相互作用關系。

為了研究渠道地表水與地下水的相互作用對渠道垂直剖面溫度分布規律的影響，本文利用上述模擬結果，選取2015年8月12日22∶00、8月24日22∶00、8月25日22∶00、8月29日22∶00，即模擬時長分別為t=16 d、t=28 d、t=29 d和t=33 d 這4個不同時刻的溫度場模擬結果，對比分析渠道地表水與地下水不同轉化關系時，渠道研究斷面的溫度分布規律，見圖8。

圖8 不同時刻渠道研究斷面的溫度分布Fig.8 Distribution of temperature at the research cross section of the cannal at different times

由圖8可見，渠道研究斷面的溫度分布由下至上大致可以分為低溫區(10～17℃)、中溫區(17～20℃)和高溫區(20～23℃)。其中，低溫區分布較為穩定，無論斷面是渠道地表水補給地下水還是地下水補給渠道地表水，低溫區的面積變化不大；高溫區分布于受氣溫、水溫影響較大的淺部非飽和層，其分布范圍較小，面積變化差異大；中溫區位于兩者之間，受到渠道地表水與地下水相互作用中熱傳導作用的影響。

當渠道地表水補給地下水時，受溫度日變異性大的渠道地表水入滲的影響，中溫區的溫度梯度整體隨埋深的增加而逐漸增大，但局部的溫度梯度方向可能相反；當地下水補給渠道地表水時，中溫區的溫度梯度方向一致，隨著埋深的增加中溫區的溫度梯度逐漸減小。這在一定程度上表征了渠道研究斷面渠道地表水與地下水的相互作用關系。

為了定量評價渠道地表水與地下水的轉化水量，本文選取渠道邊界單元，繪制出渠道研究斷面單位長度累計入滲量的變化曲線，見圖9。

圖9 渠道研究斷面單位長度累計入滲量的變化曲線Fig.9 Cumulative recharge per unit length of the research cross section of the canal

由圖9可見，0～28.6 d內渠道研究斷面單位長度累計入滲量隨時間增大，渠道地表水補給地下水，隨著渠道地表水水位的降低，渠道研究斷面單位長度累計入滲量曲線斜率減小；28.6～37 d內渠道研究斷面單位長度累計入滲量隨時間減小，渠道地表水水位低于地下水水位，轉化為地下水向渠道排泄，隨著地下水水位的降低，渠道單位長度累計入滲量曲線斜率的絕對值減小。0～28.6 d內渠道單位長度累計入滲量為35.3 m2/d；28.6～37 d內渠道研究斷面單位長度累計排泄量為9.626 m2/d；37 d內渠道單位長度累計入滲量為25.674 m2/d。

假設渠道研究斷面為典型斷面，可利用渠道研究斷面單位長度累計入滲量估算整個渠道地表水與地下水的轉化量。已知研究區斷面位于新疆石河子農八師瑪納斯河灌區頭道溝渠道的排水干渠上，干渠全長為12.6 km，南北走向，垂直等高線分布,控制灌溉面積為2萬hm2。通過分析得出:2015年7月27日22∶00至2015年8月25日12∶30，頭道溝干渠渠道地表水入滲補給地下水，地表水入滲補給水量約為4.448×105m3，2015年8月25日12∶30至9月2日22∶00，頭道溝干渠排泄地下水，地下水排泄水量約為1.213×105m3。

4 結 論

本文以新疆瑪納斯河流域炮臺鎮附近頭道溝渠道的排水干渠為研究對象，通過野外原位試驗實時監測渠道橫斷面地下水溫度、地下水水位、渠道地表水水位和不同埋深處河流沉積物的溫度變化，分析了渠道地表水與地下水相互轉化關系變化前后溫度的規律特征，并建立二維均質多孔介質變飽和土壤水熱運移數值模型，模擬分析與計算了頭道溝干渠渠道地表水與地下水的相互轉化關系和轉化量，并得到如下結論：

(1) 淺層河床沉積物(埋深為50 cm以上)的溫度曲線可示蹤表征渠道地表水與地下水的相互轉化關系。在觀測計算時段內，2015年7月27日22∶00至2015年8月25日12∶30，渠道地表水補給地下水，2015年8月25日12∶30至9月2日22∶00，地下水向渠道排泄。河床沉積物可大致分為低、中、高三個溫區，當渠道水補給地下水時，受渠道地表水入滲的影響，中溫區溫度梯度隨埋深的增加而逐漸增大，可能出現淺層沉積物溫度低于深層沉積物溫度的現象；當地下水向渠道排泄時，中溫區溫度梯度方向一致，隨著埋深的增加溫度梯度逐漸減小。河床沉積物的溫度日周期波動振幅(正弦函數擬合)隨著埋深的增加而減小，當渠道地表水補給地下水時，50 cm埋深處河床沉積物的溫度日周期波動振幅隨時間增大；當地下水向渠道排泄時，50 cm埋深處河床沉積物的溫度日周期波動振幅隨時間減小。

(2) 使用土壤水熱運移耦合模型反演得到試驗區土壤的飽和水力傳導系數為7.05 m/d。使用監測數據驗證模擬結果，溫度模擬值與實測值的RMSE為0.220，決定系數R2均大于0.95，水頭模擬值與實測值的RMSE為0.033，決定系數R2均大于0.98，模擬值和實測值之間具有較好的一致性，表明率定后的土壤水熱運移模型能較好地仿真渠道地表水與地下水相互轉化的動態變化規律。

(3) 渠道研究斷面單位長度累計入滲量為35.3 m2/d，發生在2015年8月25日12∶30時，模擬結束時渠道研究斷面單位長度累計入滲量為9.626 m2/d。2015年7月27日22∶00至2015年8月25日12∶30期間，頭道溝干渠渠道地表水入滲補給地下水，12.6 km長渠道地表水入滲補給量約為4.448×105m3；2015年8月25日12∶30至9月2日22∶00期間，頭道溝干渠排泄地下水，地下水排泄量約為1.213×105m3。