董 超，杜寧宇，駱祖江*

(1.河海大學(xué)地球科學(xué)與工程學(xué)院，江蘇南京 210098；2.河海大學(xué)河海里爾學(xué)院，江蘇南京 210098)

0 引言

地?zé)豳Y源蘊(yùn)藏豐富且無(wú)污染，是一種可再生能源[1-2]。隨著石油、煤炭等傳統(tǒng)能源逐漸枯竭，合理開發(fā)地?zé)豳Y源[3]，對(duì)緩解能源束縛和環(huán)境壓力，提高經(jīng)濟(jì)增長(zhǎng)的質(zhì)量具有重要的意義，也是實(shí)現(xiàn)大力開發(fā)新能源，實(shí)現(xiàn)能源清潔發(fā)展、安全發(fā)展、環(huán)境友好發(fā)展和可持續(xù)發(fā)展的一項(xiàng)重要的戰(zhàn)略舉措[4]。江蘇處于中國(guó)東部沿海高熱流地?zé)岙惓?，大地?zé)崃鲝?qiáng)，地溫梯度高，熱儲(chǔ)層分布廣、厚度大、埋藏條件好，具備形成中溫地?zé)豳Y源的地質(zhì)條件[5-7]。以往的地?zé)崴Y源評(píng)價(jià)，往往簡(jiǎn)化為只評(píng)價(jià)水量，將水溫看作恒定不變，或?qū)⑺亢退疁赝耆珠_來(lái)評(píng)價(jià)，和實(shí)際不相吻合。本文以江蘇省江淮生態(tài)經(jīng)濟(jì)區(qū)東南部地區(qū)地?zé)崴到y(tǒng)為例，根據(jù)研究區(qū)地質(zhì)、水文地質(zhì)及地?zé)岬刭|(zhì)特征，在概化出地?zé)崴到y(tǒng)概念模型的基礎(chǔ)上，建立了研究區(qū)地?zé)崴欠€(wěn)定滲流與熱量運(yùn)移三維耦合數(shù)值模型[8]，同時(shí)評(píng)價(jià)了現(xiàn)有地?zé)崴?0m降深條件下的可開采資源量及溫度的變化[9]，為該地區(qū)地?zé)崴目沙掷m(xù)開發(fā)利用提供了科學(xué)依據(jù)，具有較強(qiáng)的工程實(shí)際意義[10]。

1 地質(zhì)及水文地質(zhì)條件

研究區(qū)位于江淮生態(tài)經(jīng)濟(jì)區(qū)東南部的寶應(yīng)、建湖、高郵、興化4個(gè)縣(市)。地形以平原為主，總體地勢(shì)西高東低。東、南側(cè)以已存在的江淮生態(tài)經(jīng)濟(jì)區(qū)的范圍為界，西側(cè)以高郵湖、白馬湖等河流湖泊為界，北側(cè)為根據(jù)行政區(qū)域及自然邊界人為劃定的邊界，總面積約為6 900km2。

研究區(qū)地處蘇北盆地中部，地層隸屬下?lián)P子地層分區(qū)，地表均為第四系覆蓋，地層分布不全，3 000m以淺分布的地層呈現(xiàn)一老一新特點(diǎn)。即下部為震旦系—下古生界，構(gòu)成蘇北盆地基底，上部為中新生代碎屑巖，分布穩(wěn)定。

研究區(qū)所處的蘇北盆地處位于中國(guó)東部沿海高熱流地?zé)岙惓^(qū)，大地?zé)崃鲝?qiáng)，地溫梯度高，具備形成中高溫地?zé)豳Y源的地質(zhì)條件，地?zé)豳Y源成因主要為傳導(dǎo)型和中低溫對(duì)流型[11]。區(qū)內(nèi)地層深部分布有厚度巨大的碳酸鹽巖，這些地層在地質(zhì)歷史時(shí)期經(jīng)歷了強(qiáng)烈的構(gòu)造變動(dòng)和巖溶化作用，巖溶裂隙發(fā)育，為熱水型地?zé)豳Y源的賦存提供了良好的儲(chǔ)存空間。區(qū)內(nèi)斷裂構(gòu)造十分發(fā)育，但規(guī)模深淺不一，在凸起的邊部、凹陷外斜坡帶、凹陷邊緣斷階帶、隆起區(qū)的邊緣地帶與凹陷邊緣斷階帶均受區(qū)域性深斷裂控制，具有極好的水熱循環(huán)條件，也是地?zé)豳Y源富集區(qū)，主要發(fā)育NE、NNE、EW及NW向四組斷裂，NNE和NE向斷裂為大、中型新構(gòu)造斷裂構(gòu)造帶，與后期NW向張性、張扭性斷裂交匯，構(gòu)成了深部熱流和地下水循環(huán)的良好通道[12-13]。

區(qū)內(nèi)地下水根據(jù)含水介質(zhì)類型劃分為松散巖類孔隙水、碳酸鹽巖類巖溶裂隙水及碎屑巖類裂隙水。松散巖類孔隙水在區(qū)內(nèi)廣泛發(fā)育，主要賦存在鹽城組和第四紀(jì)松散堆積地層中的砂層及砂礫層中；碳酸鹽巖類巖溶裂隙水賦存于震旦系和古生界碳酸鹽巖之中；碎屑巖類裂隙水主要賦存在古近系基巖中。本次對(duì)區(qū)內(nèi)6口地?zé)峋M(jìn)行了調(diào)查(表1)。

表1 研究區(qū)地?zé)崴艣r

2 地?zé)崴到y(tǒng)概念模型

研究區(qū)熱水系統(tǒng)頂部一方面接受大氣降雨的補(bǔ)給，是一補(bǔ)給邊界，另一方面地下水又通過(guò)其蒸發(fā)，是一排泄邊界。地?zé)崴ㄟ^(guò)斷裂、裂隙的連通進(jìn)行循環(huán)，水流呈非穩(wěn)定流狀態(tài)，區(qū)內(nèi)斷裂構(gòu)造十分發(fā)育，規(guī)模深淺不一，各含水層均為非均質(zhì)各向異性，各層之間均發(fā)生水力聯(lián)系，地下水位受大氣降水和蒸發(fā)季節(jié)性變化的影響，為三維非穩(wěn)定流。由于研究區(qū)范圍較大，邊界距離地?zé)崴^遠(yuǎn)，因此模型四周均概化為第一類定水頭邊界，模型底部為隔水邊界。

研究區(qū)熱水系統(tǒng)頂部接受大氣降雨的補(bǔ)給，混入冷水，是一熱量排泄邊界；底部通過(guò)斷裂、裂隙的聯(lián)通、循環(huán)，將深層的較高溫度帶到淺層熱儲(chǔ)含水層，是一熱量補(bǔ)給邊界，根據(jù)增溫梯度推算設(shè)定為一恒溫邊界，故頂?shù)撞扛呕癁榈谝活悷崃窟吔?。?jì)算范圍內(nèi)設(shè)定地下水和含水介質(zhì)骨架的熱動(dòng)平衡是瞬時(shí)完成的，將熱儲(chǔ)含水層四周概化為第二類熱量邊界。

3 地?zé)崴欠€(wěn)定滲流與熱量運(yùn)移三維耦合數(shù)值模型

根據(jù)上述概化出的地?zé)崴到y(tǒng)概念模型，建立相應(yīng)的地?zé)崴欠€(wěn)定滲流與熱量運(yùn)移三維耦合數(shù)值模型[14]。地?zé)崴疂B流場(chǎng)與溫度場(chǎng)通過(guò)地?zé)崴鲃?dòng)方程進(jìn)行耦合計(jì)算。

3.1 地?zé)崴欠€(wěn)定滲流數(shù)值模型

取坐標(biāo)軸方向與熱儲(chǔ)含水層各向異性主滲透系數(shù)方向一致，建立如下數(shù)值模型：

(1)

式中：Kxx、Kyy、Kzz為各向異性主方向滲透系數(shù)，m/d；h為點(diǎn)(x,y,z)在t時(shí)刻的水頭值，m；W為源匯項(xiàng)(1/d)；t為時(shí)間，d；Ω為計(jì)算域；h0(x,y,z,t0)為點(diǎn)(x,y,z)處初始水頭值，m；q(x,y,z,t)為第二類邊界上單位面積的補(bǔ)給量，m/d；cos(n,x)、cos(n,y)、cos(n,z)為流量邊界外法線方向與坐標(biāo)軸方向夾角的余弦；μ為飽和差(自由面上升)或給水度(自由面下降)；Ss為儲(chǔ)水率，L/m；Γ1為第一類邊界；Γ2為第二類邊界；Γ3為自由面邊界。

3.2 地?zé)崴疅崃窟\(yùn)移數(shù)值模型

假設(shè)地?zé)崴秃橘|(zhì)骨架的熱動(dòng)平衡是瞬時(shí)完成的，并忽略由于溫度差引起水的密度不一樣而引起的上下自然對(duì)流的影響，建立地?zé)崴S熱量運(yùn)移數(shù)值模型如下：

(2)

式中：λx、λy、λz為各方向水的熱動(dòng)力彌散系數(shù)；cw為水的熱容量，J/(m3·K)；c為含水介質(zhì)的熱容量，J/(m3·K)；vx、vy、vz為地?zé)崴疂B流速度分量，m/d；T0(x,y,z)為點(diǎn)(x,y,z)處初始溫度值；T1(x,y,z,t)為第一類邊界的溫度函數(shù)；Γ1為第一類邊界；Qc為熱源匯項(xiàng)，Qc=cwW(TQ-T)；TQ為源匯項(xiàng)的溫度；n為邊界外法線向量；Q(x,y,z,t)為第二類邊界上已知的熱量或熱流函數(shù)，J/(m·d)；Γ2為第二類邊界。

3.3 地?zé)崴欠€(wěn)定滲流與熱量運(yùn)移三維耦合數(shù)值模型

地?zé)崴鬟\(yùn)動(dòng)方程：

(3)

將(1)式與(2)式通過(guò)地?zé)崴鬟\(yùn)動(dòng)方程式(3)耦合在一起，構(gòu)成研究區(qū)地?zé)崴欠€(wěn)定滲流與熱量運(yùn)移三維耦合數(shù)值模型。

4 模型求解

上述模型采用伽遼金有限元法求解，首先將地下水滲流與熱量運(yùn)移的性代數(shù)方程組代入求解方程，其插值稱為誤差函數(shù)或剩余，從某種平均意義上講希望該誤差為零。讓微分方程的近似解和精確解之差的總剩余量最小，用近似解來(lái)代替精確解。

求解滲流問題的代數(shù)方程形式：

(4)

求解熱量運(yùn)移問題的代數(shù)方程形式為：

(5)

式中：[K]為總體滲透矩陣；{h}為未知節(jié)點(diǎn)水頭列陣；[P]為總體儲(chǔ)水矩陣；[U]為總體自由面流量補(bǔ)給矩陣；{F}為已知的右端項(xiàng)；[G]為總體彌散矩陣；{T}為未知節(jié)點(diǎn)溫度列陣；[R]為總體對(duì)流矩陣；[M]為總體熱容量矩陣；{Z}為已知的右端項(xiàng)。

再對(duì)式中的時(shí)間采用隱式差分格式進(jìn)行離散，整理后可以得到：

(6)

(7)

利用上述方程，根據(jù)初始條件下的滲流場(chǎng)水頭和溫度分布，即可計(jì)算任意時(shí)間的滲流場(chǎng)水頭和溫度分布。

5 模型識(shí)別驗(yàn)證

根據(jù)研究區(qū)地下水系統(tǒng)結(jié)構(gòu)特征及幾何形狀，對(duì)研究區(qū)進(jìn)行三維剖分。平面上對(duì)地?zé)崴車椭饕膶?dǎo)水控水?dāng)嗔堰M(jìn)行網(wǎng)格加密，最終形成每層36 082個(gè)三角形單元和18 482個(gè)節(jié)點(diǎn)。垂向上考慮到地?zé)崴鏊畬游患皩雍竦挠绊?，由上往下共剖?個(gè)層位，9個(gè)計(jì)算層面，共計(jì)288 656個(gè)單元，166 338個(gè)節(jié)點(diǎn)。計(jì)算域空間單元剖分見圖1。

利用現(xiàn)有興熱1井、興熱2井、RGM1井、RGS1井、RBQ1井、寶熱1井的抽水試驗(yàn)資料對(duì)模型進(jìn)行識(shí)別、驗(yàn)證，并根據(jù)研究區(qū)各地?zé)崴跏妓?、溫度? 500m地溫分布圖及實(shí)測(cè)地溫增溫規(guī)律，確定模型初始流場(chǎng)和初始溫度場(chǎng)。圖2舉例說(shuō)明了模型第一層地?zé)崴跏剂鲌?chǎng)和初始溫度場(chǎng)。各開采井作為觀測(cè)井同時(shí)進(jìn)行水位和水溫?cái)M合。

圖1 計(jì)算域空間剖分Figure 1 Model space subdivision of computational domain

圖2 模型第一層初始流場(chǎng)及初始溫度場(chǎng)分布Figure 2 Distribution of initial flow field and initial temperature field in the first layer of modela.初始流場(chǎng)圖(m)；b.初始溫度場(chǎng)圖(℃)

全區(qū)共分108個(gè)參數(shù)分區(qū)，上部三層根據(jù)地?zé)崴植己偷貙訋r性特征進(jìn)行分區(qū)，下部地層受斷裂帶影響，根據(jù)斷裂帶劃分為斷裂影響帶和無(wú)影響區(qū)。通過(guò)反演計(jì)算，獲得了各層各參數(shù)分區(qū)的參數(shù)值，圖3舉例說(shuō)明了第八層的水文地質(zhì)參數(shù)分區(qū)，表2舉例說(shuō)明了各參數(shù)分區(qū)的參數(shù)值。各地?zé)崴樗囼?yàn)實(shí)測(cè)數(shù)據(jù)與計(jì)算數(shù)據(jù)擬合情況見圖4。從擬合結(jié)果來(lái)看，計(jì)算曲線與實(shí)測(cè)曲線擬合精度較好，總體變化趨勢(shì)一致。通對(duì)比對(duì)可以發(fā)現(xiàn)，模型計(jì)算滿足有關(guān)精度要求[15]，可用于研究區(qū)地下水滲流場(chǎng)與溫度場(chǎng)的模擬預(yù)測(cè)。

圖3 第Ⅷ層參數(shù)分區(qū)Figure 3 Model parameters partition of the eighth group

表2 第Ⅷ層各參數(shù)分區(qū)參數(shù)

續(xù)表

6 地?zé)崴砷_采資源量規(guī)劃評(píng)價(jià)

利用識(shí)別、驗(yàn)證后的地?zé)崴欠€(wěn)定滲流與熱量運(yùn)移三維耦合數(shù)值模型，對(duì)研究區(qū)現(xiàn)有6口地?zé)崴诘臒醿?chǔ)含水層2022-03-01至2032-03-01的地下水流場(chǎng)與溫度場(chǎng)變化進(jìn)行預(yù)測(cè)，每一年分為12個(gè)應(yīng)力期，共計(jì)120個(gè)應(yīng)力期。由于研究區(qū)地下水開采量少，水位、溫度變化幅度小，初始流場(chǎng)和初始溫度場(chǎng)均采用模型識(shí)別、驗(yàn)證階段的初始數(shù)據(jù)，大氣降雨量數(shù)據(jù)采用2000—2020年逐月平均值，氣溫采用年平均氣溫15℃。經(jīng)模型模擬預(yù)測(cè)，各開采井在2022年3月1日至2032年3月1日未來(lái)10a內(nèi)的水位、溫度變化歷時(shí)曲線見圖5。

圖4 研究區(qū)水位觀測(cè)值與計(jì)算值擬合Figure 4 Fitting results of the observed and calculated water level values in the study area

圖5 運(yùn)行10年水位及溫度變化歷時(shí)曲線Figure 5 Diachronic curve of water level and temperature changes during 10 years of operation

從模型運(yùn)行得出的結(jié)果可以看出，按照50m水位降深確定的開采量進(jìn)行開采，隨著開采的持續(xù)，整個(gè)研究區(qū)水位逐漸下降10a后達(dá)到平衡，地下水位處于相對(duì)穩(wěn)定狀態(tài)；各地?zé)崴疁囟茸兓怀^(guò)1℃。其中興熱1井、興熱2井、RGM1井、RGS1井和寶熱1井整體溫度場(chǎng)相對(duì)于初始狀態(tài)有所上升，RBQ1井整體溫度場(chǎng)相對(duì)于初始狀態(tài)有所下降。由于地?zé)崴Y源的持續(xù)開采，溫度在局部范圍內(nèi)將發(fā)生一定的變化，將會(huì)產(chǎn)生冷熱堆積，對(duì)周圍水環(huán)境產(chǎn)生一定的影響。

根據(jù)模擬計(jì)算結(jié)果，整個(gè)研究區(qū)的6口地?zé)崴凑詹怀^(guò)50m降深開采地?zé)崴鞯責(zé)崴砷_采資源量及水位變化統(tǒng)計(jì)見表3，合計(jì)6口地?zé)崴砷_采資源量為5 1290 m3/d。

表3 研究區(qū)各地?zé)崴A(yù)測(cè)水位變化及取水量

7 結(jié)論

1)通過(guò)建立江淮地區(qū)地?zé)崴欠€(wěn)定滲流與熱量運(yùn)移三維耦合數(shù)值模型，可以同時(shí)預(yù)測(cè)未來(lái)開采條件下水位和水溫的變化趨勢(shì)，并可根據(jù)水位和水溫兩個(gè)指標(biāo)的變化，規(guī)劃評(píng)價(jià)地?zé)崴目砷_采資源量，實(shí)現(xiàn)地?zé)崴目沙掷m(xù)開采利用。

2)江淮地區(qū)現(xiàn)有各地?zé)崴?0a內(nèi)水位降深穩(wěn)定在50m的條件下，各地?zé)崴責(zé)崴目砷_采資源量分別為興熱1井3 870 m3/d、興熱2井7 630 m3/d、RGM1井3 720 m3/d、RGS1井890 m3/d、RBQ1井1 480 m3/d、寶熱1井33 700 m3/d。全區(qū)可開采資源量為51 290 m3/d，各井的溫度變化均不超過(guò)1℃。