1 引言
青藏高原是全球气候变化的敏感区， 其动力和热力作用对大气环流、 气候变化和灾害性天气形成发展有十分重要的影响（吴统文等， 1998； 杨成等， 2020； 冯晓莉等， 2020）。高原的土壤水文过程不仅会影响土壤内部能量水分输送， 还可影响地表与大气之间的物质能量交换， 进而对高原区域的气候产生重要影响（Shukla et al， 1982； Yeh et al， 1984； Wen et al， 2003）。因此， 模式对于高原土壤水文过程的模拟能力是检验模式在青藏高原模拟效果的重要指标。
野外科考实验表明青藏高原的砾石含量非常高， 且分布广泛（Ohtsuka et al， 2008； 罗斯琼等， 2009； Arocena et al， 2012； Wu et al， 2012）。砾石是直径大于2 mm的岩石或矿物碎块， 由暴露在地表的岩石经机械风化作用形成， 常堆积在山麓或山前地带， 或由于岩石被水侵蚀破碎后， 经河流冲刷沉积产生。国内外众多学者研究表明砾石对于地表土壤水文过程有重要影响（Peck et al， 1979； Brakensiek et al， 1986； Abrahams et al， 1994； 李燕等， 2006； 何玉洁等， 2017）， 但是目前应用于区域气候模式中的土壤水热参数化方案并未考虑砾石的影响， 而忽略砾石的存在很难准确获知土壤实际水热状态（Brouwer et al， 2000）。值得关注的是潘永洁等（2015）发展了一套较适合青藏高原土壤特性的土壤参数化方案并在那曲站进行了应用与验证， 发现模拟效果有了一定的改进。然而目前针对改进土壤水热参数化方案在青藏高原的检验大多在单点开展， 并未比较高原不同区域的模拟效果（刘宜纲等， 2020）， 由于高原的砾石分布存在差异， 不同砾石含量对于地表与土壤水文过程的影响并非完全一致（Abrahams et al， 1991； 李燕等， 2006； 时忠杰等， 2008； 张志蓉等， 2008）， 此外， 考虑砾石影响后地表土壤水文过程变化的机制也有待研究。因此， 本文在耦合CLM4.5的区域气候模式RegCM4.7中分别应用原始土壤水热参数化方案与改进后的砾石参数化方案在青藏高原西部、 中部与东南部区域进行模拟， 并根据砾石分布特征在每个区域选取单点分析了两种方案模拟结果存在差异的原因。在此基础上利用中国陆面融合再分析数据（CRA-40）检验了砾石参数化方案在高原不同区域对于土壤含水量的模拟效果。
2 区域气候模式RegCM砾石参数化方案
土壤水热属性被认为是土壤中矿物属性与有机物属性的权重函数（Clapp et al， 1978； Cosby et al， 1984； Lawrence et al， 2008）。模式中矿质土壤水属性包括矿质土壤B指数 、 矿质土壤孔隙度 、 饱和矿质土壤基质势 及矿质土壤饱和导水率 ； 矿质土壤热属性包括矿质土壤固态导热率 、 矿质干土壤导热率 、 土壤容重 、 矿质土壤固态热容 和Kersten数 。原始土壤水热参数化方案仅考虑了沙石与粘土对矿质土壤水热属性的作用， 而改进的砾石参数化方案引入了砾石的影响。两种方案对比如下（表1 ）：
表1 两种参数化方案的对比Table 1 Comparison of the two parameterization schemes |
| 原始土壤水热参数化方案 | 砾石参数化方案 | |
|---|---|---|
| B指数 | ||
| 孔隙度 | ||
| 饱和基质势 | ||
| 饱和导水率 | ||
| 固态导热率 | ||
| 干导热率 | ||
| 土壤容重 | ||
| 固态热容 | ||
| Kersten数Ke ， i |
3 实验方案、 数据和评估方法
3.1 实验参数设置
选用国际理论物理中心（ICTP）开发的区域气候模式RegCM， 该模式目前最新版本为RegCM4.7， 其耦合的陆面模块为CLM4.5模式。模式中涉及的土壤水文过程均在土壤结构的前10层进行计算。利用该模式分别在高原西部、 中部与东南部区域进行模拟， 对应模拟中心选取为33°N， 82°E、 31°N， 89°E与28°N， 93°E。网格点水平分辨率设置为30 km。选用ECMWF的EIN15再分析资料作为模式初始与侧边界条件， 缓冲区大小设置为12， 边界条件每6 h更新一次， 每30 min调用一次陆面模块。根据砾石分布在高原西部、 中部与东南部分别选取了三个单点用于后续分析（简称为西部、 中部、 东南部）， 对应经纬度为32.69°N， 83.71°E、 31.49°N， 90.08°E与27.46°N， 93.03°E。模拟时段为2017年1月1日至2018年12月31日。
3.2 参照数据来源
参照数据选用同时段中国陆面融合再分析数据（CRA-40）， 根据中国气象局关于CRA-40的产品对比评估结果， CRA-40在中国区域与人工观测资料的偏差和均方根误差较小而相关系数较高， 表明CRA-40数据的精度优于其余产品。
3.3 评估方法
选用平均绝对误差（MAE）与相关系数（r）分别对砾石参数化方案与原始土壤水热参数化方案在高原不同区域模拟的土壤体积含水量进行检验， 公式如下：
式中： 为模拟结果； 为参照数据； 为样本数。
4 结果与分析
4.1 土壤水属性分析
在RegCM4.7模式中， B指数的改变会引起土壤导水率与地下水位的变化， 进而改变土壤体积含水量的变化率； 土壤孔隙度是影响水分传输的关键属性； 土壤基质势同样会影响土壤体积含水量的变化率； 土壤饱和导水率决定了土壤最大渗透能力。由于本文主要侧重于地表土壤水文过程的分析， 因此首先有必要比较两种方案在高原不同区域计算的土壤水属性的差异（图1 ）。本文将1.0 m以上的土壤深度定义为浅层土壤， 1.0~2.0 m范围内的土壤深度定义为中层土壤， 2.0 m以下定义为深层土壤。原方案在高原西部与中部计算得出的各土壤水属性在整层土壤中均为常数， 而砾石方案由于考虑了高原地区砾石含量高的实际情况， 计算得出的各土壤水属性随土壤深度而变化。
砾石方案在高原西部计算的矿质土壤B指数在整层土壤中均大于原方案， 两种方案的差异在土壤中层以下达到最大。砾石方案在高原中部计算的矿质土壤B指数在土壤浅层大于原方案， 在中层迅速回落并小于原方案， 两种方案的差异在浅层达到最大。砾石方案在高原东南部计算的矿质土壤B指数与原方案的垂直变化较为一致， 但在整层土壤中均小于原方案。砾石方案在高原西部计算的矿质土壤孔隙度在整层土壤中均小于原方案， 两种方案的差异在中层以下达到最大。砾石方案在高原中部计算的矿质土壤孔隙度在土壤浅层小于原方案， 中层以下快速上升， 两种方案的差异在浅层达到最大。砾石方案在高原东南部计算的矿质土壤孔隙度的垂直变化与高原中部较为相似。
对于饱和矿质土壤基质势， 砾石方案在高原西部的计算结果在整层土壤中均小于原方案， 两种方案的差异在中层以下达到最大。砾石方案在高原中部的计算结果在0~0.17 m与原方案保持一致， 在0.5 m降至最低， 随后迅速上升， 在土壤中层以下超过原方案。砾石方案在高原东南部的计算结果在整层土壤中均大于原方案， 两种方案的差异在土壤中层达到最大。对于矿质土壤饱和导水率， 砾石方案在高原西部的计算结果在整层土壤中均小于原方案， 两种方案的差异在土壤中层以下达到最大。砾石方案在高原中部的计算结果在土壤浅层小于原方案， 中层以下快速上升， 两种方案的差异在浅层达到最大。砾石方案在高原东南部的计算结果除土壤深层外均小于原方案。从以上对比可以发现， 两种方案在高原不同区域对于土壤水属性计算结果的差异有所不同， 由此导致在高原不同区域对于地表土壤水文过程的模拟差异同样会存在不同。
4.2 地表水文特征分析
4.2.1 地表水分输入
RegCM4.7模式中地表水分输入由到达地面的液态降水与融化的雪水组成， 随后在地表径流、 地表水存储和地表入渗之间分配。对比两种方案在高原不同区域模拟的地表水分输入强度发现（图2 ）， 两种方案模拟的地表水分输入的起止时段较为一致， 高强度时段主要集中在雨季（5 -9月）。在高原东南部模拟的地表水分输入强度明显大于其他区域。砾石方案在高原各个区域模拟的地表水分输入强度在每年的雨季较原方案均有所减小， 其中高原西部平均减小5.8×10-6 mm·s-1， 高原中部平均减小4.7×10-6 mm·s-1， 高原东南部平均减小1.2×10-5 mm·s-1。
4.2.2 地表径流
RegCM4.7模式使用的地表径流模型为SIMTOP （Niu et al， 2005）， 该模型中地表径流与地表水分输入、 饱和面积分数f sat成正比。f sat与地下水位有关， 当地下水位下降时f sat减小， 而地下水位的变化又受给水度的影响， 由于给水度是土壤水属性的函数， 从上文可知， 砾石方案模拟的地下水位较原方案必然会发生变化（图略）。由于地下水位的变化， 两种方案在高原不同区域模拟的f sat存在差异（图3 左列）。原方案在高原西部模拟的f sat除雨季个别时段外全年稳定在0.046， 砾石方案模拟的f sat的时间变化与原方案一致， 但在原方案出现波动的时段砾石方案的模拟值要大于原方案。原方案在高原中部模拟的f sat全年稳定在0.038， 而砾石方案模拟的f sat则在雨季转为波动状态， 且模拟值明显大于原方案。砾石方案在高原东南部模拟的f sat较原方案在大部分时段略有减小。同时还发现， 在高原不同区域模拟的f sat其时间变化特征存在明显的区别， 在高原西部与中部模拟的f sat通常较为稳定， 只在每年雨季才会出现短暂波动， 而在高原东南部模拟的f sat全年均处于波动状态， 原方案与砾石方案均模拟出了这种区别。
砾石方案模拟的地表径流强度对f sat较为敏感， f sat的波动会显著引起地表径流强度的变化（图3 右列）。值得注意的是， 在高原西部， 砾石方案模拟的地表径流强度即使在f sat出现波动的时段也并未大于原方案， 表明增大的f sat并未抵消地表水分输入减少所带来的影响。在高原中部， 受明显增大的f sat的影响， 砾石方案模拟的地表径流强度在雨季明显大于原方案。在高原东南部， 由于砾石方案模拟的地表水分输入与fsat均减小， 导致其模拟的地表径流强度在大部分时段小于原方案， 平均减小1.1×10-6 mm·s-1。
4.2.3 地表入渗
模式定义了土壤最大渗透能力与淹没面积分数f h2osfc用于模拟地表入渗与地表水存储。当地表入渗强度超过土壤最大渗透能力时就会产生霍顿径流补充地表水存储， 土壤最大渗透能力越小， 霍顿径流越大。土壤最大渗透能力与冰阻系数、 非饱和面积分数及饱和导水率成正比， 而冰阻系数又与土壤孔隙度成正比。图4 （右列）给出了两种方案在高原不同区域模拟的f h2osfc， 砾石方案在高原西部与中部的模拟值在原方案模拟的f h2osfc发生变化的时段偏大， 在东南部除个别时段外为0。由于砾石方案在高原不同区域计算的土壤孔隙度与饱和导水率在大部分土壤深度中均小于原方案， 而砾石方案在高原西部与中部模拟的非饱和面积分数较原方案偏小， 导致砾石方案在高原西部与中部模拟的土壤最大渗透能力低于原方案， 因此这两个区域的霍顿径流将大于原方案。而在高原东南部模拟的土壤最大渗透能力取决于土壤孔隙度、 饱和导水率与非饱和面积分数之间的竞争结果。
地表入渗是剩余地表水分输入、 非淹没面积分数与霍顿径流的函数， 当剩余地表水分输入与非淹没面积分数越大、 霍顿径流越小时地表入渗越大。由于砾石方案在高原西部与中部模拟的霍顿径流大于原方案， 且其在上述区域模拟的剩余地表水分输入与非淹没面积分数均小于原方案， 导致砾石方案在高原西部与中部模拟的地表入渗强度在雨季小于原方案（图5 ）， 其中高原西部平均减小4.1×10-6 mm·s-1， 高原中部平均减小5.9×10-6 mm·s-1。砾石方案在高原东南部模拟的地表入渗强度同样小于原方案， 平均减小1.2×10-5 mm·s-1， 由于两种方案在高原东南部模拟的非淹没面积分数除个别时段外相等， 由此推测地表水分输入的减小是导致砾石方案模拟的地表入渗强度偏小的一个主要原因。
4.3 土壤水分输送特征分析
在模式中， 地表入渗的变化会对土壤顶层的含水量造成直接影响。图6 ~8给出了两种方案在高原不同区域模拟的土壤含水量及其差异随土壤深度变化的时间剖面。在高原西部［图6 （c）］， 砾石方案模拟的浅层与深层土壤含水量较原方案有所减少， 但中层土壤含水量却高于原方案。在高原中部［图7 （c）］， 砾石方案模拟的整层土壤含水量较原方案均有所减少， 尤其在0.4~1.0 m的土壤深度范围内减少最为明显。在高原东南部［图8 （c）］， 除土壤底层外， 砾石方案模拟的土壤含水量较原方案同样偏少， 两种方案的差异在冬季的浅层土壤中最为明显。土壤水分向下传输的增大可能是砾石方案在底层模拟的土壤含水量偏高的原因。结合高原的砾石分布可以发现（即砾石含量在高原西部多集中在土壤中层以下， 在高原中部多集中在土壤浅层0.4 m附近， 在高原东南部整层均较少）， 较低的砾石含量有利于砾石方案模拟的土壤含水量向下传输， 导致该层土壤体积含水量减少。但值得注意的是， 在高砾石含量的情况下， 砾石的作用在高原西部与中部呈现相反的结果， 砾石可以阻缓或促进土壤水分的下渗。砾石方案在高原西部模拟的土壤含水量在中层偏多， 表明此区域土壤中较高的砾石含量会阻缓土壤水分的向下传输。而砾石方案在高原中部模拟的土壤含水量在0.4~1.0 m的土壤范围内减少最为明显， 表明此区域土壤中较高的砾石含量会促进土壤水分的向下传输。
图6 两种方案在青藏高原西部模拟的土壤含水量随土壤深度变化的时间剖面（单位： mm3·mm-3）Fig.6 The temporal profile of soil moisture changes with soil depth simulated by the two schemes over the western Qinghai-Xizang Plateau.unit： mm3·mm-3 |
图7 两种方案在青藏高原中部模拟的土壤含水量随土壤深度变化的时间剖面（单位： mm3·mm-3）Fig.7 The temporal profile of soil moisture changes with soil depth simulated by the two schemes over the central Qinghai-Xizang Plateau.unit： mm3·mm-3 |
4.4 模拟效果检验
图9 两种方案在青藏高原西部模拟的各层土壤含水量及对应参照数据的时间变化图9（d）为标准化处理后的结果， 其值不代表真实数据， 仅为直观展示对比结果 Fig.9 The temporal variation of soil moisture simulated by the two schemes and the corresponding reference data over the western Qinghai-Xizang Plateau.In Fig.9（d）， the result after standardization does not represent the truth values， just to visualize the comparison results |
表2 两种方案模拟的土壤含水量在高原区域的检验结果Table 2 Statistics of the reference data and the simulated soil moisture over the Qinghai-Xizang Plateau |
| 区域 | 土壤深度/m | MAE/（mm3·mm-3） | R | |||
|---|---|---|---|---|---|---|
| 原方案 | 砾石方案 | 原方案 | 砾石方案 | |||
| 高原西部 | 0.1 | 0.071 | 0.048 | 0.35 | 0.42 | |
| 0.4 | 0.095 | 0.025 | 0.57 | 0.57 | ||
| 1.0 | 0.034 | 0.030 | -0.95 | -0.25 | ||
| 2.0 | 0.056 | 0.058 | 0.52 | 0.84 | ||
| 高原中部 | 0.1 | 0.046 | 0.056 | 0.82 | 0.86 | |
| 0.4 | 0.036 | 0.049 | 0.63 | 0.69 | ||
| 1.0 | 0.051 | 0.029 | 0.45 | -0.45 | ||
| 2.0 | 0.033 | 0.014 | 0.46 | 0.42 | ||
| 高原东南部 | 0.1 | 0.066 | 0.035 | 0.78 | 0.79 | |
| 0.4 | 0.038 | 0.045 | 0.75 | 0.83 | ||
| 1.0 | 0.048 | 0.031 | 0.65 | 0.78 | ||
| 2.0 | 0.089 | 0.021 | 0.64 | 0.69 | ||
5 结论与讨论
在耦合CLM4.5的区域气候模式RegCM4.7中分别应用原始土壤水热参数化方案与改进后的砾石参数化方案在青藏高原西部、 中部与东南部区域进行模拟， 并根据砾石分布特征在每个区域选取单点分析了两种方案模拟结果存在差异的原因。在此基础上利用中国陆面融合再分析数据（CRA-40）检验了砾石参数化方案在高原不同区域对于土壤含水量的模拟效果。通过对比分析， 得到以下结论：
（1） 在应用砾石参数化方案的RegCM4.7模式中， 地表水分输入较其他变量是影响地表土壤水文过程的主要因素。
（2） 在应用砾石参数化方案的RegCM4.7模式中， 较低的砾石含量有利于土壤水分向下传输， 较高的砾石含量在高原不同区域对于土壤水分向下传输的影响作用有所不同： 在本次实验中， 在高原西部土壤中较高的砾石含量会阻缓土壤水分的向下传输， 而高原中部土壤中较高的砾石含量会促进土壤水分的向下传输。
（3） 砾石参数化方案在高原不同区域对于土壤水分的模拟效果存在差异， 在高原东南部的模拟效果提升最为明显。
本文还存在一些不足， 首先只比较了两种方案对于地表土壤水文过程模拟结果的差异， 其次在高原选取的点位过少。后续研究中将更加全面的比较两种方案的模拟结果， 并选取更密集的高原观测数据进行检验。