使用HYDRUS-1D与PEST耦合模型估算细粒土壤的水力参数

《Agricultural Water Management》：Estimation of soil hydraulic parameters for fine-textured soils using HYDRUS-1D coupled with PEST

农业水资源管理需要可靠评估土壤水力参数。本研究利用HYDRUS-1D耦合PEST工具，对曼尼托巴省细粒土在2016-2019年四个生长季的水力参数进行反演，发现表层饱和水力传导率年际波动显著（120 cm/d至10 cm/d），深层稳定（<20 cm/d），模型性能因气候波动而受影响，需考虑土壤-植物-气候相互作用优化管理策略。

　　本研究探讨了在曼尼托巴省卡曼-埃尔姆克里克地区，对一种细粒土壤进行农业水资源管理策略开发所需土壤水力参数的稳健评估。通过将HYDRUS-1D与参数估计（PEST）软件结合，研究者对四个生长季节（2016-2019年）的土壤水力参数进行了估算。为了进行校准，每个生长季节生成了100组初始参数集，利用现场测量的土壤体积含水量数据，分别在0-10厘米、10-30厘米、30-70厘米和70-130厘米的土壤层中进行分析。HYDRUS-1D与PEST的结合揭示了解决方案的非唯一性，因为多个参数集能够产生相似的模型表现。通过随机化初始参数集进行连续校准以及基于性能的筛选，研究者识别出了稳定且最大似然的参数范围。模型表现基于纳什-苏特克利夫效率系数（NSE）（0.65-0.92）、百分比偏差（PBIAS）（-9.31%-1.48%）、均方根误差与标准差的比率（RSR）（0.29-0.59）以及决定系数（R2）（0.68-0.96），在所有季节中均表现良好，除了2018年的验证期，由于该季节校准和验证期间的高变异性，模型表现未达到预期。饱和水力导度在表层（0-10厘米）显示出显著的年际变化（分别为2016年、2017年、2018年和2019年的120、10、21和120厘米/天），而在更深的土壤层中则相对稳定（<20厘米/天）。这些发现强调了在表征土壤水力参数时考虑土壤-植物-天气相互作用的重要性，以提高农业水资源管理策略的效率，并支持寒冷气候种植系统中的适应性。

曼尼托巴省以其高质量的谷物、油料作物和农业食品产品闻名，这些产品满足了不断增长的市场需求。农业和农业食品行业对省内的生产总值（GDP）贡献了大约7%。考虑到该行业对曼尼托巴省经济和本地及全球粮食安全的支撑作用，必须采取策略确保土壤能够有效支持可持续的作物生产，以及其他生态系统服务。曼尼托巴省的近70%（约470万公顷）的农田用于作物生产。这些土壤的生产力面临威胁，因为极端气候条件和集约农业实践导致的土壤结构变形，会抑制所有支持作物生产和其他生态系统服务的土壤功能和过程。

土壤结构变形导致了在管理极端水分条件时的重大挑战。过多或过少的水分含量在该省造成了超过一半的作物损失。在2022年，由于冬季积雪过多和春季降雨量大，大约有348,000公顷的耕地未能播种，造成了估计的过剩水分保险赔付金额为5350万美元。相反，在严重干旱条件下，曼尼托巴省的作物保险赔付金额显著增加，从2020年的6200万美元增长到2021年的46900万美元。这些损失突显了改进水资源管理策略的紧迫性，以帮助生产者应对洪水和干旱。

应对这些水资源挑战需要更好地理解土壤结构和水力参数及其对地表和地下土壤水分动态的影响。传统的水力参数测定方法通常繁琐、耗时、昂贵，并且往往无法代表田间条件。作为替代，水文模型被广泛用于模拟土壤水分动态和使用反演建模估算水力参数。已经开发了几种不同复杂度的模型，包括土壤-大气-植物（SWAP）模型、土壤-植被-大气传输（SVAT）模型、农业管理系统地下水负荷效应（GLEAMS）模型，以及HYDRUS包。

HYDRUS包是模拟实际条件下土壤水分动态最常用的水文模型之一，因为它提供了对理查兹方程的稳健、物理基础的解决方案。理查兹方程通过将非饱和流的达西定律与质量守恒的连续性方程结合，描述了土壤水分在变湿条件下流动的机制。HYDRUS需要输入数据，包括气象条件、模型域的边界条件、土壤物理条件、土地利用（植被类型、根系分布）以及土壤水力参数（饱和水力导度、土壤水分保持曲线及其影响该曲线的属性）。尽管需要输入参数，但物理参数、初始参数和边界条件之间的相互作用仍然需要模型校准，以确保模拟能够代表田间和水文条件。然而，由于水文条件的非线性，一个以上的参数集可以再现观测数据，从而导致解决方案的非唯一性。通常的做法是选择参数集，使其具有最佳的统计拟合（例如，最小化目标函数，最大化统计评估参数的值）。最佳参数集并不总是能准确代表田间条件，因为多个具有完全不同参数集的模型结果可以满足统计参数集的性能评估阈值。

为了更准确地解释模型性能，采用概率方法（如最大似然方法）呈现结果是一个更好的选择。HYDRUS-1D/2D/3D可以利用反演建模技术，通过已知的初始参数估计土壤水力参数。在本研究中，HYDRUS-1D的反演建模使用了ROSETTA土壤传输函数，该函数基于土壤质地、土壤密度和特定基质势下的水分含量等容易获得的土壤属性预测van Genuchten-Mualem参数。在本研究中，ROSETTA被用于其分层模式，结合了质地数据（沙、粉沙、黏土百分比）和土壤密度，生成了后续反演优化的初始参数估计。反演建模的步骤如（表3）所示。步骤1-5在反演建模过程中通过图形用户界面（GUI）完成。该过程定义了目标函数（φ），量化了模拟数据与观测数据之间的差异，并使用了Levenberg-Marquardt非线性最小化方法来最小化φ。平均土壤体积含水量（表2.2.1）被提供为输入数据，以校准模拟。初始土壤水分含量剖面（t=0）使用了模拟期间前一天（t=1）的传感器观测值定义。通过统计性能指标评估了每个层的模拟数据与观测数据之间的吻合度，如（表2.4.4）所述。在反演建模中，最多可以同时优化15个参数。然而，在本研究中，需要优化24个参数（4层×6水力参数），这使得反演建模方法具有局限性。因此，参数估计工具PEST被用作外部校准工具，以支持校准过程。

PEST利用高斯-马夸尔特-莱文伯格优化方法。创建参数估计过程和使用HYDRUS与PEST结合生成参数集所需的步骤如（表3）所示。水力参数的上下限来自文献，并根据模型收敛性进行了调整。现场测量的土壤水力导度通常超过实验室测量值，特别是在黏土质地的土壤中。因此，表层（0-10厘米）的Ks上限被设定为120厘米/天，这与黏土土壤的报告值（大约110厘米/天）一致，高于文献中报告的典型平均值，确保了Ks的潜在田间变化得到充分捕捉。

HYDRUS与PEST的耦合能够优化大量参数。生成的关联矩阵有助于识别具有强相关性的参数，从而支持决定哪些参数需要固定。本节简要介绍了HYDRUS-PEST的耦合方法。Doherty等（2010）提供了对模型-PEST接口的详细描述。按照（表3）中描述的方法，创建了连接可执行HYDRUS和PEST的文件。为了运行，PEST需要三个基本文件：模板、指令和输入控制文件。Selector.in（HYDRUS模型输入文件之一）被用作创建模板文件的基础，因为这是向PEST介绍需要优化的参数的文件。指令文件是基于模型输出文件创建的，该文件向PEST传达了需要读取并比较的输出值。这些观测初始值通过输入控制文件提供给PEST。

该控制文件还包含模板、指令、模型输入和输出文件的名称，模型执行命令，校准终止和参数转换标准。基于参数转换标准和参数界限，PEST反复运行模型，直到优化完成或终止标准被满足。优化后的参数、模拟值与观测值之间的相关系数、目标函数以及每次优化的模型运行次数被记录在由PEST生成的记录（.rec）文件中。经过100次成功的优化后，使用（表2.4.4）中的参数评估了模型性能。模型运行存在不可接受的值和收敛问题，这些模型被拒绝。剩余的校准模型参数集被分析以识别每个参数的最可能范围。对于每个参数，初始搜索范围被划分为更小的区间，并计算校准值落在每个区间内的相对频率。频率被解释为参数在该子范围内的概率。显示最高频率（即最大似然）的区间被视为最可能的范围，并被采用为最终参数估计的优化边界。

例如，表层（0-10厘米）的Ks初始范围（0-120厘米/天）被划分为六个区间：0-10、10-30、30-50、50-70、70-90和90-120厘米/天。计算了每个区间内接受的校准运行比例，并使用显示最大似然的区间来建立最终参数估计的边界。表层Ks值的初始范围是根据观察到的Ks值进行调整的，而深层Ks值则反映限制性亚土壤条件。

在本研究中，通过将HYDRUS与PEST结合，对2016年数据进行了参数估计。校准期间和验证期间的模型性能评估结果如（表8）所示。然而，优化参数集之间的差异表明解决方案的非唯一性（多个参数组合可以在相同条件下模拟相似结果），这引发了对校准参数代表实际田间条件可靠性的担忧。为了进一步探索水力参数之间的补偿性相互作用，生成了关联矩阵。仅比较了θr、θs和Ks，因为它们是物理可解释且与田间应用相关的参数。形状参数α和n被排除，因为它们是经验性的，对相互作用效应高度敏感，倾向于在拟合过程中进行补偿。

表9显示了校准的van Genuchten-Mualem水力参数之间的相关系数。在四个土层中的三个层中，θr与n之间表现出强相关性（≈0.78-0.95）。这些相关性表明土壤水分含量和土壤水分保持曲线形状参数之间存在显著的补偿性相互作用。为了减少内部补偿并提高稳定性，手动固定了θr，允许剩余参数的更可靠和物理意义的校准。θr的值使用了由科罗拉多大学开发的土壤水力性质计算器进行固定。与永久萎蔫点不同，θr是一个模型参数，反映了在无限高吸力下土壤孔隙中保留的水分量，并且通常低于基于植物定义的萎蔫阈值。根据砂、粉沙和黏土百分比作为输入，将土壤水分含量与无限高基质势（-50000米）下的最低水分含量设置为θr。最终计算的θr值如（表10-13）所示。

表10-13显示了最终参数估计运行中基于最大似然范围的参数控制设置。这些表格展示了不同深度下参数的最终范围，包括θr、θs、α、n、Ks和l。通过将参数空间缩小到高概率范围，最终估计的参数变得更加稳定和具有物理意义。深层的低Ks值（<10厘米/天）与之前对缓慢排水和可能的滞留条件的观察一致。

在2016年、2017年和2019年，模型性能在校准和验证期间均表现出色，但在2018年的验证期间表现较差，这可能是因为校准和验证期间的高季节内变异性。2018年的验证期间性能较差，可能与校准期间的天气条件变化有关。例如，在2018年，降水强度和持续时间的剧烈变化可能改变了表层的条件和土壤的物理响应，从而偏离了模型最初校准时的强迫条件。表层Ks的显著年际变化（2016年和2019年约为120厘米/天，而2017年和2018年则降低至10厘米/天）表明，天气、耕作模式和表面条件在控制入渗动态方面起着关键作用。深层的Ks值则相对稳定（<20厘米/天）。

本研究的结果表明，土壤-植物-天气相互作用对土壤水力参数的表征具有重要意义。土壤水分的动态变化受到这些相互作用的强烈影响，因此在制定农业水资源管理策略时必须加以考虑。研究结果强调了在表征土壤水力参数时考虑土壤-植物-天气相互作用的重要性，以提高农业水资源管理策略的效率，并支持寒冷气候种植系统中的适应性。

在进行连续校准后，研究者采用了基于性能的筛选方法，以识别稳定且最大似然的参数范围。通过使用100组随机初始参数集进行重复模型运行，结合蒙特卡洛采样方法，进一步优化了校准模型参数。通过将参数范围划分为更小的区间，并计算每个区间内接受的校准运行比例，研究者确定了每个参数的最可能范围。这些参数范围的确定为最终参数估计提供了数据驱动的基础。例如，表层Ks值的初始范围（0-120厘米/天）被划分为六个区间，计算了每个区间内接受的校准运行比例，并使用显示最大似然的区间来建立最终参数估计的边界。

最终参数估计的运行使用了HYDRUS-PEST模型，以提高参数估计的可靠性。通过使用重复运行中最可能的参数值作为初始参数集，研究者进一步优化了模型。参数范围的缩小使得最终估计的参数更加稳定和具有物理意义。深层的低Ks值（<10厘米/天）与之前对缓慢排水和可能的滞留条件的观察一致。此外，研究还指出，未来的研究应考虑使用多孔性模型（如双孔性模型）来更好地代表这些动态，以捕捉宏观孔隙流与土壤基质之间的相互作用。

本研究的结果不仅揭示了土壤水力参数的非唯一性，还展示了这些参数在不同年份和不同土壤层中的变化趋势。这种变化趋势对于农业水资源管理策略的制定具有重要意义，因为土壤水力参数的年际和层际变化能够影响水的流动和储存。因此，研究强调了在寒冷气候条件下，土壤水力参数的季节性和区域特定性的重要性。通过使用现场数据和模型的结合，可以更准确地评估土壤水力参数，并制定更有效的水资源管理策略。这些策略不仅需要考虑土壤的物理特性，还需要结合气候和作物管理实践的变化。