南嶺山地森林流域退水規(guī)律及影響因素

李澤華，周 平，黃遠(yuǎn)洋，徐 衛(wèi)，譚兆偉

（1. 廣東省科學(xué)院廣州地理研究所廣東南嶺森林生態(tài)系統(tǒng)國(guó)家野外科學(xué)觀測(cè)研究站/廣東省地理空間信息技術(shù)與應(yīng)用公共實(shí)驗(yàn)室/廣東省遙感與地理信息系統(tǒng)應(yīng)用重點(diǎn)實(shí)驗(yàn)室，廣州 510070；2. 重慶地質(zhì)礦產(chǎn)研究院自然資源部重慶典型礦產(chǎn)生態(tài)修復(fù)野外科學(xué)觀測(cè)研究站/重慶市萬盛礦區(qū)生態(tài)環(huán)境保護(hù)修復(fù)野外科學(xué)觀測(cè)研究站，重慶 401120）

對(duì)于包括廣東省的中國(guó)大部分地區(qū)而言，降水的年內(nèi)分配不均導(dǎo)致年徑流中絕大多數(shù)在豐水期以洪水形式入海而難以利用（黃楚珩等，2019）。與豐水期徑流相比，枯水期徑流因?yàn)樵诰S持河流生態(tài)系統(tǒng)穩(wěn)定和保障下游供水安全等方面起到更重要的作用，而受到越來越多的關(guān)注（Bruijnzeel, 2004；曾松青等，2010）。在枯水期內(nèi)由于降水持續(xù)偏少甚至無降水的時(shí)間較長(zhǎng)，徑流表現(xiàn)為逐漸消退的流域退水過程，該過程中徑流成分以地下水出流為主，主要受流域特征和水文地質(zhì)特性影響（穆文相，2020）。Brutsaert 等（1977）通過對(duì)流域概化成具有不透水底板的含水層，根據(jù)Dupuit-Boussinesq 方程的解析解和枯水期徑流消退特征，將流域水文地質(zhì)參數(shù)與退水過程建立聯(lián)系，其中提出的利用多事件的蓄泄關(guān)系推求消退系數(shù)的方法得到廣泛應(yīng)用。隨后有學(xué)者發(fā)現(xiàn)利用多事件的蓄泄關(guān)系識(shí)別存在較大的不確定性，并轉(zhuǎn)而采用基于獨(dú)立事件的蓄泄關(guān)系分析流域退水過程的變化規(guī)律及影響因素（Shaw et al.,2012；高滿等，2013）。然而，在以往基于獨(dú)立事件的蓄泄關(guān)系研究中，針對(duì)退水過程的影響因素，不同學(xué)者得出的結(jié)論存在較大差異。有學(xué)者將主要影響因素歸結(jié)為蒸散發(fā)（Federer,1973;Cadol et al., 2012;Tashie et al., 2020），也有學(xué)者認(rèn)為土壤水分起著關(guān)鍵作用（Shaw et al., 2012;2013），還有學(xué)者提出地形特征和河網(wǎng)結(jié)構(gòu)等主導(dǎo)了退水過程的變化（Biswal et al., 2010, 2014; Patnaik et al.,2015）?？梢?，有關(guān)退水規(guī)律和影響因素的結(jié)論多有不同，對(duì)于地形、植被、土壤和地質(zhì)結(jié)構(gòu)等下墊面條件更趨復(fù)雜的山地森林流域而言，其退水過程在目前更是缺乏完整的機(jī)理解釋和系統(tǒng)的實(shí)證研究，有關(guān)森林是否具有“削峰補(bǔ)枯”作用的爭(zhēng)論也一直在持續(xù)（陳軍鋒等，2001；Wang et al.,2016; Ellison et al., 2017; Li et al., 2021; Zhao et al.,2022）。

南嶺保存有地球同緯度帶上最完整的亞熱帶常綠闊葉林，是中國(guó)南方重要的生態(tài)屏障帶，也是北江和東江等眾多河流的源頭區(qū)，因?yàn)楸环Q為“廣東水塔”，對(duì)粵港澳大灣區(qū)的供水安全起到至關(guān)重要的作用。南嶺山地的地形地貌復(fù)雜性和生態(tài)系統(tǒng)多樣性為代表水文學(xué)發(fā)展方向的生態(tài)水文學(xué)提供了重要的研究對(duì)象（周國(guó)逸等，2018）。根據(jù)廣東南嶺森林生態(tài)系統(tǒng)國(guó)家野外科學(xué)觀測(cè)研究站降水觀測(cè)數(shù)據(jù)顯示，廣東南嶺國(guó)家級(jí)自然保護(hù)區(qū)2021年降水量1 231 mm，與2019年降水量2 592 mm 和2020年降水量2 333 mm 相比減少1 000 mm 以上。在此背景下，本研究基于廣東南嶺森林生態(tài)系統(tǒng)國(guó)家野外科學(xué)觀測(cè)研究站獲取的一手觀測(cè)數(shù)據(jù)，采用基于獨(dú)立事件的蓄泄關(guān)系分析方法，揭示南嶺典型的山地森林流域退水規(guī)律，識(shí)別影響退水過程的主要因素，以期為南嶺山地流域的枯水期徑流模擬及預(yù)測(cè)、水源涵養(yǎng)功能評(píng)估、河流生態(tài)系統(tǒng)保護(hù)與修復(fù)等提供科學(xué)依據(jù)。

1 材料與方法

1.1 研究區(qū)概況

研究區(qū)位于廣東南嶺國(guó)家級(jí)自然保護(hù)區(qū)，屬亞熱帶濕潤(rùn)季風(fēng)氣候，集中分布著中國(guó)最典型的中亞熱帶常綠闊葉林，土壤類型以黃壤土為主。依托廣東南嶺森林生態(tài)系統(tǒng)國(guó)家野外科學(xué)觀測(cè)研究站，位于保護(hù)區(qū)乳陽片區(qū)內(nèi)的實(shí)驗(yàn)流域（經(jīng)緯度范圍為24°55'12"—24°56'06"N 和112°59'46"—113°-0'32"E，集水面積約0.86 km2，海拔范圍為1 115～1 775 m）布設(shè)有測(cè)流堰（圖1），采用6541型浮子水位計(jì)監(jiān)測(cè)流域出口斷面水位變化，在實(shí)驗(yàn)流域內(nèi)采用基于TDR原理的CS655型土壤水分傳感器監(jiān)測(cè)不同深度 （5、10、15、20、40、60、100 cm） 土壤水分體積含量變化，采用CS451型壓力式水位計(jì)監(jiān)測(cè)地下水位變化；采用IRGASON-NM-BB-IC 型開路式渦度相關(guān)系統(tǒng)監(jiān)測(cè)實(shí)際蒸散發(fā)變化，并同步監(jiān)測(cè)降水、氣溫等氣象要素變化；以上數(shù)據(jù)監(jiān)測(cè)頻率均設(shè)定為10～30 min。

圖1 研究區(qū)示意圖Fig.1 Map of the study area

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

利用在測(cè)流堰采集的水位數(shù)據(jù)，通過堰流水力學(xué)公式計(jì)算獲得實(shí)驗(yàn)流域出口河道斷面流量數(shù)據(jù)；對(duì)不同深度的土壤水分進(jìn)行加權(quán)平均，獲得0～100 cm深度層的土壤水分體積含量數(shù)據(jù)；根據(jù)地下水監(jiān)測(cè)井壓力式水位計(jì)的布設(shè)埋深（14.0 m），將監(jiān)測(cè)的地下水位數(shù)據(jù)轉(zhuǎn)換為地下水埋深數(shù)據(jù)；將以上數(shù)據(jù)及同步監(jiān)測(cè)的氣象數(shù)據(jù)按時(shí)間頻率24 h 統(tǒng)一整理，最終獲得逐日平均流量、土壤水分、地下水埋深、潛在蒸散發(fā)、實(shí)際蒸散發(fā)等的時(shí)間序列。南嶺實(shí)驗(yàn)流域監(jiān)測(cè)數(shù)據(jù)基本信息見表1所示。

表1 南嶺實(shí)驗(yàn)流域監(jiān)測(cè)數(shù)據(jù)基本信息Table 1 Description of the data observed at the Nanling experimental catchment

1.3 研究方法

天然流域河道徑流消退過程中，流量Q作為時(shí)間t的函數(shù)，其導(dǎo)數(shù)（即隨t的變化率） 與流域的調(diào)蓄特性有關(guān)：

式中：Q為流量；t為時(shí)間；f為反映流域調(diào)蓄特性的函數(shù)，亦稱流域退水斜率曲線。當(dāng)監(jiān)測(cè)時(shí)間頻率為Δt，連續(xù)觀測(cè)的流量數(shù)據(jù)為Qi和Qi+1，式（1）可近似為：

對(duì)于具有不透水底板的含水層，在忽略源匯項(xiàng)情況下，潛水水流運(yùn)動(dòng)可通過Dupuit-Boussinesq方程描述：

式中：Sy為給水度；K為飽和滲透系數(shù)；x為橫坐標(biāo)，表示在含水層剖面的水平位置；h為地下水深度。Brutsaert等 （1977） 基于對(duì) Dupuit-Boussinesq方程的推導(dǎo)發(fā)現(xiàn)流域退水斜率曲線呈現(xiàn)冪律函數(shù)形式，可表示為：

式中：a和b為消退系數(shù)，分別對(duì)應(yīng)于流域退水斜率曲線在雙對(duì)數(shù)圖中的截距和斜率。其中，截距a主要反映徑流的消退速率，a值越大，消退速率越快；斜率b則是曲線非線性程度的量度，b值增大，曲線凹度越大，意味著高流量時(shí)退水加快而低流量時(shí)退水趨緩。

以均方根誤差RMSE（Root Mean Square Error）作為流域退水斜率曲線的擬合準(zhǔn)則：

式中：（-dQ/dt）i為觀測(cè)的流域退水速率；f（Qi）為基于觀測(cè)的流量Qi和擬合的流域退水斜率曲線計(jì)算的流域退水速率；n為流域退水事件的數(shù)據(jù)序列長(zhǎng)度。基于公式（4）和（5）逐一對(duì)退水事件的流域退水斜率曲線進(jìn)行擬合，分為不對(duì)消退系數(shù)a和b設(shè)置限制（擬合方案1）和對(duì)b設(shè)置限制（擬合方案2）兩種方案。

2 結(jié)果與分析

2.1 退水事件提取

采用2019—2021 年逐日流量時(shí)間序列篩選退水事件，考慮到本研究的流域集水面積較小（0.86 km2），由降水形成的地表徑流可在1～2 d 內(nèi)快速流出，為排除地表徑流對(duì)退水過程分析的影響，退水事件以流量拐點(diǎn)（洪峰）后第3 d 作為退水過程的起始點(diǎn)（Brutsaert et al., 1998），同時(shí)退水過程中滿足流量Q及退水速率-dQ/dt連續(xù)遞減的要求，即排除當(dāng)日流量大于前日流量、當(dāng)日退水速率快于前日退水速率的事件。根據(jù)以上原則共篩選出25 個(gè)退水事件（圖2 和表2），其中2019、2020 和 2021 年分別篩選出 10、6 和 9 個(gè)退水事件，這些退水事件中持續(xù)日數(shù)最短僅4 d，最長(zhǎng)達(dá)15 d，初始流量Q0大小介于0.002～0.092 m3/s。

圖2 2019—2021年退水事件流量過程線Fig 2 Recession events selected during 2019-2021

2.2 退水特征變化

通過在log-log雙對(duì)數(shù)坐標(biāo)空間中繪制-dQ/dt-Q關(guān)系（圖3）發(fā)現(xiàn)，-dQ/dt與Q基本上呈正比例關(guān)系，即較快的退水速率-dQ/dt對(duì)應(yīng)于較大的流量Q，但數(shù)據(jù)點(diǎn)整體尤其在中至尾部（-dQ/dt較小時(shí)）顯得較為散亂。通過不同顏色和形狀區(qū)分不同年份的數(shù)據(jù)點(diǎn)，可以發(fā)現(xiàn)2021年數(shù)據(jù)點(diǎn)與2019和2020年的數(shù)據(jù)點(diǎn)相比明顯向左偏移，由此分別擬合的流域退水斜率曲線斜率b基本保持不變，但截距a將相應(yīng)增大，說明相同流量情況下，2021年的退水速率比2019和2020年更快。

圖3 退水速率-流量關(guān)系（-dQ/dt-Q）Fig.3 Scatter points of the estimated-dQ/dt-Q

擬合方案1不對(duì)消退系數(shù)a和b設(shè)置限制（擬合結(jié)果見表2 和圖4），可以發(fā)現(xiàn)擬合的斜率b值大小介于0.6～3.0，均值為2.0，且大部分集中分布在1.5～2.5區(qū)間；在此基礎(chǔ)上，為便于比較不同退水事件的流域退水斜率曲線，擬合方案2假定各流域退水斜率曲線具有相同的斜率b（本研究假定b=2.0），可以發(fā)現(xiàn)擬合的截距a值大小介于0.6～36.3，均值為14.4（擬合結(jié)果見表2）。

圖4 斜率b的擬合值及概率密度分布Fig.4 Fitted values and their probability density distribution of slope b

表2 退水事件特征及參數(shù)Table 2 Recession events with their characteristics and estimated parameters

通過對(duì)比可以發(fā)現(xiàn)（圖5-a、c），2019年的截距a值在8月前仍然保持在6以下水平，8月后才開始出現(xiàn)增加趨勢(shì)，最大值為14.2；相比而言，2021年的截距a值增加時(shí)間提前至6月，且增幅加大，最大值達(dá)到36.3。相應(yīng)地，2019和2021年的潛在蒸散發(fā)PET 分別為2 046 和2 135 mm，差別不大；而2019和2021年的降水量分別為2 592和1 223 mm，2021年降水量不足2019年降水量的一半，可見-dQ/dt與Q的關(guān)系即流域退水斜率曲線與來水條件關(guān)系密切。

圖5 2019—2021年截距a值及逐日累積降水/潛在蒸散發(fā)變化Fig.5 Temporal variation of intercept a and daily accumulated precipitation/potential evapotranspiration during 2019-2021

2.3 退水影響因素分析

通過比較發(fā)現(xiàn)（圖6），土壤水和地下水對(duì)降水的響應(yīng)總體較好，地下水埋深較土壤水分含量變化略有延遲，說明土壤入滲能力較強(qiáng)，降水能迅速補(bǔ)充土壤水和地下水。7 月、8 月中旬至9 月中旬2個(gè)時(shí)間段連續(xù)無雨日數(shù)增多，同時(shí)潛在蒸散發(fā)PET 也相對(duì)較大，在此作用下土壤水分顯著下降，與之相比由于地下水埋深較深（5.0 m 以上），受蒸散發(fā)影響較小，因此地下水下降幅度相對(duì)緩和。

圖6 2021年逐日降水和潛在蒸散發(fā)（a）土壤水分含量和地下水埋深（b）變化Fig.6 Temporal variation of daily precipitation and potential evapotranspiration(a)soil moisture volumetric content and groundwater depth(b)in 2021

為識(shí)別退水過程的主要影響因素，將擬合方案2下各退水事件擬合的截距a值分別與退水事件的同期潛在蒸散發(fā)PET、實(shí)際蒸散發(fā)ET、地下水埋深及土壤水分含量作相關(guān)分析。從圖7可以發(fā)現(xiàn)，與截距a值相關(guān)程度最高的為地下水埋深（r=0.93*），為正相關(guān)；其次分別為土壤水分含量（r=-0.76*）、實(shí)際蒸散發(fā)ET（r=-0.61）和潛在蒸散發(fā)PET（r=-0.32），均呈負(fù)相關(guān)關(guān)系。其中，截距a值與地下水埋深及土壤水分含量的相關(guān)系數(shù)r通過了99%置信水平下的顯著性檢驗(yàn)?？梢姡叵滤裆钤诹饔蛲怂^程中起最主要的控制作用，地下水埋深越深，退水速率越快。

圖7 截距a值與退水事件同期潛在蒸散發(fā)（a）實(shí)際蒸散發(fā)（b）地下水埋深（c）土壤水分含量（d）關(guān)系散點(diǎn)圖（*99%置信水平下顯著相關(guān)）Fig.7 Scatterplots of intercept a and cocurrent potential evapotranspiration(a)evapotranspiration(b)groundwater depth(c)soil moisture volumetric content(d)(*Significant at p=0.01)

3 討論

有林地包氣帶土層比無林地疏松，土壤入滲能力較強(qiáng)，降水和地表徑流更容易補(bǔ)充土壤水和地下水，從而增加流域的蓄水量，削減洪峰流量和洪水總量，已有研究對(duì)森林“蓄水削峰”作用的認(rèn)識(shí)基本一致（Bruijnzeel, 2004; Ellison et al., 2017），在本研究中也得到了實(shí)驗(yàn)流域水文氣象長(zhǎng)期定位觀測(cè)數(shù)據(jù)的進(jìn)一步證實(shí)。然而，枯水期由于降水減少，河川徑流需要由流域蓄水（主要為地下水）補(bǔ)給時(shí)，對(duì)于森林能否保持徑流均勻和穩(wěn)定，即所謂的“放水補(bǔ)枯”作用，學(xué)者們的認(rèn)識(shí)存在較大分歧（陳軍鋒等，2001；Ellison et al., 2017; Zhao et al.,2022）。Li 等（2021）基于氣候情景的模擬研究發(fā)現(xiàn)，近幾十年來中國(guó)東江流域上游的大規(guī)模植樹造林對(duì)枯水期徑流有較為明顯的促進(jìn)作用。而Zhao等（2022）則發(fā)現(xiàn)位于東南亞季風(fēng)區(qū)的泰國(guó)湄南河流域在經(jīng)歷森林砍伐后，枯水期徑流呈現(xiàn)增加趨勢(shì)。森林覆蓋率的增加或減少對(duì)流域枯水期徑流乃至總徑流造成怎樣的影響在目前仍未有定論。對(duì)于一個(gè)流域而言，其所處的自然氣候帶、地質(zhì)地貌條件及受人類活動(dòng)的干擾等都是影響森林與水的關(guān)系的重要因素。不同的自然氣候帶有著不同的降水、入滲、蒸發(fā)和產(chǎn)匯流規(guī)律。研究區(qū)域和尺度的差異性，以及森林生態(tài)系統(tǒng)本身的復(fù)雜性都會(huì)對(duì)森林水文效應(yīng)的評(píng)價(jià)造成較大的影響，因此，森林水文效應(yīng)的區(qū)域性對(duì)比和跨尺度綜合研究是未來亟需加強(qiáng)的方向（白晉華等，2008；Ellison et al.,2017）。

對(duì)于本研究所依托的實(shí)驗(yàn)流域，流域出口處地下水位長(zhǎng)期低于河床位置，可以判斷該河段為滲失河段，河道水長(zhǎng)期向地下水滲失（地下水補(bǔ)給區(qū)）。與以往研究認(rèn)為影響退水過程的主要因素是蒸散發(fā)和土壤水分等結(jié)論（Cadol et al., 2012; Shaw et al.,2013;Tashie et al.,2020）不同，本研究表明地下水埋深是影響山地森林流域出口退水過程的最主要因素，由于地下水位低于河床位置，流域出口河段處于地下水補(bǔ)給區(qū)，因此地下水埋深越深，河道水向地下水滲失越快，從而導(dǎo)致流域出口退水速率加快。對(duì)于山地森林流域而言，河岸帶森林植被對(duì)水陸生態(tài)系統(tǒng)間的物質(zhì)、能量、信息和生物流動(dòng)具有明顯的邊緣效應(yīng)（白晉華等，2008），由此決定了地表水和地下水之間必然存在緊密而復(fù)雜的交換作用，這種作用對(duì)山地森林流域的枯水期徑流的影響目前仍然知之甚少。在森林的生態(tài)水文效應(yīng)下，山地流域河岸帶的水文地質(zhì)條件、河流-含水層之間的水力梯度、包氣帶厚度及土壤特性等是如何控制河流與含水層系統(tǒng)的水熱交互過程（王平，2018），以及由此對(duì)河流水文情勢(shì)的具體影響機(jī)制亟需進(jìn)一步揭示，從而為山地森林流域的枯水期徑流模擬及預(yù)測(cè)、水源涵養(yǎng)功能評(píng)估、河流生態(tài)系統(tǒng)保護(hù)與修復(fù)等提供科學(xué)依據(jù)。研究區(qū)屬于山地森林流域，集水面積雖然僅有0.86 km2，但海拔跨度比較大（1 115～1 775 m），劇烈的地形起伏決定了下墊面具有顯著的非均質(zhì)性和各向異性，未來需要根據(jù)流域的水文地質(zhì)條件和土壤特性開展更大范圍的長(zhǎng)期定位觀測(cè)，結(jié)合溫度示蹤和水化學(xué)方法（馬瑞等，2013；朱金峰等，2017），量化分析不同河段地表水與地下水的相互作用，從水陸界面的物質(zhì)通量角度進(jìn)一步揭示山地森林流域退水過程的潛在機(jī)理。

4 結(jié)論

采用廣東南嶺森林生態(tài)系統(tǒng)國(guó)家野外科學(xué)觀測(cè)研究站實(shí)驗(yàn)流域2019—2021年逐日平均流量、土壤水分、地下水埋深、潛在蒸散發(fā)、實(shí)際蒸散發(fā)等實(shí)時(shí)監(jiān)測(cè)數(shù)據(jù)，基于由Dupuit-Boussinesq方程推導(dǎo)的流域退水斜率曲線，對(duì)南嶺山地森林流域的退水特征進(jìn)行分析和比較，識(shí)別影響流域退水過程的主要因素，得到以下結(jié)論：

1）退水速率與流量的關(guān)系即流域退水斜率曲線在不同來水條件下變化明顯。流域退水斜率曲線在斜率b保持不變的情況下，截距a在枯水年2021年較大，在豐水年2019 年和平水年2020 年較小，表現(xiàn)出流域退水斜率曲線受豐枯變化影響明顯，相同流量情況下退水速率在枯水年要比在豐水年更快。

2）流域土壤入滲能力較強(qiáng)，降水能迅速補(bǔ)給土壤水和地下水。0～100 cm土壤水分體積含量和地下水埋深對(duì)降水過程有較好響應(yīng)，相比而言，地下水較土壤水變化略有延遲，說明土壤入滲能力較強(qiáng)，降水能迅速入滲補(bǔ)給土壤水和地下水；相對(duì)土壤水而言，地下水由于埋深較深（5.0 m以上），受蒸散發(fā)影響較小。

3）地下水埋深是影響流域退水斜率曲線變化的最主要因素。退水事件擬合的截距a值與事件同期潛在蒸散發(fā)PET、實(shí)際蒸散發(fā)ET、土壤水分和地下水埋深的相關(guān)分析結(jié)果表明，地下水埋深是影響流域退水斜率曲線變化的最主要因素，兩者呈顯著正相關(guān)關(guān)系（r=0.95），說明地下水埋深越深，流域退水速率越快。