河流潜流带氮素迁移转化数值模拟研究进展

蔡 奕，邢婧文，阮西科，周念清，黄若尧，乙东泽

(1.同济大学土木工程学院，上海 200092； 2.长江水环境教育部重点实验室，上海 200092)

随着工农业生产的快速发展，大量含氮生活污水、工业废水以及农业面源污染等通过多种途径进入河流，破坏了原有的营养盐平衡，导致水质下降，从而威胁河流生态系统健康；而且河流中超负荷承载的氮素会随水流向下游输移，进入河口后会诱发近海水体富营养化，严重时会导致大量鱼类和贝类等生物的死亡[1-2]。如何有效降低或去除河流中超标的氮素以恢复生态平衡，已成为当前水环境保护亟须解决的问题[3]。

河流潜流带(以下简称“潜流带”)是河床以下并延伸至两侧岸滩区的水分饱和沉积物含水层，地表水和地下水在此混合发生转化，其内部环境条件复杂多变[4]。潜流带中水体交换引发了碳、氮、磷、氧等物质交换和能量传递，即潜流交换，为潜流带生物群落提供了适宜生存的条件。潜流带生物地球化学反应强烈，对调节河流的生态具有重要作用[5]。潜流带沉积物中碳源物质、氧化还原环境、生物群落等因素对氮素赋存形态和循环模式产生重要影响[6]。因此，要解决河流中氮素污染问题，必须对潜流带氮素迁移转化机理和运移规律进行深入研究。

研究潜流带中氮素迁移转化机理和规律的方法主要有现场监测、室内试验和数值模拟等方法[7-9]。相比于现场监测和室内试验，数值模拟不受点位布设和尺寸效应的影响，可定量描述氮素在潜流带地下水中迁移转化的连续行为，还可对其进行反演和预测，因而得到广泛应用。由于潜流带氮素迁移转化过程复杂，影响因素多，精准量化存在一定的困难[10-11]。本文系统总结了潜流带氮素迁移转化模拟的研究现状与存在的不足，并提出今后的发展方向，可为河流污染防治和水资源保护提供参考。

1 潜流带氮素迁移转化模拟常用数学模型

1.1 潜流带水动力模型

1.2 潜流带溶质运移模型

潜流带溶质运移模拟常用对流-弥散模型。该模型基于连续性方程和质量守恒定律建立，包含对流、弥散及其他源(汇)项，通过水流速度与水动力方程耦合求解[19]。常见的源(汇)项有吸附、解吸等物理作用项及自然降解、微生物催化等生物地球化学反应项[20-21]。

①有氧呼吸；②硝化作用；③反硝化作用；④微生物同化作用；⑤植物同化作用；⑥氨化作用；⑦厌氧氨氧化作用图1 潜流带中氮素循环示意图Fig.1 Schematic diagram of nitrogen cycle in hyporheic zone

1.3 潜流带微生物生长模型

Hampton等[31-32]研究表明，好氧菌、硝化菌、反硝化菌等微生物的生长与死亡会影响潜流带中氮素迁移转化。潜流带中微生物生长是一个自限性的过程：随着微生物量(生物膜)的增长，沉积物孔隙空间逐渐被填充，生物堵塞导致底物扩散通量和沉积物渗透系数降低，使得物质和能量供给下降，造成微生物的死亡，而随着生物堵塞程度的降低，营养物质得以补充，微生物又开始新一轮的生长繁殖[33-34]。Ping等[35-37]建立了微生物生物量与沉积物孔隙度(或渗透系数)之间的定量关系，涉及微生物生长速率、死亡速率、细胞密度、初始孔隙度等参数。将计算获得的微生物量及孔隙度(或渗透系数)代入溶质运移模型和水动力模型，可实现潜流带渗流场、化学场和生物场的信息互馈。

2 耦合模型主要参数及其影响因素

2.1 水动力参数

渗透系数和水动力弥散系数是反映氮素在地下水中对流和弥散运动的重要参数。渗透系数的大小主要受沉积物组成、沉积条件及饱和度等因素的影响。潜流带沉积物是在流水中以机械方式沉积的碎屑物，通常具有沉积韵律性和流水成因的沉积构造特征，分选性较好，成层性较清晰，非均质性特征显著。潜流带既是河水和地下水的交换区，又是大多数无脊椎动物和微生物的生境，所以潜流带沉积物孔隙结构会因水流冲淤、生物行为、颗粒物质溶解与沉淀等作用影响而呈现动态变化[38-39]。此外，河岸潜流带是变饱和区，渗透系数与饱和度和含水量之间存在一定的关系。潜流带的非均质性和动态性导致了沉积物渗透系数的时空差异性，其差异可达到几个数量级[40-41]。由于潜流带水-土-生物环境复杂且多变，渗透系数变化随机性大，目前难以建立定量关系准确表征渗透系数的空间分布特征和动态变化特征。因此，用于模拟计算的渗透系数主要是通过室内试验、野外现场测定或经验估算确定的。作为水文地质的重要参数，渗透系数的取值会直接关系到渗流速度的计算并影响氮素迁移转化模拟结果[42]。可见，渗透系数的准确测定是潜流带氮素迁移转化模拟的重要基础。

水动力弥散系数包括机械弥散系数和分子扩散系数，通常采用室内或现场的弥散试验确定。Godoy等[43]研究表明，室内试验存在尺寸效应，其测定值往往比野外测量值小，两者之间可能有数量级上的差异。大量的研究成果显示影响水动力弥散系数的因素有很多，如孔隙水流速、介质粒径大小、溶质性质、温度、迁移距离等。在低流速条件下，分子扩散作用不可忽略，水动力弥散系数与流速之间并非正比关系[44]。孔隙介质粒径的增大会导致水动力弥散系数的增大[45]。饱和土的水动力弥散系数与孔隙水流速存在显著的相关性，非饱和土的水动力弥散系数与含水量有关且与浓度无关，无论是饱和土还是非饱和土，吸附性溶质的水动力弥散系数通常要大于非吸附性溶质[46]。此外，温度的升高和迁移距离的增大也会导致水动力弥散系数的增大[47]。受河水、地下水和大气温度差异的影响，潜流带的温度场时刻处于变化状态。可见，潜流带的非均质性和动态性会导致水动力弥散系数的时空差异性。在已有的相关数值模拟研究中，水动力弥散系数通常采用定值，其时空和溶质差异性并没有得到很好的体现，这在一定程度上也影响了氮素迁移转化模拟的精度。

2.2 氨氮吸附-解吸行为参数

2.3 氮素生物地球化学反应参数

潜流带中氮素生物地球化学作用受多环境因子控制，包括氧化还原电位、pH值、温度、浓度梯度、水流条件、微生物群落等[53]。Pescimoro等[54-56]研究表明，DO和DOC是潜流带氮素生化反应类型和进程的重要控制因素。在沉积物环境下，反硝化细菌酶在25～35℃范围内活性高，反硝化反应速率快，氮素去除效率高[57]。当氧化还原电位范围处于100～350 mV时，反硝化速率随着氧化还原电位增大而减小[58]。异养反硝化菌适宜的pH值为5.5～8.0，此时的反硝化速率相对较高，pH值小于5.0的酸性环境会抑制反硝化作用[59]。以硝化作用为主的生化反应会随着DO的消耗而转变为以反硝化作用为主，常通过建立DO消耗与河水滞留时间之间的函数关系实现对这两类反应的动态识别[60-61]。可见，与水动力参数和吸附-解吸参数相似，潜流带氮素生化反应参数并非为常数，而是会随着渗流运动和环境变化而发生改变。不过，在目前的模拟研究中并未考虑生化反应参数与环境因子和水流条件之间的关系。

3 数值模拟的不确定性问题

因受观测数据和控制方程的限制，潜流带氮素迁移转化模拟结果往往与真实值存在一定的偏差。建模时应充分考虑信息数据、数学模型及模型参数对模拟结果的影响，以降低过程模拟的不确定性。

3.1 数据信息获取

潜流带是地表以下的黑暗生境，其渗流场、化学场、氧化还原环境、生物菌类活性等信息获取很容易受人为干扰而失真，所以潜流带原生环境信息的准确获取是高精度数值模拟的前提。目前，针对潜流带的研究主要借助于观测井，通过在观测井中放置水位、温度、电导率、pH值、氧化还原电位等传感器或采集水、土样进行分析，以获取数值模拟所需数据[62-63]。然而，潜流带孔隙水水质、细菌分布特征等相关数据会随采样方式或条件的不同而发生变化。例如，观测井与大气连通会影响水中DO含量，可能会导致原先的还原环境转变为氧化环境，使得氮素存在形态发生改变；氧气、光照、温度等采样条件控制不当会造成沉积物土样中微生物生物量测定失真；采集频率低可能会导致多场信息变化幅度难以捕捉。因此，需要结合高分辨率高灵敏度的原位传感监测、非扰动测试及情景模拟试验等技术获取原生环境信息数据，这对氮素迁移转化模拟研究具有重要意义。

3.2 耦合模型构建

潜流带氮素迁移转化行为十分复杂，是渗流场、温度场、化学场、生物场等多场相互作用的结果，故其数值模拟的关键在于描述多场时空演变过程及其相互之间的耦合联动关系。目前的数值模拟主要关注水动力、溶质运移、微生物生长等过程的耦合计算，其中溶质运移计算通常考虑对流项、弥散项、反应项。由于土颗粒对氨氮具有一定的吸附作用，氨氮运移模型应增加吸附项，否则获得的氨氮时空分布规律与真实情况不符。朱静思等[64-65]研究表明，温度变化会影响氮素的水动力弥散系数、吸附速率和平衡常数及微生物催化反应速率等，进而控制氮素的循环过程。潜流带温度场常伴随着水量交换发生显著的动态变化[66-67]，但在现有的耦合模型中却极少考虑温度的时空演变。为了对潜流带氮素迁移转化有更全面的认识，有必要建立模型参数与温度、pH值等环境因子之间的定量关系，通过水动力、热传导、溶质运移及微生物生长模型的耦合计算，实现多过程之间的动态联动，如图2所示。此外，目前大多数的耦合模型采用一、二维数学模型[68-69]，然而，潜流带具有复杂的三维空间非均质各向异性结构，采用三维模型更符合潜流带氮素迁移转化过程的定量化表达。

图2 潜流带氮素迁移转化模型中不同过程耦合示意图Fig.2 Schematic diagram of coupling different processes involved in the model of nitrogen migration and transformation in hyporheic zones of rivers

3.3 模型参数识别

潜流带氮素迁移转化过程耦合模型结构复杂，涉及参数多，如渗透系数、水动力弥散系数、反应物最大生化反应速率等，其中一些参数在维度、位置或时间上有差异，参数不确定问题是客观存在的。基于经验估算或观测数据优化所获取的参数通常不能保证模型在各种可能情景下的模拟精度。而且，因不同算法在参数收敛轨迹上存在差异，优化结果往往也不尽相同，甚至无法判断其是否达到全局最优[70]。鉴于优化方法在复杂模型的参数识别上有一定的局限性，参数敏感性研究越来越受到关注。目前，潜流带氮素迁移转化模型的参数敏感性研究主要采用基于“最佳”估计值的扰动分析方法，其研究结果表明敏感参数主要为渗透系数、水动力弥散系数、最大生化反应速率及半饱和常数等[71-73]。此外，多参数识别通常要比单参数要复杂得多，这是因为复杂模型结构中不同参数之间可能会存在某种相关性，单独考虑各参数的“最佳”估计值未必达到最优的模型效率[74]。可见，潜流带氮素迁移转化模型的参数识别需要考虑参数时空分布的复杂性、灵敏度及各参数之间的相关性，而且合理的算法应充分体现参数不确定性而不是参数最优化。然而，这些参数识别中所要关注的问题现有研究成果还较少，还有待进一步研究。

4 结论与展望

b.深入研究潜流带氮素迁移转化中多过程间的耦合联动关系。通过试验研究探讨潜流带渗流场、温度场、化学场、生物场的时空演变规律，厘清潜流带水动力、热传导、溶质迁移及微生物活动等多过程间的耦合联动关系，建立主要模型参数与温度、pH值等环境因子间的定量关系，在此基础上对现有耦合模型进一步完善，如耦合热传导过程，考虑氨氮吸附-解吸行为和其他氮素生化反应等。

c.开展多维度多过程耦合模型优化计算研究。深入研究潜流带氮素迁移转化模型参数的复杂性、敏感性及相关性，提出基于不确定性分析的参数识别算法。分析不同过程的持续时间和涉及范围，合理设置各过程的时间步长和计算网格，以实现多过程间的同步与异步模拟。此外，三维的多过程耦合会消耗大量计算资源，需加强模拟软件的开发与集成，优化空间网格，强化并行计算，节省运行时间，以提升模拟精度和效率。