人工湿地植物水质净化作用的数值模拟研究

杜彦良，张双虎，王利军，毛战坡，马晓宇

（1.中国水利水电科学研究院，北京 100038；2.北京市水利规划设计研究院，北京 100038；3.中国电建集团有限公司，北京 100038）

1 研究背景

人工湿地被广泛用于生态修复及环境改善中，湿地对水体的净化作用是湿地的重要功能之一。关于湿地植物对水体氮磷的去除效果及机理，有大量的实验研究成果［1］。不同类型植物净化效果存在差异［2-3］，合理的湿地植物配置可提高净化效果［4］，植物密度及分布也会影响污染物去除率［5］。总结植物对水体的净化机理主要表现在4 个方面。首先是物理上的阻流作用，导致颗粒沉降、溶质输移快慢及路径发生变化等。其次，植物可改变水体溶解氧（DO）含量，沉水植物的光合作用能直接提高水中DO［6］，挺水植物通过根茎通气组织将叶片产生的氧气输运至根部［7］，通过高精仪器可监测到根部附近的泌氧量［8］。不同植物类型泌氧强度不同，从大至小为：芦苇＞香蒲＞荇菜＞水芹菜＞荷花［9］，受植物光合作用和呼吸作用影响，水体中DO 浓度呈昼夜周期性变化，在根系一定深度的底泥中，形成好氧及厌氧交替变化的环境，有利于硝化及反硝化作用的进行，去除湿地中的氮［10］。第三，植物生长提供了基质微生物的代谢的能量“碳源”［11］，根系环境成为大量微生物及底栖生物的栖息场所，被认为地球上最富有活性界面之一［10］。植物根部附近的好氧环境，促使微生物及底栖生物的生长和繁殖，加快营养物质的有机降解及消耗［12］。 最后，植物的生长吸收水体中的营养盐，有研究指出湿地植物对氮磷的吸收量小于微生物作用对氮磷的去除量［13］。综上，植物的作用可在水动力过程DO 过程以及考虑微生物摄取及植物吸收的氮磷反应过程中体现。

受植物生长的需求，人工表流湿地中适合种植区域的水深通常在0～2 m 左右，有研究认为湿地植物对水质净化效率有限［1］。为提高水体净化能力，优化工程设计，在营造适合植物生长的地形基础上，需对植物类型的配置、种植区域分布及面积大小等进行合理规划，在此过程中数学模型是必不可少的工具之一。目前有不少经典成熟的模型，如MIKE、EFDC、Delft3D 等可实现多种生化反应过程的水质模拟。植物类型及分布对湿地水动力影响的模拟研究相对较多［4，14］。Wang 等［15］开发了一个基于隔间的水动力学和水质数值模型，以研究水力梯度驱动下佛罗里达州北部沼泽湿地硫酸盐和TP的去除。JI 等［16］采用LOEM-CW 模型中考虑微生物和植物生长对磷的吸收，模拟湿地中磷的变化。Galanopoulos 等［17］采用描述了基本的物理化学和生物化学过程的AQUASIM 模型，预测湿地对氮磷去除率。李红艳等［18］等用实验率定湿地的氨氮的硝化速率、反硝化速率、有机氮磷的矿化速率等参数，在一维水动力模型耦合WASP 模型，实现对湿地水质的模拟。

本文应用MIKE 模型，参考EcoLab 中WQ 模版中水质指标的反应关系，建立反应湿地植物特性及作用的“WQ+Veg”模块，计算变量为DO、BOD、NH3-N、硝酸盐氮（NO3-N）、亚硝酸盐氮（NO2-N）和TP。考虑到植物生长区域较为稳定，对计算空间上的植物区和非植物区分别进行数值化，基于湿地植物对水质的净化作用机理，在以下反应过程中添加湿地植物的影响：湿地植物阻碍水体流动的物理过程；湿地种植区植物根部泌氧对DO 补充；湿地植物作用下微生物对氮磷的摄取；植物生长对氮磷的吸收作用。模型对河北省永年洼大量湿地植物种植前后的水质进行模拟，结果可反应植物类型、生长区域、面积等对湿地水体的氮磷响应，为湿地生态修复提供模型方法。

2 研究区域概况

永年洼位于子牙河一级支流滏阳河流域，是继白洋淀、衡水湖之后的华北区域第三大洼淀，是河北省南部唯一的内陆淡水湿地生态系统。受上游滏阳河流域用水及污染影响，永年洼湿地面积萎缩，湿地退化，水体水质恶化［19］。

为恢复湿地功能，改善洼淀内水质，2014年开始对湿地进行了大面积的生态修复工程，首先征地退田，对地形进行改造，营造水生植物的适宜水深，然后蓄水种植芦苇、菖蒲、荷花等水生植物，通过植物的自我繁殖逐渐形成区域分布。洼地积水区4.5 km2，水深1.5～3 m。2017年8月由遥感解译得到的植物分布见图1，芦苇和菖蒲面积15.73 hm2，荷花面积5.03 hm2，约占水域面积的22%。项目组于2016年和2017年分别3 次进行加密水质监测，监测断面在入口处为1#断面，2#—8#断面分别设定3～6 个监测点位不等，合计41 个点位，监测断面及永年洼的进水出水的具体位置见图1。

3 湿地水质监测

图1 永年洼湿地植物分布及监测断面

永年洼入水主要通过两个渠道的闸门控制进入，水源主要来之滏阳河和广府古城的生活污水，滏阳河承纳了邯郸市区的工业及生活废水，永年洼入流水质基本为V 类或劣V 类水体。在入流附近区域，主要种植有荷花和芦苇，为评估湿地的作用，在入流附近设定了6 个监测断面，合计31 个点位，分别监测水质、微生物等指标。2017年9月对2#—6#断面中的DO 进行监测，分别为开阔水面区、荷花种植区和芦苇种植区，监测结果见图2。

图2 永年洼2#—6#断面的监测DO 浓度比较

水体的DO 来自于大气复氧、藻类光合作用及湿地植物的根部泌氧，水体中的耗氧有化学耗氧、生物耗氧及底泥耗氧。永年洼2#—6#监测断面平均溶氧开阔水域为5.53 mg/L，荷花水域为5.76 mg/L，芦苇水域为6.15 mg/L，芦苇水域的溶解氧高于荷花水域，高于无植物的开阔水域。由图2 可知，湿地植物种植后，除了2#断面芦苇水域的DO 低于开阔水域的DO，4#断面的荷花水域DO 低于开阔水域的DO，其他断面上均呈现出有植物水域的DO 高于无植物水域的DO，其中6#芦苇水域DO 远高于无植物水域及荷花水域。2#断面靠近河流的入口，周边有植物种植，被区域水流速加快，大气的复氧水平高。而芦苇种植区域水流缓慢，底部的耗氧水平可能较高，产生DO 差异。

2016年永年洼开始种植水生植物，2017年从面积与生长水平上都较为稳定。将2016年设定为湿地植物种植前情景，2017年设定为湿地植物种植后的情景，用两年的水质数据进行对比。2016年水质监测的月份为5月、8月和10月，2017年水质监测的月份为6月、8月、11月。将断面上各点的水质进行平均，统一为8 个断面。2016年对8 个点面及出口水质进行了监测，2017年为强化示范区监测，仅对6 个断面进行了监测。为比较2016年和2017年的水质变化，统一对前6 个断面的NH3-N、TN 和TP 的平均值进行比较，分别见图3、图4 和图5。

图3 永年洼8 个断面2016、2017年均NH3-N 监测结果比较

2016年和2017年各断面平均水质监测结果表明，从入流至出流的不同断面NH3-N、TN 和TP 浓度呈现出整体降低态势，符合对数衰减模式。2016年断面NH3-N 浓度范围为0.75～1.28 mg/L，平均浓度为0.93 mg/L；2017年断面NH3-N 浓度为0.62～1.25 mg/L，平均浓度为0.77mg/L。2016年断面TN浓度范围为2.78～4.09 mg/L，平均浓度为4.26 mg/L；2017年断面TN 浓度为2.75～4.00 mg/L，平均浓度为3.02 mg/L。2016年断面TP 浓度范围为0.14～0.34 mg/L，平均浓度为0.19mg/ L；2017年断面TP浓度为0.11～0.33 mg/L，平均浓度为0.16 mg/L。水质浓度的数值2017年总体上比2016年下降，在入口处变化不大，主要是2#—6#断面的浓度下降较为明显，2016年NH3-N、TN 和TP 浓度在2#—6#断面上平均浓度下降了20.7%、8.8%和17.2%。湿地植物的种植对水质有一定的改善作用。

图4 永年洼8个断面2016、2017年均TN监测结果比较

图5 永年洼6个断面2016、2017年均TP监测结果比较

4 湿地模型及模拟

4.1 基本方程采用MIKE21 建立水动力模型，耦合EcoLab 模块来对永年洼的水质情况进行模拟。水动力模块采用二维水深平均的浅水方程：

式中： t 为时间，s；x、y 为笛卡尔坐标系；η 为水位，m；h 为水深，m； g 为重力加速度，m/s2；τsx、 τsy、 τbx、 τby为x 和y 方向的表面风应力和底部切应力，植物的作用反应在底部切应力中； Sƒ为源项，kg/（m3·s）；us、 vs为源项水流流速，m/s；Txx、Txy、Tyy为侧向应力。

水质模块通过EcoLab 的自定义及边界功能，参考其他模块中的水质反应方程，建立湿地水质的“WQ+Veg”模块，主要计算及求解变量为DO、BOD、NH3-N、NO3-N、NO2-N 和TP，二次计算变量总氮（TN）为NH3-N、NO3-N 和NO2-N 之和。各水质求解变量通用的方程为：

式中：c 为“WQ+Veg”中计算的状态变量的浓度，mg/L；Dx、Dy为x， y 方向上的扩散系数；Sc为源汇项；Pc为生化反应项。

式中：i 为特定状态变量的进程数；Processi为水体中物质浓度变化的各种生化反应过程。

DO、BOD、NH3-N、NO3-N、NO2-N 和TP 的过程方程如下：

式中：P 为植物的产氧量，mg/（m2·d）；R 为呼吸作用耗氧，mg/（m2·d）。

P 昼夜波动，由午时最大产氧量计算：

式中：Pmax为午时最大产氧量，由植物分布图给出，mg/（m2·d），同样植物呼吸作用耗氧也由空间分布给出；τ为光合作用最大产氧时刻与正午（12：00）的偏差小时数；α为相对昼长；θ1为植物光合作用和呼吸作用与温度相关的阿伦尼斯温度系数；T 为水温，℃。

湿地中生物降解过程通常在好氧环境中进行，微生物数量及活性增强，促进了水体中NH3-N 和磷酸盐（PO4）消减，微生物摄取NH3-N 及TP 的量分别由UM_N 和UM_P 表示，计算采用公式如下：

式中：UNm为微生物吸收氨氮的系数，mgN/mgBOD；UPm为微生物吸收磷的系数，mgP/mgBOD；K3为BOD 降解率，1/d；θ3为阿伦尼斯温度系数； HS-NH3为微生物摄取氨氮的半饱和浓度，mgN/L；HS-TP 为微生物摄取磷的半饱和浓度，mgP/L； HS-Pmax为有机物降解的半饱和由植物泌氧产生的DO 浓度，mgO2/L。

植物对氮磷的摄取量分别为UP_N 和UP_P，计算公式如下：

式中：UNp为植物吸收氨氮的系数，mgN/mgO2；UPP为植物吸收磷的系数，mgP/mgO2； F (N,P)为植物光合作用的营养盐限制函数。

构建一个EcoLab 的“WQ+Veg”模板分以下步骤：（1）先定义计算变量（state variables），即文中计算的6 个变量；（2）设定计算中参数（constants），包括不变的数值和空间上的变化量，例如，区域维度、BOD 降解系数K3、午时最大产氧量Pmax等；（3）设定外部驱动参量（forcing），可以是随时间变化的量，如温度、风等；（4）设定辅助量（auxiliaries），如一定气温下的饱和溶氧浓度、日变化的太阳辐照度因子等；（5）最后根据反应方程，写入计算公式，如式（12）—式（15）等。

4.2 模型参数永年洼模型计算的情景有两个，2016年湿地植物种植前和2017年湿地植物种植后的水质进行比较。2016年不考虑湿地植物的生长，洼内各参数分布均匀。 2017年考虑湿地中种植植物的种类，同时考虑种植的密度，分区设定植物作用下的曼宁数m，（m=1/n，n 为曼宁系数）、午时最大氧气生成量（mg/（m2·d））和植物呼吸率（mg/（m2·d）），见图6。

图6 模拟2017年设定的曼宁数m 及植物午时最大氧气产生量分布

模型计算中微生物摄取氨氮的系数UNm为0.012gN/gBOD，微生物细菌摄取磷的系数UPm为0.008gP/gBOD。植物摄取氨氮的系数UNp 为0.011gN/gDO，植物摄取磷的系数UPp 为0.008gP/gDO。

4.3 湿地水质模拟考虑植物生长周期，模拟时段为3月1日至7月30日，共5 个月，水动力耦合水质进行模拟。永年洼水域体积大，换水周期较长，计算前45 d 为消除水动力及水质的初始条件。根据时间的永年洼出入流及水位数据，2016年及2017年变化不大，入口处流量为1 m3/s，出口处水位按照永年洼2016年3—7月实测水位变化设置。除了出入流，风场是永年洼流场的主要驱动力之一，根据资料用多年平均风速风向设定风场。水质中的生化反应与温度相关，根据监测数据设定水温过程。

湿地植物对水流产生作用，模型入流水质边界根据入口处监测的水质数据给出，出水为出流边界。分别在计算区域中取3 个点，比较有无植物对水流的影响，见表1，其中3 个点位置见图1。

在断面2#的设置开阔水面及植物区的流速比较中可见，种植植物有束窄水流的作用，开阔水面a 点位流速增加，b 点位的流速减小，大气的复氧水平与饱和DO 和流速、水深、风速相关，流速的增加有助于大气的复氧，实测中2#断面的开阔水面DO 较高的部分原因。

模型计算5 个月，计算在4月后，淀区水质浓度变化不大，趋于稳定，取4、5月的计算数据与实测值做比较。取模型计算的从入口至出口的断面，见图1 中点画线，为与实测数据对比，在实测2#—8#断面上，取多个点位计算值平均，模型沿程计算水质数据同实测的入口至出口的1#—8#断面的NH3-N、TN 及TP 浓度比较见图7—图9。

永年洼计算结果显示从入口至出口，NH3-N、TN 和TP 浓度在沿程呈指数形式下降。实测2#—6#断面的点位大多在植物分布区域，计算的沿程结果略高于该区域的数值。2016年及2017年数据进行比较，2017年永年洼断面1#处入流水质浓度略有下降， 2017年在2#—6#断面上模型计算NH3-N、TN和TP 浓度平均下降了14.9%、7.33%和14.8%，实测数据下降比例为20.7%、8.8%和17.2%。有湿地植物2017年情景与无植物的2016年出水口处的NH3-N、TN 和TP 相比，浓度下降了26.4%、11.5%和27.6%。作为华北第三大洼淀，永年洼湿地对水体有较强的净化作用。比较入口与出口的浓度，可看出2016年NH3-N、TN 和TP 去除率为70.6%，46.5%和67.1%，2017年大量种植植物后，去除率为77.9%、51.8%和76.2%，分别提高了7.3%、5.3%和9.1%。植物种植区域大多靠近岸边，水流交换不频繁，湿地植物种植后整体的氮磷去除率有提高，但未超过10%，TN 去除率小于NH3-N。

表1 无植物（2016年）及有植物（2017年）3个点位流速比较

图7 模型2016年和2017年输出入口至出口断面沿程NH3-N 与实测值比较

图8 模型2016年和2017年输出入口至出口断面沿程TN 浓度与实测值比较

图9 模型2016年和2017年输出入口至出口断面沿程TP 浓度与实测值比较

与实测数据相比，模型计算平均误差，NH3-N 在2016年和2017年分别为8.5%和15.9%；TN 在2016年和2017年分别为6.4%和9.3%；TP 在2016年和2017年分别为21.7%和35.2%。TN 模拟较好，TP 模拟值有些差异，可能是因为未对底泥做监测，底泥中释放量未能分区精确给定。计算各指标在永年洼空间分布见图10。

图10 2016年和2017年有无植物NH3-N、TN 及TP 空间分布

模型计算结果表明，构建的EcoLab 湿地植物作用下的“WQ+Veg”模块，其中物理及生化反应的概念清晰，数值计算结果与实测的数据相比较，能反映湿地植物对氮磷的作用，模型的结果可靠。

5 结论

研究区域永年洼湿地生态系统退化、水污染严重，在地形改造，退耕还湿、水文调控的基础上，进行生态修复及水质改善。本文通过原位监测，构建湿地植物模型，模拟研究湿地植物对水体水质的影响。

（1）对永年洼湿地恢复工程建设前后2016年和2017年的水质监测表明，平均条件下永年洼中种植芦苇的水域溶解氧略高于种植荷花的水域溶解氧，高于未种植植物区域。但在个别区域，开阔水体的DO 略高，可能是因为植物种植，束窄了水流，水体流动的增加有利于水面的复氧效果。对永年洼进行水质监测表明，在进水水质差别不大的条件下，随着湿地植物的恢复，NH3-N、TN、TP 等水质指标都有不同程度的下降，在入口水质浓度变化不大的条件下，2#—6#加密监测断面上，NH3-N、TN 和TP 浓度分别下降了20.7%、8.8%和17.2%，且符合对数衰减模式。

（2）根据湿地植物对水质的作用机理，基于湿地植物的泌氧、细菌活性增加和植物吸收作用下水化学反应过程，构建“WQ+Veg”水质模块，耦合MIKE21 的水动力模块进行计算。水质计算中对相关湿地植物参数进行设定，模型模拟了永年洼2016年及2017年人工湿地种植前后湖区流速及NH3-N、TN、TP 等水质指标浓度的变化。与实测结果相比，NH3-N、TN 和TP 的平均误差为12.2%、7.8%和28.4%，模型结果可信度高，永年洼湿地对水体有较强的净化作用。比较入口与出口的浓度， 2016年湿地对NH3-N、TN 和TP 去除率为70.6%、46.5%和67.1%， 2017年有大量湿地植物后，湿地对各指标的去除率为77.9%、51.8%和76.2%，分别提高了7.3%、5.3%和9.1%。植物的种植区域，大多靠近岸边，湖岸与湖心水体交换不频繁，湿地植物提高了湿地对氮磷的去除率，植物作用对NH3-N 的去除效果高于对TN 的去除。构建的“WQ+Veg”模块，物理及生化反应概念清晰，能体现植物的类型、空间分布条件下的湿地水质变化，可为大尺度湿地设计及湿地水质模型提供技术支撑。