有机碳分布对氮素迁移转化的影响数值模拟

刘 琨，朱 焱，王丽影，程先军，杨金忠，史良胜

(1.武汉大学 水资源与水电工程科学国家重点实验室，武汉 430072；2.中水珠江规划勘测设计有限公司，广州 510610；3.中国水利水电科学研究院水利研究所，北京 100048)

1 研究背景

灌溉污水具有供给稳定，氮/磷营养物质丰富等诸多优点，在世界上很多国家得到了广泛应用[1,2]，对缓解农业水资源紧缺发挥了重要作用。但长期不合理的污水灌溉也会造成污水中富集的氮素经运移进入地下水和地表水，增加地下水和地表水污染风险。土壤水分运动受降雨/灌溉入渗、土壤蒸发、作物蒸腾等多因素影响，为复杂的非线性问题[3]。而氮素转化过程则更为复杂，一方面受到土壤水分运动的影响，另一方面，又涉及不同氮素形态之间的物理化学反应过程，影响因素众多。因此精确衡量氮素各迁移转化过程十分困难，而数值模拟则是一种十分有效的研究工具[4-8]。

经过几十年的发展，国内外关于土壤氮素转化和运移的数值模型已有数十种，如SoilN、ArcNLET、DRAINMOD-N、LEACHMN、DAISY、NITROGEN-2D、SPWS、ANIMO、RZWQM等[9-11]。土壤氮素模型主要考虑点尺度或田间尺度水氮动态过程，一般包含土壤水分模块和氮素迁移转化部分，并考虑土壤水分、温度、pH等因素对氮素转化过程的影响。但对有机碳含量的影响关注则较少。研究表明，有机碳含量对土壤中有机氮与无机氮之间的矿化/固持、无机氮的硝化、反硝化等反应速率影响十分明显[12-14]。因此，在污水灌溉等含有外源性有机碳输入情况下，土壤本底有机碳分布将发生改变，需要明确不同有机碳分布对氮素转化过程的影响,有必要对有机碳影响下的氮素迁移转化数值模型进行改进。本文基于田间尺度土壤水分运动和氮素迁移转化数值模型NPTTM，考虑外源输入性有机碳在土壤剖面的分布，增加有机碳对氮素反应过程的影响函数，对模型氮素转化过程影响因素进行改进，并将该改进的模型应用于含有机碳的再生水灌溉实验土壤水氮动态模拟。本文进一步采用该验证算例进行数值实验，分析了四种有机碳分布对土壤氮素动态变化过程的影响。

2 模型原理

2.1 水流运动模型

本文模型的水流运动过程采用Richards方程进行描述。

(1)

式中：θ为土壤体积含水率，[L3L-3]；h为土壤水分负压，[L]；S为根系吸水项或其他源汇项，[T-1]；xi为空间坐标；t是时间,[T]；KAij为各向异性无量纲张量的分量；K是非饱和水力传导度，可由下式表示：

K(h,x,z)=Ks(x,z)Kr(h,x,z)

(2)

式中：Ks为土壤饱和渗透系数,[L T-1]；Kr为相对非饱和水力传导度。

本文采用van Genuchten方程的改进模型来表示土壤体积含水率θ与土壤负压h的关系[15]：

(3)

式中：θr为残余体积含水率(即最大分子持水率)；θs为饱和体积含水率；θa和θm是土壤含水率和土壤负压关系曲线上两个假定值，且置θa=θr，而α和n都是经验常数。

2.2 氮素转化模型

模型中的氮素主要分为有机氮和无机氮。其中，有机氮根据有机物分解速率的不同，又分为快速反应有机氮(有机物、新鲜绿肥或者作物残枝落叶中的氮素)和慢速反应有机氮(腐殖质中的氮素)；无机形态的氮素主要为硝态氮和铵态氮，一般将土壤中的亚硝态氮作为反应中间的过渡态，模型中不单独考虑。本文模型考虑的土壤氮素形态及其主要转化过程如图1所示，主要氮素转化过程包括有机氮向无机氮(铵态氮)的转化，无机氮向有机氮的固持，铵态氮向硝态氮的硝化，铵态氮的挥发，硝态氮的反硝化、根系吸收及深层渗漏等。

图1 本文模型考虑的土壤氮素迁移转化过程

(1)有机氮的矿化/固持过程。关于有机氮的转化，模型有两个假设：快速反应有机物中氮素的流动方向和转化速率由该有机物中碳素的流动方向和转化速率以及土壤中有机物部分的碳氮比决定，慢速反应有机物的分解产物只有铵态氮和二氧化碳，不再转化为其他形式有机物；假设微生物中的碳氮比与土壤慢速反应有机物的碳氮比相同。

有机氮转化的动力学表达式如下：

(5)

(6)

式中：CC,l为快速反应有机物中碳的含量,[ML-3]；Kl为速率常数,[T-1]；fT、fθ为土壤温度和土壤含水量对分解速率的影响函数；CN,l为土壤快速反应有机物中氮素的含量,[ML-3]；rl是土壤中快速反应有机物中的碳氮比；CN为土壤中无机氮的含量,[ML-3]；fe为综合效率因子(快速反应有机物的分解过程中转化成微生物和腐殖质中的碳素的比例)；r0为微生物作用的合成物质与慢速反应有机物中的碳氮比。

模拟中铵态氮的硝化、挥发等过程采用的是一阶动力反应方程，影响函数和有机氮转化过程相同，考虑了土壤温度和土壤含水量的影响[16]。

(2)改进的反硝化反应速率。反硝化过程是指土壤中的硝态氮因为反硝化反应而变成气体进入大气的过程，主要受氧气、水分、有机碳含量等因素的影响。本文模拟采用一阶动力反应方程，同时加入了有机碳的影响函数，改进后的反硝化反应方程如下：

(7)

式中：CN,3是土壤溶液中硝态氮的浓度,[ML-3]；Kdn是一阶反硝化系数,[T-1]；ρ为土壤干容重,[ML-3]；Kd是土壤对溶质的吸附系数,[M-1L3]；fC,l为有机碳影响函数；fθ、fT分别为含水量和温度对反硝化速率的影响函数[8]。

本文模型考虑有机碳对反硝化的直接影响，采用如下公式计算有机碳的影响函数：

fC,l=kC,lCC,l

(8)

式中：KC,l为有机碳影响函数系数,[M-1L3]。

2.3 氮素运移模型

因土壤对铵态氮的吸附作用较强，铵态氮的运移方程可用如下方程表示：

(9)

式中：CN,4为土壤溶液中铵态氮浓度[ML-3]；θDij为饱和/非饱和水动力弥散系数，[L2T-1]；s为吸附在土壤颗粒上的铵态氮浓度，[MM-1]，采用等温吸附模型的形式，即s=KdCN,4；R4为铵态氮各种源汇项之和，[ML-3]，具体包括铵态氮的矿化/固持量、铵态氮的硝化量、根系吸收量和挥发量。

硝态氮一般不易被土壤吸附，因此在硝态氮的运移方程中不考虑土壤的吸附作用，采用如下方程表示：

(10)

式中：CN,3为土壤溶液中铵态氮浓度,[ML-3]；R3为硝态氮各种源汇项之和,[ML-3]；具体包括硝态氮的固持量、铵态氮的硝化量、硝态氮的反硝化量和根系吸收量。

2.4 温度模型

模型中采用如下经验公式计算土壤剖面各深度的温度随时间的变化：

T(z,t)=Ta+A0exp(-z/Dn) cos(ωt+φ-z/Dm) (11)

式中：T为土壤温度，℃；Ta为年平均温度,℃；A0为温度的年变幅，℃；z为计算点埋深,[L]；Dm为衰减深度,[L]；ω为温度波的频率,[T-1]；t为计算时间,[T]；φ为相位移。

3 含有机碳污水灌溉土壤水氮动态过程数值模拟

本次试验为土柱模拟再生水灌溉过程，试验所用再生水取自北京市黄村污水处理厂。试验共设置3个平行试验，在整个试验期间，3个土柱独立运行，互不干扰[17]。图2为试验装置示意图，试验所用有机玻璃柱高120 cm，直径为30 cm，距玻璃柱底部3 cm处开有直径为2 cm的排水孔，并且在排水孔中安装有可控制排水的球阀。试验土壤取自田间，混合均匀，取土时沿垂直方向每10 cm测量一次田间土壤容重并记录。试验期间，土柱放置在可移动的塑料大棚内，在有降雨时用塑料大棚进行挡雨，以免破坏试验条件；在气温较低时，用大棚进行保温，其余时间撤去大棚。灌水周期为7～10 d，实际灌水量和灌水时间根据当时耗水情况确定。每个灌水周期中，在灌水24 h后提取10、40和70 cm三个深度处的土壤溶液作为观测；在灌水72 h后打开排水阀门开始排水。灌水次数及每次灌水量、灌水中硝态氮、铵态氮、有机碳浓度及蒸发量见图3。

图2 试验装置示意图

图3 每日灌溉量、蒸发量及灌溉水中有机碳、硝态氮和铵态氮浓度

图4 温度拟合结果

通过对比累计排水量得到土壤水分运动参数，通过土壤剖面硝态氮和铵态氮浓度对比得到氮素迁移转化参数。

(1)累计排水量拟合。图5为模拟时段累计排水量模拟值与实测值对比，由图5中可知，土柱底部排水累计量实测值和模型模拟值吻合程度较高，可认为模拟的水流过程与实际相符。此时采用的土壤水分运动参数为：θs=0.42，α=4.4 m-1，n=1.23，Ks= 1.3 m/d，θr=0.02。

图5 计算过程中累计排水量模拟值与实测值对比

(2)铵态氮、硝态氮各观测时间剖面上的浓度分布。图6和图7分别表示不同观测时刻，土壤剖面铵态氮、硝态氮浓度模拟值与实测值对比结果。从图中可以看出，各观测时间铵态氮模拟值与实测值吻合较好；硝态氮模拟值在表层10 cm位置与实测值偏离较大，源于表层的硝化作用较为强烈，导致模型模拟值偏大，硝态氮模拟值在40cm以下位置与观测值一致性较强。此时采用的土壤氮素运移参数为：土壤容重ρ=1.57×103kg/m3，纵向弥散度αL=0.07 m，横向弥散度αT=0.01 m。土壤对溶质的吸附系数(铵态氮吸附系数)Kd=4.0 m3/kg，一阶快速反应有机物集合体速率常数Kl=0.004 5 d-1，一阶慢速反应有机物集合体速率常数Kh=0.000 06 d-1，慢速反应有机物的碳氮比r0=10，一阶反硝化速率常数Kdn=0.004 d-1，一阶硝化速率常数Kn=0.06 d-1。

图6 各观测时间土壤剖面铵态氮浓度分布模拟值与实测值对比

图7 各观测时间土壤剖面硝态氮浓度分布模拟值与实测值对比

(3)不同深度处铵态氮、硝态氮浓度与实测值随时间变化对比结果。图8为10、40和70 cm深处铵态氮浓度模拟值与实测值对比。从图8中可以看出，不同深度铵态氮浓度随时间变化趋势与试验值基本一致。埋深10 cm的铵态氮模拟值在26～42 d内出现了明显的峰值，这主要是因为这三次的灌水量较大，且第26、第32、第42 d灌水中铵态氮浓度很高，且其中所含的铵态氮浓度明显高于试验期间灌水中铵态氮浓度的平均值。灌水中的铵态氮含量对土壤表层的铵态氮浓度影响明显，但由于铵态氮吸附作用较强，对土壤20 cm深度以下铵态氮浓度影响较小。由于受到灌水中外源性铵态氮的影响，表层10 cm处模拟值出现明显波动，与灌溉过程一致。而在40 cm及以下各层中，模拟值变化幅度较小，这与模拟中设置的铵态氮吸附系数较大不易移动有关，而铵态氮试验值波动较大，这可能是由于铵态氮本身较易被氧化，而反应过程同时受到多种因素的影响。

图9为模拟时段10、40和70 cm深处硝态氮浓度模拟值与实测值对比。由图9可以看出，表层10 cm处硝态氮浓度仍呈现明显的峰值，且与灌水过程一致，导致此现象的原因是外源性硝态氮的输入，同时铵态氮的硝化作用较强。40 cm埋深以下硝态氮模拟值与实测值吻合较好。土壤10 cm深与40 cm深处，在第26～50 d，模拟值明显比其他时刻值大，导致此现象的原因是：由于本文模型考虑了有机碳对反硝化过程的影响，因此，土壤剖面不同深度处硝态氮浓度变化受土壤中有机碳浓度影响较大，而由于灌溉污水中有机碳的输入，土壤中的有机碳浓度总体表现为先增加后趋于稳定，反硝化速率变化基本与有机碳含量变化一致，所以模拟的浅层硝态氮浓度先增加后减小并趋于稳定。由于吸附作用，有机碳一般分布在浅层土壤中，深层土壤有机碳分布较少，因此仅在10和40 cm深度处出现了这种趋势，而在70 cm深度处，没有此变化趋势。

图8 埋深为10、40和70 cm处铵态氮浓度模拟值与实测值随时间变化对比

图9 埋深为10、40和70 cm处硝态氮浓度模拟值与实测值随时间变化对比

总体来说，本文模型模拟的硝态氮和铵态氮浓度与实测值吻合程度较高，能反应土壤氮素的动态变化过程。

(4)平衡分析。图10(a)为模拟过程中硝态氮质量平衡柱状图，由图中可以看出，硝化与反硝化作用对硝态氮含量影响最大。模拟结束时，土柱中硝态氮含量比初始时有所增加，整个模拟过程中，硝态氮的绝对误差为292.74 mg，相对误差为4.58%。图10(b)为铵态氮模拟平衡柱状图，由图中可以看出灌水与有机氮的矿化是铵态氮的主要来源，硝化作用是铵态氮的主要减少途径，排水中几乎不含铵态氮。模拟结束时，铵态氮比初始时含量略有增加，整个模拟过程中，铵态氮的绝对误差为7.27 mg，相对误差为0.14%。

(图中硝态氮(铵态氮)质量增加为+，含量减少为-)图10 氮素平衡分析

4 不同外源性有机碳分布对氮素转化过程的影响数值实验

本文采用改进后的NPTTM模型对再生水灌溉土柱试验中的氮素迁移转化过程进行了模拟，验证了模型的适用性。为了进一步探究有机碳分布对氮素过程的具体影响，本文在保证所有土壤参数、溶质运移参数与各氮素转化参数不变的前提下，计算有机碳在土壤深度上不同分配情况下土壤中硝态氮与铵态氮的变化情况。

由于输入的有机碳是按照一定的分配方式加入土壤剖面，因此，不同的分配方式会对计算结果产生影响。本文设置了4组有机物沿土壤深度分配形式，类型A为模拟采用的有机碳分配，主要是考虑一般情况下，灌水中有机物集中在土壤表层，所以人为加大表层有机碳的权重；类型B为按深度均匀分配；类型C为按深度指数衰减分配；类型D在表层50 cm以内有机碳分配量按深度线性减小，50 cm以下没有机碳的补充。具体分配方式见图11。

图11 4种有机碳沿深度分配图

图12 采用不同输入有机碳分配方式时在10、40和70 cm深处铵态氮浓度随时间变化图

图13 采用不同输入有机碳分配方式时在10、40和70 cm深处硝态氮浓度随时间变化图

图12和图13列出了不同有机碳分布影响下不同深度处铵态氮、硝态氮浓度模拟值与实测值对比。类型A、类型B、类型C、类型D表层40 cm内有机碳含量占整个剖面有机碳含量的百分比分别为85.95%、36.36%、72.66%、96%。4种不同有机碳分布所对应的土壤10、40、70 cm深度最后一次灌水前(145 d)铵态氮浓度绝对最大变幅分别为0.131、0.01、0.048 mg/L，硝态氮浓度绝对最大变幅分别为20.6、6.4、2.13 mg/L。结果表明：不同的有机碳分布对土壤铵态氮的影响较小，但对硝态氮的影响较大。除类型B外，其余不同类型有机碳分布中，有机碳含量均随深度衰减，且土壤中铵态氮、硝态氮多集中在地表，所以表层铵态氮浓度、硝态氮浓度变化受有机碳分布影响都较下层更为显著。表层40 cm内有机碳含量占整个剖面有机碳含量的百分比从36.36%(类型B)变化到96%(类型D)时，土壤10 cm深度最后一次灌水前铵态氮浓度从0.437 mg/L变为0.568 mg/L，硝态氮浓度从49.1 mg/L变为28.5 mg/L，土壤70 cm深度最后一次灌水前铵态氮浓度从0.192 mg/L变为0.144 mg/L，硝态氮浓度7.01 mg/L变为9.14 mg/L。从趋势上来说，土壤中较高的有机碳集中会增加相应位置有机氮的矿化量从而造成铵态氮浓度的增加。对于硝态氮而言，有机碳随深度平均分配还是衰减型分配对结果影响很大。类型B为随深度平均分配，表层有机碳浓度较小，反硝化作用较弱，因而导致硝态氮浓度偏高，而在深层则由于有机碳浓度偏高，导致反硝化作用较强，硝态氮浓度偏低。有类似的研究也发现在人工湿地中添加简单的有机物(如葡萄糖)或者以轧棉工业的废弃棉花作为有机碳源直接覆盖在湿地表层能够大大提高反硝化速率[18,19]。

从以上分析可以看出，有机碳的分布对土壤氮素转化过程影响较大，因而得到合理的有机碳输入模式，对污水灌溉条件下土壤氮素动态变化过程的分析十分重要。

5 结 语

本文在土壤水分和氮素运移转化模型NPTTM基础上，考虑不同外源性有机碳分配方式，并耦合加入有机碳对反硝化的影响函数，对模型进行改进。该改进模型能很好地对含外源性有机碳的再生水灌溉条件下土壤水氮动态过程进行模拟，反应土壤氮素的运移规律。用验证后的模型计算4种不同有机碳随深度分配情况，结果表明：不同的有机碳分布对土壤铵态氮的影响较小，但对硝态氮的影响较大。表层40 cm内有机碳含量占整个剖面有机碳含量的百分比从36.36%(类型B)变化到96%(类型D)时，土壤10 cm深度最后一次灌水前铵态氮浓度从0.437mg/L变为0.568mg/L，硝态氮浓度从49.1 mg/L变为28.5 mg/L，土壤70 cm深度最后一次灌水前铵态氮浓度从0.192 mg/L变为0.144 mg/L，硝态氮浓度7.01 mg/L变为9.14 mg/L，相互之间差异很大。因此，采用合理有机碳输入模式，对污水灌溉条件下土壤氮素动态变化过程的模拟分析十分重要。

□

[1] Bixio D, Thoeye C, De Koning J, et al. Wastewater reuse in europe[J]. Desalination, 2006,187(1-3):89-101.

[2] Mackie A, Woszczynski M, Farmer H, et al. Water reclamation and reuse[J]. Water Environment Research, 2009,81(10):1 406-1 418.

[3] Zhu Yan, Shi Liangsheng, Lin Lin, et al. A fully coupled numerical modeling for regional unsaturated-saturated water flow [J]. Journal of Hydrology, 2012,475:188-203.

[4] Nielsen S N, Anastácio P M, Frias A F, et al. CRISP-crayfish rice integrated system of production. 5. Simulation of nitrogen dynamics [J]. Ecological Modelling, 1999,123(1):41-52.

[5] Arheimer B, Wittgren H B. Modelling nitrogen removal in potential wetlands at the catchment scale [J]. Ecological Engineering, 2002,19(1):63-80.

[6] Jamieson P D, Semenov M A. Modelling nitrogen uptake and redistribution in wheat [J]. Field Crops Research, 2000,68(1):21-29.

[7] Bouraoui F, Grizzetti B. Modelling mitigation options to reduce diffuse nitrogen water pollution from agriculture [J]. Science of The Total Environment, 2014,468-469:1 267-1 277.

[8] Franko U, Oelschlägel B, Schenk S. Simulation of temperature, water-and nitrogen dynamics using the model CANDY [J]. Ecological Modelling, 1995,81(1-3):213-222.

[9] 陆垂裕, 杨金忠, Jayawardane N, 等. 污水灌溉系统中氮素转化运移的数值模拟分析 [J].水利学报, 2004,(5):83-88,93.

[10] 陈 研, 胡克林, 冯 凌,等. 基于土壤-作物系统模拟模型的冬小麦田间水氮优化管理 [J]. 农业工程学报, 2007,(6):55-60.

[11] Zhu Yan, Ye Ming, Roeder Eberhard, et al. Estimating ammonium and nitrate load from septic systems to surface water bodies within ArcGIS environments [J]. Journal of Hydrology, 2016,532:177-192,

[12] 崔 敏, 冉 炜, 沈其荣. 水溶性有机质对土壤硝化作用过程的影响 [J]. 生态与农村环境学报, 2006,(3):45-50.

[13] Ge Shijian, Peng Yongzhen, Wang Shuying, et al. Nitrite accumulation under constant temperature in anoxic denitrification process: The effects of carbon sources and COD/NO3-N [J]. Bioresource Technology, 2012,114:137-143.

[14] 张 乐, 何红波, 章建新,等. 不同用量葡萄糖对土壤氮素转化的影响 [J]. 土壤通报, 2008, 39(4): 775-778.

[15] van Genuchten Th M. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils [J]. Soil Science Society of America Journal, 1980,44:892-898.

[16] 王丽影, 杨金忠, 伍靖伟, 等. 再生水灌溉条件下氮磷运移转化实验与数值模拟 [J]. 地球科学(中国地质大学学报), 2008,(2):266-272.

[17] 程先军, 许 迪. 碳含量对再生水灌溉土壤氮素迁移转化规律的影响[J]. 农业工程学报, 2012,28(14):85-90.

[18] Garcia-Montiel D C, Melillo J M, Steudler P A, et al. Carbon limitations to nitrous oxide emissions in a humid tropical forest of the Brazilian Amazon [J]. Biology and Fertility of Soils, 2003,38(5):267-272.

[19] Ullah S, Faulkner S P. Use of cotton gin trash to enhance denitrification in restored forested wetlands [J]. Forest Ecology and Management, 2006,237(1-3):557-563.