张法伟, 仪律北, 郭小伟, 杨永胜, 李杰霞, 曹广民, 李英年
(1. 中国科学院 三江源国家公园研究院,青海 西宁 810008; 2. 中国科学院 西北高原生物研究所 高原生物适应与进化重点实验室,青海 西宁 810008; 3. 青海省林业和草原局 林业碳汇服务中心,青海 西宁 810008)
土壤水是土壤-植被-大气连续体的主要组成,也是四水(地表水、地下水、大气水和土壤水)转换的纽带及土壤系统内物质和能量交换的载体[1],对区域水资源的形成、保持与转换具有重要意义[2]。祁连山国家公园地处青藏高原东北部,总面积5.02万km2,是河西走廊主要内陆河的发源地和黄河流域重要的水源产流地[3]。高寒草甸是祁连山国家公园青海片区重要的植被类型和区域水文调节功能发挥的载体[4-6],频繁发生的夏季融化、冬季冻结的冻融循环是其地表显著的物理特征之一[7-8]。冻融过程导致土壤水热特性的变化会直接影响到下垫面水源涵养功能的质量和状况[9-10],是高寒地区地表水循环的关键过程之一[11-12],对植被群落演替及生态系统稳定性维持具有重要的影响[13-14]。因此,深入研究冻融过程对土壤水分储量及交换的影响机制是量化高寒草甸水源涵养功能的基础和认知区域土壤水分植被承载力的关键[15-17],也是评估中下游的水资源安全的科学支撑[18-19],具有重要的理论价值和现实意义。
土壤水分储量及通量是土壤水势对水分运移的综合作用结果[2],受到降水、蒸散发和土壤渗漏及冻融过程的强烈影响[20-21]。基于涡度相关技术的研究表明高寒草甸蒸散发的季节变异主要受辐射有效能控制,和降水量基本持平[22-23]。深层(100 cm)土壤渗漏量约相当于年降水的5%[24]。高寒草甸植被生长季土壤水分储量呈现出双峰型季节变化,密切关联于土壤温度[25-26],与降水的关系存在时间尺度效应[12,21]。由于自动观测技术无法测定土壤中固态水含量,加之冻融过程的复杂性,冻融期土壤水分交换过程是目前高寒草地生态水文学研究的难点和热点之一[13,19]。由于冻土具有较高的导热性和较低的导水性[21,24],在温度梯度及水分相变所导致的土壤水势变化的驱动下,土壤水分向冻结峰或融化峰汇聚[27]。同时,冻融期中土壤水势的温度效应十分显著[28],土壤水热传输强烈耦合,进一步增加了高寒草甸土壤水分储量及通量对冻融过程响应研究的难度[19,27-29]。基于水热传输的模型模拟成为研究冻融期水热交换过程的重要手段[8,27,30]。相关研究表明,SHAW 模型综合考虑了冻融期土壤水热物理传输过程,对冻融阶段土壤水分的模拟结果更接近于实际观测[30],能够较好地模拟高寒草地土壤水热交换[8,26,29]。但早期研究多偏重于模型的评估和改进,少有应用模型探讨高寒草地土壤水分储量及通量的变化特征。因此,以祁连山东段冷龙岭南麓的高寒草甸为对象,在SHAW 模型参数本地化的基础上,研究冻融期0~100 cm 土壤水分储量及通量的季节变化特征和环境调控过程,以期为祁连山国家公园高寒草甸土壤水源涵养功能评估提供数据支撑和理论依据。
研究区位于青海海北高寒草地生态系统国家野外科学观测研究站(简称海北站,37°37′ N,101°19′ E,海拔3 200 m)。海北站地处青藏高原东北隅祁连山东支冷龙岭南麓,为典型的高原大陆性气候,年均气温和降水分别为-1.7 ℃和580 mm,雨热同季[31]。植被类型为禾草-矮嵩草草甸,优势种为矮生嵩草(Kobresia humilis)、垂穗披碱草(Elymus nutans)、异针茅(Stipa aliena)、麻花艽(Gentiana straminea)和黄花棘豆(Oxytropis ochrocephala)等。植被地上净初级生产力和0~40 cm 地下根系现存量分别约为400 g·m-2和3 200 g·m-2,其中约90% 的根系集中在0~20 cm 土层中,形成草毡层[31]。土壤类型为高山草甸土,发育相对年轻,粗骨性强。由于地表草毡层具有较高的导水性和较低的导热性,土壤水热特性具有显著的垂向分层特异性[12]。
研究区地势平缓开阔,系统水热状况由相关仪器自动监测(表1)。土壤水热数据观测频率为30 min(CR800,Campbell,美国)。同时,基于涡度相关技术对下垫面与大气间的水热交换进行连续观测。涡度相关系统距离土壤水热观测点约100 m,采样频率为10 Hz。利用EddyPro 7.06(Li-Cor,美国)对10 Hz 高频数据进行二次坐标选择、除趋势和WPL 密度校正,输出潜热通量、显热通量、风速、风向及空气温湿等30 min平均值。同时利用CM11观测1.5 m 高度处的太阳短波辐射。群落叶面积指数(LAI)来源于MODIS 发布的陆地植被产品(MOD15A2)。LAI以通量塔为中心,空间和时间分辨率分别为250 m×250 m 和8 d,来自于美国橡树岭国家实验室的分布式主动存档中心(Distributed Active Archive Center,Oak Ridge National Laboratory,http://daac. ornl. gov/MODIS/modis. html),用于提供SHAW 模型中植被生长的LAI 信息。土壤水热、潜热通量、显热通量及其他气象因子也集成为8 d的时间序列以进行后续分析。
表1 研究区系统水热状况的观测信息Table 1 The observation information of hydrothermal feature in the study site
SHAW 模型为Flerchinger 和Saxton 于1989 年开发的一维水热耦合模型,该模型基于表面边界层的能量平衡原理及在土壤-植被-大气连续体中的水热物理传输过程,采用隐式有限差分公式表述土壤的水热交换,然后通过Newton-Raphson 迭代求解[20],已成为冻土区陆面过程研究的常用手段之一[8,26,29]。由于观测仪器无法测定土壤中固态水含量,只有生长季土壤体积含水量的观测值才能代表土壤全部含水量,方可满足SHAW 模型初始和结束时段约束数据的模拟要求,因此本研究利用2017年8 月1 日至9 月30 日和2018 年6 月1 日至7 月31 日的土壤水热观测数据作为模型的初始和结束数据,以2017 年8 月1 日至2018 年7 月31 日的逐时气象数据[气温、相对湿度、风速、降水、新雪密度(设定为0,模型依据气温计算)和太阳辐射]作为强迫资料,进行SHAW 模型模拟研究。根据研究区特征,将坡度、坡向、地表粗糙度分别设置为0°、0°和1.0 cm。0~100 cm 土壤被分为11 层,分别为0、5、10、15、20、30、40、50、60、80 cm 和100 cm,土壤水热参数设置详见表2。其中砂粒、粉粒、黏粒、有机质含量和容重等指标来自于站点的长期监测数据[32],饱和导水率、空气进入势、饱和含水量和孔径指数首先根据文献资料设定初始值[8,29-30],然后根据SHAW 模型土壤水热的拟合值与观测值的接近程度进行参数优化。干土和植被反照率分别为0.35和0.25[22]。植被最小气孔阻抗设置为380 s·m-1[23]。由于是冬季牧场,没有考虑植被残留层,生长季的植被高度、地上生物量和LAI等参数来自观测数据,叶片宽度、丛生参数和有效根深则为经验值,分别设置为1.0、0.9 cm 和45 cm。其余部分参数选取模型提供的建议值。采用SHAW 3.0.2(https://www. ars. usda. gov/pacific-west-area/boise-id/northwest-watershed-research-center/docs/shaw-model/)进行模拟研究。土壤水分储量(SWS)和水分通量(SWF)分别表示土壤中液态水与固态水的水分总量和液态水与气态水的累计交换量[20],为参数优化后SHAW 模型的输出结果。其中土壤固态水含量是土壤水势和温度的函数,通过设定冻结温度依据能量守恒计算而来。SWS0-100为0~100 cm 土壤剖面SWS 与土壤深度的积分,SWF0-100为土壤剖面SWF的累计值,以表征水分的交换强度[20]。依据高寒草甸植物根系和土壤水热分层特性,将0~20 cm、20~60 cm和60~100 cm依次划分为浅层、中层和深层土壤[25],分别用0、5、10、15、20 cm 和30、40、50、60 cm及80、100 cm 的平均值表示浅层和中层及深层土壤温湿性状。浅层土壤水分储量(SWS0-20)、中层土壤水分储量(SWS20-60)和深层土壤水分储量(SWS60-100)分别是浅层、中层、深层土壤水分储量之和,浅层土壤水分通量(SWF0-20)、中层土壤水分通量(SWF20-60)及深层土壤水分通量(SWF60-100)分别浅层、中层、深层土壤的土壤水分通量之和。
表2 SHAW模型各层土壤节点参数Table 2 The soil parameters of the SHAW model at different soil nodes
利用统计学参数,如水热观测值与模拟值的决定系数(R2)、均方根误差(RMSE)和平均绝对误差(MAE)评价SHAW模型模拟效果[7,29]。在8 d时间尺度上,利用结构方程模型评估环境因子对0~100 cm土壤水分储量(SWS0-100)及水分通量(SWF0-100)的驱动过程及相对强度。结构方程模型是在对整体模型拟合和判断的基础上,阐释多因子之间相互关系的一种多元统计分析方法,在近年来的相关研究中应用广泛[23,33]。根据相关研究结果[25],SWS0-100的主要影响因子包括气温(Ta)、降水(Rain)、蒸散发(LE)、LAI、SWS0-20、SWS20-60及土壤温度(Ts20-60)、SWS60-100及土壤温度(Ts60-100);SWF0-100的主要影响因子包括Rain、Ta与浅层土壤的温度梯度[Δ(Ta-Ts0)]、浅层与中层土壤温度梯度[Δ(Ts0-Ts20)]、中层与深层土壤温度梯度[Δ(Ts20-Ts60)]、SWF0-20、SWF20-60及SWF60-100。基于结构方程模型中的标准作用系数和作用路径表征环境因子对SWS0-100和SWF0-100的相对强度和驱动过程。结构方程模型在R 4.0.3[34]平台上利用piecewiseSEM[35]软件包实现。缩写SWS0-20、SWS20-60、SWS60-100和SWS0-100分别是0~20、20~60、60~100 cm 和0~100 cm 土壤水分储量;LAI为群落叶面积指数;Rain 为降水量;Ta为气温;LE 为蒸散发;Ts20-60为20~60 cm 土壤温度;Ts60-100为60~100 cm 土壤温度;SWF0-20、SWF20-60、SWF60-100和SWF0-100分别是0~20、20~60、60~100 cm 和0~100 cm土壤水分通量;Δ(Ta-Ts0)为Ta与Ts0-20差值;Δ(Ts0-20-Ts20-60)为Ts0-20与Ts20-60差值;Δ(Ts20-Ts60)为Ts20-60与Ts60-100差值。
研究时段的日均Ta和累计Rain分别为-1.30 ℃和494.4 mm。日均太阳短波辐射(Swin)和相对空气湿度(RH)分别为194.24 W·m-2和63.6%,均表现出植被生长季(5—10 月)较高而非生长季(11月—翌年4 月)较低的季节动态(图1)。浅层、中层和深层的日均Ts分别为2.85、2.61 ℃和2.38 ℃,Ts表现出随土壤深度(x)增加而指数下降的特征(Ts=3.76×(1-x)-0.10,R2=0.63,P<0.01)。浅层、中层和深层的日均土壤液态水含量(SWC)分别为0.22、0.25和0.17 cm3·cm-3,除在15 cm[(0.15±0.08) cm3·cm-3]和100 cm [(0.09±0.03) cm3·cm-3]处较低外,其余各层稳定在0.22~0.27 cm3·cm-3之间。Ts和SWC的变异系数(CV)分别平均为40.1%和200.0%,均随着土壤深度增加而指数下降(R2>0.79,P<0.001)。值得说明的是,相邻深度的Ts和SWC的相关系数均大于0.95,表明土壤水热具有显著的层次相关性。
图1 2017年8月1日—2018年7月31日空气温度和降水(a)、相对湿度和太阳短波辐射(b)、土壤体积含水量(c)和土壤温度(d)的季节特征Fig. 1 Seasonal variations of air temperature (Ta) and precipitation (Rain) (a), relative humidity (RH) and solar radiation(Swin) (b), volumetric soil water content (SWC) (c), and soil temperature (Ts) (d) from August 1, 2017 to July 31, 2018
逐日SWC 和逐日Ts的模拟值和实测值的对比结果表明,模型对SWC 的模拟效果较好,R2在0.80以上,RMSE 和MAE 分别小于0.05 cm3·cm-3和0.03 cm3·cm-3(表3)。模型对Ts的模拟效果随土壤深度增加而降低,R2平均为0.71,RMSE 和MAE 分别平均为2.78 ℃和2.11 ℃。SHAW 模型模拟的日均显热通量和潜热通量分别为(29.8±17.3) W·m-2和(45.2±42.9) W·m-2,较涡度相关技术观测值分别高约4.2%和1.3%。模型对潜热通量的模拟效果好于显热通量,R2分别为0.81 和0.14(图2),RMSE分别为18.67 W·m-2和15.97 W·m-2。
图2 热量通量的涡度相关技术观测值与SHAW模型模拟值的关系:潜热通量(a),显热通量(b)Fig. 2 Relationships between the simulated heat flux with SHAW model and the observed heat flux with eddy covariance techniques: Latet heat flux (a), sensible heat flux (b)
表3 SHAW模型对各节点土壤温度和土壤水分的模拟效果的统计参数Table 3 Statistical parameters of the simulation performance of the estimated soil temperature and volumetric soil water content at each node with the SHAW model
日均SWS0-100为(274.99±19.57) mm,表现出生长季低(264.51 mm)而非生长季高(286.42 mm)的季节特征,最低值(224.20 mm)和最高值(298.67 mm)分别出现在6 月初和1 月初(图3),变异系数(CV)为7.1%。日均SWS0-20和SWS20-60分别为(59.39±9.07) mm 和(134.25±12.37) mm,也表现出生长季低而非生长季高的趋势。SWS60-100平均为(81.35±5.72) mm,但表现出生长季高(84.22 mm)而非生长季低(78.22 mm)的特征。SWS0-20、SWS20-60和SWS60-100的CV 分别为15.3%、9.2%和7.0%,随土壤深度增加而降低。SWS0-20和SWS20-60的CV 在生长季和非生长季差别较小,分别平均为10.1%和7.4%,但生长季SWS60-100的CV 为7.3%,远高于非生长季的3.8%。
图3 高寒草甸不同层次土壤水分储量的季节变化Fig. 3 Seasonal variations of soil water storage of different layers in the alpine meadow
SWS0-100的结构方程模型表明,其季节变异主要受SWS20-60和SWS0-20调控,其标准作用系数分别为0.63 和0.46[图4(a)],二者可解释SWS0-100季节变异的91.4%。SWS0-20与LAI 正相关,与Ta和LE 负相关,但LAI、Ta和LE 的作用强度基本相当。SWS20-60主要受SWS0-20、Ta和Ts20-60显著影响,其中Ta的正效应(0.80)大于Ts20-60的负效应(-0.66)。Ts60-100和SWS20-60为SWS60-100的主要影响因素,其中Ts60-100正效应的标准作用系数为0.93,远大于SWS20-60的正效应。整体上看,LE 和LAI 对SWS0-100的综合效应最大,分别为-0.67 和0.59,LE 主要通过对SWS0-20(-0.29)和SWS20-60(-0.38)的负效应而LAI 则通过对SWS0-20(0.25)和SWS20-60的正效应(0.34)影响SWS0-100。Ta对SWS0-100的整体效应为-0.12。值得注意的是,在8 d 时间尺度上,研究时段内降水对SWS的季节变异无显著影响。
图4 高寒草甸土壤水分储量(a)和土壤水分通量(b)的结构方程模型Fig. 4 Piecewise SEM models for soil water storage (SWS) (a) and soil water flux (SWF) (b) in the alpine meadow
日均SWF0-20、SWF20-60和SWF60-100分别为(0.92±4.88)、(-0.34±4.09) mm·d-1和(-0.19±1.81) mm·d-1,都表现出生长季水分通量为正而非生长季为负的季节特征,即土壤水分在生长季向下运移,非生长季向上汇聚(图5)。生长季中,SWF0-20(2.85±5.88 mm·d-1)远高于SWF20-60(0.48±5.08 mm·d-1)和SWF60-100(0.30±4.34 mm·d-1)。非生长季中,SWF0-20、SWF20-60和SWF60-100差异较小,分别为(-1.18±2.10)、(-1.24±2.44) mm·d-1和(-0.73±2.24) mm·d-1,三者平均为(-1.05±2.26) mm·d-1。冻结期各层土壤水分通量基本为0 mm·d-1,但不同层次存在时间滞后,滞后速率约为0.56 d·cm-1。日均SWF0-100为(0.16±9.52) mm·d-1,在生长季和非生长季中分别为(3.27±11.74) mm·d-1和(-3.23±4.50) mm·d-1。
图5 高寒草甸不同层次土壤水分通量的季节变化Fig. 5 Seasonal variations of soil water flux of different layers in the alpine meadow
SWF0-100的结构方程模型表明,SWF0-20、SWF2060-和SWF60-100对SWF0-100的标准作用系数分别为0.43、0.36 和0.32,三者作用强度基本相当[图4(b)]。SWF0-20对SWF20-60、SWF20-60对SWF60-100具有显著的正效应。SWF0-20的季节变异主要受Rain驱动,其标准作用系数高达0.95。SWF20-60主要受浅层和深层温度梯度影响,二者标准作用系数分别为0.56 和-0.46。SWF60-100主要受中层、深层土壤温度梯度影响,二者标准作用系数分别为0.62 和-0.40。降水对SWF0-100的总体作用系数为0.99,主要通过对SWF0-20(0.40)和SWF20-60(0.35)的间接效应来影响SWF0-100。浅层、中层和深层温度梯度的标准作用系数分别为0.38、0.20和-0.44。即SWF0-100季节变异的主要影响因子是降水,其次是深层土壤温度梯度。
本研究SHAW 模型对高寒草甸土壤水热变化的模拟效果相对较好,与该模型在其他相关研究中表现(R2>0.70)较为一致[7,26,29],表明SHAW 模型适合于高寒草甸0~100 cm 土壤水热传输的模拟,模拟结果能够用以评估土壤水分状况。然而,本研究SHAW 模型对SWC 的模拟效果略好于Ts的模拟效果,这与郭林茂等[8]研究结果不同。一方面,土壤水分参数的不确定性是陆面过程模型的重要限制之一[26,29],本研究的一部分土壤水分参数(如容重、机械组成和有机质)是基于地面实测数据而来,另一部分难以直接测定的参数(如孔径指数和空气进入势等)是通过降低模型模拟结果和观测数据之间的残差优化获取的,有效地提高了SHAW 模型对SWC 的模拟效果[7]。另一方面,Ts模拟效果略差主要是由于地表(0 cm)Ts是按照气象观测规范所测定的裸露地表Ts,并非有植被覆盖的真实地表Ts,尤其在植被生长季中二者的差异较大[36]。而SHAW 模型是按照土壤节点顺序进行土壤水热模拟,在保障表层Ts模拟效果的基础上,对深层Ts进行线性插补[20],导致深层Ts的模拟值与观测值差别较大,即Ts的模拟效果随着土壤深度增加而下降,这和唐古拉地区SHAW 模型对土壤剖面Ts的模拟效果相似[26]。Ts的模拟误差也影响了地气温度梯度及显热通量的模拟效果,这与西藏那曲[29]和唐古拉地区[7]的相关研究结论一致。
由于研究区100 cm 下土壤渗漏量很少[24],SWS0-100主要取决于Rain 和LE 的平衡。SWS0-100在生长季低而非生长季较高,表明土壤冻结有利于土壤系统水分的保持[10,19],这由于冻土具有较高的隔水特性,限制了非生长季LE 的过度损失[8],加之冻结过程对下层土壤水分具有汇聚作用[12]。因为研究区在6 月初的太阳辐射较强而LAI 较小[22],系统获取的辐射有效能较大,加之冻土消融提高了表层土壤水分含量,刺激了LE 的迅速增加[23],导致SWS0-100出现低值。
SWS0-100的季节变异主要受SWS0-20和SWS20-60的影响,这主要与浅层和中层水分变化较为活跃,深层水分含量相对稳定有关[25]。SWS0-20受LE、Ta和LAI 的综合影响[图4(a)]。LE 对SWS0-20为负效应,表明LE是SWS0-20散失的主要途径[22-23]。Ta的负效应主要由于Ta可以表征大气蒸发需求,且与LE在季节变化上也具有显著的正相关[22]。LAI 越大,植被蒸腾在LE 中占比越高,系统可通过调节植被气孔行为控制水分散失[22],加之高寒植被反照率较湿润裸地高,可以降低系统的辐射有效能及水分的过度消耗[23],有利于浅层土壤水分保持。SWS20-60主要受SWS0-20的正效应影响,因为二者水分储量的季节变化具有一致性(R2=0.50,P<0.01)。由于SWS60-100的季节变异主要受控于生长季SWS60-100的变异,而后者与生长季Ts60-100显著正相关(R2=0.73,P<0.001),进而导致SWS60-100的季节变异与Ts60-100正相关,这也与相邻区域高寒灌丛草甸的研究结果一致[21]。降水对SWS0-100在8 d 时间尺度上的季节变异无显著作用,一方面由于年降水的85%集中在植被生长季,且以小降水(<3.0 mm·d-1)事件为主,加之高原太阳辐射较强,极易被LE 迅速消耗[21-22]。另一方面,强降水事件[(31.4 mm)·(48 h)-1]最多只能作用到高寒草甸40 cm 的SWC[21,24],也一定程度上削弱了降水对SWS0-100的影响强度。
SWF 的实质是各种势能综合作用下的水分质量的转移,主要包括温度梯度驱动的未冻水迁移、重力作用下的自由水下渗和土壤毛细作用下的毛细水迁移[19]。生长季中,SWF 主要驱动力是由降水、蒸散发等过程引起的土壤含水量增减而导致的土壤剖面基质势梯度的变化,土壤水迁移以自由水和毛细水为主[27],降水对SWF的影响也随着土壤深度增加而下降[8],本研究也发现SWF20-60和SWF60-100的变化主要来源于降水的间接作用[图4(b)]。这主要由于高寒草甸植物根系集中在浅层土壤中形成了一系列大孔隙网络,提高了土壤的导水率[12]。非生长季中,土壤与大气间的水气交换被冻结的表层土壤所阻隔,土壤水分迁移与冻土融冻导致的温度变化密切相关[10,27],即温度梯度是SWF 的主要影响因素[图4(b)]。一方面,冻结过程中,土壤水分在温度梯度的作用下向冻结峰汇集,土壤水分从深层向浅层聚集[27],从而导致了Δ(Ts20-60-Ts60-100)对SWF20-60和SWF60-100的负效应。消融过程中,由于季节冻土的双向消融,加之浅层融化速率较快[10,12],在浅层土壤中形成一个水分的高值区,在温度梯度、重力势和基质势的综合作用下,土壤水分向下运移[10],导致Δ(Ta-Ts0-20)和Δ(Ts0-20-Ts20-60)对SWF20-60和SWF60-100产生了正效应[图4(b)]。另一方面,温度是非生长季中土壤水分相变、未冻水含量及土水势的重要影响因素,还可通过改变土壤的导水率和水黏滞系数[28]来影响土壤水分通量。
利用祁连山南麓高寒草甸2017 年8 月1 日—2018 年7 月31 观测资料作为SHAW 模型的强迫资料,进行了冻融过程中土壤水热交换的模拟研究,并应用模型结果分析了季节冻结层SWS 及SWF 的变化特征及环境驱动,得到如下主要结论:
(1)日均SWS0-100为274.9 mm,约50%集中在20~60 cm 中层土壤,生长季和非生长季的日均SWS0-100分别为264.5 和286.4 mm。SWS0-100最低值和最高值分别出现在6 月初和1 月初,季节变异较小。蒸散发和群落叶面积指数通过影响SWS20-60和SWS0-20驱动SWS0-100的季节变化。
(2)日均SWF0-100为0.16 mm·d-1,生长季和非生长季的分别为3.27 和-3.23 mm·d-1,表现为土壤水分在生长季向下传输而非生长季向上汇聚。生长季浅层水分交换活跃,SWF0-20为2.85 mm·d-1,而非生长季中浅层、中层和深层差异较小,SWF 平均为-1.05 mm·d-1。降水和深层土壤温度梯度通过影响SWF0-20、SWF20-60和SWF60-100驱动SWF0-100的季节变化。