張將偉，盧文喜，陳 末，林 琳，李久輝，崔尚進(jìn)

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

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

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

1 假想例簡(jiǎn)介

建立了一個(gè)南北長(zhǎng)500 m，東西寬500 m的河谷假想例作為研究區(qū)，如圖1與圖1.2所示。研究區(qū)中部存在一河流，流向?yàn)楸敝聊狭鲃?dòng)，經(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)驗(yàn)值如表1。地表水系統(tǒng)源匯項(xiàng)主要有河口流入量、大氣降水補(bǔ)給、地表水體蒸發(fā)、地表水下滲補(bǔ)量以及河口流出量。

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

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

分區(qū)曼寧粗糙系數(shù)洼地儲(chǔ)量/m消減儲(chǔ)量/m耦合長(zhǎng)度/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)源匯項(xiàng)有北部邊界流入量、大氣降水、地表水下滲量、抽水井抽水量及南部邊界流出量。

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

研究區(qū)南部存在露天垃圾堆放場(chǎng)，污染物以氮元素污染為主。設(shè)置一東西向監(jiān)測(cè)線橫穿垃圾場(chǎng)，監(jiān)測(cè)污染物下滲情況(如圖2中A-A'所示)。場(chǎng)地內(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é)模擬模型

地表水的運(yùn)動(dòng)采用二維平均深度圣維南(Saint Venant)方程來(lái)描述，如模型(1)：

(1)

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

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

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

(2)

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

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

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

(3)

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

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

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

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

(4)

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

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

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

(5)

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

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

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

(6)

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

3 模型的離散

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

3.1 時(shí)間間離散

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

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

(7)

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

3.2 空間離散

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

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

4 結(jié)果分析

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

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

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

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

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

5 結(jié)論與建議

5.1 結(jié) 論

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

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

5.2 建 議

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

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

(3)本文在運(yùn)用地下水-地表水聯(lián)合模擬的方法，研究地表水與地下水的水質(zhì)水量的變化。但由于技術(shù)上的不成熟，僅考慮了溶質(zhì)的對(duì)流彌散作用，未考慮可能發(fā)生的化學(xué)反應(yīng)及生物化學(xué)反應(yīng)，致使模型與實(shí)際系統(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é)進(jìn)展,2011,(11):1345-1353.

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

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

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

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

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