中国西北干旱半干区生境恶劣，常年干旱少雨，植被稀疏，风沙严重。国家为了改善生态环境，近几十年，陆续开展了包括退耕还林还草、风沙源及荒漠化治理等一系列造林工程^［1］，以此达到减少风蚀沙化土地及增加森林覆盖率等目的。工程实施后，大面积植被恢复是否会造成土壤干燥化及地下水的过度消耗等一系列问题产生了巨大的争议，已有研究发现大模工程的实施导致了土壤水分亏缺^［2］。土壤水分的变化一直是干旱半干旱区生态系统水循环过程中的研究重点，成为土壤特征研究中重要的指标^［3］，土壤水分在干旱半干旱区是限制植被生长及影响群落结构稳定性的主要因子。因此，了解该区土壤水分变化过程，探究土壤水分盈亏状况，对干旱半干旱区生态系统的稳定与生态修复具有重要意义。

目前，在干旱半干旱区对土壤水分稳定性^［4-6］、土壤水分补给^［7-8］和土壤水分利用^［9］等方面已取得一定的研究成果。植被恢复后需要掌握其水文动态变化过程，主要是掌握水循环过程中各变量之间的分配方式，连续并全面地对每个水文过程进行动态监测难度较大，因此通过数学模型模拟水文过程、定量评价水文系统中土壤水分运移转化关系已经成为生态保护的主要内容^［10-12］。许多研究者已通过不同的土壤水分运移模型对干旱半干旱区的水分运移过程进行了系统地研究，Hydrus、Philip、SWAT等数学模型已在研究土壤水分分布变化中被广泛应用^［13-15］。

Hydrus-1D是基于Richards原理，同步将大气降水、植被吸水、地下水变化及土壤理化性质考虑在内的垂直水分循环过程，适用于全部边界条件模拟^［16-18］，编辑过程相对简单且应用过程较为灵活，同时具备参数优化功能^［19-21］，使得Hydrus-1D模型在不同的土地类型及不同植被覆盖中的土壤水分变化、溶质运移及土壤水分深层渗漏等方面进行应用^［22-23］。王宇祥等^［18］ 利用Hydrus-1D模型对沙区不同类型的沙丘与草甸区的土壤水分变化进行了模拟并进行了模型的适用性评价，以此计算出土壤深层渗漏对降雨的响应。王国帅等^［24］研究了荒漠绿洲的水盐动态过程，揭示了土壤水分的亏缺状态，以此评价其农业生态状况。李冰冰等^［22］对渭北旱塬农地及苹果园深度10 m的土壤水分时空变化进行了模拟，为研究深层土壤水分变化找到了可靠方法。毛乌素沙地腹地以沙丘为主，主要成分为细沙和粉沙粒，因此土壤水分运动过程中基本不受到团聚体的影响，更有利于通过模型准确描述各项水土保持工程实施后土壤水分变化、水分深层渗漏及地下水补给情况。

本文选取毛乌素沙地腹地典型植被恢复区杨柴（Hedysarum laeve）人工林为研究对象，基于2018年5—9月土壤深度10、30、50、70、90、110 cm处水分含量的连续实测值、气象连续观测数值、土壤水分特征曲线参数等数据，利用Hydrus-1D模型对土壤水分进行模拟，并对其适用性进行评价，以期为沙区生态保护、植被修复、水资源管理与应用等方面研究提供科学的研究方法及依据。

研究区位于鄂尔多斯市乌审旗乌兰陶勒盖治沙站内（39°08′17″N，109°31′36″E），地处毛乌素沙地的南缘，属于典型的温带大陆性气候，冬长夏短，四季分明，年降水量270—350 mm，每年7—9月降水量较大，占年降水量的60%—70%，年蒸发量2 100—2 600 mm，平均气温6—7 ℃，常年风沙较多。风沙土是研究区分布最广泛的土壤，其次是栗钙土、草甸土等。植被类型多样，主要的沙生灌木有杨柴、沙柳（Salix psammophila）、油蒿（Artemisia ordosica）等，草本层植物主要为虫实（Corispermum hyssopifolium）、沙米（Agriophyllum squarrosum）、雾冰藜（Bassia dasyphylla）等^［25］。

自2015年开始，在毛乌素沙地2002年进行飞播造林的典型杨柴人工林灌丛下布设2台土壤水分自动监测系统（Item 3345WD2 Watch Dog 2800 Weather Station，美国，图1），每台监测仪配备6个土壤水分含量自动监测传感器，分别布设在杨柴林内土壤深度为10、30、50、70、90、110 cm处，每30 min分别对不同深度土壤水分含量进行自动监测、记录并存储数据，本文选取2018年数据对其进行研究分析。同年在样地内安装1套气象监测系统（HOBO U30便携式气象站，美国，Onset公司），对环境中大气温度、太阳辐射、相对湿度和降水量等进行监测，每隔30 min监测系统自动记录并存储1次数据。

2018年8月，利用环刀采样法在杨柴人工林样地内挖取6个100 cm×100 cm×120 cm土壤剖面，在深度分别为0—20、20—40、40—60、60—80、80—100、100—120 cm土层取样，每层取样3个重复。土壤取样后带回实验室通过压力膜仪对土壤水分特征曲线进行测定，利用浸水实验对土壤饱和含水率和土壤容重进行测定，利用马氏瓶渗透仪对土壤饱和导水率进行测定，利用CLY-2000系列激光粒径分布仪（CLY-2000，国产）测定土壤粒级组成。土壤粒径分级标准采用美国农业部USDA制，黏粒（<0.002 mm）、粉粒（0.002—0.050 mm）、砂粒（0.05—2.00 mm）。测量数据平均值均用于构建模型（表1）。杨柴灌木林土壤质地以砂粒为主，每层砂粒含量都在94.41%以上。不同土层黏粒和粉粒比例不同，30—50 cm和90—110 cm土层黏粒较其他土层多，分别为4.95%和5.7%，50—70 cm土层粉粒较比其他土层占比高，为1.02%。采用压力膜仪测定出0—20、20—40、40—60、60—80、80—100、100—120 cm深度（分别对应10、30、50、70、90、110 cm深度）土层不同压力下对应的土壤含水量，导入vanfit土壤持水曲线拟合程序进行土壤水分特征参数拟合，利用实测的土壤颗粒粒径百分比和土壤干容重通过Hydrus-1D模型参数求解功能模拟各层van Genuchten特征指标，之后利用实测的土壤含水率、饱和含水率和导水率对参数进行校正，最终获取杨柴样地土壤的基本物理性质及各水力学参数（表1）。

潜在蒸发量采用Penman-Monteith公式计算，通过气象站记录的日气象数据，带入到公式中计算出日参考潜在蒸散量（ET₀）。

在2018年5—9月，每10天进行1次植被调查。釆用了SunScan冠层分析系统（Delta-T Devices SunScan，英国）对杨柴人工林叶面积指数（LAI）进行调查。采用挖掘法对杨柴的根系密度进行调查，挖取1 m×1 m×1 m的土壤剖面，分别对0—20、20—40、40—60、60—80、80—100 cm深度土层处进行根系取样，带回实验室处理得出。叶面积指数及根密度均为模型所需参数^［17］。

本文采用经典的Richards方程^［26］描述土壤水分运移过程，仅考虑垂向运动。

式中：t为时间（h）；θ为土壤体积含水率（cm³·cm^-3）；h为负压水头（cm）；z为垂向空间坐标（cm）；K（h）为非饱和导水率（cm·h^-1）；q₀（0，t）为上边界水分通量（cm·h^-1）；h（B，t）为下边界负压值（cm）；s（z，t）为根系吸水速率（cm·h^-1）。

土壤水分特征曲线和非饱和导水率利用van Genuchten方程拟合^［27］：

式中：θ_s为土壤饱和含水率（cm³·cm^-3）；θ_r为土壤残余含水率（cm³·cm^-3）；K_s为饱和导水率（cm·h^-1）；S_e为无量纲的相对含水量；α、m、n、l均为拟合参数。

通过Hydrus-1D模型模拟杨柴根系在不同土层z处吸水的值，表达式如下^［28］：

式中：α（h，t）为水分胁迫函数；β（z）为归一化根系吸水的分布函数；T_p为植被日潜在蒸发量，cm。

选取2018年5月18日至9月16日为研究期，把杨柴灌木林模拟的土壤剖面设置为6层，采用时间离散化的方式进行以小时为单元进行数值模拟，共模拟2 928 h，依据收敛的迭代次数逐渐调整时间离散的间隔^［22］。本模拟中起始时间设置为0，模拟结束时间为2 928 h。初始时间步长为0.01 h，最小时间步长为0.01 h，最大时间步长为1 h。

本文中模拟的杨柴样地地下水埋深均较大，且沙土渗透性强，所以下边界统一设置为自由排水边界，由于样地内均无径流产生，上边界都设置为无径流的大气边界。

初始土壤水力参数中残余含水率、饱和含水率、饱和导水率、α、n、l、根系吸水等模型参数是将土壤干容重及砂粒、粉粒、黏粒含量代入Rosetta公式进行预测而得（表1），其中饱和含水率、土壤干容重、饱和导水率属于实测值。水分胁迫函数是根据实测值叶面积指数、根系密度和ET₀代入模型中自带的Feddes函数模拟根系吸水过程获取^［28］。

为了确定模拟杨柴人工林土壤水分动态过程的准确性，利用决定系数R²、平均误差ME、均方根误差RMSE和平均相对误差MRE对模拟精度和误差进行评价，其中，R²的大小决定了模拟值与实测值的密切程度，ME决定了偏离程度，RMSE和MRE趋向于0即模拟值越接近于实测值，计算公式如下：

式中：yip为土壤含水率的实测值；yim为土壤含水率的模拟值；N为样本点个数。

利用Excel进行作图，土壤水分实测值和模拟值的拟合及模型各个参数的准确性采用SAS 9.4软件进行验证。

通过优化参数后Hydrus-1D模型模拟水分总体变化（图2），杨柴180 cm层土壤水分动态变化在模拟阶段的决定系数为0.587，在验证阶段为0.547。平均误差在模拟阶段和验证阶段分别为-0.0053和0.0087，均方根误差在模拟阶段和验证阶段分别为0.00017和0.00038（表2）。

选取6—8月每月首日没有降雨影响下的土壤垂直水分来评估模型的模拟能力（图3）。可以看出杨柴群落土壤水分含量的垂直方向呈先增大后减小再增大再减小的趋势，呈“M”状变化。7月1日土壤含水量先在30 cm处增大到0.038 cm³·cm^-3左右后减小，在90 cm处又增大到0.055 cm³·cm^-3左右后减小。6月1日和8月1日土壤水分在校准期与验证期间的实测值与模拟值吻合度较高，但二者之间在不同深度土层存在着一定的偏差，主要在10、50、110 cm处。

利用2018年实测的不同深度土壤水分数值对利用优化参数后Hydrus-1D模型结果进行对比（图4）。可以看出不同深度土壤水分含量的模拟值与实测值吻合度较高，模拟的结果总体上说可以较准确地反映出不同深度土壤水分的变化。10 cm和110 cm土层深度土壤水分含量模拟值稍高于实测值，30、50、70、90 cm土层深度土壤水分含量模拟值与实测值几乎重合，总体上看Hydrus-1D模型的结果与实测值拟合效果较好，相关性强。

采用不同深度土壤水分的实测值数据对所建Hydrus-1D模型进行验证（图5），验证结果中不同深度土壤水分的实测值与模拟值总体分布较集中，分布于1∶1线上或两侧，总体模拟效果较理想。110 cm深度土壤水分建模和检验阶段拟合曲线平行偏离1∶1线，土壤水分含量模拟值被高估。ME在70 cm深度建模阶段和90 cm深度检验阶段数值均大于0，土壤水分含量略小于实测值，其他深度土壤水分含量在建模和检验阶段数值均小于0。建模阶段RMSE在0.0040%—0.0217%，检验阶段RMSE在0.0055%—0.0267%。建模阶段除了110 cm深度土壤，其他深度土壤MRE在0.1490%—0.5225%，检验阶段除110 cm深度土壤，其他深度土壤MRE在0.1518%—0.6655%（表3）。除10 cm和110 cm深度的模拟阶段R²小于0.5外，其他的建模和检验阶段R²均大于0.5。总体而言，率定后的Hydrus-1D模型对毛乌素沙地杨柴灌木林不同深度土壤水分变化模拟较理想。

在降水补给杨柴人工林土壤过程中，实测值与模拟值变化具有异同（表4）。10、30、50、70 cm土层土壤水分入渗开始到达峰值时，土壤水分最大增量的差异较小，模拟效果较好，但90 cm深度模拟值小于实测值，110 cm处模拟值大于实测值。土壤水分入渗的平均速率模拟值总体小于实测值，模拟入渗过程所需时间高于实测值。从入渗深度看，模拟值与实测值入渗深度相同，都随着降水量增大，入渗深度增加。

降水是干旱半干旱地区的主要水分来源，直接影响地表水及地下水分的补给，从而决定植物群落的稳定性，因此掌握土壤水分动态变化及沙地水循环过程中土壤水分的深层渗漏量尤为重要。Hydrus-1D作为一套用于模拟饱和多孔介质中水分、能量、溶质运移的数值模型，经过改进和完善，已经在多地进行模拟与验证并得到了广泛的认可与应用^{［13，18，29-30］}，在一定程度上可以很好地模拟土壤水分及地下水分量随时间的变化趋势，为研究植物群落中土壤水分变化及地下水补给提供了有效方法。为此本实验选用优化参数后的Hydrus-1D模型对毛乌素沙地杨柴样地土壤水分运移进行了研究，模拟毛乌素0—110 cm剖面土壤的水分变化情况，探讨通过该模型与方法模拟土壤水分变化的适用性。基于2018年5—9月的实测数据进行比较，发现该模型适用于毛乌素沙地飞播杨柴灌木林，能较好地模拟土壤剖面水分的时间及空间变化特征，为进一步研究毛乌素沙地整体林水关系提供了重要的方法。

不同土层水分含量在垂直方向上有着不同的分布特征，表层土壤波动最为频繁，30—90 cm土壤波动强烈，110 cm土层水分含量较稳定，通过参数输入和优化后的水力学参数，30—90 cm土层模拟出的结果在校准期及验证期与实测值总体吻合度较好，但土壤表层10 cm和底层110 cm的实测数据与模拟值存在较大偏差。上层土壤的环境更复杂，受到蒸发、降水入渗和土壤理化性状的影响土壤水分变化较大^［17-18］，模型难以完全把各项影响因素都考虑其中，因此预测结果略有偏差。杨柴群落土壤110 cm处属于干沙层^［31］，自身水分含量较低，对降水入渗响应较小，在降水量大于20 mm时才有较小的波动，而且马上又恢复到一个较小且稳定数值，而Hydrus-1D在模拟过程中未能考虑到沙地水分分布规律的特殊性，而主要依靠土壤水力学参数和根系分布特征来进行模拟预测^［32］，因此该土层的偏差较大。可以看出在校准期不同时间段模拟的效果也存在偏差，6月拟合的效果优于7月，因为水分入渗土壤的过程相对复杂，主要受到土壤物理结构、养分含量及土壤水分初始含量等因素的影响，同时在植被群落中分布的空间异质性及瞬时气象变化都会影响模型的精准性，因此不会出现完全拟合或者瞬时拟合不紧密的结果^［23］。对降水入渗土壤过程进行模拟可以发现，入渗增量差值及深度差异性不明显，但是在入渗速率上模拟值明显小于实测值，模拟入渗过程比实际入渗过程更缓慢，模拟过程研究效果不理想，但最终降水对土壤水分补给的增量及补给深度模拟效果较为理想，可以较准确计算出土壤水分补给量。90 cm深度的降水补给模拟值小于实测值，该样地前期已做过相关研究，该处受到壤中流的影响^［31］，在降水可以入渗到该深度时水分会侧向运移到该土层。110 cm深度降水入渗差值的模拟值明显小于实测值，因为该深度土层属于干沙层，土壤水分数值总体偏小且对水分补给的响应过程较弱，而Hydrus-1D并未将该因素考虑在内，因此模拟入渗过程偏差较大。李冰冰等^［22］利用Hydrus-1D模型模拟渭北旱塬土壤水分垂直分布时得出模型决定系数为0.65—0.85，均方根误差0.01%—0.02%，段良霞等^［10］在半干旱区对不同地貌的土壤水分变化进行模拟，得出均方根误差为0.018%—0.029%，本文的研究结果与其相似，说明Hydrus-1D模型对研究区土壤水分动态能较好地预测。

利用优化土壤水力学参数后的Hydrus-1D模型对毛乌素沙地杨柴样地2018年5—9月的土壤水分进行模拟，除表层和底部拟合效果较差外，其他土层决定系数0.55—0.78，整体上模型可以较好地描述研究区土壤水分动态变化规律。