郜會彩，肖玉福，胡云進(jìn)，陳柳安，周如杰

（1.紹興文理學(xué)院土木工程學(xué)院,浙江 紹興 312000；2.浙江省巖石力學(xué)與地質(zhì)災(zāi)害重點實驗室,浙江 紹興 312000；3.浙江省山體地質(zhì)災(zāi)害防治協(xié)同創(chuàng)新中心,浙江 紹興 312000）

濕地生態(tài)系統(tǒng)通過其特殊的水文物理特性，以水循環(huán)為紐帶，削減地表徑流、降低流速和補(bǔ)給地下水發(fā)揮其水文調(diào)蓄功能[1]。濕地地表水與地下水之間的交換轉(zhuǎn)化，一直是水文地質(zhì)研究中的熱點問題[2]。以往國內(nèi)外對河流的地表水-地下水交換轉(zhuǎn)化研究較多，而對于水體相對靜止或流動較緩滯的湖泊和濕地研究相對偏少[3]。目前研究濕地地表水與地下水耦合作用的方法主要有實地調(diào)查法（直接測量法、示蹤劑法、水溫-熱量傳遞法、地球物理探測法）和數(shù)值模擬法[4]。其中濕地地表水與地下水耦合作用的數(shù)值模擬有利于理解濕地水文過程機(jī)理和預(yù)測水文情勢變化，目前濕地地表水與地下水耦合作用模型的研究大多借助已有的地表水與地下水耦合模型[5]。根據(jù)模型耦合方式不同，可分為松散耦合（“四水”轉(zhuǎn)化模型）、半松散耦合（SWATMD、GSFLOW 等）和緊密耦合（MIKE-SHE、IGSM、LL-II 等）[6]。一般濕地內(nèi)地形起伏較小，由于濕地植物的作用，造成濕地中水流速度緩慢，水力持留時間長[7]。濕地表層流復(fù)雜，周邊高地或其它水體以坡面漫流、片流、明渠流、漫灘流等多種形式進(jìn)入濕地，很難建立精確的數(shù)學(xué)模型去模擬[8]。馮夏清[9]在SWAT 模型基礎(chǔ)上重新構(gòu)建了濕地模塊來模擬扎龍濕地水文過程，實現(xiàn)了流域模型與濕地水文模型的松散耦合。Rahman 等[10]采用SWATrw 模型探討了Barak-Kushiyara 河洪泛濕地地表水與地下水之間的水量交換，結(jié)果表明水文特征是影響濕地地表水補(bǔ)給的主要因素?，F(xiàn)有模型中GSFLOW 可以將濕地作為獨特的水文響應(yīng)單元，重點考慮濕地水循環(huán)過程，模擬精度高、適用性強(qiáng)[11]。從已有研究來看，GSFLOW 已被應(yīng)用于流域地表水與地下水耦合作用研究和流域水資源配置研究[12-14]；為幫助用戶一站式構(gòu)建模型，田勇等[15-17]基于GSFLOW 開發(fā)了具有國際先進(jìn)水平的生態(tài)水文集成建模軟件系統(tǒng)Visual HEIFLOW（VHF），并以黑河流域和韓國Miho 流域為例，驗證了在大流域尺度下的適用性；但未見其在濕地地表水與地下水耦合模擬中的應(yīng)用。

濕地中有大量的表生植物和沉水植物，阻礙濕地地表水流動，使得濕地的物理結(jié)構(gòu)和化學(xué)環(huán)境更加復(fù)雜、異質(zhì)和多變[4]；并且濕地土壤是固相和液相組成的疏松多孔體系，多年植物生長沉積形成了獨特的海綿狀土壤層。因此，本文針對目前研究的不足，考慮濕地海綿層作用，對濕地明渠流和片流同時進(jìn)行模擬，利用GSFLOW 模型構(gòu)建濕地地表水與地下水耦合模型，并將其應(yīng)用于鏡湖濕地實例中，檢驗?zāi)Ｐ涂煽啃?，以期為?zhǔn)確評估濕地地表水與地下水耦合作用提供理論依據(jù)。

1 GSFLOW 模型原理

GSFLOW 模型耦合了PRMS 和MODFLOW-2005兩個模型，分別進(jìn)行濕地地表水文模擬和濕地地下水流模擬（圖1）。地表部分將濕地劃分成許多個水文響應(yīng)單元（HRU），地下部分則離散成許多個有限差分單元，在地表水與地下水相連的河段，建立重力水庫（GRV）將地表部分的HRU 與其對應(yīng)的地下部分有限差分單元連接起來[14]，利用達(dá)西定律，根據(jù)水頭差計算地表水與含水層的水交換量。

圖1 GSFLOW 結(jié)構(gòu)模型示意圖[13]Fig.1 Structural framework of the GSFLOW model [13]

1.1 地表產(chǎn)流模擬

對于內(nèi)陸濕地，降水首先被冠層截留，到達(dá)地面后一部分下滲到土壤成為地下水，其余的水沿著斜坡形成漫流，即其產(chǎn)流模式一般為超滲產(chǎn)流。GSFLOW模型將濕地分為透水區(qū)和不透水區(qū)分別進(jìn)行計算[13]，每個HRU 的透水面積和不透水面積比例需事先給定。在HRU 透水區(qū)域采用超滲產(chǎn)流機(jī)制計算地面徑流及水分下滲過程，其超滲產(chǎn)流利用貢獻(xiàn)區(qū)域的概念計算產(chǎn)流面積[14,18]。在不透水區(qū)域主要計算地表徑流量和蒸發(fā)量。

1.2 地表徑流模擬

由于絕大部分濕地地勢平坦，水流在平面大范圍的自由表面流動，水流在水平方向的運動尺度遠(yuǎn)大于垂向，所以假定水流在垂向均勻分布，屬于淺水流動。GSFLOW 中對地表水流運動的描述采用擴(kuò)散波近似的圣維南方程，其中用來計算水流在矩形或棱柱形河槽中的運動方程如下：

qstr——單位長度河道總的側(cè)向入流或出流量/（m2·d-1）；

xchannel——沿河道方向的河道長度/m；

t——時間/d。

對于明渠流，GSFLOW 采用曼寧公式計算流量：

式中：Q——河流流量/（m3·d-1）；

A——過水?dāng)嗝婷娣e/m2；

Cm——轉(zhuǎn)換常數(shù)，國際單位制中值為1；

R——水力半徑/m；

S——明渠的坡度；

n——糙率。

選用濕地植物綜合糙率公式[19]確定n，如下式所示：

式中：H——植物高度/cm。

其中沉水植物是濕地生態(tài)系統(tǒng)的重要組成部分，其生長特征和生物量積累隨濕地水位變化。通過對紹興鏡湖濕地現(xiàn)場調(diào)研測量，濕地沉水植物平均高度為3 cm，故n為0.066。

對濕地而言還普遍存在片流形式，對片流的模擬詳見下文2.1 節(jié)。

1.3 蒸散發(fā)模擬

蒸散發(fā)是濕地水循環(huán)的重要環(huán)節(jié)，也是區(qū)別于河流、湖泊水均衡的重要特點[5]。濕地地表蒸散發(fā)量是影響濕地水熱平衡的主要因素，也是能量和水分的主要消耗途徑。濕地植被供水良好，可視為植被實際蒸散發(fā)等于其潛在蒸散發(fā)。本文采用1998年世界糧農(nóng)組織（FAO）推薦的Penman-Monteith 公式計算潛在蒸散發(fā)[20]，該公式受季節(jié)、站點情況影響小，模擬效果好。公式中作物系數(shù)基于假定的參考作物，各項系數(shù)和變量結(jié)合利用氣象數(shù)據(jù)和文獻(xiàn)[21]計算得到。

1.4 地下水流模擬

濕地地表水與地下水交互滲流分為非飽和流與飽和流[5]。非飽和滲流常以多孔介質(zhì)流體達(dá)西定律與質(zhì)量守恒定律為基礎(chǔ)形成二相滲流運動的微分方程進(jìn)行模擬[22]。在GSFLOW 中非飽和滲流使用垂向一維Richards 方程運動波近似形式，公式如下：

式中：θ——含水率；

Lz——垂直運移距離/m；

K(θ)——水力滲透系數(shù)/（m·d-1）；

i——地下水蒸發(fā)速率/（m·d-1）。

飽和滲流過程將濕地地表水與地下水視為連續(xù)的整體，符合多孔介質(zhì)中潛水滲流運動的Boussinesq方程，如式（5）所示：

式中：Kxx——x軸方向水力滲透系數(shù)/（m·d-1）；

Kyy——y軸方向水力滲透系數(shù)/（m·d-1）；

Kzz——z軸方向水力滲透系數(shù)/（m·d-1）；

x、y、z——空間坐標(biāo)；

Ss——儲水率/m-1；

W——源匯項/d-1。

1.5 地下水的補(bǔ)給與排泄

濕地地表常年有積水，地表水是地下水的穩(wěn)定補(bǔ)給來源。濕地地表水與地下水的交互作用主要受濕地植被、地形坡度和濕地土壤的含水量的影響。在PRMS 模型中，土壤層分為毛細(xì)水庫、重力水庫和優(yōu)先流三個水庫[14]，這些水庫反映不同的土壤含水量（圖2）。地下水補(bǔ)給與排泄的水量通過重力水庫轉(zhuǎn)換。當(dāng)滿足時，地下水排泄到重力水庫，反之亦然。計算方程如式（6）所示：

圖2 PRMS 模型土壤層各水庫的輸入和輸出關(guān)系圖Fig.2 Relationship between the input and output of each reservoir in the soil layer of the PRMS model

CNDsz——重力水庫區(qū)域的導(dǎo)水系數(shù)/（m2·d-1）；

celtop——有限差分單元格的頂板海拔/m；

Dusz——土壤帶底部的地形起伏變化深度/m。

2 濕地片流和濕地海綿層模擬

2.1 濕地片流模擬

降水到達(dá)濕地表面后，部分水流充填洼地和地表下滲，大部分直接形成濕地地表徑流。明渠流和流過濃密植被的片流是濕地主要的地表徑流形式[23]（圖3）。在濕地植被茂密的區(qū)域，地表水流的水平速度較低。當(dāng)水位繼續(xù)升高，濕地表面完全被水淹沒，在一定的水力坡度下沿溝渠方向產(chǎn)生表面流，其速度一般比片流快。片流的方向和速度受地形坡度、水深、植被類型、降水、蒸散發(fā)等作用，其水文過程的模擬和計算相對復(fù)雜。本文采用圣維南二維擴(kuò)散方程[24]模擬濕地片流，見式（7）（8）：

圖3 濕地主要地表水流成分示意圖[24]Fig.3 Schematic map of the main surface water flow components in wetland[24]

式中：Sfx、Sfy——分別為x、y方向的摩擦坡度；

qx、qy——分別為x、y方向的片流通量/（m2·s-1）；

P——降水率/（m·s-1）；

I——入滲速率/（m·s-1）；

ET——片流蒸發(fā)速率/（m·s-1）；

Vw——其它水源或匯的流速/（m·s-1）。

2.2 濕地海綿層模擬

濕地植物根系與腐朽植物的分解很困難，容易在潮濕缺氧的條件下形成泥炭[25]。在水淹的條件下，濕地地表的泥炭蘚、植物凋落物不易腐爛，在濕地底部長期積累堆積形成一層柔軟疏松結(jié)構(gòu)，仿佛巨大厚重的海綿[26]，與泥炭層一起構(gòu)成了濕地海綿層（圖4、圖5）。該層突出特點是有較大的飽和含水量和持水量，孔隙較大，孔隙度72%～93%，既能蓄水又能透水[27]，在濕地地下水補(bǔ)給和排泄過程中具有重要意義。本文通過GSFLOW 的優(yōu)先流水庫模擬濕地海綿層，在設(shè)置優(yōu)先流閾值時，輸入濕地海綿層對應(yīng)的孔隙度。對于濕地海綿層滲透系數(shù)，本文通過野外現(xiàn)場豎管法（圖6）來進(jìn)行測定[28]，該方法簡便易行，精度高。在選取的試驗點上將直徑6 cm 的PVC 豎管垂直插入濕地，盡量不破壞濕地海綿層的原有結(jié)構(gòu)。濕地海綿層垂向滲透系數(shù)計算公式為：

圖4 濕地地表水與地下水相互作用地點及濕地海綿層位置 [4]Fig.4 Location of wetland surface water and groundwater interaction and location of the wetland spongy layer [4]

圖5 現(xiàn)場采集濕地海綿層樣本Fig.5 In situ collection of wetland spongy layer samples

圖6 豎管法測試濕地海綿層滲透系數(shù)示意圖Fig.6 Measurement of the hydraulic conductivity in the wetland spongy layer using straight standpipes

式中：Kv——垂向滲透系數(shù)/（m·d-1）；

Lv——管中濕地海綿層厚度/m；

h1、h2——分別為測管中t1、t2時刻水頭值/m。

測量結(jié)果如表1所示，經(jīng)過計算得到濕地海綿層滲透系數(shù)的范圍為15～45 m/d。

表1 濕地海綿層垂向滲透系數(shù)Table 1 Vertical hydraulic conductivities in the wetland spongy layer

3 應(yīng)用實例

3.1 研究區(qū)概況

研究區(qū)（紹興市鏡湖國家城市濕地公園）位于浙江省紹興市鏡湖新區(qū)，面積為5.13 km2。研究區(qū)數(shù)字高程模型（DEM）如圖7所示。區(qū)內(nèi)河流屬蕭紹曹運河水系，河湖密布，水網(wǎng)交織。研究區(qū)處于江山—紹興深斷裂東北端兩側(cè)，晚燕山早期白堊紀(jì)構(gòu)造沉積盆地內(nèi)[29]，區(qū)內(nèi)構(gòu)造被大面積的第四系松散堆積層所覆蓋，未見大的斷裂和皺褶構(gòu)造，巖石出露區(qū)主要發(fā)育節(jié)理裂隙。根據(jù)浙江省水文地質(zhì)勘察資料，將研究區(qū)含水層等厚劃分為三層：上層為濕地海綿層，厚度為0.5 m；中間為潛水含水層，厚度為3 m；下部為穩(wěn)定分布的粉質(zhì)黏土，厚度為45 m。

圖7 研究區(qū)DEM 圖Fig.7 DEM map of the study area

3.2 數(shù)據(jù)收集與整理

GSFLOW 模型的主要輸入數(shù)據(jù)包括：DEM 數(shù)據(jù)、土地利用數(shù)據(jù)、土壤類型、氣象驅(qū)動數(shù)據(jù)、水文數(shù)據(jù)及土壤屬性數(shù)據(jù)，各數(shù)據(jù)來源如表2所示。

表2 GSFLOW 模型所需數(shù)據(jù)類型Table 2 Data types for the GSFLOW model

借助ArcGIS10.2，基于研究區(qū)DEM 數(shù)據(jù)，劃分流域。將土地利用類型重分類，根據(jù)研究區(qū)現(xiàn)場情況概化為建筑用地、綠地、水域和山林地等四種；對土地利用類型進(jìn)行編碼賦值（圖8）。由于研究區(qū)域范圍較小，為便于區(qū)分將土壤類型概化為黃棕壤、高鈣土、褐土、潮土和其他等五種類型（圖9）。

圖8 研究區(qū)土地利用類型分布圖Fig.8 Zonation of the land use types in the study area

圖9 研究區(qū)土壤類型分布圖Fig.9 Zonation of soil types in the study area

3.3 數(shù)值模型構(gòu)建

根據(jù)研究區(qū)的地形、地貌、水文地質(zhì)概況進(jìn)行空間離散。按照DEM 圖（12.5 m×12.5 m）和河網(wǎng)水系進(jìn)行HRU 劃分，一共生成4 632 個HRU。本次模擬的地下部分主要是第四系潛水含水層，先給定一個初始滲透系數(shù)和補(bǔ)給率，經(jīng)過率定得到最終滲透系數(shù)和補(bǔ)給率。濕地地下水模擬采用MODFLOW-2005，對整個研究區(qū)進(jìn)行網(wǎng)格剖分，共分為5 599 個單元格。

研究區(qū)上邊界為地表面，接受降水入滲補(bǔ)給、河流補(bǔ)給、蒸發(fā)、排泄等水量交換，在GSFLOW 中處理為可自由移動界面。下邊界是凝灰質(zhì)砂巖基巖，處理為隔水邊界。對于模型的四周邊界，高低起伏較大的山坡處取至分水嶺，作為隔水邊界；多數(shù)邊界附近地勢平坦，為減少計算工作量，根據(jù)研究區(qū)DEM 圖、水域分布圖及實測水文情況，做如下處理：河道橫斷面較寬的河流，常年水位較恒定，作為定水頭邊界（圖10的Γ1）；河道橫斷面較窄的河流，常年流量較穩(wěn)定，作為定流量邊界（圖10 的Γ2）；其它周界由于水量交換較少，按隔水邊界處理（圖10 的Γ3）。

圖10 研究區(qū)邊界條件示意圖Fig.10 Schematic diagram of the boundary conditions in the study area

盡管研究區(qū)面積較小，但由于滲透系數(shù)是影響濕地地表水與地下水交互作用的重要因素，因此根據(jù)前人研究結(jié)果及《浙江省紹興市環(huán)境地質(zhì)調(diào)查評價報告》等，結(jié)合研究區(qū)水文地質(zhì)條件、鉆孔資料和工程地質(zhì)勘察資料，對地下水模型中水平滲透系數(shù)Kh、垂直滲透系數(shù)Kv和儲水率Ss進(jìn)行詳細(xì)分區(qū)（圖11）。表3 為各分區(qū)水文地質(zhì)參數(shù)取值范圍。

表3 研究區(qū)各分區(qū)水文地質(zhì)參數(shù)取值范圍Table 3 Range of hydrogeological parameters in each subdomain of the study area

圖11 研究區(qū)水文地質(zhì)參數(shù)分區(qū)示意圖Fig.11 Zonation of hydrogeological parameters in the study area

3.4 模型參數(shù)率定及驗證

以2005—2019年逐日降水量、最低和最高氣溫、蒸發(fā)量等氣象數(shù)據(jù)作為驅(qū)動數(shù)據(jù)，對模型參數(shù)進(jìn)行率定和驗證?？紤]到模型初始值對模擬結(jié)果的影響以及濕地現(xiàn)場水文資料的限制，將模擬期（2005年1月1日—2019年12月31日）劃分為兩個階段：階段1 從2005年1月1日—2014年12月31日，為濕地地表水與地下水耦合模型參數(shù)的率定期；階段2 從2015年1月1日—2019年12月31日，為模型參數(shù)的驗證期。在模型率定過程中，采用相關(guān)系數(shù)評價模型率定結(jié)果。

在模型率定過程中，受濕地獨特水力特征和氣候條件的影響，土壤表層飽和含水率、快速壤中流非線性系數(shù)、線性貢獻(xiàn)面積算法指數(shù)、HRU 不透水區(qū)最大滯留儲量等參數(shù)為濕地地表水與地下水耦合模型的敏感參數(shù)，在參數(shù)變動范圍內(nèi)調(diào)整參數(shù)值，來率定模型使其達(dá)到良好的模擬效果。利用率定好的參數(shù)，選取研究區(qū)內(nèi)地下水位對模型模擬結(jié)果進(jìn)行驗證（圖12）。驗證期濕地地表水與地下水耦合模型平均水位與模擬區(qū)水文站觀測地下水位的相關(guān)系數(shù)R2=0.671。將模擬結(jié)果與實測濕地河道水位進(jìn)行比對（圖13），采用均方根誤差（RMSE）來衡量模擬水位和實測水位兩個數(shù)據(jù)序列之間的誤差，計算可得RMSE=0.086，表明模擬水位與實測水位誤差很小。

圖12 濕地地下水水位對比圖Fig.12 Comparison between the observed and calculated groundwater levels in the wetland

圖13 濕地河道2019年模擬水位與實測水位對比圖Fig.13 Comparison between the observed and calculated surface water levels in the wetland

從以上率定和驗證結(jié)果可知，濕地植被茂盛、地勢起伏較小，片流和明渠流是濕地地表徑流的主要形式，考慮濕地多年植物生長沉積形成的獨特海綿層的調(diào)蓄和滲透作用，將濕地片流與明渠流同時進(jìn)行模擬是合理可靠的。

3.5 研究區(qū)水均衡分析

為進(jìn)一步揭示濕地地表水與地下水交互作用的動態(tài)特征，利用校正后的模型模擬長時間尺度下濕地的地表水與地下水水量交互過程，并對研究區(qū)進(jìn)行水均衡分析（圖14、表4）。

表4 研究區(qū)2005—2019年的年平均水均衡計算表Table 4 Average annual water balance of the study area from 2005 to 2019

圖14 研究區(qū)2005—2019年水均衡模擬示意圖Fig.14 Simulated water budgets of the study area from 2005 to 2019

由圖14 可知，研究區(qū)域水量輸入有降水量、地表水流入、地下水定水頭邊界入流；輸出項包括地表水流出、蒸散發(fā)和地下水定水頭邊界出流。研究區(qū)水量主要輸入途徑是降水和地表水流入，分別占總水量輸入的41.12%和52.71%；主要輸出途徑是蒸散發(fā)和地表水流出，分別占總水量輸出的50.98%和39.72%。研究區(qū)總蓄水量反映水量的盈余和虧缺，經(jīng)計算總蓄水量為正，代表濕地水儲量增加。濕地地表水儲量變化不大，表明濕地海綿層蓄水能力較強(qiáng)；但濕地地表水與地下水交換頻繁，說明濕地海綿層透水性好、下滲能力強(qiáng)，因此在濕地地表水與地下水耦合模擬中考慮濕地海綿層非常必要。另外，研究區(qū)濕地地表水與地下水交互作用顯著且復(fù)雜，在區(qū)域水循環(huán)中發(fā)揮了一定作用。

3.6 濕地地表水與地下水轉(zhuǎn)化關(guān)系

經(jīng)過率定驗證后的濕地地表水與地下水耦合模擬模型可較好地模擬濕地地下水位的動態(tài)變化。因此，可用來研究濕地地表水與地下水的轉(zhuǎn)化關(guān)系，評價流域內(nèi)2005—2019年的濕地地下水與河水水量補(bǔ)排關(guān)系（圖15），圖中正值表示濕地地表水對地下水的補(bǔ)給量，負(fù)值表示地下水向地表水的排泄量（即基流量）。流域內(nèi)濕地地下水與河流的補(bǔ)排關(guān)系以地下水向河流排泄為主，但在每年6—9月份雨季期間，大量降水使得濕地地表水水位升高，造成濕地地表水對地下水的補(bǔ)給。

圖15 研究區(qū)2005—2019年濕地地表水與地下水補(bǔ)排關(guān)系Fig.15 Recharge and discharge between wetland surface water and groundwater in the study area from 2005 to 2019

4 結(jié)論

（1）濕地地表水與地下水耦合模擬需要考慮地形、土壤、植被等要素對濕地水文過程的影響，并對濕地海綿層、濕地明渠流和片流等同時進(jìn)行模擬。

（2）本文所建立的模型充分考慮了濕地水文特征，加入濕地海綿層，并通過二維擴(kuò)散方程模擬濕地片流、一維圣維南方程模擬濕地明渠流，率定和驗證結(jié)果表明模型具有較高精度，能較好地反映研究區(qū)水量的變化狀況。

（3）所構(gòu)建的耦合模型可定量計算濕地地表水與地下水交互作用關(guān)系和濕地水量收支平衡：鏡湖濕地水資源通過地表徑流和降水得到補(bǔ)給，分別占總水量輸入的41.12%和52.71%；以蒸散發(fā)和地表徑流排泄，分別占總水量輸出的50.98%和39.72%。