唐繼張， 周維博， 安寶軍， 楊 浩， 李 穎

(1.長(zhǎng)安大學(xué) 水利與環(huán)境學(xué)院， 陜西 西安 710064； 2.長(zhǎng)安大學(xué) 旱區(qū)地下水文與生態(tài)效應(yīng)教育部重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710064； 3.陜西省西咸新區(qū)灃東新城斗門(mén)水庫(kù)建設(shè)管理中心， 陜西 西安 710086)

1 研究背景

城市化進(jìn)程的快速推進(jìn)使得現(xiàn)代社會(huì)對(duì)自然的改造力度越來(lái)越大，城市建設(shè)發(fā)展過(guò)程中對(duì)水資源的要求也在不斷提高[1-4]。近年來(lái)，人工湖泊在城市的規(guī)劃定位中突顯了其重要價(jià)值[5-8]：不僅具有旅游價(jià)值，推動(dòng)當(dāng)?shù)丶爸苓叺纳虡I(yè)發(fā)展、提高居民的生活質(zhì)量，而且一些湖泊還兼具調(diào)洪蓄水、保障城市供水安全、改善局地小氣候等生態(tài)價(jià)值。但與天然湖泊相比，人工湖泊較為封閉，水環(huán)境單一，且水體流動(dòng)性較差，納污及自?xún)裟芰Σ蛔鉡9-10]，加之周邊人類(lèi)活動(dòng)的影響，容易引起湖泊富營(yíng)養(yǎng)化、水質(zhì)惡化等問(wèn)題。若不及時(shí)有效地進(jìn)行預(yù)防治理，會(huì)對(duì)城市的建設(shè)發(fā)展起到阻礙作用。明確人工湖泊的水流場(chǎng)特征及水體交換能力，對(duì)預(yù)防湖泊水質(zhì)惡化與污染治理、促進(jìn)城市水環(huán)境健康發(fā)展意義重大[11-13]。

許多學(xué)者針對(duì)湖泊水體交換能力問(wèn)題進(jìn)行了大量研究。例如，許莉萍等[14]運(yùn)用數(shù)值方法建立水動(dòng)力和對(duì)流擴(kuò)散模型，研究了人工湖的水體更新時(shí)間和水體交換率。Delhezé等[15]建立了一維河流和二維河口的耦合模型對(duì)斯凱爾特河河口對(duì)流擴(kuò)散進(jìn)行研究，得到了河口水體滯留、輸運(yùn)及影響時(shí)間。王冼民等[16]利用EcoTaihu模型，在不同蓄水量方案下研究太湖的氮、磷濃度最低值的換水周期。李云良等[17]運(yùn)用數(shù)值模擬方法，定量對(duì)比了不同季節(jié)時(shí)鄱陽(yáng)湖的換水周期，得出了鄱陽(yáng)湖在四季的換水能力大小。Umgiesser等[18]建立了地中海瀉湖三維數(shù)學(xué)模型，得到了潮流流場(chǎng)、水體交換率、水體輸運(yùn)時(shí)間和水體混合等情況。李大鳴等[19]基于物質(zhì)遷移和對(duì)流擴(kuò)散理論，研究了洋河水庫(kù)中TN、TP、NH3—N質(zhì)量濃度的變化和分布。姜恒志等[20]針對(duì)普蘭店灣海域在潮流作用下的水交換能力，將海灣分為不同區(qū)域建立潮流數(shù)值模型，得到各區(qū)域的水體交換率。鑒于數(shù)值模擬方法在湖泊研究中的成熟運(yùn)用[21-22]，本文采用MIKE軟件，在水動(dòng)力基礎(chǔ)上耦合染色劑模型，從時(shí)間及空間尺度上，對(duì)昆明池(試驗(yàn)段)進(jìn)行模擬研究，結(jié)合湖區(qū)流場(chǎng)特征，定量分析整個(gè)湖區(qū)的換水能力，探究湖區(qū)在正常運(yùn)行條件下最可能出現(xiàn)水質(zhì)惡化問(wèn)題的水域及原因。研究結(jié)果為昆明池(試驗(yàn)段)水體的污染預(yù)防與治理提供了科學(xué)指導(dǎo)，也為后期昆明池?cái)U(kuò)建工程提供了理論依據(jù)，還對(duì)于一般人工湖泊的水環(huán)境綜合管理及調(diào)度運(yùn)行具有重要理論意義和參考價(jià)值。

2 資料來(lái)源與研究方法

2.1 研究區(qū)概況

西安昆明池(另稱(chēng)斗門(mén)水庫(kù))位于東經(jīng)108°45′～108°48′，北緯34°11′～34°13′之間。地處關(guān)中平原中部，秦嶺以北，屬暖溫帶半濕潤(rùn)大陸性季風(fēng)氣候區(qū)，其利用原昆明池遺址進(jìn)行擴(kuò)建，分為南湖和北湖，是一座以調(diào)蓄“引漢濟(jì)渭”工程來(lái)水、兼顧城市生活供水、灃河防洪及改善生態(tài)環(huán)境的綜合性平原水庫(kù)工程。

本次研究對(duì)象為昆明池試驗(yàn)段工程(以下簡(jiǎn)稱(chēng)為昆明池)，湖區(qū)設(shè)有5座人工島嶼，水域面積0.47 km2，庫(kù)底高程396.8 m，正常蓄水位400.47 m，蓄水深度3.67 m，蓄水量為155×104m3。湖區(qū)位于灃河與太平河之間，由灃惠渠引水，再退水至太平河，水系較為豐富。但灃河及灃惠渠水質(zhì)較差，且圍繞昆明池打造的休閑公園和周邊建筑群以及一系列的人為活動(dòng)使得湖區(qū)現(xiàn)狀水質(zhì)較差，局部水域出現(xiàn)富營(yíng)養(yǎng)化狀態(tài)。研究概況見(jiàn)圖1。

圖1 研究區(qū)概況

2.2 數(shù)據(jù)來(lái)源

本文所用昆明池地形及水位數(shù)據(jù)來(lái)自《陜西省斗門(mén)水庫(kù)工程環(huán)境影響報(bào)告書(shū)》及《西咸新區(qū)灃東新城斗門(mén)調(diào)節(jié)水庫(kù)巖土工程勘察報(bào)告》；降雨蒸發(fā)、風(fēng)速風(fēng)向資料來(lái)自于課題組于研究區(qū)設(shè)置的小型氣象站自動(dòng)監(jiān)測(cè)采集。

2.3 研究方法

2.3.1 MIKE21 FM水動(dòng)力模型

(1)模型基本方程。MIKE21 FM為平面二維水流數(shù)學(xué)模型，該模型遵循Navier-Stokes方程，且服從Boussinesq假定和流體靜定假設(shè)。其二維水流連續(xù)方程為：

(1)

二維水流動(dòng)量方程為：

(2)

(3)

(4)

(5)

考慮流體黏性影響，水流的渦黏項(xiàng)用Tij表示，其值按下式估算：

(6)

(7)

(8)

式中：A為黏性系數(shù)，m2/s。

(2)方程組離散求解。模型采用有限體積法進(jìn)行研究區(qū)域的離散求解，離散為若干個(gè)不規(guī)則的三角形網(wǎng)格，可采用三角形與四邊形相結(jié)合的混合網(wǎng)格對(duì)局部區(qū)域進(jìn)行加密，以便更好地反映復(fù)雜邊界，有限體積法保證了水量和動(dòng)量在計(jì)算域內(nèi)的守恒。本次對(duì)昆明池地形離散共計(jì)生成4 178個(gè)網(wǎng)格，2 307個(gè)節(jié)點(diǎn)，最終剖分網(wǎng)格見(jiàn)圖2(圖中已標(biāo)注進(jìn)水口、出水口及龍頭監(jiān)測(cè)點(diǎn)位置)。

(3)模型模塊的選擇。水動(dòng)力模塊(HD)是模型模擬的基礎(chǔ)，主要通過(guò)對(duì)研究區(qū)域的地形網(wǎng)格化處理、控制網(wǎng)格干濕條件、選取合適的求解方法并考慮模型的輸入輸出作用，從而計(jì)算出模擬的研究區(qū)水深和流場(chǎng)分布規(guī)律。對(duì)流擴(kuò)散(AD)模塊通過(guò)設(shè)置不同類(lèi)型的擴(kuò)散系數(shù)模擬物質(zhì)在水體中經(jīng)過(guò)對(duì)流和擴(kuò)散過(guò)程而發(fā)生的擴(kuò)散現(xiàn)象特征，其優(yōu)點(diǎn)在于忽略了污染物質(zhì)在水中運(yùn)移時(shí)發(fā)生的復(fù)雜的物理、化學(xué)、生物及生態(tài)過(guò)程，突出了區(qū)域水體流場(chǎng)的作用?？紤]湖區(qū)的水體環(huán)境及研究目的，本次模擬選用水動(dòng)力模塊和對(duì)流擴(kuò)散模塊作為研究手段。

圖2 研究區(qū)網(wǎng)格剖分

(4)定解條件的確定。湖岸及湖心島邊界面的法相流速為零，設(shè)為零法向流速的陸地邊界；進(jìn)水口為流量，設(shè)定為實(shí)測(cè)的隨時(shí)間變化的流量邊界。出水口為水位，設(shè)定為實(shí)測(cè)的水位邊界；模型輸入項(xiàng)主要還包括湖區(qū)風(fēng)速風(fēng)向及降雨蒸發(fā)，皆設(shè)定為在空間上恒定，隨時(shí)間變化。

2.3.2 水體交換研究方法

(1)換水周期。關(guān)于湖泊水體交換能力的討論，許多學(xué)者從不同角度進(jìn)行了定量的研究，提出了水齡、水力滯留時(shí)間、更新時(shí)間、曝光時(shí)間和換水周期等概念。根據(jù)研究對(duì)象及模型特點(diǎn)，本文選用換水周期來(lái)定量研究湖區(qū)的換水能力，其中換水周期為可溶性保守物質(zhì)初始的單位濃度經(jīng)過(guò)水流對(duì)流擴(kuò)散作用，衰減到某一程度所需要的時(shí)間，表征湖泊對(duì)于湖區(qū)水體自我凈化能力的大小[17,23]?；诤^(qū)內(nèi)保守物質(zhì)的最終濃度幾乎不可能為零這一事實(shí)，本文換水周期選用濃度變化的指數(shù)衰減函數(shù)表示，具體公式如下：

Ct=C0·e-t/Tf

(9)

式中:t為時(shí)間，d;C0為保守物質(zhì)初始濃度值,g/L;Ct為保守物質(zhì)在t時(shí)刻的濃度值,g/L。

由公式(9)知，當(dāng)t=Tf(V/Q)時(shí)，濃度衰減至初始濃度的e-1，即將剩余濃度降低至初始濃度的37%時(shí)所需要的時(shí)間定義為換水周期。

(2)頻率分布曲線。頻率分布曲線[24]能夠直觀反映整個(gè)湖泊的換水周期面積分布，基于頻率曲線的定義，通過(guò)分析變異性，進(jìn)一步探討湖泊的換水能力，評(píng)價(jià)換水周期的優(yōu)劣。曲線越陡峭，說(shuō)明湖泊換水能力變異性越弱，系統(tǒng)換水越快；曲線越平緩，則說(shuō)明湖泊換水能力變異性越強(qiáng)，系統(tǒng)換水越慢。

2.3.3 模擬方案及模型設(shè)置 本文選定MIKE21 FM的水動(dòng)力和對(duì)流擴(kuò)散模塊，在湖泊水動(dòng)力模型基礎(chǔ)上通過(guò)耦合染色劑模型，模擬計(jì)算出昆明池的流場(chǎng)特征和換水周期分布。根據(jù)上文方法定義，本文將昆明池水體設(shè)置為單位濃度為1的保守型示蹤劑作為初始濃度場(chǎng)，湖區(qū)進(jìn)水口設(shè)為濃度為0的流量邊界，湖口開(kāi)邊界設(shè)為自由出流邊界，隨著模型的運(yùn)行，通過(guò)監(jiān)測(cè)研究區(qū)內(nèi)每一個(gè)網(wǎng)格單元的剩余濃度變化，得到昆明池的換水周期空間分布。

考慮到湖區(qū)在夏季降水及蒸發(fā)量較大，且周邊公園日開(kāi)放時(shí)常較長(zhǎng)，人流較多，這一時(shí)期湖水水質(zhì)遭受外界因素的影響較大，故結(jié)合現(xiàn)有數(shù)據(jù)資料，設(shè)定本次模擬時(shí)段為夏季，即2018年6月1日至9月1日。

2.3.4 模型率定與驗(yàn)證 根據(jù)實(shí)測(cè)水位數(shù)據(jù)對(duì)模型進(jìn)行反復(fù)率定，調(diào)參得到的模型主要參數(shù)為：干濕網(wǎng)格臨界水深為0.05 m、Smagorinsky系數(shù)為0.28、曼寧糙率系數(shù)為28 m1/3/s。再利用研究區(qū)中龍頭監(jiān)測(cè)點(diǎn)7月份實(shí)測(cè)水位數(shù)據(jù)對(duì)模型進(jìn)行驗(yàn)證，驗(yàn)證結(jié)果如圖3所示，模擬過(guò)程水位與實(shí)測(cè)基本一致。

圖3 龍頭監(jiān)測(cè)點(diǎn)水位驗(yàn)證

(10)

依據(jù)公式(10)計(jì)算得到該模型的Skill值為0.94，大于標(biāo)準(zhǔn)值0.65，表明該模型擬合精度極好，可用于對(duì)昆明池的模擬計(jì)算。

3 結(jié)果與分析

3.1 水動(dòng)力場(chǎng)模擬

運(yùn)行模型對(duì)昆明池進(jìn)行水動(dòng)力模擬，得到夏季昆明池流場(chǎng)模擬結(jié)果，在模擬期內(nèi)，各時(shí)段流場(chǎng)變化較為一致，故選擇模擬期中間時(shí)段某一天(2018年7月26日)作為代表性流場(chǎng)進(jìn)行分析，結(jié)果見(jiàn)圖4。由圖4可見(jiàn)，昆明池在人工島附近均形成多個(gè)較為明顯的環(huán)流，方向呈逆時(shí)針流動(dòng)，其成因或?yàn)楹^(qū)在東側(cè)進(jìn)水口引入水流，西南側(cè)出水口排水，人工島在水流運(yùn)動(dòng)路徑上與其發(fā)生碰撞，改變了水流方向形成環(huán)流。在沿岸也出現(xiàn)許多較小環(huán)流，這可能是由于昆明池湖岸線起伏，特別是北岸在碼頭、觀水臺(tái)等曲折較大的臨水地帶，水流受其擾動(dòng)作用，形成范圍較小的環(huán)流。由于整個(gè)湖區(qū)水域在南北方向上呈“中間窄、兩邊大”的特點(diǎn)，水流在東西方向交換過(guò)程中，經(jīng)過(guò)湖中間地帶時(shí)受到地形阻礙，改變了水流方向，故在湖東西部區(qū)域分別形成較大的環(huán)流。

圖4 2018-07-26昆明池流場(chǎng)模擬結(jié)果

受環(huán)流等影響，昆明池流速也具有明顯分區(qū)域的特征，流速分布模擬結(jié)果見(jiàn)圖5。

由圖5可看出，貼近人工島的環(huán)流流速較大，約為0.04～0.06 m/s，但湖東北部人工島附近，因?yàn)樗蚱x水流主方向，故流速較低，約為0.002～0.005 m/s。沿湖岸線受地形因素形成的各小環(huán)流的近岸流速較大，約為0.04～0.07 m/s，個(gè)別區(qū)域流速達(dá)到0.09 m/s。而在湖中心區(qū)域，水流較為緩慢，流速基本低于0.004 m/s，且距離湖岸較遠(yuǎn)區(qū)域流速約小于0.002 m/s。

3.2 換水周期模擬

對(duì)昆明池進(jìn)行換水周期模擬，得到夏季昆明池各區(qū)域的換水天數(shù)。其中，龍頭、出水口處的濃度-時(shí)間曲線見(jiàn)圖6，分析圖6得出，龍頭處的換水周期為27 d，出口處換水周期為36 d，即使新水至東側(cè)進(jìn)水口引入，也分別需要約27、36 d才能對(duì)龍頭和出水口處水體進(jìn)行較完全的置換。這可能是由于龍頭位置處于地形夾角處，水流循環(huán)受阻，而出口位置距進(jìn)水口較遠(yuǎn)，水流路徑較長(zhǎng)，故更新緩慢，換水周期較長(zhǎng)。各“時(shí)間-濃度”為指數(shù)衰減形式，進(jìn)一步表明了本文所用的換水周期計(jì)算方法在本次模擬中具有良好的適用性，體現(xiàn)了模擬結(jié)果的可靠、合理性。

圖5 2018-07-26昆明池流速分布模擬結(jié)果

圖6 昆明池龍頭和出口站點(diǎn)濃度隨時(shí)間變化曲線

進(jìn)一步對(duì)整個(gè)湖區(qū)換水周期分析(見(jiàn)圖7)，得出湖區(qū)換水周期具有明顯的空間異質(zhì)性，根據(jù)其換水周期在東西方向上階梯遞增的特征，大致將湖區(qū)劃分成4個(gè)區(qū)域(Ⅰ、Ⅱ、Ⅲ、Ⅳ)。其中：靠近進(jìn)水口的區(qū)域Ⅰ，其水流路徑最短，故換水周期最小，約為1～10 d，但東北角區(qū)域換水周期有所增長(zhǎng)，約為10～17 d，或是因?yàn)樵搮^(qū)域偏離進(jìn)水口主流方向，且水流與該處的人工島發(fā)生碰撞，阻礙了其向東北角擴(kuò)散，故降低了東北角水流交換能力，造成換水周期延長(zhǎng)；位于龍頭前端的區(qū)域Ⅱ，其換水周期約為17～28 d，且與區(qū)域Ⅲ在湖岸地形收縮處形成一明顯分界線。結(jié)合湖區(qū)流場(chǎng)特征，分析可能是由于區(qū)域Ⅱ存在的較大環(huán)流造成的，且環(huán)流邊界恰位于分界線附近，故在此處，部分水流前進(jìn)方向改變被迫回流至區(qū)域Ⅱ內(nèi)，形成一明顯的換水周期分界線，且隨著水流路徑的增長(zhǎng)，換水周期較區(qū)域Ⅰ也有所增大；穿過(guò)東部湖區(qū)，水流到達(dá)西部湖區(qū)；區(qū)域Ⅲ換水周期約為28～32 d，在小環(huán)流，特別是靠近湖南岸一帶的小環(huán)流作用下，水流流動(dòng)方向改變，做循環(huán)流動(dòng)，換水周期普遍較大，并且在西側(cè)與區(qū)域Ⅳ形成一明顯界限；區(qū)域Ⅳ處于昆明池最西側(cè)，其進(jìn)水水流路徑最長(zhǎng)，普遍換水周期約為32～35 d，在人工島附近及靠近湖岸一帶，部分區(qū)域換水周期超過(guò)35 d，這同樣是受人工島與湖岸形成的環(huán)流影響所致。

圖7 昆明池?fù)Q水周期空間分布

3.3 頻率分布曲線

根據(jù)昆明池?fù)Q水周期模擬結(jié)果，做出其換水周期頻率分布曲線，見(jiàn)圖8。分析圖8可知：曲線上存在3個(gè)明顯拐點(diǎn)(分別對(duì)應(yīng)換水周期為10、17、34 d)，將曲線分成4段，呈現(xiàn)出“平緩-陡峭-平緩-陡峭”的特點(diǎn)，表明昆明池?fù)Q水能力變異性不具備明顯的強(qiáng)弱特性，即昆明池不能簡(jiǎn)單地定義為一個(gè)快速換水或慢速換水的湖泊系統(tǒng)。其中，約10%的昆明池湖區(qū)換水周期≤10 d，約40%的昆明池湖區(qū)換水周期≤17 d，約70%的昆明池?fù)Q水周期≤34 d，還有不足10%的昆明池湖區(qū)換水周期超過(guò)了35 d。

圖8 昆明池?fù)Q水能力變異性分析頻率分布曲線

整體上看，昆明池流場(chǎng)及換水周期分布均主要受湖中人工島及湖岸地形因素的影響。在水動(dòng)力交換弱和受地形影響形成環(huán)流的區(qū)域，水流更新時(shí)間較長(zhǎng)，換水周期較大，這種換水周期分布特征與水動(dòng)力場(chǎng)的密切相關(guān)性，表明了模型模擬結(jié)果的合理性。通過(guò)現(xiàn)場(chǎng)調(diào)研，昆明池在人工島附近及湖西部區(qū)域，更容易出現(xiàn)水體富營(yíng)養(yǎng)化等問(wèn)題，這也與昆明池模擬結(jié)果相印證，進(jìn)一步說(shuō)明了本次模型模擬的科學(xué)性與結(jié)果的合理性。通過(guò)對(duì)昆明池流場(chǎng)特征及換水周期分布的探究，能夠明確昆明池湖區(qū)范圍內(nèi)易受水體交換更新不利影響而發(fā)生水質(zhì)惡化問(wèn)題的區(qū)域，有利于對(duì)湖區(qū)水質(zhì)問(wèn)題作出針對(duì)性的預(yù)防治理措施，維持湖區(qū)水環(huán)境的健康可持續(xù)運(yùn)行。

4 結(jié) 論

(1)昆明池流場(chǎng)受人工島及湖岸地形影響較大。在各人工島附近，均形成多個(gè)逆時(shí)針流向的小環(huán)流，近島處流速較大，平均可達(dá)0.04～0.06 m/s；靠近湖岸，存在許多明顯小環(huán)流，流速約為0.04～0.07 m/s，最大流速可達(dá)0.09 m/s；在湖東西部均存在較大環(huán)流，流速基本低于0.004 m/s。

(2)昆明池?fù)Q水周期具有明顯空間異質(zhì)性，在東西方向呈階梯遞增：區(qū)域Ⅰ約1～17 d，區(qū)域Ⅱ約17～28 d，區(qū)域Ⅲ約28～32 d，區(qū)域Ⅳ約32～35 d，個(gè)別區(qū)域超過(guò)35 d。

(3)昆明池屬于快速換水與慢速換水并存的湖泊系統(tǒng)，約10%的水域換水周期≤10 d， 40%的水域換水周期≤17 d，70%的水域換水周期≤34 d，不足10%的水域換水周期超過(guò)35 d。

(4)加強(qiáng)對(duì)湖區(qū)西部的監(jiān)測(cè)，提高該水域的水體交換能力，是進(jìn)一步改善昆明池水質(zhì)，預(yù)防湖區(qū)出現(xiàn)富營(yíng)養(yǎng)化問(wèn)題的關(guān)鍵。