孫 潔

（1.中煤科工集團(tuán)西安研究院有限公司，陜西 西安710054；2.陜西省煤礦水害防治技術(shù)重點(diǎn)實(shí)驗室，陜西 西安710077）

降水入滲補(bǔ)給既是水循環(huán)的關(guān)鍵環(huán)節(jié)，又是水均衡的關(guān)鍵要素，準(zhǔn)確估算降水入滲補(bǔ)給量關(guān)系著地下水資源的開發(fā)和利用[1]。影響降水入滲補(bǔ)給的因素很多，主要包括氣候變化（降水和蒸發(fā)）[2]、灌溉[3]、土地利用類型[4]、植被狀況[5]及地下水埋深（包氣帶厚度）[6-8]等，不同因素組合下降水入滲系數(shù)的估算較為困難。我國東部草原區(qū)分布有大量的露天礦，大型露天礦開采導(dǎo)致地下水位下降，甚至第四系地下水出現(xiàn)完全干涸的情況；另一方面，露天礦外排土場是在原始地表上堆積形成的，掩蓋了第四系水文觀測孔[7]。以上2種情況導(dǎo)致無法正常的捕捉第四系地下水動態(tài)，使得傳統(tǒng)的地下水補(bǔ)給量計算方法如地下水位波動法無法適用于高強(qiáng)度開采露天礦區(qū)。

目前國內(nèi)外對降水入滲補(bǔ)給地下水的研究較多，Huo S[8]等人根據(jù)地中滲透儀和數(shù)值模擬的方法研究了華北平原地下水位下降背景下的降水入滲系數(shù)變化規(guī)律；束龍倉[9]等利用遙感技術(shù)對城市化影響下不同土地利用類型的大氣降水入滲補(bǔ)給量進(jìn)行分析；尹德超[10]等提出以流量衰減分析為基礎(chǔ)的巖溶區(qū)次降水入滲補(bǔ)給系數(shù)計算方法。但是，以上研究大多關(guān)注農(nóng)灌區(qū)及特殊地貌區(qū)，對地下水位動態(tài)捕捉難度大的露天礦區(qū)降水入滲研究較少[11]。針對東部草原區(qū)大型露天煤礦，采煤地下水位持續(xù)下降和土壤大面積重構(gòu)形成了厚度更大、非均質(zhì)性更強(qiáng)的包氣帶，極大的影響著礦區(qū)降水-土壤水-地下水的相互轉(zhuǎn)化關(guān)系，當(dāng)入滲量較低時，地下水資源量相對較小，土壤儲水量較大，繼而有利于依靠土壤水生長的草原植被。相反，不利于草原植被的生長。因此，弄清降水入滲對露天礦區(qū)地下水位下降和土壤重構(gòu)的響應(yīng)機(jī)制成為了草原露天礦區(qū)亟待解決的科學(xué)問題。

為此以位于呼倫貝爾草原的寶日希勒露天煤礦為研究區(qū)，利用野外采集的巖土樣品在實(shí)驗室進(jìn)行顆粒分析及干密度測試，構(gòu)建研究區(qū)典型包氣帶巖性結(jié)構(gòu)剖面。根據(jù)當(dāng)?shù)貧v年降水、蒸發(fā)資料，利用數(shù)值模擬技術(shù)[12]探討露天礦區(qū)采煤地下水位下降和土壤重構(gòu)對降水入滲的影響。通過研究，為露天礦區(qū)地下水埋深較大條件下評價降水對地下水的補(bǔ)給提供簡易有效的方法，并指導(dǎo)區(qū)域地下水資源開發(fā)和排土場土壤重構(gòu)的實(shí)踐。

1 研究區(qū)概況

寶日希勒露天煤礦位于呼倫貝爾市陳旗境內(nèi)，屬于陳旗煤田的重要組成部分，西北邊界為莫爾格勒河，南邊界為海拉爾河，北部及東北部與低山丘陵相接，其宏觀地貌顯示為略有起伏的高平原，多年平均降水量為327.76 mm，多年平均潛在蒸散發(fā)量為737.67 mm。區(qū)內(nèi)地下水主要接受大氣降水的補(bǔ)給，東北部邊緣地帶接受山區(qū)基巖裂隙水的側(cè)向徑流補(bǔ)給；地下水流向基本與地形保持一致，總體呈現(xiàn)由東北向西南徑流的趨勢，最終排泄于海拉爾河、莫勒格爾河及其下游湖裙帶，受露天煤礦開采的影響，地下水局部向露天坑徑流。區(qū)內(nèi)地下水的主要排泄方式為蒸發(fā)及露天礦疏干排水。

根據(jù)研究區(qū)歷次地質(zhì)勘察資料及野外民井調(diào)查，得到研究區(qū)55口井的地下水位資料，研究區(qū)地下水流場如圖1。

圖1 研究區(qū)地下水流場Fig.1 The groundwater table in research area

從圖1可以看出，研究區(qū)露天煤礦周邊出現(xiàn)地下水降落漏斗，周邊地下水位埋深普遍較大（部分第四系水井干涸），其余地下水埋深基本小于10 m，55口井的平均地下水埋深為10 m，因此，本次研究擬設(shè)置的地下水位埋深為10 m，該情景下研究降水入滲對淺層地下水的補(bǔ)給具有實(shí)際背景意義。

2 模型構(gòu)建

2.1 典型包氣帶剖面構(gòu)建

研究區(qū)包氣帶典型剖面示意圖如圖2。針對草原區(qū)，按照2個層段采樣進(jìn)行，分別為黑色腐殖土（0～50 cm）和黃土（＞50 cm），每層重復(fù)取樣3個，共計12個典型剖面；針對排土場，排土場主要由大量的剝離物和排棄物堆積而成，上層0～50 cm為腐殖土；50 cm以下為沙礫石、亞沙土、亞黏土等組成的多層結(jié)構(gòu)，土壤組成復(fù)雜。結(jié)合寶日希勒露天煤礦排土場現(xiàn)場施工工藝，將研究區(qū)排土場的土層概化為“腐殖土（0～50 cm）+黏土（50～70 cm）+中砂（>70 cm）”的理想模型，探求排土場重構(gòu)條件下的土壤水運(yùn)移及入滲系數(shù)。

圖2 研究區(qū)典型包氣帶剖面示意圖Fig.2 The typical vadose zone profile in research area

2.2 控制方程

無植被條件下，包氣帶土壤水分運(yùn)移方程可用下式表示：

式中：C（h）為容水度，cm－1；h為土壤水壓力水頭，cm；t為時間，d；z為垂向坐標(biāo)，cm；K（h）為包氣帶的滲透系數(shù)，cm/d。

V an-Genuchten-Mualem模型可以較好地描述包氣帶土壤水分特征和滲透系數(shù)曲線[13-14]：

式中：θs為飽和含水量，cm3/cm3；θr為殘余含水量，cm3/cm3；Ks為飽和滲透系數(shù)，cm/d；m、n、α為相關(guān)土壤參數(shù)，其中m=1-1/n，n＞1；l為彎曲度參數(shù)，通常取0.5；Se為中間變量。

將草原區(qū)12個剖面的包氣帶巖土樣品進(jìn)行顆粒分析，草原區(qū)土壤顆粒分析曲線如圖3。

由圖3可以看出，草原區(qū)腐殖土和下覆黃土的顆粒分布特征差別不大，根據(jù)國際制土壤質(zhì)地分類標(biāo)準(zhǔn)，都可命名為砂質(zhì)黏壤土。由于根據(jù)土壤粒徑分布、密度等指標(biāo)采用物理-經(jīng)驗?zāi)Ｐ皖A(yù)測土水特征曲線參數(shù)基本上能夠滿足模型對精度的要求，而且已成為流域尺度研究中的一種方便可行的方法[15]，為此利用ROSSTA軟件根據(jù)砂、粉、黏粒含量生成草原區(qū)腐殖土和黃土所對應(yīng)的土壤水力參數(shù)?？紤]到排土場重構(gòu)土壤質(zhì)地的高度不均質(zhì)性，黏土和中砂的土壤水力學(xué)參數(shù)采用Carsel RF and Parrish RS的研究成果[16]，土水特征曲線參數(shù)見表1。

圖3 草原區(qū)土壤顆粒分析曲線Fig.3 The sample grain size analytical curves in the undisturbed steppe

表1 土水特征曲線參數(shù)Table 1 Hydraulic parameters of the vadose zone

2.3 初始及邊界條件

假設(shè)模擬的包氣帶剖面最初處于干燥狀態(tài)，設(shè)定0.07 cm3/cm3，即約等于殘余含水率。上邊界的氣象資料根據(jù)水文頻率計算方法，對1950—2016年的降水量進(jìn)行分析，選擇平水年1970年的降水、蒸發(fā)作為上邊界條件進(jìn)行計算；同時，對計算出來的包氣帶含水率進(jìn)行反復(fù)迭代，直至含水率出現(xiàn)穩(wěn)定，此穩(wěn)定含水率分布作為本次模擬計算的初始含水率。模型的上邊界選取大氣邊界。選取本次研究深度（10 m處）的自由排水作為模擬計算的下邊界，定義下邊界處的土壤水流通量（滲漏量）為降水對地下水的“潛在”補(bǔ)給量[17]。因此，本次土壤水運(yùn)動模擬的邊界條件可總結(jié)為：

式中：θ0為初始含水率，cm3/cm3；q0（0，t）為凈入滲速率，cm/d。

2.4 模型建立及設(shè)置

為減小初始條件對模擬結(jié)果的影響，針對1990—2016年長達(dá)27年的包氣帶土壤水運(yùn)移及降水入滲進(jìn)行模擬計算，將上文中的土壤水力學(xué)參數(shù)、1990—2016年的氣象數(shù)據(jù)及土壤初始含水率輸入模型中。其中，模擬的時間步長選取1 d，模擬的空間步長為5 cm。剖面上輸出節(jié)點(diǎn)設(shè)置遵循以下原則：0～1 m水量交換最為頻繁，每隔10 cm設(shè)置1個觀測點(diǎn)；1～3 m水量交換較為頻繁，每隔50 cm設(shè)置1個觀測點(diǎn)；3 m以下水量交換不頻繁，每隔100 cm設(shè)置1個觀測點(diǎn)。

3 結(jié)果討論

3.1 利用土壤水通量估算降水入滲系數(shù)

研究區(qū)降雨及蒸發(fā)量與模型下邊界10 m處土壤水流通量隨時間的變化特征如圖4。

圖4 降雨及蒸發(fā)量與模型下邊界10 m處土壤水流通量隨時間的變化特征Fig.4 Variation of annual precipitation,evaporation,and deep drainage over time

由圖4可知，研究區(qū)降雨及蒸發(fā)量與土壤水流通量三者的多年平均值分別為20.38、357.36、203.85 mm。其中年際降水量從232.5～619.1 mm（變異系數(shù)為24.86%），土面蒸發(fā)量從145.82～323.57 mm（變異系數(shù)為21%），與之對應(yīng)的底部土壤水通量從14.58～27.85 mm（變異系數(shù)為19%）。同時，可以看出土面蒸發(fā)量和降水量呈現(xiàn)出明顯的正相關(guān)關(guān)系，但是，10 m處土壤水通量的峰值并不和降水量的峰值保持一致，滯后于降水量高達(dá)3年左右，表明草原區(qū)包氣帶巖性（砂質(zhì)黏壤土）不利于降水的入滲。李娜等人研究了北京市大興區(qū)厚包氣帶土壤水流通量變化特征，提出了利用土壤水流通量估算降水入滲系數(shù)的方法[17]，定義為：

式中：α為降水入滲系數(shù)；Fs為多年平均土壤水通量，m/a；P為多年平均降水量，m/a。

由于野外調(diào)查研究區(qū)55口井的平均地下水埋深為10 m左右，因此考慮將模型底部10 m處的土壤水通量視為降水對地下水的潛在補(bǔ)給量，則研究區(qū)多年平均土壤水通量除以降水量即可得到草原區(qū)的降水入滲系數(shù)為0.06，其值遠(yuǎn)小于西部荒漠區(qū)[18]。

為驗證利用土壤水?dāng)?shù)值模擬方法計算降水入滲系數(shù)的準(zhǔn)確性，利用露天礦區(qū)煤礦疏放水資料確定區(qū)域計算降水入滲系數(shù)。寶日希勒露天礦及周邊煤礦已開采多年，形成了穩(wěn)定的地下水降落漏斗。根據(jù)地下水動力學(xué)原理，地下水漏斗內(nèi)排水量與該區(qū)域降水入滲量近似相等[19]，該區(qū)域地下水漏斗內(nèi)排水量與降水量的比值即為高平原區(qū)的降水入滲系數(shù)。

式中：α為降水入滲系數(shù)；F為入滲面積，m2；X為研究區(qū)多年平均降雨量，m/a，QS為穩(wěn)定疏排水量，m3。

根據(jù)內(nèi)蒙古東部高平原降落漏斗的面積、多年平均降水量、穩(wěn)定疏干排水量等參數(shù)計算出研究區(qū)的降水入滲系數(shù)為0.065，該計算結(jié)果和陳旗煤田水文地質(zhì)調(diào)查報告研究結(jié)果基本一致，也印證了利用土壤水?dāng)?shù)值模擬方法計算降水入滲系數(shù)的準(zhǔn)確性。高平原降落漏斗內(nèi)穩(wěn)定疏干排水量統(tǒng)計表見表2。

3.2 降水入滲對礦區(qū)地下水位下降的響應(yīng)

草原區(qū)降水入滲系數(shù)隨埋深的變化如圖5。隨著地下水埋深由0.5 m增加到2 m，降水入滲系數(shù)由0.11減小至0.06；當(dāng)?shù)叵滤裆畲笥? m時，降水入滲系數(shù)基本保持穩(wěn)定。因此利用深度2 m處的土壤水流通量估算研究區(qū)的降水入滲系數(shù)具有較好的可行性。露天礦區(qū)大規(guī)模煤炭開采不可避免導(dǎo)致區(qū)域地下水位下降，部分觀點(diǎn)認(rèn)為入滲系數(shù)隨著地下水埋深的增大呈減小趨勢，即采煤地下水位下降一定程度上減少了降水對地下水的補(bǔ)給。由圖5可知，地下水位下降對降水入滲補(bǔ)給地下水的影響是存在的，但是影響也是有限度的。入滲水分達(dá)到潛水極限蒸發(fā)深度以下，入滲系數(shù)與地下水埋深基本無關(guān)。地下水位的下降只是引發(fā)土壤水的入滲途徑增長，從而導(dǎo)致降水通過包氣帶運(yùn)移到潛水面的時間增大，但是對于地下水資源潛在的補(bǔ)給量影響較小。寶日希勒礦區(qū)采前地下水位埋深普遍大于2 m，因此可推斷采煤水位下降對降水入滲補(bǔ)給地下水影響較小。

表2 高平原降落漏斗內(nèi)穩(wěn)定疏干排水量統(tǒng)計表Table 2 The stable dewatering quantity of groundwater depression cone

圖5 草原區(qū)降水入滲系數(shù)隨埋深的變化Fig.5 Precipitation infiltration coefficient of the soil profile in undisturbed steppe

3.3 排土場土壤重構(gòu)對降水入滲的影響

當(dāng)?shù)叵滤宦裆畲笥? m時，降水入滲和地下水位埋深基本無關(guān)。因此，排土場水分運(yùn)移模型下邊界設(shè)置為3 m，大于極限蒸發(fā)深度。排土場穩(wěn)定含水率和入滲系數(shù)隨土壤深度的變化規(guī)律如圖6。圖6（a）表示排土場土壤剖面深度和穩(wěn)定含水率的關(guān)系，排土場50 cm以上含水率隨土壤剖面深度的增加呈增大趨勢，由于50 cm以下的土壤結(jié)構(gòu)發(fā)生改變，含水率急劇減小并保持穩(wěn)定。排土場土層結(jié)構(gòu)為“腐殖土+黏土+中砂”時，50 cm以上的土壤含水率平均為0.208 cm3/cm3，遠(yuǎn)高于草原區(qū)，說明現(xiàn)有土壤重構(gòu)模式對50 cm以上的含水率具有一定的促進(jìn)作用；圖6（b）反映了排土場降水入滲系數(shù)隨土壤埋深的變化規(guī)律，排土場結(jié)構(gòu)為“腐殖土+黏土+中砂”時，入滲系數(shù)接近0，遠(yuǎn)小于草原區(qū)的入滲系數(shù)0.06，說明寶日希勒露天礦排土場土壤重構(gòu)效應(yīng)縮減了降水對地下水的有效補(bǔ)給量。

圖6 排土場穩(wěn)定含水率和入滲系數(shù)隨土壤深度的變化規(guī)律Fig.6 The soil water content and precipitation infiltration coefficient of the soil profile in refuse dump

露天礦區(qū)排土場重構(gòu)導(dǎo)致草原局部發(fā)生地貌變化，并引發(fā)一系列的水文生態(tài)效應(yīng)，研究旱區(qū)露天煤礦復(fù)墾土壤的水分運(yùn)動特征是構(gòu)造理想新土體的前提[20]。研究結(jié)果認(rèn)為目前排土場“腐殖土+黏土+中砂”的土壤組合結(jié)構(gòu)對50 cm以上腐殖土的水分具有一定的促進(jìn)作用，但是對地下水的補(bǔ)給影響較大。在“腐殖土+黏土+中砂”的條件下，由于黏土的滲透性能較差，土壤水在黏土中向下運(yùn)移受阻，導(dǎo)致黏土上部含水率較高，即黏土的“隔水效應(yīng)”對土壤水分的運(yùn)移起主控作用。

4結(jié)論

1）根據(jù)研究區(qū)地下水埋深，結(jié)合1990—2016年的氣象資料，對降雨入滲條件下草原區(qū)包氣帶的土壤水流通量進(jìn)行了數(shù)值模擬計算。模擬結(jié)果顯示：自然降雨條件下10 m處土壤水分通量平均值為20.38 mm，降雨入滲系數(shù)為0.06，結(jié)果和利用礦區(qū)疏放水資料計算的降水入滲系數(shù)較為一致，驗證了利用土壤水?dāng)?shù)值模擬估算降水入滲系數(shù)的準(zhǔn)確性。

2）隨著地下水埋深由0.5 m增加到2 m，降水入滲系數(shù)由0.11減小至0.06，當(dāng)?shù)叵滤裆畲笥? m時，降水入滲系數(shù)基本保持穩(wěn)定，因此利用深度2 m處的土壤水流通量估算研究區(qū)的降水入滲系數(shù)具有較好的可行性；寶日希勒礦區(qū)采前地下水位埋深普遍大于2 m，因此可推斷采煤水位下降對降水入滲補(bǔ)給地下水影響較小。

3）相比草原區(qū)，目前露天礦排土場“腐殖土+黏土+中砂”的模式，50 cm以上的土壤含水率平均為0.208 cm3/cm3，遠(yuǎn)高于原生草原區(qū)，說明現(xiàn)有土壤重構(gòu)模式大大提高了50 cm以上土壤的含水率，對露天礦排土場的植被恢復(fù)有積極作用，但一定程度縮減了降水入滲系數(shù)，不利于地下水接受大氣降水的入滲補(bǔ)給。