地下水源熱泵運行期間地下水流場和溫度場的變化規律研究

徐連三，華 杉，劉紅衛，柯 立，李永輝，王富強

(1.湖北省地質局 武漢水文地質工程地質大隊,湖北 武漢 430051； 2.湖北省地熱能研究推廣中心,湖北 武漢 430051)

地下水源熱泵系統將地下水作為熱能載體,利用其溫度常年穩定,不易受外界環境干擾的特性,夏季制冷,冬季取暖[1]。在國家宏觀政策引導和地方政策支持下,地下水源熱泵系統工程發展迅速,同時也暴露出一些問題,如抽—灌井布局、回灌井回灌、熱貫通及系統運行對地下水環境的影響等。國內外一些學者針對如此多的問題開展了大量研究。如徐貴來等[2]通過實地對地下水源熱泵系統運行期間地溫變化的監測,分析研究巖土層溫度變化規律,提出了提高系統能源的措施和建議；王佳樂等[3]應用有限元模擬軟件模擬了巖溶含水層地下水源熱泵運行對地下水溫度場的影響；Park B等[4]通過試驗觀測發現注水井回灌導致地下水熱彌散系數變化,從而影響含水層滲流場和溫度場；張淑秘等[5]應用水熱耦合模型模擬了不同布井模式下含水層溫度變化和熱貫通特征。地下水源熱泵系統抽灌井布局,模擬預測地下水源熱泵受水文地質條件因素影響,各區域地下水源熱泵系統運行對地下水環境影響不盡相同,且運行模式不同影響亦不同。鑒于此,本研究以江漢油田供熱移交改造項目——向中能源站地下水源熱泵系統工程為例,利用MODFLOW軟件的MT3D模塊建立地下水流—熱運移耦合數值模擬模型,模擬地下水源熱泵運行期間地下水流場和溫度場的變化規律,旨在為后期地下水源熱泵系統優化設計,減輕系統運行對地下水環境的影響提供依據。

1 研究區概況

1.1 研究區水文地質條件

研究區位于湖北省潛江市廣華區(見圖1),距潛江市區約23 km,地貌上屬典型的平原湖區地貌,地勢西北部略高、東南部略低,海拔24～40 m。屬亞熱帶濕潤季風氣候區,四季分明,多年日平均氣溫18.3 ℃。

圖1 研究區地理位置圖Fig.1 Geographical location map of the study area

整體屬于漢江平原水文地質單元,200 m以淺主要為第四系松散巖類孔隙水,可劃分為兩個承壓含水層組：第一層承壓含水層埋深約13.0 m,厚度約為29.4 m,主要物質成分為粉砂、粉細砂,單井涌水量500～1 000 m3/d；第二層承壓含水層埋深約46.0 m,厚度約為78.3 m,主要物質成分為圓礫、中砂和細砂,單井涌水量2 000～3 200 m3/d。

區內地下水位埋深隨季節波動,枯水期地下水位埋深約為2.0 m,豐水期地下水位埋深0.5 m左右。地下水流向近似由北向南,根據水位標高初步測算區內天然水力坡度約為0.23‰。根據對工作區內地下水溫度的監測結果顯示,在150 m深度范圍內,地下水取水口溫度普遍在18.5～19.5 ℃。區內第四系地下水化學類型主要為HCO3-Ca·Na及HCO3-Ca·Mg,礦化度<1 g/L,屬于低礦化度淡水。

1.2 地下水源熱泵工程布置

區內布置有16眼抽水井和38眼回灌井,其中抽水井出水量120 m3/h,回灌井回灌量為50.53 m3/h,抽水和回灌均在第二層承壓含水層中進行。抽水井和回灌井呈交叉分布,具體布局如圖2所示。

圖2 研究區抽灌井分布Fig.2 Distribution of pumping and filling wells in the study area

該系統工程只取暖,不制冷。運行模式分為運行期、蓄熱期和恢復期,運行期時間為每年12月1日—次年2月28日,共3個月(90 d),所有抽灌井正常運行。蓄熱期為每年的6月1日—9月30日,共4個月(120 d)。為了保證在取熱運行期間被抽取熱源后盡快得到恢復,避免次年取熱受到影響,在系統停運后每年6—9月,與系統取熱運行期間一致,所有抽灌井正常運行,將抽取的地下水充分利用夏季氣溫高的特點進行儲存加熱,回灌水補熱溫度為25 ℃。其余時間段(3、4、5、10和11月)均為自然恢復期。

2 水—熱耦合數值模擬模型

2.1 水文地質概念模型

(1) 模擬平面范圍。在系統分析區域地下水滲流場的基礎上,參考能源站設計的抽灌井分布情況,以該區域為中心垂直和平行于水流方向各延伸1 km,以降低邊界對大規模開采引起的地下水流場變化的影響,構成2 km×2 km的矩形區域作為本次模擬范圍。

(2) 模擬垂向分層。根據前期勘查情況結合區域地質資料,本次模擬深度為地表至150 m深度范圍內。模型垂直方向上共分為5層：第1層為人工雜填土、粉質粘土和淤泥質粉質粘土,概化為潛水含水層,厚度為13.2 m；第2層為粉砂、細砂等,概化為淺層承壓含水層,層厚29.4 m；第3層為粉質粘土層,概化為第一弱透水層,厚度3.4 m；第4層為圓礫、中砂和細砂等,夾有薄層粉砂質粘土,整體概化為深層承壓含水層,厚度78.3 m；第5層為粘土層,概化為第二弱透水層,厚度25.7 m。根據地表高程數據及各層厚度,對模型各層的頂底板高程進行賦值,構建研究區三維結構模型。

(3) 邊界條件界定。區內地下水整體上由北向南徑流,由此平行于地下水流向的東西兩側邊界條件取為隔水邊界,垂直于水流方向的北側邊界為側向補給邊界,南部邊界為側向排泄邊界,兩者均取為流量邊界。頂部邊界為開放邊界,接受降雨入滲補給、農田灌溉水入滲補給、人工開采、蒸發排泄等,概化為流量邊界。根據江漢平原地下水均衡計算結果可知,地下水補給量與排泄量相對均衡差為2%左右；江漢平原區地表水豐富,地表水與地下水之間交互頻繁,降水量及其入滲補給量均比較大,區內地下水系統整體上保持均衡狀態。因此將研究區內第四系含水系統視為動態平衡系統,在模型中將入滲補給量等總補給量與地下水開采量和蒸發量等總排泄量視為均衡,而且模型的范圍比抽水回灌區大得多,還存在回灌補給,開采對邊界流量的影響可忽略,仍然保持原來的進出平衡,因此在模型中取為零流量邊界。模型底邊界為粉質粘土層,概化為隔水邊界。

研究區范圍較小,大地熱流的空間變異小,不同位置的地溫梯度變化小,地下水以水平運動為主,水流對地溫的影響較小；回灌情況下地下水開采對地下水流場的影響范圍較小,對邊界的影響可以忽略,因此可將側邊界取為定溫度邊界。上邊界為該區的恒溫帶,設置為定溫度邊界,近似取值為當地多年平均氣溫。同樣的方式,下邊界亦取為定溫度熱流邊界。

2.2 數學模型

2.2.1地下水流動數學模型

地下水—熱耦合模型涉及地下水流動和熱運移過程。根據滲流連續性方程和達西定律,地源熱泵系統地下水流動可用偏微分方程及其定解條件表示：

(1)

式中：Kx、Ky和Kz分別為x、y、z方向的滲透系數,m/d；H為含水層的水頭值,m；M為含水層厚度,m；ε為源匯項,m/d；S為給水度或比彈性釋水系數,潛水含水層取重力給水度,承壓含水層取彈性釋水系數；Ω為模擬范圍；n為邊界外法向方向單位向量；Γ為側邊界；B為底邊界；H0為初始水頭,m。

2.2.2熱量運移模型

地下熱水主要以液態的形式存在,在多孔介質中地下熱水熱量運移模式主要包括熱跟隨水流傳輸(對流)及熱量在水和多孔介質中的熱傳導,其方程可表示為[6]：

(2)

式中：n表示孔隙度；ρm表示含水介質的密度,kg/m3；ρw表示液體的密度,kg/m3；α表示縱向彌散度,m；cw表示液體的比熱容,J/(m3·K)；cm表示含水介質的比熱容,J/(m3·K)；λm表示含水介質熱傳導系數,J/(m·s·K)；qh表示熱量的輸出或輸入,J/s。

由于地下熱水系統中的熱量傳輸機制以對流、傳導為主,其數學表達式與地下水系統中溶質運移的數學模型一致。而且地下水溶質運移的數值模擬程序更為成熟,有學者使用地下水溶質運移模擬程序開展地下熱水的運移和熱傳輸模擬,模擬結果顯示MT3DMS的模擬結果與FEFLOW模擬結果相近。通用地下水數值模擬軟件(GMS v10.1)教程中也包含使用MT3DMS進行地下熱水流和熱傳輸模擬。

溶質運移的數學方程為：

(3)

式中：ρb表示含水介質的密度，kg/m3;Kd表示有效分布系數,m3/kg；Dm表示有效分子擴散系數,m2/d；va表示地下水流速,m/d；qCk表示源匯項,kg/(m3·s)。

從方程式(2)和(3)可以看出,MT3D模塊求解溶質運移與熱量運移方程的形式是一致的,MT3D模塊中溶質運移與熱量運移參數之間的關系如下式所示：

(4)

(5)

(6)

2.3 數值模型及參數設置

2.3.1離散化處理

在水文地質概念模型基礎上,利用MODFLOW數值模擬軟件的MT3D模塊建立數值模擬模型。本次模擬范圍為2 000 m×2 000 m×150 m。為了提高模型的計算速度,同時保證模型的計算精度,平面上采用有限差分法進行等間距剖分,網格間距為50 m。水平方向上,在抽水井和回灌井群區范圍內,對剖分網格進行加密處理,網格間距為10 m。模型垂直方向上共分為潛水含水層、淺層承壓含水層、第一弱透水層、深層承壓含水層和第二弱透水層等5層,共56 160個單元。

本次模擬以能源站熱源井剛開始運行的時間(2019年12月1日)作為初始時刻,以天為單位,對能源站地源熱泵系統運行20年后的地下水流場和溫度場特征進行預測。地下水三維水—熱模型計算過程中均采用非穩定流模擬,模型計算初始步長設置為1 d。

2.3.2邊界條件

(1) 水流邊界。根據研究區能源站鉆孔抽水試驗計算結果可知,抽水量為50 m3/h、80 m3/h和130 m3/h時,計算的影響半徑R分別為56、100和300 m。研究區側向邊界范圍遠遠大于抽水井影響半徑,因此將側向邊界作為隔水邊界進行處理,取零流量。模型底邊界為隔水邊界,取零流量。

(2) 溫度邊界。將模型第一層設置為恒溫帶,溫度取潛江市的多年日平均氣溫,取值18.3 ℃。

側邊界取為定溫度邊界,邊界溫度根據地溫梯度進行計算,賦入模型中。參考《中國大陸地區大地熱流數據匯編(第四版)》中江漢盆地的地溫梯度,取值范圍19.7～27.9 ℃/km,平均值為22.9 ℃/km。底邊界取定溫度邊界,溫度按照地溫梯度計算而來。

2.3.3初始條件

本次模擬主要根據研究區的地表高程數據,結合向中能源站鉆孔時監測的水位埋深,對模型中各剖分單元的初始水位進行賦值。研究區熱儲層(深層承壓含水層)地下水初始水位高程取值30 m。

依據地表多年日平均氣溫數據,結合研究區地溫梯度,按照各層的厚度計算其溫度,然后作為初始溫度分層賦入模型。研究區熱儲層(深層承壓含水層)初始地下水溫度取值19.5 ℃。

2.3.4模型參數的選用

構建三維水—熱耦合運移模型時,需要對水文地質參數和含水層骨架的熱物理性質參數進行賦值。根據前述水文地質概念模型設計,本次模擬需要對滲透系數(Kxx=Kyy、Kzz)、給水度、彈性釋水系數、孔隙度、流體熱傳導系數、固體熱傳導系數、流體比熱容和固體比熱容等參數進行賦值。根據野外試驗和室內熱物理參數試驗,結合相關文獻資料[7],初步確定了模型所需參數,不同層位(潛水含水層、淺層承壓含水層、第一弱透水層、深層承壓含水層和第二弱透水層)參數數值如表1所示。

表1 模擬區各層參數統計表Table 1 Statistical table of parameters of each layer in simulation area

利用MT3D模塊進行溶質運移和熱量運移模擬時,主要利用Advection、Dispersion和Chemical reaction程序包進行方程求解。其中Advection程序包采用軟件默認的參數進行,Dispersion和Chemical reaction程序包中各層參數賦值如表2所示。

表2 模擬區Dispersion和Chemical reaction程序包中各層參數表Table 2 Parameter table of each layer in the simulation area Dispersion and Chemical reaction package

2.4 抽灌井運行模擬結果分析

在驗證模型邊界條件、初始條件、參數分區和源匯項的基礎上,按區內地下水源熱泵工程實際運行模式,建立預測模型,模擬熱源井運行20年的地下水流場和溫度場的空間分布。

(1) 地下水流場特征。地源熱泵系統供暖期地下水流場分布如圖3所示,熱源井進行供暖期間在抽水井分布比較集中的區域形成小型水位降落漏斗,與初始地下水位(30 m)相比,此時水位落差接近4.5 m。實際抽水和回灌模式下熱源井進行完全回灌,水源熱泵系統運行1年和20年后的流場空間分布趨勢相同,并均與研究區初始地下水流場空間分布保持一致。

圖3 熱源井運行90天的地下水位等值線圖(圖中紅色圓點代表抽水井,黑色圓點代表回灌井)Fig.3 Contour map of groundwater level in 90 days operation of heat source well

從圖4可以看出,地源熱泵系統剛開始運行時,抽水井處水位迅速下降,此后水位保持較穩定狀態,停止運行后,水位迅速恢復至初始狀態,隨著夏季蓄熱期地熱井的運行,抽水井處水位同樣迅速下降后保持穩定,恢復期水位迅速恢復。回灌井處的水位隨著各個時期的變化呈現先升高、在保持穩定后下降的趨勢,如此循環反復(圖5)。

圖4 不同位置水位動態隨時間變化曲線(運行1年)Fig.4 Water level dynamic curve with time at different positions

圖5 不同位置水位動態隨時間變化曲線(運行20年)Fig.5 Water level dynamic curve with time at different positions

(2) 溫度場特征。地源熱泵系統供暖期的地下水溫度場分布如圖6所示。地源熱泵系統在運行期間,通過井群回灌向地下含水層注入冷水,此時注入的冷水堆積形成以回灌井為中心的冷水體,然后向回灌井周邊延伸,導致回灌井周邊溫度降低。系統運行期間距離回灌井10 m和20 m處地下水溫度受回灌冷水體的影響呈現迅速下降的趨勢,溫度分別降低至9 ℃和12 ℃；夏季蓄熱期(6—9月)回灌水體溫度較高(25 ℃),此時溫度呈現迅速上升的趨勢,溫度升高至24 ℃和22 ℃,恢復期溫度維持在較穩定狀態。距離回灌井40 m處和抽水井附近受回灌冷水體影響相對較小,溫度呈現緩慢下降趨勢,與初始溫度(19.5 ℃)相比,溫度下降2 ℃左右；夏季蓄熱期，地下水繼續受回灌冷水體影響呈現下降趨勢,之后恢復期受蓄熱期回灌水體的影響,溫度呈上升趨勢,并升高至18 ℃左右。

圖6 熱源井運行90天的溫度等值線圖Fig.6 Temperature contour map of heat source well in 90 days operation

抽水井和回灌井呈交叉分布,而且存在夏季蓄熱期(6—9月),此時地源熱泵系統運行1年和20年的溫度場如圖7和圖8所示。結合不同位置溫度隨時間變化動態曲線(圖9和圖10)可以看出,該模式下地源熱泵系統在回補(地下熱蓄)期末抽水井處溫度回升至接近20 ℃,回灌井處溫度穩定在24 ℃左右。由此可以看出,當抽水井和回灌井交叉分布,且存在蓄熱期時,抽水井和回灌井處溫度可達20～24 ℃,這種模式對提高系統能效和保證冬季供暖效果更為有利。

圖7 熱源井運行1年的溫度等值線圖Fig.7 Temperature contour map of heat source well in one year operation

圖8 熱源井運行20年的溫度等值線圖Fig.8 Temperature contour map of heat source well in operation for 20 years

圖9 不同位置溫度隨時間變化曲線(運行1年)Fig.9 Time varying curve of temperature at different positions(1 year operation)

圖10 不同位置溫度隨時間變化曲線(運行20年)Fig.10 Time dependence curve of temperature at different positions(20 yeasr operation)

3 結論與建議

(1) 本文針對研究區實際布置運行的地下水源熱泵系統工程進行地下水—熱耦合模擬,結果表明,在抽灌井交叉布局且存在蓄熱期回補熱量的情況下,地下水源熱泵運行期結束后抽灌井處地下水位、溫度能較快地得到恢復,次年—若干年后取水、取熱不受影響,地下水源熱泵系統工程能持續有效地運行。

(2) 模擬過程中采用概化地層、綜合參數進行取值,由于地層(含水層)滲透性的差異,實際抽水回灌過程中可能出現強滲透性層(圓礫層)流量過大導致的滲流、換熱集中的現象,從而帶來偏差。具體工程實施中,應盡量采取措施調節換熱場區含水層段流量及熱量交換,以維持充分換熱及熱平衡。