淺層地?zé)崮懿煌_發(fā)方案對(duì)深層地下淡水咸化效應(yīng)研究

李 兆，駱祖江，楊 璐，胡 穎

（1．河海大學(xué) 地球科學(xué)與工程學(xué)院，江蘇 南京 211100；2．江蘇省地質(zhì)工程勘察院，江蘇 南京 211102）

1 研究背景

淡水是人類生產(chǎn)、生活不可缺少的資源，我國(guó)655個(gè)城市中超過400個(gè)以地下水作為飲用水源，然而大量研究顯示濱海地區(qū)地下淡水資源有咸化的趨勢(shì)［1-2］，地下淡水資源正在受到威脅，為了進(jìn)一步保護(hù)地下淡水資源，咸水入侵防控成為了當(dāng)今地下水科學(xué)研究的熱點(diǎn)。

研究認(rèn)為地下水開采導(dǎo)致地下水位下降，咸、淡水分界面向內(nèi)陸侵襲是導(dǎo)致地下水咸化的主要原因［3-4］，在濱海地區(qū)開展諸多關(guān)于海水入侵的研究［5-7］。然而很多濱海地區(qū)地下水受沉積環(huán)境影響，天然狀態(tài)下水質(zhì)較差，特別是在孔隙地下水常呈現(xiàn) “上咸下淡”的格局，即淺部含水層貯存咸水，深部含水層貯存淡水［8］，周慧芳等［9］就南通地區(qū)地下水咸化進(jìn)行了系統(tǒng)研究，提出淡水咸化原因有海水入侵、水—巖作用和淺部咸水越流入侵。Ma等［10］針對(duì)淺部咸水越流入侵深部淡水這一問題，提出適量開采淺層地下咸水，減小與深層淡水之間的水力梯度，以緩解深部淡水咸化，為咸水入侵防控提供新的思路，但是開采淺部咸水如何利用成為了新的問題。

近年來我國(guó)大力發(fā)展淺層地?zé)崮芤詫?shí)現(xiàn)清潔取暖、節(jié)能減排，淺層地?zé)崮艿拈_發(fā)利用可以有效解決淺層咸水的利用問題。地下水源熱泵系統(tǒng)因其造價(jià)低、換熱功率高而被廣泛使用［11］，過去針對(duì)地下水源熱泵系統(tǒng)的研究大多是針對(duì)熱量運(yùn)移［12］，而關(guān)于地下水源熱泵系統(tǒng)運(yùn)行對(duì)地下水水質(zhì)影響的研究較少，De Keuleneer等［13］曾就濱海地區(qū)利用地下水源熱泵系統(tǒng)控制海水入侵，并就不同抽、灌井布井方式進(jìn)行研究，但是沒探討地下水源熱泵系統(tǒng)應(yīng)用對(duì)淺部咸水入侵深部淡水的影響。

本次研究基于地下水滲流、熱量運(yùn)移和溶質(zhì)運(yùn)移理論，以江蘇省地質(zhì)工程勘察院（南通分院）科研實(shí)驗(yàn)樓地下水源熱泵系統(tǒng)為研究對(duì)象，建立地下水滲流—熱量運(yùn)移—溶質(zhì)運(yùn)移多場(chǎng)耦合數(shù)值模型，在考慮合理科學(xué)開發(fā)淺層地?zé)崮艿幕A(chǔ)上，進(jìn)一步評(píng)價(jià)淺層地?zé)崮荛_發(fā)對(duì)深層地下淡水咸化控制效應(yīng)。

2 工程概況

科研實(shí)驗(yàn)樓位于南通市崇川區(qū)，分為主樓（6F）和附樓（2F），采用地下水源熱泵系統(tǒng)進(jìn)行制冷和供暖，根據(jù)設(shè)計(jì)指標(biāo)，科研實(shí)驗(yàn)樓冬季最大負(fù)荷345 kW，夏季最大負(fù)荷465 k W。地下水源熱泵系統(tǒng)工程布置了1口抽水井和1回灌井，受場(chǎng)地限制，抽水井和回灌井之間的距離為58 m。場(chǎng)地布置示意圖見圖1。

圖1 江蘇省地質(zhì)工程勘查院（南通分院）地下水源熱泵場(chǎng)地布置示意

2．1 沉積地層 南通市第四系厚度一般可達(dá)300 m，其水文地質(zhì)體結(jié)構(gòu)是砂土與黏性土互層組成的復(fù)合含水層系統(tǒng)。區(qū)內(nèi)下更新統(tǒng)（Q1）埋深在190 m以下，厚度60 m左右。由中粗砂、粉細(xì)砂構(gòu)成，為河湖相沉積，含水砂層構(gòu)成區(qū)第Ⅲ承壓含水層組；江蘇省地質(zhì)工程勘察院（南通分院）一帶缺失中更新統(tǒng)（Q2）砂層，即研究區(qū)缺失第Ⅱ承壓含水層；區(qū)內(nèi)上更新統(tǒng)（Q3）埋深在40～120 m之間，厚度80～100 m。由中粗砂、粉細(xì)砂構(gòu)成，受兩次海浸影響，形成海陸交互相沉積，含水砂層構(gòu)成區(qū)內(nèi)第Ⅰ承壓含水層組；區(qū)內(nèi)全新統(tǒng)（Q4）厚度40～50 m，由粉細(xì)砂、亞砂土構(gòu)成，含淤泥質(zhì)，為三角洲海陸交互相沉積。上部多為長(zhǎng)江沖積土層，屬近現(xiàn)代沉積物［14］。

2．2 水文地質(zhì)條件 工程場(chǎng)地的地下水以松散孔隙水為主，主要賦存在砂層中，根據(jù)沉積環(huán)境可劃分為潛水、第Ⅰ、第Ⅲ承壓含水層，各含水層的富水性好，但水質(zhì)差異較大。含水層之間夾黏土層，分別第Ⅰ、第Ⅲ黏性土弱透水層，水文地質(zhì)剖面圖見圖2。潛水主要賦存于全新統(tǒng)地層（Q4）中，含水層下段具有一定微承壓特性，水位埋深1．5～2．0 m，水質(zhì)微咸，TDS在800～1600 mg/L之間，地下水化學(xué)類型為HCO3·Cl-Ca·Mg和HCO3·Cl·SO4-Ca·Na，潛水與地表水有著密切的水力聯(lián)系。第Ⅰ承壓水主要賦存于上更新統(tǒng)（Q3）中，水位埋深在3 m左右，水質(zhì)差，TDS在1100～1800 mg/L之間，水化學(xué)類型為Cl·HCO3-Ca·Na。第Ⅲ承壓水主要賦存于下更新統(tǒng)（Q1）中，水質(zhì)好，TDS在700 mg/L左右，水化學(xué)類型為HCO3·Cl-Na·Ca。

3 數(shù)值模型

3．1 概念模型 以工程場(chǎng)地為中心向四周各擴(kuò)展1000 m作為計(jì)算邊界，模型底部以第Ⅲ承壓含水層底部作為底部邊界。根據(jù)土體的性質(zhì)，模型垂向上共剖分5層，分別為潛水含水層、第Ⅰ弱透水層、第Ⅰ承壓含水層、第Ⅲ弱透水層和第Ⅲ承壓含水層。由于計(jì)算邊界與抽、灌井距離遠(yuǎn)，結(jié)合南通地區(qū)水文地質(zhì)條件、滲透系數(shù)參考值和影響半徑經(jīng)驗(yàn)公式，地下水源熱泵系統(tǒng)的運(yùn)行過程中抽、灌地下水影響范圍小于1000 m，地下水源熱泵系統(tǒng)運(yùn)行對(duì)數(shù)值模型邊界影響很小。因此，四周邊界處理為Dirichlet邊界，即水力學(xué)上處理為定水頭邊界、熱力學(xué)上處理為恒溫邊界、水化學(xué)上處理為恒定濃度邊界；根據(jù)計(jì)算域水力學(xué)特征，模型頂面處理為自由面邊界，底部處理為Neumann邊界，即水力學(xué)上處理為隔水邊界、熱力學(xué)上處理為熱通量邊界、水化學(xué)上處理為零通量邊界。

本項(xiàng)目利用第Ⅰ承壓含水層的淺層地?zé)崮?，需要抽取、回灌第Ⅰ承壓含水層的地下水，由于各層存在水位差，各含水層之間存在水力聯(lián)系，地下水流特征概化為三維非穩(wěn)定流。地下水源熱泵系統(tǒng)運(yùn)行過程中，假定每個(gè)節(jié)點(diǎn)的溫度同時(shí)代表地下水和土顆粒溫度，忽略地下水密度隨溫度發(fā)生的變化。3．2 數(shù)學(xué)模型 各含水層之間的黏性土層作為弱透水層和含水層一起直接參與計(jì)算，地下水流為三維非穩(wěn)定流。取坐標(biāo)軸方向與各向異性介質(zhì)主方向一致，地下水系統(tǒng)滲流控制方程如下［15］：

式中：Kxx、Kyy、Kzz為各方向水的滲透系數(shù)，m/d；h為地下水位，m；W為源匯項(xiàng)，1/d；Ss為含水層貯水率，1/m。

假設(shè)地下水和含水介質(zhì)骨架的熱動(dòng)平衡是瞬時(shí)完成的，即含水介質(zhì)骨架與周圍地下水具有相同的溫度，并忽略由于溫度差引起水的密度不一樣而引起的上下自然對(duì)流的影響。地下水系統(tǒng)熱量運(yùn)移控制方程如下［16-18］：

式中：λx、λy、λz為各方向水的熱動(dòng)力彌散系數(shù)，W/（m·℃）；T為巖土體溫度，℃；cw、cs分別為水和土體的體積熱容量，J/（m3·℃）；u為多孔介質(zhì)中地下水平均流速，m/d；Qc為熱的源匯項(xiàng)，J/（m3·d）；βL、βT分別為縱向和橫向熱彌散度，m；δij為Kronecker記號(hào)；ui、uj為地下水平均流速分量，m/d；λ′ij、λp分別為熱機(jī)械彌散系數(shù)和熱傳導(dǎo)系數(shù)，W/（m·℃）。

根據(jù)多孔介質(zhì)流體動(dòng)力學(xué)，溶質(zhì)遷移的過程中滿足質(zhì)量守恒定律，本次研究針對(duì)稀物質(zhì)在水中遷移過程，可假定流體的密度為定值。取坐標(biāo)軸方向與主滲透方向一致，則地下水溶質(zhì)運(yùn)移控制方程為［19-21］：

式中：Dxx、Dyy、Dzz為各方向的水動(dòng)力彌散系數(shù)分量，m2/d；c為地下水中溶質(zhì)的質(zhì)量濃度，mg/L；n為孔隙度；I為溶質(zhì)源匯項(xiàng)，mg/（L·d）；αL、αT分別為縱向和橫向溶質(zhì)彌散度，m；D′ij為機(jī)械彌散系數(shù)，m2/d；Dd為分子有效擴(kuò)散系數(shù)，m2/d。

研究表明：當(dāng)?shù)叵滤魉佥^快時(shí)，機(jī)械彌散系數(shù)遠(yuǎn)大于分子有效擴(kuò)散系數(shù)，因此忽略分子擴(kuò)散系數(shù)的影響［22］。

地下水平均流速分量分別表示為：

根據(jù)概念模型，建立與之相應(yīng)的初始條件和邊界條件，所建立的數(shù)值模型如圖3—4所示。

圖3 模型剖分平面示意圖

圖4 三維數(shù)值模型示意圖

3．3 模型識(shí)別、驗(yàn)證 以2017年6月作為模型計(jì)算的初始時(shí)刻，由于模擬區(qū)相對(duì)較小，同層初始水位、地下水化學(xué)濃度、地溫梯度變化較小，模型潛水、第Ⅰ、第Ⅲ承壓含水層的初始條件根據(jù)實(shí)測(cè)水位、地下水化學(xué)濃度、巖土體溫度進(jìn)行賦值，相鄰黏性土弱透水層初始條件通過插值獲取，模型各層地下水位初始值從上往下依次為0．87，0．8，0．73，-9．14及-19．01 m；Cl-濃度初始值從上往下依次為2570，2100，1630，879及128 mg/L；平均溫度初始值從上往下依次為17．74，18．03，19．74，22．19及24．40℃。

在上述模型的基礎(chǔ)上，利用抽水試驗(yàn)反演模型的水力學(xué)參數(shù)，利用2018年11月16日至2019年3月3日地下水源熱泵系統(tǒng)試運(yùn)行監(jiān)測(cè)數(shù)據(jù)反演巖土體熱物性參數(shù)，計(jì)算域土力學(xué)參數(shù)參考土工試驗(yàn)結(jié)果，水動(dòng)力彌散系數(shù)參考以往經(jīng)驗(yàn)系數(shù)［23-24］。

采用試錯(cuò)法來反演模型的水力學(xué)參數(shù)和熱物性參數(shù)，將不同參數(shù)下計(jì)算的地下水位與巖土體溫度與實(shí)測(cè)值對(duì)比，獲取擬合精度最好的參數(shù)。對(duì)比結(jié)果顯示：地下水位和巖土體溫度的實(shí)測(cè)值與趨勢(shì)相同且計(jì)算誤差小，所建立的模型可信。各層土體參數(shù)見表1，地下水位、巖土體溫度擬合結(jié)果見圖5—8。

圖5 抽水試驗(yàn)地下水位擬合對(duì)比圖

表1 各層土體參數(shù)一覽表

3．4 參數(shù)敏感性分析 系統(tǒng)特性X主要由n個(gè)因素x＝（x1，x2，…，xn）所決定，分別令各因素在其各自可能范圍內(nèi)變動(dòng)，分析系統(tǒng)特性X偏離基準(zhǔn)狀態(tài)X*的趨勢(shì)和程度，以確定這些因素的變化對(duì)系統(tǒng)特征的影響程度。若分析xi對(duì)特性X的影響，可令其余各參數(shù)基準(zhǔn)值固定不變，而令xi在其可能的范圍內(nèi)變動(dòng)，這時(shí)系統(tǒng)特性X表現(xiàn)為：

圖6 G2觀測(cè)井埋深55 m處溫度擬合對(duì)比圖

圖7 G2觀測(cè)井埋深90 m處溫度擬合對(duì)比圖

圖8 G2觀測(cè)井埋深125 m處溫度擬合對(duì)比圖

若xi的微小變化引起X的較大變化，表明X對(duì)xi敏感，此時(shí)xi是高敏感性參數(shù)；若xi的變化范圍較大而X變化微小時(shí)，表明X對(duì)xi不敏感，此時(shí)xi是低敏感性參數(shù)。本次研究的參數(shù)靈敏度應(yīng)為抽水井溫度和第Ⅲ承壓含水層Cl-濃度，由于本次研究建立的數(shù)值模型涉及到的參數(shù)較多，且具有不同的量綱。因此，在敏感性分析時(shí)定義參數(shù)的相對(duì)靈敏度來表征參數(shù)的敏感性［25］，公式表示如下：

表2 參數(shù)靈敏度計(jì)算結(jié)果一覽表

由計(jì)算結(jié)果可知：對(duì)于抽水井熱貫通影響縱向熱彌散度最大，熱容量次之。對(duì)于第Ⅲ承壓含水層Cl-濃度影響縱向彌散度最大，孔隙度次之。本次研究縱向熱彌散度、熱容量參數(shù)通過現(xiàn)場(chǎng)試驗(yàn)反演獲得，孔隙度根據(jù)室內(nèi)土工試驗(yàn)獲得，確保了模型的可靠性，而由于模型識(shí)別驗(yàn)證階段是時(shí)間較短，第Ⅲ承壓含水層Cl-濃度反應(yīng)不明顯，因此各層溶質(zhì)彌散度參考前期研究成果，以確保模型能客觀反映含水層Cl-濃度的變化特征。

4 模擬結(jié)果

4．1 原設(shè)計(jì)方案預(yù)測(cè)結(jié)果 根據(jù)地下水源熱泵系統(tǒng)的設(shè)計(jì)運(yùn)行方案，設(shè)計(jì)擬采用1組地下水源熱泵換熱機(jī)組，設(shè)計(jì)負(fù)荷見表3，按照應(yīng)用需求，提供整個(gè)科研實(shí)驗(yàn)樓所需的夏天制冷、冬天制熱的能量，確定工程夏季6—9月開機(jī)，共運(yùn)行122 d；冬季12、1、2月開機(jī)，共運(yùn)行90 d。

表3 科研實(shí)驗(yàn)樓工程冷熱負(fù)荷一覽表

科研實(shí)驗(yàn)樓設(shè)計(jì)利用溫差為5℃，每天運(yùn)行10 h，經(jīng)式（15）和式（16）計(jì)算得到夏季循環(huán)水量為943．88 m3/d，冬季循環(huán)水量為403．88 m3/d，抽汲的地下水以100%進(jìn)行回灌，利用識(shí)別、驗(yàn)證后的數(shù)值模型來模擬地下水源熱泵系統(tǒng)運(yùn)行10年第Ⅰ承壓含水層溫度和第Ⅲ承壓含水層濃度變化趨勢(shì)。計(jì)算結(jié)果見圖9—12。

圖9 夏季第Ⅰ承壓含水層水位等值線圖（單位：m）

圖11 第Ⅰ含水層平均溫度等值線圖（單位：℃）

圖12 第Ⅲ承壓含水層Cl-濃度等值線圖（單位：mg/L）

式中：Qcool、Qheat分別為夏季和冬季循環(huán)水量，m3/d；Pcool、Pheat分別為建筑的冷、熱負(fù)荷，kW；Pinput-cool、Pinput-heat分別為地源熱泵機(jī)組制冷和熱輸入功率，kW；ΔT為利用溫差，℃；toperation為每天運(yùn)行時(shí)間，h；cw為水的體積熱容量，J/（m3·℃）。

若按照設(shè)計(jì)方案運(yùn)行，10年后第Ⅰ承壓含水層中抽水井處水體溫度為20．59℃，上升了0．85℃；第Ⅲ承壓含水層Cl-濃度最大為156．11 mg/L，增大了28．11 mg/L。計(jì)算結(jié)果顯示，位于長(zhǎng)江三角洲地區(qū)的南通市由于夏季炎熱，冬季溫和，夏季的制冷周期和制冷負(fù)荷高于冬季制熱周期和制熱負(fù)荷，按照既有方案運(yùn)行地下水源熱泵系統(tǒng)會(huì)出現(xiàn)熱失衡問題，即出現(xiàn)輕度熱貫通現(xiàn)象［26］；深部淡水水質(zhì)咸化，即第Ⅲ承壓淡水Cl-濃度升高，而由于第Ⅲ承壓含水層上部弱透水層厚度較厚，Cl-濃度增大幅度有限。

4．2 改進(jìn)方案預(yù)測(cè)結(jié)果 為了科學(xué)合理開發(fā)淺層地?zé)崮?，控制地下水系統(tǒng)環(huán)境惡化，對(duì)江蘇省地質(zhì)工程勘察院（南通分院）科研實(shí)驗(yàn)樓地下水源熱泵系統(tǒng)運(yùn)行方案進(jìn)行比選，即在抽、灌井位置和建筑冷、熱負(fù)荷不變的基礎(chǔ)上，通過調(diào)整利用溫差和夏季灌采比來控制熱貫通和深層地下淡水咸化，利用已建立的數(shù)值模型，模擬不同運(yùn)行方案下的地下水位、巖土體溫度、第Ⅲ承壓含水層咸化發(fā)展變化趨勢(shì)，分析改進(jìn)方案對(duì)地下水環(huán)境惡化的控制效應(yīng)。

以往研究顯示，增大利用溫差和減小灌采比有利于控制熱貫通的發(fā)生［27］，本次研究采用增大利用溫差和減小夏季灌采比的改進(jìn)方案進(jìn)行分析，將夏季灌采比調(diào)整為80%和60%，冬季灌采比仍維持100%，以盡量減小夏季和冬季注入含水層的熱量平衡；同時(shí)將利用溫差調(diào)整為7和9℃，以減小含水層循環(huán)水量，弱化抽、灌井之間的相互影響，評(píng)價(jià)地下水源熱泵系統(tǒng)運(yùn)行10年地下水環(huán)境的變化，具體方案及模擬預(yù)測(cè)結(jié)果見表4、圖13和圖14。

表4 改進(jìn)運(yùn)行方案和模擬預(yù)測(cè)結(jié)果一覽表

圖13 抽水井溫度變化曲線

圖14 第Ⅲ承壓含水層Cl-濃度變化曲線

計(jì)算結(jié)果顯示采用方案8（利用溫差9℃、夏季灌采比60%）的方案抽水井溫度和第Ⅲ承壓含水層Cl-濃度上升幅度最小，在井間距、冷熱負(fù)荷不變的條件下，將利用溫差由5℃分別提高至7和9℃運(yùn)行地下水源熱泵系統(tǒng)10年，抽水井溫度分別降低了0．06和0．13℃，第Ⅲ承壓含水層Cl-濃度最大值分別降低了5．58和9．09 mg/L；將夏季灌采比分別從100%降低至80%、60%，運(yùn)行地下水源熱泵系統(tǒng)10年，與原運(yùn)行方案相比較抽水井溫度分別降低了0．29和0．53℃，第Ⅲ承壓含水層Cl-濃度最大值分別降低了2．73和5．20 mg/L。

5 結(jié)論

以江蘇省地質(zhì)工程勘察院（南通分院）地下水源熱泵系統(tǒng)為例，建立孔隙地下水滲流—熱量運(yùn)移—溶質(zhì)運(yùn)移多場(chǎng)耦合模型，研究地下水源熱泵系統(tǒng)運(yùn)行過程中地下水位、Cl-濃度和巖土體溫度的變化趨勢(shì)。獲得的主要結(jié)論如下：（1）按照原設(shè)計(jì)方案運(yùn)行10年，由于夏季冷負(fù)荷高于冬季熱負(fù)荷，第Ⅰ承壓含水層出現(xiàn)輕度熱貫通現(xiàn)象；同時(shí)在水動(dòng)力作用下，第Ⅰ承壓含水層的咸水入侵深部第Ⅲ承壓含水層淡水，Cl-濃度升高相對(duì)較小。（2）增大利用溫差和減小夏季灌采比能有效緩解熱貫通和深層淡水咸化的發(fā)展趨勢(shì)，增大利用溫差對(duì)深層地下淡水咸化的控制效應(yīng)更為顯著，而減小夏季灌采比對(duì)熱貫通的控制效應(yīng)更明顯。（3）通過模擬計(jì)算可以發(fā)現(xiàn)，調(diào)整利用溫差和夏季灌采比可實(shí)現(xiàn)地下水源熱泵系統(tǒng)運(yùn)行方案改進(jìn)，緩解熱貫通，對(duì)第Ⅲ承壓含水層咸化也存在影響，改進(jìn)運(yùn)行方案的本質(zhì)是通過調(diào)整循環(huán)水量改變地下水動(dòng)力場(chǎng)，進(jìn)而控制地下水系統(tǒng)環(huán)境的惡化。