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

譚志容，高宗軍，趙季初

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

?

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

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

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

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

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

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

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

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

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

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

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

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

1.2 館陶組熱儲特征

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

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

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

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

2.1 概念模型

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

圖3 熱儲概念模型圖

2.2 無限含水層邊界的概化

2.2.1 邊界條件的分類

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

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

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

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

2.2.2 模型邊界的選擇

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

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

2.3 模型的建立與識別

2.3.1 單元剖分

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

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

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

2.3.2 抽水井數(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é)峋募虚_采時間段為11月15日～至翌年3月15日，約120 d，單井開采量約80 m3/h；其余時間段開采量極小，在模型中可以忽略。

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

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

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

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

T=K·M

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

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

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

μe=Ss·M

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

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

2)邊界的輸入

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

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

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

2.4 模型的校正

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

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

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

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

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

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

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

圖4 模型校正前計算值與實測值校正圖

圖5 DR14-1、DR34-1兩眼觀測孔水位計算值與觀測值擬合曲線

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

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

4 結(jié)語

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

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

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

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

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

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

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

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

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

[6]童彥釗, 吳平, 韓強強. 地下水數(shù)值模擬軟件簡述[J]. 科技信息.2012, (29):48-48.

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

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

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

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

[11]趙季初,譚志容,王成明.德城區(qū)館陶組熱儲水文地質(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-)，女，湖北荊州人，高級工程師，主要從事水文地質(zhì)方面工作。

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

P314.1

A

1004-1184(2017)03-0019-03