一种结合空间显式水文模型模拟土壤氮淋失的方法

1.本发明属于水文水化学领域，具体涉及一种结合空间显式水文模型模拟土壤氮淋失的方法。


背景技术：

2.地表淡水是人类淡水资源的重要组成部分，是人类生活、农业与工业用水的重要来源。氮素(主要亚硝酸盐氮，硝酸盐氮和氨氮)是地表水的主要污染物之一。农业生产过程中大量使用化肥是水体中氮素的重要来源之一，残存在土壤中的氮素肥力从农业集水区的土壤释放(也称淋失，下文不作区分)到地下水中，并进一步排放到地表水体，使得地表水质富营养化，直接威胁人类饮用水安全。有效地模拟该释放过程，对于理解地下地表水氮污染过程以及农业管理具有十分重要的意义。
3.在研究地下水流场和氮排放方式的过程中，研究人员常常针对地下水进行数学建模。基于物理过程的、空间显式的数值模型(如hydrogeosphere，mt3dms，parflow等)能模拟三维空间内连续水流和溶质运移过程，提供流场和溶质分布。但利用此类模型对氮素的运移、排放进行模拟的有效性需借助准确描述氮源从土壤运移到地下水的过程来实现。
4.氮随着降雨进入浅层土壤中并经过离子化、反硝化作用、植被吸收的过程后释放到地下水中。许多科学工作者已经提出了一些算法来模拟这一释放过程，根据特点分为两类：一阶衰减法和扩散法。一阶衰减法将淋失过程描述为一阶衰减过程，需要构造衰减常数和限制变量，常见的限制变量包含土壤的水饱和度或降雨量。扩散法将淋失过程描述为扩散过程，需要构造扩散系数和计算土壤和地下水之间浓度差。这些方法能在一定程度上描述淋失过程，并在实际案例中得到运用，但与上述数值模型相结合时会面临一些困难，如深层单元格无法确定降雨量、土壤和地下水之间的浓度差计算困难等。这导致模拟的结果具有较高的系统性误差。


技术实现要素：

5.发明目的：本发明针对现有技术中存在的不足，公开一种结合空间显式水文模型模拟土壤氮淋失的方法，该方法能模拟土壤无机氮的淋失过程，可以更好地与空间显式数值模型结合，从而模拟氮的输入、运移和排放过程。
6.技术方案：本发明所公开的结合空间显式水文模型模拟土壤氮淋失的方法，其特征在于，包括以下步骤：
7.s1、基于物理过程建立空间显式水文数值模型，结合水文数值模型在研究区域里计算浅层多孔介质中每一个剖分单元格(共n个)的无机氮载量sin；
8.s2、对于第j个单元格(j为单元格编号，1≤j≤n)，获取多孔介质参数及模拟的水文变量；
9.s3、对于第j个单元格和特定时间步长δt，求土壤无机氮sin的淋失速率lea，淋失的无机氮进入地下水形成溶解无机氮din；
10.s4、对于第j个单元格和特定时间步长δt，淋失速率lea作为源汇项，参与水文数值模型的求解，得到研究区域din浓度的空间分布，及边界处的din排放浓度；
11.s5、将淋失速率lea在整个研究区域和模拟时间积分，求得全区域、全时间跨度的土壤氮淋失速率。
12.其中，s1中，所述水文数值模型包含溶质运移模块，所述的无机氮载量sin在空间和时间上是可变的。
13.进一步的，s2中，所述水文变量为孔隙水的达西流速q。
14.进一步的，s3具体为：
15.在第j个单元格空间内和特定时间步长δt内，基于达西流速q、孔隙率θ、垂向土壤厚度d、以及淋失系数a，更新该步长时间内土壤无机氮淋失速率lea,
[0016][0017]
从上式可知，a为唯一参数，且决定淋失速率lea随时间变化的关键变量为孔隙水的达西流速q。
[0018]
进一步的，s3-5中，所述淋失速率lea为单元格上时空可变的变量。
[0019]
有益效果：与现有技术相比，本发明优点为：首先，本发明提出的淋失算法采用了局部变量达西流速，易在水文数值模型中获取，因而对空间显式数值模型具有较好的适应性；其次，本发明模拟氮的输入、运移和排放过程的有效性高；最后，本发明提出氮淋失为达西流速的函数，是一种改进的一阶衰减法算法，相比常规衰减算法具备较高的精度和较小的误差。
附图说明
[0020]
图1为本发明实施例中的研究区域，左图为空间显式水文模型网格，右图为该区域含水层土壤类型；
[0021]
图2为氮在土壤中的转化路径概念图；
[0022]
图3为实施例中分别使用本发明与对照算法得到的2001年至2010年排放浓度(集水区出口处)与测量值的比较；
[0023]
图4为实施例中分别使用本发明与对照算法模拟得到的氮载量与通量与真实值的比较。
具体实施方式
[0024]
下面结合附图和实施例对本发明做进一步的解释说明。
[0025]
一种结合空间显式水文模型模拟土壤氮淋失的方法，包括以下步骤：
[0026]
s1、基于物理过程建立空间显式水文数值模型(有限元法)，结合水文数值模型在研究区域里计算每一个剖分浅层多孔介质(土壤)单元格的无机氮载量sin，所述水文数值模型包含溶质运移模块，所述无机氮载量sin在空间和时间上是可变的。
[0027]
s2、对于第j个单元格，获取多孔介质参数及模拟的水文变量，所述水文变量为孔隙水的达西流速q。
[0028]
s3、对于第j个单元格和特定时间步长δt，采用有效的算法求土壤无机氮sin的淋失速率lea，所述淋失速率lea为单元格上时空可变的变量，淋失的无机氮进入地下水形成溶解无机氮din。具体为：在第j个单元格空间内和特定时间步长δt内，基于达西流速q、孔隙率θ、垂向土壤厚度d、以及淋失系数a，更新该步长时间内土壤无机氮淋失速率lea,
[0029][0030]
从上式可知，a为唯一参数，且决定淋失速率lea随时间变化的关键变量为孔隙水的达西流速q。
[0031]
s4、对于第j个单元格和特定时间步长δt，淋失速率lea作为源汇项，参与水文数值模型(溶质运移模块)的求解，得到研究区域din浓度的空间分布，及边界处的din排放浓度；
[0032]
s5、将淋失速率lea在整个研究区域和模拟时间积分，求得全区域、全时间跨度的土壤氮淋失速率。
[0033]
在一个小型农业集水区中验证本发明是能结合空间显式数值模型有效地模拟土壤氮淋失的过程。
[0034]
如图1所示，该区域面积约1.4平方公里，含水层厚度约2～5米，地表种植农作物为小麦，该区域两边为坡地，地下水从中间较低处汇集进入河流并排放至下游。首先针对该区域建立三维空间显式模型，离散为13860个单元格(棱柱)，利用hydrogeosphere(或其它软件如modflow，mt3dms，parflow等)求解地下水流方程和溶解无机氮din运移方程。
[0035]
在浅层土壤(0.5m厚度)单元格中构建外部氮输入
→
有机氮son
→
无机氮sin的简化的转化路径(如图2所示)，根据资料设定该路径的一系列参数，用于计算sin的载量，即潜在的最大淋失速率。
[0036]
在每一个时间步长内，对于每一个获取浅层土壤单元格，获取孔隙率，土壤厚度(0.5m)，达西流速三个参数或变量值。达西流速q可由模拟得三方向流速求得：
[0037][0038]
在每一个时间步长内，对于每一个获取的浅层土壤单元格，根据s2中获取的参数或变量值，及sin载量，计算淋失速率：
[0039][0040]
上式中，a为常数，在模型校正中进行调整。为了与本发明进行对照比较，另用一种普遍采用的、以土壤湿度为限制变量的淋失算法计算lea[yang et al.,2021]：
[0041]
对照算法：lea＝fw·
λ
·
sin，其中fw＝(s-sr)/(1-sr)
[0042]
上式中fw为土壤湿度限制因子，λ为衰减系数(在模型校正中进行调整)，s为土壤水饱和度，sr为土壤残余饱和度。
[0043]
在每一个时间步长内，对于每一个获取的浅层土壤单元格，将淋失量作为源汇项添加到s1所述水文模型的din(溶解无机氮)运移方程中，并进一步利用hydrogeosphere求解，可得到din浓度的空间分布，以及该集水区出口处的din浓度。将该浓度与2001年至2011
年间实际测量值进行对比，并以此来校正模型参数，以达到最佳吻合度和最小误差。经过对参数的校正，通过本发明算法和对照算法所能达到的最佳模拟结果如图3所示。结果表明，本发明能较好的模拟氮排放浓度(nse值达到0.52，表1)。与对照算法相比，本发明与实际测量值更吻合(nse较大，表1)，误差更低(rmse较小，表1)。因此可见，本发明模拟结果具备较高的精度和较小的误差。
[0044]
表1发明与对照算法模拟得到的地表水氮浓度与测量值的吻合度与误差
[0045] 对照算法本发明理想值吻合度nse0.49[-]0.52[-]1[-]误差rmse1.84mg/l1.80mg/l0mg/l
[0046]
分别将本发明和对照算法所模拟得到的氮载量(sin等)和氮通量(lea等)，在整个研究区域和模拟时间(2001年至2011年)积分，求得区域年平均值，如图4和表2所示。本发明计算的各项值与真实值(相对可信)基本吻合，说明本发明能较好地模拟氮淋失及其它过程。与对照算法相比，本发明模拟结果的误差与其相似或较低。
[0047]
表2本发明与对照算法模拟得到的各项氮载量与通量值
[0048][0049]
通过以上结合水文数值模型的运用与分析可知，本发明使用结合空间显式水文模型模拟土壤氮淋失的方法能够有效的模拟土壤氮淋失，并准确地重现河道中无机氮浓度、氮载量及通量。与其它技术相比，该方法针对空间显式模型适应性高，计算的精度更高。