基于HYDRUS-1D模型模拟关中盆地氮在包气带中的迁移转化规律

2015-10-13王小丹王文科曹玉清夏雨波

华北地质 2015年4期

王小丹，凤　蔚，王文科，曹玉清，夏雨波

（1.中国地质调查局天津地质调查中心，天津300170；2.中国地质环境监测院，北京100081；3.长安大学环境科学与工程学院，西安710064；4.吉林大学环境与资源学院，长春130026）

王小丹1，凤蔚2，王文科3，曹玉清4，夏雨波1

以陕西关中盆地为研究区，利用HYDRUS-1D软件构建土壤水分运移模型、作物根系吸水模型和溶质运移模型，模拟了“三氮”在包气带中的迁移转化过程。结果表明：1）作物根系吸水吸氮规律一致，且主要吸收氨氮形式的氮素，吸收率为35%；2）亚硝态氮和硝态氮更容易被淋失到地下水中，且主要以硝态氮为主；3）不同包气带岩性对“三氮”向下迁移的速率和迁移量影响很大。

三氮；包气带；迁移转化；HYDRUS-1D；关中盆地

进入21世纪，我国氮肥的施用量超过作物地上部吸收量，氮肥对环境的负面影响开始凸显［1］，最突出的表现之一是“三氮”经由包气带淋失到浅层地下水从而对地下水造成污染，因此，氮素在包气带中的迁移转化规律及氮污染等问题成为国内学者近年来关注的热点。关中盆地作为陕西省重要的农业基地，长期以来受农业施肥和灌溉的影响，NO3-的淋失作用比较强烈，淋失量较大，使得关中盆地浅层地下水普遍受到NO3-污染的威胁［2-4］。自2002年开始，王文科率领的研究小组一直对关中盆地的“三氮”污染进行研究，分析“三氮”在包气带中的迁移转化规律、空间分布特征、地下水脆弱性评价、地下水氮污染、地下水氮污染健康风险评价以及与“三氮”有关的微生物分布及其反硝化能力等，取得了一系列研究成果［5-13］。

本文尝试在这些研究成果的基础上，在关中盆地首次利用美国盐土实验室（USSalinity laboratory）开发的HYDRUS-1D软件，将土壤水分运移模型、作物根系吸水模型和溶质运移模型相耦合，模拟研究区典型包气带结构中“三氮”的迁移转化过程，研究不同包气带岩性、结构对“三氮”迁移转化的影响。

1　研究区概况

关中盆地位于陕西省中部，地处黄土高原南缘，是一个三面环山，西高东低，西窄东宽的盆地，渭河自西向东穿过盆地中部。从南北山前到盆地中心，呈阶梯状依次分布有山前洪积扇、黄土台塬、冲积平原等地貌类型（图1）。第四系松散层孔隙水、孔隙裂隙水是关中盆地潜水的主要类型，其补给来源主要是降水和地表水的入渗等［2］。

2　模型建立

2.1数学模型

2.1.1土壤水分运移模型

土壤水分运移模型是用来描述水分在土壤中的运移过程的。本文采用经典的Richards方程来描述土壤水分运移过程，忽略土壤水平和侧向水流运动，仅考虑一维垂向运移。

式中：C（h）为土壤容水度（1/cm）；K（h）为非饱和渗透系数（cm/d）；h为压力水头或基质势（cm）；S（z，t）为单位时间单位体积土壤中根系吸水率（d-1）；t为时间变量（d）；Z为空间变量（cm），地表为原点，向下为正；L为潜水埋深（cm）；h0（z）为初始压力水头分布（cm）；ε（t）为入渗强度（cm/d）。

图1　关中盆地地貌图Fig.1 Geomorpho logy o f Guanzhong Basin

2.1.2作物根系吸水模型

根系吸水率表示由于根系吸水而在单位时间内从单位体积土壤中流失的水分体积，本文采用以水势差为基础的Feddes模型，即：

式中：α（h，z）表示水分胁迫反应方程；β（z）为根系吸水分布函数（cm-1）描述了根系吸水的空间变异；Tp为作物潜在蒸腾率，（cm·d-1）。

2.1.3溶质运移模型

在模拟中，考虑氮肥施用后在土壤中的运移和转化过程，采用传统的对流-弥散方程（Convection-DispersionEquation，CDE）来描述溶质运移过程。

对于NH4+，考虑吸附、硝化、和根系吸收作用，CDE方程为：

2.2初始条件及边界条件

对于水分运移的上边界，选择通量已知的大气边界条件（Atmospheric Boundary Condition w ithSurface Layer），在作物生育期内逐日加入通过上边界的变量值，包括降水量、灌溉量、棵间潜在蒸发率和作物潜在蒸腾率；其中，叶面拦截雨量忽略不计。由于表土层导水率较大，即使有较高强度降水或灌溉发生，也会很快入渗，因此，地面径流忽略不计。下边界选在潜水面处，使用定水头边界（Constant Pressure Head）。初始条件通过含水率给定，根据实验室测定的不同岩性中含水率的经验值给定，下边界处为饱和含水率。

对于溶质运移的上边界，选择溶质通量边界（Concentration Flux Boundary Condition），施肥量通过上边界水流通量与其浓度的乘积体现。下边界选择浓度零梯度边界条件（ZeroGradient），计算土壤水渗漏量和硝态氮淋失量。由于氮的转化过程十分复

2.3主要参数的确定

本次模型重点模拟“三氮”在包气带的迁移转化过程及去向，因此，模型参数选用研究小组试验测定的参数及研究区统计条件下的氮肥施用量和灌溉量等。

2.3.1土壤水分特征参数及溶质运移参数

土壤水分特征参数和溶质运移参数见表1、表2，主要通过原位及实验室试验测定获得。

表1　土壤水分特征参数表Tab le 1 Parameters o f the soilmoisture

2.3.2灌溉量和施肥量

研究区域种植的农作物主要是冬小麦和夏玉米，因此本次模拟选定这两种作物。根据冬小麦-夏玉米的需水量和关中盆地农作物长期以来的灌溉制度，在每个冬小麦-夏玉米的轮作期内共灌溉五次，分别为冬小麦的播前灌、返青灌、拔节灌和夏玉米的播前灌和拔节灌，总灌溉量为280m3/亩，其中冬小麦150m3/亩和夏玉米130m3/亩。另外，取各行政区多年施肥量的平均值分别作为该区内的施肥量，折算成施氮量加入到模型中，见表3。

表2　溶质运移模型参数Tab le 2　Parameters o f so lute transportmode l

表3　　不同地区农田施氮量（kg/hm2·a）Tab le 3　Fa rm land nitrogen amoun t in thedifferenta reas

2.3.3作物根系模型参数

研究区域种植的农作物主要是冬小麦和夏玉米，因此本次模拟选定这两种作物，模型参数利用软件中提供的经验参数，并忽略溶质对水分的胁迫作用。

由于玉米根系主要集中在土壤0～40 cm的深度，占整个根量的95%左右，所以玉米吸收水分和养分主要依靠这部分根［14］；同时，根据赵秉强等［15］的研究，玉米最大扎根深度可达160 cm，但99%以上的根量集中在120 cm土层以上。小麦根系深度一般在80～210 cm之间［16］，并且主要分布在70 cm土层内［17］。因此，综合考虑小麦和玉米的根系深度和特征，在模拟时将最大扎根深度取为100 cm。

3　模型运行结果及分析

在上述建模与确定参数的基础上，经过数值模拟得到根系吸水吸氮及三氮在不同岩性包气带中的迁移转化规律。

3.1根系吸水吸氮规律

以西安市氮肥施用量为例，模拟冬小麦从播种到收割期间根系吸收水分及“三氮”规律，模拟结果如下（图2、表4）：

（1）由图2可以看出，根系吸水吸氮规律一致。小麦根系对水分的吸收主要集中在播种后的前两个月和返青直至收割的后三个月，对“三氮”的吸收也同样集中在这段时间。关中盆地冬小麦一般在10月份播种，第二年6月收割，在冬季的三个月（12月到2月）里，由于小麦处于休眠期，基本不生长，因此，对水分和“三氮”的吸收量较小。

（2）比较小麦根系对“三氮”的吸收可以看出，小麦根系主要吸收形式的氮素，对亚硝态氮和硝态氮的吸收量很少，这也是和容易受到淋滤作用影响的原因。

图2　小麦根系对水分、氨氮、亚硝态氮和硝态氮的吸收Fig.2 Wheat root absorp tion o f the water，NH4+，NO2-and NO3-

表4　模拟计算区冬小麦生长季氮素平衡各项比较（以N计）（单位：kg/hm2）Tab le 4 Nitrogen ba lance o f simu lation area during the w inter wheatgrow ing season

（3）由模拟结果表4可以看出，当施氮量为大约254.8 kg/hm2时，根系对“三氮”的吸收量约为89.3 kg/ hm2，吸收率为35%，而大约有65%的氮素残留在土壤中，一部分被土壤颗粒吸附，一部分随水分向下运移而被淋滤到地下水中，从而对地下水污染构成威胁。

（4）模型运行结果显示，在冬小麦生长季氮的淋失量并不大，只有0.133 kg/hm2，这可能是因为模型计算时间较短，只计算了冬小麦的生长时间（约250天），大量的硝态氮还尚未运移到潜水面，所以模型结果显示的淋失量偏小。

3.2单一岩性中氮的迁移规律

3.2.1模拟结果

选择中砂作为单一介质进行模拟，以春玉米-冬小麦为种植作物，计算时间选择五年，模拟结果如图3-图5所示。

由图2-图4可见，在单层中砂中，氨氮、亚硝态氮和硝态氮都在随水流向下迁移，但是迁移规律不同。主要表现在：

图3　单层中砂中不同时间氨氮剖面分布图Fig.3 NH4+con ten t p ro file o f the m ono laye r sand in d iffe ren t tim e

图4　单层中砂中不同时间亚硝态氮剖面分布图Fig.4 NO2-content p ro file o f the mono layer sand in different time

图5　单层中砂中不同时间硝态氮剖面分布图Fig.5 NO3-con ten t p ro file o f the m ono laye r sand in d iffe ren t tim e

（1）氨氮进入土层后，随水分向下迁移量迅速减少，迁移速度较慢。在前五年，氨氮主要累积在大约1m以上的土层。第一年，氨氮迁移到大约50 cm的土层，在0～50 cm土层，氨氮浓度从0.043mg/cm3减少到0；第三年，氨氮在土壤表层出现累积现象，最大累积量出现在大约25 cm的土层，浓度达到0.047 mg/cm3，迁移至约85 cm的土层；第五年，最大累积量出现在大约30 cm的土层，达到0.05mg/cm3，迁移至约120 cm的土层。由此推算，氨氮随水分向下迁移速度约为15～20 cm/a。

（2）亚硝态氮随水分向下迁移速度较快，在土壤表层同样出现了累积现象。第一年，最大累积量大约出现在20 cm处土层，最大累积量约为0.01mg/cm3，亚硝态氮迁移到了400 cm以下的土层；第三年，最大累积量大约出现在30 cm处的土层，最大累积量约为0.018mg/cm3，亚硝态氮已经迁移到潜水面；第五年，最大累积量大约出现在50 cm处的土层，最大累积量约为0.023mg/cm3，进入潜水面的亚硝态氮通量大约为0.000 7mg/cm3。

（3）硝态氮随水分向下迁移速度快，在土壤表层也有累积现象。第一年，硝态氮就已经随水流迁移到了潜水面，最大累积量出现在大约25 cm的土层，最大累积量约为0.002 5mg/cm3；第三年，最大累积量出现在大约40 cm的土层，最大累积量约为0.018 mg/cm3，进入潜水面的硝态氮通量大约为0.002 2 mg/cm3；第五年，最大累积量出现在大约60 cm的土层，最大累积量约为0.023mg/cm3，进入潜水面的硝态氮通量大约为0.000 4mg/cm3。

3.2.2形成原因

分析其原因，主要有以下几点：

（1）由根系吸水吸氮模型结果可以知道，氨氮进入土壤后，有一部分由于植物固氮作用而被吸收，有一部分被土壤颗粒吸附而固定，还有一部分通过硝化作用转化成了亚硝态氮。因此，氨氮进入土层后浓度迅速减少，向下迁移速度很慢，并在土壤表层出现累积现象。

（2）由于氨氮的硝化作用，亚硝态氮也在土壤剖面中出现了累积；由于氨氮主要出现在大约1m以上的土层，亚硝态氮的最大累积量的峰值点也主要出现在土壤表层。同时，由于土壤颗粒对亚硝态氮几乎不吸附，而且中砂岩性颗粒较粗，透水性好，水分向下迁移速度快，从而导致亚硝态氮随水分向下迁移的速度也较快，在第二年已经迁移到了潜水面，并导致进入到潜水面的亚硝态氮通量随着时间的推移不断增大，对浅层地下水水质构成了威胁。

（3）经由硝化作用，亚硝态氮转化成了硝态氮，导致土壤剖面中也出现了硝态氮的累计现象，累积量峰值点也主要出现在1m以上的土层，且随时间变化向下不断推移。同时，由于土壤颗粒对硝态氮几乎不吸附，而且中砂岩性颗粒粗，透水性好，使得硝态氮迁移速度很快，淋失量也较大。

3.3不同岩性中氮的迁移转化规律

3.3.1模拟结果

图6　不同岩性结构中氨氮的迁移规律Fig.6 M ig ra tion o f the NH4+in d iffe rent litho logy

图7　不同岩性结构中亚硝态氮的迁移规律Fig.7 Migration o f the NO2-in different litho logy

以春玉米-冬小麦为种植作物，计算时间选择二十年，包气带厚度为5m，以关中盆地包气带两种比较典型的岩性中砂和黄土为例，进行不同岩性中氮素迁移转化规律研究。模型计算结果如图6-图8所示：

图8　不同岩性结构中硝态氮的迁移规律Fig.8 M ig ra tion o f the NO3-in d iffe ren t litho logy

由图6-图8可以看出，不同岩性条件下，氨氮、亚硝态氮和硝态氮随水流向下迁移的规律不同，主要表现在：

（1）20年后，氨氮在中砂中已经迁移到了300 cm的土层中，并且在大约40 cm处出现了累积峰值点，峰值浓度约为0.065mg/cm3；在同样的时间内，氨氮在黄土中只迁移到了大约25 cm处，并且浓度值随土层深度急剧减少，从土壤表层的0.06mg/cm3变为0。

（2）亚硝态氮在中砂中已经迁移到了潜水面，进入潜水面的通量大约为0.003mg/cm3，并且在土层大约80 cm处也出现了累积峰值点，峰值浓度约为0.022mg/cm3；而在黄土中亚硝态氮只迁移到了80 cm深的土层中，并且在0～30 cm处亚硝态氮浓度迅速从0.042mg/cm3减少到0.005mg/cm3，而在30～80 cm处其浓度从0.005mg/cm3减少到0。

（3）硝态氮在中砂中也迁移到了潜水面，进入潜水面的通量大约为0.011mg/cm3，并且在大约90 cm处也出现了累积现象，最大浓度值为0.012 5 mg/ cm3；而在黄土中亚硝态氮只迁移到了150 cm深的土层中，并且在土壤表层的0～30 cm处其浓度迅速从0.02mg/cm3减少到0.005mg/cm3，而在30～120 cm处其浓度从0.005mg/cm3减少到0。

3.3.2形成原因

“三氮”在中砂中随水分向下迁移的能力比在黄土中强，迁移速度快。这主要是由于粗颗粒岩性土层孔隙大，透水性好，土壤颗粒对氨氮吸附性弱，因此，氨氮易随水流向下迁移；而细颗粒土层岩性结构致密，透水性差，土壤颗粒比表面积大，对氨氮吸附性强，大量的氨氮由于吸附作用而被固定到土壤颗粒中，由氨氮转化来的亚硝态氮和硝态氮的量也相应变少，因此，氮素随水流向下迁移速度慢，大量的氮主要积累在了土壤表层。

4　讨论

（2）根据模型模拟结果，当施氮量为254.8 kg/ hm2时，根系对“三氮”的吸收量为89.3 kg/hm2，吸收率仅为35%，这与朱兆良等人关于我国氮肥利用的研究成果相符，即“当季作物（对氮肥）的吸收利用率约为35%”［18，19］，而剩余65%的氮素经由吸附、硝化—反硝化、氨挥发、淋滤等作用而损失。因此，改变传统的施肥灌溉模式，适当控制施肥量，不仅可以提高氮肥的利用率，而且对防治地下水氮污染也有重要的作用。

（3）氨氮进入土层后，由于植物吸收、土壤颗粒吸附和硝化作用而迅速减少，随水分向下迁移速度很慢，在中砂中其速度大约为15～20 cm/a，第五年时，氨氮迁移到120 cm处，最大累积量为0.05mg/cm3，出现在30 cm处。亚硝态氮和硝态氮由于根系吸收、硝化、反硝化及吸附作用等也迅速减少，但是随水分向下迁移速度很快，第一年都已经迁移值五米处的潜水面，到第五年时，进入潜水面的亚硝态氮和硝态氮通量分别为0.000 7mg/cm3和0.004mg/cm3，进一步证明，氮素主要以硝态氮形式进入地下水，亚硝态氮次之。

（4）不同包气带岩性对“三氮”向下迁移影响较大，表现为粗颗粒岩性中“三氮”随水流向下迁移速度较快，而细颗粒岩性中迁移速度很慢。在模型模拟的二十年时间里，当氨氮在中砂中已经迁移到大约300 cm处时，在黄土中只迁移到大约25 cm处；亚硝态氮和硝态氮在中砂中都已经迁移至五米处的潜水面，进入潜水面的通量分别为0.003mg/L和0.011 mg/L，而在黄土中亚硝态氮和硝态氮分别迁移至80 cm和150 cm深处。可见，“三氮”在包气带中的迁移速率主要受包气带岩性影响，粗颗粒岩性中迁移速度快，细颗粒岩性中迁移速度慢。

5　结论

经过对比，本次模型模拟结果与前人已有研究成果相符，因此，利用HYDRUS-1D软件可以较好地反映出“三氮”在包气带中的迁移转化规律，研究方法可行，根据模拟结果，主要得出以下结论：

（1）作物根系吸水吸氮规律一致，且主要吸收NH+4形式的氮素。在我国目前氮肥施用水平上，作物对氮素的吸收率较低，仅为35%左右，剩余NH+4形式的氮素由于吸附、硝化—反硝化、挥发、淋滤等作用而迅速减少，因此，随水流向下迁移速率较慢。

（2）作物对亚硝态氮和硝态氮的吸收量很少，同时，包气带对其吸附能力相对较弱，因而，向下迁移的量相对较大，而NO2-不稳定，大部分经过硝化过程转化为相对稳定的NO3-，因此，硝态氮受淋滤作用影响最明显，淋失量最大。

（3）不同包气带岩性对“三氮”向下迁移影响较大，表现为粗颗粒岩性中“三氮”随水流向下迁移速度较快，而细颗粒岩性中迁移速度很慢。

总之，模拟结果进一步证明，氮素主要以亚硝态氮和硝态氮的形式进入地下水中，并以硝态氮为主，其迁移速度和迁移量受包气带岩性影响很大。

［1］张卫峰，马林，黄高强，等.中国氮肥发展、贡献和挑战［J］.中国农业科学，2013，46（15）:3161-3171.

［2］姜桂华.关中盆地地下水脆弱性研究［D］.西安:长安大学，2002.

［3］姜桂华，王文科，乔小英，等.关中盆地地下水特殊脆弱性及其评价［J］.吉林大学学报（自然科学版），2009，39（6）: 1106-1110.

［4］王小丹.关中盆地氮对浅层地下水污染的数值模拟与预警研究［D］.长安大学，2008.

［5］姜桂华，王文科，杨晓婷，等.关中盆地地下水硝酸盐污染分析及防治对策［J］.水资源保护，2002，68（2）:6-8.

［6］李翔.关中盆地马兰黄土中三氮转化及相关微生物关系实验研究［D］.西安:长安大学，2003.

［7］李翔，王文科，杨胜科，等.马兰黄土中与三氮有关的微生物的分布及其反硝化能力［J］.环境污染与防治，2006，28 （6）:411-414.

［8］姜桂华，王文科.关中盆地包气带氮迁移转化数值模拟及预测［J］.西北大学学报（自然科学版），2007，37（5）: 825-829.

［9］邓林，王文科，杨晓婷，等.关中盆地地下水硝酸盐含量的空间变异特征［J］.干旱区资源与环境，2008，22（10）: 152-155.

［10］段磊.关中盆地地下环境氮污染机理与地下水质安全评价［D］.西安:长安大学，2010.

［11］段磊，王文科，孙亚乔等.关中盆地浅层地下水氮污染的健康风险评价［J］.水文地质工程地质，2011，38（3）:92-97.

［12］段磊，王文科，杨晓婷，等.关中盆地浅层地下水氮污染的时空变化规律及其防治措施［J］.干旱区资源与环境，2011，25（8）:133-137.

［13］李慧，王文科，段磊，等.关中盆地包气带中“三氮”的分布特征及影响因素分析［J］.安徽农业科学，2013，41（10）: 4567-4570，4595.

［14］王树安.作物栽培学各论［M］.北京：中国农业出版社，1995.136，169-170.

［15］赵秉强，张福锁，李增嘉，等.间套作条件下作物根系数量与活性的空间分布及变化规律研究.Ⅱ.间作早春玉米根系数量与活性的空间分布及变化规律［J］.作物学报，2001，27（6）：974-979.

［16］马元喜，王晨阳，周继泽.小麦根系主要生态效应的研究［J］.河南农业大学学报，1994，28（1）:12-18.

［17］马福生，刘洪禄，吴文勇，等.再生水灌溉对冬小麦根冠发育及产量的影响［J］.农业工程学报，2008，24（2）：57-63.

［18］Zhu Z L，Chen D L.Nitrogen fertilizer use in China -Contributions to food production，impacts on the environment andbestmanagement strategies［J］.Nutrient Cycling inAgroecosystems，2002，63:117-127.

［19］朱兆良，金继运.保障我国粮食安全的肥料问题［J］.植物营养与肥料学报，2013，19（2）:259-273.

［20］徐军祥，邢立亭.济南泉域岩溶水数值预报与供水保泉对策［J］.地质调查与研究，2008，30（3）；209-213.

P333.9

1672－4135（2015）04－0291-08

2015-09-15

陕西省环保局环境科研项目“陕西省地下水污染现状调查（Z0709009）”

王小丹（1982-），女，硕士，工程师，主要从事环境地质水文地质研究，E-mail：71573377@qq.com。