一种农田农业水文生态系统模拟方法

本发明涉及水文生态系统模拟，特别涉及一种农田农业水文生态系统模拟方法。

背景技术：

1、作物模型是一种基于作物生长发育规律的数学模型，能够模拟作物在不同环境条件下的生长和发育过程。随着计算机技术和数据采集技术的发展，作物模型在农业生产和研究领域中得到了广泛应用，可以帮助研究者更好地理解作物的生长规律，预测作物产量和品质，优化农业生产管理措施，提高农业生产效率。

2、作物模型不受时间和地域限制，能够准确模拟作物的生长发育和产量形成等生理过程。根据模拟的作物种类，可以分为通用模型和专用模型。通用模型如aquacrop、epic和cropsyst，可以描述水稻、小麦、玉米等多种作物的生理生态多过程。专用模型如水稻模型oryza和ceres-rice，能够对水稻品种属性和稻田条件进行更细致和准确的描述。

3、目前国际上广泛应用的水稻生长模型有oryza系列模型、ceres-rice模型、apsim-rice模型、dssat模型、warm模型等。但这些模型一般都只侧重于水稻模拟的某一方面，如涉及水稻生理过程的有ceres-rice和wofost、土壤氮素转化过程的有dndc和daisy、管理措施表现的如whcns和cropsyst等。其中，oryza、aquacrop和ceres-rice是应用广泛的水稻稻田模拟代表模型。oryza模型是一个详细描述了水稻生长和非生物因素的典型水稻模型，已经在多个国家、多种水稻生长环境下广泛应用。

4、目前的水稻生长模拟研究多集中于热带/亚热带地区，相关的参数集也都是热带/亚热带地区应用得出的。对于寒地水稻生长的系统参数化和模拟研究还较少。在不同的土壤和气候种植环境中，水稻的生长状况变化很大，特别是在我国拥有肥沃黑土地的东北寒地地区，其无论在气候条件、自然环境还是土壤性质上，均与热带/亚热带地区存在较大差异。而了解黑土地的氮素利用效率，对于保护黑土地的质量具有重大意义。

5、现有技术中缺少适用于寒地水稻生长-水文生态-产量形成过程模拟的、能够调控生育期灌溉控制阈值、定量模拟土壤氮素归趋方-氮动态转化/利用的水稻模型。oryza2000模型是由荷兰瓦赫宁根大学和国际水稻研究所irri联合开发的模拟水稻生长和产量形成过程的水稻专用模型，广泛应用于水稻生产和农业研究领域。虽然该模型能够模拟水稻的生育期的逐日生物量累积、叶面积指数、光合产物分配、籽粒产量和氮素分配等关键过程，以评估不同农业管理措施对水稻生产的影响。但该模型在土壤有机质丰富且气候寒冷的东北黑土区应用时，其模拟效果差，且对于稻田节水灌溉制度调控上存在灌溉触发机制及灌水量单一，不足以满足控制灌溉需求的相关问题。因此，本发明旨在完善稻区的农田水文-生态系统过程模拟机理，基于完备的作物模块生物量分配框架，提出寒地稻区农田水文-生态多过程模拟普适性模型和模拟方法。

技术实现思路

1、本发明的主要目的在于提供一种农田农业水文生态系统模拟方法，可以有效解决背景技术中的问题。

2、为实现上述目的，本发明采取的技术方案为，

3、一种农田农业水文生态系统模拟方法，所述方法通过构建水稻耦合模型对农田水文生态多过程进行模拟，所述水稻耦合模型包括：用于计算逐日潜在蒸散发etp、潜在蒸发量ep、潜在作物蒸腾量tp的et模块；用于根据生物量累积和分配方式模拟水稻生长和产量形成的作物生长模块；所述水稻耦合模型还包括：

4、用于对土壤剖面中的水分再分配、排水和入渗进行模拟的土壤水模块；

5、用于模拟土壤各层的温度，并将逐日逐层的土壤温度用于计算氮素各转化过程中的温度修正因子的土壤热模块；

6、用于描述稻田土壤中的微生物对氮素矿化/固持影响过程的土壤碳氮模块；

7、用于描述nh4-n和no3-n在土壤介质中的迁移转化，涉及对流、弥散和各种反应过程的土壤氮模块。

8、所述土壤水模块的计算过程包括：

9、确定应用于稻田土壤表面的积水层和土壤剖面中的每一层水均衡的计算方式，对给定时间步长的所有水均衡分量、田面积水深度和逐层土壤体积含水量的状态进行模拟，其中，低洼稻田的水均衡计算方法如下:

10、δw＝δws+δhpond＝(i+r+c-e-t-s-p-d)×δt  (62)

11、式中，δw表示稻田土壤剖面和积水层的总存储水变化，δws表示土壤水分存储量的变化(mm)，δhpond表示稻田表面的积水深度变化，i表示灌溉供水(mm d-1)，r表示降雨(mmd-1)，c表示从地下水向土壤中的毛细上升(mm d-1)，e表示土壤蒸发(mm d-1)，t表示作物蒸腾(mm d-1)，s是田块的侧向渗流损失(mm d-1)，p是垂直渗流损失(mm d-1)，d是地表排水/产生的径流(超过田坎高度溢出形成的排水/径流)(mm d-1)；

12、设定灌溉、降雨和毛管上升水是土壤水分补给的来源，只有当土壤水张力大于相应的田间持水量时，从地下水到土壤剖面的毛管上升才会发生，在darcy方程的基础上，毛管上升量作为渗透系数、土壤水张力和地下水埋深距离的函数进行估计，其中，毛管上升只能随着与地下水埋深距离的增加而减小，而且是从离地下水位更近的土壤分层开始计算；如果计算出的土壤剖面上层的毛细管上升量比下层大，则毛细管上升被重置为离地下水位较近的土层的数值；

13、采用penman-monteith方法计算eto后，进一步将其分为辐射项etrd和空气动力项etae，用于估算稻田潜在的蒸发速率和蒸腾速率ep和tp，在考虑土壤水分胁迫的条件下，确定了实际蒸发速率和蒸腾速率ea和ta；

14、根据darcy方程对梯田的侧向渗流进行估算，在两个相邻田块之间没有高差时，侧向渗流量为0；

15、利用newton-raphson迭代方法根据犁底层的水力传导特性、土壤水梯度和土壤厚度计算犁底层的垂直渗透量，超过稻田容量的水被排出土层，其最大排出速率等于犁底层的饱和导水率，将蒸散发和渗漏损失计算在内，超过田坎高度的积水为田面径流量；

16、设定根系直接从根区所在土层的土壤水中吸收水分，在一定深度的潜在蒸腾速率sp(z)通过根区土层的潜在蒸腾速率(tp)计算得出；对于高地/旱地水稻，根系密度分布对根系吸水有很大影响；因此，sp(z)通过在任意形状的根区引入非均匀分布的tp来确定，是指定根长密度的函数，如下所示：

17、

18、式中，zr是根系深度(cm)，br(z)是相对根系吸水分布，通常与根系长度的密度分布(-)有关，最终的sp(z)根据bouman等的水分胁迫计算结果等比例削减；

19、采用如下排水计算方式模拟土壤剖面中的水分再分配、排水和入渗，具体为：

20、

21、式中，δθ/δt是单位时间土壤体积含水量的减少量，δt(cm3cm-3d-1)，τ是排水特性(-)，θsat是饱和时的土壤体积含水量(cm3cm-3)，θfc是田间持水量(cm3cm-3)。

22、所述土壤热模块的计算过程包括：

23、通过傅里叶定律对热通量密度(jt，w m-2)与土壤温度(t，℃)之间的关系进行描述：

24、

25、式中，λ是土壤表观导热系数(j s-1m-1℃-1)，z是各土层的深度(m)；式(1)与连续性方程相结合得到随时间变化的一维(1-d)热扩散的微分方程：

26、

27、式中，ch是土壤的体积比热(j m-3℃-1)，t是时间(s)，t是土壤温度(℃)；ch和λ的值是根据土壤成分和土壤实际体积含水量θ(cm3cm-3)估算的，具体计算如下：

28、

29、λ＝(a+bθ-(a-d)exp(-(cθ)4)  (4)

30、a＝0.65-0.78ρ+0.60ρ2  (5)

31、b＝1.06ρ  (6)

32、c＝1.0+2.6clay-0.5  (7)

33、d＝0.3+0.1ρ2  (8)

34、式中，ρ是土壤干容重(g cm-3)，clay是土壤中粘粒含量的比例(g g-1)。顶部边界由空气温度(tah，℃)和空气与第一个土层之间的热导率(λ0，w m-1℃-1)驱动；底部边界条件设定为土壤底部的观测土壤温度；

35、如果导热系数和土壤热容量不随深度变化，则可将λ从偏导中取出并与ch相结合，得到

36、

37、式中，dh为热扩散率(m2s-1)；

38、设定土层间传递的热量可以储存在土层中，也可以从土层原本储存的热量中提取；采用中心有限差分格式在时间和空间维度求解式(2)；对于给定的土层i，式(2)可写成式(10)：

39、

40、δt＝tj+1-tj   (11)

41、δzu＝zi-1-zi  (12)

42、δzl＝zi-zi+1  (13)

43、式中，δt是时间增量，上标j和j+1分别表示t和t+δt时的值，是加权平均温度(℃)，δzi是第i个土层的厚度；可以通过tj和tj+1时刻的温度计算得到：

44、

45、式中，η是一个权重系数，范围在0到1.0之间，一般取η值为0.6；结合每个节点的式(10)和式(14)，将各方程组写成以下矩阵格式：

46、

47、对于土层i，j+1时刻的新温度为tni，j时刻的温度为ti，而di为边界温度(℃)；对于从2到n-1的节点：

48、

49、式中，

50、

51、

52、

53、

54、对顶部节点，即节点＝1：

55、

56、a1＝0  (22)

57、

58、

59、

60、式中，tnair为已知的j+1时刻的空气温度(℃)；

61、对底部节点，即节点＝n：

62、

63、式中，

64、

65、

66、cn＝0  (29)

67、

68、式中，tnb为已知的j+1时刻的底部温度(℃)，该数值采用土壤温度的实测值；

69、其中，土壤热模块计算的结果是时间步长为1小时的土壤分层温度，而传递给其他模块的每层土壤的日平均温度(td，℃)是该模拟日的24小时内逐小时土壤温度的平均值。

70、所述土壤碳氮模块的计算过程包括：

71、设定溶质nh4-n或no3-n随着土壤水的流动而被输送到稻田积水、土壤剖面各层等不同位置，其中，qnra和qnir分别为降雨和灌溉水中的输入氮量(mg m-2)；i为剖面中土层的顺序；qn0和qni分别为田面积水和第i层土壤中的氮量(mg m-2)；c0和ci分别为田面积水和第i层土壤中的氮浓度(mg l-1)；qnrunf为径流中的氮损失量(mg m-2)；qncrop,i、qnl,i和qngw,i分别是作物吸收的氮量、氮的侧向渗漏量和地下水补给的氮量(mg m-2)；gwl是地下水位(m)；

72、溶质的向下运动与田间水分的垂向渗漏有关，而来自地下水的毛管水上升则会导致溶质向上运移；基于土壤水均衡方程，得到田面积水和每个土壤层的溶质平衡方程；

73、对于田面积水：

74、

75、式中，sw是田面积水中的溶质总量(mg m-2)，ir和ra分别是灌溉和降雨(mm d-1)，cir和cra分别是灌溉水和降雨中的溶质浓度(mg l-1)，q0是田面积水层和第一个土层之间的水通量(mm d-1)，ru是田面积水超过田坎高度后溢出田坎形成的地表径流(mm d-1)，c0是在一个模拟时间步长结束时更新的田面积水溶质浓度(mg l-1)，上标j和j+1分别表示t和t+δt时刻的值，δt是计算时间步长(＝1天)；

76、对于土壤各层：

77、

78、式中，下标i是指从表层开始向下的土壤层数，si是第i层的储存溶质总量(mg m-2)，qi是第i层和(i+1)层之间的水通量(mm d-1)，ci是第i层的溶质浓度(mg l-1)，cr是地下水的毛管上升量(mm d-1)，zi是第i层土壤的厚度(mm)，sc是氮素转化过程的每日源汇项(mgl-1d-1)；

79、其中，nh4-n和no3-n的每日源汇项(sc)计算如下：

80、scnh4＝rhyurea+rad+rmin-rvot-rnit-rup  (33)

81、scno3＝rnit-rden-rup  (34)

82、式中，scnh4和scno3分别为nh4-n和no3-n浓度(mg l-1d-1)，rhyurea、rad、rmin、rvot、rnit、rden和rup分别为尿素水解、固体吸附、矿化、氨的挥发、硝化、反硝化和作物吸收氮的速率(mgl-1d-1)；

83、各源汇项的估算假设和计算方法如下：

84、a、尿素水解

85、假定地表肥料，即尿素，均匀地融入表土层，参考apsim模型，随后开始尿素的水解反应；

86、则每层土壤的实际水解率(rhyurea，mg n l-1d-1)估算如下：

87、rhyurea＝min[1,fureamin(fw_urea,ft_urea)]urea  (35)

88、式中，urea为每个时间步长中尿素的最大可利用量(mg n l-1d-1)；furea为潜在的尿素水解速率；fw_urea和ft_urea分别为尿素水解的土壤水分和土壤温度降低系数；furea、fw_urea和ft_ureaf的计算公式如下：

89、furea＝min(1.0,1.12+1.31oc+0.203ph-0.155oc×ph)  (36)

90、

91、ft_urea＝0.027t+0.19  (38)

92、式中，oc，ph和t分别为土壤各层的有机碳含量(％)、氢离子浓度(-)和土壤温度(℃)；

93、b、固态吸附

94、溶质每天在土壤颗粒表面的吸附和释放的速率rad的计算方法如下：

95、

96、式中，qd是土壤颗粒的溶质吸附能力(mg kg-1)，kd是freundl ich系数(cm3g-1)，cref是用来使n无量纲的溶质浓度的参考值(mg l-1)，c是该层土壤溶质浓度(mg l-1)，n(-)是freundlich指数；

97、c、氨挥发

98、氨挥发是氨向大气的损失，土壤中氨挥发的损失受土壤温度(t，℃)和深度(d，cm)的影响，氨挥发率和校正因子可通过下式计算得出：

99、

100、ft_vot＝0.25exp(0.0693t)  (41)

101、fd_vot＝exp(-0.05d)  (42)

102、式中，rvot为潜在氨挥发速率(mg n l-1d-1)；c(nh4)为各土层中氨态氮(nh4-n)浓度(mg n l-1)；ft_vot和fd_vot分别为土壤温度和土壤深度对nh3挥发的校正系数(0-1)；

103、d、硝化作用

104、硝化过程是一个好氧过程，与各土层的nh4-n浓度(mg n l-1)、土壤深度(d，m)、土壤温度(t,℃)和土壤体积含水量(θ,cm3cm-3)有关；硝化速率(rnit,mg n l-1d-1)用下式进行估算，其中，日最大硝化量限制在总可利用nh4-n的10％；

105、rnit＝max[0,min(0.1,ft_nitfw_nitfd_nit)]c(nh4)  (43)

106、由于硝化需要氧气，因此土层中的氧气交换会影响硝化率，在土壤层中，通气性随着土壤深度的增加而降低，因此氧气供应也随之减少；假定硝化率每增加30cm减少55％，则与土壤深度相关的校正因子fd_nit计算如下：

107、

108、与硝化作用相关的土壤温度ft_nit和水分校正因子fw_nit计算如下:

109、

110、

111、式中，tavg为土壤层的平均温度；硝化产物为硝态氮(no3-n)和一氧化二氮(n2o)，比例分别为98％和2％；

112、e、反硝化作用

113、反硝化是一个厌氧过程，假设只有当土壤水饱和度高于0.80时才会发生反硝化；反硝化速率(rden,mg n l-1d-1)采用一阶方程计算，与土壤温度(t,℃)和有机碳含量(oc,％)的函数关系如下：

114、rden＝min[0.05,1-exp(-1.4ft_denoc)]c(no3)  (47)

115、

116、式中，ft_den为反硝化土壤温度校正因子(0-1)，c(no3)为各土层硝态氮(no3-n，mgn l-1)浓度；

117、在饱和和非饱和条件下，利用下式计算反硝化过程中的n2o排放量，反硝化产生的n2o计算如下:

118、

119、fw_den＝exp(-23.77+23.77wfps)  (50)

120、式中，fw_den是土壤湿度因子(-)，wfps是土壤孔隙中充满水的比例(0-1)，αden是总反硝化率rden中以n2o形式出现的最大比例(-)，默认值为0.5；其中，在硝化过程中损失的n2o部分范围为0.001到0.05；

121、f、作物吸氮速率

122、作物吸氮速率rup的计算公式如下：

123、rup＝kr·ta·c  (51)

124、式中，kr是根系吸收溶质校正因子(-)，ta是根区层的实际作物蒸腾量(mm d-1)。

125、所述土壤氮模块的计算过程包括：

126、设定土壤有机物被分为三个氮库，具体为：

127、代表更易腐烂的土壤微生物生物量和微生物产物的微生物生物量库bom，其中，土壤中的bom通常占活性腐殖质的4％以下，对整体影响较小；

128、包括土壤有机物的其余部分的腐殖质库hum；

129、包含植物残体、农家肥和泥浆的有机物库fom；

130、三个氮库的所有分解率均采用式37至39的一阶反应动力学进行计算：

131、

132、

133、

134、式中，ci为按fom、bom和hum的顺序排列的第i个碳库的碳含量；ki为一阶分解率系数，该系数根据所考虑的非生物因素进行修正；αi为碳周转胁迫系数；fi为碳流分配系数；ei为基质利用效率。

135、土壤温度和湿度是影响上述碳氮周转的两个主要的非生物因素，而土壤粘粒含量是影响hum的一个特殊因素，通过式下式(55)-(56)计算得到影响因子；假定不同因子之间不存在交互作用，综合效应为乘法效应；允许有机质的分解过程中碳在不同子库之间流动；则有：

136、α1＝α2＝fm(t)fm(θ)  (55)

137、α3＝fm(t)fm(θ)fm(clay)  (56)

138、胁迫方程fm的计算方式如下：

139、

140、

141、

142、式中，θ和θs分别为土层中实际体积含水量和饱和体积含水量(cm3cm-3)；θwp，θfc和θ25分别为作物枯萎时的土壤体积含水量、田间持水量和θwp+0.25(θfc-θwp)；cc为土层中粘土含量(g g-1)；cc′设为0.25kg kg-1，a取2.0；

143、不同碳氮库之间的流量以碳为单位计算，相应的氮通量取决于接收库的碳氮比；单一土层中的每个有机物库都有一个特定的无量纲的碳氮比cnri值，氮的净矿化方程计算如下：

144、

145、式中，dci是第i层土壤中每个库的碳分解含量(mg l-1)。

146、所述土壤水模块中灌溉制度的确定原则为：

147、灌溉事件是通过控制灌溉时间和灌溉深度来确定；

148、灌溉时间由不同的触发条件决定；

149、灌溉事件的触发机制包括：临界田面积水深度、临界土壤水张力、临界土壤含水量和田面积水层消失的天数；

150、灌水的触发条件由用户自主设置，随着作物生长发育阶段变化；

151、灌溉深度(irad，cm)根据作物发育阶段进行计算确定，计算公式为：

152、

153、式中，i是土壤剖面中从表层至底层的土层序号，n是土壤剖面中位于犁底层上方的土层数，δd是土层的厚度(cm)，θ是实际土壤体积含水量(cm3cm-3)，θs是饱和土壤体积含水量(cm3cm-3)，t是模拟运行的天数(d)，x是用户自定义的灌溉后积水深度，即积水深度的上限(cm)。

154、本发明具有如下有益效果，

155、与现有技术相比，本发明技术方案提出的适用于寒区粳稻农田的农业水文-生态多过程耦合模拟模型及方法，开发了适用于农田水热运移-碳氮转化-作物生长-产量形成多过程耦合的水文-生态系统模型，并构建了适用于不同灌溉模式的灌溉制度决策模块，为水氮高效管理提供了精准的定量工具；

156、与现有技术相比，本发明技术方案阐明了复杂灌排环境下稻田土壤-作物系统的水分利用途径和氮素归趋，明晰了寒区水稻潜在产量的影响因素及稻田节水潜力。