柯賢敏,歐安鋒,王 瑋,牛富俊,高澤永

(1. 長安大學(xué)水利與環(huán)境學(xué)院,陜西 西安 710054；2. 長安大學(xué)旱區(qū)地下水文與生態(tài)效應(yīng)教育部重點(diǎn)實驗室,陜西 西安 710054；3. 中國科學(xué)院西北生態(tài)環(huán)境資源研究院凍土工程國家重點(diǎn)實驗室,甘肅 蘭州 730000)

青藏高原是中國多年凍土的代表,影響著中國乃至全球的氣候變化[1]。近年來,在人為活動和氣候變暖的影響下,青藏高原多年凍土逐步退化[2- 3]。多年凍土退化最直接的體現(xiàn)是凍土活動層厚度增大和熱喀斯特湖數(shù)量及面積的增大[4- 5]。熱喀斯特湖作為重要的熱源,對區(qū)域水文循環(huán)、碳循環(huán)以及生態(tài)環(huán)境[6- 7]具有重要影響。目前,國外開展的研究主要集中于熱喀斯特湖對地貌、區(qū)域環(huán)境及碳循環(huán)的影響[8- 9]。青藏高原熱喀斯特湖的研究主要集中于水熱過程、分布特征、湖水水文地球化學(xué)特征及其對工程穩(wěn)定性影響方面[10- 13],關(guān)于熱喀斯特湖與多年凍土相互作用的研究相對較少。深入了解熱喀斯特湖與多年凍土間的相互作用,有助于更好地理解多年凍土區(qū)的生態(tài)水文過程。

長期以來,學(xué)者們多通過遙感影像和現(xiàn)場監(jiān)測手段分析熱喀斯特湖的發(fā)展過程及其與凍土間的水熱過程[14]。然而,熱喀斯特湖與周圍多年凍土的長期相互作用十分復(fù)雜,涉及熱傳導(dǎo)、水分遷移、冰水相變等過程,短期的觀測數(shù)據(jù)無法揭示和預(yù)測熱喀斯特湖演化過程及其與凍土間的相互作用[15]。因此,數(shù)值模擬技術(shù)成為研究凍土水熱變化、熱喀斯特湖發(fā)育的重要手段。國外針對多年凍土區(qū)的水熱耦合模型主要有2個研究方向:一個是只考慮凍土內(nèi)部的水熱變化[16- 17],另一個是考慮在大氣- 植被- 土壤整個系統(tǒng)中的水熱變化,再加上降水和融雪等過程[18- 19]。目前為止,針對青藏高原熱喀斯特湖開展的模擬研究相對較少。林戰(zhàn)舉等[20]通過構(gòu)建熱傳導(dǎo)數(shù)值模型,預(yù)測了熱喀斯特湖對多年凍土的影響,結(jié)果表明湖的溫度對湖周圍凍土影響顯著；尹國安等[21]通過構(gòu)建一維熱傳導(dǎo)模型,模擬了不同深度的湖在氣候變暖背景下對多年凍土的影響；令鋒等[22]利用熱傳導(dǎo)模型計算了熱喀斯特湖不同擴(kuò)張速率對多年凍土退化的影響,結(jié)果表明湖的擴(kuò)張速率越大,對凍土退化影響越大；紀(jì)舒婷[23]利用數(shù)值模型模擬了凍土退化對湖泊儲量的影響,指出凍土退化增加了地下水入湖量；楊振等[24]利用移動網(wǎng)格技術(shù)構(gòu)建了熱喀斯特湖演化的熱傳導(dǎo)數(shù)值模型,分析了熱喀斯特湖深度和半徑對凍土的影響,并得出恒定的熱喀斯特湖邊界會高估對湖底凍土的影響、低估對湖附近凍土的影響。

上述研究多是基于熱傳導(dǎo)理論開展的,忽略了熱對流過程對凍土退化的貢獻(xiàn),這些方法只適用于地下水流速很小的環(huán)境。青藏高原多年凍土分布格局和動態(tài)變化與地下水密切相關(guān)[25],僅考慮熱傳導(dǎo)過程對凍土的影響是不全面的。為此,本文以北麓河盆地典型熱喀斯特湖區(qū)域為研究對象,結(jié)合野外監(jiān)測與數(shù)值模擬手段,對熱喀斯特湖與多年凍土間的相互作用展開研究,分析熱喀斯特湖的水均衡變化特征,探討地質(zhì)環(huán)境和氣候變暖對熱喀斯特湖和生態(tài)環(huán)境發(fā)展方向的影響。

1 研究區(qū)與研究方法

1.1 研究區(qū)概況

北麓河研究區(qū)位于青藏高原北麓河盆地(圖1),屬于青藏高原亞寒帶半干旱氣候。研究區(qū)年平均氣溫、降水量、蒸發(fā)量分別為-4 ℃、274.70 mm和1 421.40 mm。降水通常以雨、雪和冰雹的形式發(fā)生于6—8月。研究區(qū)內(nèi)地層起伏不大,自上而下依次為第四系全新統(tǒng)沖洪積層、上第三系湖相沉積層。全區(qū)分布厚度約為30～100 m的多年凍土,活動層厚度一般為1.60～3.40 m。地下水的排泄方式主要有蒸發(fā)、以泉的形式溢流于地表和向熱喀斯特湖及凍土層下水補(bǔ)給等。本文選取的熱喀斯特湖(BLH- B)位于北麓河盆地內(nèi)部,面積約為7 hm2,最大深度為2 m(圖1),湖周圍下墊面類型為高寒草地。

圖1 BLH- B湖溫度監(jiān)測點(diǎn)布設(shè)Fig.1 Distribution of temperature monitoring sites in BLH- B

1.2 野外試驗和數(shù)據(jù)

采用MRDT605型溫度測量儀對湖中和湖附近區(qū)域溫度進(jìn)行監(jiān)測。湖東北、西南和東南側(cè)各設(shè)置了3個監(jiān)測點(diǎn)(圖1),監(jiān)測深度為地面之下0.10 m、0.25 m、0.50 m、1.00 m,監(jiān)測頻率為4 h/次；在湖中布設(shè)了3個溫度監(jiān)測點(diǎn)(圖1),用于監(jiān)測湖底之下的地溫。湖底溫度和湖周圍地溫(地面之下0.10 m)分別用于模型湖底邊界和地表邊界的賦值。此外,湖周圍地溫監(jiān)測數(shù)據(jù)也用于對計算模型進(jìn)行驗證。

1.3 水- 冰- 熱耦合模型

1.3.1 主要控制方程

SUTRA能夠解決飽和非飽和條件下的溶質(zhì)運(yùn)移和地下水流動等問題。然而,SUTRA原始程序并不能對多年凍土區(qū)地下水的凍結(jié)和融化進(jìn)行刻畫。為此,McKenzie等[26]對SUTRA程序進(jìn)行了改進(jìn),得到能夠模擬飽和條件下地下水凍結(jié)與融化過程的SUTRA- ICE程序,該程序同時考慮了熱傳導(dǎo)和熱對流過程。SUTRA- ICE程序的地下水流動控制方程為

(1)

SUTRA- ICE程序中能量傳遞的實質(zhì)是冰、孔隙介質(zhì)和液體中的能量平衡,其方程可以表示為

(2)

1.3.2 模型結(jié)構(gòu)和邊界條件

根據(jù)模擬區(qū)實際地形坡度(約0.87°),利用水- 冰- 熱耦合計算程序,采用柱坐標(biāo)系建立二維剖面模型。模型橫向由外流域邊界到湖泊中心(X=0)的500 m徑向距離確定,縱向由深度確定,模型徑向截面如圖2所示。模型在垂向上共劃分了3大層,自上而下地層滲透率依次為5×10-13m2、1×10-13m2和1×10-14m2。模型中間層和底層共劃分了49小層,單層厚度為10 m。為了更好地展示活動層的變化,將頂層進(jìn)一步劃分為100層,單層厚度為0.10 m。此外,為了分析熱喀斯特湖周圍凍土的變化情況,將湖及湖岸附近20 m區(qū)域進(jìn)一步劃分為35列,每列寬2 m；其他區(qū)域劃分為43列,列寬為10 m。模型厚度為H=2πX。模型的頂部設(shè)定為空氣熱力學(xué)邊界,其溫度為大氣溫度,壓強(qiáng)為0 Pa；模型底部設(shè)定為地?zé)嵬窟吔?0.06 W/m2)；兩側(cè)為零通量邊界,不與外界發(fā)生水熱交換。模型西北角設(shè)定1個寬50 m、深2 m的湖,利用溫度和壓強(qiáng)邊界來刻畫湖泊。

圖2 模型結(jié)構(gòu)、網(wǎng)格剖分及邊界條件概化Fig.2 Figure showing the structure,grid,and boundary condition of model

2 結(jié)果及分析

2.1 模型驗證

構(gòu)建多年平均氣溫(-3 ℃)和地?zé)嵬抗餐饔孟碌姆€(wěn)定流模型,得到穩(wěn)定的多年凍土結(jié)構(gòu)。穩(wěn)定流模型的溫度和壓強(qiáng)計算結(jié)果作為非穩(wěn)定流模型的初始條件?；诜€(wěn)定流模型,在模型西北角增加1個深2 m、寬50 m的湖泊,根據(jù)2018—2019年熱喀斯特湖和地表溫度監(jiān)測數(shù)據(jù)分別對模型地表邊界和湖底邊界進(jìn)行賦值,并利用溫度監(jiān)測數(shù)據(jù)對計算模型進(jìn)行驗證。實測溫度與計算溫度組成的坐標(biāo)點(diǎn)均位于1∶1直線附近,相關(guān)系數(shù)(R2)為0.97～0.99,均方根誤差(ERMS)為0.48～1.31 ℃(圖3),計算結(jié)果與實測結(jié)果擬合較好。因此,該模型可用于模擬熱喀斯特湖的排泄量及其與多年凍土間的相互作用。模型參數(shù)如表1所示。

圖3 實測溫度與計算溫度的對比Fig.3 Comparison of the measured and calculated temperature

表1 模型最終參數(shù)表

2.2 多年凍土及活動層變化特征

為了揭示熱喀斯特湖長期作用下多年凍土和活動層的變化特征,在穩(wěn)定流模型的基礎(chǔ)上構(gòu)建了長達(dá)1 000 a的非穩(wěn)定流數(shù)值模型。該模型設(shè)置了12 000個時間步長,每個時間步長代表1個月。地表和湖底邊界溫度根據(jù)多年平均溫度進(jìn)行循環(huán)(圖4)。當(dāng)?shù)貙颖柡投刃∮?.05時(無顏色覆蓋區(qū)域),認(rèn)為地層未發(fā)生凍結(jié)；模型中的剩余液態(tài)水飽和度設(shè)定為0.05,即凍結(jié)的凍土中也有少量水分存在,故代表多年凍土存在的冰飽和度范圍為0.05～0.95。冰飽和度的增加(由紅色區(qū)域向藍(lán)色區(qū)域變化)代表了地層凍結(jié)程度增強(qiáng)(圖5)。

圖4 熱喀斯特湖和大氣溫度變化Fig.4 Temperature variation of thermokarst lake and air

湖水持續(xù)向多年凍土傳遞熱量,引起湖周圍凍土逐步退化,多年凍土面積逐漸減小(表2)。第200 a,湖底凍土已經(jīng)開始退化但尚未融穿；而熱喀斯特湖持續(xù)作用400 a,湖底形成了貫穿融區(qū)。貫穿融區(qū)形成前,地下冰限制了地下水的流動,地下水流速緩慢,熱傳導(dǎo)占主導(dǎo)作用，此時熱量主要在垂向上進(jìn)行傳遞,垂向上的凍土退化速度快于水平方向；貫穿融區(qū)形成后,熱對流占主導(dǎo)作用,熱量伴隨著地下水的流動進(jìn)行傳遞。在湖附近區(qū)域,多年凍土底界面因不斷受到地下水流的熱侵蝕而抬升，同時,地下水流速和湖的排泄量有所增大,而湖獲得的補(bǔ)給量減小。地表的排泄量及其獲得的補(bǔ)給量在貫穿融區(qū)形成前后基本穩(wěn)定,在模擬期末出現(xiàn)增長，由于熱喀斯特湖的影響范圍有限,距離湖較遠(yuǎn)的區(qū)域多年凍土保持穩(wěn)定；貫穿融區(qū)形成后,地下水流速增大,攜帶熱量的地下水對多年凍土的熱侵蝕作用增強(qiáng),活動層厚度增大了0.23 m。整個模擬時段內(nèi),活動層地下水保持近水平流動。在多年凍土的影響下,凍土層下水、凍土層上水、凍土層間水以及湖水間的水力聯(lián)系微弱,凍土層下水在靜水壓力作用下自東南向西北流動,流速自下而上逐漸增大。對比第200 a和第1 000 a計算結(jié)果可以得出,貫穿融區(qū)形成后,水流通道恢復(fù),凍土層下水的流向發(fā)生了改變,各水體間的水力聯(lián)系增強(qiáng),形成了新且穩(wěn)定的地下水循環(huán)系統(tǒng)。

圖5 不同時期多年凍土、流速矢量及活動層分布計算結(jié)果Fig.5 Distribution of permafrost,active layer,and velocity vector in different periods

表2 不同時段模型計算結(jié)果統(tǒng)計表

圖6截取了模型第400～420 a的計算結(jié)果,來反映湖邊觀測孔處(X=60 m)的活動層凍融循環(huán)過程。地表長期處于升溫、降溫循環(huán)過程中,故形成了厚度為3.35 m的活動層。11月至次年5月,地表溫度低于0 ℃,活動層處于凍結(jié)狀態(tài)；而6—10月,地表溫度高于0 ℃,活動層處于融化狀態(tài)。整個模擬時段內(nèi)年際間地表溫度未升高,故活動層厚度保持穩(wěn)定。地表熱量向深部傳遞的過程中逐漸被消耗,故活動層深部處于融化狀態(tài)的時間比淺部短,而處于凍結(jié)狀態(tài)的時間比淺部有所延長。此外,隨著深度的增加,深部多年凍土對地表溫度的響應(yīng)也表現(xiàn)出一定的滯后性。

圖6 湖邊觀測孔處活動層厚度計算結(jié)果Fig.6 Active layer thickness of lakeside observation well

2.3 熱喀斯特湖水均衡計算結(jié)果分析

圖7 熱喀斯特湖水均衡計算結(jié)果Fig.7 Results of thermokarst lake water balance

多年凍土和熱喀斯特湖的水熱過程勢必會引起地表邊界和湖底邊界的水均衡發(fā)生變化。湖底邊界獲得的補(bǔ)給量(進(jìn)入湖底邊界的水量)及其自身的排泄量(湖底邊界流出的水量)決定了湖底邊界的水均衡狀況。如圖7(a)所示,熱喀斯特湖獲得的補(bǔ)給量與其自身的排泄量同步變化。隨著湖底溫度的升高,流入湖底的凍土層上水水量增加,湖與凍土層下水間的水頭差減小,湖的排泄量減小，為了維持湖水位保持不變,湖獲得的補(bǔ)給量也同步減小；反之,湖水溫度降低過程中,湖的排泄量及其獲得的補(bǔ)給量均增加。在1個溫度周期內(nèi),湖表現(xiàn)為負(fù)均衡,即獲得的補(bǔ)給量小于自身的排泄量。湖向周圍凍土傳遞熱量需要一定的時間,從而導(dǎo)致湖的補(bǔ)給量和排泄量滯后于其溫度的變化。熱喀斯特湖在整個模擬時段內(nèi)也表現(xiàn)為負(fù)均衡,即一直處于排水狀態(tài)(圖7(b))。模型第285～388 a計算結(jié)果表明:隨著凍土含冰量的降低,地下水流動通道恢復(fù),熱喀斯特湖的排泄量逐漸增大,其獲得的補(bǔ)給量減小；熱喀斯特湖持續(xù)作用388 a時,湖底形成貫穿融區(qū),熱喀斯特湖的凈排泄量達(dá)到最大；貫穿融區(qū)形成后,凍土層上水、凍土層下水以及湖水間的水力聯(lián)系增強(qiáng),形成新的地下水循環(huán)系統(tǒng)；隨后,由于湖底邊界的壓強(qiáng)是恒定的,故熱喀斯特湖的排泄量開始減小并逐漸恢復(fù)至初始狀態(tài)。

3 討 論

3.1 不同地質(zhì)環(huán)境下的熱喀斯特湖發(fā)展方向

地層滲透率影響著地下水的徑流速度和湖的排泄量。圖8(a)顯示了不同地層滲透率(方案1—方案4的地層滲透率分別為1×10-12m2、1×10-13m2、1×10-14m2和1×10-40m2)地層環(huán)境下熱喀斯特湖凈排泄量的變化特征。隨著地層滲透率的增大,地下水流動受到的阻滯作用減小,熱喀斯特湖的排泄量逐漸增大。地層滲透率很小時(1×10-40m2),地下水近似處于靜止?fàn)顟B(tài),湖的排泄量接近于0 kg/s。湖底形成貫穿融區(qū)時,湖水將持續(xù)向深層地下水排泄,湖的排泄量有所增大,引起湖水進(jìn)一步流失。當(dāng)湖獲得的補(bǔ)給量得不到保證時,熱喀斯特湖將逐漸排干,最終形成新的凍土。北極地區(qū)地層巖性以砂卵石為主,地層滲透率大,熱喀斯特湖底凍土融穿后,湖的排泄量增大,而其獲得的補(bǔ)給量有限,故一些熱喀斯特湖逐漸干涸消失[27]。本次模擬的青藏高原熱喀斯特湖區(qū)域,地層巖性以砂巖和泥巖為主,地層滲透性弱,湖的排泄量相對較小。雖然湖底形成貫穿融區(qū)后,湖的排泄量有所增大,但增大程度有限。此外,氣候變暖加速凍土退化,湖獲得來自地下水的補(bǔ)給量增大,降水量也能得到保證,故熱喀斯特湖逐漸擴(kuò)張[28]。

圖8 不同滲透率地層環(huán)境下湖和地表排泄量計算結(jié)果Fig.8 Results of the lake and ground- surface discharge in different permeability environments

3.2 不同地質(zhì)環(huán)境下的生態(tài)環(huán)境效應(yīng)

天然條件下,降水量、蒸發(fā)量、地表徑流量與地下水徑流量處于動態(tài)均衡狀態(tài)。當(dāng)降水入滲補(bǔ)給量大于維持地表水位所需補(bǔ)給量時,地下水完全飽和處于蓄滿產(chǎn)流狀態(tài),形成地表徑流；而當(dāng)降水入滲補(bǔ)給量小于維持地表水位所需補(bǔ)給量時,潛水面將會下降。不同的地質(zhì)環(huán)境下,凍土活動層厚度及地表的排泄量(地表邊界流出的水量)也有所不同。隨著地層滲透率的增大,地表的排泄量也逐步增大(圖8(b)),從而需要獲取更多的補(bǔ)給量才能維持地表水位；當(dāng)?shù)貙訚B透率較小時(1×10-14m2和1×10-40m2),地表的排泄量接近于0 kg/s，此時,很小的補(bǔ)給量便能維持地表水位。

地表排泄量會直接影響地表生態(tài)環(huán)境的發(fā)展。地質(zhì)環(huán)境對生態(tài)環(huán)境的影響主要表現(xiàn)在以下2個方面:① 滲透率較大時,地表排泄量大,用于維持地表水位的補(bǔ)給量也越大。當(dāng)?shù)乇慝@得的補(bǔ)給量小于該值時,地下水位將不斷下降，植被生長與凍土水熱過程的平衡受到干擾,本就十分脆弱的植被生態(tài)系統(tǒng)被破壞,不僅使地表侵蝕容忍量驟降,而且生態(tài)環(huán)境將逐漸向荒漠化的趨勢發(fā)展。② 滲透率較小時,地表排泄量較小,用于維持地表水位的補(bǔ)給量也較小。當(dāng)降水入滲補(bǔ)給量大于該值時,地下完全飽和,處于蓄滿產(chǎn)流狀態(tài)并形成地表徑流,生態(tài)環(huán)境將向著沼澤濕地方向發(fā)展。對于青藏高原而言,地表朝著沼澤濕地的趨勢發(fā)展,有利于發(fā)揮高原濕地的生態(tài)蓄水、水源補(bǔ)給和氣候調(diào)節(jié)等重要生態(tài)功能,使生態(tài)系統(tǒng)良性循環(huán)。

3.3 氣候變暖條件下熱喀斯特湖與多年凍土相互作用

在氣候變暖的影響下,多年凍土不斷退化,對凍土水文過程及高原生態(tài)環(huán)境產(chǎn)生了重要影響[29]。為探究氣候變暖背景下熱喀斯特湖與多年凍土的相互作用,在前述模型基礎(chǔ)上構(gòu)建了湖和地表升溫方案。升溫方案中,湖和地表邊界溫度在多年平均溫度(圖4)的基礎(chǔ)上,以0.03 ℃/a升溫直到167 a后保持不變,溫度整體升高5 ℃。由圖9(a)和圖9(b)可以看出,第100 a,模型范圍內(nèi)存在高含冰量凍土；而到第400 a,僅模型中部殘留了小范圍的低含冰量(冰飽和度小于0.06)凍土,表明地表和湖溫度的升高加速了多年凍土的退化。對比第100 a和第400 a計算結(jié)果可知,凍土退化使水流通道得到恢復(fù),地下水流速增大；同時,地下水流向也發(fā)生了變化。模擬中后期,深部凍土完全退化,地表僅在冬季發(fā)生凍結(jié),此時形成了季節(jié)凍土(圖9(c))。模型第50～70 a間,在地表升溫的作用下,活動層厚度呈現(xiàn)逐漸增加的變化趨勢,由3.81 m增加至4.50 m(圖9(c))。由于受到湖的側(cè)向熱侵蝕,活動層之下的凍土逐步退化,融區(qū)上界面不斷上移；最終,深部凍土全部退化,地表形成季節(jié)凍土。氣候變暖加速凍土的退化,將加速釋放儲存在凍土中的碳,從而加劇全球氣候變暖,進(jìn)一步加速凍土的退化。此外,在氣候變暖的背景下,凍土退化可能導(dǎo)致潛水面下降,青藏高原生態(tài)環(huán)境更容易向荒漠化方向發(fā)展。

圖9 升溫方案模擬結(jié)果Fig.9 Simulated results of the temperature rise scheme

4 結(jié) 論

本文以青藏高原北麓河盆地典型熱喀斯特湖區(qū)域為研究對象,建立了考慮熱傳導(dǎo)和熱對流過程的“水- 冰- 熱”耦合數(shù)值模型,并對其進(jìn)行了驗證。模擬了熱喀斯特湖與多年凍土的相互作用關(guān)系及熱喀斯特湖的排泄量,并對不同地質(zhì)和氣候變暖條件下的熱喀斯特湖發(fā)育方向及生態(tài)效應(yīng)進(jìn)行了探討。主要結(jié)論如下:

(1) 熱喀斯特湖持續(xù)向多年凍土傳遞熱量,導(dǎo)致多年凍土逐步退化。貫穿融區(qū)形成前后,熱傳導(dǎo)和熱對流在多年凍土退化過程中分別占主導(dǎo)作用。在地表溫度作用下,形成的活動層厚度為3.35 m。

(2) 貫穿融區(qū)形成后,凍土層下水、凍土層上水以及湖水間的水力聯(lián)系增強(qiáng),形成新的地下水循環(huán)系統(tǒng)。整個模擬時段內(nèi),熱喀斯特湖表現(xiàn)為負(fù)均衡,其排泄量在285～388 a間顯著增加。

(3) 地層滲透率大時,地表和湖的排泄量大,熱喀斯特湖易于疏干,生態(tài)環(huán)境將逐漸向荒漠化方向發(fā)展；反之,地表和湖的排泄量小,熱喀斯特湖呈現(xiàn)擴(kuò)張趨勢,生態(tài)環(huán)境向沼澤濕地方向發(fā)展。

(4) 氣候變暖加速了凍土退化和熱喀斯特湖的發(fā)展,多年凍土逐漸轉(zhuǎn)變?yōu)榧竟?jié)凍土,嚴(yán)重影響寒區(qū)生態(tài)水文過程。