万家寨水库水环境特性分析

高玮岐，胡 煜，任华堂，Tomasz TYMISKI

(1.中央民族大学生命与环境科学学院，北京 100081；2.弗罗茨瓦夫环境与生命科学大学环境工程学院，弗罗茨瓦夫 50-363)

0 引 言

万家寨水库位于黄河北干流上段托克托至龙口河段峡谷内，是太原、大同、朔州3市的工业和生活用水主要水源地。近年来上游取水量增加，且上游工业排污持续增加，水环境污染呈现出日益严重之势[1]，严重威胁山西人民的生产生活。开展万家寨水库水环境纳污特性的研究，对于引水区域的供水安全具有重要的现实意义。

对于万家寨水库的研究，现有成果主要集中于凌汛期水动力特性的研究。冀鸿兰[2]等根据历史实测冰情资料，统计分析得出上游河段冰情特征规律，给出上游河段防凌调度原则；张宝森[3]等以水泥厂为观测点，对凌汛期库尾水位变幅增大受库区水位因子及环境因子影响的比重进行研究发现，岔河口形成冰塞、冰坝是上游水位突升的主要影响因子；梁贵生[4]等研究发现，库区冰塞头部形成位置取决于河道条件、水库水位和河道流量。近年，对于万家寨水库水环境污染的研究日渐重视，漆强强[5]选取具代表性断面，进行pH值、高锰酸盐、化学需氧量、总氮等污染因子的监测发现，大部分断面水质指标均有不同程度的超标现象；吴喜军[6]等应用改进的内梅罗污染指数法、灰色聚类法和单因子评价法对库区水质进行评价，得出黄河干流整体水质状况较好的结论；王若仪[7]应用营养分级评分法和综合营养状态指数法进行总磷、总氮等5项污染因子的评价发现，万家寨水库的叶绿素a、透明度已达到轻度富营养水平，总氮达到中度富营养水平。

目前，对于研究区域水环境污染的研究成果多为水质总体评价，对于库区水环境特性精细的时空分布研究相对匮乏。为此，本文应用水环境数学模型，对万家寨水库库区水环境进行模拟，分析水动力特性和水环境特性，并计算典型工况下区域内COD最大排放浓度分布及水环境容量。

1 数学模型

1.1 控制方程

采用二维水流水质模型进行求解[8]，即

(1)

(2)

(3)

式中，η为水位高度；D为单元的全水深；u、v分别为流速矢量U在x、y方向的分量；g为重力加速度；Cb为拖曳力系数，与水面的粗糙度有关；ρ为水体密度；AM为涡粘系数；f=2Ωsinφ为科氏力参数，Ω为地球自转角速度，φ为地理纬度。

采用对流扩散方程模拟区域污染物输移扩散规律，即

(4)

式中，C为污染物垂向平均浓度；AH为水平扩散系数；S为源项。

应用有限体积法(FVM)对于上述物质连续方程、动量方程以及浓度方程进行离散，主要包括对对流及扩散通量做体积积分和时间积分处理，根据对流通量、扩散通量、源(汇)项的表达式即可建立各单元的物质守恒方程，将所有单元的浓度守恒方程进行联立求解，进而得知全场污染物的浓度分布特性。

1.2 网格划分及计算条件

计算区域为39°34′4″N～40°13′05″N，111°10′23″E～111°27′06″E，上游边界取在拐上断面，下游边界取在坝上断面。研究区域网格划分见图1。采用二维水动力模型进行计算，库区平面应用非结构三角形网格进行贴体划分，其划分网格10 476个，三角形网格节点数为6 109个，三角形界面数16 584个。沿河宽方向网格间距为33.0～66.5 m，平均间距为51.1 m；沿河道方向网格间距为38.7～62.2 m，平均间距为53.6 m。

图1 万家寨库区网格划分

1.3 计算条件

根据已掌握的万家寨水库长期监测数据，采用上游入库流量、下游坝前水位边界进行控制，将拐上断面多年月均监测资料作为模型入库的水质边界条件。模型计算时间步长为0.1 s，底部拖曳力系数为0.025。参考研究库区相似河段及相关水库水质模型中有关各污染物衰减系数的研究成果[9-10]，根据季节的差异，非冰封期化学需氧量(COD)的降解系数设为0.240 d-1，氨氮设置为0.105 d-1。

2 模型验证

从水动力和水质2个方面进行验证。采用水库库容日变化过程对水动力结果进行验证，见图2。从图2可知，实测值和计算值平均相对误差为1.08%，具有较好的一致性。

图2 模拟库容与实测库容验证

选取COD与氨氮浓度进行不同水文年条件下水质模型坝前污染物浓度验证，见图3、4。利用平均对数比偏差Dr计算计算值与实测值的差异性，公式如下

(5)

图3 不同水文年条件下水质模型COD浓度验证

图4 不同水文年条件下水质模型氨氮浓度验证

式中，DLm为实测值；DLp为计算值；n为样本个数。

经计算，各种条件下2种水质指标计算值与实测值之间的平均对数比偏差均在0.1～0.3范围内，最大为0.297，满足工程应用精度要求。

图5 不同水文年下夏季坝前30 km范围内流速分布状态

3 结果及分析

3.1 水动力场分析

图5给出了典型丰水年、平水年及枯水年万家寨水库坝前30 km区域内夏季流速分布状态。从图5可知，不同水文年流速差异较大，丰水年流速>平水年流速>枯水年流速，这是由于不同水文条件下的流量差异所致。此外，由于紊流作用，在曲率较大的岸边界处会形成漩涡，由于离心惯性力和重力导致局部水体出现掺混现象，有利于水体交换。

3.2 水力停留时间

水力停留时间(HRT)是反映水体交换能力的重要参数，对水体自净有着重要影响[11]，计算公式如下

(6)

式中，HRT为水力停留时间；V为水库库容；Q为入流流量。

库区水力停留时间计算统计见表1。从表1可知，在典型丰水年、平水年、枯水年3种水文年条件下，水库段水力停留时间在1.76～11.78 d范围内，万家寨水库具备典型的河道水库型水动力特征。具体而言，夏季水力停留时间小于冬季水力停留时间，最短水力停留时间出现在丰水年夏季，为1.76 d；最长水力停留时间出现在枯水年冬季，为11.78 d。产生此结果的原因是水库在丰水年夏季具有较常年各季节偏高的径流量，水力停留时间较其他时期短。枯水年冬季具有较常年各季节偏低的径流量，水力停留时间较其他时期长。

表1 库区水力停留时间计算统计

3.3 COD的最大排放浓度及水环境容量计算

3.3.1计算工况

为了探究水库水环境容量，假设万家寨水库上游入流断面处以固定流量连续排放一定浓度的有机污水，设计3种不同的计算工况，见表2。

根据库区HJ 338—2018《饮用水水源保护区划分技术规范》，饮用水源一级保护区的水质应满足二类水质标准，其中水质指标COD应小于15 mg/L。在此要求下，本文采用控制断面法，将下游饮水水源一级保护区的3个监测断面作为控制断面(见图6)，根据控制断面的水质要求确定水环境容量。

表2 连续源污染排放计算设计工况

图6 饮用水源保护区控制断面

3.3.2污染物浓度空间分布

在满足水库下游3条控制断面水质指标COD浓度均达到二类水质标准(15 mg/L)的条件下，利用数学模型计算了枯水期、平水期、丰水期3种水文条件下的最大允许排放浓度，结果分别为426、310 mg/L与252 mg/L。在该污染负荷条件下，COD浓度的空间分布见图7。从图7可知，总体而言，污染物浓度随着水流方向纵向浓度逐渐降低，库尾处污染物COD浓度最大，入库后污染物随水流向下游迁移，但随着污染物运动过程中迁移扩散及降解作用，沿水流输运方向污染物浓度会逐渐减小。

图7 不同水文条件下区域COD浓度空间分布

丰水期、平水期、枯水期污染带上下游浓度空间差异分别为230.59、279.65 mg/L和289.48 mg/L。丰水期入库流量较大，由于推移作用强，纵向迁移距离较长；枯水期入库流量较小，推动作用弱，污染物拉伸范围较小，浓度空间分布更不均匀，浓度空间差异较大。

3.3.3水环境容量

(1)计算区域。选择入库处(拐上断面)到坝前共73 km为研究水功能区域。水环境容量W采用下式计算

(7)

式中，Q为设计流量；C0为段首污染物浓度；(t2-t1)为时间间隔。

(2)不同水文条件下水环境容量分析。设计工况下水环境容量在不同水文条件下的计算结果见表3。其中，入库COD浓度分别为在枯水期、平水期、丰水期3种水文条件下，连续污染源COD的最大临界排放浓度。从表3可知，在枯水期、平水期、丰水期设计流量下，下游3个控制断面COD均能达到二类水质标准(15 mg/L)时，入库污染物COD日均最大负荷总量分别为14 060、15 959 t和18 137 t。在等长河段内，水质目标相同的前提下，不同水文条件下COD的水环境容量差异明显，丰水期>平水期>枯水期。产生此差异的原因是污染物的迁移受流速的影响，流量越大，流速越大，污染物的对流扩散及迁移能力越强；枯水期流量小，水流流速较小，污染物迁移速度较为缓慢，容易聚集在库区内；丰水期水流流速较大，污染物扩散迅速，水力停留时间短，在库区聚集的可能性小，所以丰水期污染物COD日均负荷总量最大。该水环境容量仅根据下游断面的水质要求计算得到，在水库管理中还需根据整个库区的水质要求确定可容纳的污染负荷。

表3 库区水环境容量

4 结 语

本文采用水环境数学模型，对万家寨水库库区水动力特性和水环境特性进行了分析，计算出典型工况下区域内污染物浓度空间分布和水环境容量。在固定源连续排放下，本文设定典型的水文条件，利用数学模型计算，得到连续污染源的COD最大允许排放浓度与相应的COD负荷水环境容量，为库区污染物总量消减提供参考依据。