張?jiān)骆?，徐明德，原文?/p>

(1.中國(guó)輻射防護(hù)研究院，太原 030006；2.太原理工大學(xué)，太原 030024)

前 言

汾河水庫(kù)是引黃入晉工程的調(diào)蓄水庫(kù)，是山西省最大的飲用水水庫(kù)[1]。上世紀(jì)初，國(guó)內(nèi)外就已開展了泥沙規(guī)律的研究工作，但研究?jī)?nèi)容多局限于泥沙自身運(yùn)動(dòng)分析，例如泥沙的運(yùn)輸與淤積[2-3]，對(duì)于泥沙對(duì)水體水質(zhì)的影響規(guī)律研究卻很少涉及，而對(duì)于北方干旱地區(qū)水庫(kù)的相關(guān)研究更是鮮有報(bào)道。本文選取汾河水庫(kù)為研究對(duì)象，以汾河水庫(kù)水動(dòng)力[4]—泥沙轉(zhuǎn)移—富營(yíng)養(yǎng)化耦合模型[5]進(jìn)行模擬分析，研究高含量SSC對(duì)汾河水庫(kù)庫(kù)區(qū)水質(zhì)產(chǎn)生的影響，從而揭示SSC對(duì)庫(kù)區(qū)水質(zhì)的影響規(guī)律。

1 研究區(qū)域

1.1 區(qū)域概況

汾河水庫(kù)位于汾河干流上游段，是太原地區(qū)主要水源地。南北長(zhǎng)15 km，東西寬5 km，其總面積約達(dá)32 km2。汾河水庫(kù)的流域總面積達(dá)到5 268 km2，平均流量每秒可達(dá)21.9m3，水庫(kù)總流量為7.22 億m3，泥沙總庫(kù)容為3.45億m3。

1.2 水質(zhì)監(jiān)測(cè)

1.2.1 監(jiān)測(cè)點(diǎn)位：考慮汾河水庫(kù)自然特點(diǎn)，設(shè)置A、B、C、D、E、F6個(gè)監(jiān)測(cè)點(diǎn)，其中A點(diǎn)和B點(diǎn)分別位于汾河入庫(kù)處和澗河入庫(kù)處，C點(diǎn)、D點(diǎn)和E點(diǎn)分別位于庫(kù)區(qū)上段、中段和后段，F(xiàn)點(diǎn)位于壩址處，各監(jiān)測(cè)點(diǎn)的位置如圖1所示。

圖1 汾河水庫(kù)監(jiān)測(cè)點(diǎn)位Fig.1 Monitoring sites of Fenhe reservoir

1.2.2 監(jiān)測(cè)對(duì)象：表層水樣、沉積物樣品。

1.2.3 監(jiān)測(cè)因子及方法：包括SSC、Chla、TP、TN、固態(tài)氮等監(jiān)測(cè)因子，監(jiān)測(cè)方法選用《水質(zhì)懸浮物的測(cè)定重量法》(GB/T 11901-1989)中的重量法、冷凍浸提法、分光光度法、紫外分光光度法[6]、超聲萃取法[7]。

1.2.4 監(jiān)測(cè)時(shí)間：2017年6月11日～2017年9月23日。

2 汾河水庫(kù)模型的構(gòu)建

2.1 水動(dòng)力模型

本模型的理論基礎(chǔ)為靜水壓力和布辛尼克斯(Boussinesq)[8]假設(shè)，納維葉-斯托克斯(Navier-Stokes)[9]方程對(duì)水動(dòng)力方程進(jìn)行控制。二維水動(dòng)力模型如下：

(1)

(2)

(3)

其中：φ(x,y,t)為表面高程；p、q為在x、y方向上的流量分密度；h(x,y,t)為水深均值；C為固定值，謝才系數(shù)；g為重力加速度；τs(V)為風(fēng)應(yīng)力；Vx、Vy分別為風(fēng)應(yīng)力沿x、y方向的分量；Ω(x,y)為科里奧利參數(shù)；Pa(x,y,t)為大氣壓；ρw為水密度；S為源質(zhì)量；Six，Siy分別為在x、y方向上源質(zhì)量的動(dòng)量分量；τxx、τxy、τyy為在xx、xy、yy方向上的應(yīng)力有效力。

2.2 泥沙轉(zhuǎn)移模型

通過水深均一的二維模型(模型中非粘性泥沙顆粒的直徑要求達(dá)到0.063 mm以上)來模擬泥沙轉(zhuǎn)移。該模型通常用以分析湖泊水庫(kù)中非粘性泥沙的絮凝，同時(shí)也可用以分析泥沙轉(zhuǎn)移對(duì)水底高程的改變。模型如下：

(4)

其中：ρ為垂直方向上的泥沙平均質(zhì)量密度；u、v為水流沿x、y方向上的流速；Dx、Dy為泥 沙的擴(kuò)散系數(shù)；Ti為源匯項(xiàng)。

2.3 富營(yíng)養(yǎng)化模型

模型公式如下：

(5)

其中：C1為富營(yíng)養(yǎng)化垂向平均濃度；Dx、Dy為x、y方向的擴(kuò)散系數(shù)；kp為線性衰減率；CS為源濃度；S為點(diǎn)源排放量；Pc為富營(yíng)養(yǎng)化過程組。

(6)

其中：ci為富營(yíng)養(yǎng)化濃度；m為狀態(tài)變量過程數(shù)。

3 汾河水庫(kù)含沙水質(zhì)模擬

3.1 水動(dòng)力模擬

3.1.1 模擬范圍

選取全庫(kù)水域?yàn)槟M范圍，運(yùn)用GPS和Erdas工具對(duì)水庫(kù)進(jìn)行現(xiàn)場(chǎng)定位，同時(shí)依靠水庫(kù)周邊的地形特點(diǎn)，確定水庫(kù)的邊緣點(diǎn)坐標(biāo)，通過地理信息處理技術(shù)將坐標(biāo)進(jìn)行數(shù)字化處理，確定汾河水庫(kù)的水庫(kù)邊界?？紤]水庫(kù)周邊環(huán)境存在一定復(fù)雜性，依據(jù)非結(jié)構(gòu)化三角形網(wǎng)絡(luò)系統(tǒng)將整個(gè)庫(kù)區(qū)劃分為3 160個(gè)三角網(wǎng)格，垂直方向無分層。水庫(kù)模擬邊界與地形見圖1。

3.1.2 定解條件

3.1.2.1 初始條件：假定全部模擬區(qū)域均以靜止?fàn)顟B(tài)開始，即u=v=0；水位初始值為沿水流方向取得的邊界兩端的平均值，即1 126.2 m。

3.1.2.2 邊界條件：通過岸壁法處理閉邊界，其處理方法為設(shè)置法線方向上的濃度和速度為零；汾河開邊界，采用“干濕點(diǎn)判別法”的方式進(jìn)行分析，增水時(shí)的水深不低于0.1 m，此時(shí)為“濕點(diǎn)”，按水域進(jìn)行處理，退水時(shí)的水深為0.005 m，水速為零，視為 “干點(diǎn)”，按陸域來處理。

3.1.3 模擬參數(shù)率定和模型驗(yàn)證

3.1.3.1 模擬時(shí)間步長(zhǎng)。模擬時(shí)間為2017年2月1日至2017年12月1日。依據(jù)線性方程，得到最終時(shí)間步長(zhǎng)為3600 s，其滿足CFL條件。

3.1.3.2 渦黏系數(shù)。用Smagorinsky[10]公式確定水平渦黏系數(shù)，其值為0.8。

3.1.3.3 摩擦力。通過曼寧公式確定摩擦力，曼寧系數(shù)值取32 m1/3/s。

3.1.3.4 降雨量與蒸發(fā)量。查閱汾河水庫(kù)庫(kù)區(qū)相關(guān)資料，獲得在時(shí)間尺度上的降雨量和蒸發(fā)量。

3.1.3.5 源與匯。取進(jìn)入汾河水庫(kù)的水流為源匯項(xiàng)，北部汾河入水口和西部澗河入水口為汾河水庫(kù)的主要進(jìn)水口，汾河水庫(kù)壩址處為出水口。每年7～8月，汾河水庫(kù)處于枯水期，而在1、5、6、9～12月，汾河水庫(kù)處于平水期，2～4月，汾河水庫(kù)處于枯水期。豐水期，入庫(kù)流量最大為60 m3/s，枯水期入庫(kù)流量最小為4 m3/s。汾河水庫(kù)的主要目的為下游提供灌溉用水和日常生活飲用水[11]。每年3月中旬至4月中旬、8月中旬至9月初，下游開閘放水灌溉，出口流量可達(dá)20～40 m3/s，其余時(shí)間則主要為古交礦區(qū)和太原市提供日常生活飲用水，此時(shí)水庫(kù)的出水流量約為3.12 m3/s。依據(jù)不同時(shí)間段水量的變化，形成相應(yīng)的時(shí)間序列。

3.1.3.6 參數(shù)率定和模型驗(yàn)證。以汾河水庫(kù)的表面水位值為驗(yàn)證資料，經(jīng)過模型參數(shù)調(diào)整，確定計(jì)算參數(shù)。模擬輸出了2017年7月～2017年9月水庫(kù)表面水位值時(shí)間序列，將輸出水位值與2017年《山西水文月報(bào)》中的實(shí)際水位值進(jìn)行對(duì)比，結(jié)果表明6月11日、6月20日、7月11日三天汾河水庫(kù)的實(shí)際測(cè)量水位為：1 126.03、1 126.18、1 126.50 mm，通過模擬，所得到的水位為：1 126.00、1 126.15、1 126.48 mm。模擬表面水位值與水位實(shí)測(cè)值絕對(duì)誤差小于0.03，相對(duì)誤差小于0.003%，因此參數(shù)設(shè)置合理，模型可用性較高。

3.2 水動(dòng)力—泥沙輸移—富營(yíng)養(yǎng)化模型

3.2.1 模擬范圍

模擬范圍為全庫(kù)水域。

3.2.2 模型參數(shù)

3.2.2.1 模擬時(shí)間步長(zhǎng)、渦黏系數(shù)、摩擦力、降雨量與蒸發(fā)量參數(shù)取值同水動(dòng)力模擬一致。

3.2.2.2 源與匯。與水動(dòng)力模擬所建立的時(shí)間序列一致之外，還需要對(duì)TN、TP以及Chla等進(jìn)行檢測(cè)并建立時(shí)間序列文件。例選取北部上游Chla濃度并建立時(shí)間序列，濃度介于0.021～0.048 mg/L。

3.2.3 耦合模型的參數(shù)率定與驗(yàn)證

3.2.3.1 水動(dòng)力—泥沙耦合模型

選取2017年6月11日～2017年6月30日B點(diǎn)位置處的數(shù)據(jù)進(jìn)行驗(yàn)證，通過實(shí)際數(shù)據(jù)的輸入和調(diào)整，從而確定準(zhǔn)確數(shù)據(jù)，數(shù)據(jù)類型主要包括：孔隙率系數(shù)為0.2，相對(duì)密度系數(shù)1.2，曼寧系數(shù)32 m1/3/s，擴(kuò)散系數(shù)1，模型模擬輸出結(jié)果與實(shí)際監(jiān)測(cè)值對(duì)比結(jié)果如圖2(a)所示。

3.2.3.2 水動(dòng)力—富營(yíng)養(yǎng)化耦合模型

選用2017年8月1日～2017年8月10日實(shí)地測(cè)量B監(jiān)測(cè)點(diǎn)的Chla、TN、TP、DO等，通過模型對(duì)數(shù)值進(jìn)行驗(yàn)證。因篇幅限制，本研究?jī)H以Chla為例，在誤差可接受范圍內(nèi)，對(duì)水質(zhì)模型相關(guān)參數(shù)進(jìn)行驗(yàn)證，結(jié)果見圖2(b)。

以B點(diǎn)泥沙濃度為例，模擬值與實(shí)測(cè)值基本吻合。最大絕對(duì)值偏差7.5 g/m3，相對(duì)誤差1.05%。因此本文構(gòu)建的耦合模型具有可行性。

圖2 B點(diǎn)模擬值與實(shí)測(cè)值對(duì)比Fig.2 Comparison of B the simalated value and measured value at point B

3.3 方案的設(shè)定

設(shè)置不同的模型方案，對(duì)各個(gè)方案進(jìn)行對(duì)比分析，選出具有代表性的控制點(diǎn)分析汾河水庫(kù)水質(zhì)受到泥沙的影響。從空間距離上，依據(jù)庫(kù)區(qū)入水口離壩址的遠(yuǎn)近順序選取了6個(gè)代表性控制點(diǎn)，如圖1所示。設(shè)定三種方案，方案1：水動(dòng)力—泥沙耦合模型 ，為①；方案2：水動(dòng)力—富營(yíng)養(yǎng)化耦合模型，為②～④；方案3：水動(dòng)力—泥沙—富營(yíng)養(yǎng)化耦合模型，為⑤～⑦，方案設(shè)置分類見表1。其中，降水量、蒸散量、風(fēng)力等，采用實(shí)際測(cè)量數(shù)據(jù)。模擬時(shí)間：2017年6月11日至2017年9月23日，共計(jì)104天。

3.4 模擬結(jié)果分析

3.4.1 模擬結(jié)果

在水動(dòng)力模擬基礎(chǔ)上，Chla、TP、TN、DO濃度模擬分布圖，見圖3。

圖3 Chla、TP、TN、DO濃度分布Fig.3 Distribution of Chla、TP、TN、DO concentration

3.4.2 結(jié)果分析

3.4.2.1 通過水動(dòng)力—泥沙耦合模型模擬，分析泥沙濃度變化情況。以下為模擬結(jié)果：選擇6個(gè)模擬點(diǎn)，以進(jìn)水口位置為基準(zhǔn)點(diǎn)，A、B較近，其次為C、D、E、F點(diǎn)，結(jié)合位置分為兩組，具體仿真結(jié)果見圖4。

圖4 各監(jiān)測(cè)點(diǎn)泥沙濃度變化圖Fig.4 Sediment concentration variation of monitoring points

汛期SSC的時(shí)空分布：A點(diǎn)泥沙濃度的急劇增加主要來自汾河干流來水的直接補(bǔ)給；B點(diǎn)泥沙濃度的變化較為緩和。A、B兩監(jiān)測(cè)點(diǎn)的研究結(jié)果表明：2016年7月25日前后，兩監(jiān)測(cè)點(diǎn)的SSC濃度出現(xiàn)了劇烈增加，之后泥沙濃度達(dá)到峰值。實(shí)時(shí)監(jiān)測(cè)結(jié)果顯示，7月22日至7月28日，汾河干流的進(jìn)水量為30～60 m3/s，澗河處，其進(jìn)水量達(dá)到了本年最大值，為1.13～2.26 m3/s，這直接影響了SSC濃度，使其變化十分明顯。A點(diǎn)與進(jìn)水口距離較近，水體中泥沙的含量變化大，相反，C、D、E、F點(diǎn)與進(jìn)水口距離遠(yuǎn)，泥沙的濃度變化可以忽略不計(jì)。綜上所述：含沙水在某種程度上會(huì)對(duì)水庫(kù)泥沙濃度產(chǎn)生影響，同時(shí)這種影響的范圍也是有限的。

3.4.2.2 通過水動(dòng)力—富營(yíng)養(yǎng)化耦合模型和水動(dòng)力—泥沙—富營(yíng)養(yǎng)化耦合模型模擬，并選取具有代表性的監(jiān)測(cè)點(diǎn)位，模擬不同時(shí)段下的Chla濃度、TP濃度、TN濃度、DO濃度，見圖5。

圖5 各代表性監(jiān)測(cè)點(diǎn)Chla、TP、TN、DO濃度變化圖Fig.5 Variatiion of Chla、TP、TN、DO concentrations at representative monitoring points

(1)SSC對(duì)水庫(kù)Chla的影響：汛期前期，Chla含量受到SSC的影響基本可以忽略不計(jì)；A、B點(diǎn)Chla含量在SSC大量存在前后差異不大，并趨于吻合。此外，在汛期前期，如果水庫(kù)的水量進(jìn)出發(fā)生變化的話，會(huì)對(duì)Chla的含量產(chǎn)生直接影響。汛期期間，B點(diǎn)高含量SSC存在時(shí)Chla含量比SSC不存在時(shí)越來越小。汛期過后，B點(diǎn)Chla濃度較非含沙水高。另外，由于E、F監(jiān)測(cè)點(diǎn)距離進(jìn)水口較遠(yuǎn)，這就導(dǎo)致水體中泥沙的含量對(duì)Chla的含量并無明顯影響。綜上，整體來說，在整個(gè)汾河水庫(kù)體系中，水庫(kù)進(jìn)水對(duì)于Chla的含量影響并不是很大。

(2)SSC對(duì)水庫(kù)TP的影響：SSC大量存在和SSC不存在兩種情況下，各監(jiān)測(cè)點(diǎn)位的TP含量總是前者高于后者。

(3)SSC對(duì)水庫(kù)TN的影響：C點(diǎn)TN含量差異不大，D和E點(diǎn)TN的含量差距明顯，當(dāng)TN的濃度上升后，其濃度就很難下降，這種現(xiàn)象與TP的變化規(guī)律存在明顯區(qū)別，當(dāng)某些監(jiān)測(cè)點(diǎn)受SSC的干擾作用較弱時(shí)，這種現(xiàn)象與水體中TP含量的時(shí)空分布規(guī)律是一致的，當(dāng)進(jìn)入水庫(kù)的水不含SSC時(shí)，水體中TN的含量就會(huì)遠(yuǎn)大于進(jìn)水存在大量泥沙的狀況。

(4)SSC對(duì)水庫(kù)DO 的影響：研究結(jié)果表明，SSC幾乎對(duì)DO的含量不會(huì)產(chǎn)生影響，即使當(dāng)水體中含有大量的SSC時(shí)，水體中的DO含量也幾乎不會(huì)發(fā)生任何變化。因此泥沙存在與否，二者對(duì)DO含量的影響差異性不大。

3.5 底泥污染物釋放量計(jì)算

根據(jù)經(jīng)驗(yàn)公式進(jìn)行底泥污染物釋放量估算。

M=R×S×t×10-3

(6)

其中：M為底泥污染物釋放量，t；R為污染物釋放通量，mg/m2.d；S為庫(kù)區(qū)面積，km2；t為釋放時(shí)間，d，取t=104 d。

M’=C×S×h+C×Q×10-6

(7)

其中：M’為水體污染物總負(fù)荷，t；C為水體污染物平均濃度，mg/l；h為平均水深，m；Q為年出庫(kù)水量。

由底泥污染物釋放量計(jì)算得出，汛期后期，有2.5 t的污染物來自于底泥釋放的TN，有4.16 t的污染物來自于底泥釋放的TP，水庫(kù)庫(kù)區(qū)TN、TP的總負(fù)荷分別為615.33 t和16.47 t，此兩者的貢獻(xiàn)率分別為0.406%和25.25%，這一結(jié)果說明水庫(kù)水體磷負(fù)荷以及富營(yíng)養(yǎng)化是造成底泥內(nèi)源磷釋放的重要因素，需積極開展污染治理。

3.6 污染防治對(duì)策

(1)預(yù)防性措施：在汾河水庫(kù)上游采取攔沙方式，減少泥沙入庫(kù)，是控制懸浮泥沙入庫(kù)的最根本方法。其中水土保持，防止流失，對(duì)減少水庫(kù)淤積有著關(guān)鍵作用，同時(shí)還可以減少面源污染物，改善生態(tài)環(huán)境質(zhì)量。(2)清淤措施：采用“蓄清排渾”的運(yùn)行方式，在豐水豐沙年內(nèi)，利用高含沙量的水流，在保持一定低水位的條件下，進(jìn)行異重流排沙，通過及時(shí)打開排沙閘門的方式使泥沙高效排出，有利于改變下游水生生物的生長(zhǎng)與繁殖。

4 結(jié) 論

4.1 基于MIKE21，構(gòu)建了水動(dòng)力-泥沙轉(zhuǎn)移-富營(yíng)養(yǎng)化模型，研究污染物的轉(zhuǎn)移與擴(kuò)散過程、SSC的轉(zhuǎn)移規(guī)律等，并分析了水庫(kù)的水質(zhì)情況。

4.2 根據(jù)模型結(jié)果，對(duì)水庫(kù)的TN、TP、Chla及DO四項(xiàng)指標(biāo)綜合分析，同時(shí)結(jié)合SSC在時(shí)間、空間上的分布規(guī)律，依靠對(duì)比分析方式分析水庫(kù)水質(zhì)受到含沙水和不含沙水的影響，闡明了泥沙對(duì)水質(zhì)產(chǎn)生的影響。

4.3 水庫(kù)水動(dòng)力—泥沙轉(zhuǎn)移—富營(yíng)養(yǎng)化耦合模型模擬結(jié)果表明，SSC對(duì)TN、TP的影響最大，對(duì)Chla、DO濃度的分布影響很小。

4.4 通過計(jì)算底泥的釋放量，得出TP對(duì)水庫(kù)的影響最大，針對(duì)汾河水庫(kù)的特征，分別從預(yù)防性措施和清淤措施兩方面對(duì)汾河水庫(kù)提出針對(duì)性的污染防治對(duì)策，其中保持水土，防止流失，對(duì)控制懸浮泥沙入庫(kù)起著關(guān)鍵作用。