基于水文集成模型ParFlow的黑河流域下游地下水-地表水相互作用模擬研究

陸 崢,胡錦華,張 圓,李宗超，楊 晨,楊曉帆*

(1.北京師范大學地理科學學部地表過程與資源生態(tài)國家重點實驗室，北京 100875；2.北京師范大學地理科學學部自然資源學院，北京 100875；3.中國環(huán)境監(jiān)測總站，北京 100012；4.普林斯頓大學高草甸環(huán)境研究所土木與環(huán)境工程系，美國 普林斯頓 08544)

地下水和地表水是區(qū)域和全球水文循環(huán)的基本要素，它們在溫度和化學組分上存在著普遍性的差異；且兩者具有緊密耦合的水力聯(lián)系和頻繁的轉(zhuǎn)化關系，共同構成了一個完整、復雜、密不可分的系統(tǒng)。因此，定量研究地下水與地表水的相互作用規(guī)律對流域水文和生態(tài)環(huán)境科學具有重要意義。水文集成模型耦合了地下水與地表水模型，能夠較為準確地描述自然界中水循環(huán)變化的復雜驅(qū)動因素，從而描述各種驅(qū)動因素對地下水與地表水水文循環(huán)的影響機理。然而，由于自然系統(tǒng)存在高度的空間異質(zhì)性，同時地下水-地表水相互轉(zhuǎn)換的同時還伴隨著能量轉(zhuǎn)化和溶質(zhì)運移，因而定量描述其物理過程及與之耦合的生物地球化學過程具有極大的挑戰(zhàn)性。另外，針對不同尺度所采用的不同求解方案也會帶來難以量化的不確定性。為了解決上述問題，需要建立和發(fā)展基于物理過程的、高時空分辨率的水文集成模型。但是，當前的水文集成模型依舊受到耦合方式、網(wǎng)格和離散化方案、驅(qū)動和校驗數(shù)據(jù)質(zhì)量、計算效率等因素的制約，因此評估地下水-地表水耦合模型在上述幾方面的表現(xiàn)尤為重要。通常，在解決具體的流域?qū)嶋H問題前需要使用簡單的二維剖面算例對模型進行測試、驗證和評價。典型的水文地質(zhì)剖面數(shù)值模擬研究可以定量描述局部地表徑流和地下水的側(cè)向及垂向運動，從而獲得精確的二維地下水動力場，對于揭示區(qū)域地下水與地表水之間的轉(zhuǎn)換關系和水文循環(huán)規(guī)律特征具有重要意義。在二維剖面數(shù)值模擬算例中，常常需將實際問題進行簡化、抽象，即建立空間結構更加簡單的幾何概念模型(如圖1所示，將真實流域水文循環(huán)概念模型簡化為理想化的“V”字形雙傾斜坡面流域模型)，并簡化、提煉相應的物理過程[如圖1(a)中的降水、蒸散發(fā)、植被冠層截留、下滲、地下水側(cè)向流等]，最終利用解析解或?qū)嶒炗^測對目標模型模擬的結果進行驗證和評估。尤其在評估地下水-地表水耦合模型時，需基于簡化的水文地質(zhì)剖面結構，驗證其數(shù)值解的準確性并量化其不確定性，評價其水文響應模擬結果的時空分布演變特征和統(tǒng)計特征。為此，本文使用基于物理過程的、地下水-地表水緊密耦合的水文集成模型ParFlow，依據(jù)不同案例設置，定量描述了地表水與地下水的相互作用，模擬計算了不同情境下地下水-地表水交互的水文響應。

圖1 真實流域水文循環(huán)的簡化建模方案示意圖Fig.1 Diagram of simplified model of the hydrologic cycles in a reale-world watershed

1 研究方法

1.1 ParFlow模型

本研究中所采用的水文集成模型ParFlow (https://www.parflow.org/) 是一款開源的、可并行計算的、面向?qū)ο蟮娜S地下水-地表水耦合模型。ParFlow模型基于三維地下水流動Richards方程和二維地表水流動Saint-Venant方程的二階動力波(Kinematic wave)或擴散波(Diffusive wave)近似解，將飽水帶、包氣帶和地表水視為一個完整的水文連續(xù)體，以實現(xiàn)全面刻畫不同尺度(例如坡面、子流域、流域、國家/區(qū)域、大陸和全球尺度等)下的水文循環(huán)過程。ParFlow模型的主要特征見表1，更詳細的求解處理方案和功能參見相關文獻。

表1 ParFlow模型的主要特征Table 1 Summary and specification of ParFlow

圖2為ParFlow (Parallel Flow)模型在美國俄克拉荷馬州Little Washita流域的一個應用案例。在Little Washita流域的一系列研究中，首先進行了簡單的基準測試算例研究，隨后循序漸進地擴展到二維剖面的各類研究，最終實現(xiàn)了不同尺度下的流域水文循環(huán)三維模擬。秉承上述理念和邏輯遞進關系，本研究選擇黑河流域下游的典型二維剖面，建立了高精度的地下水-地表水耦合模型，驗證了模型的模擬精度，并定量分析了水文響應特征，為后續(xù)的流域三維地下水-地表水耦合模擬及其他相關工作奠定堅實的基礎。

圖2 ParFlow模型在美國Little Washita流域的地下水-地表水交互模擬應用案例示意圖Fig.2 Diagram of application of ParFlow to groundwater-surface water simulation case over the Little Washita,USA

ParFlow模型可跨系統(tǒng)適用于Linux/Unix/OSX，包含諸多不同的求解器，其中IMPES求解器適用于單相、完全飽和的情況，RICHARDS求解器適用于地下水變飽和的情況，用戶可以自主選擇是否耦合地表流動過程，因此ParFlow模型可用于求解飽和帶-中間帶-包氣帶-地表水流偏微分方程定解問題。本文主要介紹變飽和水流和地下水-地表水耦合的求解器，分別利用三維Richards方程和二維Saint-Venant方程來求解地下水流動過程和地表水流動過程，地下與地表的耦合主要是通過在地下水與地表水交界處的邊界條件來實現(xiàn)，使地下水的三維Richards方程和地表水的二維Saint-Venant方程中的壓力水頭保持一致，從而保證了地下水-地表水交界面上壓力的連續(xù)性。這種直接耦合的方法，避免了定義界面導水率(interface conductance)，并且降低了數(shù)值求解耦合系統(tǒng)過程中的難度。其中，地下水-地表水交界處的邊界條件可根據(jù)有無積水及積水深度在第一類邊界條件(Dirichlet型)和第二類邊界條件(Neumann型)間進行切換。對于變飽和地下水流動方程和地表水流動方程分別采用中心差分格式的有限差分法和標準迎風格式的有限體積法進行空間離散，時間離散方法均為隱式后向歐拉法。此外，利用高效穩(wěn)定的Newton-Krylov非線性求解器進行矩陣求解，所采用的網(wǎng)格為結構化規(guī)則網(wǎng)格，可選正交網(wǎng)格或地形網(wǎng)格。

1.2 控制方程

1.2.1 地下水流動方程

依據(jù)質(zhì)量、動量守恒定律，地下水流動(包括包氣帶、中間帶的壤中流和飽和帶基流)水量平衡方程采用的是三維Richards方程，可表示如下：

(

1

)

式中：

h

為壓力水頭(L)；

t

為時間(T)；

S

為比容(L)；

S

(h

)

為土壤相對飽和度(無量綱)；φ為土壤孔隙度(無量綱)；

z

為垂直坐標(向上為正)；

q

和

q

為源匯項(L/T，具體的介紹見下文)；為達西通量(L/T)，可表示為

(

2

)

其中：為土壤水力傳導率(L/T)；

K

為土壤相對滲透率(無量綱)。

土壤相對飽和度與土壤相對滲透率之間的關系可通過van Genuchten模型獲得，其計算公式分別如下：

(

3

)

(

4

)

式中：

α

為孔隙進氣值參數(shù)

(

L

);n

為土壤孔隙分布參數(shù)(無量綱)；

S

為土壤相對飽和含水量(無量綱)；

S

為土壤相對殘余飽和度(無量綱)。

1.2.2 地表水流動方程

依據(jù)質(zhì)量守恒定律，地表淺水流動水量平衡方程采用的是二維Saint-Venant方程，可表示如下：

(5)

式中:

h

為積水深度(L);為地表流速(L/T)

;q

為與地下水-地表水的交換通量(L/T)

;q

為降水通量(L/T)。

在ParFlow模型中，根據(jù)地表水動力波近似，若忽略動量平衡方程中的局部慣性項、對流慣性項和壓力差項,則二維Saint-Venant方程可簡化為

g(S

-

S

)=

0

(

6

)

式中：

g

為重力加速度(L/T)；

S

為坡面斜率(無量綱)，其中

i

為

x

或

y

軸方向；

S

為摩擦斜率(無量綱)，可表示為

(

7

)

此外，根據(jù)Manning-Strickler方程可以建立地表流速與積水深度

h

之間的關系式為

(

8

)

式中：

N

為曼寧粗糙度

(

T/L

)

。

1.2.3 地下水與地表水的耦合方法

地下水與地表水的耦合通過自由表面邊界條件來實現(xiàn)。由于需要保證地下水-地表水邊界處壓力水頭和通量連續(xù)，則地下水的壓力水頭

h

和地表水的積水深度

h

均等于地下水-地表水邊界處的垂直平均壓力

p

，即

:

p=h

=h

(

9

)

地下水的上邊界流動通量與地下水-地表水邊界處的交換通量相同，即

:

q

=q

(

10

)

根據(jù)地下水動量方程

(

2

)

可知，地下水-地表水邊界處Neumann型邊界條件

q

[

若積水深度大于0，則地下水-地表水邊界處由第二類邊界條件(Neumann型)可改為第一類邊界條件(Dirichlet型)

]

為

(

11

)

地下水-地表水邊界處的交換通量為

q

=τ(h

-

h

)

(

12

)

式中：

τ

為邊界比例系數(shù)。

q

與公式(5)中地下水-地表水的交換通量相同，且可改寫為

(

13

)

式中

:

‖

p,

0‖算子結果返回值為

p

和0中的最大值。將上式結果代入公式

(

11

)

中，可得到：

(

14

)

將公式

(

14

)

左側(cè)項重新代回公式

(

1

)

中，可得到：

(

15

)

當

p

小于0時，地下水-地表水邊界處為第二類邊界條件(Neumann型)，地表水流動邊界方程關閉，公式(15)為標準的地下水三維Richards方程；當

p

大于0，也就是存在地表積水時，地下水-地表水邊界處由第二類邊界條件(Neumann型)轉(zhuǎn)化為第一類邊界條件(Dirichlet型)，地表水二維Saint-Venant方程被激活，則利用公式(14)計算。

2 研究區(qū)概況與模擬方法

2.1 研究區(qū)概況

黑河流域位于97°6′～102°0′E，37°30′～42°42′N間，東西、南北各橫跨約5°，發(fā)源于祁連山中段，北至中蒙邊境，東與石羊河流域接壤，西與疏勒河流域毗鄰，總面積為14.3×10km(見圖3)。該地區(qū)高程范圍為869～5544 m，氣候主要受中高緯度西風帶環(huán)流的控制和極地冷氣團的影響，氣候干燥，降水稀少而集中，多大風，日照充足，太陽輻射強烈，晝夜溫差大，年平均氣溫約-1℃，年均降水量為500 mm。按照氣候區(qū)劃，黑河流域總體上屬于半干旱高寒氣候。黑河流域上游地處高原亞寒帶亞干旱區(qū)，中下游為中溫帶干旱型氣候，獨特的氣候水文條件使之成為開展流域水文過程研究的理想?yún)^(qū)域。北京師范大學和中國科學院西北生態(tài)環(huán)境資源研究院于2012年9月起在黑河中游和下游地區(qū)共同布設了水文氣象觀測網(wǎng)，其中額濟納旗核心綠洲二道橋東至七道橋典型河岸林區(qū)域為黑河流域下游核心觀測區(qū)，區(qū)內(nèi)布設了5個水文氣象觀測站點。

圖3 黑河流域下游巴牙吉呼(BJH)地區(qū)二維剖面設置Fig.3 Configurations of the BJH case

黑河流域下游的巴牙吉呼(Bajajihu,BJH)地區(qū)在行政區(qū)劃上隸屬內(nèi)蒙古自治區(qū)阿拉善盟額濟納旗，其氣候類型屬于中溫帶極干旱區(qū)，陸表下墊面多為裸地、草地。下游額濟納綠洲是天然的綠洲生態(tài)系統(tǒng)，結構簡單且極度脆弱，屬于極干旱氣候區(qū)，多年平均降水量不足45 mm，地帶性植被多為溫帶小灌木、半灌木荒漠植被。該地區(qū)因降水稀少，植被生長主要依靠地下水，因此獲取區(qū)域水文循環(huán)特征與量化各要素之間的轉(zhuǎn)化關系極為重要。另外，雖然該地區(qū)地形較為平坦，但含水層水文地質(zhì)情況極其復雜，地下水流向(尤其側(cè)向流)尚不明確，因而在此處進行二維剖面的模擬意義深遠。同時在前期的研究中，干旱荒漠區(qū)地表水與地下水的轉(zhuǎn)化機制及轉(zhuǎn)化過程、淺層地下水的形成機理及運移方式等科學問題尚未有較完善的解決方案。

因此，本次在兩個水文氣象觀測站點之間構建了長為850 m的BJH地區(qū)典型二維剖面，旨在驗證ParFlow模型預測精度的同時，概算地表水-地下水交互過程及其中的水文響應。西北的胡楊林站海拔為876 m，東南的混合林站海拔為874 m，兩站點之間的海拔呈現(xiàn)單調(diào)遞減的趨勢，地表覆蓋為裸土并有塊狀的稀疏胡楊和檉柳。

2.2 模擬方法

2.2.1 二維剖面概述及模擬設置

BJH地區(qū)二維剖面的計算區(qū)域如圖3(a)所示，長為850 m，高為4 m，水平分辨率(即

x

方向)為10 m，垂直分辨率(即

z

方向)為0.01 m。本次模擬總時長設置為24 h(2015年9月3日)，時間步長設為10 min，與降水觀測的時間間隔相同；地下水水位原始數(shù)據(jù)時間間隔為30 min，采用線性內(nèi)插方法獲取時間間隔為10 min的地下水水位值；蒸散發(fā)值為日平均觀測值，包括裸地蒸發(fā)量和植被蒸騰量。ParFlow模型可以輸出多類地表水和地下水變量，諸如地表徑流、地表水儲量、包氣帶水儲量、飽和帶水儲量、網(wǎng)格中心水位等。取二維剖面兩側(cè)(對應胡楊林站表層土壤2 cm和4 cm)的土壤飽和度換算成土壤水分，并與時域反射儀(Time Domain Reflectometry,TDR)觀測值進行對比，同時分析了6個典型時刻下BJH地區(qū)二維剖面的土壤水分分布和4個不同水平位置處的土壤垂向飽和度分布曲線。

2.2.2 初始條件和邊界條件設置

根據(jù)BJH地區(qū)二維剖面兩端站點的觀測水位設置初始水位，采用線性內(nèi)插方法獲得剖面內(nèi)部每個網(wǎng)格的初始水頭。模擬區(qū)域的底部、前側(cè)和后側(cè)均設置為零通量邊界，左側(cè)和右側(cè)采用定水頭的邊界條件，其值取自剖面兩端的胡楊林站和混合林站的觀測水位。BJH剖面2015年9月3日的降水和地下水水位時間序列變化曲線，見圖4。

圖4 黑河流域下游BJH地區(qū)剖面2015年9月3日 的降水和地下水水位時間序列變化曲線Fig.4 Time-series of precipitation and water table depths on 3rd Sept.,2015 for the BJH transect

2.2.3 模擬參數(shù)設置

孔隙度、飽和導水率、曼寧粗糙度、van Genuchten模型參數(shù)均設置為均質(zhì)，其中孔隙度、飽和導水率和van Genuchten模型參數(shù)取自中國土壤數(shù)據(jù)集，曼寧粗糙度修改自文獻[32]，具體參數(shù)設置詳見表2。

表2 黑河流域下游BJH地區(qū)案例中的參數(shù)設置Table 2 Parameters in the BJH cases

3 結果與討論

ParFlow模型模擬的胡楊林站表層土壤水分(2 cm和4 cm)時間序列曲線與TDR觀測結果的對比，見圖5。

圖5 ParFlow模型模擬的胡楊林站表層土壤水分與 TDR觀測值的對比圖Fig.5 Surface soil moisture time-series comparison between ParFlow simulation and TDR observation at the Populus euphratica site

由圖5可見，利用ParFlow模型模擬得到的胡楊林站表層土壤水分模擬值與觀測結果基本一致，2 cm和4 cm深度的模擬值與觀測值的Pearson相關系數(shù)(

R

)分別為0.94和0.93，且均方根誤差(RMSD)較低(小于0.008 m/m，即小于整體動態(tài)變化范圍的10%)；降雨前后表層土壤水分ParFlow模型的模擬值與觀測值有較高的一致性，但在降雨期間，表層土壤水分ParFlow模型的模擬時間序列曲線與觀測值之間的增速和增幅有所差異：降雨開始后不久，表層土壤水分ParFlow模型的模擬值迅速上升到峰值，而觀測值上升速度相對緩慢，推測有兩個原因?qū)е铝诉@個現(xiàn)象，即土壤水分的“記憶效應”和胡楊林的冠層截留?？傮w而言，表層土壤水分的ParFlow模型模擬值與觀測值顯示出很好的契合度和相關性。

6個典型時刻下黑河流域下游巴牙吉呼(BJH)地區(qū)二維剖面的土壤水分分布，見圖6。此6個時刻分別為初始狀態(tài)(00∶00)、降水事件的起始點(06∶00)、降水強度最大時刻(11∶30)、降水事件結束(14∶30)、排水的中間時刻(17∶30)、模擬的最后一個時間步長(24∶00)。

圖6 6個典型時刻下黑河流域下游巴牙吉呼(BJH)地區(qū)二維剖面的土壤水分分布圖Fig.6 Sectional distribution plot of soil moisture in BJH transect at six typical moments in the downstream of the Heihe River Basin

由圖6可見，00∶00 與06∶00時刻相比，剖面區(qū)域頂部土壤水分差異明顯，說明降水開始后表層土壤水分模擬值迅速增加，至06∶00時下滲面與地表平行；隨著降水的持續(xù)進行，地表積水因重力作用沿斜坡向下滲透，并隨著降雨強度的持續(xù)增加而下滲加快，剖面頂部逐漸飽和，至11∶30時飽和深度達到1 m以上，此時模擬區(qū)域的最下層土壤和地表淺層土壤均呈現(xiàn)飽和狀態(tài)，中間1～3 m處呈現(xiàn)出一個垂直分布的非飽和過渡帶；臨近降雨結束，降雨速率下降，導致地表淺層土壤水分的飽和度減小,隨著降水事件的結束(即11∶30后)，中間飽和部分的水分繼續(xù)下滲，但因地形作用地下水開始橫向流動，側(cè)向流使得地下水水位不再與地表保持平行；14∶30后地下水水位開始從高海拔處向低海拔處傾斜，由于持續(xù)的蒸發(fā)作用，表層土壤開始干燥，至17∶30以后整個區(qū)域達到相對穩(wěn)定狀態(tài)，土壤水分剖面圖趨于穩(wěn)定。這6個典型時刻的土壤水分變化表明，ParFlow模型模擬能夠高時空精度地刻畫出BJH剖面上土壤水分運移的動力學過程，可較為精準地表征降水、蒸散發(fā)驅(qū)動下剖面內(nèi)的水文循環(huán)過程。

此外，為了進一步檢驗ParFlow模型對各時間步長下土壤水分遷移的模擬情況，本次還繪制了4個不同水平位置(

x

=1 m、28 m、56 m、85 m)處土壤垂向飽和度時間變化廓線(見圖7)，時間間隔為30 min，與模擬所用的時間步長相同，并用不同顏色表示從00∶00(紅色)到24∶00(紫色)時刻的土壤瞬時垂向飽和度分布。從不同水平位置處土壤垂向飽和度分布曲線的變化情況可以看出，土壤垂向飽和度時間變化廓線受降雨速率、蒸發(fā)、地下水水位變化和地形的綜合影響，在各個因素的共同作用下不同水平位置處的土壤垂向飽和度產(chǎn)生了不同的響應。在降雨事件之前，4個水平位置處的土壤垂向飽和度時間變化廓線對蒸發(fā)的響應模式幾乎是相同的；降雨事件期間的土壤垂向飽和度時間變化廓線受降雨強度控制，在降雨事件停止后(從17∶30到24∶00時刻，見圖7中紫色的廓線)，土壤垂向飽和度時間變化廓線的分布開始穩(wěn)定；降雨事件后不久，由于蒸發(fā)作用，近地表處飽和消退，土壤垂向飽和度時間變化廓線的分布開始趨于穩(wěn)定。此外，在

x

=28 m和

x

=56 m處的土壤垂向飽和度時間變化廓線在各個時間步長上的表現(xiàn)近似，但在降水停止后的幾個時間步長內(nèi)差異較大，這表明除了降雨事件開始后的一小段時間外，這兩個位置處側(cè)向的排泄和補給幾乎相同。然而，由于地下水水位的變化，剖面邊緣

x

=1 m和

x

=85 m處的土壤垂向飽和度時間變化廓線表現(xiàn)出不同的分布模式，緣于地下水水位變化和地形的綜合影響。從以上4組土壤垂向飽和度時間變化廓線可見，ParFlow模型模擬的土壤垂向飽和度分布時間序列具有一定的準確性和敏感性，受地形、蒸散發(fā)等驅(qū)動因素控制而產(chǎn)生不同的水文響應，可以較好地表征實際水文動力學過程。

圖7 不同水平位置處土壤垂向飽和度分布曲線的對比Fig.7 Comparison of vertical saturation distribution curves of soil at different horizontal positions in different time steps

4 結論與展望

本文使用基于物理過程的、地下-地表緊密耦合的水文集成模型ParFlow，依據(jù)不同案例設置，定量描述了地表水與地下水的相互作用，模擬計算了不同情境下地下水-地表水交互的水文響應。從簡單到復雜的基準測試案例，再到簡化的實際問題，ParFlow模型均可較好地解決非線性的水文動力學問題，展現(xiàn)出對驅(qū)動因素較為精準的水文響應信息。具體的，在黑河流域下游的巴牙吉呼地區(qū)(BJH)剖面展開案例研究，獲取了小時內(nèi)時間步長的水文響應。通過對土壤水分和土壤飽和度模擬結果的驗證以及時空分布分析，表明ParFlow模型在模擬地下水-地表水交互方面具備良好的時空準確度。

然而，目前的ParFlow模型仍存在許多缺陷和不確定性，例如水文響應信號依然受到模型驅(qū)動數(shù)據(jù)和空間異質(zhì)性的影響，同時尺度效應也會影響模擬結果的代表性。此外，ParFlow亦可以耦合陸面過程模型(CLM)和大氣模型(WRF或ARPS)，目前均已有較為成熟的版本，但其應用仍有待開發(fā)和利用。未來水文模型應與生物地球化學循環(huán)、植被動力學模型、溶質(zhì)運移模型相結合，以更加全面地刻畫自然界中真實的水文過程。