土壤的碳、氮、磷元素作为土壤养分的重要组成部分, 在森林生态系统中具有重要的生态价值, 对改良土壤性质, 促进植被发育、土壤酶活性、土壤微生物活性与生态系统能量循环等方面发挥着显著的作用^[1-5]。土壤碳、氮、磷元素在生态系统中参与动态的循环过程, 不同循环过程中的元素会对生态环境存在不同的影响^[6]。作为衡量土壤质量的重要指标, 土壤碳、氮、磷元素与土地生产力联系紧密, 人为地增加碳氮磷元素能够在一定程度上促进生态系统的恢复及植株的生长^[7]。但土壤中的碳、氮、磷等元素容易随着击溅、冲刷等土壤侵蚀过程产生运移, 造成水体富营养化、土壤养分流失等环境问题；另一方面, 土壤的碳储量的减少会促进有机碳的排放, 加剧温室效应^[8]。

冀北山地地处华北滦河上游, 为国家划定的重要的生态防线之一, 担负护卫京津冀生态安全、阻挡风沙的重任^[9]。20世纪60年代, 当地营造了以华北落叶松为优势树种的人工林, 90年代以来, 出于保护天然林资源的政策的原因, 天然林林分蓄积的开发受到严格限制, 当地对已营造的人工混交林采取了以封为主, 人工辅育为辅的措施, 以加快成林速度并提高林分质量^[10]。高敏等^[11]对该地区落叶松人工林的土壤动物群落多样性进行研究发现, 经人工辅育的落叶松混交林通过对地上植被组成与结构的调控, 能够改变土壤微生物种类、改善土壤结构及理化性质, 提高了林分环境的异质性, 促进了土壤动物群落的多样性的提高。但迄今为止, 关于次生林立地因子与植被对不同土壤深度层次土壤养分性状之间的关系的研究尚不多见。

本文意在从林分、地形等因子的角度出发, 探究不同土壤养分指标之间的计量关系与影响土壤养分含量及空间分布的因素。为揭示次生林土壤养分变化以及合理抚育山区次生林, 提高次生林生产力提供必要的依据。

1 研究区域概况

研究区域位于河北省承德市木兰围场自然保护区内北沟林场, 该区域位于坝上高原与坝下山地的交界地带, 北接御道口牧场与河北省塞罕坝机械林场, 东临滦河最大支流小滦河, 地理坐标41°47′—42°06′ N, 116°51′—117°45′ E。海拔1010—1939.9 m, 整体地势呈现西北高, 东南低的局面, 属暖温带亚湿润针叶林区向中温带半干旱森林草原区的过渡带, 为大陆性季风性山地气候, 区域年降水量380—560 mm, 集中于6、7、8月, 全年降水量曲线呈“单峰”型变化。多年平均气温-1.4—4.7℃, 日照充足, 雨热同季, 冬长夏短, 四季分明, 昼夜温差较大。

木兰围场自然保护区区域内植被类型丰富, 呈现明显的地带性分布。分布有原始森林生态群落退化后形成的阔叶灌丛与草甸, 其中草甸以地榆-细叶苔草草甸与披碱草-委陵菜草甸为主要类型的杂草类草甸。研究样地为经火灾干扰后次生演替, 于火烧迹地上重新生长起的先锋群落。20世纪四五十年代, 研究区域曾经历一次对部分大径乔木的择伐, 后被封育进入自然恢复状态。目前已演替为以山杨(Populus davidiana Dode)、白桦(Betula platyphylla Suk)及落叶松(Larix principis-rupprechtii Mayr.)为主要优势树种的天然次生林。研究样地林分为区域地带性植被在自然恢复、演替过程中的过渡阶段。除上述优势树种外, 研究样地还分布有数量相当的杜松(Juniperus rigida Sieb. et Zucc.)、油松(Pinus tabulaeformis Carr.)、红松(Pinus koraiensis Siebold et Zuccarini.)、樟子松(Pinus sylvestris Linn. var. mongolica Litv.)、旱柳(Salix matsudana Koidz.)及栎类(Quercus)等常绿及落叶乔木。

2 研究方法 2.1 样方设置及调查

试验地选择林场内面积为200 m×200 m的固定样地一处, 海拔高度1330—1440 m, 使用全站仪将样地划分为100个面积为20 m×20 m的样方。于2017年5月份对所有样方内的植被、土壤及地形等因子进行全面细致的调查。研究区域设置的100组样方中, 设定编号数尾数为0的样方位于图 1最上方。由图 2可见, 样地海拔高度整体呈西北-东南方向递减, 各组样地间坡向不一。位于最上方的10组样地中, 中部的40、50与60号样地的海拔高度相对左右两侧相对较低。

2.2 植被调查

基于20 m×20 m面积的样方, 对样方内所有胸径≥1 cm的木本植物进行调查, 调查内容包括种名、树高、胸径、冠幅及定位坐标；在每个20 m×20 m面积样方内的四角及中心点设置1 m×1 m面积的小样方用于草本植物的调查, 内容包括种名、盖度、高度、生长状况及分布情况等。对高度不足1.5 m及胸径不足1 cm的灌木, 统计平均胸径、平均树高及株数。

2.3 样方内叶面积指数测定

将每块样地划分为10 m×10 m的小样方4块, 于每块小样方的中心点作为叶面积指数(leaf area index, LAI)取样点。为避免太阳光照对镜头摄像效果的影响, 试验选择于光照较弱的清晨或日落时分进行观测。观测时, 使用CI-110叶面积指数测定仪于1.5 m高度处通过鱼眼镜头向上拍摄, 测定时仪器须保持水平。由仪器自带软件根据拍摄结果进行自动计算结果值。计算4块小样方所测LAI值的平均值作为样地的LAI值。

2.4 土壤因子测定 2.4.1 土壤水分测定

以10 m×10 m的小样方中心点为测点使用便携式时域反射仪TDR测定0—20 cm、20—40 cm与40—60 cm三个土壤深度层次的含水量, 测定时间选择于测定前连续3日无降雨的阴天进行。

2.4.2 土壤养分测定

因实验需求, 测定土壤养分指标需获取较大量的土壤样品, 实验基于20 m×20 m的样方面积, 于样方中心点1 m范围内设置3—5个采样点, 使用土壤采样器分层采样后混合, 得到样方中心范围内土壤的复合样品。因实验设计样方为方正的矩阵式排列, 20 m×20 m面积的样方于中心布设1处土壤采样点, 所得结果能够体现同一指标在不同样地间的变异。土壤养分测定指标包括：pH值、土壤有机碳(Soil organic carbon, SOC)、全氮(Total nitrogen, TN)、全磷(Total phosphorus, TP)、全钾(Total potassium, TK)、速效磷(Available phosphorus, AP)、速效钾(Available potassium, AK)。其中, pH值使用电位法测定, SOC、TN、TP含量需对土壤研磨并过直径0.25 mm筛后, 分别利用重铬酸钾外加热氧化法、半微量凯氏定氮法、NaOH熔融-钼锑抗比色法测定^[12], AP采用NaHCO₃浸提-钼锑抗比色法测定^[13], AK、TK采用浓度为2 mol/L的HNO₃溶液振荡过滤测定^[14]。

2.5 不同分解程度的枯落物观测及采集

以20 m×20 m面积样方的中心点及对角线1/4、3/4距离处布设测点, 每处测点使用直尺测定枯落物未分解层、半分解层、分解层各自的厚度。在样方的中心点处设50 cm×50 cm的小样方分层采集并带回室内作烘干处理, 推算出各层枯落物的量。

2.6 数据处理

根据试验测定结果, 计算每组样地土壤C:N、C:P与N:P值。绘制据计算结果所绘制各养分指标含量折线图。

使用SPSS 22.0软件对所有试验所观测、获取的数据进行正态分布检验, 符合正态分布的数据可直接用于下一步的分析；不符合正态分布的数据, 进行对数转换处理后用于Pearson相关性检验。

研究选定叶面积指数(Leaf area index)、乔木胸高断面积(Arbor basal area)、乔木密度(Arbor density)、冠幅(Canopy)、草本盖度(Herb coverage)、灌木盖度(Shrub coverage)、灌木胸高断面积(Shrub basal area)、灌木密度(Shrub density)8组林分因子, 土壤含水率(Soil moisture)、未分解层(Undecomposed)、半分解层(Semi-decomposed)、分解层(Decomposed)各层枯落物厚度4组土壤因子与海拔高度(Altitude)、曲率(Convex)、坡向(Slope direction)、坡位(Slope position)、坡度(Slope degree)5组地形因子共17组因子作为解释变量。坡位、坡向3种因子, 作为参与冗余分析(Redundancy analysis, RDA)的解释变量的一部分, 为更好地解释土壤养分的变异, 对坡向、坡位因子进行定量化赋值。对坡向而言, 将正北、东北、西北、正东、东南、正西、西南及正南共8个坡向分别赋值1、2、3、4、5、6、7、8^[15-19]；计算样地距坡顶距离占整个坡面长度的百分比作为定量化的坡位值。

使用RDA分析解释变量与土壤养分关系之前, 需要对解释变量进行去趋势对应分析(De-trend correspondence analysis, DCA)。结果显示, 排序轴的梯度长分别为1.096、0.496、0.600、0.354, 低于RDA分析要求的梯度长不高于3的上限, 因此选用RDA进行排序分析优于选用单峰的典范对应分析(Canonical correspondence analysis, CCA)。为分析各解释变量之间的相互作用及因子自身对响应变量的解释能力, 在进行RDA分析后, 基于其结果进行变差分解(Variance part analysis, VPA)处理。数据处理、分析及绘图全过程于R 3.6.1中完成, 其余绘图于Origin 2018中完成。

3 结果及分析 3.1 天然次生林植被群落组成

天然次生林样地内共调查乔木计16种, 3126株, 乔木个体数较多的有山杨(Populus davidiana)、白桦(Betula platyphylla)、华北落叶松(Larix principis-rupprechtii)、华北五角枫(Acer mono Maxim)、蒙古栎(Quercus mongolica)等；灌木计18种, 其中重要值大于1的共计9种, 按从重要值从大到小的排序为毛榛(Corylus mandshurica)、锦带花(Weigela florida)、迎红杜鹃(Rhododendron mucronulatum)、沙梾(Swida bretchneideri)、土庄绣线菊(Spiraea pubescens)、金花忍冬(Lonicera chrysantha)、照山白(Rhododendron micranthum)、红瑞木(Cornus alba)、丁香(Syzygium aromaticum)。其中毛榛为灌木层中的绝对优势种, 重要值达到66.63；位于第二位的锦带花重要值仅9.52；草本计71种, 其中重要值≥1的草本21种, 从大到小排序依次为, 细针苔草(Carex rigescens)、小红菊(Dendranthema chanetii)、风毛菊(Saussurea japonica)、糙苏(Phlomis umbrosa)、舞鹊草(Maianthemum bifolium)、华北耧斗菜(Aquilegia yabeana)。细针苔草为草本层中的绝对优势种(表 1)。

研究统计所有样地间优势乔木山杨(Populus davidiana)、白桦(Betula platyphylla)与优势灌木毛榛(Corylus mandshurica)在所有100组样地中的分布及差异, 统计发现, 山杨与白桦主要分布于1—50号, 研究区域分布的3126株乔木中, 1—50号样地乔木数量为2448株, 占比78.31%。其中位于图 1左下方的25组样地分布数量为1612株, 占总数51.57%。分布于1—50号样地的山杨数量1746株, 占山杨总数81.82%；相同样地白桦数量为1028株, 占白桦总数58.41%。毛榛数量分布与山杨、白桦相似, 1—50号样地分布毛榛9229株, 占总数77.65%。其中, 13、14、15、29与46号样地分布毛榛数量均大于300株。

3.2 土壤养分含量分布与指标之间的相关性

在样地内, 土壤pH值的变化波动较小, 整个研究区域内的土壤pH值最高时不超过6.8, 研究区域土壤属于酸性土壤。其余各项指标在不同的土壤深度层次、不同的样地之间的分布有显著的差异(图 3)。土壤有机碳在不同深度土层中分布差异显著, 表层土壤SOC值在55—60号、66—70号与76—80号样地含量最高, 下层土壤SOC值较高的样地主要集中于位于图 1左下方25组样地内, 其中40 cm以下深度土壤SOC于1、2号样地分布最集中。所有样地内全氮含量分布基本均匀, 13号与14号样地为TN含量最高的样地。45—50号样地与55—60号样地为全磷含量较高的样地, 相同样地全磷含量随土壤深度增加呈逐步减小趋势, 有效磷由全磷中的部分成分转化而来, 但分布规律与全磷不同, AP含量较高的样地位于TP高分布样地右侧的65—70号样地, 除此之外, 随着土壤深度加深, 位于65号样地下方的63号样地及61号样地相继出现较高的AP含量, 含量值与70号样地相当。相同土壤深度土壤速效钾分布规律与土壤有机碳分布相近。

不同土壤深度各指标间的相关关系的显著性存在差异, 但部分指标与其他指标的相关关系的显著性随土壤深度的变化始终存在。TN与SOC之间的相关性由浅层及深层土壤保持极显著相关(P < 0.001), 说明样地间的C:N比波动平缓, 能够稳定于一定范围内。与之相似, SOC与AK值之间存在显著的相关关系(P < 0.01), 根据图 3中不同样地SOC与AK值的分布可见, 两者呈现极显著的正相关。出现这一现象的原因可能为由部分土壤全钾转化形成的AK能够按照植被根系所需土壤有机碳的含量, 及时形成并满足各组样地中植被根系的需要。由SOC及TN值所计算得出的各样地C:N值与AK保持显著相关同样能够说明这一点(表 2)。

3.3 林分、海拔及其他土壤因子与土壤养分指标之间的关系

在17组解释变量的条件下, RDA分析能够分别解释0—20 cm、20—40 cm与40—60 cm三个深度层次土壤养分指标的48.64%、53.45%及49.23%的差异性(P < 0.001, 表 3, 表 4)。其中以SOC、AK两项土壤养分指标被解释变量解释比例最大(表 3)。在0—20 cm深度土壤中, 17组参与环境因子解释的变量有12组对SOC产生解释, 其中, 灌木密度、冠幅、坡向、地表凹凸度、灌木盖度、叶面积指数、枯落物层变量对SOC解释明显。而解释AK含量分布的指标主要有坡度与灌木盖度、乔木密度等与植被生长有关的因子。20—40 cm深度土壤中, 解释SOC的因子发生了明显改变, 影响SOC指标的因素数量不及表层土壤, 土壤水分、乔灌胸高断面积因素对该层的土壤养分值有一定的影响。40—60 cm土壤中, 解释SOC含量及分布的因子主要为与灌木有关的因子。不同土壤深度土层解释AK的因子中, 部分因子始终存在, 灌木盖度HC, 坡度SLD与乔木密度AD均对所有土壤层次的AK作出解释(图 4)。

VPA结果显示, 基于前述的RDA分析结果, 在没有海拔高度、坡度、凹凸度、坡位、坡向等地形因子介入时, 林分因子与土壤因子的交互作用不能对土壤养分的差异性进行解释。由海拔高度、凹凸度、坡度、坡位、坡向组成的地形因子组, 其单独对各指标的解释率较低。除对40—60 cm较深层土壤养分尚有7.1%的解释比率外, 单独对40 cm以上的浅层土壤的解释能力很微弱, 但与林分、土壤因子相交互后, 解释能力明显增强。纵观3个土壤深度层次, 地形因子与其他因子相交互后, 解释比率占各自整体的比重分别达到40.69%、60.80%、71.25%, 分别为其单独解释比例的6.45、8.06、3.21倍(图 5)。

4 讨论与结论 4.1 讨论 4.1.1 不同土壤深度层次土壤养分含量分布差异及关系

研究区域所属的天然次生林封禁较早, 经多年次生演替与群落迭代, 林分生长已进入近熟林至成熟林阶段。其林下土壤养分含量与分布格局相比于研究区域附近林龄较小的次生纯林存在较大差异, 林场内其他纯林林分浅层土壤有机碳含量平均值34.41 g/kg, 最高值不超过58.85 g/kg^[20], 而研究区域各样地浅层土壤有机碳含量大多在80 g/kg以上, 较高的14、57与84号样地, 得益于较大量的枯落物存在, SOC含量超过120 g/kg。研究发现, 区域内1—50号样地因乔灌木分布密集, 林下产生较厚的枯落物层, 经微生物分解后形成SOC含量较高的样地, SOC的增加同时提高了AK的含量。这与Mercik等^[21]研究得出的向农田施用有机肥以提高其有机质含量进而提高土壤钾的有效性的结论一致。

TP作为衡量土壤各种形态磷元素的一个综合指标, 在一定程度上衡量着土壤的有效肥力的高低^{[20, 22]}。20 cm以下TP与其他土壤性质指标之间的相关性却很低。造成这一现象的主要原因可能为, 作为组成TP的主要成分有效P、有机P, 其主要来源于土壤微生物对枯落物及其他主要植物残体的分解^[23], 有研究表明土壤可溶性碳的含量可以影响土壤P的含量、组分及有效性^[24], 随着土壤深度的增加, 凋落物分解形成的养分难以大量地进入较深的土壤, 多年来该次生林分长期处于封育状态, 不同于同一地区的人工林进行定期择伐后会保留相当数量的采伐剩余物, 以增加造林地有机碳的含量为土壤微生物提供充足的养分与能量, 处于封育状态的次生林, 表层土壤有机碳经过淋溶后到达深层土壤的数量更少, 使深层土壤微生物得到的养分来源有限, 会使土壤以NaOH-Po形式保留土壤TP的能力降低, 加之缺少采伐物的覆盖, 因淋溶造成的土壤侵蚀也会导致一定的土壤P的流失^[25], 使同一地点采集的20 cm以下土壤TP、AP的含量明显低于表层土壤。AP含量较高的60—70号样地方位特殊, 65—70号样地位于左右两侧样地的下方位处, 61—64号样地位于65—70号样地的下方位处, 土壤P随降雨的淋溶入渗至较深土壤中形成壤中流, 随坡向向下坡位样地集中, 研究发现61、63号样地20 cm以下深度土壤AP含量明显高于浅层土壤, 其较高的AP可能来自于原分布于上坡位表层土壤中的AP。但调查发现, 尽管AP为重要的土壤养分来源, 富集AP的样地植被密度在研究区域内相对偏低, 远不及AP含量相对较低的1—50号样地。

涉及土壤碳、氮、磷养分的固持与矿化的过程, 土壤的C:N、C:P与N:P值反映土壤养分的转化^[26]。分别表征土壤有机碳的矿化程度, 土壤有机碳质对土壤中P的释放与吸收能力以及土壤氮素含量是否受限。经实验测定并计算, 研究区域多数样地C:N比接近或超过20, 全国平均水平为10.9^[27]；多数样地土壤C:P平均值自浅层至深层分别为92.14、74.61与69.25, 高于全国平均水平60.0^[28], 其中部分样地浅层与中层土壤C:P接近或超过200, 有研究表明^[29-30], C:P比超过200的土壤足以使P发生净固持, 降低土壤速效磷的转化效率。与C:N及C:P不同, 所有样地N:P平均值自浅层及深层仅6.44、4.58与4.75, 与全国平均水平相当, 但远小于邻近区域天然纯林^[20](图 3)。根据曹娟等人^[31]的研究结论, 林下土壤N:P值变化与林龄有很大关系, 在幼龄林阶段, N:P会随林龄的增加而明显增加, 在中龄林阶段达到最大值后, 又逐渐减小, 与研究区域所得出的结果一致。对各深度层土壤N:P值与土壤TN、TP值进行相关性分析发现, TN值与N:P的相关性更加显著, 说明影响N:P值的主要因素为氮素量的限制。

4.1.2 林分因子与地形因子对土壤养分含量的影响

研究选择的8组林分因子对SOC、AK两组土壤养分变量产生解释, 对其他养分变量解释时, 即使增加与其他因子的交替仍然相对微弱(表 3、图 4、图 5)。表明影响土壤养分的林分因子涉及植被种类, 形态等其他因子。研究发现部分土壤指标, 如TN的含量明显地偏低于同时期所测的人工林分。有研究表明, 主要优势乔木树种山杨、白桦的根系在土壤低氮的状态下, 山杨通过降低吸收根中皮层细胞的直径, 降低皮层厚度；白桦通过增加皮层细胞直径、减少皮层细胞层数降低皮层厚度。皮层厚度降低后将促使吸收根吸取更多的土壤N, 进而形成恶性循环^[32]。因研究方法存在的局限性, 用于RDA分析的数据矩阵仅能以样地为单位输入参数, 未能考虑树种的因素, 可能是造成林分因子对作为相应变量的土壤养分解释率相对较低的原因。

表层土壤由于接收了大量地表枯落物分解后形成的各种有机碳及其他土壤养分, 致使影响其含量、分布的主要因子始终为枯落物分解量, 弱化了其他因子的影响。这与淑敏等^[33]、王凯等^[34]所研究的结果完全一致。但20 cm以下的土壤接收的地表枯落物分解后的养分量远不及表层土壤, 地形等因素对土壤养分的影响随之逐渐凸显。海拔、坡度、坡向、坡位及曲率等地形因子显著影响20 cm深度以下的各项土壤养分指标(图 4)。其中, SOC、TK与各项地形因子呈较显著的负相关关系, 其他土壤养分指标与地形因子呈显著的正相关关系。随着土壤深度的增加, 地形因子对土壤养分含量、分布的解释能力逐渐增加(图 5)。魏亚伟等^[35]、苏日古嘎等^[36]、宋爱云等^[37]、沈泽昊等^[38]均认为地形因子对物种分布、林分结构、土壤养分的影响作用较为显著。地形可以通过对光照强度、温度、降水、土壤侵蚀等环境因子的再分配间接地影响土壤养分分布及植被生长状况^[39-40]。据试验设计时对样方号的标定, 样方号趋于整十数的样方趋于低坡位位置, 从上坡位到下坡位每10块样地的0—20 cm深度土壤SOC、TN含量均有较为明显的增加趋势(图 3), 即下坡位的SOC、TN含量明显高于上、中坡位。与卫玮等^[19]的研究结果一致。主要原因在于, 上、中坡位的土壤在一定的降雨击溅及重力因素的侵蚀作用下, 部分土壤随着时间的推移逐渐地被搬运并堆积至下坡位处, SOC、TN等土壤养分也随着土壤颗粒向下坡位运移, 使下坡位的土层较上、中坡位偏厚, 养分容易富集。除此之外, 海拔因子也是影响山区植物群落、土壤养分分布的主导因子之一, 它综合反映了林地降水、温度、相对湿度等诸多的气象因素, 在各种大小的研究尺度上都对物种的分布、多样性有重要影响^[41]。刘志理等^[42]在总海拔高差仅有89 m的小兴安岭红松林地样地内, 仍发现海拔高度作为研究所选定的地形因子之一, 能够显著地影响叶面积指数的空间分布特征, 主要表现为叶面积指数值随海拔的变化呈显著的增加趋势。本研究区域海拔高差最大仅有102 m, 仍然表现了与20 cm以下深度土壤各项养分的显著相关甚至极显著相关。但较大空间尺度的海拔高差是否会对土壤养分、植被群落组成带来更显著的差异还需进一步地研究。

4.2 结论

本文探究基于次生演替并封育多年的天然次生林下土壤生态化学计量特征及差异, 结果表明：封育多年的天然次生林下土壤中C:N与C:P大于同区域内其他林分, 也大于全国平均水平, N:P值相较于同区域内其他经营较晚并处于幼龄林与中龄林的林分较低, 与全国平均水平相当, 次生林生长主要受N控制。封育多年的次生林可以进行适当地采伐剩余物, 施加N素等抚育措施以保持天然次生林内土壤C、N、P含量更加均衡。多数样地土壤不同深度C、N、P含量维持在较高水平, 植物对土壤养分保持着较强的吸收能力, 土壤磷、钾的固持效果明显, 能够转化为足量的速效磷与速效钾供给植物生长。研究范围内土壤生态化学计量特征不仅与植被类型、植被密度紧密相关, 同时与林地地形联系密切。在未来人工林经营与天然植被恢复过程中, 应当注重林地实际状况, 因地制宜, 选择合适的播种位置, 促进植物吸收土壤养分的同时提高养分在生态系统中的循环能力。