黃熠麗,崔婧嫄,陳思緒,張海平

(同濟大學(xué)環(huán)境科學(xué)與工程學(xué)院,上海 200092)

pH是水溶液的一個重要特性參數(shù),是影響水生生物繁殖,反映水質(zhì)狀況的一個綜合指標(biāo)。pH過低可能引發(fā)水體酸化等環(huán)境問題,過高則可能伴隨水體富營養(yǎng)化等風(fēng)險[1-2]。在飲用水處理和輸送過程中,pH被認(rèn)為是水處理化學(xué)中最重要的指標(biāo)之一[3],其值高低和波動程度對管網(wǎng)輸送、原水預(yù)處理、混凝沉淀、飲用水消毒等過程都有很大的影響。因此,建立湖庫pH機理性模型,有助于深刻認(rèn)識和描述pH與水質(zhì)因子間的關(guān)系,探討湖庫pH的變化特征與影響因素;也有助于預(yù)測湖庫pH變化趨勢,指導(dǎo)后續(xù)水處理工藝中混凝劑的投加量,以節(jié)約經(jīng)濟成本;在湖庫實測點位有限的情況下,還有助于進一步了解pH的空間分布,對于湖庫特別是飲用水水源型水庫的管理有著重要意義。

然而,由于pH的非保守性質(zhì)[4],通過模擬水域生物、化學(xué)及物理過程以模擬pH的變化并不容易,目前應(yīng)用于自然水體pH值建模與預(yù)測的成熟模型還相對較少。已有的模型可歸納為經(jīng)驗?zāi)Ｐ蚚5-7]、pH控制模型[8-10]、生物地球化學(xué)模型[11-14]3類。經(jīng)驗?zāi)Ｐ投嗷趐H與其他水質(zhì)指標(biāo)間的定性與定量關(guān)系建立經(jīng)驗公式,未能揭示pH的影響機理;pH控制模型則多用于養(yǎng)殖體系和實驗室培養(yǎng)等人為控制下的水體,以監(jiān)測、調(diào)節(jié)并預(yù)測水體pH,此類模型較少考慮水動力的對流擴散影響,無法模擬較大流域中pH的時空異質(zhì)性;生物地球化學(xué)模型通常具有更復(fù)雜的模擬結(jié)構(gòu)和較高的模擬精度[15],但此類模型通常需要大量的參數(shù),在實際應(yīng)用中,部分參數(shù)(如各類菌落數(shù)、磷酸鹽濃度、鋁離子濃度等)難以獲得,使得其在一些場合不能適用,且目前的研究主要集中于海洋、酸水排放區(qū)等水域,其生態(tài)系統(tǒng)與湖庫有較大區(qū)別。

在湖庫條件下,pH的變化與藻類活動呈現(xiàn)出了顯著的相關(guān)性,被認(rèn)為是葉綠素a(Chl-a)濃度的被動因子[16]。但已有的研究[17]中,較少有進一步研究湖庫凈生產(chǎn)力并量化其對水體pH的作用,在模擬過程中或是簡化了浮游植物的生化影響,或是通過復(fù)雜的模板模擬藻類生長、死亡、排泄等過程,而藻類實際的生物活動更為復(fù)雜,參數(shù)精度受限,從而影響pH模擬的結(jié)果。鑒于以上背景,本研究基于藻類活動對湖庫pH的顯著影響,利用MIKE 21建立水動力模型,以長期實測Chl-a濃度數(shù)據(jù)為基礎(chǔ),旨在提出一套實用型的湖庫pH模擬方法,并具備以下特征:(1)模型參數(shù)易在湖庫的日常管理中獲得;(2)模型方法能夠反映湖庫存在的生物活動、水化學(xué)平衡及流體動力學(xué)過程;(3)模擬結(jié)果能夠體現(xiàn)藻類活動和湖庫流態(tài)對pH時空波動的影響。最后,將模擬方法置于華東地區(qū)兩座重要飲用水水源型水庫中進行驗證,可為湖庫的pH預(yù)警預(yù)報提供技術(shù)支撐,促進湖庫的科學(xué)管理和供水安全保障。

1 模型方法與構(gòu)建

1.1 湖庫碳循環(huán)

研究[5]發(fā)現(xiàn),湖庫水體pH的動態(tài)平衡主要受游離CO2及碳酸鹽的平衡系統(tǒng)影響,即碳循環(huán)是影響淡水湖庫pH的主要框架,水體溶解無機碳(DIC)濃度與總堿度(TA)共同決定了水體pH的大小,這一過程如圖1所示,涉及到浮游植物、浮游動物、沉積物等對象,它們通過不同的生物化學(xué)作用影響水中碳酸鹽的含量。浮游植物在水體中的生命活動主要由藻類生長、新陳代謝,以及被浮游動物攝食、沉降和死亡這幾部分組成;浮游動物的生命過程則主要包括攝食藻類、呼吸、排泄及死亡等過程;湖泊和底泥中的碎屑有機物、糞便排泄物在生物及微生物作用下會發(fā)生礦化降解[18];而在水氣表面,氣體受水體與大氣的分壓差溶解或逸出。這些過程都伴隨著CO2的產(chǎn)生與消耗,引起水體DIC濃度的變化。

圖1 湖庫無機碳循環(huán)基本框架及其對pH值的影響Fig.1 Basic Framework of Inorganic Carbon Cycle in Lakes and/or Reservoirs and the Effect on pH Values

1.2 藻類活動對pH的影響

如前所述,CO2含量直接與水中浮游生物特別是水生植物的數(shù)量、生命活動的旺盛程度有關(guān),而Chl-a濃度可用于表征水中藻類生物量[22]。因此,湖庫pH的動態(tài)變化理論上和水體Chl-a濃度具備相關(guān)性。對華東地區(qū)某水庫2020年的pH與Chl-a濃度建立回歸關(guān)系,如圖2所示?？傮w而言,pH隨Chl-a濃度的增加而升高。

圖2 華東地區(qū)某水庫輸水口pH值與Chl-a濃度相關(guān)性Fig.2 Correlation between pH Values and Chl-a Concentration at a Reservoir Outlet in East China

同時,在任一Chl-a濃度下,pH仍有較大幅度波動,表明pH同時還受其他因子影響,例如受水溫和光照強度的規(guī)律性影響及水生植物初級生產(chǎn)的綜合性影響,pH表現(xiàn)出明顯的晝夜變化、季節(jié)循環(huán)以及與許多其他因素相關(guān)的強烈不規(guī)則波動[23-25],總體呈現(xiàn)為白天升高、夜間降低,春夏季高、秋冬季低的趨勢[26-28]。水溫主要從CO2溶解度對pH產(chǎn)生影響,隨水溫升高,CO2溶解度大幅降低,逸出量減少,pH變低,且溫度升高同樣提高了藻類的生命活性;光照強度則主要通過藻類光合作用對pH產(chǎn)生影響,白天藻類光合作用占據(jù)優(yōu)勢,水體pH升高,到傍晚達(dá)到最高值,夜間水生生物呼吸產(chǎn)生CO2,pH逐漸降低,到次日早晨達(dá)到最低值。因此,本文通過模擬藻類在不同光照、溫度等情景下由于初級凈生產(chǎn)力產(chǎn)生的CO2,同時考慮水-氣交換、有機物、浮游動物等作用,模擬湖庫pH變化過程。

1.3 CO2的模擬與計算

CO2的模擬是構(gòu)建湖庫pH模型的基礎(chǔ),本文應(yīng)用MIKE 21軟件中的ECOLab模板平臺,編制CO2計算程序,并與水動力模型耦合。

CO2的模擬分為對流擴散、水-氣交換和生化反應(yīng)3個方面。對流擴散作用以水動力模型為基礎(chǔ),前期已通過對氯化物濃度的率定獲得擴散系數(shù)取值;CO2在水-氣界面的轉(zhuǎn)移量采用大氣碳通量(Cflux)進行計算;生化反應(yīng)則主要考慮了藻類呼吸作用(CChla_respt)、藻類光合作用(CChla_phtsyn)、浮游動物呼吸作用(Czoop_respt)、有機物降解(CCODMn-CO2)4個過程,上述過程可表達(dá)為式(1)。

(1)

其中:CCO2——CO2質(zhì)量濃度,mg/L;

τ——單位時間,d;

Cflux——CO2傳輸通量,mg/(L·d);

CChla_respt——藻類呼吸作用產(chǎn)生的CO2,mg/(L·d);

CChla_phtsyn——藻類光合作用消耗的CO2,mg/(L·d);

Czoop_respt——浮游動物呼吸作用產(chǎn)生的CO2,mg/(L·d);

CCODMn-CO2——有機物降解產(chǎn)生的CO2,mg/(L·d)。

1.3.1 水-氣交換

基于1992年Wanninkhof[29]理論和2014年更新的Wanninkhof[30]理論,在水-大氣接觸面的CO2傳輸通量采用式(2)計算,其中0.010 56為單位轉(zhuǎn)換因子。

Cflux=0.010 56kCO2K0(pCO2Air-pCO2Water)/d

(2)

其中:kCO2——CO2轉(zhuǎn)移速率,cm/h;

K0——CO2的溶解度,mol/(L·atm),1 atm=101 325 Pa;

pCO2Air——大氣的CO2達(dá)到平衡狀態(tài)時的分壓,設(shè)為恒定值380 μatm[31];

pCO2Water——水表面的CO2達(dá)到平衡狀態(tài)時的分壓,μatm,在給定鹽度、溫度、壓力等環(huán)境條件下,可通過輸入DIC和pH,在CO2SYS-Program程序中求解得到(返回單位為Pa);

d——網(wǎng)格水深,m。

在淡水條件下,CO2轉(zhuǎn)移速率用式(3)～式(4)計算,對于風(fēng)速為3～15 m/s有效,適用溫度為-2～40 ℃,總體上不確定性約為20%[30]。

(3)

Sc=1 923.6-125.06t+4.377 3t2-
0.085 681t3+0.000 702 8t4

(4)

其中:u10——水面以上10 m處的風(fēng)速,m/s;

Sc——施密特系數(shù);

600——淡水在20 ℃下的施密特系數(shù);

t——水溫,℃。

在淡水條件下,CO2的溶解度K0通過式(5)計算。

(5)

其中:T——水溫,K。

1.3.2 藻類呼吸作用

藻類呼吸作用產(chǎn)生的CO2通過式(6)進行計算。

(6)

其中:rresp——藻類呼吸速率,d-1;

tresp——呼吸作用的溫度系數(shù);

CChla——Chl-a質(zhì)量濃度,mg/L,作為實測值輸入到模型中;

CDO——DO質(zhì)量濃度,mg/L,作為實測值輸入到模型中;

Dresp——DO半飽和常數(shù),mg/L;

pC/pChla——C/Chl-a,取33;

pCO2/pC——CO2/C,為3.67。

研究[32]表明,水溫對控制藻類呼吸作用具有重要意義,被認(rèn)為是控制季節(jié)性生長的主要因素。水溫對于呼吸作用的影響通過引入溫度校正系數(shù)實現(xiàn),根據(jù)參考溫度為20 ℃的Arrhenius表達(dá)式返回溫度相關(guān)性。DO作為藻類新陳代謝的限制因子之一,可通過米氏方程量化其限制作用。

1.3.3 藻類光合作用

藻類光合作用消耗的CO2通過式(7)～式(8)進行計算。

(7)

(8)

其中:mave——日均藻類生長系數(shù),d-1,通過式(8)進行計算;

tphtsyn——光合作用的溫度系數(shù);

Cphtsyn——CO2半飽和常數(shù),mg/L;

mday——白天的藻類生長系數(shù),d-1;

Tsunrise——當(dāng)日日出時刻,h;

Tsunset——當(dāng)日日落時刻,h。

光照對藻類光合作用有決定性的影響。本模型對光照強度做了簡化處理,認(rèn)為白天的光照強度為一恒定值,夜間為0。水溫對于光合作用的影響同樣通過溫度校正系數(shù)實現(xiàn),并通過基于CO2半飽和常數(shù)的米氏方程計算水體CO2濃度對于藻類光合作用的限制作用。

1.3.4 浮游動物呼吸作用

浮游動物呼吸作用產(chǎn)生的CO2通過式(9)進行計算。

Czoop_respt=rzoop×tzoopt-20×z

(9)

其中:rzoop——浮游動物呼吸作用產(chǎn)生CO2的速率,d-1;

tzoop——浮游動物呼吸作用的溫度系數(shù);

z——浮游動物與細(xì)菌生物量,mg/L,設(shè)為常數(shù)。

1.3.5 有機物降解

有機物降解產(chǎn)生的CO2通過式(10)進行計算。

(10)

其中:KCODMn——CODMn降解系數(shù),d-1;

tCODMn——CODMn降解的溫度系數(shù);

CCODMn——CODMn質(zhì)量濃度,mg/L,作為實測值輸入到模型中,(CCODMn+1)表示將實測CODMn數(shù)據(jù)轉(zhuǎn)化為TOC/(mg·L-1)的經(jīng)驗值;

DCODMn——CODMn降解的半速率反應(yīng)常數(shù),mg/L。

1.4 DIC與pH的模擬與計算

在平衡狀態(tài)中,CO2進入水中的反應(yīng)可通過式(11)～式(14)進行描述。

(11)

(12)

(13)

(14)

其中:K1——碳酸的一級解離常數(shù),在淡水條件下受溫度影響;

K2——碳酸的二級解離常數(shù),在淡水條件下受溫度影響;

[H+]——H+摩爾濃度,mol/L;

NDIC——DIC摩爾濃度,mol/L。

(15)

在淡水條件下,不考慮磷堿度及硅堿度,淡水堿度可通過式(16)～式(17)。

NTA=AC+[OH-]-[H+]

(16)

(17)

其中:NTA——TA摩爾濃度,mol/L;

AC——碳堿度的摩爾濃度,mol/L;

[OH-]——OH-摩爾濃度,mol/L。

結(jié)合上述公式,可以得到TA、DIC、[H+]與之間存在的關(guān)系式,如式(18)。

(18)

其中:Kw——水的離子積常數(shù),在淡水條件下受溫度影響。

2 模型應(yīng)用

2.1 模擬水體

選取A、B兩座華東地區(qū)飲用水源水庫對本文所建立的pH模型進行驗證。兩座水庫的水文水質(zhì)特征有較大差別。A水庫總面積達(dá)66.26 km2,設(shè)計有效庫容為4.35億m3,供水規(guī)模逾800萬m3/d,庫區(qū)水力停留時間較長,約為20 d,水流緩慢,局部存在藻類過度繁殖和聚集等問題,庫內(nèi)Chl-a質(zhì)量濃度在2～41 μg/L,輸水區(qū)pH值則明顯大于上游來水,為8.0～9.0,pH偏高且波動劇烈,對于自然水體而言變化明顯。B水庫總面積為1.35 km2,設(shè)計有效庫容為953萬m3,供水能力逾200萬m3/d,水體交換時間短,水力停留時間僅為3～5 d,整體流通性強,大型藻類、沉水植物難以生長,大部分區(qū)域浮游藻類質(zhì)量濃度低,常年為1～5 μg/L,輸水區(qū)pH值為7.8～8.4,處于相對較低的水平。

利用MIKE 21軟件分別建立兩座水庫的二維水動力水質(zhì)模型。A水庫模型地形網(wǎng)格數(shù)量共102 826個,平均每個網(wǎng)格分辨率為644 m2,B水庫模型地形網(wǎng)格數(shù)量共26 042個,平均每個網(wǎng)格分辨率為52 m2。利用模型計算每個網(wǎng)格的垂向平均水動力參數(shù)(水位、流量和流速等)和水質(zhì)參數(shù)(CO2、DIC、pH等)。模型計算時間為2020年5月1日—2021年5月1日。

利用庫區(qū)實測水位資料對模型進行率定,實測數(shù)據(jù)來自5 min/次的在線數(shù)據(jù)。A水庫平均相對誤差(MRE)為4.6%,均方根誤差(RMSE)為0.16 m;B水庫的MRE為0.73%,RMSE為0.06 m,表明所構(gòu)建的水動力模型較好地模擬了水文水動力要素的主要變化特征,為后續(xù)pH建模提供了可靠基礎(chǔ)。

圖3 模擬與實測pH值比較(2020年5月—2021年5月)Fig.3 Comparison of Simulated and Measured pH Values (from May 2020 to May 2021)

2.2 模型條件設(shè)置與參數(shù)取值

根據(jù)水庫上游庫外的實測pH、TA、水溫等數(shù)據(jù),利用CO2SYS-Program程序計算上游庫外的CO2濃度作為模型邊界條件。庫區(qū)風(fēng)速、水深、水溫、Chl-a濃度、CODMn、DO濃度、大氣CO2分壓、鹽度和TA是模型輸入條件。其中,水深數(shù)據(jù)從水動力模型讀取;風(fēng)速、水溫、Chl-a濃度、CODMn、DO濃度為實測值;大氣CO2分壓設(shè)為恒定380 μatm。通過將實測日堿度輸入和平均值輸入進行測試比較發(fā)現(xiàn),模型對這兩座水庫的堿度輸入不敏感,結(jié)果相差較小,為簡化起見,采用全年平均值輸入到模型中參與計算,但這點是否適用于其他水庫,還需進一步觀察。

參數(shù)取值的合理性直接影響到模型對水體系統(tǒng)表征的可靠性與合理性[33-34],進而影響到模型精度。在本研究中,通過文獻(xiàn)查閱及過往運行經(jīng)驗得出參數(shù)取值范圍后,用實測數(shù)據(jù)對模型進行率定,得到的參數(shù)值如表1所示,與其他淡水湖庫的研究結(jié)果總體一致[35-36]。A水庫藻類活動強烈,實測得到的CODMn包含了藻類生長而產(chǎn)生的有機物,因此,在表觀上幾乎未發(fā)生降解,即綜合降解系數(shù)為0。

表1 模型參數(shù)設(shè)置Tab.1 Setup of Model Parameter

2.3 模型結(jié)果

將模擬得到的輸水口pH與實測數(shù)據(jù)進行比較,如圖3所示。實測數(shù)據(jù)的監(jiān)測頻率為每天一次。A水庫的MRE為1.02%,平均絕對誤差(MAE)為0.086,RMSE為0.111;B水庫的MRE為0.77%,MAE為0.062,RMSE為0.082。

本模型考慮了對流擴散、水-氣交換、藻類活動等過程,模擬的總體結(jié)果顯示,模型對于兩座水庫pH的高低及全年變化趨勢都能有較好地重現(xiàn)與反演,且模擬趨勢與實測趨勢具有良好一致性,表明模型的設(shè)置總體合理。

2.4 模型驗證

基于表1中的率定參數(shù),選取2022年4月1日—6月30日的數(shù)據(jù)對模型進行驗證,結(jié)果如圖4所示。A水庫的MRE為1.24%,MAE為0.106,RMSE為0.133;B水庫的MRE為0.97%,MAE為0.081,RMSE為0.106,驗證精度與率定模型相當(dāng)。

圖4 模擬與實測pH值比較(2022年4月—2022年6月)Fig.4 Comparison of Simulated and Measured pH Values (from April 2022 to June 2022)

3 結(jié)果分析與討論

3.1 結(jié)果分析

從模擬結(jié)果來看,本文所構(gòu)建的模型能夠較好地重現(xiàn)pH在自然水體中的變化趨勢,模型模擬結(jié)果與實測數(shù)據(jù)的MAE和RMSE均較低,表現(xiàn)出良好的模擬性能。所選擇的模擬對象水體A、B水庫的水文水質(zhì)情況差異較大,而該模型均能較好地對其pH進行模擬,呈現(xiàn)出一定的普適性。

A水庫的實測pH波動相對較大,模型對于pH起峰回落的過程及其所對應(yīng)的時間響應(yīng)良好,特別是pH谷值的模擬精度較高,在大多數(shù)情況下均較接近于實測值;pH峰值受藻類活性、氣候條件、光照強度等環(huán)境因子的疊加影響,實際驅(qū)動作用可能比模型更為強烈,模擬值往往相較于實測值偏低,這與Cerco等[37]在富營養(yǎng)淺水海灣的模擬結(jié)果相符。

模型對相關(guān)數(shù)據(jù)獲取困難的一些生化過程進行了簡化描述,如將水體浮游動物量設(shè)為定值、不考慮優(yōu)勢藻種的季節(jié)性變化、底部沉積物溶解等,Liang等[17]的研究同樣表明上述因素引起的pH變化并不明顯,甚至可能增加模擬難度并降低精度。盡管在兩座水庫的模擬中未發(fā)現(xiàn)由此產(chǎn)生明顯誤差,將來在獲得更多觀測資料后,可進一步檢驗?zāi)Ｐ头椒?優(yōu)化模型結(jié)構(gòu)。

3.2 參數(shù)敏感性分析

本文所構(gòu)建的模型考慮了Chl-a濃度、DO濃度、水溫、CODMn、堿度、風(fēng)速、水深、光照、CO2濃度、pCO2等實測值或計算值,利用局部敏感性分析中的擾動法,量化所選參數(shù)對模擬結(jié)果的影響,其計算過程如式(19)。

(19)

其中:SS——參數(shù)敏感性數(shù)值;

ΔS——模擬結(jié)果的變化量,文中指ΔpH;

S——模擬結(jié)果,文中指pH;

ΔE——參數(shù)變化量,變化幅度設(shè)為合理范圍內(nèi)的±20%;

E——參數(shù)值。

計算得兩座水庫的敏感性數(shù)值及其絕對值排序如表2所示。

表2 參數(shù)敏感性及排序Tab.2 Parameter Sensitivity and Ranking

篩選出敏感性數(shù)值大于1%的參數(shù),在A水庫中表現(xiàn)為mday=光照時長>tresp>tphtsyn>tzoop>rzoop=zoop>Chl-a濃度>rresp,在B水庫中表現(xiàn)為CO2濃度>tphtsyn>風(fēng)速>水深>水溫?？傮w而言,模型對Chl-a濃度、水溫較為敏感,且部分因子呈現(xiàn)出異參同效的特征。Chl-a濃度及其相關(guān)參數(shù)(mday、rresp等)是影響湖庫特別是富營養(yǎng)化湖庫pH變化的重要驅(qū)動因子,其敏感性數(shù)值受到藻類指標(biāo)濃度的高度影響[38];水溫相關(guān)參數(shù)(tphtsyn、tresp等)敏感性也較高,且隨溫度季節(jié)性變化呈現(xiàn)出時間變異性特征[39];CO2作為邊界條件,對B水庫等面積較小、水動力較強的湖庫影響較大。

3.3 結(jié)論

(1)本研究基于藻類活動對于湖庫pH的顯著影響,利用Chl-a濃度、水溫、光照、風(fēng)速、DO濃度、CODMn、CO2濃度、pCO2、TA等實測值或計算值,通過耦合水動力模型,考慮對流擴散作用、水-氣交換、藻類呼吸作用、藻類光合作用、浮游動物呼吸作用、有機物降解等生態(tài)動力學(xué)過程,建立了適用于湖庫的pH模型。

(2)所構(gòu)建的模型應(yīng)用于華東地區(qū)兩座飲用水源水庫。率定結(jié)果表明,pH模擬的MRE分別為1.02%和0.77%,MAE分別為0.086和0.062,RMSE分別為0.111和0.082,且模擬趨勢與實測趨勢具有良好一致性。模型驗證精度基本相當(dāng),模擬精度較高,可較好反映湖庫的pH動態(tài)變化過程。

(3)所模擬的兩座水庫存在較為不同的水文水質(zhì)特征,而本研究所建立的模型均能較好地對其pH進行模擬,初步表明該模型具有華東地區(qū)湖庫pH值建模的普適性。模型參數(shù)對于模型結(jié)果的影響較大,且參數(shù)敏感性隨湖庫水文水質(zhì)特征表現(xiàn)出差異。同時,模型所需輸入數(shù)據(jù)基本可從日常監(jiān)測資料獲取,模型的實用性較強。