基于HydroGeoSphere的河谷地區(qū)地表水地下水水流水質(zhì)聯(lián)合模擬

張將偉，盧文喜，陳 末，林 琳，李久輝，崔尚進

(吉林大學(xué)環(huán)境與資源學(xué)院，吉林 長春 130012)

從水循環(huán)的角度來看，地表水與地下水均是水文循環(huán)的重要一環(huán)，二者在水量與水質(zhì)上相互作用，緊密聯(lián)系。為了更加合理的反映流域水文循環(huán)實際過程，應(yīng)當(dāng)將地下水與地表水系統(tǒng)作為一個整體進行模擬分析[1]。自20世紀(jì)70年代Pinder等對耦合模型進行了研究探索[2]后，地下水地表水耦合模型迅速發(fā)展，Govindaraju等建立地表渠系水流和地下水流耦合模型；Maxwell等將地下水地表水耦合模型ParFlow與氣象模型ARPSZ結(jié)合[3]；王中根等耦合地表水模型SWAT與地下水模型MODFLOW，應(yīng)用于海河流域[4]。但上述模型僅僅考慮水流耦合模擬，而忽視了地表水地下水在水質(zhì)方面的聯(lián)系，現(xiàn)代社會水量與水質(zhì)均是水資源開發(fā)利用的核心制約[5]，進行地表水地下水水流水質(zhì)聯(lián)合模擬是發(fā)掘區(qū)域水資源系統(tǒng)潛力的必要手段。

目前已有多種模型可以進行地表水地下水聯(lián)合模擬，其中HydroGeoSphere模型從水文循環(huán)物理過程出發(fā)，通過耦合方程將地表水地下水?dāng)?shù)學(xué)模擬模型平列起來，同時求解，可以更客觀地描述水分轉(zhuǎn)化關(guān)系與溶質(zhì)運移過程[6]。本文綜合考慮了坡面漫流、河道流與包氣帶水、飽和帶水的相互聯(lián)系，利用HydroGeoSphere軟件建立地表水地下水水量水質(zhì)聯(lián)合模擬模型，藉此得到降雨量變化情況下研究區(qū)內(nèi)水流和總氮的分布特征與變化規(guī)律，為今后進行地表水地下水聯(lián)合調(diào)度、污染預(yù)報等方面研究提供科技支撐，能夠更好地發(fā)揮區(qū)域水資源的整體性功能。

1 假想例簡介

建立了一個南北長500 m，東西寬500 m的河谷假想例作為研究區(qū)，如圖1與圖1.2所示。研究區(qū)中部存在一河流，流向為北至南流動，經(jīng)入流斷面流入研究區(qū)(如圖2中Γ4所示)，由出流斷面流出(如圖2中Γ2所示)。研究區(qū)東西兩側(cè)山脊線為分水嶺(如圖2中Γ1，Γ3所示)。研究區(qū)多年平均降雨量為800 mm/a，水面蒸發(fā)量為1 000 mm/a。研究區(qū)內(nèi)土地利用類型有草地和灘地，如圖2所示，其地表水水力參數(shù)取經(jīng)驗值如表1。地表水系統(tǒng)源匯項主要有河口流入量、大氣降水補給、地表水體蒸發(fā)、地表水下滲補量以及河口流出量。

圖1 研究區(qū)三維圖Fig.1 3-D diagram of research area

圖2 研究區(qū)平面圖Fig.2 Planimetric map of research area

分區(qū)曼寧粗糙系數(shù)洼地儲量/m消減儲量/m耦合長度/m草地0.300.250.0020.01灘地0.030.250.0020.01

研究區(qū)含水層以第四系沉積物為主，水文地質(zhì)參數(shù)取值如表2所示。含水層南側(cè)和北側(cè)(如圖2中Γ2,Γ4所示)設(shè)為定水頭邊界，水頭值分別為7 m和9 m，東側(cè)和西側(cè)邊界(如圖2中Γ1,Γ3所示)為隔水邊界。研究區(qū)地下水系統(tǒng)源匯項有北部邊界流入量、大氣降水、地表水下滲量、抽水井抽水量及南部邊界流出量。

表2 研究區(qū)水文地質(zhì)參數(shù)表Tab.2 Value of hydrogeology parameter of study area

研究區(qū)南部存在露天垃圾堆放場，污染物以氮元素污染為主。設(shè)置一東西向監(jiān)測線橫穿垃圾場，監(jiān)測污染物下滲情況(如圖2中A-A'所示)。場地內(nèi)縱向彌散系數(shù)10 m，橫向彌散系數(shù)1 m。研究區(qū)東側(cè)、北側(cè)和西側(cè)(如圖2中Γ1,Γ3,Γ4所示)為零溶質(zhì)通量邊界；南側(cè)水流交換通量邊界(如圖2中Γ2所示)概化為第三類邊界條件。研究區(qū)南部擬建一抽水井，抽水方案如表3所示。

表3 抽水方案表Tab.3 Scheme of pumping

2 數(shù)學(xué)模擬模型

2.1 水流數(shù)學(xué)模型

2.1.1 地表水水流數(shù)學(xué)模擬模型

地表水的運動采用二維平均深度圣維南(Saint Venant)方程來描述，如模型(1)：

(1)

式中：do為地表徑流深度，[L]；Kox和Koy分別為X軸和Y軸方向上的等效水力傳導(dǎo)系數(shù)，[LT-1]；h0(x,y)為(x,y)點初始水位，[L]；Q為地表水的源匯項，為單位面積體積流量，[LT-1]；Γgs為地表水與地下水之間交換通量，[T-1]；φ0為等效地表孔隙度，無量綱；Γ1,Γ3為分水嶺，Γ2,Γ4入流斷面和出流斷面。

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

研究區(qū)內(nèi)地下水為三維變飽和地下水流，使用Richard方程描述其水流運動，如模型(2)：

(2)

式中：Kxx,Kyy,Kzz為滲透系數(shù)，[LT-1]；Krw為相對滲透率，無量綱；H為地下水水頭，[L]；W為地下水的源匯項，為單位體積流量，[T-1]；Γgs為地表水與地下水之間交換通量，[T-1]；φ為孔隙度，無量綱；Sw為飽和度，無量綱；Ss為儲水系數(shù)，[L-1]；qG為邊界上的水流通量，[LT-1]；H0為地下水初始水頭，[L]；Γ1,Γ3為隔水邊界， 已知流量邊界；DG為研究區(qū)地下水區(qū)域。

2.1.3 地表水-地下水水流模型耦合

本次研究采用雙重節(jié)點法對地下水地表水模型進行耦合。這種方法假設(shè)地表水與地下水之間存在一個水量交換層，地表水與地下水交換過程符合達西定律。并在達西定律的基礎(chǔ)上引入相對滲透系數(shù)，反映包氣帶對地表水地下水水量交換的影響，如(3)式所示：

(3)

式中：qgs為地下水和地表水之間單位面積體積流量[LT-1]；do為地表徑流深度[L]；Γgs為地表水與地下水之間交換通量[T-1]；krw為相對滲透率，無量綱；Kzz為地下水表層介質(zhì)滲透系數(shù)[LT-1]；H表示地下水水頭[L]；h表示地表水水頭[L]；lexch為地表水與地下水之間的耦合長度[L]。

2.2 溶質(zhì)運移數(shù)學(xué)模型

2.2.1 地表水溶質(zhì)運移數(shù)學(xué)模擬模型

地表水中溶質(zhì)二維運移過程數(shù)學(xué)模擬方程如模型(4)：

(4)

式中：qo為地表水流速，[LT-1]；c為地表水溶質(zhì)濃度，ML-3；Do為地表水水動力彌散系數(shù)，L2T-1；Ωs為地表水溶質(zhì)的源匯項，ML-3T-1；Ωgs為地下水地表水之間溶質(zhì)交換量，ML-3T-1；Ro地表水溶質(zhì)阻滯系數(shù)，無量綱；Γ1,Γ2,Γ3,Γ4為紐曼邊界。

2.2.2 地下水溶質(zhì)運移數(shù)學(xué)模擬模型

使用三維變飽和溶質(zhì)運移模型刻畫地下水溶質(zhì)變化，控制方程如模型(5)：

(5)

式中：q為地下水流速，量綱[LT-1]；C為地下水溶質(zhì)濃度，[ML-3]；D為水動力彌散系數(shù)，[L2T-1]；Ωg為地表水地下水間溶質(zhì)的交換量，[ML-3T-1]；Ωg為地下水溶質(zhì)的源匯項，[ML-3T-1]；n為孔隙度，無量綱；R為阻滯系數(shù)，無量綱；Γ1,Γ3為零通量紐曼邊界，Γ2,Γ4第三類邊界條件。

2.2.3 地表水-地下水溶質(zhì)運移模型耦合

地表水與地下水溶質(zhì)運移模型耦合方法還處于探索發(fā)展階段[7]，本文假設(shè)地表水地下水間的溶質(zhì)運移過程符合Fick定律,使用一階多項式描述地表水地下水間的溶質(zhì)運移過程如式(6)所示：

(6)

式中：Ωgs為地表水地下水的溶質(zhì)交換量，[ML-3T-1]；Γ0為上游節(jié)點流量，[T-1]；C為地下水溶質(zhì)濃度，[ML-3]；c為地表水溶質(zhì)濃度，[ML-3]。

3 模型的離散

本文在耦合地表水模型與地下水模型后，對耦合模型進行統(tǒng)一的時間空間離散。

3.1 時間間離散

地表水地下水聯(lián)合模擬的關(guān)鍵在于時間尺度[9]，采取自適應(yīng)時間步長方法對聯(lián)合模擬模型進行時間離散可以在滿足模型精度同時有效提高計算效率。

本次研究模擬期為120 d。自適應(yīng)時間步長法各項參數(shù)如下所示。初始時間0 s，初始時間步長0.5 s，單個時間步長內(nèi)允許最大水頭變化0.5 m，允許最大飽和度變化0.1，最小步長遞增倍數(shù)0.5，最大步長遞增倍數(shù)2，最大時間步長由經(jīng)驗公式[8](7)計算得到：

(7)

式中：S為貯水系數(shù)；a為流域邊長；T為導(dǎo)水系數(shù)。據(jù)此計算得最大時間步長為62 500 s。

3.2 空間離散

研究區(qū)空間離散使用三角網(wǎng)格剖分，平面上共剖分出989個節(jié)點，1 888個單元，網(wǎng)格平均分布。垂向上將研究區(qū)分為15層，為研究地表水與地下水的交互作用，將地表水地下水交換界面附近剖分得更精細(xì)。如圖3所示。

圖3 研究區(qū)空間離散圖Fig.3 Spatial discrete graph of study area

4 結(jié)果分析

抽水井以方案1方案2方案3抽水時，地表水地下水聯(lián)合模擬模型運行結(jié)果如下。

方案1:模擬期末刻，抽水井內(nèi)水位為7.04 m，污染物隨降水下滲進入地下水，下滲至地面以下4 m，垂向污染范圍約410.79 m2，如圖4(a)所示。地表水出流斷面徑流量789.54 m3/d，污染物濃度為0.0043 g/L。研究區(qū)地表水地下水補排關(guān)系為地表水補給地下水。

方案2：模擬期末刻，抽水井水位降至6.73 m。污染物下滲4.75 m，垂向污染范圍約430.29 m2，如圖4(b)所示。地表水出流斷面徑流量為747.37 m3/d，污染物濃度為0.004 g/L。研究區(qū)地表水地下水補排關(guān)系為地表水補給地下水。

方案3：模擬期末刻，抽水井水位為6.23 m，污染物下滲5.3 m，垂向污染范圍約472.52 m2，如圖4(c)所示。地表水出流斷面徑流量為705.39 m3/d，污染物濃度為0.003 8 g/L。研究區(qū)地表水地下水補排關(guān)系為地表水補給地下水。

圖4 各方案污染羽垂向分布圖Fig.4 Vertical distribution of pollution plume in each scheme

5 結(jié)論與建議

5.1 結(jié) 論

(1)模擬時段內(nèi)研究區(qū)地表水地下水補排關(guān)系為地表水補給地下水。隨抽水井抽水量增加，地下水水位持續(xù)下降，大量地表水下滲補給地下水。當(dāng)抽水井抽水量為43 m3/d時，河口地表徑流量減少5.34%，當(dāng)抽水井抽水量為86 m3/d時，河口地表徑流量減少10.66%。

(2)地下水的開采會導(dǎo)致地下水污染加重。隨抽水量增加，地表水地下水水頭差增大，地表水下滲量增大，地表水污染物通過對流不斷向地下水下滲。抽水井抽水量為43 m3/d時，地下水垂向污染范圍增大4.14%。抽水井抽水量為86 m3/d時，地下水垂向污染范圍增大15.03%。

5.2 建 議

(1)地表水-地下水聯(lián)合模擬還處于發(fā)展階段，但其應(yīng)用前景極為廣闊。可結(jié)合3S技術(shù)研究水循環(huán)過程，也可針對生態(tài)基流，水庫最小下泄流量，水資源評價中地表水地下水資源重復(fù)量等專門問題進行研究。

(2)地表水-地下水聯(lián)合模擬模型對水文地質(zhì)資料與地表水水力資料要求高。可利用靈敏度分析性分析等方法確定敏感度較大參數(shù)，簡化非敏感參數(shù)，從而達到降低資料要求同時保證模擬精度的要求。

(3)本文在運用地下水-地表水聯(lián)合模擬的方法，研究地表水與地下水的水質(zhì)水量的變化。但由于技術(shù)上的不成熟，僅考慮了溶質(zhì)的對流彌散作用，未考慮可能發(fā)生的化學(xué)反應(yīng)及生物化學(xué)反應(yīng)，致使模型與實際系統(tǒng)仍存在一定誤差。

□

[1] Maxwell R M, Miller N L. Development of a coupled land surface and groundwater model [J]. Journal of Hydrometeorology, 2005,(6):233-247.

[2] Pinder G F, Sauer S P. Numerical simulation of flood wave modification due to bank storage effects [J]. Water Resources Research, 1971, 7(1):63-70.

[3] Maxwell R M, Miller N L. Development of a coupled land surface and groundwater model [J]. Journal of Hydrometeorology, 2005,(6):233-247.

[4] 王中根,朱新軍,李 尉,等.海河流域地表水與地下水耦合模擬[J].地理科學(xué)進展,2011,(11):1345-1353.

[5] 王建華,肖偉華,王浩等.變化環(huán)境下河流水量水質(zhì)聯(lián)合模擬與評價[J].科學(xué)通報,2013,58(12):1 101-1 108.

[6] 凌敏華,陳 喜,程勤波,等.地表水與地下水耦合模型研究進展[J].水利水電科技進展,2010,(4):79-84.

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

[8] 薛禹群,謝春紅.地下水?dāng)?shù)值模擬[M].北京:科學(xué)出版社,2007:248-253.

[9] 鄧 潔,魏加華,邵景力. 河渠與地下水相互轉(zhuǎn)化耦合模型研究進展[J]. 南水北調(diào)與水利科技,2008,(2):75-79.