神東礦區(qū)降雨入滲補(bǔ)給強(qiáng)度研究*

何 淵

(中煤科工集團(tuán)西安研究院有限公司，陜西 西安710054)

0 引 言

神東礦區(qū)地處陜北黃土高原北部、鄂爾多斯高原毛烏素沙漠的東南緣。該區(qū)屬干旱半干旱溫帶大陸性氣候，降雨稀少，蒸發(fā)強(qiáng)烈，生態(tài)環(huán)境脆弱，礦區(qū)的發(fā)展帶來了一系列水資源及生態(tài)環(huán)境問題，合理開發(fā)利用地下水資源需要對(duì)區(qū)域的地下水資源做出科學(xué)的評(píng)價(jià)，已有研究表明，區(qū)域地下水循環(huán)方式主要以垂向交換為主［1］。作為地下水資源評(píng)價(jià)中補(bǔ)給量的主要構(gòu)成部分，以往的計(jì)算大多采用經(jīng)驗(yàn)公式計(jì)算的方法簡(jiǎn)單的將降雨入滲補(bǔ)給地下水量概化為降雨入滲系數(shù)，誤差較大。本次研究采用野外調(diào)查包氣帶巖性結(jié)構(gòu)及取樣，采用數(shù)值法模擬降雨入滲過程，模型中參數(shù)由實(shí)驗(yàn)室測(cè)定給出，對(duì)于土樣的實(shí)驗(yàn)室測(cè)定主要采用離心機(jī)、環(huán)刀法及馬爾文激光顆粒分析儀測(cè)定各項(xiàng)土壤水力學(xué)參數(shù)，在此基礎(chǔ)上結(jié)合當(dāng)?shù)氐臍庀筚Y料建立研究區(qū)包氣帶數(shù)值模擬模型，計(jì)算出區(qū)域不同地段降雨入滲強(qiáng)度及入滲系數(shù)，同時(shí)得出的參數(shù)也可用于類似區(qū)域的水文地質(zhì)計(jì)算［2］，為地下水資源評(píng)價(jià)提供科學(xué)依據(jù)。

1 野外調(diào)查及取樣

項(xiàng)目組于2013 年在研究區(qū)由南至北垂直烏蘭木倫河按地質(zhì)地貌特征布置了3 條剖面線進(jìn)行包氣帶巖性調(diào)查，同時(shí)取樣。取樣范圍基本控制研究區(qū)域，每條剖面布置4 個(gè)取樣點(diǎn)，沿剖面布置的每個(gè)取樣點(diǎn)取包氣帶埋深0 ～3 m 顆分樣，同時(shí)在示范礦井(補(bǔ)連塔與大柳塔煤礦)范圍及包氣帶巖性(風(fēng)積沙和黃土狀圖)典型區(qū)域取0 ～3 m 包氣帶原狀樣，每個(gè)深度原狀樣均取雙樣。

2 研究區(qū)包氣帶巖性分布

根據(jù)野外調(diào)查及結(jié)合前人已做過的地質(zhì)地貌圖，遙感解譯圖等資料分析，研究區(qū)在地貌單元上可劃分為風(fēng)沙地貌和黃土地貌2 大類。以風(fēng)沙地貌(土壤質(zhì)地分類為砂土、壤質(zhì)砂土)為主，覆蓋區(qū)內(nèi)90%以上面積，其地表形態(tài)表現(xiàn)為風(fēng)沙灘地及沙丘-沙梁，在風(fēng)沙灘地區(qū)，湖泊海子星羅棋布;黃土地貌(土壤質(zhì)地分類為砂質(zhì)壤土、壤土為主)僅在區(qū)內(nèi)東南部沿溝壑地帶分布，地形切割較深，沖溝發(fā)育，縱橫交錯(cuò)，溝谷陡峻狹窄，植被稀少，水土流失嚴(yán)重［3-5］。

圖1 水文地質(zhì)概念模型示意圖Fig.1 Hydrogeological conceptual model

3 包氣帶水分運(yùn)移數(shù)值模擬

3.1 水文地質(zhì)概念模型

根據(jù)野外典型點(diǎn)取樣的包氣帶巖性結(jié)構(gòu)概化出的水文地質(zhì)概念模型如圖1 所示。由于所計(jì)算區(qū)域包氣帶中水分運(yùn)移主要是垂向運(yùn)動(dòng)，因此，將其概化為垂向一維流［6-8］。試驗(yàn)介質(zhì)概化為2 層，z 軸坐標(biāo)原點(diǎn)選在地表，向上為正。已有研究認(rèn)為，在包氣帶厚度小于潛水蒸發(fā)極限深度的條件下，隨著包氣帶厚度的增大，降雨入滲強(qiáng)度逐漸減小，大于3 m 后，年降雨入滲補(bǔ)給強(qiáng)度將趨于穩(wěn)定，與埋深無明顯相關(guān)關(guān)系。根據(jù)本次研究收集區(qū)域內(nèi)地質(zhì)、水文地質(zhì)資料分析，神東礦區(qū)內(nèi)地下水位埋深除溝谷地段外基本均大于3 m，故本次數(shù)值模擬的包氣帶計(jì)算剖面厚度概化為3 m.

3.2 參數(shù)獲取

3.2.1 土壤水分特征曲線

此次土壤水力學(xué)參數(shù)獲取的試驗(yàn)采用日本himac 公司生產(chǎn)的CR21G 離心機(jī)測(cè)定土壤水分特征曲線，測(cè)定時(shí)溫度恒定20 ℃.對(duì)于原狀土，首先將標(biāo)準(zhǔn)環(huán)刀(面積100 cm2)取得原狀土樣浸水飽和，并測(cè)定其飽和含水量。然后將飽和的原狀土樣品放入離心機(jī)裝置中，根據(jù)設(shè)計(jì)壓力設(shè)定離心機(jī)轉(zhuǎn)速，加壓范圍為10 ～1 000 kPa，分別為10，20，40，60，80，100，200，400，600，800，1 000 kPa.在每個(gè)壓力達(dá)到平衡后，取出樣品通過稱重獲得質(zhì)量含水量。試驗(yàn)測(cè)得8 個(gè)原狀樣土壤水分特征曲線如圖2 所示。

本次計(jì)算利用RETC 軟件對(duì)8 個(gè)原狀樣利用實(shí)驗(yàn)室測(cè)得的土壤水分特征曲線進(jìn)行參數(shù)反演，得出8 個(gè)原狀樣的土壤水力學(xué)參數(shù)。其結(jié)果見表1.

表1 土壤水力學(xué)參數(shù)測(cè)定結(jié)果表Tab.1 Soilhydraulic parameter(original state samples)

圖2 原狀樣土壤水分特征曲線Fig.2 Soil water characteristic curve of the original state sample

3.2.2 土壤飽和導(dǎo)水率

土壤飽和導(dǎo)水率測(cè)定采用環(huán)刀法，用環(huán)刀取原狀土樣，浸水后，在單位水壓梯度下，根據(jù)達(dá)西定律，求得通過垂直于水流方向的單位土壤截面積的水流速度，稱為土壤的飽和導(dǎo)水率［9-10］。測(cè)得結(jié)果見表2.

表2 原狀樣土壤飽和導(dǎo)水率測(cè)定結(jié)果Tab.2 Saturated Soil hydraulic conductivity of original state sample

3.2.3 土壤顆粒分析

顆粒分析采用世界通用的馬爾文激光顆粒分析儀測(cè)定(型號(hào)為MS2000)，該分析儀是采用激光衍射法來測(cè)定土壤顆粒之間的粒度大小及粒度分級(jí)。

本次計(jì)算中顆粒分析土樣的土壤水力參數(shù)獲取是利用RETC 軟件中的土壤轉(zhuǎn)換函數(shù)功能，即根據(jù)土壤的顆粒級(jí)配中砂粒、粉粒、粘粒的百分含量以及土壤容重等土壤物理性質(zhì)數(shù)據(jù)，直接輸出Van-Genuchten 模型中的參數(shù)。

顆粒分析土樣參數(shù)反演結(jié)果見表3.

3.3 數(shù)學(xué)模型

式中 θ 為土壤體積含水率，cm3·cm-3;k 為非飽和滲透系數(shù)，cm·hr-1;t 為時(shí)間變量，hr-1;z 為空間變量，cm;h 為壓力水頭(在包氣帶中稱為負(fù)壓);地表為原點(diǎn)，向上為正。

表3 顆粒分析樣土壤水力學(xué)參數(shù)測(cè)定結(jié)果表Tab.3 Soil hydraulic parameter(particle-size analysis of soil samples)

)

3.4 模型的初始條件和邊界條件

模型的初始條件設(shè)置為計(jì)算剖面的負(fù)壓值，假定地表為負(fù)壓值為0 m，計(jì)算剖面底部(3 m)處負(fù)壓值為-3 m，線性分布。

模型的上邊界條件由氣象條件決定(本次計(jì)算采用研究區(qū)內(nèi)氣象站近20 a 的氣象資料)，HYDRUS-1D 軟件采用能量平衡方程將連續(xù)的氣象資料轉(zhuǎn)換為地表的能量通量。并計(jì)算出地表蒸發(fā)作為模型的上邊界條件。模型下邊界條件考慮為自由排泄邊界［11］。

3.5 數(shù)學(xué)模型的求解

上述數(shù)學(xué)模型求解采用美國(guó)鹽土實(shí)驗(yàn)室開發(fā)的hydrus-1D 軟件解算此模型。數(shù)值離散采用有限元法，由程序自動(dòng)完成。長(zhǎng)度單位定為cm，時(shí)間單位定為d.初始時(shí)間步長(zhǎng)為0.01 d，計(jì)算過程中，程序會(huì)根據(jù)計(jì)算收斂的情況自定調(diào)整時(shí)間步長(zhǎng)，設(shè)定最小時(shí)間步長(zhǎng)為0.000 01 d，最大時(shí)間步長(zhǎng)為1 d［12-16］。

經(jīng)數(shù)值模擬計(jì)算，研究區(qū)95%降水頻率年地下水入滲補(bǔ)給量計(jì)算的結(jié)果見表5.

設(shè)地下水獲得的降雨入滲補(bǔ)給量為pr，大氣降水量記為p，降雨入滲補(bǔ)給系數(shù)為，根據(jù)降雨入滲補(bǔ)給系數(shù)的定義，有α =pr/p. 計(jì)算得出區(qū)域各計(jì)算點(diǎn)地下水不同的入滲補(bǔ)給系數(shù)，見表4.

表4 研究區(qū)降雨入滲強(qiáng)度及入滲系數(shù)匯總表Tab.4 Rainfall infiltration intensity and coefficient in study area

4 結(jié) 論

1)研究區(qū)在地貌單元上可劃分為風(fēng)沙地貌和黃土地貌兩大類。在風(fēng)沙灘地區(qū)以風(fēng)沙地貌(土壤質(zhì)地分類為砂土、壤質(zhì)砂土)為主，覆蓋區(qū)內(nèi)90%以上面積，其地表形態(tài)表現(xiàn)為風(fēng)沙灘地及沙丘，沙梁;黃土地貌(土壤質(zhì)地分類為砂質(zhì)壤土、壤土為主)僅在區(qū)內(nèi)東南部沿溝壑地帶分布。

2)本次數(shù)值模擬的包氣帶計(jì)算剖面厚度取為3 m，結(jié)合當(dāng)?shù)氐臍庀髼l件及典型點(diǎn)巖性結(jié)構(gòu)計(jì)算研究區(qū)降雨入滲強(qiáng)度及入滲補(bǔ)給系數(shù)。從研究區(qū)降雨入滲系數(shù)計(jì)算結(jié)果可以看出，區(qū)域降雨入滲系數(shù)大致范圍在18.4% ～27.7%.

3)區(qū)域降雨入滲強(qiáng)度主要受到當(dāng)?shù)貧庀髼l件、包氣帶巖性結(jié)構(gòu)及地下水位埋深(即包氣帶厚度)的影響。由本次計(jì)算結(jié)果可以看出，在相同降雨特征和水位埋深條件下，粗顆粒巖性的降雨入滲系數(shù)大于細(xì)顆粒巖性，即砂土＞壤質(zhì)砂土＞砂質(zhì)壤土＞壤土＞粉壤土。所以在研究區(qū)內(nèi)地表廣泛分布的風(fēng)積沙對(duì)降雨入滲十分有利，對(duì)地下水資源的保護(hù)起到積極的作用。

4)從計(jì)算的結(jié)果來看，在研究區(qū)除溝谷地段只要地下水埋深大于潛水蒸發(fā)的極限埋深，潛水蒸發(fā)微弱，對(duì)潛水面來講降水補(bǔ)給是主要的，研究區(qū)廣泛分布的風(fēng)積砂對(duì)區(qū)域降雨入滲起到了積極的作用。

5)文中擬合出的研究區(qū)各巖性水文地質(zhì)參數(shù)可用于該區(qū)域及類似區(qū)域的水文地質(zhì)計(jì)算。

6)包氣帶水分運(yùn)移規(guī)律受多種因素影響，本次研究?jī)H考慮了降雨量、蒸發(fā)量、包氣帶巖性、厚度等因素的影響，對(duì)于氣、熱、植被等的影響因素未來尚需開展專門的研究工作。

References

［1］ 何 淵.鄂爾多斯盆地沙漠高原區(qū)湖泊和潛水面蒸發(fā)能力研究［D］.西安:長(zhǎng)安大學(xué)，2006.HE Yuan.Research on evaporation intensity of lakes and groundwater in desert plateau of Ordos Basin［D］. Xi’an:Chang’an University，2006.

［2］ 黃金廷.鄂爾多斯盆地沙漠高原區(qū)降雨入滲補(bǔ)給地下水研究［D］.西安:長(zhǎng)安大學(xué)，2006.HUANG Jin-ting.Research on rainfall recharge in desert plateau of Ordos Basin［D］.Xi’an:Chang’an University，2006.

［3］ 肖 暉，李建新，韓 偉，等.鄂爾多斯盆地南緣渭北隆起中新生代構(gòu)造抬升及演化［J］.西安科技大學(xué)學(xué)報(bào)，2013，33(5):359 -363.XIAO Hui，LI Jian-xin，HAN Wei，et al. The tectonic uplift time and evolution characteristics of Weibei uplift in the south edge of Ordos Basin［J］.Journal of Xi’an University of Science and Technology，2013，33(5):359-363.

［4］ 崔宏俊，劉寶平，歐陽畹町.鄂爾多斯盆地延長(zhǎng)氣田石盒子組沉積特征分析［J］. 西安科技大學(xué)學(xué)報(bào)，2012，32(6):726 -731.CUI Hong-jun，LIU Bao-ping，OUYANG Wan-ding. Analysis of the sedimentary characteristics of Shihezi formation in Yanchang gasfield，Ordos Basin［J］.Journal of Xi’an University of Science and Technology，2012，32(6):726 -731.

［5］ 王 云，梁 明，汪桂生.基于ArcGIS 的流域水文特征分析［J］.西安科技大學(xué)學(xué)報(bào)，2012，32(5):581 -585.WANG Yun，LIANG Ming，WANG Gui-sheng. Analysis of hydrological features based on ArcGIS［J］.Journal of Xi’an University of Science and Technology，2012，32(5):581 -585.

［6］ 李 翔，王金安，張少杰.復(fù)雜地質(zhì)體三維數(shù)值建模方法研究［J］.西安科技大學(xué)學(xué)報(bào)，2012，32(6):676-681.LI Xiang，WANG Jin-an，ZHANG Shao-jie.3D modeling method of complicated geological body［J］. Journal of Xi’an University of Science and Technology，2012，32(6):676 -681.

［7］ 計(jì) 宏，余學(xué)義.厚黃土覆蓋層采動(dòng)破壞數(shù)值模擬研究［J］.西安科技大學(xué)學(xué)報(bào)，2012，32(2):209 -213.JI Hong，YU Xue-yi. Numerical simulation of surface subsidence in thick loess mining seam［J］. Journal of Xi’an University of Science and Technology，2012，32(2):209 -213.

［8］ 白如紅，郝君，何龍偉，等.淺埋煤層綜采工作面礦壓顯現(xiàn)規(guī)律數(shù)值研究［J］.西安科技大學(xué)學(xué)報(bào)，2012，32(2):676 -681.BAI Ru-hong，HAO Jun，HE Long-wei，et al.Numerical study of mine pressure regularity in fully-mechanized face of shallow coal seam［J］.Journal of Xi’an University of Science and Technology，2012，32(2):676 -681.

［9］ 李志強(qiáng)，唐 旭.水射流卸壓增滲及抽采瓦斯效果的滲流力學(xué)數(shù)值解［J］. 西安科技大學(xué)學(xué)報(bào)，2012，32(4):464 -469.LI Zhi-qiang，TANG Xu. Stress relief and permeability increasing by water jet and numerical method of seepage mechanics equation for gas extraction effect［J］.Journal of Xi’an University of Science and Technology，2012，32(4):464 -469.

［10］雷志棟，楊詩秀，謝森傳. 土壤水動(dòng)力學(xué)［M］. 北京:清華大學(xué)出版社，1988.LEI Zhi-dong，YANG Shi-xiu，XIE Sen-chuan.Soil water dynamics［M］.Beijing:Tsinghua University Press，1988.

［11］唐麗云，楊更社，趙 濤，等.重力壩壩體水平變形模型研究［J］.西安科技大學(xué)學(xué)報(bào)，2012，32(2):193 -197.TANG Li-yun，YANG Geng-she，ZHAO Tao，et al.Horizontal deformation model of gravity dam［J］. Journal of Xi’an University of Science and Technology，2012，32(2):193 -197.

［12］王松陽，李樹剛，王紅勝.水滲流對(duì)薄層復(fù)合頂板的弱化影響［J］.西安科技大學(xué)學(xué)報(bào)，2012，32(6):686-690.WANG Song-yang，LI Shu-gang，WANG Hong-sheng.Weaken influence of water seepage on thin-layer compound roof［J］. Journal of Xi’an University of Science and Technology，2012，32(6):686 -690.

［13］ 周 明.利用地中滲透儀觀測(cè)資料進(jìn)行降雨入滲補(bǔ)給規(guī)律分析［J］. 地質(zhì)科技情報(bào)，2002，21(1):37 -40.ZHOU Ming.Analysis of precipitation recharge using observed data of lysimeter［J］. Geological Science and Technology Information，2002，21(1):37 -40.

［14］趙貴章.鄂爾多斯盆地風(fēng)沙灘地區(qū)包氣帶水——地下水轉(zhuǎn)化機(jī)理研究［D］.西安:長(zhǎng)安大學(xué)，2011.ZHAO Gui-zhang.Study on transformation mechanism of vadose zone water-groundwater in the wind-blown sand area of the Ordos Basin［D］.Xi’an:Chang’an University，2011.

［15］馮 平，李建柱.土壤水分下滲機(jī)制及其在半干旱區(qū)產(chǎn)流模擬中的應(yīng)用［J］.干旱區(qū)資源與環(huán)境，2008，22(6):95 -98.FENG Ping，LI Jian-zhu. Ifiltration mechanism of soil moisture and its application in run off generation modelling［J］. Journal of Arid Land Resources and Environment，2008，22(6):95 -98.

［16］杜濤濤，竇林名，藍(lán) 航.定向水力致裂防沖原理數(shù)值模擬研究［J］. 西安科技大學(xué)學(xué)報(bào)，2012，32(4):444 -449.DU Tao-tao，DOU Lin-ming，LAN Hang. Simulation on rockburst prevention by directional hydraulic fracture［J］.Journal of Xi’an University of Science and Technology，2012，32(4):444 -449.