龔發(fā)露，王裕成，蘭 佳，羅純良，羅瀲蔥，李慧赟，張如楓

(1.云南大學(xué)生態(tài)與環(huán)境學(xué)院，云南 昆明 650500； 2.杭州市生態(tài)環(huán)境局淳安分局，浙江 杭州 311700；3.昆明市滇池高原湖泊研究院，云南 昆明 650500； 4.中國科學(xué)院南京地理與湖泊研究所，江蘇 南京 210008)

溫躍層是深水湖泊中常見的物理現(xiàn)象，其變化規(guī)律和分布特征極大地影響著湖泊生態(tài)系統(tǒng)的結(jié)構(gòu)和功能。溫躍層會阻礙表層與底層水體之間的物質(zhì)和能量交換，降低底層水體的氧含量，導(dǎo)致底泥污染物釋放，致使水體營養(yǎng)鹽濃度升高。溫躍層的厚度、深度與持續(xù)時(shí)間受湖泊形態(tài)特征、湖體水質(zhì)和氣象因子等因素的影響，如湖泊深度和表面積大小是決定湖泊是否有條件形成溫躍層的重要因素[1]。在湖泊熱力分層穩(wěn)定期間，溫躍層深度與表層水溫存在顯著負(fù)相關(guān)關(guān)系，與透明度和水深存在顯著正相關(guān)關(guān)系[2]。氣溫升高可導(dǎo)致湖泊溫度分層時(shí)間延長，湖泊穩(wěn)定性增強(qiáng)。水位作為湖泊最基本的水文要素之一，對溫躍層的結(jié)構(gòu)和持續(xù)時(shí)間具有重要影響。水位下降后，中下層水體可接收更多的太陽輻射而升溫，上下水層溫差減小，分層強(qiáng)度減弱[3]。因此，如果人為大量排水致使水位下降，會在短期內(nèi)迅速降低水體穩(wěn)定度，可能使水體提前進(jìn)入混合期[4]。強(qiáng)降雨和入湖流量的異常增多，可導(dǎo)致水位迅速上漲并使穩(wěn)定的水體分層結(jié)構(gòu)遭到破壞[5]。夏季水位的大幅波動會對溫躍層結(jié)構(gòu)產(chǎn)生顯著影響[6]。在水位升降對溫躍層的影響方面已有很多研究成果，但鮮有涉及極端水位(年內(nèi)最高水位與最低水位之差)變化條件下溫躍層特征變化的研究。受氣候變化影響，從長時(shí)間尺度分析，我國長江中下游地區(qū)暴雨發(fā)生次數(shù)呈一定的增加趨勢，相應(yīng)會造成湖泊水位波動幅度和頻率的變化[7]。因此，研究極端水位變化條件下的溫躍層結(jié)構(gòu)特征，對水庫水環(huán)境綜合管理、水利工程調(diào)度和湖泊生態(tài)系統(tǒng)穩(wěn)定維持具有重要意義。

水位的監(jiān)測手段主要包括水位原位監(jiān)測、衛(wèi)星測高、遙感提取、模型模擬[8-10]等，其中水文水動力模型和數(shù)理統(tǒng)計(jì)模型等模型模擬方法運(yùn)用較為普遍。例如：李育等[11]采用瞬態(tài)氣候演變模型和湖泊水量平衡模型對東亞及中亞地區(qū)全新世以來湖泊有效水分變化進(jìn)行了模擬，分析了不同區(qū)域湖泊水位的主要影響因素；李瓊等[9]基于水位-流量轉(zhuǎn)換模型實(shí)現(xiàn)了宜昌市沮璋河洪峰水位的實(shí)時(shí)預(yù)報(bào)；華祖林等[12]從河流形態(tài)、水生生物和水質(zhì)等多維度出發(fā)，計(jì)算了秦淮河不同河段的生態(tài)水位。這些研究采用的方法為湖泊水位計(jì)算提供了良好的借鑒。水動力學(xué)模型DYRESM(dynamic reservoir simulation model)具有參數(shù)較少且率定較易的優(yōu)點(diǎn)，在湖泊水位研究中得到廣泛運(yùn)用。例如：Fadel等[13]通過DYRESM精確地模擬了Karaoun水庫的水位波動情況，模擬值與觀測值之間的誤差僅為0.4～0.7 m；在缺少完整入流數(shù)據(jù)時(shí)，Weinberger等[14]基于水量平衡法，通過DYRESM計(jì)算出了Ammersee湖1975—2007年的逐日入流量；Gu等[15]利用DYRESM模擬了海南省松濤水庫2016年3月至2017年3月的水位與水溫，發(fā)現(xiàn)其能很好地反映水庫的水位變化與溫度分層情況。此外，DYRESM還被廣泛應(yīng)用于水體溫度垂直分層的模擬，如加拿大Clearwater湖[16]、美國Mendota湖[17]和我國太湖[18]等。

目前對千島湖溫度分層及其影響因素的研究較多[19-21]，但主要集中于氣象強(qiáng)迫因子(氣溫與降水)和湖體特征(深度與透明度)對溫躍層的影響，很少考慮極端水位變化情況下的溫度分層規(guī)律。因此，本文以研究千島湖溫躍層對極端水位變化的響應(yīng)為目的，采用DYRESM模擬不同水位變化情景下的水溫分層情況，并通過湖泊分析程序Lake Analyzer(LA)計(jì)算并分析溫躍層特征量和湖泊穩(wěn)定性參數(shù)在不同水位變幅情況下的變化規(guī)律，以期為千島湖水環(huán)境保護(hù)提供參考。

1 研究區(qū)概況

千島湖(29°22′N～29°50′N，118°34′E～119°15′E)位于浙江省杭州市淳安縣境內(nèi)，是攔蓄新安江而形成的人工深水湖泊。自1960年完成蓄水至2016年，千島湖年均水位96.67 m(黃?；?，年平均水位變幅為9.91 m，其中2007年水位變幅最小(4.29 m)，1973年最大(22.57 m)。2000年以來年平均水位變幅為7.01 m，最大水位變幅出現(xiàn)在2014年，為11.49 m[22]。多年平均出入湖水量分別為91.7億m3與94.5億m3。千島湖入湖河道共有30多條，其中新安江、東源港、武強(qiáng)溪、楓林港和云源港5條為主要入湖河道(圖1)，共占千島湖總?cè)牒髁康?3.57%[23]，千島湖主要出流為下游的千島湖大壩。流域內(nèi)降水集中于4—8月，多年平均降水量為1 733 mm，年均氣溫為17.3℃[24]，年平均日照時(shí)數(shù)為1 951 h，年蒸發(fā)量為1 362 mm[22]。

圖1 千島湖流域及湖體水質(zhì)和氣象監(jiān)測點(diǎn)位Fig.1 Map of Lake Qiandao catchment and sites of water quality monitoring and weather station

2 研究數(shù)據(jù)與方法

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

研究所用數(shù)據(jù)來源及其屬性見表1，主要包括千島湖出入湖河道流量、水位和氣象數(shù)據(jù)。

表1 研究數(shù)據(jù)

2.2 研究方法

2.2.1 水動力學(xué)模型DYRESM

DYRESM是由西澳大利亞大學(xué)研發(fā)的一維水動力模型，可模擬湖庫中溫度、鹽度及密度的垂向分布，其輸入數(shù)據(jù)包括湖底高程、出入湖河道水量和氣象數(shù)據(jù)，輸出參數(shù)為不同深度水層的溫度和鹽度(適用于海洋及咸水湖)，可與水生態(tài)動力學(xué)模型(computational aquatic ecosystem dynamics model，CAEDYM)耦合，進(jìn)行水質(zhì)和生態(tài)系統(tǒng)模擬。

模型含有的文件包括：物理參數(shù)(.par)、配置(.cfg)、入流(.inf)、出流(.wdr)、氣象(.met)，其中出入流和氣象文件中均為逐日數(shù)據(jù)。物理參數(shù)文件中包含的參數(shù)見表2。氣象文件由平均氣溫、太陽輻射量、云量、日累計(jì)降水量、平均風(fēng)速和平均水汽壓組成。

表2 DYRESM主要物理參數(shù)

a.模型入流數(shù)據(jù)。因缺乏河道流量監(jiān)測數(shù)據(jù)，本文采用新安江模型[31]計(jì)算。張倚銘等[32]采用漁梁水文站2006年1月1日至2012年12月31日和2013年1月1日至12月31日兩個(gè)時(shí)間段的逐日監(jiān)測流量數(shù)據(jù)分別對新安江模型進(jìn)行了率定與驗(yàn)證，率定期與驗(yàn)證期的納什效率系數(shù)分別為0.98與0.95，水量平衡誤差百分比均為1%，說明新安江模型在漁梁站及以上子流域內(nèi)模擬精度較高。而漁梁站以上子流域的坡度、土地利用情況、氣候條件等與千島湖其他子流域非常相似，因此本文直接使用通過漁梁站監(jiān)測數(shù)據(jù)率定后的新安江模型，計(jì)算了2006年1月1日至2020年12月31日新安江、東源港、武強(qiáng)溪、楓林港和云源港5條主要入湖河道的逐日入湖流量，入流水溫則為實(shí)測值。

b.模型驗(yàn)證。為確定DYRESM對水溫和水位的模擬精度，采用大壩前國控點(diǎn)2018年1月1日至12月31日的實(shí)測水位和水溫?cái)?shù)據(jù)對模型進(jìn)行驗(yàn)證，對比分析水位和水深0.5 m、10.5 m、20.5 m的水溫模擬結(jié)果與實(shí)測結(jié)果(圖2)，模擬值和實(shí)測值之間的相關(guān)系數(shù)分別為0.99和0.99、0.99、0.92，均方根誤差分別為0.29 mm和1.08℃、1.10℃、1.47℃，說明模型對千島湖水位和水溫的模擬精度較高。

圖2 千島湖大壩前國控點(diǎn)水位和水溫模擬值與實(shí)測值對比Fig.2 Comparison of simulated and measured daily water level and water temperature at Dabaqian

2.2.2 湖泊分析程度LA

LA是計(jì)算湖泊(水庫)穩(wěn)定性參數(shù)和熱力分層指標(biāo)的程序。程序輸入不同監(jiān)測頻率的水溫、平均風(fēng)速、水深等指標(biāo)，輸出溫躍層表層深度(MetaT)、溫躍層底層深度(MetaB)、湖泊數(shù)(Ln)、韋德伯恩數(shù)、施密特穩(wěn)定度(St)、浮力頻率(N2)等參數(shù)。本文主要采用Ln與St作為水體穩(wěn)定度參數(shù)，其中St反映水體從現(xiàn)狀態(tài)轉(zhuǎn)換成完全混合狀態(tài)所需要的能量，在一定程度上可以度量溫躍層的強(qiáng)度，該值越大表明湖泊分層越穩(wěn)定，越接近0表明湖泊越趨于混合狀態(tài)；Ln表示外界因素(如風(fēng)力、出入流和人為活動)對水體分層穩(wěn)定性的干擾程度，其值越小代表水體越趨于混合。本文所用溫躍層厚度(MetaTH)為MetaT與MetaB的差值。

3 模擬結(jié)果與分析

3.1 情景設(shè)置

降雨和蒸發(fā)雖然能在一定程度上影響湖泊水位的變化，但引發(fā)千島湖水位變動的主要原因是入湖河道流量和出庫流量的波動，因此在配置模擬情景時(shí)，主要根據(jù)千島湖出入湖河道歷史流量平均值(代表多年平均情況)與水位波動幅度較大情況下的實(shí)際值(代表極端水位變化情況)進(jìn)行設(shè)定，共設(shè)置兩種情情景：情景1代表千島湖極端水位波動情況，發(fā)生于2013年7月1日至2014年7月31日，該情景中5條主要入湖河道的流量通過新安江模型計(jì)算得到；情景2中5條主要入湖河道的流量為2006—2020年的多年平均實(shí)測入湖流量。兩種情景均采用相同的初始水位(104.40 m)，模擬時(shí)段均為2013年7月1日至2014年7月31日，氣象數(shù)據(jù)均為該時(shí)段內(nèi)的實(shí)測數(shù)據(jù)，時(shí)間步長設(shè)為1 h。

情景1水位最大變幅達(dá)到11.49 m，2013年7月1日水位為104.40 m，經(jīng)歷231 d水位下降至95.01 m，其中2013年7、8月水位下降速度較快，下降速度分別達(dá)到0.09 m/d與0.07 m/d，再歷經(jīng)152 d于2014年7月18日達(dá)到最高水位106.50 m；情景2最高水位與最低水位相差5.48 m，從2013年7月1日的104.40 m下降到2014年2月16日的98.92 m，再歷經(jīng)153 d上漲3.58 m達(dá)到102.50 m。

3.2 模擬結(jié)果

圖3、圖4與圖5分別為兩種情景水溫隨水位的逐日變化、水體穩(wěn)定度參數(shù)和溫躍層范圍的模擬結(jié)果。由圖3(a)、圖4(a)和圖5(a)可知，情景1溫躍層在7—9月相對穩(wěn)定，10—11月逐漸消失，12月水體完全混合。2013年秋冬季，MetaT與MetaB逐漸增大，平均值分別為21.59 m與25.18 m；而MetaTH整體較小，平均厚度僅為2.53 m，低于同時(shí)期情景2平均厚度(4.75 m)，最大厚度出現(xiàn)于2013年8月(5.02 m)；St與Ln總體呈下降趨勢，平均值分別為5 252.00 J/m2與87.00。2014年1—6月，情景1中水位逐漸上升，而MetaT基本不變，維持在水深6 m左右，MetaB不斷增大，5月MetaTH成為該時(shí)段的最大值，為8.45 m；7月溫躍層位置明顯加深，位于水深19.8～24.5 m；2014年下半年St和Ln平均值分別為2 821.00 J/m2與25.90。

圖3 千島湖大壩前國控點(diǎn)水溫隨水位的變化Fig.3 Chang of water temperature with depth at Dabaqian

由圖3(b)、圖4(b)與圖5(b)可知， 2013年7—8月情景2溫躍層較淺，僅存在于水深12～20 m之間，9—12月MetaT與MetaB的平均值分別為25.90 m和29.19 m，溫躍層位于水深19.0～35.0 m之間；2013年7—12月，MetaTH逐漸變小，7月平均厚度最大(7.19 m)，12月最小(0.2 m)。2014年1—7月，MetaT從3.85 m變?yōu)?9.50 m，MetaB從6.10 m增至23.73 m，其中，6—7月二者變化最明顯，分別增加9.95 m和7.63 m。

圖4 千島湖大壩前國控點(diǎn)水體穩(wěn)定度參數(shù)Fig.4 Water stability parameters at Dabaqian

圖5 千島湖大壩前國控點(diǎn)溫躍層范圍Fig.5 Thermocline ranges at Dabaqian

情景1溫躍層平均水位98.86 m，St與Ln的平均值分別為4 070.26 J/m2和55.93；情景2溫躍層平均水位為101.10 m，St與Ln的平均值分別為4 189.45 J/m2和51.04。情景1與情景2相比，平均水位低2.24 m，St均值低119.19 J/m2，可見極端水位變化情況下，水體穩(wěn)定度也相應(yīng)降低。

3.3 典型時(shí)段溫度垂直分層差異

3.3.1 溫躍層持續(xù)時(shí)間差異

分析兩種情景溫躍層的形成與消失時(shí)間可知，情景1溫躍層從模擬開始時(shí)間(2013年7月1日)持續(xù)到2013年11月20日，在2014年3月23日重新出現(xiàn)，水體混合期122 d；情景2溫躍層自模擬開始時(shí)間持續(xù)到2013年12月5日，于2014年3月16日重新出現(xiàn)，水體混合期102 d。由St與Ln計(jì)算結(jié)果可知，2013年11月至2014年3月(水體混合期)情景1平均水位比情景2低4.09 m，St平均值比情景2小71.70 J/m2，說明冬季水位降低時(shí)湖泊混合相對較強(qiáng)，穩(wěn)定度較差，而處于高水位且水位波動幅度較小時(shí)穩(wěn)定度較高，溫躍層相對穩(wěn)定。

3.3.2 溫躍層結(jié)構(gòu)差異

模擬時(shí)段內(nèi)，2013年7月20—30日，情景1與情景2平均水位分別為99.76 m與103.41 m；2014年7月20—30日，情景1與情景2平均水位分別為106.35 m與101.98 m，水位差異達(dá)到最大。選擇兩情景在上述時(shí)段內(nèi)的水體分層情況進(jìn)行對比，分析水位差異最大時(shí)段內(nèi)的溫躍層結(jié)構(gòu)特征，結(jié)果如圖6所示。2013年7月20—30日(圖6(a))，情景1溫躍層出現(xiàn)在水深10～35 m，MetaT與MetaB處水溫相差18.07℃；情景2 MetaTH較小，水層的溫度變化梯度較大，從水深10～27 m水溫下降15.6℃。2014年7月20—30日(圖6(b)),情景2從表層到水深15 m的水層溫度梯度變化較小，水深15～35 m出現(xiàn)溫躍層，水體溫度下降15.2℃，而水深35 m以下水體無明顯溫度變化，此時(shí)水體為單一、穩(wěn)定的溫躍層結(jié)構(gòu)。情景1中由于大氣降水和河道徑流的匯入，水位上漲較快，使得表層水體發(fā)生了一定程度的擾動，水深15 m以上的水層水溫明顯升高，且在水深30～40 m出現(xiàn)15℃左右的“均溫層”，水深40～50 m水體溫度梯度為0.29 ℃/m，溫躍層再次出現(xiàn)，即水位上漲之后MetaB相對情景2有所增大，且溫躍層被水深30～40 m處的均溫層切斷，水體呈從表層到深層依次為“混合層、溫躍層、均溫層(30～40 m水深)、溫躍層、均溫層”的雙溫躍層結(jié)構(gòu)。

圖6 模擬水溫隨深度的變化Fig.6 Changes of simulated water temperature with depth

4 討 論

千島湖是一個(gè)典型的亞熱帶單混合水體，正常水位情況下，水體從3月下旬出現(xiàn)分層，至12月初分層消失，溫躍層持續(xù)時(shí)間較長。本文所設(shè)置的兩種情景中，水位與St、Ln均呈顯著正相關(guān)關(guān)系，情景1相關(guān)系數(shù)分別為0.81和0.41(樣本量為396)，情景2相關(guān)系數(shù)分別為0.79和0.53。因此，水位越高，MetaTH越大，湖泊水體越穩(wěn)定，分層越明顯。極端水位變化情況下，對溫躍層持續(xù)時(shí)間的影響主要體現(xiàn)在水體穩(wěn)定度降低和混合期提前，縮短了溫躍層持續(xù)時(shí)間。情景1秋季溫躍層比情景2提前15 d結(jié)束，而春季溫躍層開始日期僅相差5 d，說明水位對水體分層的影響在秋季更明顯，這與Wang等[4]對廣東省流溪河水庫的研究結(jié)果一致，進(jìn)一步說明了本文所采用的模型預(yù)測的可靠性。

在國內(nèi)外具有相同氣候背景的其他湖泊中，溫躍層持續(xù)時(shí)間具有與千島湖類似的特征。例如：韓國Yongdam水庫(最大水深70 m)，其溫度分層時(shí)間為4—11月[33]；云貴高原烏江流域的百花湖水庫(最大水深45 m)分層時(shí)間為4—10月[34]；流溪河水庫(最大水深73 m)溫躍層出現(xiàn)在4—12月[4]。但是，包括千島湖在內(nèi)的這些湖泊，受氣候變化的影響，隨著暴雨頻率和強(qiáng)度的增強(qiáng)，極端降水事件出現(xiàn)的概率會增大，意味著湖泊徑流量和水位波動會更加激烈和復(fù)雜化[7,35]，對溫躍層的沖擊將不可忽視。而在這些極端事件的發(fā)生過程中，很難依靠人工監(jiān)測獲取水位和水溫垂直分布的信息，因此依靠浮標(biāo)探頭進(jìn)行自動監(jiān)測和模型模擬將是獲取這些信息的有效手段。相比于探頭自動監(jiān)測技術(shù)，數(shù)值模型具有經(jīng)濟(jì)和高效的特點(diǎn)，且所獲數(shù)據(jù)具有時(shí)空尺度上的連續(xù)性。因此，在湖泊物理特性的研究方面(例如水溫和水量)，模型模擬技術(shù)是對人工監(jiān)測和自動監(jiān)測的有效補(bǔ)充。

模型模擬結(jié)果表明，千島湖水位較高時(shí)水體出現(xiàn)雙溫躍層結(jié)構(gòu)。通過對影響溫躍層的各種可能因子的主成分分析可知，影響MetaTH的第一主成分為入湖河道水溫和氣溫，貢獻(xiàn)率達(dá)44%，兩者與MetaTH相關(guān)系數(shù)分別為0.78和0.71(樣本量均為396)，而入湖流量、降水量和云量為第二主成分(貢獻(xiàn)率18%)?？梢娗u湖雙溫躍層結(jié)構(gòu)的形成不僅與氣候條件有關(guān)，也與水文條件有關(guān)。這種雙溫躍層結(jié)構(gòu)其實(shí)普遍存在于深水湖庫中，如海南松濤水庫[36]和二灘水電站水庫[37]均有發(fā)現(xiàn)；Liu等[21]研究表明，強(qiáng)降雨條件導(dǎo)致短期水量大量涌入，改變了表層水溫，導(dǎo)致水體出現(xiàn)暫時(shí)的雙溫躍層結(jié)構(gòu)。雙溫躍層結(jié)構(gòu)形成期間，河道徑流匯入與其密度相近的水層并發(fā)生混合和熱量交換，這在一定程度上使得深層水體能得到表層富氧水體的補(bǔ)充，為浮游動植物的生長生活提供溶解氧。同時(shí)分層強(qiáng)度的削弱也有利于水體中營養(yǎng)鹽的擴(kuò)散，降低了局部水域營養(yǎng)鹽的濃度。在千島湖，除暴雨過程外，其他氣候條件基本穩(wěn)定，雖然氣候變化可能導(dǎo)致氣溫緩慢上升，但這種效果在短期內(nèi)幾乎可以忽略。因此由暴雨事件、千島湖發(fā)電和水利調(diào)度造成的水文條件的急驟改變，進(jìn)而造成的對千島湖水體垂直溫度結(jié)構(gòu)的破壞以及對水質(zhì)的影響，值得更進(jìn)一步深入研究。

已有研究表明，重要的水動力過程(如極端水位變化)會引起上覆水和沉積物的理化性質(zhì)的變化，進(jìn)而影響甚至改變底泥污染物的釋放規(guī)律[38]。當(dāng)溫躍層存在時(shí)水體垂直交換受阻，DO的水平分層結(jié)構(gòu)穩(wěn)定，底部水體氮磷營養(yǎng)鹽、金屬離子及有機(jī)物等污染物在厭氧環(huán)境中釋放[39]，并可能在沉積物-水界面呈現(xiàn)一定程度的營養(yǎng)鹽累積效應(yīng)。極端水位變化情況下，冬季溫躍層提前消失，不穩(wěn)定的水體可能導(dǎo)致底層污染物向上層水體擴(kuò)散，上層水體中營養(yǎng)鹽富集可能加快浮游植物的生長，并增加水華暴發(fā)的可能性[40-41]，給飲用水源安全帶來影響[42]。因此，對類似于千島湖的大型深水湖庫，水文條件劇烈波動所造成的水環(huán)境的改變和對生態(tài)系統(tǒng)的沖擊不能低估。隨著自動監(jiān)測技術(shù)的愈發(fā)先進(jìn)和模擬技術(shù)的愈加成熟，這種影響和沖擊的定量評估會變得越來越容易，能給相關(guān)管理部門提供制定有關(guān)極端氣候事件應(yīng)對策略的科學(xué)依據(jù)。

5 結(jié) 論

a.水位劇烈波動條件下，溫躍層持續(xù)時(shí)間較正常水位波動情況下明顯縮短。情景1水體混合期為2013年11月20日至2014年3月23日，情景2則為2013年12月5日至2014年3月16日，情景1相對于情景2溫躍層持續(xù)時(shí)間縮短20 d。

b.高水位時(shí)，情景1湖體呈“混合層、溫躍層、均溫層(30～40 m水層)、溫躍層、均溫層”的雙溫躍層結(jié)構(gòu)，情景2則為“混合層、溫躍層、均溫層”單溫躍層結(jié)構(gòu)。

c.情景1平均水位較情景2低2.24 m，St均值低119.19 J/m2，表明極端水位變化條件下出現(xiàn)的低水位情況能顯著降低水體穩(wěn)定度。