湖泊富营养化是水环境领域面临的突出问题^[1-4].为有效预防、控制和治理湖泊富营养化问题，国内外围绕其成因、危害以及控制和管理等方面开展了大量的研究^[5-10]，其中成因分析是有效治理的基础.研究表明，氮(N)、磷(P)是湖泊生态系统必需的营养元素，正常的氮、磷平衡是维持湖泊健康的必要条件；但大量氮、磷的流入和累积会改变湖泊内营养盐的通量平衡，增加富营养化的风险^[11-13].因此对湖泊恢复而言，开展湖泊氮、磷营养盐循环和通量平衡研究就显得尤为重要^{[9, 14-15]}.具体而言，氮、磷的过量输入会导致湖泊营养盐循环发生改变，藻类大量生长、消亡，影响水体的理化条件并进一步改变氮、磷的内部循环过程，增加湖泊富营养化的风险^[16].研究表明，随着湖泊外源污染负荷的削减，湖泊底质内源释放对湖泊富营养化的贡献也逐渐突显出来^[17].氮、磷的内源释放和内部循环会使湖泊在大量削减外源负荷后，也继续维持长时间的富营养化状态，阻碍湖泊修复.研究发现，对于部分富营养化的湖泊，内源负荷可以占到总入湖负荷量的50%以上^[18].同时，一些富营养化湖泊的修复过程表明：在有效控制外源负荷的前提下，仍需要10~20年才能使湖泊的生态功能得到明显的恢复^[19].因此，通过对湖泊氮、磷循环通量的追踪，可分析湖体中所发生的氮、磷迁移转化过程，识别不同氮、磷循环过程的响应机制，加深对富营养化过程及机理的认知.机理模型是研究湖泊富营养化过程和营养盐通量平衡的理想手段^[20-21]，通过对水体中与富营养化相关机理过程的模拟，为制定更为细致和准确的治理方案提供科学依据.在模型的应用中，简单机理模型过程简单，无法实现对重要营养盐过程的刻画；而复杂的三维水动力-水质模型研究则大多关注于模型的性能及其不确定性，对氮、磷循环过程的关注较少^[22].

本文基于以EFDC为内核的三维水动力-水质模型(IWIND-LR)，开发了利用数值积分的氮、磷营养盐通量与存量计算方法，并以滇池外海为例开展应用研究.滇池(24°40′~25°02′N, 102°36′~102°47′E)是我国重点治理的“三湖”之一，湖面面积为306.3 km²，湖体体积为15.6×10⁸ m³，平均水深4.4 m，流域面积2920 km^2[23].由于流域社会经济的快速发展，滇池从1980s开始出现富营养化现象，目前呈重度富营养状态；外海水质在V~劣V类之间波动，且呈现北部偏高、南部偏低的趋势^[24].本文基于前期经过校正和验证的滇池三维水动力-水质模型^[23]，在此基础上重点开展氮、磷的内部循环过程通量研究，分析滇池水质和藻类生物量的变化过程，识别不同氮、磷循环过程对湖泊富营养化的贡献，为滇池富营养化治理提供科学依据.

1 研究方法 1.1 模型构建

模型构建基于前期经过校正和验证的滇池三维水动力-水质模型^[23]，主要包含三维水动力和温度动力学模拟以及内部耦合的水质模型系统；可模拟26个水质状态变量，并可同时模拟底泥与水体间的营养盐交换动力学以及底泥对溶解氧的动态作用，动态跟踪水龄及低溶解氧区的时空分布.其核心的水质模块主要基于物质守恒，刻画物理输送、大气交换、吸附与解吸、藻类吸收、底床-水界面交换、硝化与反硝化以及沉积成岩作用等，基本控制方程如下：

$\begin{array}{l} \frac{{\partial C}}{{\partial t}} + \frac{{\partial \left( {uC} \right)}}{{\partial x}} + \frac{{\partial \left( {vC} \right)}}{{\partial y}} + \frac{{\partial \left( {wC} \right)}}{{\partial z}} = \frac{\partial }{{\partial x}}\left( {{K_x}\frac{{\partial C}}{{\partial x}}} \right) + \frac{\partial }{{\partial y}}\\ \left( {{K_y}\frac{{\partial C}}{{\partial y}}} \right) + \frac{\partial }{{\partial z}}\left( {{K_z}\frac{{\partial C}}{{\partial z}}} \right) + {S_c} \end{array}$

(1)

式中, C为水质状态变量浓度, u、v、w分别为x、y、z方向的速度矢量, K_x、K_y、K_z分别为x、y、z方向的湍流扩散系数, S_c为每单位体积内部和外部的源和汇.在一阶动力学中，存在以下关系：

${S_c} = k \cdot C + R$

(2)

式中, k为动力学速率，R为由于外部负荷和内部反应引起的源汇项.

在前期的滇池三维水动力-水质模型^[23]构建中，将滇池外海水体划分为垂向6层，每层664个水平网格(尺度为634±148 m)(图 1)，主要模拟的状态变量为水温、溶解氧、藻类中氮磷含量、颗粒态有机氮磷浓度、溶解态有机氮磷浓度、磷酸盐浓度、硝态氮浓度、氨氮浓度等.模型的边界条件主要包括入湖支流的流量、温度、水质、沉积床初始条件以及与时间相关的气象条件和大气边界数据^[23]，其中大气湿沉降条件为磷酸盐浓度0.04 mg/L，硝态氮浓度0.70 mg/L.模型底泥模块初始条件为：活性有机磷、有机氮和有机碳底泥体积浓度分别为23±15.3、103.2±13.8和741±176.5 g/m³.

1.2 营养盐存量与通量核算

全湖氮、磷营养盐各组分的存量和通量计算模块是通过与三维水质模型方程耦合，并与水动力-水质模型采用相同离散化分辨率进行数值积分得到的.各组分存量考虑水质组分的时空分异性，各相关状态变量的存量通过对全湖所有网格进行积分获取：

$M\left( {i,t} \right) = {\rm{d}}x{\rm{d}}y{\rm{d}}z$

(3)

式中，M为i组分的存量，c为湖体中i组分的浓度.氮元素各组分存量考虑：藻类活体氮存量(algaeN)、颗粒有机氮存量(PON)、溶解态有机氮存量(DON)、氨氮存量(NH₃)、硝酸盐存量(NO₃^-)；磷元素各组分存量考虑：藻类活体磷存量(algaeP)、颗粒有机磷存量(POP)、溶解态有机磷存量(DOP)、正磷酸盐存量(PO₄).

各氮、磷循环过程通量通过对全湖所有网格进行积分获取，氮循环各个过程通量考虑以下过程：①流域输入(包括藻类活体氮入湖量、颗粒有机氮入湖量、溶解态有机氮入湖量、氨氮入湖量、硝酸盐/亚硝酸盐氮入湖量)；②大气沉降(包括硝酸盐氮沉降)；③出流(包括藻类活体氮出湖量、颗粒有机氮出湖量、溶解态有机氮出湖量、氨氮出湖量、硝酸盐/亚硝酸盐出湖量)；④内源(包括氨氮底泥交换通量、硝酸盐氮底泥交换通量)；⑤沉降(包括藻类活体氮沉降量和颗粒有机氮沉降量)；⑥反硝化过程和固氮.磷循环各个过程通量考虑以下过程：①流域输入(包括藻类活体磷入湖量、颗粒有机磷入湖量、溶解态有机磷入湖量、正磷酸盐入湖量)；②大气沉降(包括正磷酸盐沉降组分)；③出流(包括藻类活体磷出湖量、颗粒有机磷出湖量、溶解态有机磷出湖量、正磷酸盐出湖量)；④内源(包括正磷酸盐底泥交换通量)；⑤沉降(包括藻类活体磷沉降量、颗粒有机磷沉降量、正磷酸盐沉降量).上述过程的主要通量计算方程如下：

流域输入通量：

${F_{{\rm{WS}}(i)}} = \sum\limits_{t = {t_0}}^T {\sum\limits_{k = 1}^N {Q\left( {t,k} \right)} } \cdot C\left( {t,k,i} \right)$

(4)

累积大气沉降通量：

${{F}_{\text{AIR}(i)}}=\iiint{({{D}_{\text{d}}}+{{D}_{\text{w}}})}\text{d}x\text{d}y\text{d}t$

(5)

累积出流通量：

${F_{{\rm{O}}(i)}} = \sum\limits_{k = 1}^M {\int q } \left( {k,t} \right)c\left( {k,t} \right){\rm{d}}t$

(6)

累积底泥-水体通量：

${{F}_{\text{Ben}(i)}}=\iiint{{{F}_{\text{b}(i)}}}\text{d}x\text{d}y\text{d}t$

(7)

累积沉降通量：

${{F}_{\text{Set}(i)}}=\iiint{{{F}_{\text{s}(i)}}}\text{d}x\text{d}y\text{d}t$

(8)

累积反硝化通量：

${{F}_{\text{Den}}}=\iiint{\int{{{F}_{\text{dn}}}}}\text{d}x\text{d}y\text{d}z\text{d}t$

(9)

累积固氮量：

${{F}_{\text{Nfix}}}=\iiint{\int{{{F}_{\text{Nfix}}}}}\text{d}x\text{d}y\text{d}z\text{d}t$

(10)

式中, Q为入湖流量(m³/d)，C为每条入湖河流中i组分的浓度(g/m³)，D_d为干沉降(g/(m²·d))，D_w为湿沉降(g/(m²·d))，q为出流流量(m³/d)，c为出流中i组分的浓度(g/m³)，F_b为i组分的底泥释放速率((g/(m²·d))，F_s为i组分的沉降速率(g/(m²·d))，F_dn为反硝化速率(g/(m³·d))，F_Nfix为固氮速率(g/(m³·d)).需要注意的是，根据以上方程计算的氮、磷循环通量结果并不是瞬时通量结果，而是累积通量结果.为计算湖泊氮、磷循环通量和存量，模型以相同条件连续运行5年以达到相对稳定状态，取第5年结果计算湖泊氮、磷循环通量及存量.

2 研究结果 2.1 模型校验

模型校验是水质模型构建的核心步骤，目的是推算合理的参数化方案以实现模型对实际湖体过程的再现^[25].滇池模型的校验分2个步骤：水动力模型的校验、水质模型校验；其中水动力模型以水位和温度作为评判标准，水质模型校验以滇池外海8个国家常规监测点的月监测数据，主要校准溶解氧(DO)、总氮(TN)、总磷(TP)、NH₃-N、叶绿素a(Chl.a)浓度等.由于本文采用的是前期经过校正和验证的模型，限于篇幅，其校正和验证在此不再赘述，具体结果请参阅相关文献^[23].由校正和验证结果可知^[23]，TN、TP和Chl.a的模拟结果与实测值吻合较好，RMSE分别为0.619、0.064和0.040 mg/L，模型能够较为真实反映湖体物理、化学和生物过程的变化.校核后的模型参数取值如下：水体最大硝化反应速率为0.05 d^-1，参照温度是25℃，温度效应系数为0.045；反硝化速率为0.06 d^-1；有机氮矿化速率为0.05 d^-1，有机磷矿化速率为0.04 d^-1，有机碳为0.06 d^-1，矿化速率温度修正系数为0.069；颗粒营养盐沉降速率为0.1 m/d；蓝藻、硅藻和绿藻生长速率分别为1.35、2.5和1.9 d^-1.底泥通量采用底泥成岩模块计算，主要参数为：G1类有机物降解速率为0.015 d^-1，G2类为0.001 d^-1，G3类为0，温度修正系数为1.09；底泥硝化最佳反应速率为0.09 d^-1.

2.2 滇池营养盐通量存量平衡核算

在模型校准的基础上，对滇池不同形态氮、磷营养盐的存量及主要源汇过程的通量进行计算(图 2).由结果可知，氨氮、硝态氮和有机氮之间存在明显的相互转化关系和季节特征；磷酸盐是磷元素各个组分中存量最大的组分，呈现明显的夏季增加、冬季减少的趋势.氮、磷的底泥释放也存在明显的季节变化，其释放通量都表现为夏季较大、冬季较低的规律，其中氨氮和磷酸盐在冬季为负值，即底泥的吸收过程.反硝化的季节特征与其他过程不同，7月之前的反硝化通量明显小于之后的反硝化通量.

2.3 关键源汇过程识别

表 1汇总了2003年滇池湖体氮、磷循环各个过程通量的年均值，以TN和TP计，滇池TN的年总输入量为7620.92 t，总输出量为7637.31 t；TP的年总输入量为450.23 t，总输出量为429.57 t.对于氮而言，贡献最大的源过程是陆域流入过程，占总流入的67.76%，其次是底泥释放过程，氨氮和硝态氮共占29.04%，大气沉降和固氮作用共占3.2%.氮的去除过程则主要为反硝化和沉降，共占总去除量的95.5%，其中反硝化48.99%，沉降为46.51%.对于磷而言，陆域流入过程占总流入量的79.77%，其次是磷酸盐底泥释放作用，为18.3%.磷去除则主要为沉降作用，其中颗粒态磷沉降占62.77%，藻类沉降占27.87%.

3 讨论 3.1 数值模拟与基于监测的通量核算结果的差异

首先需要明确的是本文中的营养盐质量平衡关系是基于三维复杂水动力-水质模型的模拟结果，是对全湖所有网格进行数值积分所得到的.这不同于一般的营养盐质量平衡分析方法，一般意义上的营养盐质量分析方法是基于有限的监测数据计算得到的，在空间上并不能做到与三维模型模拟结果一样细致，这也就导致2种方法所计算的营养盐质量平衡会有一定的差异^[22].虽然营养盐质量平衡分析在更大的程度上是一种估算，并不要求观测数据的空间水平和垂直分布，但对于面积较大或深度较大的湖泊而言，会造成结果的明显偏差.而如果借助于三维复杂水质模型，在监测数据的基础上，对全湖各个网格进行积分后再计算营养盐质量平衡关系就会避免不必要的偏差.以滇池为例，根据对8个湖体常规监测点位的数据，计算滇池TP存量为231.4 t，TN存量为3076 t；而通过水动力-水质模型计算可知，滇池的TP和TN存量分别为350.2和2525.4 t，二者之间均有较大的差异，计算结果TP偏高、TN略低.由国家标准监测方法可知，所检测的样品均取自水面下0.5 m，其监测数据主要表征的是表层水体的水质状态.对于滇池而言，特别是磷元素，其主要来源之一是底泥中磷酸盐的释放，这会造成湖体磷酸盐形成自下而上的浓度梯度^[26-28].这种湖体垂向浓度梯度的存在会导致在只用表层水质数据估算营养盐通量平衡时产生系统性误差(图 3).除此之外，考虑到湖泊复杂的水动力条件以及监测点位的分布并不能完全代表全湖的空间差异，这也会在估算时造成一定的误差.相比而言，基于三维水质模型的营养盐通量平衡计算则会涉及到全湖所有已划分的网格，在模型已校准的情况下，比传统的估算方法更为精确，且能提供更为详细的通量结果.

3.2 滇池氮、磷营养盐存量变化

对湖体内氮、磷元素的存量分析可知，不同形态的氮或磷在时间和通量上的变化趋势各异，具体表现在：

① 氨氮和硝态氮的变化趋势不同步.在春季，溶解氧浓度相对较高，随着氨氮浓度的降低，硝态氮浓度经历了先升高后降低的过程，这主要是因为在有氧环境下氨氮通过硝化作用可以生成硝态氮.而夏季由于藻类的大量生长甚至暴发，水体中溶解氧浓度降低，硝态氮会通过反硝化作用生成氮气逸出湖体^[29-30].在7月中下旬藻类暴发后，反硝化作用显著增强(图 2e)，同时缺氧环境也抑制了硝化作用，促进了底泥释放氨氮的过程，所以硝态氮在夏季保持较低浓度，而氨氮由于底泥释放作用和藻类吸收的共同作用会出现较大的波动.秋、冬季则由于藻类生物量降低，水体溶解氧浓度升高会使水体内氨氮和硝态氮浓度逐渐升高.

② 水体中颗粒态有机氮和溶解态有机氮浓度比例在1:5左右，且变化趋势大致相同，这符合湖沼学基本定律中的描述^[31].但是有机氮和氨氮随时间的变化趋势差异非常明显，这表明有机氮的矿化过程受水体DO、pH和微生物活性的影响较大，有明显的季节变化，主要表现为春、冬季矿化速率较高，夏季矿化速率较低.夏季氨氮的来源主要为外源输入和内源底泥释放，多种因素的共同作用导致了滇池水体氨氮浓度随季节的复杂变化趋势^[32-33].

3.3 氮、磷内源释放的贡献及变化

湖体内氮、磷底泥释放和沉降通量是研究关注的重点.氨氮和磷酸盐的底泥释放过程有明显的季节特征，夏季表现为释放，冬季则会表现为吸收.造成这种现象的主要因素除风力扰动外，最为主要的是底泥氧化还原条件的改变.在缺氧条件下，底泥对氨氮和磷酸盐的吸附效率降低，从而促进了氨氮和磷酸盐的底泥释放作用^{[17, 27, 34]}.综合考虑颗粒态有机氮的沉降和氨氮、硝态氮的底泥释放作用可得出，氮的水体-底泥界面交换过程对氮循环总体的贡献表现为释放，释放量为237.6 t/a.而磷的水体-底泥界面交换过程与氮不同，该过程对磷循环总体的贡献表现为沉降，沉降量为173.9 t/a.但需要注意的是，底泥磷酸盐的释放随季节变化非常大，在藻类暴发期间表现为强烈的释放作用，日最大值已达到24.2 t/d，是湖体中磷酸盐的主要来源之一^[35-36].对比之前针对滇池底泥氮、磷内源释放的研究发现，利用模拟实验的方法计算氨氮底泥释放量约为1394.1~6542.2 t/a，磷酸盐的底泥释放量约为74.9~380.2 t/a；本文模型结果中氮、磷的底泥释放量分别为2213.05和82.39 t/a，模型结果与文献中的实验结果范围基本吻合^[37-38].氮、磷元素间的上述差异主要与其化学性质相关，磷循环属于沉积态循环，其沉降到湖底的磷会经过沉积成岩过程被永久的埋藏在湖底；而氮的化学性质相对活泼，其沉降到湖底的主要为有机残体，经微生物分解后会主要以氨氮的形式再次释放到水体，使湖体氮浓度维持在较高的水平.通过对滇池氮、磷元素化学性质以及关键源汇过程通量的分析，可以为滇池富营养化的治理提供理论和决策支持.

4 结论

以IWIND-LR为计算平台，核算了滇池外海氮、磷循环通量及其在氮、磷循环整体中的贡献.研究发现，通过三维数值模型计算湖泊营养盐存量和关键过程通量以及质量平衡关系相较于传统的质量平衡方法更为详细，且可计算反硝化、底泥释放等很难进行实验测定的过程.其中，对于2003年滇池TN的年总输入量为7620.92 t，总输出量为7637.31 t；TP的年总输入量为450.23 t，总输出量为429.57 t.沉降和反硝化则是主要的氮、磷去除途径，反硝化占总氮去除量的48.99%，沉降占总磷去除量的62.77%.必须注意的是，上述的结果并不能推广为滇池长期的趋势.今后的研究将集中在2个方面：①通过实验手段，对模型所涉及的关键过程参数进行测定和核实，减少模型结果的不确定性和可能存在的异参同效；②对滇池进行更长期的模拟，分析不同水文气象条件和流域负荷条件下的营养盐平衡和通量，为滇池富营养化控制提供更扎实的理论依据和技术支持.