江西武夷山国家级自然保护区典型林地类型土壤水力特性差异及土壤水分调控效应模拟

徐燕星, 刘 昭, 涂安国, 雷 平, 胡 皓 ,赵佳鼎,石芬芬

(1.江西水利职业学院，330013，南昌；2.江西省水利科学院，江西省土壤侵蚀与防治重点实验室，330029，南昌；3.江西武夷山国家级自然保护区管理局，334500，江西，上饶)

0 引言

土壤水是森林生态系统物质和能量循环的关键载体，但同时森林土壤水分消耗、存储和转化又受到森林植被生长的影响[1]。有关研究表明，植被活体根系穿插以及枯落物腐解能够改善土壤的理化性质[2]，不同林地类型下土壤孔隙率、入渗性能等土壤水分物理性质差异明显[3-5]。虽然目前已有大量关于森林土壤水分的研究成果,但关于不同林地类型下土壤水力特性差异及其产生的土壤水分调控效应研究还不够深入，对水源是起到了涵养作用还是消耗作用的问题亟需回答。

森林土壤水分受降雨、径流、蒸发、根系吸水、植物截留等多个水文因素共同影响，为评估不同林地类型土壤水力特性这一影响因素差异下的土壤水分调控效应，建立土壤水分运动模型是可行的办法[6]。目前描述土壤水分运动模型主要分为Richards模型、动力波模型和均衡模型，其中以达西定律和质量守恒定律推导的Richards 方程具有坚实的物理基础，对于水分通量计算具有较大的优势，在各种水文模型中得到了广泛应用[7]。但是，由于该方程是高度非线性的抛物线型偏微分方程，其数值求解方法具有一定挑战性，一直是热点研究领域[8]。其中Ross方法是一种求解Richards 方程的主变量切换技术，在模型收敛性、计算效率和计算精度上均具有优势[6, 9]，这为分析土壤水分运动过程提供了有效工具。

1 研究区概况

江西武夷山国家级自然保护区地处江西省东北部的铅山县境内(117°39′30″～117°55′47″E， 27°48′11″～28°00′35″N)，坐落于武夷山脉的主峰(黄岗山，海拔2 160.8 m)区域，是我国东南陆域最高山地的自然保护区，总面积16 007 hm2，其中核心区面积4 835 hm2，占保护区面积的30.2%；缓冲区面积2 021 hm2，占保护区面积的12.6%；实验区面积9 151 hm2，占保护区面积的57.2%。区内保存有世界同纬度现存最完整、最典型、最原始的中亚热带原生性中山森林生态系统，是珍贵树种群落和珍稀野生动植物的自然集中分布区。

叶家厂小流域为本研究中的典型研究区域，主要植被类型分布为：阔叶林52.2%，针阔混交林21.4%，针叶林10.6%，毛竹林10.4%，中山矮曲林和草甸5.4%；而在该流域河岸带范围内，毛竹林面积可达50%以上，因此毛竹林极具代表性。本研究在叶家厂小流域范围内选取了毛竹林、杉木林、檫木林和柳杉林四种林地类型作为典型林地类型。其中杉木林是针叶林的典型代表，檫木林(伴生少量杉木)是阔叶林(针阔混交林)的典型代表，柳杉林是古树林的典型代表。于2018年11月在该4种典型林地类型中各选取20 m×20 m的样方进行植被调查，毛竹、檫木、杉木和柳杉分别代表各典型植物群落中的优势树种，其主要特征如表 1所示。

表1 各样方优势树种的主要特征

2 研究方法

2.1 样品采集与指标测定

典型林地类型样方土层较薄，平均土壤厚度仅35～47 cm，10 cm以上为腐殖土层。以100 cm3环刀和压力膜仪专用环刀在各样方随机采集0～10 cm和10～40 cm的原状土样(3次重复)。采用烘干法测试土壤容重，采用饱和入渗试验测试饱和水力传导度，采用压力膜仪分别测试了各样品在10 cm、20 cm、30 cm、40 cm、60 cm、80 cm、100 cm、200 cm、300 cm、500 cm、1 000 cm、1 500 cm等11组土壤负压的土壤含水量，以获取土壤水分特征曲线。此外，在样方附近区域选取裸地、茶园、水田3片样地作为参照，同步采集和测试上述土壤物理指标。其中，裸地位于保护区内靠近叶家厂小流域，已风化成土，但无任何植被覆盖；茶园位于保护区内靠近武夷山镇簧村，有较强的人类活动；水田位于保护区外武夷山镇往铅山县方向公路旁。

2.2 数据分析方法

2.2.1 土壤水动力参数 本文采用van Genuchten-Mualem 模型描述土壤体积含水率θ(cm3/cm3)、土壤负压h(m)和土壤水力传导度K(θ) (m/d)之间的函数关系，即

(1)

(2)

式中：α为土壤介质进气值的倒数(cm-1)，一般土壤粘性越强，该值越小；m和n是与土壤颗粒大小分布相关的参数(无量纲)； 满足m=1-1/n，且0

2.2.2 土壤水运动模型 本文主要是讨论各典型植物群落样方和各参照样地的土壤水动力参数的不同条件下，土壤水分运动过程对于降雨入渗和土面蒸发的响应。在建模之前，需对一些条件进行限定和假设。为简化起见，采用一维Richards方程模拟入渗和蒸发过程。土层底板为坚硬基岩，下边界可假设为不透水层(隔水边界)。由于该区平均土壤厚度35～47 cm，因此可统一假定各处理的土层厚度为40 cm，分为0～10 cm和10～40 cm的2层。初始含水率剖面均设置为底部节点刚好饱和时的稳定剖面(各点重力势与负压之和处处相同，当无蒸发或降雨时，含水率剖面稳定)。模拟土壤水库从上边界的补排过程，设置上边界条件为虚拟的大气边界，设为前20日蒸发后降雨，参考作物腾发量强度根据经验值设为2 mm/d，净降雨强度设为30 mm/d(根据叶家厂雨量站2017―2019年观测，月内平均降雨强度为2.7～37.9 mm/d)；由于林地坡度较陡，且土壤导水性好，可认为底板处饱和后，水迅速以饱和流向下坡运移，因此以底部节点刚好饱和时，作为降雨结束时间，及模拟时段内不产生地表径流。

综上，土壤水运动控制方程表示为：

(3)

其中上边界条件可表示为

(4)

下边界可考虑为隔水边界，即

(5)

初始条件设置为

(6)

式(3)～式(6)中，z为垂向坐标(m)，向下为正；zb为土体下边界深度，本文取为0.4 m；t为模拟时间(d)；β为饱和-非饱和状态参数，在饱和时取为1，非饱和时为0；μs为弹性释水系数，本文主要考虑非饱和土体含水率的变化，该参数根据土壤性质，取经验值0.000 1 m-1；r(θ)为通过潜在蒸发量计算实际蒸发量的参数(无量纲)，与表层土壤含水率有关；ETp和Pr为单位时间内的潜在蒸发量和降雨量(m/d)，根据前文假设可表示为：

ETp=0.002 m/d， 0≤t≤tend

(7)

(8)

式中：t为计算时刻，tend为模型终止时刻，即土体底部节点刚好饱和时刻。本文对于式(3)～式(6)采用计算精度高、数值收敛性较好的Ross方法进行求解。

3 结果与分析

3.1 不同林地类型土壤水力特性差异

土壤水动力参数平均值如表2所示，各典型植物群落样方上下层土之间的参数值差异较小，土壤渗透性较好，饱和水力传导度Ks可达102cm/d；饱和含水量较大，均超过0.5。裸地和茶园上下层土之间的饱和水力传导度Ks差异很大，可达1个数量级。变异系数Cv表征了随机变量的离散程度，Cv大于0.36为高度变异，0.16～0.35为中等变异，小于0.16为小变异[10]。分析土壤水动力参数变异系数可发现，参数Ks属于中等变异至高度变异范围，参数α绝大部分属于高度变异，参数n、θr和θs为小变异范围。

表2 典型林地类型样地和参照样地土壤水动力参数统计结果

根据表2中的土壤水动力参数和式，可绘制不同典型林地类型样方和各参照样地的土壤水分特征曲线，如图1所示。可以发现各典型林地类型样方的土壤水分特征曲线靠右，茶园、水田和裸地的曲线依次靠右。说明在同样负压水平下，各典型林地类型样方的土壤比参照样地土壤的持水性更好。各典型林地类型样方内比较土壤持水性，杉木林略占优势，而毛竹林稍显弱。

图1 典型林地类型样地和参照样地土壤水分特征曲线

3.2 不同林地类型土壤水分调控效应模拟

通过模拟，可输出模型终止时刻tend(即降雨期末时刻，如表 3所示)典型林地类型样地土壤的底板开始产流时间(模型终止时刻)均比参照样地长，说明森林土壤本身(不包含植物的作用、土层厚度影响等)就具有延缓径流输出的特性。各不同典型林地类型样地中，土壤延缓径流输出的能力由大到小依次为毛竹林、柳杉林、杉木林、檫木林。

不同典型时刻(蒸发期末时刻、降雨期末时刻)的土壤含水量模拟结果如图2和图3所示。典型林地类型样地土壤的蒸发期末时刻含水量高于各参照样地，说明森林土壤土质(不包含植物的作用、土层厚度影响等)的抗旱能力更强。典型林地类型样地土壤中，抗旱能力由大到小依次为杉木林、檫木林、柳杉林和毛竹林。典型林地类型样地土壤的降雨期末时刻含水量高于各参照样地，说明森林土壤土质(不包含植物的作用、土层厚度影响等)的土壤水库蓄积能力更大。典型林地类型样地土壤中，土壤水库蓄积能力由大到小依次为杉木林、柳杉林、毛竹林和檫木林。

表3 典型林地类型样地和参照样地土壤水分模拟终止时刻

图2 蒸发期末时刻土壤含水量剖面分布

图3 降雨期末时刻土壤含水量剖面分布

在蒸发期，土壤耗水去向为大气，土壤耗水量越大说明土壤增加空气湿度的能力越强。蒸发期土壤耗水量可采用蒸发期起止时刻的剖面累积含水量之差计算，即：

ΔV(t)=V(ts)-V(te)

(9)

式中：te和ts分别表示蒸发期起、止时刻，V(t)表示t时刻的土体储水量；ΔV(t)表示土壤水分消耗。其中，土体储水量V(t)可由下式计算：

(10)

式中：N表示节点总数，本文取为20；i=1,2,...,N表示节点序号；θi(t)为t时刻节点i的土壤体积含水率(cm3/cm3)，为模型模拟结果；Δdi为节点i的控制长度(m)，本文模拟土壤剖面深度共计0.40 m，考虑各节点等距分布，故Δdi均为0.02 m。计算结果如图 4所示，森林土壤土质(不包含植物的作用、土层厚度影响等)的土壤增加空气湿度能力明显强于参照样地土壤。典型林地类型样地土壤的蒸发期土壤耗水量为裸地的4.2～7.3倍。各典型林地类型样地中，土壤增加空气湿度能力由大到小依次为杉木林、檫木林、柳杉林、毛竹林。

图4 蒸发期土壤耗水量计算结果

4 结论

本文对江西武夷山自然保护区典型林地类型样方(檫木林、柳杉林、毛竹林和杉木林)和参照样地(裸地、茶园、水田)的土壤水力性质和土壤水分进行了分析，主要取得了以下结果。

1)分析了典型林地类型样方和参照样地土壤水力性质，发现典型林地类型样方的土层较薄，上层土与下层土壤水力性质差异小，土壤透水性较好，饱和水力传导度Ks可达102cm/d，而参照样地裸地和茶园上下层土的饱和水力传导度Ks差异很大，可达1个数量级。总体上，各典型林地类型样方的土壤比参照样地土壤的持水性更好。

2)考虑典型林地类型样方中土壤水力性质的不同，以Richards方程模拟了土壤典型降雨-蒸发事件下的土壤水运动过程，从而说明土壤水分调控效应。各典型林地类型土壤比参照样地能更好地延缓径流输出、蓄积土壤水分以及增加空气湿度。各典型林地类型中，土壤延缓径流输出的能力由大到小依次为毛竹林、柳杉林、杉木林、檫木林，土壤水库蓄积能力由大到小依次为杉木林、柳杉林、毛竹林和檫木林，土壤增加空气湿度能力由大到小依次为杉木林、檫木林、柳杉林和毛竹林。