基于地表水-地下水耦合模型的海底地下淡水排放量估算

羅 璐 , 俞 烜 , 劉丙軍 , 陳曉宏

(1. 中山大學土木工程學院, 廣東 珠海 519082; 2. 南方海洋科學與工程廣東省實驗室(珠海), 廣東 珠海519082)

全球有大約40%的人口和2/3的大型城市集中在沿海地區(qū)[1], 日益加劇的人類活動使得沿海地區(qū)所承受的環(huán)境壓力越來越大, 近年來近海環(huán)境水污染事件如赤潮、有害藻華等也頻繁發(fā)生, 這與陸源營養(yǎng)物質(zhì)和污染物的輸入密切相關[2-3]。以往的研究更多的是關注通過地表徑流和大氣沉降的污染物質(zhì)輸入, 而忽略了通過海底地下水排泄(SGD)輸入的營養(yǎng)物質(zhì)與污染物對海洋生態(tài)環(huán)境的影響。海底地下水排泄(SGD)是指不考慮流體成分、來源或驅(qū)動力,通過陸架邊緣由海底進入海洋的所有流體, 主要包括陸源地下淡水排泄(submarine fresh groundwater discharge, SFGD)和再循環(huán)海水地下水排泄(submarine circulate groundwater discharge, SCGD)[4-5]。已有研究表明, SGD所攜帶的營養(yǎng)鹽、有機物及其他化學物質(zhì)在某些地區(qū)接近甚至超過了地表徑流的輸入量[6-7],SFGD作為SGD的重要組成部分, 帶來大量的陸源營養(yǎng)物質(zhì), 是陸地向海洋輸送化學物質(zhì)重要且隱蔽的通道, 與沿海生態(tài)環(huán)境密切相關[8]。同時, SFGD對于某些沿海社區(qū)來說也是重要的淡水來源[5], 研究SFGD對近岸人類生活環(huán)境和生態(tài)環(huán)境具有重要影響和意義。

目前, 研究SFGD的方法主要有三種: 直接測量法、同位素示蹤法和水平衡與水文模型法[9]。SFGD速率具有空間分布不均勻性, 直接測量法難以在空間和時間尺度上對SFGD速率進行展延[10], 無法進行流域或區(qū)域尺度上的估算。同位素示蹤法[11-13]適用于不同尺度范圍SGD的評估, 但并不能很好的解釋與其總SGD相關的陸源地下淡水成分[14], 無法對研究區(qū)域進行SFGD空間分布以及時序變化分析。水平衡與水文模型法是用于估算區(qū)域SFGD速率的常用方法[15-16], 許多研究采用水文模型法對海岸SFGD進行估算時[17-18], 由于未考慮地下水與河流的交互作用, 無法反映SFGD對河流季節(jié)變化的響應規(guī)律, 由于海岸水文地質(zhì)參數(shù)獲取較難, 將水文耦合模型擴展應用到海岸帶研究的較少[19]。隨著海岸地區(qū)經(jīng)濟文化重要性的不斷加強, 近海生態(tài)環(huán)境要求的不斷提高, 進一步了解沿海流域與河口之間的水文過程, 對沿海地區(qū)陸源地下淡水排泄進行高分辨率時空估計顯得愈加重要[9,20]。

本次研究利用賓州州立集成水文模型(Penn State Integrated Hydrologic Model)PIHM將地表水-地下水進行耦合, 構建海岸地表水-地下水耦合模型, 并應用于廣東省東南部紅海灣區(qū)域, 估計沿海流域的陸源地下淡水通量, 對SFGD速率進行時空變異分析,從而進一步了解沿海流域與河口之間相互作用的水文過程。研究基于賓州州立集成水文模型PIHM和高分辨率水文數(shù)據(jù)集, 通過模擬沿海地表、地下和海洋之間的水量交換過程, 估計流域SFGD速率, 對其空間上進行流域尺度分析、影響因素分析, 時間上進行降水響應分析, 揭示流域SFGD速率的時空變異規(guī)律。

1 研究區(qū)概況與數(shù)據(jù)來源

1.1 研究區(qū)概況

研究區(qū)域位于廣東省東南部沿海, 如圖1所示,其地理位置為22°42′21″～23°29′25″N, 114°57′47″～115°46′48″E, 北倚蓮花山脈, 南瀕南海, 從東至西依次由螺河流域片區(qū)、黃江河流域片區(qū)以及赤石河流域片區(qū)構成。區(qū)域總面積3 510 km2, 海岸線總長190 km,地勢西北高東南低, 研究區(qū)內(nèi)所有河流匯入南海。

圖1 研究區(qū)域概況Fig. 1 Study area

研究區(qū)屬南亞熱帶季風氣候區(qū), 年平均降水量為1 900～2 500 mm。降雨時空分配不均勻, 空間上從西北山區(qū)向東南沿海區(qū)遞減, 時間上夏多冬少,其中4月—9月(汛期)的降水占全年降水的85%以上。年平均氣溫為21.9 ℃[21], 全年無霜期為347 d, 年平均日照為2 179.1 h, 日照率為49%; 太陽輻射總量年平均120 kcal/cm2以上。年平均水面蒸發(fā)量為1 351.5 mm,沿海一帶水面蒸發(fā)量較大, 山丘區(qū)水面蒸發(fā)量較小。研究區(qū)濕度為78%～83%。據(jù)陸豐氣象站統(tǒng)計資料表明,年最大瞬時風速為40 m/s (風向為SE、NW), 年最大10 min平均風速為25 m/s (風向為W)。

1.2 數(shù)據(jù)來源

模型所需的基礎數(shù)據(jù)包括研究區(qū)地形數(shù)據(jù)(DEM)、土地利用類型及分布信息(LUCC)、土壤屬性數(shù)據(jù)庫信息、水文氣象資料等。

DEM數(shù)據(jù)選用30 m空間分辨率的美國國家航空航天局(NASA)的SRTM數(shù)據(jù)(http://srtm.csi.cgiar.org)。水文氣象數(shù)據(jù)來自中國氣象局, 主要包括日降水、日平均氣溫、日照、相對濕度以及風速等。模型徑流模擬結果率定與驗證的標準數(shù)據(jù)來自區(qū)域內(nèi)唯一水文站蕉坑站2001年—2002年實測日徑流量,數(shù)據(jù)系列完整連續(xù)。

本研究的土壤類型數(shù)據(jù)來源于中國科學院資源環(huán)境科學數(shù)據(jù)中心(http://www.resdc.cn)。由Pedotransfer函數(shù)[23]來估算飽和水力傳導系數(shù)以及Van Genuchten持水參數(shù)α和n(表1)。

表1 研究區(qū)域土壤類型及主要參數(shù)Tab. 1 Descriptions and essential parameters of the soil types in the study area

土地利用類型選用清華大學地球科學系統(tǒng)研究中心2017年全球土地覆蓋10 m空間分辨率數(shù)據(jù)集[24](http://data.ess.tsinghua.edu.cn/fromglc10_2017v01.html)。將土地利用類型轉換為PIHM模型自帶的土壤屬性數(shù)據(jù)庫,研究區(qū)內(nèi)有8種土地覆蓋類型, 包括常綠林(46.7%)、耕地(32.0%)、城市地面(6.4%)、開放水域(5.9%)、草地(5.1%)、灌木林(3.2%)、荒地(0.6%)和林木濕地(0.2%)。

2 研究方法

2.1 模型介紹

賓州州立集成水文模型PIHM是基于物理過程的、完全耦合的分布式水文模型[25], 以完全耦合的方式模擬截留、入滲、補給、蒸發(fā)蒸騰、地表徑流、地下徑流和河網(wǎng)匯流過程(圖2)。區(qū)域分解受觀測點、邊界條件等水文特征的約束, 根據(jù)流域的地貌或水文特征, 改變流域空間離散的分辨率。每個控制體上進行水文過程建模, 包括地表徑流、地下徑流和渠道徑流的偏微分方程, 以及葉面截留、入滲、補給和蒸散發(fā)的常微分方程。采用有限體積法將偏微分方程離散化為常微分方程, 每個網(wǎng)格都有對應的局部常微分方程系統(tǒng), 整個研究區(qū)域組合形成一個總常微分方程系統(tǒng), 使用SUNDIALS進行求解[26]。

圖2 PIHM模型產(chǎn)匯流結構Fig. 2 PennState Integrated Hydrologic Model processes

在PIHM的原始方程中, 側重于對閉合流域進行水文模擬。通過對方程進行修改, 在流域的海岸邊界處添加水頭壓力[27], 將海岸水文邊界條件整合到模型中, 使得PIHM能夠考慮海岸邊界對海岸流域的影響。

2.2 基本方程

總蒸發(fā)量為林冠截留量(ec)、植被蒸騰量(et)、土壤蒸發(fā)量(es)的總和。采用Penman Monteith方程計算潛在蒸發(fā)量, 采用Noah_LSM方程來計算實際蒸散發(fā)量:

式中,ep為潛在蒸散發(fā)量;Rn為植被地表凈輻射;G為土壤熱通量密度;εs- εa為空氣蒸汽壓差額;ρa為空氣密度;Cp為空氣比熱; Δ為飽和蒸氣壓-溫度關系曲線的斜率;γ為計量常數(shù);rs,ra分別為表面和空氣動力阻力;σf為植被覆蓋率;Wc為截冠含水量;S為最大冠層容量;Bc為冠層阻力;β為土壤水飽和度;θref為田間持水量;θw為萎蔫點;hv為植被截留量;pv為總降水量當量;pt為地表儲水的入滲量和徑流量; 下標m為空間網(wǎng)格編號, 范圍從1到三角形網(wǎng)格總數(shù)。

趙新宇等(2013)利用吉林大學公眾幸福指數(shù)課題組關于2012年中國公眾主觀幸福感問卷調(diào)查數(shù)據(jù)，運用有序概率模型考察了絕對收入、相對收入和預期對公眾主觀幸福感的影響。發(fā)現(xiàn)我國存在“幸福悖論”現(xiàn)象，相對收入對公眾主觀幸福感有顯著促進作用，其效果強于絕對收入；預期對于中、低收入群體的主觀幸福感具有顯著正向作用。

使用二維Saint Venant方程來描述坡面漫流, 使用一維Saint Venant方程描述渠道徑流, 使用Richard方程來描述地下徑流, 用半離散方式表示為:

式中,h0為地表淺水深度;為元素i到相鄰元素j的歸一化側向地表流量;pt、q+、e分別為凈降水量、入滲量、蒸發(fā)量; 下標m為空間網(wǎng)格編號, 范圍從1到三角形網(wǎng)格總數(shù);θs為含水率;hu為非飽和地下水深度;hg為地下水深度;q0為非飽和-飽和區(qū)之間的地下水通量;為從元素i到鄰近元素j的歸一化側向地下流量;hc為河道水深;p、e分別為來自河道段的降水量和蒸發(fā)量;分別為來自河道兩側的地表側向流量和與河道相互作用的地下水流量; 各河道段的上下游流量分別為和; 下標k為河道段編號,范圍從1到河道段總數(shù)。

淡水地下水通量采用以下公式進行計算[27]:

式中,qfreshSGD為淡水地下水通量, 單位m2/s;Ssat為飽和地下水儲量, 單位m;Zb為基底高程, 單位m;K為水力傳導系數(shù), 單位m/d;SL為潮位高, 單位m, 設為平均海平面高度。

2.3 模型構建

用三角網(wǎng)格剖分法對模型進行剖分, 共得到節(jié)點6 473個, 三角形單元12 206個, 河網(wǎng)剖分為231條河道段。沿海岸線單元設置為定水頭邊界。單元平均面積約為2.88×105m2, 單元邊長為500～1 000 m, 主要徑流寬度為50～100 m。模型空間離散圖如圖1(e)所示。模擬區(qū)域初始地表水設置為0 m,初始地下水位設定為地表以下5 m, 沿海邊界處的初始地下水位高程設置為0 m(平均海平面)。通過輸入1992年至2014年的氣象數(shù)據(jù)對研究區(qū)域的水文過程進行模擬。

3 結果與討論

3.1 模型率定與驗證

以1992年為預熱期, 利用區(qū)域內(nèi)唯一水文站蕉坑站2001—2002年實測日徑流量, 對模型土壤及土地覆蓋關鍵參數(shù)進行率定和驗證, 選取相關系數(shù)(R)、相對誤差(RE)、納什效率系數(shù)(NSE)[28]和Kling-Gupta系數(shù)(KGE)[29]作為評價指標。R、NSE用于評估日徑流過程模擬效果, RE用于分析年徑流總量的模擬精度, KGE是一個集成相關系數(shù)、模擬值和實測值的均值與誤差的評價指標, 用于綜合反映模擬效果。關鍵參數(shù)包括土壤孔隙率、土壤水力參數(shù)、含水層導水率、優(yōu)先流深度、河床曼寧糙率。早期模型應用提供了這些關鍵參數(shù)上下界的參考值[30-31]以及率定方法[26]。由于地表水-地下水耦合模擬的計算過程耗時高, 需要將參數(shù)分成兩組: 產(chǎn)流參數(shù)組與蒸發(fā)參數(shù)組, 進行分別調(diào)參。產(chǎn)流參數(shù)調(diào)參選取了2001年7月1日到至7月15日的降雨徑流過程, 利用CMA-ES(協(xié)方差矩陣自適應進化策略)方法進行優(yōu)化。將優(yōu)化結果代入后, 利用2001年水量平衡進行蒸發(fā)組參數(shù)的率定。驗證期為2002年1月1日至2002年12月31日。模型率定和驗證期結果見圖3。

圖3 蕉坑水文站實測日徑流量與模擬日徑流量過程Fig. 3 Observed and simulated streamflows at Jiaokeng station

圖3中可以看出, 率定期和驗證期模擬的日均流量過程與實測過程吻合。NSE在率定期和驗證期分別為0.74和0.79,R分別為0.90和0.91, 年徑流相對誤差RE低于10%, KGE分別為0.78和0.73, 四個評價指標充分說明模擬滿足精度要求。模型率定得到的主要參數(shù)結果見表1和表2。

表2 研究區(qū)域土地覆蓋類型及主要參數(shù)水Tab. 2 Land cover types and essential parameters in the study area

3.2 SFGD影響因素分析

在本次研究中, SFGD的影響因素主要為地形坡度(i)、滲透系數(shù)(k)以及含水層厚度(t)。模型離散化將研究區(qū)域沿海邊界剖分成216段, 計算每一段邊界的坡度、滲透系數(shù)、含水層厚度以及年均單位長度SFGD速率, 通過SFGD速率與各參數(shù)關系圖(圖4)可以看出,當i≤2%時, lg(SFGD)隨i的增加而增加, 當i＞2%時,lg(SFGD)逐漸趨于穩(wěn)定值。按照地形坡度是否大于2%,對SFGD速率與地形坡度(i)、滲透系數(shù)(k)以及含水層厚度(t)進行多元線性回歸分析(表3)。結果顯示, 當i≤2%時, lg(SFGD)與地形坡度呈極顯著正相關性(P＜0.01,ri= 0.617), 與滲透系數(shù)(k)無顯著相關性(P＞0.05,rk= -0.101), 與含水層厚度(t)關系不顯著(P＞0.05,rt= 0.465)。表明地形坡度是SFGD速率的主控因子, 地形坡度通過影響水力坡度從而影響流向海岸的地下水,隨著地形坡度的增加, SFGD速率增加。SFGD速率空間分布圖(圖5)顯示, 地形坡度大的沿海地區(qū)SFGD速率大, 地形平緩或低洼地區(qū)SFGD速率較小。當i＞2%時, 回歸模型相關系數(shù)R較小, 表明回歸模型沒有明顯相關性, lg(SFGD)與各個自變量均無顯著相關。

圖4 模擬SFGD速率與相關水文參數(shù)關系圖Fig. 4 Relationship between submarine fresh groundwater discharge (SFGD) flux and related parameters

圖5 SFGD速率空間分布圖Fig. 5 Spatial distribution maps of SFGD

表3 lg(SFGD)多元線性回歸分析系數(shù)Tab. 3 Coefficients of multiple linear regression analysis of lg(SFGD)

3.3 SFGD流域尺度分析

研究區(qū)總面積為3.51×103km2, 海岸線長190.32 km,從東至西按照主河道依次可劃分為三個子流域片區(qū),分別是螺河流域片區(qū)、黃江河流域片區(qū)和赤石河流域片區(qū)。通過分析研究區(qū)1993年—2014年共22 a的SFGD模擬結果(表4), 年平均SFGD速率為20.69×107m3/a (2.98 m2/d或0.016 1 cm/d), 占研究區(qū)內(nèi)河道徑流入海總量的3.28%。三個子流域中, 赤石河流域的單位長度SFGD速率最大, 為3.21 m2/d, 黃江河流域為3.05 m2/d, 與赤石河流域相近, 螺河流域的單位長度SFGD速率最小, 為0.7 m2/d, 根據(jù)4.2及4.3的分析結果, 造成這種差異的重要原因是由于三個子流域之間地形和形狀的差異, 螺河流域海岸線短, 且沿海地區(qū)幾乎都是平原, 地表坡度小, 因此海底地下水排泄量小。

表4 海底地下水排泄量結果統(tǒng)計Tab. 4 Results of submarine groundwater discharge

3.4 SFGD對降水的響應分析

統(tǒng)計三個子流域片區(qū)入?？谔?993年—2014年間主要降水事件中地表徑流峰值與SFGD通量峰值響應時間, 采用配對t檢驗進行分析比較(圖6), 螺河流域片區(qū)、黃江河流域片區(qū)及赤石河流域片區(qū)檢驗的t值分別為7.638、15.426、6.603, 其雙側檢驗的P值均小于0.05, 表明檢驗結果差異顯著。對比分析發(fā)現(xiàn)相對于地表徑流, SFGD對降雨的響應更快, 平均快1～3 d。這種差異的產(chǎn)生可能是由于SFGD與流域地表徑流的補給區(qū)不同造成的。Sawyer等[32]在對美國海岸的陸源地下淡水進行估算時, 將沿海流域依據(jù)其地貌地質(zhì)劃分為了不同類別的區(qū)域(圖7),SFGD的補給區(qū)為徑流集水區(qū)之外的向海集水區(qū)域(圖7中的白色區(qū)域), 該區(qū)域內(nèi)的地下水直接流向海洋。而地表徑流則來源于徑流集水區(qū)(圖7中的灰色區(qū)域), 在本研究中, 計算的地表徑流降水響應時間為流域入海口處的降水響應時間, 徑流集水區(qū)的流域范圍要更大, 因此其匯流時間會大于SFGD補給區(qū)的匯流時間。以上分析表明, 當發(fā)生降雨事件時,陸源地下淡水會先于地表徑流達到峰值, 因此, 在沿海地區(qū)進行地下水資源管理和SFGD觀測時, 需要加強監(jiān)管以更快進行應對。以上分析表明, 在沿海地區(qū)進行地下水資源管理和SFGD觀測時, 需要投入更多的精力以更快進行應對。

圖6 地表徑流峰值與SFGD峰值響應時間對比Fig. 6 Comparison between the peak values of runoff and SFGD flux

圖7 SFGD補給區(qū)示意圖(改自Sawyer等[32])Fig. 7 Recharge area of SFGD (modify from Sawyer et al.[32])

3.5 SFGD模擬結果比較

研究區(qū)域位于廣東省東南部紅海灣區(qū)域, SFGD速率為2.1×108m3/a, 該區(qū)域目前尚無SFGD評估結果可供直接比對。最近兩篇估算全球海岸SFGD通量的研究中, Zhou等[18]關于研究區(qū)域的估算結果為2.4×108m3/a, 與我們的估算結果相近; Luijendijk等[17]關于研究區(qū)域的估算結果為 2.66×107(2.16×106～1.48×108) m3/a, 低于我們的評估結果。在進行全球海岸SFGD評估中, 采用的流域數(shù)據(jù)庫比較粗糙, 而且由于考慮到模型計算成本, 在局部區(qū)域上會進行簡化計算。本研究采用高分辨率水文氣象數(shù)據(jù)集和詳細的水文地質(zhì)參數(shù), 考慮了地表水-地下水交互作用,對SFGD的時空分析有更高的分辨率。

將本研究的SFGD估計值與表5中其他研究區(qū)域的SFGD估計值進行比較。為方便比較, 將單位統(tǒng)一為單位時間內(nèi)單位海岸線長度排放量(m2/d), 即研究區(qū)域海岸線長度除SFGD排泄量。研究區(qū)的SFGD速率為2.78 m2/d, 在其他海岸SFGD評估結果范圍內(nèi)。沿海流域受到復雜的近海動力過程影響, 如潮汐、波浪以及海平面變化等。Kuan等[33]結合物理實驗和數(shù)值模擬法研究改變內(nèi)陸淡水輸入以及潮汐力對海底地下水排泄的影響, 發(fā)現(xiàn)潮汐會通過改變淡水排泄區(qū)域?qū)挾葟亩餝FGD時間排放規(guī)律的變化; Michael等[34]的研究表明海平面上升對SFGD的影響主要取決于沿海含水層的水文地質(zhì)屬性, 通量控制屬性(flux-controlled)的沿海含水層受到海平面上升帶來的影響較小, 而水頭控制屬性(head-controlled)的沿海含水層則會受到較大影響。本文對研究區(qū)域近海驅(qū)動力進行簡化處理, 將海岸邊界設置為定水頭邊界, 對小時尺度的SFGD估算存在一定誤差, 在日尺度及月尺度的誤差影響更為明顯, 將在以后的研究中進行進一步的考慮。

表5 其他地區(qū)SFGD速率Tab. 5 Submarine fresh groundwater discharge in different coastal zones

4 結論

本文將賓州州立集成水文模型PIHM拓展應用到沿海地區(qū), 通過開發(fā)沿海區(qū)域尺度地表水—地下水耦合模型, 估算區(qū)域SFGD通量, 深入揭示SFGD通量時空變異特征, 進一步理解沿海流域與河口之間相互作用的水文過程。模型在廣東省東南部紅海灣區(qū)域模擬效果良好。主要結論如下:

1) 當坡度i≤2%時, 地形坡度是SFGD的主控因子, 隨著地形坡度的上升, SFGD排放速率增加,地形坡度大的沿海地區(qū)SFGD速率大, 地形平緩或低洼地區(qū)SFGD速率較小; 當坡度i＞2%時, SFGD速率與各個自變量均無顯著相關。

2) 研究區(qū)域年均SFGD通量占研究區(qū)內(nèi)河道徑流年均總量的3.28%。三個子流域片區(qū)之間SFGD通量的差異主要是流域形狀和地形的差異引起的, 螺河流域海岸線短且沿海主要為平原地區(qū), 因此海底地下水排泄量小。

3) 分別對三個子流域片區(qū)的地表徑流峰值與SFGD峰值在主要降雨事件中的響應時間進行分析,結果表明SFGD比地表徑流對降雨的響應更快, 峰值出現(xiàn)時間比地表徑流平均早1～3 d。

4) 本文將地表水-地下水耦合模型應用于海岸地區(qū), 所建立的地表水-地下水耦合模型為分析SFGD通量的時空變異規(guī)律提供了一種有效的數(shù)值方法?，F(xiàn)階段模型缺少與流域內(nèi)地下水水位數(shù)據(jù)的對比, 將在以后的工作中著重考慮。