滕飛達，張 博，梁秀娟，李明乾，楊偉飛，張芷豪

（1．吉林大學(xué)地下水資源與環(huán)境教育部重點實驗室，吉林 長春130021；2．油頁巖地下原位轉(zhuǎn)化與鉆采技術(shù)國家地方聯(lián)合工程實驗室，吉林長春130021；3．吉林大學(xué)新能源與環(huán)境學(xué)院，吉林長春130021）

傍河水源地含水層通常具有富水性強、地下水位埋藏淺、水量穩(wěn)定、水質(zhì)優(yōu)良、易集中開采和便于管理等優(yōu)點，是地下水開發(fā)常見類型［1-2］。河流和地下水之間存在著密切的水力聯(lián)系和頻繁的轉(zhuǎn)化［3］，地下水資源實際上有相當(dāng)大的一部分是襲奪了河水資源，地下水可開采量的大小和河流的最枯流量、含水層厚度和導(dǎo)水性有關(guān)［4］。傍河取水能夠激發(fā)并增大河水對含水層的入滲補給［5］，以此實現(xiàn)地表水與地下水資源的聯(lián)合利用。國內(nèi)外很多學(xué)者采用地表水與地下水耦合模型，研究地下水與地表水的交互關(guān)系。Jolly等［6］指出濕地地下水和地表水交互作用受地表水水位和地下水水位之差、局部地形和地下水流場等因素的影響。丁元芳等［7］應(yīng)用MODFLOW對渾河李官堡傍河地下水源地流場進行數(shù)值模擬，將渾河概化為一類定水頭邊界，試算不同的開采方案模擬出相應(yīng)的流場。王麗雪［8］對第二松花江流域傍河水源地進行地下水?dāng)?shù)值模擬，將第二松花江及其支流飲馬河與伊通河均概化為河流邊界。張奇瑩等［9］將神木縣馬鎮(zhèn)合河村境內(nèi)的黃河概化為河流邊界、上下游邊界概化為一類定水頭邊界，模擬地下水資源傍河取水開采方案。國內(nèi)研究通常將河流概化為一類定水頭邊界或河流邊界處理，通過模擬出的流場設(shè)計開采方案，沒有充分考慮枯、豐水期河流水位與流量的差異及水均衡中地表水與地下水轉(zhuǎn)化量對開采方案的影響。為彌補以往研究的不足，筆者應(yīng)用SFR模型概化白馬河-吉利河，建立地表水與地下水耦合模型，結(jié)合地表水與地下水轉(zhuǎn)化量及流場設(shè)計開采量方案。

1 研究區(qū)概況

白馬河-吉利河流域位于青島市，屬暖溫帶沿海濕潤季風(fēng)區(qū)［10］，多年平均氣溫為12．5℃，年平均相對濕度為74%，年平均降水量為 830 mm［11］。7—9月為豐水期，4—6月、10月為平水期，11月—次年3月為枯水期，豐水期降水量占全年降水量的60%以上，河水豐、枯水期流量差別巨大，在枯水年枯水期會出現(xiàn)斷流現(xiàn)象。吉利河干流長39．8 km，流域面積285．1 km2；白馬河干流長 44．2 km，流域面積 262．1 km2。 研究區(qū)位于白馬河-吉利河流域中下游，總面積為148．5 km2。河流北、東、西面為山地，兩河河谷成Y形。

白馬河第四系松散巖類含水層組按成因可劃分為沖洪積層、砂礫石層，坡積-坡洪積黏質(zhì)砂土夾碎石層、沖積-海積砂層。地下水富水性分區(qū)見圖1，地下水補給途徑為大氣降水、灌溉入滲和河流入滲補給?；鶐r區(qū)接受大氣降水入滲補給后，地下水由兩側(cè)向河谷流動，進入平原。地下水運動方向大體與地表徑流方向一致，兩河成為主要排泄帶。地下水排泄方式主要為人工開采、潛水蒸發(fā)和徑流排泄。地下水動態(tài)季節(jié)性變化明顯，6月初地下水位下降至最低，7—9月豐水期水位回升速度較快。豐水年水位回升，枯水年水位下降。

圖1 地下水富水性分區(qū)

2 研究方法

SFR（Streamflow-Routing）模型是美國地質(zhì)勘探局（USGS）開發(fā)的河床匯流模型［12］，原理為省略加速項的擴散波方程。河流邊界模型也是由USGS開發(fā)的，通過設(shè)置上下游河流水位、河床高程及水力傳導(dǎo)系數(shù)等參數(shù)概化河流，與SFR模型的共同點都是通過河床水力傳導(dǎo)實現(xiàn)地下水與河水之間的轉(zhuǎn)換，但是它不具備匯流功能，其水位均由觀測獲得，需要大量長時間的觀測數(shù)據(jù)，不如SFR模型靈活。

白馬河-吉利河中下游地區(qū)取水來源以地下水為主、地表水為輔，本文綜合考慮河水與含水層飽和帶水的相互聯(lián)系，根據(jù)河床與地下水之間的河床水力傳導(dǎo)系數(shù)、河床厚度、水深及斷面形狀等參數(shù)，將SFR模型與MODFLOW模型進行耦合，計算豐、平、枯水期河流流量與上下游水位，使河流概化更符合實際情況、河水與地下水轉(zhuǎn)化量計算更加精確，提高水源地地下水資源量的計算精度；再分別將白馬河-吉利河概化為第一類定水頭邊界模型、河流邊界模型，計算3種河流概化情形下的地下水可開采量，并與實際勘查的地下水最大可開采量進行對比分析。

2．1 地下水?dāng)?shù)學(xué)模型

根據(jù)水文地質(zhì)概念模型，建立的地下水?dāng)?shù)學(xué)模型為

式中：h為地下水頭，m；kx、ky、kz分別為x、y、z方向的滲透系數(shù)，m／d；h1為含水層第一類邊界水頭，m；h0為含水層初始水頭，m；W為開采強度，m3／d；Σ1為含水層第一類邊界；D為研究區(qū)；μ為給水系數(shù)；t為時間，d。

2．2 SFR地表水與地下水耦合模型

根據(jù)白馬河-吉利河水文資料構(gòu)建SFR耦合模型。河水一方面接受地表匯流的補給，一方面與地下水之間發(fā)生交換。河床與地下含水層之間的匯流表示為q地表匯流和q地下匯流兩部分。 為簡化計算，省略河水的加速運動，選用曼寧公式描述河水運動，于是圣維南方程組可改寫為

式中：B為水面寬，m；q地表匯流、q地下匯流分別為地表匯流與地下匯流的流量，m3／d；S0為河床縱比降；Sf為摩阻比降；Q為上游來水流量，m3／d。

方程的源匯項q地表匯流和q地下匯流分別由地表降水徑流和地下匯流計算得出。地下水?dāng)?shù)學(xué)模型的邊界條件是河流擴散波方程提供的水位和流量，而河流擴散波方程的源匯項又是由地下水流微分方程計算得到的。

2．3 河水與地下水交換模式

河水與地下水交換分為3種模式：一種是地下水位高于河床高程，另一種是地下水位略低于河床高程，這兩種均屬于飽和滲流機制，河床垂向和側(cè)向滲流量是地下水與河水水頭差的函數(shù)；第三種是地下水位遠低于河床高程，二者之間形成非飽和滲流。無論是哪種模式，河床與含水層的水位差決定著流入和流出量。河道與地下水的非穩(wěn)定流耦合模型是基于對每一時間步長的河床矩陣與地下水矩陣進行耦合求解，從而提高了交互計算的性能。把地下水流數(shù)值方程與擴散波方程聯(lián)立求解，即得到q地下匯流。河床水力傳導(dǎo)系數(shù)計算公式為

式中：C為河床水力傳導(dǎo)系數(shù)；k為河床滲透系數(shù)，m／d；w為河床寬度，m；T為河床厚度，m；L為上下游距離，m。

2．4 所用資料

模型所用數(shù)據(jù)資料包括：①2017年（現(xiàn)狀年）流量監(jiān)測及降水資料（見圖2）；②現(xiàn)狀年實測地下水位資料，地下水位監(jiān)測點（25處）分布于研究區(qū)全境，流量監(jiān)測點位于中游（見圖3）；③上下游河床形狀概化資料（見圖 4）。

圖2 現(xiàn)狀年月降水量、流量過程線

圖3 地下水位及流量監(jiān)測點分布

圖4 河道形狀概化

3 模型構(gòu)建與處理

3．1 水文地質(zhì)概念模型構(gòu)建

根據(jù)研究區(qū)的水文地質(zhì)條件與水動力場，確定概念模型。①含水層概化：把含水層概化為3層，第一層為透水性稍弱的細砂層，第二層為透水性很強的粗砂礫石層，第三層為基巖（見圖5）；②水動力條件概化：將研究區(qū)概化為非均質(zhì)各向同性的含水系統(tǒng)，地下水運移符合達西定律，水位隨時間變化，地下水流可概化為三維達西流；③邊界條件概化：白馬河-吉利河下游連通海洋，設(shè)置為定水頭邊界，入?？谒怀Ｄ攴€(wěn)定，高程為0，概化為一類水頭邊界，河谷東部、西部和北部均為山地，將地下水分水嶺概化為零流量邊界；④源匯項概化：將河流邊界分別概化為SFR耦合模型、河流邊界模型和第一類定水頭邊界模型3種邊界。大氣降水和灌溉回滲作面狀入滲處理；通過ET模塊設(shè)定蒸發(fā)極限深度和蒸發(fā)速率，計算蒸發(fā)量；水源地開采排泄作點狀和面狀排泄處理。

圖5 研究區(qū)模型概化

3．2 河流動力學(xué)概念模型構(gòu)建

根據(jù)現(xiàn)場踏勘，將上、中、下游河床概化為變寬度的梯形河床，河床形狀和坡度按照DEM確定。將白馬河-吉利河概化為SFR耦合模型，設(shè)置現(xiàn)狀年流量、上下游水位、河床滲透系數(shù)、河床厚度、水深及斷面形狀等，使河流概化更符合實際情況，再將河流概化為第一類定水頭邊界和河流邊界，對比分析3種不同河流概化條件下模擬出的地表水與地下水轉(zhuǎn)化量、允許開采的地下水量的差異，以判定各河流概化方法的優(yōu)劣。

3．3 模型識別與驗證

利用SFR耦合模型模擬研究區(qū)枯水期地下水位，并與實際枯水期水位進行對比。選用現(xiàn)狀年2017年1—3月、7—9月兩個時段分別進行模型的識別與驗證（見圖6、圖7）。模型識別期計算水位與實測水位擬合誤差小于0．5 m的數(shù)據(jù)點占已知水位數(shù)據(jù)點數(shù)量的75%以上，模型驗證期豐水期流場與實測流場非常接近，表明模型具有較高的精度和可信度。

圖6 模型識別期實測水位與計算流場

圖7 模型驗證期實測水位與計算流場

3．4 傍河井點開采量設(shè)置

沿河兩岸布置130口開采井，當(dāng)模型初始地下水位下降含水層2／3厚度時即為該研究區(qū)地下水可開采量，對3種模型設(shè)置相同位置的開采井，計算SFR耦合模型、河流邊界模型及第一類定水頭邊界模型條件下的地下水開采量，對以上3種情況模擬地下水位下降含水層2／3厚度時的地下水最大允許開采量。通過計算，SFR耦合模型、河流邊界模型、第一類定水頭邊界模型模擬的豐、平、枯水期每口開采井開采量分別為700～900、700～800、900～1000 m3／d。

4 討 論

采用傳統(tǒng)的地下水?dāng)?shù)值模型進行傍河地下水開采量評價，可能會出現(xiàn)一定的偏差。與傳統(tǒng)的地下水?dāng)?shù)值模型相比，SFR耦合模型考慮了河流豐、枯水期流量變化對地下水開采的影響。為了驗證SFR耦合模型的優(yōu)越性，采用河流邊界模型和第一類定水頭邊界模型建立研究區(qū)傍河開采地下水?dāng)?shù)值模型，與SFR耦合模型的水均衡計算結(jié)果進行對比（見表1）。

表1 水均衡計算結(jié)果萬 m3／d

4．1 SFR耦合模型與定水頭邊界模型對比分析

4．1．1 地表水與地下水轉(zhuǎn)化量

從表1可知：豐水期SFR耦合模型、第一類定水頭邊界模型計算的河水補給地下水量分別為3．08萬、7．55萬 m3／d，地下水轉(zhuǎn)化為河水的量分別為 2．66 萬、5．67萬m3／d；平水期SFR耦合模型、第一類定水頭邊界模型計算的河水補給地下水量分別為 2．39萬、9．46萬 m3／d，地下水轉(zhuǎn)化為河水的量分別為 2．30 萬、5．78萬m3／d；枯水期SFR耦合模型、第一類定水頭邊界模型計算的河水補給地下水量分別為2．27萬、9．99萬m3／d，地下水轉(zhuǎn)化為河水的量分別為 2．40 萬、4．72 萬m3／d。

由表1可知，第一類定水頭邊界模型計算的河流補給量是最大的，且枯水期河水對地下水的補給量最大，豐、平、枯水期都是河流補給量大于地下水轉(zhuǎn)化為河水的量，即地表水補給地下水。SFR耦合模型計算的地表水與地下水轉(zhuǎn)化較為均衡，豐水期河流補給量略大于地下水轉(zhuǎn)化為河水的量，為地表水補給地下水；平水期河流補給量與地下水轉(zhuǎn)化為河水的量基本持平；枯水期河流補給量略小于地下水轉(zhuǎn)化為河水的量，為地下水補給地表水。定水頭模型計算的枯水期河水補給地下水的量大于豐水期河水補給地下水的量，是“無限供水的邊界”，不符合白馬河-吉利河季節(jié)性河流特征。綜上，SFR耦合模型相對于定水頭邊界模型更能準(zhǔn)確計算出地表水與地下水的轉(zhuǎn)化量，計算結(jié)果更為合理。

4．1．2 模擬開采量

由表1可知，豐水期SFR耦合模型、第一類定水頭邊界模型計算的可開采量分別為8．39萬、10．82萬m3／d，平水期SFR耦合模型、第一類定水頭邊界模型計算的可開采量分別為 7．4 萬、9．24 萬 m3／d，枯水期SFR耦合模型、第一類定水頭邊界模型計算的可開采量分別為 7．05 萬、8．8 萬 m3／d。 由表 1 可知，第一類定水頭邊界模型模擬的傍河開采量比SFR耦合模型的大，第一類定水頭邊界模型計算的最大允許開采量為3 511萬m3／a，SFR耦合模型計算的最大開采量為2 779萬m3／a，參照《青島市膠南縣吉利河-白馬河地區(qū)供水水文地質(zhì)勘察報告》及《白馬河-吉利河地下水模型及脆弱性安全性評價》等以往研究資料，研究區(qū)平水年地下水可開采量為2 870萬m3／a，SFR耦合模型比第一類定水頭邊界模型計算結(jié)果更符合實際。

4．2 SFR耦合模型與河流邊界模型對比分析

4．2．1 地表水與地下水轉(zhuǎn)化量

由表1可知，豐水期SFR耦合模型、河流邊界模型計算的河水補給地下水量分別為3．08萬、2．83萬m3／d；平水期SFR耦合模型、河流邊界模型計算的河流補給地下水量分別為 2．39 萬、1．77 萬 m3／d，地下水轉(zhuǎn)化為河水的量分別為 2．3 萬、3．57 萬 m3／d；枯水期SFR耦合模型、河流邊界模型河計算的水補給地下水量分別為 2．27 萬、1．42 萬 m3／d，河水補給地下水量分別為 2．40 萬、2．85萬 m3／d。 河流邊界模型計算的豐、平、枯水期地下水轉(zhuǎn)化為河水的量與SFR耦合模型的相當(dāng)，但地下水轉(zhuǎn)化為河水的量較大，即地下水補給河水。與河流邊界模型相比，SFR耦合模型具備匯流功能，地表水與地下水轉(zhuǎn)化量更為均衡，與研究區(qū)豐、平、枯水期地表水與地下水實際轉(zhuǎn)化關(guān)系更為相符。

4．2．2 模擬開采量

由表1可知，豐水期SFR耦合模型、河流邊界模型計算的可開采量分別為 8．39 萬、7．40 萬 m3／d，平水期SFR耦合模型、河流邊界模型計算的可開采量分別為 7．40 萬、7．00 萬 m3／d，枯水期 SFR 耦合模型、河流邊界模型計算的可開采量分別為 7．05萬、6．77萬m3／d。綜上，SFR耦合模型計算的開采量在豐、平、枯水期均略大于河流邊界模型的，SFR耦合模型計算的開采量為2 779萬m3／a，河流邊界模型計算的開采量為2 576萬m3／a（遠小于研究區(qū)平水年地下水可開采量2 870萬m3／a）。與河流邊界模型相比，SFR耦合模型計算結(jié)果更符合實際。

5 結(jié) 論

本文綜合考慮白馬河-吉利河在豐、平、枯水期的流量差異，應(yīng)用SFR地表水與地下水耦合模型，使河流概化更符合實際情況，地表水與地下水的轉(zhuǎn)化量計算更加精確。SFR地表水與地下水耦合模型比河流邊界模型更能精細概化河道特征，計算的開采量更接近實際情況。在枯水期部分河段斷流時，無河水補給地下水，而第一類定水頭邊界模型為無限補給，使計算的地下水資源量明顯偏大，與實際不符。綜上，SFR地表水地下水耦合模型比傳統(tǒng)的地下水?dāng)?shù)值模型更適合于傍河地下水水源地地表水與地下水轉(zhuǎn)化量計算及開采量評價。