基于EFDC的潘家口水库富营养化模拟

徐婉珍

摘要 基于EFDC建立潘家口水库水动力与富营养化模型，模拟了2006～2007年潘家口水库内的水动力、营养物质循环以及藻类生长代谢，并对模型进行了校准与验证。结果表明，TN的模拟效果最好，最大相对误差为4.97%；TP的模拟效果次之，最大相对误差为-15.77%；而Chla的模拟效果最差，最大相对误差18.29%。模拟与实测结果对比表明，基于EFDC建立的潘家口水库富营养化模型，可以比较真实地反映库区内营养物质的循环以及藻类的生长代谢过程。

关键词水动力；水质；富营养化；潘家口水库；EFDC

中图分类号S181文献标识码A文章编号0517-6611（2015）34-115-04

潘家口水库位于海河流域滦河水系的滦河干流上，地跨河北省的承德、唐山两市，是天津市、唐山市重要的水源地，承担着提供两地工业、农业、城市用水的重任。截至2013年，潘家口水库已经向天津市、唐山市以及滦下灌区提供优质水资源近366亿m3，其中向天津市供水167亿m3，向唐山以及滦下灌区供水199亿m3，成为天津市和唐山市可靠的供水水源，为天津市和唐山市的可持续发展做出了重要贡献[1]。

随着区域经济的发展，人民生活水平的提高，潘家口水库所在流域的点源污染与面源污染的输入越来越被关注与重视。上游地区城市生活污水对水源地造成严重污染。农作物使用的大量农药、化肥，随着上游的来水直接进入库区，造成库区氮、磷含量明显增加，富营养化程度升高。水库周边大量的工厂也导致了大量污染物的输入。同时，网箱养鱼成为重要的内源，造成的威胁不容小觑，在经济利益的驱使下，潘家口水网箱养鱼从2003年至今呈现快速扩大的趋势，养殖密度几乎超过了水库的承载能力。投放的饵料以及产生的鱼类粪便中含有大量的氮、磷，加速了水库的富营养化进程[2]。

水库作为调节、蓄水和控制水资源的人造系统，是实现水资源合理开发利用的有效途径，而水库相对封闭，水体停留时间较长，营养盐循环受到干扰，大量的氮磷元素排入水体，引起水体中藻类、浮游生物的迅速繁殖，从而导致水体富营养化，严重影响居民的饮水安全[3]。为明晰潘家口水库内的富营养化过程，笔者利用EFDC模型建立潘家口水库富营养化模型，对潘家口水库内的水动力、营养物质的循环以及藻类的生长代谢做了较系统的模拟。

1模型简介

EFDC（The Environmental Fluid Dynamics Code）模型是美国国家环保署（USEPA）推荐的三维地表水水动力模型[4-5]，可实现河流、湖泊、水库、湿地系统、河口和海洋等水体的水动力学和水质模拟，是一个多参数有限差分模型。EFDC模型采用MellorYamada 2.5阶紊流闭合方程，根据需要可以分别进行一维、二维和三维计算。模型包括水动力、水质、有毒物质、底质、风浪和泥沙模块，用于模拟水系统一维、二维和三维流场、物质输运（包括水温、盐分、粘性和非粘性泥沙的输运）、生态过程及淡水入流，可以通过控制输入文件进行不同模块的模拟。模型在水平方向采用直角坐标或正交曲线坐标，垂直方向采用σ坐标变换，可以较好地拟合固定岸边界和底部地形。在水动力计算方面，动力学方程采用有限差分法求解，水平方向采用交错网格离散，时间积分采用二阶精度的有限差分法，以及内外模式分裂技术，即采用剪切应力或斜压力的内部模块和自由表面重力波或正压力的外模块分开计算。外模块采用半隐式三层时间格式计算，因传播速度快，所以允许较小的时间步长。内模块采用考虑了垂直扩散的隐式格式，传播速度慢，允许较大的时间步长，其在干湿交替带区域采用干湿网格技术。该模型提供源程序，可根据需要对源程序进行修改，从而达到最佳的模拟效果。

水动力学方程是基于三维不可压缩的、变密度紊流边界层方程组，为了便于处理由密度差而引起的浮升力项，常常采用Boussinesq假设。在水平方向上采用曲线正交坐标变换，在垂直方向上采用sigma坐标变换后，得到以下控制方程：

动量方程为：

（mHu）t+（myHuu）x+（mxHvu）y+

（mwu）z-（mf+vmyx-umxy）Hv=

-myH（gζ+p）x-my（hx-zHx）pz+z（m1HAvuz）+Qu（1）

（mHv）t+（myHuv）x+（mxHvv）y+

（mwv）z+（mf+vmyx-umxy）Hu=

-mxH（gζ+p）y-mx（hy-zHy）pz+z（m1HAuvz）+Qv（2）

pz=-gHρ-ρ0ρ0=-gHb（3）

連续方程为：

（mζ）t+（myHu）x+（mxHv）y+（mw）z=0

（4）

（mζ）t+（myH∫10udz）x+（mxH∫10vdz）y=0

（5）

ρ=ρ（P，S，T）（6）

EFDC模型中的富营养化模块模拟了21种物质的状态变量，包括3种藻类中碳、氮、磷组分、硅循环、溶解氧以及粪便大肠杆菌。水体中的污染物控制方程由溶于水中的污染物、吸附于水中可溶性物质上的污染物以及吸附于水中悬浮颗粒物上的污染物共同构成。溶于水中的污染物迁移控制方程如下：

t（mxmyHCw）+x（myHuCw）+y（mxHvCw）+z（mxmywCw）

=z（mxmyAbHzCw）+mxmyH（ikidSSiχis+jKjdDDjχjD）-

mxmyH[i（KiaSSi）（ψwCw）（iS+χiS）+j（KjaDDj）（ψwCw）（jD-χjD）+γCw]

（7）

式中，Cw为溶解于水中污染物的单位体积质量浓度；χS为吸附于i类悬浮物的污染物单位体积质量浓度；χD为吸附于j类溶解物质的单位体积质量浓度；为孔隙率；ψw为可吸附溶解于水中污染物的水所占比例；Ka为吸附速率；Kd为解吸速率；γ为净线性衰减系数；S为水中悬浮物单位体积质量浓度；D为水中溶解物质单位体积质量浓度。

EFDC的富营养化模块包含3种藻类：蓝藻、绿藻和硅藻，描述这些藻类动力的过程大致相同，只是在不同藻类的方程中改变了参数值。模块中的源汇项主要是藻类的生长、新陈代谢、捕食、沉积和外来负荷，方程如下：

tBc=（Pc-BMc-PRc）Bc+Z（WScBc）+WBcV（8）

式中，Bc为藻类生物量（g C/m3）；t为时间（d）；Pc为藻类的生产率（d-1）；BMc为藻类新陈代谢速率（d-1）；PRc为藻类捕食率（d-1）；WSc为沉积速率（m/d）；WBc为藻类的外源输入（g C/d）；V为网格单元体积（m3）。

2模型建立与数据来源

2.1模型建立

清河口断面至潘家口水库大坝区域水面开阔、水深较大，污染物浓度在水平及垂直方向上均有较大的差别，为了能充分揭示污染物在潘家口水库的迁移扩散规律，需要进行三维模拟。其中，水平方向上划分为637个单元格，垂直方向上平均划分为4层，共计2 916个单元格（图1）。模型共有5处入流，分别为清河口、瀑河口、石泉浩、沟门以及菜子沟门，瀑河口与菜子沟门有水质点源输入，除此外的入流根据水质监测数据输入水质组分的浓度时间序列。校核点为燕子峪与潘家口两处。水库出流为潘家口水库大坝，采用上层出流。

图1潘家口水库EFDC模型边界条件与校核点示意图

2.2数据来源

除2006～2007年潘家口水库库区的气象资料来源于气象资料分享网站和流域年鉴外，

潘家口水库库底地形图、水文测站基本信息、2006～2007年潘家口水库水文监测资料（水位、流量）、水质监测数据和点源入流的数据均来源于水利部海河水利委员会。

3模型的校核与验证

水动力的模拟是潘家口富营养化模型的基础。水动力的模拟决定了物质运输的路径以及物质混合的方式，这是模型首先要评价的初步内容。因为缺乏流速资料，仅以潘家口水库大坝坝上水位进行水动力的校验。结果表明，模拟值与实测值吻合良好（图2），平均误差为0.056 m。

EFDC的富营养化模拟了蓝藻、绿藻和硅藻在太阳辐射、水温、营养盐水平影响下的生长、代谢、死亡。每种藻类都有其自己的生长、演替规律，藻动力模型通过改变与生产、代谢、布施、沉积速率等相关的关键参数来模拟不同的过程。细胞个数和叶绿素浓度这两个指标可以描述水体中的藻类生物量。细胞计数、染色剂浓度都是计量藻类生物量的通用方法。直接的细胞数计量与体积计量十分费力，需要专门的生物分类以及样品保存。单独某种藻类的叶绿素浓度测算也十分复杂。但是，水体中的叶绿素浓度可以直接指示生态系统中的生物量与富营养化水平，因此模型采用叶绿素浓度代表总藻类生物量。

富营养化模型的校准关键在于藻类初级生产力、有机物分解、营养盐浓度的模拟。模型校准的目标是选择与藻类生长、营养盐转化、呼吸过程相关的最合适的参数组合，减小模拟值与实际值的误差。对浮游生物体以及营养盐输移转化过程的理解，有助于提高富营养化过程模拟的准确性与精确性。模型初步运行应用相关文献推荐的经验参数，逐渐调整个别参数，让模拟值与实测值间接近吻合。模型的验证是用另外一套独立的数据，在不同的外界条件下检验模型反映水环境系统变化的能力。藻类动力过程的关键参数包括生长速率、新陈代谢速率以及捕食率，这3个参数描述了藻类生物量的增加与减少。极端的高温与低温会严重影响藻类浓度。此外，营养盐浓度也会影响藻类的生长速率。富营养化模型中的各个状态变量也取决于硝化过程、反硝化过程以及有机物的水解过程。

在模型校验之后进行敏感性测试，可以加深对模型表现

的理解，并檢验模型设置的关键参数。结果表明，藻类生产率是最敏感的内部参数，磷的半饱和浓度次之。敏感性分析显示，10%藻类生产率的变化可以导致大于43.5%的叶绿素浓度变化。藻类生产浓度作为内部参数，取决于藻的种类及其特征，而磷的半饱和浓度反映了外界营养盐负荷对藻生物量的影响。在模拟中，TP负荷发生10%的扰动，会导致叶绿素浓度至少35.7%的变幅；TN发生10%的扰动，导致叶绿素浓度大约13.2%的变幅，所以可以认为在模型中磷负荷的输入是影响藻类生长最敏感的外因。

在模型的校准和验证中，将模拟值与实测值相比较，计算平均误差（ME）与平均相对误差（MRE），从而评价模拟的效果。具体计算公式为：

ME=1nni=1（Oi-Pi）（9）

MRE=1nni=1Oi-PiOi×100%（10）

式中，Oi为实测数据；Pi为模拟数据；n为校准与验证数据的数目。校核与验证结果表明，TN的模拟效果最好，TP次之，Chla的模拟偏差最大。对于Chla，在校准时期燕子峪的最大相对误差达18.29%，在验证时期潘家口的最大相对误差达18.13%；TP模拟的最大相对误差在校准时期与验证时期分别达-15.77%、15.52%；TN模拟效果最好，最大相对误差出现在验证时期的燕子峪，达4.97%（表1）。TN、TP、Chla模拟的相对误差都控制在20%以内。模拟结果显示，该模型可以描述潘家口水库表层的藻类生长。在校准、验证中，营养盐循环的模拟有着更好的吻合度。

由图3～6可知，2006和2007年Chla模拟值与实测值呈现较好的一致性，在120 d左右，低温环境逐渐改变，藻类生命活动逐渐活跃起来，随着气温回暖，生长速率、代谢速率、繁殖速率都逐渐上升，180～240 d Chla浓度达到峰值。Chla浓度的变化过程主要依赖于水温的变化，180～240 d是北方的夏季，温度为25～30 ℃，正好蓝藻、绿藻、硅藻的最适温度分别为27.5、20、25 ℃左右；同时，充足的太阳辐射为藻类的光合作用提供了足够的能量。该段时间内，Chla浓度呈现波动变化，原因是7、8月份正是北方雨水充沛的时节，流域上游的大量来水一方面会在短时间内稀释Chla，但又因为大量营养物的输入，最终又促进了藻类的生长，同时大坝放水，也对库内的Chla浓度造成了影响。此外，库内该段时期内浮游动物与鱼类的大量捕食，也是造成Chla浓度波动的原因之一。在240 d之后，温度下降是造成Chla浓度下降的主要原因。

营养物质的浓度也受到藻类生长的影响，TN、TP浓度会在120 d左右出现明显下降，因为此时温度升高，藻类生长，开始吸收营养物质，导致在该时间段会出现氮、磷浓度的明显下降。但随着温度逐渐升高，藻类生长逐渐旺盛，TN、TP浓度却并没有一直下降，这是因为上游来水也逐渐增大，同时带来流域内大量的营养物质，随着雨季的到来，氮、磷的浓度在7、8月份达到最高峰，因为此段时间常伴随着水库阶段性大量的放水，导致了库内TN、TP浓度的波动。240 d后雨季结束，氮、磷营养物质的输入大幅度下降，导致240 d左右TN、TP的浓度又一次明显的下降，而300 d后又普遍呈上升的趋势，这是由浮游植物的死亡与分解引起的。而在0～60与330～360 d大约为期3月的时间内，库区冰封，既没有明显影响库内营养物质浓度的水量、营养物质输入，也没有浮游植物的生长代谢，也没有水库的开闸放水，处于非常平稳的状态，变化微小，呈缓慢下降的趋势，可能是因为吸附在颗粒物上的营养物质随着颗粒的沉淀而进入了水底沉积物中，脱离了水体。