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

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

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

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

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

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

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

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

1 研究區(qū)概況

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

2 熱儲特征

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

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

分布，儲集條件良好，厚度發(fā)育穩(wě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)組年開采量69×104m3，水位58～62 m，年降幅1 m左右[21]。

3 概念模型

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

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

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

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

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

5 模型數(shù)值模擬

5.1 時空離散

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

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

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

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

5.3 模型的識別驗證

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

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

圖2 長觀井水位埋深觀測值與計算值曲線對比圖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 模型識別后水文地質(zhì)參數(shù)取值Table 1 Values of hydrogeological parameter model identification

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

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

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

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

6.2 開采預(yù)測

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

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

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

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

7 結(jié)論

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

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

[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]郭曉東,田輝,張梅桂，等.我國地下水?dāng)?shù)值模擬軟件應(yīng)用進展[J].地下水,2010，32(4)：5-7.

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

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

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

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

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

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

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

[10]田光輝，曾梅香，程萬慶，等.灰色系統(tǒng)GM（1，1）模型在天津地下熱水水位預(yù)測中的應(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]張志忠，武強.河水與地下水耦合模型的建立與應(yīng)用[J].遼寧工程技術(shù)大學(xué)學(xué)報，2004，23(4)：449-452.

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

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

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

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

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

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

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

[22]姬亞東，柴學(xué)周，劉其聲，等.大區(qū)域地下水流數(shù)值模擬研究現(xià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

天津市國土資源與房屋管理局項目：天津濱海新區(qū)橋沽—看財莊地?zé)岙惓^(qū)普查：（國土房任[2008]10號）

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