徐連三，華 杉，劉紅衛(wèi)，柯 立，李永輝，王富強(qiáng)

(1.湖北省地質(zhì)局 武漢水文地質(zhì)工程地質(zhì)大隊(duì),湖北 武漢 430051； 2.湖北省地?zé)崮苎芯客茝V中心,湖北 武漢 430051)

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

1 研究區(qū)概況

1.1 研究區(qū)水文地質(zhì)條件

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

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

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

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

1.2 地下水源熱泵工程布置

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

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

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

2 水—熱耦合數(shù)值模擬模型

2.1 水文地質(zhì)概念模型

(1) 模擬平面范圍。在系統(tǒng)分析區(qū)域地下水滲流場的基礎(chǔ)上,參考能源站設(shè)計(jì)的抽灌井分布情況,以該區(qū)域?yàn)橹行拇怪焙推叫杏谒鞣较蚋餮由? km,以降低邊界對大規(guī)模開采引起的地下水流場變化的影響,構(gòu)成2 km×2 km的矩形區(qū)域作為本次模擬范圍。

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

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

研究區(qū)范圍較小,大地?zé)崃鞯目臻g變異小,不同位置的地溫梯度變化小,地下水以水平運(yùn)動為主,水流對地溫的影響較??；回灌情況下地下水開采對地下水流場的影響范圍較小,對邊界的影響可以忽略,因此可將側(cè)邊界取為定溫度邊界。上邊界為該區(qū)的恒溫帶,設(shè)置為定溫度邊界,近似取值為當(dāng)?shù)囟嗄昶骄鶜鉁?。同樣的方?下邊界亦取為定溫度熱流邊界。

2.2 數(shù)學(xué)模型

2.2.1地下水流動數(shù)學(xué)模型

地下水—熱耦合模型涉及地下水流動和熱運(yùn)移過程。根據(jù)滲流連續(xù)性方程和達(dá)西定律,地源熱泵系統(tǒng)地下水流動可用偏微分方程及其定解條件表示：

(1)

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

2.2.2熱量運(yùn)移模型

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

(2)

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

由于地下熱水系統(tǒng)中的熱量傳輸機(jī)制以對流、傳導(dǎo)為主,其數(shù)學(xué)表達(dá)式與地下水系統(tǒng)中溶質(zhì)運(yùn)移的數(shù)學(xué)模型一致。而且地下水溶質(zhì)運(yùn)移的數(shù)值模擬程序更為成熟,有學(xué)者使用地下水溶質(zhì)運(yùn)移模擬程序開展地下熱水的運(yùn)移和熱傳輸模擬,模擬結(jié)果顯示MT3DMS的模擬結(jié)果與FEFLOW模擬結(jié)果相近。通用地下水?dāng)?shù)值模擬軟件(GMS v10.1)教程中也包含使用MT3DMS進(jìn)行地下熱水流和熱傳輸模擬。

溶質(zhì)運(yùn)移的數(shù)學(xué)方程為：

(3)

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

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

(4)

(5)

(6)

2.3 數(shù)值模型及參數(shù)設(shè)置

2.3.1離散化處理

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

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

2.3.2邊界條件

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

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

側(cè)邊界取為定溫度邊界,邊界溫度根據(jù)地溫梯度進(jìn)行計(jì)算,賦入模型中。參考《中國大陸地區(qū)大地?zé)崃鲾?shù)據(jù)匯編(第四版)》中江漢盆地的地溫梯度,取值范圍19.7～27.9 ℃/km,平均值為22.9 ℃/km。底邊界取定溫度邊界,溫度按照地溫梯度計(jì)算而來。

2.3.3初始條件

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

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

2.3.4模型參數(shù)的選用

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

表1 模擬區(qū)各層參數(shù)統(tǒng)計(jì)表Table 1 Statistical table of parameters of each layer in simulation area

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

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

2.4 抽灌井運(yùn)行模擬結(jié)果分析

在驗(yàn)證模型邊界條件、初始條件、參數(shù)分區(qū)和源匯項(xiàng)的基礎(chǔ)上,按區(qū)內(nèi)地下水源熱泵工程實(shí)際運(yùn)行模式,建立預(yù)測模型,模擬熱源井運(yùn)行20年的地下水流場和溫度場的空間分布。

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

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

從圖4可以看出,地源熱泵系統(tǒng)剛開始運(yùn)行時(shí),抽水井處水位迅速下降,此后水位保持較穩(wěn)定狀態(tài),停止運(yùn)行后,水位迅速恢復(fù)至初始狀態(tài),隨著夏季蓄熱期地?zé)峋倪\(yùn)行,抽水井處水位同樣迅速下降后保持穩(wěn)定,恢復(fù)期水位迅速恢復(fù)。回灌井處的水位隨著各個(gè)時(shí)期的變化呈現(xiàn)先升高、在保持穩(wěn)定后下降的趨勢,如此循環(huán)反復(fù)(圖5)。

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

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

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

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

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

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

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

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

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

3 結(jié)論與建議

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

(2) 模擬過程中采用概化地層、綜合參數(shù)進(jìn)行取值,由于地層(含水層)滲透性的差異,實(shí)際抽水回灌過程中可能出現(xiàn)強(qiáng)滲透性層(圓礫層)流量過大導(dǎo)致的滲流、換熱集中的現(xiàn)象,從而帶來偏差。具體工程實(shí)施中,應(yīng)盡量采取措施調(diào)節(jié)換熱場區(qū)含水層段流量及熱量交換,以維持充分換熱及熱平衡。