趙志宏,劉桂宏,王佳鋮,徐浩然

(清華大學(xué) 土木工程系,北京 100084)

如期實(shí)現(xiàn)碳達(dá)峰、碳中和,是推動(dòng)我國(guó)經(jīng)濟(jì)高質(zhì)量發(fā)展和生態(tài)環(huán)境高水平保護(hù)的內(nèi)在要求,而建設(shè)清潔低碳安全高效的能源體系,是實(shí)現(xiàn)碳達(dá)峰、碳中和的必由之路?！吨腥A人民共和國(guó)國(guó)民經(jīng)濟(jì)和社會(huì)發(fā)展第十四個(gè)五年規(guī)劃和2035年遠(yuǎn)景目標(biāo)綱要》指出要因地制宜開(kāi)發(fā)利用地?zé)崮堋５責(zé)崮苁且环N清潔低碳、分布廣泛、資源豐富、安全穩(wěn)定的優(yōu)質(zhì)可再生能源[1-4]。根據(jù)埋藏深度可將地?zé)崮芊譃?類(lèi):200 m以淺為淺層地?zé)崮?200～3 000 m為中深層地?zé)崮?3 000 m以下為深層地?zé)崮?。目?已實(shí)現(xiàn)規(guī)模化開(kāi)發(fā)利用的水熱型地?zé)豳Y源在中深層和深層熱儲(chǔ)中均有廣泛分布。主要針對(duì)中深層和深層水熱型地?zé)豳Y源開(kāi)展研究,為描述簡(jiǎn)便,以下統(tǒng)稱(chēng)為深層地?zé)崮堋?/p>

地?zé)崮荛_(kāi)采過(guò)程涉及復(fù)雜的溫度場(chǎng)、滲流場(chǎng)、應(yīng)力場(chǎng)、化學(xué)場(chǎng)耦合效應(yīng)[5-7]。溫度場(chǎng)與滲流場(chǎng)的耦合是地?zé)崮荛_(kāi)采的核心過(guò)程,主要是回灌的尾水在熱儲(chǔ)層中滲流換熱重新恢復(fù)高溫。當(dāng)需要對(duì)熱儲(chǔ)層進(jìn)行增產(chǎn)穩(wěn)產(chǎn)改造時(shí),或當(dāng)考慮儲(chǔ)層變形、回灌堵塞等問(wèn)題時(shí),則需要在考慮溫度場(chǎng)與滲流場(chǎng)耦合的基礎(chǔ)上,進(jìn)一步考慮應(yīng)力場(chǎng)、化學(xué)場(chǎng)的耦合作用。數(shù)值模擬已成為研究評(píng)價(jià)深層地?zé)崮荛_(kāi)采多場(chǎng)耦合效應(yīng)的重要方法,比如TOUGH-FLAC、OpenGeoSys、COMSOL等多場(chǎng)耦合模擬軟件已廣泛應(yīng)用于深層資源評(píng)價(jià)與開(kāi)采方案優(yōu)化設(shè)計(jì)[8-26]。雖然多場(chǎng)耦合數(shù)值模型可以精細(xì)化地預(yù)測(cè)熱儲(chǔ)溫度場(chǎng)、滲流場(chǎng)、應(yīng)力場(chǎng)、化學(xué)場(chǎng)的動(dòng)態(tài)響應(yīng)過(guò)程,但計(jì)算效率通常不高,特別是對(duì)于復(fù)雜的井-儲(chǔ)模型。為此,將三維地?zé)峋?jiǎn)化為一維線(xiàn)單元,實(shí)現(xiàn)了地?zé)崽锍叨葟?fù)雜群井系統(tǒng)多場(chǎng)耦合效應(yīng)的高效模擬[7,10,23,27-28]。盡管考慮多場(chǎng)耦合效應(yīng)的深層熱儲(chǔ)數(shù)值模型與算法已日臻成熟,但在深層地?zé)崮荛_(kāi)采過(guò)程中的多場(chǎng)耦合災(zāi)變效應(yīng)仍未受到足夠的重視,尚未形成統(tǒng)一的標(biāo)準(zhǔn)體系對(duì)地?zé)衢_(kāi)采多場(chǎng)耦合災(zāi)變效應(yīng)進(jìn)行評(píng)價(jià)與調(diào)控。

本文總結(jié)筆者團(tuán)隊(duì)近年來(lái)針對(duì)水熱型地?zé)豳Y源開(kāi)發(fā)利用多場(chǎng)耦合效應(yīng)開(kāi)展的數(shù)值模擬研究成果,提出地?zé)崮芸沙掷m(xù)開(kāi)采的多指標(biāo)評(píng)價(jià)體系,充分發(fā)揮地?zé)峋?儲(chǔ)系統(tǒng)多場(chǎng)耦合數(shù)值模型的優(yōu)勢(shì),建立面向城市深層地?zé)崮芸沙掷m(xù)開(kāi)采的資源量評(píng)價(jià)與采灌方案優(yōu)化設(shè)計(jì)方法,旨在促進(jìn)城市深層地?zé)崮芨笠?guī)模和更高質(zhì)量的開(kāi)發(fā)利用。

1 城市深層地?zé)崮芸沙掷m(xù)開(kāi)采的多指標(biāo)評(píng)價(jià)體系

國(guó)家能源行業(yè)標(biāo)準(zhǔn)《地?zé)崮苄g(shù)語(yǔ)》(NB/T 10097—2018)中,有3個(gè)術(shù)語(yǔ)的定義均提到了“可持續(xù)”的標(biāo)準(zhǔn)[29]:① “可開(kāi)采量”,在地?zé)崽锟辈椤㈤_(kāi)采和監(jiān)測(cè)的基礎(chǔ)上,考慮到可持續(xù)開(kāi)發(fā),經(jīng)擬合計(jì)算允許每年合理開(kāi)采的地?zé)崃黧w量和熱量;② “動(dòng)態(tài)監(jiān)測(cè)”,地?zé)豳Y源在勘探、開(kāi)采及停采階段,連續(xù)記錄水位、井口溫度、井口壓力、開(kāi)采量、回灌量和蒸汽比例等,并定時(shí)分析地?zé)崃黧w化學(xué)組分和同位素的過(guò)程,基于此來(lái)判斷熱儲(chǔ)溫度、壓力、流體化學(xué)組份含量及資源量的動(dòng)態(tài)變化,為地?zé)豳Y源的可持續(xù)利用與管理提供依據(jù);③ “地?zé)崮軆?yōu)化開(kāi)采”,采用最優(yōu)化的方法開(kāi)采熱儲(chǔ),在可持續(xù)且不會(huì)帶來(lái)環(huán)境危害的基礎(chǔ)上獲得最大采熱量?？梢?jiàn),可持續(xù)已成為地?zé)豳Y源勘探、評(píng)價(jià)、開(kāi)發(fā)的重要前提,然而對(duì)于可持續(xù)的評(píng)價(jià)指標(biāo)及其標(biāo)準(zhǔn)仍未形成共識(shí)。以下從溫度、水位、變形3個(gè)方面給出深層地?zé)崮芸沙掷m(xù)開(kāi)采的評(píng)價(jià)標(biāo)準(zhǔn):

(1)熱突破時(shí)間。開(kāi)采井井口溫度下降不超過(guò)允許值,本文暫定溫降2 ℃對(duì)應(yīng)的時(shí)間為熱突破時(shí)間,本文暫定最短熱突破時(shí)間為100 a。

(2)水位降深。開(kāi)采井水位降深不超過(guò)允許值,本文暫定允許水位降深為-150 m。

(3)地表變形。長(zhǎng)期采灌引起的地面變形不超過(guò)允許值,本文暫定地表變形允許值為20 mm。

開(kāi)采井需同時(shí)滿(mǎn)足上述3個(gè)環(huán)境約束條件才可運(yùn)行,否則應(yīng)予以關(guān)停,即地?zé)峋畨勖?tlife)由上述3個(gè)約束條件共同決定。因此,可采地?zé)崮蹺prod的計(jì)算式為

Eprod=qprotlifeρfCp,f(Tpro-Tin)

(1)

其中,qpro為開(kāi)采流量;ρf為流體密度;Cp,f為流體比熱容;Tpro為開(kāi)采溫度;Tin為回灌溫度。前人的研究中,tlife在30～100 a間變化,本文按tlife= 100 a考慮,是對(duì)深層地?zé)崮芸沙掷m(xù)開(kāi)采較為嚴(yán)格的限制。因此,需因地制宜選擇地?zé)衢_(kāi)采模式,并對(duì)采灌井布局、井深、采灌量、回灌溫度以及回灌方式等參數(shù)進(jìn)行優(yōu)化設(shè)計(jì),在可持續(xù)開(kāi)采條件下最大限度地發(fā)揮地?zé)峋淖匀划a(chǎn)能。

需要特別說(shuō)明的是,深層地?zé)崮荛_(kāi)發(fā)利用與地表沉降的相關(guān)性研究仍在進(jìn)行當(dāng)中,本文將地表變形作為限制性條件之一,而且允許值極為嚴(yán)苛,主要目的是完全杜絕深層地?zé)崮荛_(kāi)采對(duì)地表變形的影響。根據(jù)干涉合成孔徑雷達(dá)(InSAR)數(shù)據(jù),雄安新區(qū)在2015—2019年的地表沉降速率為37.9～104.8 mm/a,有學(xué)者推斷地表沉降的原因主要是地?zé)豳Y源的開(kāi)發(fā)利用[30]。基于多孔介質(zhì)有效應(yīng)力原理,有學(xué)者認(rèn)為雄安新區(qū)地表沉降與第四系地下水的超采有關(guān),其中地下水超采是地面沉降的主要因素,砂巖熱儲(chǔ)層的無(wú)序且沒(méi)有采灌均衡的開(kāi)發(fā)成為部分地區(qū)誘發(fā)沉降的次要因素[31]。根據(jù)水準(zhǔn)測(cè)量數(shù)據(jù),2015—2019年山東德城區(qū)500 m以淺地層的地表沉降速率為37.9 mm/a,500 m以深地層的地表沉降速率為11.1 mm/a,其中500～800 m層段土體壓縮量約為8.2 mm/a,800 m以深地層壓縮量約為2.9 mm/a。因此,在多場(chǎng)耦合作用下,深層地?zé)衢_(kāi)采確有可能導(dǎo)致儲(chǔ)層變形,但對(duì)地表沉降的影響有限[7],后續(xù)應(yīng)結(jié)合實(shí)測(cè)數(shù)據(jù),并建立數(shù)值模型精細(xì)評(píng)估深層地?zé)衢_(kāi)采對(duì)地表沉降的影響規(guī)律。

2 城市深層地?zé)峋?儲(chǔ)系統(tǒng)多場(chǎng)耦合模擬方法

深層地?zé)嵯到y(tǒng)主要包括兩部分,一是蓋層、儲(chǔ)層、斷層等地質(zhì)結(jié)構(gòu),二是回灌井、開(kāi)采井等工程結(jié)構(gòu)。地面換熱系統(tǒng)不在本文考慮范圍之內(nèi),相關(guān)內(nèi)容見(jiàn)文獻(xiàn)[31]。以下簡(jiǎn)要介紹深層地?zé)峋?儲(chǔ)系統(tǒng)多場(chǎng)耦合模型及非均質(zhì)儲(chǔ)層構(gòu)造方法[7,26]。

2.1 井-儲(chǔ)系統(tǒng)多場(chǎng)耦合數(shù)值模型

圖1給出了井-儲(chǔ)系統(tǒng)多場(chǎng)耦合模型原理。將熱儲(chǔ)層按飽和孔隙介質(zhì)考慮,采用多孔介質(zhì)熱彈性理論描述熱儲(chǔ)層應(yīng)力、滲流、傳熱三場(chǎng)耦合過(guò)程。蓋層與儲(chǔ)層類(lèi)似,但不考慮其中的滲流過(guò)程。將斷層按飽和的光滑平板模型來(lái)考慮,分別采用簡(jiǎn)化的線(xiàn)性模型描述其法向、切向變形過(guò)程,同時(shí)考慮剪脹對(duì)裂隙變形的影響;采用立方定律描述裂隙滲流過(guò)程;傳熱過(guò)程包括熱對(duì)流、熱傳導(dǎo)、熱交換等。為了提高計(jì)算效率而又不降低精度,引入一維線(xiàn)單元對(duì)三維井筒結(jié)構(gòu)進(jìn)行簡(jiǎn)化,考慮沿井筒軸向的對(duì)流傳熱過(guò)程,井筒內(nèi)流體與圍巖的熱交換通過(guò)等效換熱系數(shù)來(lái)近似考慮。具體的控制方程可見(jiàn)文獻(xiàn)[7,26],本文不再贅述。熱儲(chǔ)層和熱蓋層采用四面體單元進(jìn)行剖分,地?zé)峋捎镁€(xiàn)單元進(jìn)行剖分,為兼顧計(jì)算精度與效率,地?zé)峋車(chē)爸饕獰醿?chǔ)層的網(wǎng)格尺寸應(yīng)相對(duì)細(xì)化,而在遠(yuǎn)離地?zé)峋吧w層區(qū)域則可使用較粗化網(wǎng)格,有助于減少網(wǎng)格數(shù)量,節(jié)省計(jì)算時(shí)間。建議地?zé)峋車(chē)木W(wǎng)格尺寸不大于井徑的10倍。

圖1 井-儲(chǔ)系統(tǒng)多場(chǎng)耦合數(shù)值模型原理[26]Fig.1 Schematic diagram of coupled multi-field numerical model of well-reservoir system[26]

2.2 非均質(zhì)儲(chǔ)層構(gòu)建方法

熱儲(chǔ)層的滲透率、比熱容等參數(shù)具有較強(qiáng)的非均質(zhì)性,其空間非均質(zhì)性不僅決定了流體的流動(dòng)過(guò)程,還影響了傳熱和變形過(guò)程[21,24,32]。轉(zhuǎn)向帶法[33-34]、傅里葉積分法[35]等各類(lèi)地質(zhì)統(tǒng)計(jì)學(xué)方法已被應(yīng)用于生成非均質(zhì)熱儲(chǔ)模型。轉(zhuǎn)向帶法的主要特點(diǎn)是將高維的隨機(jī)過(guò)程簡(jiǎn)化為一維情況,而傅里葉積分法可基于多種協(xié)方差函數(shù)構(gòu)建非均質(zhì)模型。無(wú)論采用何種方法生成隨機(jī)模型,其結(jié)果均具有較好的一致性。決定熱儲(chǔ)參數(shù)非均勻分布特征的控制參數(shù)主要有平均值、方差、相關(guān)長(zhǎng)度等。生成非均質(zhì)熱儲(chǔ)模型的通常做法是(圖2中,σ為標(biāo)準(zhǔn)差;x為隨機(jī)變量;μ為數(shù)學(xué)期望):① 選定合適的熱儲(chǔ)參數(shù)非均勻分布形式(正態(tài)分布、對(duì)數(shù)正態(tài)分布等)及其控制參數(shù)(平均值、方差、相關(guān)長(zhǎng)度等);② 采用地質(zhì)統(tǒng)計(jì)學(xué)方法生成隨機(jī)參數(shù)場(chǎng);③ 將熱儲(chǔ)層離散為一定數(shù)量的代表性體積單元,賦予隨機(jī)參數(shù)場(chǎng);④ 基于熱儲(chǔ)參數(shù)間的物理關(guān)系,生成其他熱儲(chǔ)參數(shù)的隨機(jī)分布場(chǎng)。

圖2 非均質(zhì)熱儲(chǔ)模型生成過(guò)程示意Fig.2 Schematic diagram of generating heterogeneous geothermal reservoir models

3 基于數(shù)值模擬的可采資源量評(píng)價(jià)

3.1 可采資源量評(píng)價(jià)方法

國(guó)家標(biāo)準(zhǔn)《地?zé)豳Y源地質(zhì)勘查規(guī)范》(GB/T 11615—2010)[36]中給出了地?zé)釤崃髁糠?、熱?chǔ)法、解析模型法、統(tǒng)計(jì)分析法、數(shù)值模型法、比擬法等地?zé)豳Y源/儲(chǔ)量計(jì)算方法。在地?zé)崽锏目辈槌潭缺容^高,并且具有一定時(shí)期的開(kāi)采歷史和比較齊全的監(jiān)測(cè)資料時(shí),尤其是城市地區(qū)的地?zé)崽?應(yīng)建立地?zé)崽锏臄?shù)值模擬模型,用以計(jì)算/評(píng)價(jià)地?zé)醿?chǔ)量,并作為地?zé)崽锕芾淼墓ぞ?。因?筆者基于深層地?zé)峋?儲(chǔ)系統(tǒng)多場(chǎng)耦合模擬方法,計(jì)算滿(mǎn)足可持續(xù)開(kāi)采標(biāo)準(zhǔn)的地?zé)豳Y源可采量,具體流程如下:

(1)地質(zhì)模型構(gòu)建。根據(jù)地?zé)崽镞吔缗c評(píng)價(jià)區(qū)域確定研究區(qū)范圍,逐層導(dǎo)入各地層頂?shù)装?、斷層面幾何?shù)據(jù),構(gòu)建所有地層界面、斷層面,并通過(guò)地層界面、斷層面、模型實(shí)體間的相互分割,刪除多余部分,形成研究區(qū)三維地質(zhì)模型。

(2)井-儲(chǔ)模型構(gòu)建。在地質(zhì)模型的基礎(chǔ)上,按照地?zé)峋鴺?biāo)插入一維線(xiàn)單元代表地?zé)峋?線(xiàn)單元的長(zhǎng)度即井深。根據(jù)研究區(qū)斷層帶位置和熱儲(chǔ)層的富水性差異,對(duì)熱儲(chǔ)參數(shù)進(jìn)行分區(qū),并賦予合適的模型參數(shù)。設(shè)置合適的溫度場(chǎng)、滲流場(chǎng)、地應(yīng)力場(chǎng)初始條件,并根據(jù)研究區(qū)熱儲(chǔ)層地?zé)岬刭|(zhì)構(gòu)造特征、地?zé)嵝纬蓹C(jī)理、水文地質(zhì)條件,確定整個(gè)模型熱儲(chǔ)系統(tǒng)的熱水力邊界條件。根據(jù)采灌量、回灌溫度等數(shù)據(jù),為每一眼地?zé)峋O(shè)置溫度和流體的流入流出邊界條件。對(duì)于數(shù)值模型的校檢,可以對(duì)地?zé)崽镞^(guò)去某時(shí)期的采灌活動(dòng)進(jìn)行模擬計(jì)算,將模擬結(jié)果與現(xiàn)場(chǎng)監(jiān)測(cè)數(shù)據(jù)相對(duì)比,以此檢驗(yàn)?zāi)Ｐ徒⒌暮侠硇?。主要通過(guò)校正熱儲(chǔ)層參數(shù)使模擬結(jié)果與監(jiān)測(cè)數(shù)據(jù)相吻合,并最終確定模型參數(shù)。建模的滲流過(guò)程中,可調(diào)整滲透率、孔隙率和儲(chǔ)水系數(shù)使水位模擬值與水位監(jiān)測(cè)數(shù)據(jù)相吻合。而在傳熱過(guò)程中,可調(diào)整地溫梯度、比熱容及導(dǎo)熱系數(shù)等,使溫度模擬值與測(cè)溫?cái)?shù)據(jù)相吻合。

(3)開(kāi)采方案優(yōu)選。在采灌開(kāi)始前需先進(jìn)行地應(yīng)力平衡,計(jì)算模型在初始熱水力多場(chǎng)耦合作用下的平衡狀態(tài),然后再運(yùn)行采灌條件下的數(shù)值模型,根據(jù)生產(chǎn)需求確定計(jì)算時(shí)長(zhǎng),考慮到采灌周期以月為單位,時(shí)步設(shè)置應(yīng)不大于1個(gè)月。為數(shù)值模型設(shè)置時(shí)間周期函數(shù),將每年離散成供暖季和非供暖季2個(gè)時(shí)間段。計(jì)算結(jié)果以開(kāi)采井水位、溫度和沉降量變化為主。

借鑒地?zé)衢_(kāi)發(fā)誘發(fā)地震風(fēng)險(xiǎn)管控的“紅綠燈系統(tǒng)”[37],針對(duì)不同開(kāi)采方案下開(kāi)采井的可持續(xù)開(kāi)采風(fēng)險(xiǎn)進(jìn)行分級(jí)評(píng)價(jià)(表1)。

表1 地?zé)峋沙掷m(xù)開(kāi)采風(fēng)險(xiǎn)等級(jí)評(píng)價(jià)標(biāo)準(zhǔn)

3.2 工程案例

以雄安新區(qū)容東安置區(qū)(面積約13 km2)為案例,計(jì)算該區(qū)域的地?zé)豳Y源可采量。從容城地?zé)崽锏刭|(zhì)模型中截取邊長(zhǎng)為10 km、深為5 km的區(qū)域建立井-儲(chǔ)數(shù)值模型(圖3),具體參數(shù)取值見(jiàn)表2。模型內(nèi)現(xiàn)有12口開(kāi)采井和10口回灌井,主要位于原容城縣城,另在容東安置區(qū)規(guī)劃建設(shè)20口地?zé)峋?。將該?儲(chǔ)模型共剖分為四面體單元約26萬(wàn)個(gè),一維線(xiàn)單元2 000個(gè)。根據(jù)水位監(jiān)測(cè)結(jié)果,以地表作為壓力基準(zhǔn),通過(guò)水柱高度換算構(gòu)建初始水壓場(chǎng)。根據(jù)氣象及有關(guān)地?zé)岬刭|(zhì)資料,以恒溫層帶溫度為基準(zhǔn),通過(guò)地溫梯度獲得熱儲(chǔ)溫度,并插值得到初始溫度場(chǎng)。聯(lián)合深孔水壓致裂法與壓磁電感法測(cè)量熱儲(chǔ)地應(yīng)力,并通過(guò)施加邊界條件形成初始地應(yīng)力場(chǎng)。對(duì)于滲流過(guò)程,側(cè)邊界為恒壓約束,壓力分布與深度呈線(xiàn)性相關(guān);上、下邊界不透水。對(duì)于傳熱,側(cè)邊界為開(kāi)邊界,溫度由地表溫度和插值得到的地溫梯度場(chǎng)根據(jù)深度計(jì)算得到,模型的頂部和底部邊界為熱絕緣邊界。對(duì)于力學(xué),側(cè)邊界為應(yīng)力約束,大小主應(yīng)力分別為(-0.033 1z+7.529 6)MPa和(-0.020 8z+6.092 3)MPa,其中z為深度;上邊界是自由的,但下邊界是固定的。對(duì)于地?zé)峋木?回灌井設(shè)置為回灌流量和回灌溫度,開(kāi)采井設(shè)置為開(kāi)采流量和熱流出。對(duì)于地?zé)峋阊鄱?回灌井設(shè)置為質(zhì)量流量和熱流出,開(kāi)采井設(shè)置為質(zhì)量流量和溫度。通過(guò)監(jiān)測(cè)井水位變化和測(cè)溫?cái)?shù)據(jù),對(duì)模型進(jìn)行校準(zhǔn)。校準(zhǔn)后的數(shù)值模型能夠較為準(zhǔn)確地反映熱儲(chǔ)中的熱-水-力多場(chǎng)耦合過(guò)程,可用于下一步的計(jì)算。

圖3 雄安新區(qū)容東安置區(qū)井-儲(chǔ)系統(tǒng)多場(chǎng)耦合數(shù)值模型[26]Fig.3 Coupled multi-field numerical model of well-reservoir system in the Rongdong resettlement area in Xiongan New Area[26]

為評(píng)價(jià)容東安置區(qū)深層地?zé)崮芸刹闪?擬布置地?zé)峋?0口,井深2 000 m,井距500 m。在供暖季,回灌井和開(kāi)采井以相同的流量200 m3/h運(yùn)行,而在非供暖季地?zé)峋贿\(yùn)行。在采灌均衡條件下,考慮3種布井方案,包括“軌道式”、“棋盤(pán)式”和“集中采灌式”(圖4)。

圖4 雄安新區(qū)容東安置區(qū)熱儲(chǔ)初始溫度場(chǎng)、流場(chǎng)與布井方案[26]Fig.4 Initial temperature field,flow field and well layout scheme in the Rongdong resettlement area in Xiongan New Area[26]

考慮到熱儲(chǔ)流體運(yùn)動(dòng)方向,開(kāi)采井布置在上游區(qū)域,防止開(kāi)采井過(guò)早發(fā)生熱突破。

圖5對(duì)3種布井方案的水位場(chǎng)、溫度場(chǎng)和位移場(chǎng)進(jìn)行了比較。圖6給出了3種布井方案中10口開(kāi)采井的開(kāi)采風(fēng)險(xiǎn)等級(jí)分布。在“軌道式”布井方案中,開(kāi)采井0～5號(hào)和9號(hào)只有熱突破風(fēng)險(xiǎn),開(kāi)采風(fēng)險(xiǎn)等級(jí)較低;開(kāi)采井6～8號(hào)同時(shí)面臨熱突破和垂直位移風(fēng)險(xiǎn),開(kāi)采風(fēng)險(xiǎn)等級(jí)為較高。在“棋盤(pán)式”布井方案中,開(kāi)采井0～2、4、7和9號(hào)的開(kāi)采風(fēng)險(xiǎn)等級(jí)為低,均滿(mǎn)足地?zé)豳Y源可持續(xù)開(kāi)發(fā)利用的評(píng)價(jià)標(biāo)準(zhǔn);開(kāi)采井3、5、6和8號(hào)只有垂直位移風(fēng)險(xiǎn),開(kāi)采風(fēng)險(xiǎn)等級(jí)為較低。在“集中采灌式”布井方案中,開(kāi)采井0、1和4號(hào)只有水位風(fēng)險(xiǎn),開(kāi)采風(fēng)險(xiǎn)等級(jí)為較低;開(kāi)采井3、6、8和9號(hào)同時(shí)面臨熱突破和垂直位移風(fēng)險(xiǎn),開(kāi)采井7號(hào)同時(shí)面臨水位和垂直位移風(fēng)險(xiǎn),這5口開(kāi)采井開(kāi)采風(fēng)險(xiǎn)等級(jí)為較高;開(kāi)采井2和5號(hào)完全不滿(mǎn)足地?zé)豳Y源可持續(xù)開(kāi)發(fā)利用的評(píng)價(jià)標(biāo)準(zhǔn),開(kāi)采風(fēng)險(xiǎn)等級(jí)為高。以上結(jié)果表明:在均質(zhì)模型中,“棋盤(pán)式”布井方案開(kāi)采風(fēng)險(xiǎn)最低,其次是“軌道式”布井方案,而“集中采灌式”布井方案的開(kāi)采風(fēng)險(xiǎn)最高?；诘?zé)豳Y源可持續(xù)開(kāi)發(fā)利用的評(píng)價(jià)標(biāo)準(zhǔn),可以得到每口開(kāi)采井的壽命,根據(jù)壽命計(jì)算了3種布井方案在100 a內(nèi)的可采地?zé)崮芸偭?如圖7所示。

此外,假設(shè)容東安置區(qū)熱儲(chǔ)的滲透率服從均值為4.32×10-12m2、方差為1的正態(tài)分布,相關(guān)長(zhǎng)度為300 m,采用第2.2節(jié)的方法生成非均質(zhì)儲(chǔ)層模型。設(shè)置與均質(zhì)模型相同的初始和邊界條件,對(duì)非均質(zhì)熱儲(chǔ)模型“軌道式”、“棋盤(pán)式”和“集中采灌式”3種布井方案進(jìn)行對(duì)比分析。與均質(zhì)熱儲(chǔ)模型中的近似橢圓形或圓形冷鋒面相比,非均質(zhì)熱儲(chǔ)模型中的冷鋒面由于優(yōu)勢(shì)滲流通道而變得不規(guī)則,由回灌引起的垂直位移在非均質(zhì)熱儲(chǔ)模型中不如均質(zhì)熱儲(chǔ)模型顯著(圖5)。與均質(zhì)熱儲(chǔ)模型相比,水位下降到-150 m以下的地?zé)峋當(dāng)?shù)量明顯減少,而熱突破或超過(guò)臨界垂直位移的地?zé)峋當(dāng)?shù)量也有減少,這是因?yàn)椴灰?guī)則的流動(dòng)通道阻礙了熱突破,降低了開(kāi)采井的熱收縮。因此,部分開(kāi)采井從由均質(zhì)模型中的高風(fēng)險(xiǎn)變?yōu)榉蔷|(zhì)模型中的低風(fēng)險(xiǎn)。結(jié)合圖6和圖7可以得出“集中采灌式”是非均質(zhì)熱儲(chǔ)層的最優(yōu)群井開(kāi)發(fā)模式,在實(shí)際工程中應(yīng)充分考慮熱儲(chǔ)非均質(zhì)特性的影響,通過(guò)優(yōu)化井位降低開(kāi)采井風(fēng)險(xiǎn)。

圖5 3種布井方案的水位場(chǎng)、溫度場(chǎng)和位移場(chǎng)[26]Fig.5 Water level field,temperature field and displacement field of three well layout schemes[26]

圖6 3種布井方案下10口開(kāi)采井的風(fēng)險(xiǎn)等級(jí)對(duì)比(彩色圓圈代表開(kāi)采井,白色圓圈代表回灌井)Fig.6 Comparison of risk rating of 10 production wells under three well layout schemes (Colorful circles are production wells,white circles are injection wells)

圖7 容東安置區(qū)不同采灌方案100 a可采地?zé)崮蹻ig.7 Recoverable geothermal energy of 100 years under different exploitation schemes in the Rongdong resettlement area

4 基于代理模型的地?zé)崮荛_(kāi)采方案優(yōu)化設(shè)計(jì)

深層地?zé)豳Y源開(kāi)發(fā)利用中,開(kāi)采參數(shù)的優(yōu)化設(shè)計(jì)對(duì)提高地?zé)岵墒章示哂兄匾绊?但目前還僅停留在基于敏感性的優(yōu)化設(shè)計(jì)。通過(guò)優(yōu)化開(kāi)采參數(shù),在最大限度可持續(xù)開(kāi)采地?zé)崮艿耐瑫r(shí),避免誘發(fā)地下水位下降、熱突破和地面沉降等環(huán)境災(zāi)害。根據(jù)地?zé)衢_(kāi)發(fā)方案的需求,基于地?zé)峋?儲(chǔ)系統(tǒng)多場(chǎng)耦合數(shù)值模型,結(jié)合實(shí)驗(yàn)設(shè)計(jì)、代理模型以及多響應(yīng)優(yōu)化分析確定地?zé)峥刹蓞^(qū)的優(yōu)化對(duì)井開(kāi)發(fā)模式,以實(shí)現(xiàn)最大限度地可持續(xù)開(kāi)發(fā)利用地?zé)豳Y源,但儲(chǔ)層參數(shù)的不確定性對(duì)可采地?zé)崮艿挠绊懖辉诒疚难芯康姆秶鷥?nèi)。

4.1 開(kāi)采參數(shù)優(yōu)化設(shè)計(jì)方法

地?zé)衢_(kāi)采參數(shù)優(yōu)化工作流程如圖8所示。首先,根據(jù)地?zé)衢_(kāi)發(fā)方案的需求,收集地?zé)崽锏牡責(zé)岬刭|(zhì)資料,在此基礎(chǔ)上建立三維熱儲(chǔ)地質(zhì)模型,并利用監(jiān)測(cè)數(shù)據(jù)對(duì)模型進(jìn)行校準(zhǔn)。然后,根據(jù)實(shí)驗(yàn)設(shè)計(jì)矩陣,設(shè)定人為控制參數(shù)后運(yùn)行數(shù)值模型,并進(jìn)行統(tǒng)計(jì)分析,得到熱儲(chǔ)響應(yīng)對(duì)人為控制參數(shù)的敏感性大小。最后,建立熱儲(chǔ)響應(yīng)和人為控制參數(shù)之間的代理模型,通過(guò)數(shù)值模型驗(yàn)證代理模型的準(zhǔn)確性,并利用多響應(yīng)優(yōu)化方法,得到人為控制參數(shù)的最優(yōu)組合來(lái)滿(mǎn)足地?zé)衢_(kāi)發(fā)方案中在地?zé)衢_(kāi)采區(qū)最大限度地開(kāi)采地?zé)崮艿男枨?。具體流程如下:

(1)選定擬優(yōu)化的開(kāi)采參數(shù)并確定其上下限,指定熱儲(chǔ)響應(yīng)。本文選取采灌量、回灌溫度、采灌井間距、井深和回灌井方位等參數(shù)作為主控因子,根據(jù)地質(zhì)模型和地?zé)峋\(yùn)行狀態(tài)確定每個(gè)因子的上下限,將熱突破時(shí)間、水位降深、地表變形和可采地?zé)崮茏鳛闊醿?chǔ)響應(yīng)。

(2)基于Plackett-Burman設(shè)計(jì)構(gòu)建熱儲(chǔ)響應(yīng)代理模型。Plackett-Burman是二水平的部分試驗(yàn)設(shè)計(jì),對(duì)每個(gè)因子分別取高低兩水平,即上限“+”和下限“-”,通過(guò)比較各個(gè)因子兩水平之間的差異來(lái)確定因子的顯著性,最終能夠快速、有效地完成影響因素的篩選。Plackett-Burman設(shè)計(jì)矩陣隨機(jī)生成,每列所包含的高低水平數(shù)相等,即在N次試驗(yàn)中,每個(gè)因子的高低水平各出現(xiàn)N/2次。使用第2.1節(jié)中的城市深層地?zé)峋?儲(chǔ)系統(tǒng)多場(chǎng)耦合數(shù)值模型單獨(dú)運(yùn)行每次實(shí)驗(yàn)。采用多元回歸方法擬合熱儲(chǔ)響應(yīng)與人為控制參數(shù)之間的多項(xiàng)式代理模型。

(3)開(kāi)采參數(shù)優(yōu)化。多響應(yīng)優(yōu)化的目的是解決響應(yīng)之間的沖突,并確定聯(lián)合優(yōu)化單個(gè)響應(yīng)或多個(gè)響應(yīng)的參數(shù)設(shè)置組合。解決多響應(yīng)之間沖突的一種主流策略是降維,該策略將多響應(yīng)模型降為具有單個(gè)合意性函數(shù)的模型,然后將其作為單目標(biāo)優(yōu)化問(wèn)題進(jìn)行求解。合意性函數(shù)法已被廣泛用于解決各種多響應(yīng)優(yōu)化問(wèn)題,它將每個(gè)響應(yīng)轉(zhuǎn)換為單個(gè)合意性函數(shù),然后聚合為復(fù)合合意性,其定義為所有單個(gè)合意性函數(shù)的加權(quán)幾何均值。本文通過(guò)尋找一組人為控制參數(shù)取值使地?zé)衢_(kāi)采區(qū)的可采地?zé)崮茏畲蠡?而不超過(guò)第2.1節(jié)中定義的水位埋深、熱突破時(shí)間和垂直位移的臨界值,按這3個(gè)環(huán)境約束條件同等重要考慮,即權(quán)重相同,暫未考慮賦予3個(gè)環(huán)境約束條件不同的權(quán)重。首先獲得每個(gè)響應(yīng)的單個(gè)合意性(d),然后將單個(gè)合意性加權(quán)以形成復(fù)合合意性(D)。采用帝國(guó)主義競(jìng)爭(zhēng)算法(ICA)使復(fù)合合意性(D)最大化。

圖8 地?zé)衢_(kāi)采方案優(yōu)化設(shè)計(jì)流程[39]Fig.8 Flow chart of optimization design of geothermal exploitation schemes[39]

4.2 工程案例

以北京城市副中心及周邊區(qū)域作為案例,針對(duì)地?zé)釋?duì)井系統(tǒng)開(kāi)采參數(shù)進(jìn)行優(yōu)化設(shè)計(jì),實(shí)現(xiàn)地?zé)豳Y源可開(kāi)采量的最大化。從北京通州區(qū)三維地質(zhì)模型[38]中截取邊長(zhǎng)為10 km、深度為5 km的區(qū)域建立數(shù)值模型,在提高計(jì)算效率的同時(shí)且有較好的代表性,具體參數(shù)取值見(jiàn)表3。該區(qū)域現(xiàn)有7口開(kāi)采井和4口回灌井,井深為2 203～3 002 m,開(kāi)采溫度為45～51 ℃,采灌量為1 603～2 573 m3/d。所建立的三維井-儲(chǔ)數(shù)值模型如圖9所示,共剖分四面單元約20萬(wàn)個(gè),三角形網(wǎng)格約3萬(wàn)個(gè),一維線(xiàn)單元約500個(gè)。模型的初始溫度場(chǎng)可根據(jù)研究區(qū)內(nèi)各地層的地溫梯度來(lái)確定,初始?jí)毫?chǎng)根據(jù)研究區(qū)水位分布監(jiān)測(cè)結(jié)果,以地表作為壓力基準(zhǔn)面,進(jìn)行換算來(lái)確定。初始地應(yīng)力通過(guò)水力壓裂法對(duì)北京地區(qū)的5個(gè)鉆孔進(jìn)行原位應(yīng)力測(cè)量得到。對(duì)于滲流過(guò)程,側(cè)邊界為恒壓約束,壓力分布與深度呈線(xiàn)性相關(guān);上、下邊界不透水。對(duì)于傳熱,側(cè)邊界為開(kāi)邊界,溫度場(chǎng)由插值得到,模型的頂部和底部邊界為熱絕緣邊界。對(duì)于力學(xué),側(cè)邊界為應(yīng)力約束,大小主應(yīng)力分別為(-0.032 8z+2.5)MPa和(-0.022 1z+2)MPa,上邊界是自由的,但下邊界是固定的。對(duì)于地?zé)峋木?回灌井設(shè)置為回灌流量和回灌溫度,開(kāi)采井設(shè)置為開(kāi)采流量和熱流出。對(duì)于地?zé)峋阊鄱?回灌井設(shè)置為質(zhì)量流量和熱流出,開(kāi)采井設(shè)置為質(zhì)量流量和溫度。

擬在副中心東北、正南2個(gè)位置各布置一對(duì)地?zé)峋?需對(duì)采灌方案進(jìn)行優(yōu)化設(shè)計(jì)。首先確定擬優(yōu)化的開(kāi)采參數(shù)并確定其上下限,見(jiàn)表4。采用Plackett-Burman設(shè)計(jì)給出敏感性計(jì)算方案(表5),需運(yùn)行模型共12次,其中開(kāi)采參數(shù)在上、下限水平下各運(yùn)行6次。在回灌開(kāi)始前先進(jìn)行地應(yīng)力平衡,計(jì)算模型在初始熱-水-力多場(chǎng)耦合作用下的平衡狀態(tài),然后將采灌井按照Plackett-Burman設(shè)計(jì)中的參數(shù)水平設(shè)置并運(yùn)行模型,所有模型的初始和邊界條件都相同。計(jì)算時(shí)長(zhǎng)為100 a,并將每年分為供暖季4個(gè)月,非供暖季8個(gè)月,時(shí)步取1個(gè)月,在此期間,連續(xù)監(jiān)測(cè)儲(chǔ)層模型中的水位、溫度和位移的演變。使用1臺(tái)CPU為i7-4790 K、內(nèi)存為32 GB的普通電腦,每次模擬平均需要12 h完成。

表3 通州地?zé)崽锬Ｐ蛥?shù)[39]

圖9 通州地?zé)崽锶S地質(zhì)模型示意[39]Fig.9 Schematic diagram of three-dimensional geological model of Tongzhou geothermal field[39]

表4 北京城市副中心地?zé)釋?duì)井系統(tǒng)開(kāi)采參數(shù)[39]

表5 Plackett-Burman設(shè)計(jì)方案[39]

根據(jù)對(duì)井系統(tǒng)I的計(jì)算結(jié)果(表5),通過(guò)多元回歸得到熱儲(chǔ)響應(yīng)與開(kāi)采參數(shù)之間的多項(xiàng)式代理模型:

ER=8.59+0.001 172Q-0.293 5T+0.004 689D+ 0.000 243d-0.012 74α(R2=0.93)

(2)

(3)

(4)

(5)

式中,ER為可采地?zé)崮?LW為水位;DV為垂直位移;TBT為熱突破時(shí)間。

決定系數(shù)R2在0.88～0.99,說(shuō)明所建立的代理模型可以較好地預(yù)測(cè)熱儲(chǔ)響應(yīng)與開(kāi)采參數(shù)之間的對(duì)應(yīng)關(guān)系。圖10給出了對(duì)井系統(tǒng)I的最優(yōu)開(kāi)采參數(shù)。提高流量可增大可采地?zé)崮?、垂直位移和水位埋?但會(huì)縮短熱突破時(shí)間,故最優(yōu)流量(2 181 m3/d)在取值范圍中間。降低回灌溫度均會(huì)使4個(gè)熱儲(chǔ)響應(yīng)參數(shù)增大,但與可采地?zé)崮芎痛怪蔽灰葡啾?回灌溫度對(duì)熱突破時(shí)間和水位埋深的影響不太顯著,故最優(yōu)回灌溫度為下限25 ℃。增大井距可增大可采地?zé)崮?、熱突破時(shí)間和水位埋深,但會(huì)降低垂直位移,故最優(yōu)井距(677 m)在取值范圍中間。增大井深會(huì)略微增大可采地?zé)崮?、垂直位移和水位埋?但會(huì)縮短熱突破時(shí)間?？紤]到增加井深會(huì)增加鉆井成本,故最優(yōu)井深為下限2 357 m。增大回灌井旋轉(zhuǎn)角可減小可采地?zé)崮?、熱突破時(shí)間和水位埋深,但會(huì)增加垂直位移。不過(guò),與可采地?zé)崮芎蜔嵬黄茣r(shí)間相比,增大回灌井旋轉(zhuǎn)角對(duì)垂直位移和水位埋深的影響較小,故回灌井最優(yōu)方位應(yīng)選在滲流方向的下游。

將上述最優(yōu)開(kāi)采參數(shù)作為城市副中心井-儲(chǔ)數(shù)值模型的輸入?yún)?shù),采用COMSOL計(jì)算熱儲(chǔ)響應(yīng),進(jìn)一步檢驗(yàn)代理模型的準(zhǔn)確性。結(jié)果表明,由于采灌不均衡,通州副中心出現(xiàn)了大范圍的抽水漏斗,最大水位埋深超過(guò)了200 m,無(wú)法滿(mǎn)足可持續(xù)開(kāi)發(fā)利用的要求。而地?zé)釋?duì)井系統(tǒng)I的水位、熱突破和垂直位移均滿(mǎn)足可持續(xù)開(kāi)發(fā)利用的評(píng)價(jià)標(biāo)準(zhǔn)。代理模型與數(shù)值模型的對(duì)比結(jié)果表明,4個(gè)熱儲(chǔ)響應(yīng)的最大計(jì)算誤差為11.6%,說(shuō)明代理模型能夠較準(zhǔn)確地表征熱儲(chǔ)響應(yīng)和開(kāi)采參數(shù)之間的關(guān)系。與Plackett-Burman設(shè)計(jì)的12次運(yùn)行結(jié)果相比,采用優(yōu)化后的開(kāi)采參數(shù)得到了滿(mǎn)足可持續(xù)開(kāi)采標(biāo)準(zhǔn)的最大可采地?zé)崮?圖11)。

采用類(lèi)似的流程,得到了對(duì)井系統(tǒng)II的代理模型及其最優(yōu)開(kāi)采參數(shù)(圖10)。相比于對(duì)井系統(tǒng)I,對(duì)井系統(tǒng)II的最優(yōu)流量和井深較小,而最優(yōu)井距較大?？赡艿脑蚴菍?duì)井系統(tǒng)II的平均滲透系數(shù)(0.38 m/d)小于對(duì)井系統(tǒng)I(0.5 m/d),需采用較小的流量確保水位埋深不超過(guò)150 m。因此,由于儲(chǔ)層參數(shù)的非均質(zhì)性和地下流場(chǎng)的各向異性,同一地?zé)崽镏胁煌責(zé)衢_(kāi)采區(qū)內(nèi)對(duì)井系統(tǒng)的最優(yōu)人為控制參數(shù)可能會(huì)有明顯的不同。

圖11 可采地?zé)崮芘c開(kāi)采風(fēng)險(xiǎn)等級(jí)[39]Fig.11 Recoverable geothermal energy and production risk rating[39]

5 結(jié) 論

(1)從溫度、水位、變形3個(gè)方面給出了深層地?zé)崮芸沙掷m(xù)開(kāi)采的多指標(biāo)評(píng)價(jià)體系,以及基于數(shù)值模型的可采地?zé)崮苡?jì)算公式。但是,各評(píng)價(jià)指標(biāo)的允許值應(yīng)根據(jù)具體情況合理制定。

(2)基于深層地?zé)峋?儲(chǔ)系統(tǒng)多場(chǎng)耦合數(shù)值模型給出了滿(mǎn)足可持續(xù)開(kāi)采標(biāo)準(zhǔn)的可采地?zé)崮茉u(píng)價(jià)方法,提出了開(kāi)采井風(fēng)險(xiǎn)評(píng)價(jià)的“紅綠燈系統(tǒng)”。

(3)基于深層地?zé)峋?儲(chǔ)系統(tǒng)多場(chǎng)耦合數(shù)值模型給出了滿(mǎn)足可持續(xù)開(kāi)采標(biāo)準(zhǔn)的采灌方案優(yōu)化設(shè)計(jì)方法,使用Plackett-Burman設(shè)計(jì)構(gòu)建了代理模型,通過(guò)多響應(yīng)優(yōu)化實(shí)現(xiàn)了參數(shù)的快速優(yōu)化。

致謝研究過(guò)程中得到了中國(guó)地質(zhì)科學(xué)院水文地質(zhì)環(huán)境地質(zhì)研究所王貴玲研究員和馬峰研究員級(jí)高工、北京市地?zé)嵴{(diào)查研究所張進(jìn)平正高工等專(zhuān)家的指導(dǎo)和幫助,在此衷心感謝!