MODFLOW在天津?yàn)I海新區(qū)地?zé)崃黧w數(shù)值模擬中的應(yīng)用

張芬娜，賈志，2，高亮，2，劉九龍

（1.天津地?zé)峥辈殚_(kāi)發(fā)設(shè)計(jì)院，天津 300250；2.中國(guó)地質(zhì)大學(xué)（北京），北京 100083）

MODFLOW在天津?yàn)I海新區(qū)地?zé)崃黧w數(shù)值模擬中的應(yīng)用

張芬娜1，賈志1，2，高亮1，2，劉九龍1

（1.天津地?zé)峥辈殚_(kāi)發(fā)設(shè)計(jì)院，天津 300250；2.中國(guó)地質(zhì)大學(xué)（北京），北京 100083）

為了合理開(kāi)發(fā)利用地?zé)豳Y源，本文以天津?yàn)I海新區(qū)橋沽-看財(cái)莊新近系明化鎮(zhèn)組地?zé)崃黧w為例，應(yīng)用Visual MODFLOW軟件對(duì)明化鎮(zhèn)組地?zé)崃黧w進(jìn)行數(shù)值模擬。對(duì)模型的識(shí)別與檢驗(yàn)表明，所建模型能夠較好的反映明化鎮(zhèn)組地?zé)崃黧w的地下水位動(dòng)態(tài)特征。計(jì)算得出明化鎮(zhèn)組100年水位接近150 m的地?zé)崃黧w可采資源量為10.7×106m3/a。通過(guò)預(yù)測(cè)按設(shè)計(jì)開(kāi)采量進(jìn)行開(kāi)采后明化鎮(zhèn)組地?zé)崃黧w地下水位的變化趨勢(shì)，每天新增1000 m3開(kāi)采量是可行的。

天津?yàn)I海新區(qū)；地?zé)崃黧w；數(shù)值模擬；MODFLOW

近幾十年來(lái)，地下水模型已被廣泛應(yīng)用于地下水資源評(píng)價(jià)、預(yù)報(bào)和管理之中[1]。目前國(guó)內(nèi)常用的模擬地下水位動(dòng)態(tài)預(yù)測(cè)的軟件有Visual MODFLOW、FEFLOW和GMS[2]等。有許多學(xué)者建立了MODFLOW或者FEFLOW地下水流模型，用水文觀測(cè)井的地下水動(dòng)態(tài)觀測(cè)資料，對(duì)數(shù)值模型進(jìn)行了識(shí)別和校驗(yàn)后，用于預(yù)測(cè)不同開(kāi)采方案的地下水動(dòng)態(tài)變化[3-8]。天津市地?zé)崃黧w地下水位模擬研究成果有胡燕建立的天津市濱海新區(qū)新近系館陶組地?zé)崃黧w流場(chǎng)三維數(shù)值模型[9]和田光輝建立的灰色系統(tǒng)GM(1,1)模型[10]等。但在地?zé)崃黧w地下水位動(dòng)態(tài)預(yù)測(cè)的模擬中利用Visual MODFLOW和FEFLOW軟件的較少，Visual MODFLOW軟件在計(jì)算模型的建立、計(jì)算結(jié)果的可視化、數(shù)據(jù)前后處理以及與其他軟件數(shù)據(jù)信息交互等方面能力較好[11-14]。因此，本文運(yùn)用Visual MODFLOW建立明化鎮(zhèn)組地?zé)崃黧w三維水流模擬模型，計(jì)算明化鎮(zhèn)組的可采資源量，并通過(guò)模擬地?zé)崃黧w地下水位動(dòng)態(tài)，來(lái)判斷增加設(shè)計(jì)開(kāi)采量是否可行。

1 研究區(qū)概況

研究區(qū)地處天津?yàn)I海新區(qū)北部橋沽－看財(cái)莊一帶，位于Ⅲ級(jí)構(gòu)造單元黃驊坳陷的北部，橫跨寧河凸起和北塘凹陷兩個(gè)Ⅳ級(jí)構(gòu)造單元[15-16]，區(qū)內(nèi)及其附近主要發(fā)育的斷裂有：漢沽斷裂、寧河斷裂、大神堂斷裂、茶淀斷裂和北塘東斷裂（圖1）。在當(dāng)今經(jīng)濟(jì)技術(shù)條件下區(qū)內(nèi)可開(kāi)采的熱儲(chǔ)層有新近系明化鎮(zhèn)組、館陶組和古近系東營(yíng)組[17]。本文主要以新近系明化鎮(zhèn)組作為研究目的層。

2 熱儲(chǔ)特征

圖1 研究區(qū)構(gòu)造簡(jiǎn)圖Fig.1 Sketch map of the structures

明化鎮(zhèn)組作為一套曲流相碎屑巖，在全區(qū)均有

分布，儲(chǔ)集條件良好，厚度發(fā)育穩(wěn)定，具有北部寧河凸起淺，南部北塘凹陷深的特點(diǎn)[18]。頂板埋深376～506 m，砂巖占總厚度的25.7%～29%，單井涌水量54～85 m3/h，液面溫度25～53℃，總礦化度439.1～584.3 mg/L，pH值8.15～8.69，水質(zhì)類型以HCO3-Na型為主[19]，局部地區(qū)為HCO3·SO4-Na型水[20]。2012年明化鎮(zhèn)組年開(kāi)采量69×104m3，水位58～62 m，年降幅1 m左右[21]。

3 概念模型

計(jì)算范圍：為了在地下水流數(shù)值模擬過(guò)程中更加準(zhǔn)確地給定邊界水位（或流量），在確定模擬區(qū)范圍時(shí)，將模擬邊界適當(dāng)向外擴(kuò)展到具有長(zhǎng)觀孔游樂(lè)港1（觀3）的區(qū)域，將東南部邊界擴(kuò)展到海岸線，整個(gè)數(shù)值模擬區(qū)域總面積為423.4 km2（圖1）。

邊界條件設(shè)定：把明化鎮(zhèn)組頂?shù)装宥楦羲吔纭Ｍㄟ^(guò)本區(qū)的地?zé)岬刭|(zhì)、水文地質(zhì)條件和地?zé)崃黧w開(kāi)發(fā)利用特點(diǎn)、同位素分析結(jié)果可知，計(jì)算區(qū)新近系明化鎮(zhèn)組熱儲(chǔ)層為層狀結(jié)構(gòu)，分布范圍廣闊，不是完整的水文地質(zhì)單元。因此，側(cè)向邊界設(shè)為自然流量邊界（通用水頭邊界）。

4 數(shù)學(xué)模型的建立

模擬區(qū)內(nèi)水文地質(zhì)參數(shù)隨空間變化，參數(shù)在水平方向沒(méi)有明顯的方向性，因此地?zé)崃黧w模型可概化為三維、非均質(zhì)水平各向同性的非穩(wěn)定熱流體系統(tǒng)，在笛卡爾坐標(biāo)系下可描述為[22]：

式中：Ω為滲流區(qū)域；h為熱流體的水位標(biāo)高(m)；Kx、Ky、Kz分別為熱流體儲(chǔ)積層中x、y、z方向的滲透系數(shù)(m/d)；Kn為邊界面法向方向熱流體儲(chǔ)積層滲透系數(shù)(m/d)；Kz為邊界內(nèi)側(cè)熱流體儲(chǔ)積層垂向滲透系數(shù)(m/d)；μs為單位儲(chǔ)水率（1/m）；σ為通用水頭邊界的阻力系數(shù)，σ=L/K；L為模型邊界到通用水頭邊界的水平距離(m)；K為模型邊界到通用水頭邊界之間熱流體儲(chǔ)積層平均滲透系數(shù)(m/d)；ε為源匯項(xiàng)(m/ d)；h0為初始水位分布(m)，h0=h0(x,y,z)；hb為通用水頭邊界上的水位分布(m)；n為邊界面的法線方向；z→為熱儲(chǔ)層頂、底邊界面的法線方向；Γ1為通用水頭邊界；Γ2為滲流區(qū)域頂、底邊界。

5 模型數(shù)值模擬

5.1 時(shí)空離散

模擬區(qū)面積約514 km2。利用Visual Modflow程序?qū)γ骰?zhèn)組水位動(dòng)態(tài)進(jìn)行模擬，采用矩形等距剖分，剖分成80×60個(gè)單元格，其中活動(dòng)單元格3955個(gè)，非活動(dòng)單元格845個(gè)。有效活動(dòng)單元計(jì)算區(qū)面積423.4 km2，有效活動(dòng)單元計(jì)算區(qū)即為本次工作的模擬計(jì)算區(qū)。

5.2 數(shù)據(jù)處理

以2012年4月明化鎮(zhèn)組熱儲(chǔ)層動(dòng)態(tài)長(zhǎng)觀資料為基礎(chǔ)，對(duì)地?zé)崃黧w的溫度、密度統(tǒng)一校正（熱儲(chǔ)層溫度統(tǒng)一校正到40℃），利用Kriging優(yōu)化插值法取得初始水位，通過(guò)模擬運(yùn)算，修改流場(chǎng)，得到該層的初始流場(chǎng)。

源匯項(xiàng)包括側(cè)向徑流補(bǔ)給、層間越流補(bǔ)給、人工回灌量、人工開(kāi)采等。將源匯項(xiàng)輸入時(shí)，有具體井位的按坐標(biāo)輸入，其它單元的源匯項(xiàng)概化為井的形式，平均賦于區(qū)內(nèi)的各個(gè)網(wǎng)格點(diǎn)中。其中層間越流量由Visual MODFLOW模型軟件根據(jù)所賦相關(guān)參數(shù)自動(dòng)進(jìn)行處理。

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

本次模型識(shí)別驗(yàn)證采用試估－校正法，即反復(fù)調(diào)整參數(shù)，達(dá)到較為理想的擬合結(jié)果。本次模型識(shí)別時(shí)段選擇2012年4月1日－2013年10月1日，明化鎮(zhèn)組參與擬合的長(zhǎng)觀井共4眼（圖1）?？蛇M(jìn)行地?zé)崃黧w水位歷時(shí)曲線和區(qū)域流場(chǎng)的擬合（圖2），計(jì)算出來(lái)的長(zhǎng)觀井地?zé)崃黧w水位動(dòng)態(tài)變化與實(shí)際變化趨勢(shì)一致，熱流體流場(chǎng)擬合等值線與實(shí)測(cè)流場(chǎng)吻合較好（圖3）。根據(jù)明化鎮(zhèn)組熱儲(chǔ)層的埋藏分布和控?zé)針?gòu)造分布特征分為3區(qū)，各區(qū)的水文地質(zhì)參數(shù)主要有：滲透系數(shù)K和彈性釋水系數(shù)μ*。各區(qū)的水文地質(zhì)參數(shù)初值主要參考

相關(guān)的報(bào)告①胡燕，林黎，孫寶成，等.《天津市濱海新區(qū)地?zé)豳Y源可持續(xù)開(kāi)發(fā)潛力評(píng)價(jià)》天津：天津地?zé)峥辈殚_(kāi)發(fā)設(shè)計(jì)院，2008.。運(yùn)行Visual MODFLOW，可得到在給定水文地質(zhì)參數(shù)和源匯項(xiàng)條件下的地下水位時(shí)空分布，通過(guò)擬合同時(shí)期熱儲(chǔ)層的流場(chǎng)和典型孔的動(dòng)態(tài)過(guò)程曲線，識(shí)別熱儲(chǔ)層的參數(shù)，識(shí)別后的水文地質(zhì)參數(shù)（表1）符合實(shí)際的水文地質(zhì)條件。

圖2 長(zhǎng)觀井水位埋深觀測(cè)值與計(jì)算值曲線對(duì)比圖Fig.2 Comparison between the water-level-burial-depth curve of the long-term-observation well and the calculated curve

圖3 明化鎮(zhèn)組熱流體水位擬合圖Fig.3 Minghuazhen geothermal fluid level fitting chart

表1 模型識(shí)別后水文地質(zhì)參數(shù)取值Table 1 Values of hydrogeological parameter model identification

6 模擬模型的應(yīng)用

6.1 地?zé)峥刹少Y源量的計(jì)算

計(jì)算熱流體可采資源量采用方法的不同，可導(dǎo)致結(jié)果差異較大。熱流體可開(kāi)采資源量多少取決于開(kāi)采條件下的補(bǔ)給量，也與開(kāi)采的技術(shù)經(jīng)濟(jì)條件和開(kāi)采方案有關(guān)。在約束限制條件一定的情況下，開(kāi)采井布局的合理性成為求取可采量大小的關(guān)鍵。虛擬井采用梅花型布局，集中布置在地?zé)岬刭|(zhì)條件較好的地區(qū)，在地?zé)岬刭|(zhì)條件較差的地區(qū)適當(dāng)布置幾眼虛擬井，代入模型進(jìn)行計(jì)算。時(shí)間以2012年10月的地?zé)崃黧w流場(chǎng)作為開(kāi)始時(shí)刻，2112年10月作為結(jié)束時(shí)間；開(kāi)采量在計(jì)算中調(diào)節(jié)，求出100年地?zé)崃黧w水位最接近150 m的地?zé)崃黧w采量作為該熱儲(chǔ)層的可采資源量，結(jié)果見(jiàn)表2。

表2 明化鎮(zhèn)組地?zé)崃黧w可采量表Table 2 Minghuazhen group geothermal fluid recoverable scale

6.2 開(kāi)采預(yù)測(cè)

鑒于明化鎮(zhèn)組熱儲(chǔ)層為模擬區(qū)主要開(kāi)采層段，開(kāi)采量較大，考慮到該熱儲(chǔ)層回灌效果不好的現(xiàn)狀，本次預(yù)測(cè)暫不考慮回灌問(wèn)題，故本次在現(xiàn)有開(kāi)采井的基礎(chǔ)上，在地?zé)衢_(kāi)采較少的地段增加少量開(kāi)采虛擬井，本次布置4眼，開(kāi)采量按1000 m3/d代入模型進(jìn)行運(yùn)算，以2012年10月為初始時(shí)刻，初始水位按60 m（20℃）計(jì)算，分別預(yù)測(cè)10年和30年地?zé)崃黧w水位的

變化情況（圖4、圖5）。從預(yù)測(cè)結(jié)果可以看出，30年后模擬區(qū)明化鎮(zhèn)漏斗中心（漢沽城區(qū)）水位將達(dá)到96 m左右，年降幅約1.2 m/a。由此可見(jiàn)，在現(xiàn)有開(kāi)采量的基礎(chǔ)上，合理布設(shè)井位，每天新增1000 m3的開(kāi)采量是可行的。

圖4 明化鎮(zhèn)組10年預(yù)測(cè)地?zé)崃黧w水位等值線圖Fig.4 Minghuazhen group 10 year forecast geothermal fluid level contour map

圖5 明化鎮(zhèn)組30年預(yù)測(cè)地?zé)崃黧w水位等值線圖Fig.5 Minghuazhen group 30 year forecast geothermal fluid level contour map

7 結(jié)論

（1）本文在水文地質(zhì)條件和地?zé)岬刭|(zhì)條件分析的基礎(chǔ)上，運(yùn)用Visual Modflow三維模型，計(jì)算出研究區(qū)100年明化鎮(zhèn)組地?zé)崃黧w水位接近150 m的地?zé)崃黧w可采資源量為10.7×108m3/100 a。

（2）綜合考慮數(shù)值模擬的預(yù)報(bào)結(jié)果以及區(qū)域以往動(dòng)態(tài)資料認(rèn)為：本區(qū)在原有開(kāi)采量的基礎(chǔ)上，合理布設(shè)開(kāi)采井，每天新增1000 m3的開(kāi)采量是可行的。

[1]Anderson M P,Woessner W W.Applied groundwater modeling simulation of flow and advective transport[M]. San Diego：Academic Press,1992,215-286.

[2]郭曉東,田輝,張梅桂，等.我國(guó)地下水?dāng)?shù)值模擬軟件應(yīng)用進(jìn)展[J].地下水,2010，32(4)：5-7.

[3]李樹(shù)文，王文靜，崔新玲.基于Visual Modflow的成安縣地下水動(dòng)態(tài)研究[J].河北工程大學(xué)學(xué)報(bào)（自然科學(xué)版），2008，（4）：49-52.

[4]任印國(guó)，柳華武，李明良，等.石家莊市東部平原FEFLOW地下水?dāng)?shù)值模擬與研究[J].水文，2009，29（5）：59-62.

[5]賈金生，田冰，劉昌明.Visual MODFLOW在地下水模擬中的應(yīng)用—以河北省欒城縣為例.[J].河北農(nóng)業(yè)大學(xué)學(xué)報(bào)，2003，26(2)：71-78.

[6]高衛(wèi)東，孟磊，張海榮，等.FEFLOW軟件在地下水動(dòng)態(tài)預(yù)測(cè)中的應(yīng)用[J].上海地質(zhì)，2008，4：10-13.

[7]邵景力，崔亞莉，趙云章，等.黃河下游影響帶（河南段）三維地下水流數(shù)值模擬模型及其應(yīng)用[J].吉林大學(xué)學(xué)報(bào)（地球科學(xué)版），2003，33(1)：53-55.

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

[9]胡燕，沈鍵，劉杰，等.天津市濱海新區(qū)新近系館陶組地?zé)崃黧w流場(chǎng)三維數(shù)值模擬方法[C].第十三屆中國(guó)科協(xié)年會(huì)第14分會(huì)場(chǎng)-地?zé)崮荛_(kāi)發(fā)利用與低碳經(jīng)濟(jì)研討會(huì)論文集[M].北京:中國(guó)科學(xué)技術(shù)協(xié)會(huì)，2011：106-113.

[10]田光輝，曾梅香，程萬(wàn)慶，等.灰色系統(tǒng)GM（1，1）模型在天津地下熱水水位預(yù)測(cè)中的應(yīng)用[J].地下水，2009，31(6)：8-9.

[11]John Marler,Shemin Ge The Permeability of the Elkhom Fault Zone，South Park，Colorado[J]Groundwater, 2003，41(3):223-233.

[12]Wu Qiang,Yin Zhongmin,Wu Xun.Groundwater Study in the Weerselo Area,Overijssel,the Netherlands[M]. Beijing：Petroleum Industry Press,2001,47-74．

[13]Moench A F.Importance of the Vadose Zone in Analyses of Unconfined Aquifer Tests[J].Groundwater,2004,42 (2):223-233.

[14]張志忠，武強(qiáng).河水與地下水耦合模型的建立與應(yīng)用[J].遼寧工程技術(shù)大學(xué)學(xué)報(bào)，2004，23(4)：449-452.

[15]胡燕，林黎，林建旺，等.天津市地?zé)豳Y源可持續(xù)開(kāi)發(fā)潛力評(píng)價(jià)報(bào)告[R].天津：天津地?zé)峥辈殚_(kāi)發(fā)設(shè)計(jì)院，2007.

[16]劉杰，宋美鈺，田光輝.天津地?zé)豳Y源開(kāi)發(fā)利用現(xiàn)狀及可持續(xù)開(kāi)發(fā)利用建議[J].地質(zhì)調(diào)查與研究，2012，35(1):67-73.

[17]胡燕，高寶珠，靳寶珍，等.天津地?zé)崃黧w水化學(xué)分布特征及形成機(jī)理[J].地質(zhì)調(diào)查與研究，2007，30(3):213-218.

[18]楊玉新，高亮，張芬娜，等.天津?yàn)I海新區(qū)橋沽-看財(cái)莊地?zé)岙惓^(qū)普查報(bào)告[R].天津：天津地?zé)峥辈殚_(kāi)發(fā)設(shè)計(jì)院，2009.

[19]潘家明，眭彪，李耀義.天津?yàn)I海新區(qū)中新生態(tài)城地?zé)豳Y源賦存條件探討[J].地質(zhì)調(diào)查與研究，2010，33(3):185-189.

[20]阮傳俠，于彥，高寶珠，等.天津?yàn)I海新區(qū)地?zé)崃黧w水化學(xué)特征分析[J].地下水，2010，32(1)：51-53.

[21]田光輝，宋美鈺，劉東林，等.2012年天津地?zé)豳Y源開(kāi)發(fā)利用動(dòng)態(tài)監(jiān)測(cè)年報(bào)[R].天津：天津地?zé)峥辈殚_(kāi)發(fā)設(shè)計(jì)院，2012.

[22]姬亞?wèn)|，柴學(xué)周，劉其聲，等.大區(qū)域地下水流數(shù)值模擬研究現(xiàn)狀及存在問(wèn)題[J].煤田地質(zhì)與勘探，2009，37（5）：32-36.

Abstract:We carry out the deep-penetrating geochemical prospecting test for thick loess coverage area in the Yuncheng basin.This work adopts four methods simultaneously,including clay absorption phase metal ions(Cu Co Ni Zn etc.),thermal released carbonates(ΔC),conductivity(η%)and fluoride ion selective electrode.The test objects are two copper-nickel-polymetallic sulfide deposits buried in 200～300 meters deep under the loess.In the experiment,the four methods all captured the geochemical information from 200-meter-deep metal deposits covered by the thick loess in different degrees.It proves that the methods are effective and feasible preliminarily. There is obvious positive correlation between the anomaly intensity of the clay absorption phase metal ions(Cu Co Ni Zn)and the concentrating strength in the mineral ore.The indirect prospecting information including thermal released carbonates(ΔC)and conductivity(η%)shows the same tendency.Furthermore,we superimpose the anomaly information of the four methods and form the assemblage anomaly,which overcome the defect of single information weak and unstable,indicate the buried polymetallic ore effectively,and reflect the mineral ore anomaly more obviously.The cost of production and test of the four methods can be accepted by grass-roots unit,because the field sampling techniques are simple and easy in operation,consistent with conventional soil geochemical exploration methods.The development of the deep-penetrating geochemical prospecting method has the popularization value in similar soil coverage area.

Key words:thick loess;deep-penetrating geochemical method;clay absorption phase metal ions;thermal released carbonates(ΔC);conductivity;fluoride ion selective electrode;assemblage information

MODFLOW Application on the Numerical Simulation of Geothermal Fluid in Tianjin Binhai New Area

ZHANG Fen-na1,JIAZhi1,2,GAO Liang1,LIU Jiu-long1
（1.Tianjin geothermal exploration and designing institute,Tianjin,300250,China; 2.China University of Geosciences,Beijing,100083,China）

In order to rational development and utilization of geothermal resources in Tianjin Binhai NewArea, we take Minghuazhen geothermal fluid in the Qiaogu-Kancaizhuang Neogene Minghuazhen Formation as anexample,do numerical simulation with Visual MODFLOW software.The result Indicate that the calculation model can reflect the dynamic characteristics of groundwater level in the geothermal fluid.It is suggested that in 100 years,geothermal fluid geothermal fluid level close to 150 m recoverable reserves is about 10.7×106 m3/a. Accoding to the prediction of changes in the design of mining Minghuazhen geothermal fluid of underground water level,to add 1000 m3exploitation of dayis feasible.

Tianjin Binhai NewArea；geothermal fluid；numerical simulation；MODFLOW;calculation;Minghuazhen Formation

Study on Prospecting Test of Deep-Penetrating Geochemical Multi-Mathod in Thick Loess Coverage Area of the Zhongtiao Mountain

JIN Zhi-bin,ZHOU Xin-peng,ZHANG Shuang-kui,ZHANG Lu
(Geophysical and Geochemical Exploration Institute of Shanxi Province,Yuncheng,Shanxi,044000,China)

P314.1

A

1672－4135（2014）03－0237-05

2014-05-06

天津市國(guó)土資源與房屋管理局項(xiàng)目：天津?yàn)I海新區(qū)橋沽—看財(cái)莊地?zé)岙惓^(qū)普查：（國(guó)土房任[2008]10號(hào)）

張芬娜（1981-），女，工程師，2005年畢業(yè)于中國(guó)地質(zhì)大學(xué)（北京）水文與水資源工程，現(xiàn)主要從事地?zé)岱矫娴难芯?，E-mail:ying-ying22@163.com。