?

閩江入海斷面溶解無機氮長時間序列分析及入海通量估算

2024-01-05 13:14陳克亮吳燁飛陳文花林云杉任保衛
中國環境監測 2023年6期
關鍵詞:閩江入海徑流量

王 顥,陳克亮,吳燁飛,陳文花,林云杉,白 亮,任保衛

1.福建省環境監測中心站,福建 福州 350003 2.自然資源部第三海洋研究所,福建 廈門 361005 3.廈門海洋職業技術學院,福建 廈門 361100 4.自然資源部海洋生態保護與修復重點實驗室,福建 廈門 361005 5.福建省漁業資源監測中心,福建 福州 350003

溶解無機氮(DIN)包括亞硝酸鹽氮(NO2-N)、硝酸鹽氮(NO3-N)和氨氮(NH3-N),是海洋初級生產者維持生長和繁殖所必需的營養物質。過量的DIN可導致富營養化,并在一定條件下引發赤潮等海洋生態環境問題[1]。河流輸送是DIN進入海洋的主要途徑之一,河流入海DIN濃度、組分、通量及其同其他生源要素之間關系的改變可導致河口及鄰近海域浮游植物種群結構發生變化,進而對近海生態系統造成結構性影響[2-4]。閩江是福建省流域面積最大的河流,干流全長為210.7 km,在下游于閩江口注入東海。閩江口是中國重要河口之一,近年來其生態系統始終處于亞健康狀態,富營養化程度較為嚴重[5]。由于閩江口DIN分布不僅受到閩江徑流輸入影響,還受到臺灣海峽流和浙閩沿岸流的季節性作用[6],故對閩江DIN入海濃度、組分及通量變化進行分析對于深入研究閩江口富營養化機制十分必要。

目前,國內對于入海通量的研究主要通過模型和實測2種方法進行。模型法利用產、排污模式結合GIS對通量進行估算[7-8],其優點是易于工況比較,但計算過程需使用大量經驗系數,在未經驗證的情況下偏差較大。實測法基于現場實地監測結果,早期人工成本較高。近年來隨著在線自動監測技術的發展,實測法成本降低的同時精細度不斷提高,逐漸成為入海通量監測的發展趨勢。許多學者利用實測結果對包括閩江在內的中國江河污染物入海通量進行了研究[9-11],但同水文、氣象領域相比,環境類參數監測起步較晚,積累數據相對有限,導致污染物入海通量的研究報道多集中于5~10年甚至更短時間段,10年以上長時間序列分析報道尚不多見。此外,入海通量的估算需要用到大量水文數據,但部分區域受資料獲取限制,在估算過程中多使用年均濃度或年均徑流量,導致結果過于粗略。由于改革開放以來中國沿海地區社會經濟發展非常迅速,污染排放出現較大變化,短時間序列分析對于污染排放與環境響應之間的關系判斷不夠準確。因此,該研究利用1985—2021年閩江水質監測資料對DIN入海濃度的長時間序列變化進行分析,并結合逐月徑流量資料對DIN入海通量進行估算,從而獲得相對精細的結果。在此基礎上,結合統計年鑒等資料對37年來影響閩江DIN入海濃度及通量變化的主要因素進行分析和討論,為進一步深入研究閩江口富營養化演變及成因機制,推進污染物入??偭靠刂萍傲饔?海域協同治理提供技術參考。

1 數據來源與研究方法

1.1 研究區域介紹

閩江位于福建省中北部(圖1),流域面積為60 992 km2,多年平均降水量為1 717.3 mm(福建部分),水資源十分豐富。閩江梯級電站開發程度較高,位于干流的水口電站是華東地區最大的水電站。電站于1986年開工,1993年4月完成蓄水,1996年正式投入運行,屬不完全季調節電站[12]。電站下游58 km處為福建省省會福州市,閩江在此穿城而過后折向東北進入河口區。位于馬尾區的閩安至長樂區金剛腿為傳統河海分界線,環保部門在此設有水質監測斷面。閩江感潮河段長度約為72.5 km,覆蓋整個福州市區。潮區界上游的竹岐水文站為入??刂普?其觀測的流量通常作為閩江干流入?;鶞柿髁?。

1.2 數據來源

1985—2021年閩江入海斷面(閩安)水質監測數據來源于福建省地表水環境監測網。其中1985—2002年的監測頻率較低,為豐水期(5月)、平水期(9月)、枯水期(11月)各1次,個別月份有加測,全年3~4次。2003—2008年監測頻率提高為逢單月監測1次,全年共6次。2008年以后又提高至每月1次,全年共12次。采樣選擇在小潮低平潮以確保水樣鹽度在2以下,并在8 h內完成實驗室分析。其中NO2-N采用磺胺-萘乙二胺比色法,NO3-N采用酚二磺酸比色法,NH3-N采用納氏試劑比色法,上述分析方法在37年間未發生過變化,3種氮鹽濃度的代數和即為DIN濃度。徑流量觀測數據(竹岐水文站)的逐月值來自水利部門發布的水情年報等資料[13-15]。用于關聯分析的相關數據來自歷年福建省及閩江流域所屬設區市(縣)發布的統計年鑒。

1.3 研究方法

季節性肯達爾檢驗(SK檢驗)由美國地質調查局(USGS) HIRSCH等[16]提出,廣泛用于環境科學領域趨勢性分析與評估[17-19]。SK檢驗的優勢在于僅要求年度水質數據遵從相同的概率分布形式,不要求數據一定為正態分布,部分時段數據缺失或未檢出對結果也不構成影響。由于水質數據變化容易受到徑流量變化的影響,個別洪旱事件可導致水質結果出現極值,從而對SK檢驗結果造成干擾。為剔除這一干擾,采用CLEVELAND等[20-21]提出的局部加權回歸散點平滑法(LOWESS),以徑流量作為輔助變量對水質數據進行處理后獲得殘差(也稱流量調節濃度),進一步對殘差進行SK檢驗獲得變化趨勢。同時,為考察徑流量的影響程度,運算時一并給出不經LOWESS處理的SK檢驗結果供對照分析。此外,由于數據時間跨度較大,受水文水質等條件影響,DIN及其各組分的中短期趨勢和長期趨勢可能存在差異,故在討論趨勢時采用整體分析和分段分析相結合的方式。其中整體分析使用37年的數據直接進行SK檢驗,結果反映的是長期變化趨勢。分段分析則結合閩江中下游發生的重大事件將其分為4個時段,結果反映的是中短期變化趨勢。4個時段具體為1985—1992年,代表水口電站大壩蓄水前的時段;1993—1999年,代表水口電站大壩蓄水后的時段;2000—2011年,代表福州市城市布局“東擴南進”,工業外遷的時段;2012—2021年,代表后工業時代城市建設加速,人口大量涌入的時段。SK檢驗使用USGS提供的Kendall計算程序在DOS界面下進行,具體運算方法參見USGS 2005-5275號報告[22]。

入海通量的計算公式為

(1)

式中:Qi為第i個時段的平均徑流量,Ci為第i個時段的平均濃度。當水文和水質觀測不同頻時,采用代表濃度Cri代替平均濃度Ci。

(2)

通常情況下,Qi觀測頻率較高而Cri觀測頻率較低。由于觀測頻率越低,通量估算產生的偏差越大,故有部分學者試圖采用模型內插法使Qi和Cri同頻以便于提高通量計算精度[23]。但對于有梯級水電開發的河流,其流量-濃度關系較為復雜,模型選擇不當時可能會獲得相反效果,故該研究不對Cri做差分處理,即直接進行通量計算。由于1985—2002年監測頻率只有3~4次/年,故筆者在計算通量時根據閩江流域的水文特征,按照豐水期(4—8月)、平水期(9—10月)、枯水期(11月至次年3月)進行劃分,先計算各水期的入海通量,再將3個水期的通量進行加和得到全年的通量。為評估較低監測頻率對通量估算偏差的影響程度,選取2008年以后各年數據按公式(3)計算相對偏差。

(3)

式中:各項含義與公式(1)和公式(2)相同,結果顯示,2種代表時段的估算結果除個別洪澇年份相對偏差為25%左右外,大部分年份中相對偏差為10%以內,整體偏差比較小。

2 結果分析

2.1 入海徑流量年際變化情況

竹岐水文站的觀測結果見圖2。

圖2 1985—2021年閩江入海年徑流量變化Fig.2 Annual runoff variation of Minjiang River during 1985-2021

由圖2可知,閩江入海徑流量年際變化比較劇烈。按照1950—2020年多年平均徑流量(5.40×1010m3)為基準進行統計,1985—2021年閩江入海徑流量距平百分率超過10%的有33年,其中負距平有19年,正距平有14年。距平百分率超過20%的有20年,其中負距平為12年,正距平為8年。整體上負距平年數明顯多于正距平年數,與文獻[24]中閩江徑流量在近年總體呈下降趨勢的結論一致,但距平百分率超過50%的均為正距平年,分別為1998年(59%)、2010年(57%)和2016年(75%)。此外,同2000年前相比,2000年以后出現較大距平的年份較多,表明進入21世紀后,閩江入海徑流量波動明顯增強。

2.2 入海斷面DIN及其各組分濃度年際變化

入海斷面DIN及其各組分濃度數據、年際變化見表1和圖3,其中NO2-N濃度整體呈先升后降的特征。1985—2008年NO2-N濃度整體有所上升但波動較為劇烈,年際間差別較大。2009年以后,NO2-N濃度逐漸回落并穩定于0.02 mg/L附近。NO3-N濃度整體呈“升-降-升”的特征,其中1985—1991年NO3-N濃度不斷升高并達到最大值(2.722 mg/L),1991年以后NO3-N濃度逐漸回落,至1999年時已接近1985年的水平。2000年以后,NO3-N濃度呈階梯性上升(上升-回落-波動-再次上升),其中,2002—2005年及2019—2021年是NO3-N上升明顯的2個時段。NH3-N濃度整體呈先升后降的特征。其中1985—1992年NH3-N濃度不斷升高并達到最大值(0.643 mg/L)。1992年以后,NH3-N濃度經歷了較長穩定期后開始波動下降,其中2008—2015年和2016—2020年下降較為明顯,并在2021年出現了監測以來的最小值。DIN濃度年際變化主要受NO3-N濃度和NH3-N濃度的共同影響。2000年以前,NO3-N和NH3-N濃度的變化趨勢較為一致,故DIN濃度也表現出先升后降的特征。2000年以后,NO3-N濃度同NH3-N濃度的變化趨勢發生背離,因NO3-N濃度較大且變化相對劇烈,故DIN濃度變化趨勢同NO3-N濃度較為接近,大部分時段呈小幅波動上升狀態。

表1 1985—2021年閩安斷面DIN及其各組分濃度特征Table 1 Concentration characteristics of DIN and its components in Min’an section during 1985-2021

圖3 閩安斷面DIN及其各組分的年際變化Fig.3 Interannual variations of DIN and its components in Min’an section

2.3 SK檢驗結果

SK檢驗結果見表2。由表2可知,NO2-N濃度在1985—2021年整體呈顯著下降趨勢,但在劃分的4個時段中Sen斜率依次由正轉負,表明分時段NO2-N濃度存在先升后降的趨勢。Sen絕對值極小且P值始終大于0.1,表明在波動狀態下的NO2-N濃度變化較為緩慢,各分時段趨勢均未達到顯著水平。此外,有LOWESS的檢驗結果與無LOWESS的檢驗結果一致,表明徑流變化對NO2-N濃度變化影響不顯著。

表2 閩安斷面DIN及其各組分時間變化趨勢分析結果Table 2 Temporal variation trend analysis results of DIN and its components in Min’an section

NO3-N濃度在1985—2021年整體呈極顯著上升趨勢,其中1985—1992年和2012—2021年NO3-N濃度呈顯著上升趨勢,1993—1999年NO3-N濃度呈顯著下降趨勢,2000—2011年NO3-N濃度無顯著升降趨勢。此外,除2000—2011年外的3個時段中,有LOWESS的檢驗結果與無LOWESS的檢驗結果不一致,表明徑流變化對NO3-N濃度變化存在顯著影響。1985—1992年徑流變化對NO3-N濃度變化起到弱化作用,而1993—1999、2012—2021年徑流變化對NO3-N濃度變化則起到增強作用。

NH3-N濃度在1985—2021年整體呈極顯著下降趨勢,但分時段條件下,其僅在2012—2021年呈極顯著下降趨勢,其余時段均呈上升趨勢但不顯著。此外,有LOWESS的檢驗結果與無LOWESS的檢驗結果一致,表明徑流對NH3-N濃度變化影響不顯著。

DIN濃度在1985—2021年整體呈上升趨勢但不顯著,其中1985—1992年和1993—1999年DIN濃度分別呈顯著上升和顯著下降趨勢,其余時段升降趨勢不顯著。此外,同NO3-N類似,DIN存在有LOWESS的檢驗結果與無LOWESS的檢驗結果在1985—1992年和1993—1999年不一致的現象,表明徑流變化同樣對DIN濃度變化存在影響。

2.4 DIN各組分所占比重年際變化

1985—2021年閩安斷面DIN各組分所占比重變化較明顯(表3和圖4)。其中NO2-N在DIN中所占比重最小,大部分時段比重均低于2%,僅在個別年份(1987、1998、2004、2008年)比重略高但均未超過5%,整體呈波動狀態,這同天然水體中NO2-N熱力學性質不穩定,容易進一步轉化為NH3-N或NO3-N有關[26]。NO3-N在DIN中所占比重最大,1985—1998年為61.9%~76.6%。1998年以后NO3-N比重逐漸增大,至2021年已達90.4%,表明NO3-N已成為DIN的主要組成部分。NH3-N在DIN中所占比重較小,1985—1996年為21.8%~35.0%。1999年起NH3-N比重逐漸減小,至2021年降至8.7%,為監測以來的最小值。

表3 1985—2021年閩安斷面DIN各組分比重的特征值Table 3 Characteristic percentage of DIN components in Min’an section during 1985-2021

圖4 1985—2021年閩安斷面DIN組分年際變化Fig.4 Secular variation of DIN components in Min’an section during 1985-2021

2.5 DIN及其各組分的入海通量

1985—2021年閩江DIN及其各組分的入海通量見表4和圖5,其中DIN入海通量除在個別極端洪旱年份波動較大外整體呈緩慢增加的特征。各組分中NO2-N入海通量以2005年為拐點呈先增加后減少的特征。NO3-N入海通量整體呈上升特征,尤其是2019—2021年NO3-N入海濃度快速升高導致年徑流量明顯下降時,入海通量也未出現明顯減少。NH3-N入海通量以1998年為拐點呈先增加后減少的特征,并在2010年以后降幅擴大。由于入海通量是濃度和徑流量共同作用的結果,這些因素的不同組合使入海通量的年際變化比較復雜[27]。但整體上,DIN及其各組分入海通量變化在短期內同徑流量變化比較接近,而在長期則同濃度變化更為類似。這主要是由于入海通量的變化主要取決于各因素中變幅較大者,短期內,徑流量的年際變化同濃度變化相比較為劇烈,而在長期來看,濃度變化幅度則遠大于年徑流量變化幅度(即距平百分率)。

表4 1985—2021年閩江DIN及其各組分入海通量特征值Table 4 Characteristic flux data of DIN and its components in Minjiang River during 1985-2021

圖5 1985—2021年閩江DIN及其各組分入海通量變化情況Fig.5 Flux variations of DIN and its components in Minjiang River during 1985-2021

3 討論

3.1 DIN及其各組分同源排放的關系

閩江流域DIN的來源較復雜,能夠對入海斷面構成影響的主要有工業源、農業源、生活源和大氣源。其中大氣源對流域氮輸入貢獻率較小(約為16.27%),農業源氮肥施用貢獻率最大(約為55.2%)[28]。此外,工業源和生活源具有影響局部性和集中性的特點,20世紀90年代,閩江流域大規模的水利開發建設基本完成后,下游水質同中上游工業和生活污染排放的直接響應關系不再顯著。因此該研究重點討論全流域農業源(以含氮化肥施用情況為例)以及下游福州市(以下如無特別指明,均指福州市在閩江流域范圍內的行政區域)工業源和生活源的變化對入海斷面DIN濃度及入海通量演變的影響。

閩江流域含氮化肥施用量年際變化情況及福州市人口與城市污水處理率變化情況見表5和圖6、圖7。

表5 福州市閩江流域人口及城市污水處理率變化情況Table 5 Variations of population and percentage of municipal sewage treatment of Fuzhou in Minjiang River basin

注:含氮化肥數據中,復合肥按照含氮33%的系數進行折純計算。圖6 閩江流域含氮化肥施用量年際變化情況Fig.6 Variation of nitrogenous fertilizer application in Minjiang River basin

圖7 福州市閩江流域人口及污水處理率變化情況Fig.7 Variations of population and treatment rate of domestic sewage of Fuzhou in Minjiang River basin

整體上,閩江入海斷面DIN及其各組分濃度演變同流域環境變化及下游福州市含氮污水排放存在較好的時間關聯性。具體體現在1985—1992年,閩江流域含氮化肥施用量快速增長,短短8年時間增幅近1倍。同時,下游福州市在改革開放初期工業迅猛發展,而污水處理設施建設基本處于空白,絕大多數醫藥、食品、造紙和化學工業企業產生的高濃度含氮工業廢水伴隨生活污水直排或通過城市內河排入閩江[29],導致該時期DIN所有組分濃度均呈上升特征。

1993—1999年,閩江流域含氮化肥施用量增長放緩并在1997年后穩定為16×104~17×104t。這一時期,閩江流域水環境綜合整治起步,1996年福州市第一座(也是全流域第一座)城市污水處理廠投入運行[29]。此外,隨著水口電站等梯級水利設施建成,水庫蓄水在改變污染物-徑流之間關系的同時,對中上游排放的污染物起到一定的滯留及緩沖作用[30-32]。在這些因素的共同作用下,該時期入海斷面NH3-N濃度上升的態勢得到遏制,NO3-N濃度甚至有所下降。但該趨勢有時也會受到短期事件的影響,如1998年8月福州開始實施引閩江水沖刷內河工程以改善城市水質[33]。引水初期大量內河黑臭水體排入閩江,導致當年入海斷面NO2-N和NH3-N濃度分別出現監測以來的最大值和次大值,并在1999年基本回落至原有水平。

2000—2011年,閩江流域含氮化肥施用量繼續保持穩定。同時下游福州市開始實施“東擴南進”發展戰略。原先位于城市建成區的工業企業根據要求陸續搬遷至流域外的濱海港口工業集中區,至2011年重污染企業全部遷出,導致工業源NH3-N劇減。而隨著外來人口大量涌入及城鎮化水平不斷提高,生活源在含氮廢水排放量中的比重大大增加。另一方面,福州市城市污水處理率在此時期大幅提高,但早期大部分污水處理廠只對NH3-N提出排放控制要求。其采用的好氧處理工藝僅能將NH3-N轉化為NO3-N,而總氮(TN)的去除率只有10%~20%[34]。這就導致12年間隨著城市污水處理量不斷增加,NO3-N排放逐漸增多,從而引起入海斷面NO3-N濃度緩慢升高,呈現出同NH3-N濃度變化背離的現象。

2011—2021年,受農業結構調整以及化肥使用量零增長減量化行動[35]的影響,全流域含氮化肥施用量自2016年起迅速下降。由于閩江流域降雨量大,土壤對肥料滯留能力弱,氮素流失較快[36],導致減量化行動實施后進入流域水體的氮通量迅速減少。與此同時,福州市外來人口快速增長,并在2015年突破100萬,導致生活源含氮廢水排放量進一步增加。由于該時期NH3-N已被列入國家生態環境保護規劃的主要污染物排放總量控制指標,福州市在城市污水處理率提高較為有限的情況下對污水處理廠進行提標改造。目前大部分污水處理廠出水可達《城鎮污水處理廠污染物排放標準》(GB 18918—2002)規定的一級A標準,個別新建廠出水甚至能達到《地表水環境質量標準》(GB 3838—2002)中規定的Ⅳ類標準。然而,根據對不同地區污水處理廠氮排放的報道[37-38],NH3-N平均去除率達97.8%,但TN去除率普遍不足80%,水廠出水NO3-N在TN中比重接近90%[37],表明除NH3-N效率高,除NO3-N能力不足的情況普遍存在。受這些因素的影響,在2000—2011年形成的NO3-N上升、NH3-N下降的趨勢在2011—2021年被進一步擴大,當2020、2021年閩江處于枯水年份時,這種現象表現得尤為明顯。因此,含氮化肥施用量減少及對NH3-N的強化治理是引起2011—2021年NH3-N比重顯著下降、NO3-N比重顯著上升的主要原因。

3.2 DIN及其各組分同徑流量之間的關系

SK檢驗結果表明,DIN各組分中NO2-N和NH3-N濃度變化在4個時段受徑流影響均不顯著,而NO3-N濃度變化在1985—1992、1993—1999、2011—2021年均顯著受到徑流影響??紤]到NO2-N和NH3-N濃度相對較低,且在河流、湖庫輸送過程中參與了較多生化過程,因此其同徑流量關系不明顯。而NO3-N是DIN組分中最穩定的形態,其濃度變化可能受徑流影響較大。為驗證徑流影響的方式,采用SK檢驗附帶的流量校正方程對NO3-N濃度和徑流量關系進行擬合,結果(表6)表明,1985—1999年擬合的方程為面源代表方程,而2011—2021年擬合的方程為點源代表方程,即存在從面源向點源轉變的現象,其原因可能是閩江下游福州市自21世紀以來環境治理設施不斷完善,越來越多的面源污染被收集至集中處理設施進行處理,導致污染排放方式發生轉變。

表6 不同時段NO3-N同徑流量的擬合結果Table 6 Regression results of NO3-N and flux in different time sections

3.3 DIN及各組分同其他突發環境事件的關系

雖然趨勢分析表明DIN及其各組分同污染排放關聯性整體較好,但研究發現部分突發環境事件對閩江DIN入海濃度及通量在短期同樣存在較大影響,導致個別數據出現“異常值”。通過對37年間閩江各類突發事件的梳理和分析,發現水口庫區生態災害和臺風對DIN短期影響比較突出。水口電站建成后,雖然減少了氮磷營養鹽向下游和河口的輸送,但同時也造成了其在庫區積累[30],為入侵水生植物生長創造了良好的條件。21世紀開始,水口庫區水葫蘆暴發時有發生,其中2003—2004年尤為嚴重。由于當時打撈能力有限,高度聚集的水葫蘆在庫區自生自滅并引起下游水質急劇惡化,致使入海斷面NH3-N在這2年中均呈現較高濃度[39-40]。

福建處于臺風多發區。臺風帶來的暴雨影響時間雖短,但其對河流物質輸送存在顯著影響[41]。2016年閩江下游因臺風降雨影響出現了較大規模的洪災,致使全年入海徑流量超過1998年閩江洪水時的水平,進而導致NO3-N和DIN入海通量達到自1985年監測以來的最大值。臺風暴雨等短期因素帶來DIN的大量輸入,可導致近海營養鹽局部失衡,繼而引發赤潮等海洋生態災害[42],因此對于極端條件下DIN的短期變化需要給予格外關注。

4 結論與建議

1) 1985—2021年的水質監測結合SK檢驗表明,閩江入海斷面DIN濃度整體呈上升趨勢但不顯著。DIN各組分中NO2-N和NH3-N濃度分別呈顯著和極顯著下降趨勢,而NO3-N濃度整體呈極顯著上升趨勢。37年間NH3-N在DIN中比重大幅降低,而NO3-N的比重則相反,目前已成為DIN的主要組成部分。

2) DIN及其各組分入海通量的長期變化與其各自濃度變化趨勢較接近。其中DIN入海通量為3.59×104~14.85×104t,在37年間緩慢增加。DIN各組分中,NO2-N和NH3-N入海通量呈先增加后減少的特征,而NO3-N入海通量整體呈上升的特征。從長期來看,DIN及其各組分濃度與流域環境演變及下游福州市污染排放存在較好的時間關聯性,但短期易受臺風等突發環境事件影響,需要格外關注。

3) 針對NO3-N濃度及其在DIN中比重明顯增大的情況,建議不斷完善城市污水處理廠處理工藝,逐步提高廢水脫硝處理效率,并進一步提升城市污水處理廠出水綜合利用水平。同時結合資源環境承載力適度控制城市人口規模,以遏制氮排放進一步增加。此外,鑒于臺風等短期事件對入海通量影響較大,建議進一步完善水質自動站等在線監測手段,以提高江河污染物入海通量的監測精度。

猜你喜歡
閩江入海徑流量
閩江雨情
金橋(2021年11期)2021-11-20
少兒美術(快樂歷史地理)(2020年5期)2020-09-11
沈葆楨題閩江仰止亭
高中地理校本課程的開發與實施——以閩江環境保護校本開發為例
上天入海我主沉浮
水文比擬法在計算河川徑流量時的修正
SCS模型在紅壤土坡地降雨徑流量估算中的應用
資江流域徑流量演變規律研究
上天入海 與夢同在
91香蕉高清国产线观看免费-97夜夜澡人人爽人人喊a-99久久久无码国产精品9-国产亚洲日韩欧美综合