孫驍磊，田軍倉，朱 磊(1.寧夏大學土木與水利工程學院,銀川 750021;2.寧夏節(jié)水灌溉與水資源調控工程技術研究中心，銀川 750021;3.旱區(qū)現(xiàn)代農業(yè)水資源高效利用教育部工程研究中心，銀川 750021)

井渠結合灌區(qū)具有充分利用地表水地下水、循環(huán)利用灌溉及渠系滲漏水、適時調控灌區(qū)地下水位、防治土壤發(fā)生次生鹽堿化、能夠實現(xiàn)適時灌溉、提高水資源利用率等突出特點，是我國西北地區(qū)灌區(qū)實現(xiàn)農業(yè)高效用水不可或缺的組成部分[1]。

2008年，寧夏水利廳確定了惠農渠、唐徠渠、漢渠3個典型引黃自流灌區(qū)試點，實施井渠結合灌溉綜合規(guī)劃，總灌溉面積達到0.262萬hm2，機井分布達223眼。2010年，井渠結合的灌溉面積，平均已達到1萬hm2規(guī)模，機井布置達到600眼左右[2]。但因當?shù)剞r戶反應井水涼、地下水礦化度較高使土壤積鹽等問題，人工開采進行井灌積極性不高。尤其在作物需水量較大的六七月份，黃河來水較充沛時期，機井大多處于閑置或者棄用的狀態(tài)，甚少用井水，破壞了地下水的采補平衡。只有在黃河水枯竭和斷流時期進行井水補灌。人工開采數(shù)量相對較少，使得地下水埋深超過適宜的埋深臨界。因此分析地下水-地表水水分相互轉化關系，確定適宜的地下水-地表水配水比是迫切需要解決的問題。

針對地下水-地表水的聯(lián)合耦合模擬，國內外不少學者開展了大量的工作。張浩佳，吳劍鋒等，通過降雨徑流的模擬系統(tǒng)PRMS和三維有線差分模型MODFLOW提出地下水-地表水GSFLOW模型，提出不同水文參數(shù)對于流域水資源影響提供理論和技術支撐[3]。Kim等針對SWAT模型在地下水及MODFLOW模型在地表水上模擬的不足，對SWAT模型優(yōu)先進行輸參，得到的結果由MODFLOW運行，輸出后的結果最后由SWAT進行計算。得到的模型能夠較全面的模擬地下水的補排規(guī)律以及地表水地下水的相互作用[4]。劉路廣，崔遠來等，通過將SWAT和MODFLOW中的差分網(wǎng)格進行相互對應，對改進的SWAT模型每個HRU地下水補給計算數(shù)值，再通過HRU-cells界面作為MODFLOW地下水模型補給板塊的輸入項，完成了試驗灌區(qū)地下地表水分布式模型耦合[5]。曾獻奎針對凌海市大、小凌河扇地內的水文地質條件，采用HydroGeoSphere軟件建立試驗灌區(qū)地下水-地表水水流運動及溶質運移耦合的模擬模型，選擇總氮作為模擬對象，對研究區(qū)的水流及總氮濃度進行模擬分析與預測，并對當?shù)厮Y源管理與規(guī)劃提出建議[6]。

地下水地表水耦合模擬已經趨于完善，但是目前國內利用HydroGeoSphere軟件對灌區(qū)地下水地表水進行耦合模擬的研究較少，并且通過耦合模擬對引黃灌區(qū)適宜地下水埋深進行探討的相關研究就更少。因此本文通過建立地下水-地表水耦合模擬模型，運用HydroGeoSphere 對模型進行計算求解，并通過建立的地下水模擬模型，調整用水比例，對研究區(qū)內的地下水動態(tài)變化進行預測分析，得到在生態(tài)安全埋深下，渠井用水最優(yōu)比例。

1 材料與方法

1.1 研究區(qū)概況

本文研究以寧夏銀北引黃灌區(qū)為背景，試驗灌區(qū)位于寧夏回族自治區(qū)惠農區(qū)李崗村試驗區(qū)，東經106°～107°E，北緯38°～40°N之間， 海拔1 091.0 m，地處寧夏引黃灌區(qū)，東臨黃河，西倚賀蘭山，總控制面積212.4 hm2，耕地面積143 hm2，試驗灌區(qū)的主要種植作物為小麥和玉米，以及西瓜、芹菜、番茄等經濟作物，作物生育期總灌溉水量42.373 7 萬m3，總抽水量10.307 2 萬m3，現(xiàn)狀年下渠首引水量水和抽井水量的用水比例接近4∶1。

試驗區(qū)砂礫石分布廣泛，灌區(qū)土壤的質地主要為沙壤土和壤土兩種，耕作層的土壤干密度為1.45 g/cm2，土壤全鹽量0.48～0.93 g/kg，適合農業(yè)耕作。試驗區(qū)東部邊界為官泗渠，西、南、北分別為第五排水溝、蔡家溝、艾家溝。試驗區(qū)共有抽水井19眼，農渠39條，黃河來水在每年的4月底-8月底時期。

1.2 研究內容與方法

本次研究采用田間勘測與數(shù)值計算結合的方式進行開展，研究寧夏銀北井渠結合灌區(qū)地表水與地下水水分運移規(guī)律。試驗灌區(qū)內共設置6眼地下水位的觀測井，試驗觀測內容主要為：大氣觀測、田間灌溉水量、灌區(qū)排水量觀測以及地下水的動態(tài)變化觀測。觀測井具體布置見圖1。

圖1 試驗灌區(qū)觀測井的平面布置圖Fig.1 Test area observation wells

2 地下水-地表水耦合模擬

2.1 水文地質概念模型

試驗灌區(qū)屬賀蘭山山前坳陷，地下水埋深較淺，地下水開采主要以潛水為主，含水層為單一的潛水層，潛水層以下的亞沙土、亞黏土、黏土部分作為相對隔水底板[7]。承壓水層一般的埋深在60 m以下，為遠源的補給水，本灌區(qū)不作為灌溉用水水源。因本試驗區(qū)地下水的開采未深入到承壓水含水層，所以本次模擬以潛水為主。

盡管試驗灌區(qū)地下水運動以垂向為主，水平方向上幾乎無流動，試驗區(qū)農溝較多，一部分地下水先通過各農溝匯入南北兩側斗溝(艾家溝和蔡家溝)，再由南北的斗溝排泄到西部的第五排水支溝排出。所以在概化潛水含水層的側向邊界條件時，將南北邊界，即艾家溝和蔡家溝，其主要作用為側向排水，設定為第二類邊界。西部邊界即第五排水支溝，其主要作用為側向排水，設定為第二類邊界。試驗區(qū)東部為官泗渠，概化為河流邊界。研究區(qū)內的地表水系統(tǒng)的源匯項主要有大氣降雨補給、灌溉入滲補給、蒸發(fā)蒸騰、人工開采等。

2.2 數(shù)學模型

2.2.1地下水水流數(shù)學模擬模型

根據(jù)研究區(qū)水文地質模型，采用三維Richards方程定解條件及初始條件來描述研究區(qū)飽和、非飽和地下水的非穩(wěn)定運動，建立相應的數(shù)學模型如下[8,9]：

(1)

q=-Kkr▽(ψ+z)

(2)

式中：q為地下水的滲流量；QGS為地表水入滲量；QG為地下水流的源匯項；Sw為飽和度，無量綱；Ss為貯水系數(shù)；ψ為壓力水頭；θs為飽和含水率(近似等于孔隙度)，無量綱；Ho為邊界上水頭；Γ1、Γ2和Γ3分別為給定水頭邊界、零通量邊界和給定流量邊界；qG為邊界上的水流通量；Ω為試驗區(qū)的地下水區(qū)域；K為滲透系數(shù)張量；kr為相對滲透率，無量綱；z為位置水頭。

通過van Genuchten公式得到相對滲透率kr：

(4)

式中：Se為有效含水量，無量綱；Swr為殘余含水率，無量綱；α為空氣負壓；β為孔隙大小的分布指數(shù)，無量綱；lp為孔隙連通性參數(shù)，無量綱。

2.2.2地表水水流數(shù)學模擬模型

研究區(qū)內的地表徑流受到微地形的影響，其運動屬于緩變不穩(wěn)定的波動。所以采用忽略動量項的二維圣維南方程組(二維淺水方程)描述地表徑流運動[6,8,9]。

(5)

式中：φo為地表孔隙度，無量綱；ho為水面高程；Kox、Koy為x、y方向上的地表傳導系數(shù)；Qo為地表水的源匯項；do為地表徑流深度；Do為地表水范圍；h(x,y,t)為給定水深第一類邊界條件；qo1(x,y,t)為給定通量的第二類邊界；Kn為邊界法向的地表傳導系數(shù)張量；qo2(x,y,t)為臨界深度邊界條件；qo3(x,y,t)為零深梯度邊界條件；So為地表坡度，無量綱；h(x,y,0)為初始條件。

Kox,Koy可由下式求出：

(6)

式中：nx和ny分別是x和y方向上的曼寧糙率系數(shù)。

2.2.3地下水-地表水水流數(shù)學模型的耦合

本次研究所采用的軟件HydroGeoSphere采用雙重節(jié)點法，從物理機制上將地下水和地表水進行耦合處理，進行模型的耦合時，首先將二維平面地表水模型重疊在地下水模型的頂部，對地下水、地表水模型進行相同的空間和時間的離散。將試驗區(qū)剖分為31 080個節(jié)點和54 648個單位，并選擇2015年3月27日-7月4日(共100 d)作為耦合模擬的率定期，2015年7月5日-10月12日(共計100 d)作為耦合模擬的驗證期。本次模擬的時間步長輸出為1 d。

試驗區(qū)表層的地表水模型的節(jié)點與地下水模型的頂部節(jié)點獲得相同的空間坐標，也就是說耦合模型表層的節(jié)點兼有地表水和地下水的雙重屬性，將上下層節(jié)點同時進行水力聯(lián)系。在模型中地表水與地下水體間的水頭保持完全相同，在模擬的過程中，以方程整合的方式對地表水與地下水量以及兩者之間的交換進行水力計算[6,8,9]。

在這樣的耦合方式下，并不假設相互系統(tǒng)間水頭的連續(xù)性，水流交換的描述是通過達西流關系來實現(xiàn)：

dQGS=kroKso(h-hs)

(7)

式中：d為地表水與地下水的耦合的長度；QGS為地下水和地表水的通量交換;kro為上游節(jié)點的相對滲透率，無量綱；Kso為地下水的表層介質滲透系數(shù);h為地下水水頭；hs為地表水水頭。

2.3 模型的率定與檢驗

2.3.1模型的率定

本文將試驗區(qū)在豎直方向上劃分為23層，根據(jù)抽水井的實際分布情況，以19眼抽水井平均分配的原則將試驗區(qū)劃分為六塊平面區(qū)域建立模型(見圖2)，利用GRIDBUILDER對試驗區(qū)二維平面區(qū)域圖進行前期處理，再進行網(wǎng)格剖分，在邊界處及邊界排水處，采用柵格加密進行描述(見圖3)。

圖2 試驗灌區(qū)平面區(qū)域劃分圖Fig.2 Flat test area zoning

圖3 試驗灌區(qū)模型柵格剖分圖Fig.3 Experimental zone model meshing figure

試驗灌區(qū)的水文地質參數(shù)的初始值主要通過前人的研究資料并結合實際問題得出[10,11](見表1)。選擇2015年3月27日-7月4日(共100 d)，作為模擬的率定期，確定模型的時間步長輸出為1 d。為驗證所選取的水文地質參數(shù)的真實性和所建立的數(shù)學模型的可靠性，在研究區(qū)內對已布置的6眼潛水觀測井進行觀測，利用實際觀測得到的地下水動態(tài)變化資料來對模型進行率定。試驗灌區(qū)空間離散圖，如圖4所示。

表1 主要水文地質參數(shù)的取值Tab.1 The main values of hydrogeological parameters

2.3.2模型的驗證

為了進一步驗證上述建立的模型的可靠性，選取2015年7月5日-10月12日(共計100 d)的地下水位動態(tài)觀測資料對模型進行驗證，試驗區(qū)率定期和驗證期各觀測井地下水埋深變化實測值與計算值如圖5所示。1號井、2號井、4號井、6號井種植作物為玉米，玉米全生育期灌溉3次，所列觀測井渠井灌水時間較統(tǒng)一，即5月中后旬到9月中旬為地表徑流與地下水水分相互交換的激烈期，明顯看到有3次上升和3次下降。5月初到6月中旬，灌溉的補給是導致地下水位的上升的主要原因，在這個階段，地下水量的消耗相對較少，地下水位上升約0.9 m。6月初到9月初，受到田間灌溉、大氣降雨以及潛水蒸發(fā)等因素的影響，地下水水位呈波動性變化，但總體呈現(xiàn)出緩慢下降的趨勢，下降約0.6 m。這說明，地下水的消耗不僅僅有灌區(qū)人工開采的原因，也表明作物在生育期內要消耗了一部分地下水，地下水是作物需水的不可或缺的一部分。

圖4 試驗灌區(qū)模型空間離散圖Fig.4 Discrete model space test area

圖5 率定期與驗證期內地下水埋深的實測值與模擬值對比Fig.5 Comparison of measured and simulated values of groundwater level at regular and verification stages

從圖5可以看出，經率定之后的模型在驗證期內計算出的地下水埋深與實測值的變化趨勢較為統(tǒng)一，模型計算得到的地下水埋深與地實際觀測埋深的絕對誤差小于0.2 m的時間天數(shù)占總觀測天數(shù)的90%以上。表明模型所選擇的含水層結構、概化得到的邊界條件以及水文地質參數(shù)的選取合理。同時試驗灌區(qū)的地下水系統(tǒng)特征能夠通過所建立的數(shù)學模型真實地刻畫得到，模型用來對試驗灌區(qū)的地下水埋深進行預測是真實可靠的。

圖6為模型進行計算后試驗區(qū)不同時期飽和度的分布圖，其中云圖表示飽和度的大小。圖6表示了在各個時期土壤水流飽和度的發(fā)展。從圖6可以看出，在模擬初期及3月末4月初，試驗區(qū)未進行灌溉，部分區(qū)域只以單井進行少量的灌溉，試驗區(qū)水位較低，各觀測井地下水埋深達到2 m以下，飽和度在0.69～0.76之間不斷變化。在七八月試驗區(qū)機井基本處于閑置狀態(tài)，但渠灌溉水量為最大時期，試驗區(qū)地下水位不斷升高，飽和度不斷變大。在灌水期水量達到頂峰時期，即地下水埋深小于1 m不斷靠近地表面時，飽和度垂向變化劇烈，灌水后飽和度最高達到0.9～0.96。進入九十月份，灌溉停止，土壤水流飽和度不斷降低到0.68以下，水位不斷下降至最低點達到2.2 m以下，驗證期末時段地下水的三維水位圖見圖7。

圖6 試驗灌區(qū)不同時期飽和度分布圖Fig.6 The saturation distribution of the test area in different periods

圖7 驗證期末時刻(2015-10-12)試驗區(qū)地下水三維水位圖Fig.7 Verification of the final time (2015-10-12) test area groundwater three-dimensional water level

3 基于耦合模型的最優(yōu)渠井用水比例計算分析

試驗灌區(qū)地下水與地表水相互聯(lián)系密切，灌區(qū)內渠引水量和井抽水量及渠井配水比對灌區(qū)內地下水埋深的動態(tài)變化影響大，且試驗灌區(qū)人工開采靈活，井灌設施配套齊全，通過改變渠井用水比，得到適宜的地下水埋深是切實可行的。

在不超過試驗灌區(qū)渠引水量和地下水安全開采量的條件下，渠井的用水比例不宜過大，過大會引起地下水位上升，潛水無效蒸發(fā)增加，水資源重復利用率低，引起土壤次生鹽堿化；渠井用水比例也不宜過小，過小會造成地下水嚴重超采，導致地下水漏斗[12,13]。

基于上述提出的地下水模擬，針對現(xiàn)狀年渠井比例，進行調整，對試驗區(qū)地下水埋深進行模擬預測。在現(xiàn)狀年4∶1下，不斷減小渠井用水比例，試驗區(qū)各模擬情景及計算結果見表2。

表2 試驗區(qū)各模擬情景設置下計算結果Tab.2 Calculation results under different simulation scenarios in the test area

由上述不同情景設置下計算結果表明：①在保證總用水不變的情況下，渠井用水比例不斷減小對緩解地下水水位上升有明顯的效果。通過情景1和情景5的對照分析，渠井用水比為3∶1下的渠首有效引水量比渠井用水比為2∶1下的渠首有效引水量減小了11.4%，地下水位下降近22.5%，地下水位平均下降約0.38 m。情景5與現(xiàn)狀年相比，試驗區(qū)減少渠首引黃量7.678 280 萬m3，每公頃每年減少引黃量約361.35 m3/hm2，地下水水位下降0.4 m。②隨著渠井用水比例的不斷減小，由3∶1減小到2∶1的過程中可以看出，觀測井地下水埋深的最小值也在不斷增大(見圖8)，由渠井用水比3∶1下的1.134 m減小到渠井用水比2∶1下的1.808 m，地下水埋深下降達到60%，下降0.674 m。得到在增加地下水的開采比例后，對地下水位的下降有一定的促進作用。③通過上述提出的符合寧夏引黃灌區(qū)生態(tài)安全的適宜地下水埋深，并結合上述情景計算結果，可以得到，在渠井比為2：1下的地下水埋深時，模型計算結果表明地下水最小埋深為1.808，平均埋深為2.074，符合上文提出的地下水埋深維持在1.8～3.0 m的范圍。該情景下，即2∶1的渠井用水比可以適度緩解地下水補排失衡，模擬期間，地下水位變幅變化較穩(wěn)定，人工調蓄有效措施控制地下水位持續(xù)上升，能夠避免危害性的地質災害現(xiàn)象發(fā)生。

圖8為不同的渠井用水比例下，利用模型進行預測分別計算不同渠井用水比例下地下水埋設變化情況。

圖8 不同渠井用水比例下地下水埋深變化Fig.8 Variation of groundwater depth under different water use ratio

由圖8可得適當減小渠井用水比例，以井補渠，可以得出黃河引水量減少明顯，同時適當加大開采地下水量，能夠保證地下水仍處在適宜埋深。

4 結 語

(1)在分析概化寧夏惠農李崗村試驗區(qū)水文地質的基礎上，采用雙重節(jié)點方式對地表水模型與地下水模型進行耦合，并通過HydroGeoSphere軟件對耦合模型進行計算計算求解。全面分析了試驗區(qū)灌水期地下水系統(tǒng)輸入輸出的轉化關系，在三四月以井灌為主時期，地下水開采量較大，并且處于集中開采期，地下水水位的下降幅度變化較大，地下水平均埋深在2 m左右，飽和度在0.69～0.76之間變化。七八月份為渠灌水量最大時期，地下水水位明顯升高，地下水埋深最高達到1 m，超過生態(tài)安全水位，部分水位較高地區(qū)飽和度上升到0.94～0.96。當九十月份渠首停止引水后，地下水位又重新下降，最低達到2.2 m左右，飽和度下降到0.68以下。

試驗區(qū)灌水期地下水地表水水分變化規(guī)律為當?shù)氐叵滤Y源的合理開發(fā)和科學管理提供了依據(jù)。

(2)結合地下水地表水耦合模擬模型，得出在不同情景設置下，即不同渠井用水比例下地下水埋深的變化情況。計算結果表明，引黃灌區(qū)目前渠井灌水比例在4∶1左右，地下水平均埋深接近1.5 m，越過生態(tài)安全水位。建議試驗灌區(qū)渠引水量與井抽水量比例控制在2∶1，試驗區(qū)減少渠首引黃量7.678 280 萬m3，每公頃每年減少引黃量約361.35 m3/hm2，節(jié)水效果明顯，且地下水埋深能夠維持在1.8～2.6 m的適宜安全埋深范圍內。

[1] 史海濱，田軍倉，劉慶華.灌溉排水工程學[M]. 北京：中國水利水電出版社，2006：210-216.

[2] 馬德仁，劉學軍. 寧夏井渠結合灌區(qū)管理技術研究與實踐[J].人民長江，2014，45(14):107-111.

[3] 張浩佳，吳劍鋒，林 錦，等. 基于GSFLOW 的地下水-地表水耦合模擬與分析[J].工程勘察，2015，43(5):34-38.

[4] Kim N W，Chung I M，Won Y S, et al. Development and application of the integrated SWAT-MODFLOW model [J]. Journal of Hydrology, 2008,356:1-16.

[5] 劉路廣，崔遠來.灌區(qū)地表水-地下水耦合模擬的構建[J].水利學報，2012，43(7):826-833.

[6] 曾獻奎. 基于HydroGeoSphere的凌海市大、小凌河扇地地下水-地表水耦合數(shù)值模擬研究[D]. 長春：吉林大學，2009.

[7] 席文娟.寧夏石嘴山地區(qū)地下水水質模擬與變化趨勢分析[D].西安：長安大學,2013.

[8] 朱 磊，田軍倉，孫驍磊.基于全耦合的地表徑流與土壤水分運動數(shù)值模擬[J].水科學進展，2015，26(3):322-330.

[9] R Therrien，R G McLaren, E A Sudicky, et al. HydroGeoSphere: a three-dimensional numerical model describing fully-integrated subsurface and surface flow and solute transport[R]. Waterloo: University of Waterloo, 2009.

[10] 劉玉珍. 平原區(qū)地下水系統(tǒng)模擬與決策支持研究[D]. 遼寧大連：大連理工大學, 2008.

[11] Marcel G S, Feike J L. Improved prediction of unsaturated hydraulic conductivity with themualem-van Genuchten model [J]. Soil Science Society of American Journal, 2000,64:843-851.

[12] 周維博，李佩成.井渠結合灌區(qū)節(jié)水灌溉的有效途徑[J].沈陽農業(yè)大學學報, 2004,35(5-6):473-475.

[13] 樊自立，陳亞寧，李和平，等．中國西北干旱區(qū)生態(tài)地下水埋深適宜深度的確定[J].干旱區(qū)資源與環(huán)境，2008，22(2):1-5.