譚志容，高宗軍，趙季初

(1.山東省魯北地質(zhì)工程勘察院，山東 德州 253072；2.山東科技大學(xué)地科學(xué)院，山東 青島 266590)

?

基于Modflow軟件的德州市城區(qū)館陶組熱儲(chǔ)數(shù)值模擬研究

譚志容1，高宗軍2，趙季初1

(1.山東省魯北地質(zhì)工程勘察院，山東 德州 253072；2.山東科技大學(xué)地科學(xué)院，山東 青島 266590)

德州市位于山東省西北部，區(qū)內(nèi)地?zé)豳Y源豐富，城區(qū)已有地?zé)峋?5眼，地?zé)豳Y源大量開發(fā)，對(duì)地?zé)豳Y源的管理提出了新的要求，本文以德州市城區(qū)為研究區(qū)，充分收集資料，采用Modflow軟件建立德州市城區(qū)館陶組熱儲(chǔ)模型，對(duì)館陶組熱儲(chǔ)的埋藏分布規(guī)律及城區(qū)地?zé)崽锏牡責(zé)崴畠?chǔ)量進(jìn)行計(jì)算評(píng)價(jià)，將館陶組熱儲(chǔ)概化為平面上無限延伸且頂、底板為隔水層的均質(zhì)含水層。采用非穩(wěn)定流抽水試驗(yàn)求取的水文地質(zhì)參數(shù)對(duì)模型進(jìn)行識(shí)別，利用兩眼長(zhǎng)期監(jiān)測(cè)孔的水位動(dòng)態(tài)觀測(cè)數(shù)據(jù)及歷年來地?zé)豳Y源的開采情況對(duì)模型進(jìn)行校正。校正過程中發(fā)現(xiàn)與熱儲(chǔ)砂巖互層分布的泥巖壓密釋水是熱儲(chǔ)地?zé)崴闹饕獊碓粗弧Ｒ缘責(zé)崴唤捣? m/a為限制條件，德城區(qū)館陶組熱儲(chǔ)地?zé)崴脑试S開采量為532萬 m3/a。

德州市城區(qū)；館陶組熱儲(chǔ)；數(shù)值模擬；通用水頭邊界；允許開采量

德州市位于山東省西北部，在大地構(gòu)造單元上處于華北平原的東南部。自1996年在德州市城區(qū)開展地?zé)豳Y源普查、施工了兩眼探采結(jié)合孔以來，區(qū)內(nèi)地?zé)豳Y源開發(fā)蓬勃發(fā)展，并帶動(dòng)了周邊地區(qū)地?zé)豳Y源的開發(fā)。目前德州市各縣級(jí)城區(qū)均有數(shù)眼至數(shù)十眼地?zé)峋坏?，德州市城區(qū)已有地?zé)峋?2眼。

地?zé)豳Y源的大量開發(fā)，對(duì)地?zé)豳Y源的管理工作提出了新的要求，其中地?zé)崴刹少Y源量的確定是整個(gè)地?zé)豳Y源管理工作的核心。根據(jù)規(guī)范要求：“在地?zé)崽锏目辈槌潭缺容^高，并且具有一定時(shí)間的開采歷史，具有比較齊全的監(jiān)測(cè)資料時(shí)，應(yīng)建立地?zé)崽锏臄?shù)值模擬模型，用以計(jì)算/評(píng)價(jià)地?zé)醿?chǔ)量，并作為地?zé)崽锕芾淼墓ぞ摺盵1]。

本文在充分收集整理以往勘探成果的基礎(chǔ)上，借鑒前人在地下水?dāng)?shù)值模擬方面取得的成功經(jīng)驗(yàn)[2-6]，采用Modflow軟件建立德州市城區(qū)館陶組熱儲(chǔ)的數(shù)值模型，對(duì)德州市城區(qū)地?zé)崽锏牡責(zé)崴畠?chǔ)量進(jìn)行計(jì)算與評(píng)價(jià)。

1 地?zé)岬刭|(zhì)條件

1.1 地質(zhì)構(gòu)造

德州市城區(qū)在大地構(gòu)造單元上位于臨清坳陷的東北部，東鄰埕寧隆起、西接滄縣隆起，北與黃驊坳陷相通，為兩隆夾一坳的地塹型構(gòu)造(圖1)。其內(nèi)沉積了巨厚的新生代地層，是地下熱水的賦存場(chǎng)所，其中館陶組下段砂巖是區(qū)內(nèi)地?zé)豳Y源開發(fā)的目標(biāo)熱儲(chǔ)[7-8]。

圖1 模擬區(qū)構(gòu)造位置綱要圖

1.2 館陶組熱儲(chǔ)特征

館陶組熱儲(chǔ)在整個(gè)華北平原內(nèi)皆有分布[9]，具有下粗上細(xì)的正旋回沉積特征，底部砂礫巖層分布較穩(wěn)定。平面上熱儲(chǔ)呈層狀連續(xù)分布，垂向上熱儲(chǔ)與泥巖呈互層分布。熱儲(chǔ)巖性主要為河流相、沖積扇相的細(xì)砂巖、粗砂巖、含礫砂巖與砂礫巖，礫石呈半圓狀，磨圓度中等，礦物成分主要為石英與長(zhǎng)石。

模擬區(qū)館陶組熱儲(chǔ)的厚度與埋深受基底構(gòu)造的影響，底板埋深1 300～1 600 m，熱儲(chǔ)厚120～200 m(圖2)，單井涌水量80～100 m3/h，水化學(xué)類型為Cl—Na型，礦化度5 g/L左右。

圖2 模擬區(qū)館陶組熱儲(chǔ)厚度與底板埋深圖

2 基于Modflow軟件的數(shù)值模型

2.1 概念模型

根據(jù)區(qū)內(nèi)館陶組熱儲(chǔ)特征，為降低模型的復(fù)雜程度，本次將整個(gè)研究區(qū)概化為一個(gè)均質(zhì)各向同性的穩(wěn)定含水層。熱儲(chǔ)頂?shù)装寰植加泻穸染薮蟮哪鄮r層，在模擬期內(nèi)可視為隔水層。熱儲(chǔ)在整個(gè)華北平原內(nèi)均有分布，相對(duì)于整個(gè)華北平原而言，可將熱儲(chǔ)概化為平面上無限延伸的含水層，其概念模型如圖3所示。

圖3 熱儲(chǔ)概念模型圖

2.2 無限含水層邊界的概化

2.2.1 邊界條件的分類

地下水流動(dòng)問題的邊界條件主要分為給定水頭邊界條件(第一類邊界)、給定流量邊界條件(第二類邊界)與潛水面邊界(陳崇希等.1999)[10]。

Modeflow 軟件將模擬區(qū)劃分為活動(dòng)單元(Active Cells)與非活動(dòng)單元(Inactive Cells)，非活動(dòng)單元與活動(dòng)單元之間的邊界可認(rèn)為是零流量邊界(特殊的第二類邊界)。此外，該軟件還給出了8種不同的含水層邊界條件，其中適宜本概念模型邊界條件的有兩種，即定水頭邊界(Constant Head)與通用水頭邊界(General Head)。

定水頭邊界(第一類邊界)可通過邊界附近的觀測(cè)井監(jiān)測(cè)資料，給出不同時(shí)間段的初始時(shí)水位與結(jié)束時(shí)水位，軟件自動(dòng)插值確定各模擬期邊界的水位，當(dāng)模擬區(qū)邊界附近有較多的監(jiān)測(cè)井時(shí)，這類邊界條件能較真實(shí)地反映模擬區(qū)與外界的水量交換。

通用水頭邊界(第二類邊界)通過計(jì)算該邊界處活動(dòng)單元的水位與模擬區(qū)外一定距離處給定的水位之間的差值，計(jì)算流入或流出模擬區(qū)的水量。

2.2.2 模型邊界的選擇

根據(jù)概念模型，熱儲(chǔ)概化為四周向外無限延伸的均質(zhì)各向同性含水層，最為適宜的邊界條件應(yīng)該采用給定水頭邊界，但由于模擬區(qū)邊界附近無監(jiān)測(cè)井分布，沒有辦法給定各模擬時(shí)間段熱儲(chǔ)的水位。若采用零流量邊界來定義模擬區(qū)四周的邊界，則模型的范圍要外擴(kuò)，至少達(dá)到模擬區(qū)抽水井的影響范圍之外，采用這種方法不能重點(diǎn)突出模擬區(qū)的水位變動(dòng)情況。

因此，本次采用通用水頭邊界來刻畫這種無限延伸的含水層邊界條件，即假定離模擬區(qū)邊界一定距離處的水位保持不變，這種假定也與實(shí)際情況基本相符。

2.3 模型的建立與識(shí)別

2.3.1 單元剖分

本次以德州市城區(qū)的地理底圖作為模型建立的底圖，將模擬區(qū)概化為一個(gè)南北長(zhǎng)37 453 m，東西寬35758m的矩形區(qū)域。

根據(jù)地?zé)峋你@探與成井資料，垂向上將整個(gè)模擬區(qū)概化為深度0～1 300 m、1 300～1 550 m及1 550～2 000 m三個(gè)地質(zhì)層。其中1 300～1 550 m為區(qū)內(nèi)地?zé)峋娜∷畬?，即熱?chǔ)層，屬館陶組下段，地層巖性為砂礫巖、粗砂巖及中粗砂巖與泥巖的互層，其中砂巖占地層厚度的比例60%左右。0～1 300 m為熱儲(chǔ)蓋層，從上至下為第四系、明化鎮(zhèn)組及館陶組上段；1 550～2 000 m為熱儲(chǔ)下伏層，從上至下為東營(yíng)組與沙河街組。

平面上將整個(gè)模擬區(qū)劃分為100行×100列的均勻網(wǎng)。

2.3.2 抽水井?dāng)?shù)據(jù)處理

德州市城區(qū)地?zé)豳Y源開采始于1998年，截止目前共有地?zé)峋?2眼，其中1998年投入使用2眼，2002年新增3眼；2003年新增7眼，2004年新增4眼；2005年新增11眼；2006年新增24眼；2007年新增3眼；2008年新增8眼；2009年新增6眼；2010年新增1眼；2011新增3眼。

各地?zé)峋募虚_采時(shí)間段為11月15日～至翌年3月15日，約120 d，單井開采量約80 m3/h；其余時(shí)間段開采量極小，在模型中可以忽略。

因此，各井自投入使用開始，每年開采應(yīng)力期概化為兩段：即每年的11月15日～至翌年3月15日，單井開采量1 920 m3/d；每年的3月16日～11月14日開采量為零。以1998年1月1日為模擬的初始時(shí)間，以2016年2月21日為模擬的結(jié)束時(shí)間，總時(shí)長(zhǎng)為6 625 d。

2.3.3 水文地質(zhì)參數(shù)及邊界

1) 滲透系數(shù)(K)與貯水率(Ss)

根據(jù)非穩(wěn)定流抽水試驗(yàn)結(jié)果，區(qū)內(nèi)館陶熱儲(chǔ)的導(dǎo)水系數(shù)為944.5 m2/d，貯水系數(shù)(彈性釋水系數(shù))為2.48×10-4[11]。本次模型的熱儲(chǔ)厚度取值250 m，根據(jù)導(dǎo)水系數(shù)公式：

T=K·M

式中：T為導(dǎo)水系數(shù)(m2/d)；K為滲透系數(shù)(m/d)；M為熱儲(chǔ)厚度(m)。

經(jīng)計(jì)算可得，熱儲(chǔ)的平均滲透系數(shù)為3.778 m/d。

根據(jù)貯水率(Ss)的定義，當(dāng)含水層水頭變化一個(gè)單位時(shí)，從單位體積含水層中，由于水體積膨脹(或壓縮)以及介質(zhì)骨架的壓縮(或伸長(zhǎng))而釋放(或貯存)的彈性水量[12]。根據(jù)貯水系數(shù)公式：

μe=Ss·M

式中：μe為貯水系數(shù)(無量綱)；Ss為貯水率(1/m)。

經(jīng)計(jì)算，熱儲(chǔ)的平均貯水率為9.992×10-7(1/m)。

2)邊界的輸入

本次假定在抽水影響范圍以外的熱儲(chǔ)保持定水頭不變，其關(guān)鍵是影響范圍的確定，采用裘布依公式進(jìn)行估算[13]：

式中ΔS為距開采中心R處與開采中心的熱儲(chǔ)水位差(m)；Q為地?zé)崴_采量，本次取72眼井的共同開采量(m3/d)；R為影響半徑(m)；r為72眼井圈定的假想圓半徑(m)；

根據(jù)模擬區(qū)兩觀測(cè)孔的最大水位降約68 m左右，推算的R值為135 km左右，即將通用水頭在離模擬區(qū)外135 km處的水位設(shè)為初始值“零”。

2.4 模型的校正

利用模擬區(qū)兩眼觀測(cè)孔長(zhǎng)期水位動(dòng)態(tài)監(jiān)測(cè)數(shù)據(jù)對(duì)模型進(jìn)行校正。兩眼監(jiān)測(cè)孔的編號(hào)分別為DR14-1與DR34-1。運(yùn)行模型后，結(jié)果顯示，計(jì)算值與觀測(cè)值的擬合度差，模型參數(shù)不合理，需要校正。

通過對(duì)各水文地質(zhì)參數(shù)進(jìn)行分析可知，由于抽水試驗(yàn)取得的彈性系數(shù)主要反應(yīng)含水層的貯水/釋水能力，而實(shí)際情況是熱儲(chǔ)為砂、泥巖互層，在長(zhǎng)期抽水作用下，泥巖的壓密釋水不容忽視，因此貯水率值要增大；滲透系數(shù)主要是反映熱儲(chǔ)的滲透能力，區(qū)內(nèi)熱儲(chǔ)較厚，泥巖對(duì)熱儲(chǔ)滲透性能的影響程度略輕，因此熱儲(chǔ)平均滲透系數(shù)也要適當(dāng)增大；另外，受復(fù)雜的地質(zhì)條件影響，側(cè)向逕流補(bǔ)給也沒有想象中那么大，在通用水頭邊界賦值時(shí)，在不改變距離R值時(shí)應(yīng)將模擬區(qū)外熱儲(chǔ)的滲透系數(shù)減小(相當(dāng)于增加了R的距離)。

根據(jù)以上分析，經(jīng)反復(fù)調(diào)試，最終確定校正前后的水文地質(zhì)參數(shù)如下表所示(表1)：

表1 校正前后水文地質(zhì)參數(shù)對(duì)比表

注：K’為通用水頭邊界賦值的模擬區(qū)外熱儲(chǔ)滲透系數(shù)

經(jīng)參數(shù)校正后，檢查校正結(jié)果圖(Calibration graph)發(fā)現(xiàn)，計(jì)算值=觀測(cè)值的直線位于95%置信區(qū)間的中心，絕大部分觀測(cè)值位于95%區(qū)間內(nèi)(圖4)，說明模型已校正。

另外，兩個(gè)觀測(cè)井的水位降深觀測(cè)值與計(jì)算值的擬合度較高(圖5)，符合《地下水資源管理模型工作要求》(GB/T 14497-93)中關(guān)于“水位擬合時(shí)相對(duì)誤差小于10%的節(jié)點(diǎn)必須占已知水位節(jié)點(diǎn)的70%以上”[14]的規(guī)定。

圖4 模型校正前計(jì)算值與實(shí)測(cè)值校正圖

圖5 DR14-1、DR34-1兩眼觀測(cè)孔水位計(jì)算值與觀測(cè)值擬合曲線

3 地?zé)崴试S開采量計(jì)算結(jié)果

經(jīng)校正后，所建立的數(shù)值計(jì)算模型具有較好的準(zhǔn)確度。利用該模型既可用來預(yù)測(cè)不同開采條件下的水位變動(dòng)情況，亦可通過控制熱儲(chǔ)水位的年降幅，求得該允許的水位年降幅條件下的開采量作為地?zé)崴目刹少Y源量。參照規(guī)范[1]，以年壓力下降速度不大于0.02 Mpa(相當(dāng)于2 m/a)作為限制條件，運(yùn)行該模型，求得德州市城區(qū)內(nèi)館陶組熱儲(chǔ)地?zé)崴目刹少Y源量為532萬 m3/a。

4 結(jié)語(yǔ)

(1)德州市城區(qū)館陶組熱儲(chǔ)地?zé)崴_發(fā)始于1998年，目前已有地?zé)崴?2眼。地?zé)崴畽C(jī)制開采期是每年11月15日至來年的3月15日，單井開采量為80 m3/h。

(2)由區(qū)域地質(zhì)條件，德州市城區(qū)館陶組熱儲(chǔ)可概化為厚度250 m、頂?shù)装鍨楦羲畬拥钠矫鏌o限延伸含水層，模型采用了通用水頭邊界。

(3)分析研究發(fā)現(xiàn)，開采條件下熱儲(chǔ)的貯水率大于單井抽水試驗(yàn)時(shí)求取的貯水率，說明泥巖壓密釋水是熱儲(chǔ)地?zé)崴闹匾獊碓?，故貯水率上調(diào)了1個(gè)數(shù)量級(jí)。由此可見，長(zhǎng)期地?zé)崴_采，勢(shì)必造成不可恢復(fù)的地面沉降，應(yīng)引起重視。

(4)以地?zé)崴唤捣? m/a為限制條件，德州市城區(qū)館陶組熱儲(chǔ)地?zé)崴刹少Y源量為532萬 m3/a。

[1]GB/T 11615-2010.地?zé)豳Y源地質(zhì)勘查規(guī)范[S].

[2]薛禹群.中國(guó)地下水?dāng)?shù)值模擬的現(xiàn)狀與展望[J].高校地質(zhì)學(xué)報(bào).2010,16(1):1-6.

[3]孫從軍,韓振波,趙振，等.地下水?dāng)?shù)值模擬的研究與應(yīng)用進(jìn)展[J].環(huán)境工程.2013,31(5):9-13.

[4]武強(qiáng)，朱斌，徐華，等.MODFLOW 在淮北地下水?dāng)?shù)值模擬中的應(yīng)用[J].遼寧工程技術(shù)大學(xué)學(xué)報(bào).2005(4):503.

[5]沈媛媛,蔣云鐘,雷曉輝，等.地下水?dāng)?shù)值模擬中人為邊界的處理方法研究[J].水文地質(zhì)工程地質(zhì).2008(6).

[6]童彥釗, 吳平, 韓強(qiáng)強(qiáng). 地下水?dāng)?shù)值模擬軟件簡(jiǎn)述[J]. 科技信息.2012, (29):48-48.

[7]顏世強(qiáng),潘懋,鄒祖光，等.山東德州凹陷地下熱水地球化學(xué)特征及成因[J].中國(guó)地質(zhì).2007,34(1):149-152.

[8]趙季初.山東省德州市城區(qū)館陶組熱儲(chǔ)地?zé)豳Y源評(píng)價(jià)[D].北京:中國(guó)地質(zhì)大學(xué)(北京).2007.

[9]陳墨香.華北地?zé)醄M].北京：科學(xué)出版社.1988，90-91.

[10]陳崇希,林敏.地下水動(dòng)力學(xué)[M].武漢:中國(guó)地質(zhì)大學(xué)出版社(武漢市洪山區(qū)魯磨路31號(hào)).1999:32.

[11]趙季初,譚志容,王成明.德城區(qū)館陶組熱儲(chǔ)水文地質(zhì)參數(shù)探討[J].地?zé)崮?2008(2):16-18.

[12]張惠昌. 關(guān)于貯水系數(shù)與貯水率意義的商榷[J]. 勘察科學(xué)技術(shù).1987(4).

[13]JGJ/T111-98，建筑與市政降水工程技術(shù)規(guī)范[S].23-24.

[14]GB/T 14497-1993，地下水資源管理模型工作要求[S].

Study on Numerical Simulation of thermal reservoir of Guantao formation in Dezhou based on Modflow software

TAN Zhi-rong1，GAO Zong-Jun2，Zhao Ji-chu1

(1.Shandong Lubei Geological Engineering Investigation Institute, Dezhou, Shandong, 253072, China；2.Academy of Sciences of Shandong University of Science and Technology, Qingdao, Shandong, 266590, China)

Shandong province is located in the northwest region, rich in geothermal resources, the city has 75 geothermal wells, geothermal resources development, put forward new requirements for the management of geothermal resources, this paper takes Dezhou city as the study area, to collect sufficient data, using Modflow software to establish the Dezhou City geothermal reservoir of Guantao formation model. Water distribution and urban buried geothermal field of reservoir of Guantao formation in the evaluation calculation, the reservoir of Guantao formation is generalized to plane and infinite top floor aquifuge homogeneous aquifer. The unsteady flow pumping test of hydro geological parameters obtained from the test is used to identify the model, with both eyes the long-term monitoring of water level dynamic observation data of holes and over the years, exploitation of geothermal resources of the model. The sandstone and reservoir interblended mudstone compaction released water is the main reservoir of geothermal water found in the process of correction One of the sources is to reduce the geothermal water level of 2m/a as the limiting condition, the allowable exploitation of geothermal water in the Guantao formation of area is 5 million 320 thousand

Dezhou city; guantao reservoir; numerical simulation; general head boundary; allowable exploitation

2017-02-20

譚志容(1975-)，女，湖北荊州人，高級(jí)工程師，主要從事水文地質(zhì)方面工作。

高宗軍(1964-)，男，山東泰安人，教授，研究方向：水文地質(zhì)。

P314.1

A

1004-1184(2017)03-0019-03