趙 振，陳惠娟，秦光雄，祁澤學(xué)

(青海省環(huán)境地質(zhì)勘查局青海省環(huán)境地質(zhì)重點實驗室，青海 西寧 810001)

西北干旱區(qū)是我國缺水最為嚴重的地區(qū)之一[1]。年平均降水量很少，飲用水的主要來源是地下水，因此，地下水資源對當?shù)貒窠?jīng)濟發(fā)展具有重要影響[2]。地下水資源數(shù)學(xué)模型是對地下水環(huán)境和地下水資源變化進行定量研究的常用方法。目前，國內(nèi)外專家學(xué)者已開發(fā)出多種地下水資源數(shù)學(xué)模型，Visual Modflow便是其中應(yīng)用最為廣泛的一種，該模型由Modflow、Modpath和MT 3D等模塊構(gòu)成，采用有限差分法進行計算，具有數(shù)據(jù)輸入及操作過程便捷、計算精度高等特點[3]。

青海省烏蘭盆地是典型的西北干旱地區(qū)，隨著烏蘭盆地社會經(jīng)濟的快速發(fā)展，對水資源的需求日益增大，已達到以水定產(chǎn)的階段。該地區(qū)的水資源包括少量地表水資源和大量的地下水資源，供水主要依賴地下水資源[4]。烏蘭盆地水文地質(zhì)條件復(fù)雜，地下水與地表水之間交換頻繁，對該地區(qū)的地下水資源進行綜合評價，進一步分析其可開采潛力，具有重要的理論意義和實用價值[5-8]。

1 研究區(qū)概況

烏蘭盆地隸屬于青海省海西蒙古族藏族自治州，位于青海省中部，海西蒙古族藏族自治州和柴達木盆地東部。地理坐標為東經(jīng)97°01′～99°27′、北緯36°19′～37°20′，是柴達木盆地的次級盆地，三面環(huán)山，東西走向，東寬西窄，盆地內(nèi)平坦開闊。氣候類型屬典型的內(nèi)陸干旱氣候，年平均氣溫4.0 ℃，年平均降水量215.9 mm，年平均蒸發(fā)量1 898.3 mm。都蘭河和賽什克河為區(qū)內(nèi)主要河流，都蘭河年均流量1.435 m3/s，賽什克河年均流量1.250 m3/s[8]。

烏蘭盆地自成一個獨立完整的水文地質(zhì)單元。區(qū)內(nèi)的第四系砂礫卵石及含泥砂礫卵石層沉積厚度巨大，自盆地中心向周邊山區(qū)逐漸變薄，出山口僅30～50 m，盆地中心松散層厚度大于300 m，且地層巖性顆粒大小亦呈由粗變細，泥質(zhì)含量逐漸增大，含水層單一的潛水含水層逐漸演變成多層結(jié)構(gòu)的承壓-自流含水層。

烏蘭盆地地下水主要排泄項有地下水開采、潛水蒸散發(fā)和排湖量等，地下水的補給來源主要為大氣降水入滲補給、農(nóng)田灌溉回歸補給、地表水滲漏補給、山前暴雨洪流入滲補給和出山口河流潛流量等[9]。流動方向呈現(xiàn)出由盆地周邊向中心匯聚，且盆地周邊徑流速率快，滲透性好，水力坡度大，滲透系數(shù)約30～100 m/d，水位埋深一般為30～50 m。中心地下水徑流速率較周邊緩慢，滲透系數(shù)約 2.2～10 m/d，水位埋深一般為1～5 m，局部地區(qū)具有承壓-自流特征，主要通過地面及湖水蒸發(fā)和少量人工開采方式排泄。烏蘭盆地地下水富水性分布極不均勻，沖洪積扇中部強于盆地中心，沖洪積扇中部單井計算涌水量可達5 359～8 885 m3/d，而盆地中心單井計算涌水量僅為629.2～728.7 m3/d。地下水水化學(xué)類型亦復(fù)雜多樣，主要為Cl·SO4-Na·Mg·Ca。

2 烏蘭盆地地下水資源現(xiàn)狀

2.1 地下水補給量

a.大氣降水入滲補給量。本文采用日降水數(shù)據(jù)計算大氣降水入滲補給量，在模型中大氣降水入滲補給量可用入滲系數(shù)法近似估算。采用Surfer V8軟件用克里格平面插值方法得到研究區(qū)多年平均月降水量等值線圖。根據(jù)烏蘭盆地地貌圖，結(jié)合不同地貌單元大氣降水入滲系數(shù)的取值，確定出模擬區(qū)大氣降水入滲補給系數(shù)平面分布(圖1)。然后將區(qū)內(nèi)近年每月降水量等值線圖與大氣降水入滲補給系數(shù)α平面分布圖進行套合，計算出每個小區(qū)域的大氣降水補給強度，α=0.10和α=0.13分布面積分別為4.28×108m2和3.17×108m2，有效降水量僅在α=0.13的區(qū)域分布，為46.04 mm/a，最終計算出的大氣降水入滲補給量為0.52萬m3/d。

b.農(nóng)田灌溉回歸補給量。研究區(qū)農(nóng)田灌溉回歸補給量包括地表水和地下水灌溉回歸補給量，其計算公式為

Q回=β灌Q灌

(1)

圖1 大氣降水入滲補給系數(shù)分區(qū)

式中：Q回為農(nóng)田灌溉回歸補給量；β灌為灌溉回滲補給地下水系數(shù)；Q灌為地表水或地下水灌溉水量。經(jīng)計算，研究區(qū)農(nóng)田灌溉回歸補給量為5.54萬m3/d。

c.河流與洪流滲漏補給量。研究區(qū)內(nèi)主要河流有都蘭河和賽什克河，兩條河流均補給地下水，根據(jù)河流徑流量與地下水的補給關(guān)系，得出河流對地下水的滲漏補給量約為2.94萬m3/d。洪流補給量按照出現(xiàn)洪流的月份分配到年內(nèi)。

d.側(cè)向徑流補給量。根據(jù)烏蘭盆地水文地質(zhì)概念模型，研究區(qū)南部、北部及東部邊界含水層存在地下水側(cè)向徑流補給，屬二類邊界。這些邊界流量可通過達西定律[10-11]求得，計算公式為

(2)

式中：Q側(cè)補為邊界側(cè)向徑流補給量；Kij為各邊界滲透系數(shù)；Iij為各邊界水力坡度；Aij為各邊界過水斷面面積。

邊界處水力坡度及過水斷面寬度在研究區(qū)地下水等水位線圖上量取(圖2)，含水層厚度取斷面處平均厚度。經(jīng)計算研究區(qū)邊界側(cè)向徑流補給總量為7.54萬m3/d。得到邊界側(cè)向徑流量的計算結(jié)果后，可進一步根據(jù)補給邊界的位置、邊界的長度將這些補給量以井的形式加入模型，并約定補給為正。

圖2 地下水側(cè)向徑流補給量計算斷面

2.2 地下水排泄量

a.地下水開采。地下水開采主要發(fā)生在盆地北部，開采主要用于農(nóng)田灌溉、人畜飲水和工業(yè)用水，根據(jù)實際調(diào)查資料，開采量較小，為1.5萬m3/d，以抽水井的形式加入模型中。

b.潛水蒸發(fā)。研究區(qū)內(nèi)地下水蒸發(fā)量根據(jù)阿維揚諾夫斯基公式采用離散化方法進行計算，計算公式為

(3)

式中：Q蒸為潛水蒸發(fā)蒸騰排泄量；E為各離散單元大水面蒸發(fā)度；Si為第i單元潛水水位埋深；ΔS為潛水蒸發(fā)極限深度，根據(jù)烏蘭縣地下水勘察成果，本文ΔS取4 m；Ai為第i單元面積；n為指數(shù)，取值為1。

根據(jù)收集的烏蘭縣氣象站2006—2013年的氣象資料，本次計算所需的日大面積水域蒸發(fā)度采用彭曼公式求得，即：

(4)

式中：ET0為應(yīng)用彭曼公式計算的日大水面蒸發(fā)度；Δ為飽和水汽壓與溫度曲線的斜率；Rn為凈輻射；G為土壤熱通量；T為2 m高度處的日平均氣溫；u2為2 m高度處的風速；es為飽和水汽壓；ea為實際水汽壓；γ為干濕表常數(shù)。由此計算的大水面多年平均蒸發(fā)度為1 206.1 mm。

根據(jù)模擬區(qū)的地形等高線及地下水流場，對水位埋深小于4 m的區(qū)域進行分區(qū)。地下水埋深小于1 m的沼澤濕地，地下水蒸發(fā)量為1 206.1 mm/a，面積為18.58 km2；1～4 m的淺埋帶，地下水蒸發(fā)量為46.19 mm/a，面積為367.79 km2。研究區(qū)總地下水蒸發(fā)量約為9.37萬m3/d。

c.湖水蒸發(fā)。研究區(qū)內(nèi)都蘭湖、柯柯鹽湖、鹽湖主要接受沖洪積扇地下水的補給，其次為大氣降水和盆地5處較大構(gòu)造上升泉的補給。2010年9月至2013年12月，湖面面積基本保持穩(wěn)定，因此，地下水年均排湖量約為年均湖面蒸發(fā)量減去年均大氣降水的直接補給量以及盆地5處較大構(gòu)造上升泉的直接入湖補給量。據(jù)此求得的地下水湖水蒸發(fā)量為6.25萬m3/d。

綜上，根據(jù)補給量總和計算得到烏蘭盆地地下水資源量為17.03萬m3/d，根據(jù)排泄量總和計算得到烏蘭盆地地下水資源量為17.11萬m3/d，根據(jù)斷面凈流量計算得到烏蘭盆地地下水資源量為17.75萬m3/d，3個計算結(jié)果接近。由于補給量總和計算考慮因素最全面，故本文烏蘭盆地地下水水資源量取值為17.03萬m3/d。

3 地下水資源數(shù)學(xué)模型

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

根據(jù)地下水資源評價原則，本次模擬以整個烏蘭盆地作為獨立完整的水文地質(zhì)單元，南以布赫特山為界，北以牦牛山為界，南北寬約15 km，東西寬約 50 km，面積約為755 km2。研究區(qū)邊界條件充分考慮了自然邊界以及地下水循環(huán)模式，反映了盆地地下水流的特征，可較好地模擬實際情況[9]。研究區(qū)屬山間斷陷盆地，含水層空間結(jié)構(gòu)具備西北干旱區(qū)特有山前沖洪積扇主要特征，具有明顯分帶。根據(jù)水文地質(zhì)鉆探揭露，沖洪積扇后緣，地下水形成單一結(jié)構(gòu)潛水，具有含水層厚度大、巖性顆粒粗、地下水位埋藏深等特點。沖洪積扇中前緣含水層過渡為沙與黏土互層的多層結(jié)構(gòu)，形成多層結(jié)構(gòu)的上層潛水及下層承壓水。研究區(qū)的地面標高采用數(shù)字高程模型表示，運用Mapgis對模擬范圍內(nèi)1∶10 000數(shù)字化地形圖進行處理，獲取研究區(qū)含水層空間結(jié)構(gòu)分布[10]。

研究區(qū)西部邊界是地表水分水嶺，垂直于地下水等水位線，可概化為隔水邊界；東北部邊界幾乎垂直于等水位線，可概化為隔水邊界；其余邊界與地下水等水位線斜交，均接受了側(cè)向的補給，可概化為第二類流量邊界；區(qū)內(nèi)湖泊排泄地下水，為了準確估算地下水開采對湖泊面積的影響，因此概化為第三類邊界。潛水面是上部邊界，大氣降水入滲、潛水蒸發(fā)等垂向水交換作用發(fā)生在該面上，因此概化為潛水面邊界[11]。下部邊界為基巖，概化為隔水邊界(圖3)。

圖3 水文地質(zhì)概念模型示意圖

3.2 地下水流數(shù)學(xué)模型

研究區(qū)地下水流數(shù)學(xué)模型采用Visual Modflow模型[12-13]，公式為

(5)

H(x,y,z,t)|t=0=H0(x,y,z) ((x,y,z)∈D)

(6)

(7)

(8)

(9)

(10)

式中：D為模擬范圍；H為地下水位；K為滲透系數(shù)；Ss為儲水率；μ為潛水重力給水度；H0(x,y,z)為地下水位初值；β為湖泊三類邊界；AB、CD為研究區(qū)西部與東北部邊界(隔水邊界)；AD、BC為研究區(qū)南、北邊界(流入邊界)；q(x,y,z,t)為二類邊界單寬流量；n為滲流區(qū)邊界法線方向；Qi為第i口井抽水量；δi為第i個單元的狄拉克函數(shù)。

據(jù)研究區(qū)含水層和地下水流場特征，將計算區(qū)域劃分成200 m×250 m的網(wǎng)格，將整個模擬區(qū)在平面上從南到北剖解成108×270個網(wǎng)格，含水層分為兩層，共有58 320個單元格，時間步長為1 d。識別模型與檢驗過程遵循以下原則：地下水流場、地下水動態(tài)過程、地下水均衡變化的模擬和實際情況基本相同；水文地質(zhì)參數(shù)方面，識別的參數(shù)符合現(xiàn)實水文地質(zhì)條件[14]。

根據(jù)2012年12月31日野外調(diào)查繪制流場為模型校正依據(jù)，流場和觀測孔分布見圖4，再次利用區(qū)域模型進行檢驗，檢驗區(qū)段為：2012年12月31日至2013年6月30日。模型識別時，采用有限差分法(Modflow軟件包)分析數(shù)學(xué)模型的數(shù)值，采用分區(qū)法反求參數(shù)[15]，即假定每一個參數(shù)區(qū)都是均質(zhì)的，把擬定的參數(shù)初值代入到有限差分法數(shù)值模型當中，計算各個時段各個節(jié)點的水位值，接著把計算的水位值與實測值進行比較，不斷調(diào)試參數(shù)值，當二者之間的誤差滿足要求時，即認定為該假定參數(shù)即代表含水層參數(shù)[16]。本研究對滲透系數(shù)、大氣降水入滲補給系數(shù)、地下水蒸發(fā)量、地表水體(河流、湖泊等)滲漏補給量及地下水開采量的處理方法主要參考文獻[17]。

模型的識別及驗證的結(jié)果主要包括以下3個方面：①地下水動態(tài)觀測水位；②模擬流場與觀測流場(地下水位空間分布采用克里金插值方法)；③模擬計算水量與長期檢測水量。擬合過程中主要對模型中的滲透系數(shù)、給水度與蒸發(fā)強度做調(diào)整[18]。圖5為代表性觀測孔水位計算值與觀測值對比，可見計算值擬合曲線與觀測值擬合曲線趨勢較一致，滿足誤差精度的要求。圖6為模型驗證期末時刻地下水水位情況，可見模擬結(jié)果正確地反映了地下水的補徑排特征。總體來看，所建地下水模型與實際水文地質(zhì)條件符合程度較高，參數(shù)分區(qū)較為合理、仿真度較高，能夠真實地反映研究區(qū)地下水流場變化規(guī)律。

(a) 初始水位

(b) 觀測孔分布

圖5 代表性監(jiān)測孔水位擬合曲線

3.3 地下水可開采資源量評價

根據(jù)上述模型以及當?shù)亟?jīng)濟發(fā)展規(guī)劃和用戶需求，都蘭湖泊面積變化、綠洲帶面積和地下水位等約束條件下，制定了相應(yīng)地下水開采方案。擬建阿木內(nèi)格溝水源地、賽什克河水源地和都蘭河水源地3處C級精度水源地，開采量分別為0.5萬m3/d、3.0萬m3/d和3.5萬m3/d。地下水總可采量為8.5萬m3/d時，研究區(qū)多年平均水均衡計算結(jié)果為：補給項包括河流滲漏3.0萬m3/d、河流潛流7.21萬m3/d，渠系、田間滲漏5.58萬m3/d、山區(qū)洪流量0.71萬m3/d、大氣降水量0.65萬m3/d；排泄項包括地下水蒸發(fā)量4.91萬m3/d、湖水排泄3.75萬m3/d。補排均衡差為-0.01萬m3/d。模型計算末時刻(20年后)地下水流場模擬見圖7，開采區(qū)降深分布見圖8。可見，烏蘭盆地地下水允許開采量為8.5萬m3/d的條件下，水源地的開采井群最大降深分別為7 m、5 m和1 m，都蘭湖和濕地沼澤退化面積百分數(shù)分別為40%和45.3%。井群地下水位降深和湖泊濕地面積減少量基本符合該區(qū)地下水可持續(xù)開發(fā)利用的生態(tài)約束條件，該開采方案是合理的。

圖6 模型驗證期末時刻地下水水位模擬結(jié)果

圖7 開采20年后地下水等水位線

圖8 開采20年后地下水水位降深

4 地下水可開采潛力分析

開展烏蘭盆地地下水可開采潛力分析，對烏蘭循環(huán)經(jīng)濟區(qū)社會經(jīng)濟發(fā)展和水資源優(yōu)化配置利用具有重要意義[19-20]。烏蘭盆地地下水可開采潛力分析，是在烏蘭盆地地下水資源評價的基礎(chǔ)上進一步評價，主要依據(jù)自然資源部發(fā)布DZ/T 0282-2015《水文地質(zhì)調(diào)查規(guī)范(1∶50 000)》進行[21]，地下水可開采潛力和開采程度計算公式為

Q潛=Q可-Q開

(11)

ν=Q開/Q可

(12)

式中：Q潛為被評價區(qū)域地下水可開采潛力；Q可為被評價區(qū)域地下水可開采資源量；Q開為被評價區(qū)域地下水實際開采量；ν為被評價區(qū)域地下水開采系數(shù)，小于0.4即為開采程度低。

目前，烏蘭盆地開采量為1.5萬m3/d，烏蘭盆地地下水可開采資源量為8.5萬m3/d，研究區(qū)地下水可開采潛力計算結(jié)果為：平水年時，Q潛為 7.0萬m3/d，ν為0.176，小于0.4，說明烏蘭盆地現(xiàn)在條件下，烏蘭盆地地下水開采程度低，開采潛力巨大。

根據(jù)評價結(jié)果，擬建水源地可開采量為 8.5萬m3/d，主要來自襲奪46.4%的潛水無效蒸發(fā)量和41.3%的排湖量。在該方案實施下，都蘭湖面積退化了13.46 km2，約為湖泊總面積的40%，沼澤濕地面積退化了8.42 km2，約為沼澤濕地總面積的45.3%，基本能夠維持湖泊及濕地周圍生態(tài)環(huán)境的良性循環(huán)。而對于沖洪積平原前緣鹽漬土的改良、抑制下游沼澤區(qū)樹木、草場及農(nóng)作物漬害有一定的正效應(yīng)。擬建水源地開采對環(huán)境的影響程度在可調(diào)節(jié)的良性范圍之內(nèi)，不會產(chǎn)生危害性的地質(zhì)災(zāi)害現(xiàn)象。

5 結(jié) 論

a.烏蘭盆地是一個獨立完整的水文地質(zhì)單元。在現(xiàn)狀條件下，盆地內(nèi)無較大取水工程，地下水多年平均補給量與排泄量基本持平，烏蘭盆地地下水資源量為17.03萬m3/d。

b.利用Visual Modflow模型模擬烏蘭盆地地下水流情況，經(jīng)驗證，該模型模擬結(jié)果與實際水文地質(zhì)條件符合程度較高，參數(shù)分區(qū)較為合理，能夠準確地模擬研究區(qū)地下水流場變化規(guī)律。

c.利用地下水數(shù)值模擬方法計算得到烏蘭盆地地下水可開采量為8.5萬m3/d，占全盆地多年平均補給量的49.9%，地下水開采程度低，潛力巨大。擬建設(shè)水源地投入開采后，該地區(qū)的濕地面積及都蘭湖水位均將有一定減少，但仍然基本能維持盆地內(nèi)湖泊及濕地周圍生態(tài)環(huán)境的良性循環(huán)。

致謝：感謝長安大學(xué)王文科教授、段磊副教授、李倩睿碩士在烏蘭盆地地下水數(shù)值模擬方面給予的指導(dǎo)。