2023, Vol. 55
2. 武汉大学水资源与水电工程科学国家重点实验室, 武汉 430072
冻土是广泛存在于我国北方寒冷地区的一种特殊水土体系，其明显区别于非冻土的水热状况和迁移特征[1]。一方面，冻融过程中的水势是以温度势和基质势为主导，冰的存在导致土壤导水率的降低从而影响入渗过程；另一方面，由于冻土中水的运移引起了土壤热特性参数在时空上的变化，从而影响土壤中热流和温度分布。土壤水的相变和对流传热，使土壤水、热迁移具有强烈的耦合性，二者在土壤冻融过程中的分布和迁移密不可分，共同构成了受多种因素综合影响的复杂系统[2]。寒区冻土的水热状况影响作物的生长发育，对于越冬作物(如冬小麦)，越冬期不同土壤水热状况将直接影响作物的出苗率，并间接影响作物幼苗发育和物质积累，对作物生长起到了承上启下的作用[3]；而对于非越冬作物(如水稻、大豆、玉米等)来说，越冬末期的土壤水热状况影响作物种植的初始土壤环境条件，从而影响作物的种植时间、种子萌发和根系生长等重要生长过程，对作物生长起到了引导作用。因此，开展寒区冻融期土壤热迁移数值模拟研究，掌握冻融过程中的土壤热状况和变化规律，分析模型输入参数的不确定性，识别关键因素，具有重要实践意义。
由美国农业部Flerchinger博士等[4–5]建立的SHAW(Simultaneous heat and water)模型被广泛用于冻融水热特征及陆面生态环境变化等的研究。国内外学者对该模型在寒区水热过程模拟的适用性进行了探究。DeGaetano等[6]将该模型用于寒区最大冻深模拟，与实测值对比，表明模型可以较好地模拟研究区的最大冻深；Lu等[7]基于SHAW模型模拟了盐渍地不同覆盖条件下冻融土壤水热运移规律，结果表明模型可用于土壤温度模拟，且土壤温度的模拟效果均好于土壤含水率；Chen等[8]基于SHAW模型探究了不同地下水位、不同土壤质地下冻融土壤水热过程。多数学者仅对模型进行应用，但是冻融期土壤水热迁移问题的复杂性、影响因素的多变性以及强烈的耦合性，使得模型参数不确定性对模拟值精度影响作用较大，导致模拟结果的不确定性增加，严重影响模拟的稳定性与准确度。可见，定量评价模型输入参数不确定性是SHAW模型合理预测的基础和前提。目前用于不确定性分析的方法有OAT(One factor at a time)[9]、FOEA(First-order error analysis)[10]、蒙特卡罗法(Monte carlo)[11]等。许多学者已将这些方法广泛运用于不同陆面过程模型的校正和不确定性分析中。对于SHAW模型，Flerchinger等[12]首次采用OAT法分析了SHAW模型中土壤水力学参数对土壤含水率输出结果的敏感性影响。本文在SHAW模型研究中，引入拉丁超立方取样(Latin hypercube sampling，LHS)方法，该方法作为定额抽样法的延伸，相比其他抽样方法，其优点在于能够确保采样点覆盖整个输入随机变量的分布区间，同时可克服如Monte carlo法计算成本过高的缺点[13]。
本研究以松嫩平原黑土区为研究区域，考虑初始积雪覆盖条件，利用LHS方法对SHAW模型输入参数进行取样并对土壤冻结深度和温度的输出进行不确定性分析，同时通过标准秩逐步回归分析量化各参数在不同深度土壤剖面的不确定性贡献率。研究成果有效减少模型的率定、验证工作，为模型参数间相互作用、相互影响等规律的探索提供参考资料。
1 材料与方法 1.1 研究区概况试验区位于松嫩平原黑土区东北农业大学综合试验场，地理位置为126°45'32″E、45°44'41″N，属中温带大陆性季风气候，冬季寒冷干燥，受极地大陆气团控制，冷气团长期滞留。年平均气温3.6 ℃，冬季多年平均气温为–14.2 ℃，最冷月平均气温–19.6 ℃。多年平均积雪覆盖期为110 d，降雪期集中在每年11月至次年1月，历年平均降雪量为23.6 mm，年际变化介于6.6 ~ 35.4 mm，最大冻深为180 cm。试验区位于平原腹地，地势平坦，无植被覆盖，经人工取样分析可知，该地区0 ~ 60 cm土层为黑色壤土，60 cm以下为黑色黏土，土壤理化性质见表 1。
|
表 1 不同土层的物理性质 Table 1 Physical properties of different soil layers |
在东北农业大学试验场设置了土壤温度观测点，测定设备为土壤温湿度传感器(ET100，北京东方生态，中国)，测定深度为0 ~ 20、20 ~ 40、40 ~ 60、60 ~ 100、100 ~ 140、140 ~ 180 cm土层深度的温度。除此之外，地块还设置冻土器用来监测土壤冻结深度变化情况，通过人工测量的方式来获取具体数据。
1.3 模型描述SHAW模型为一维水热耦合模型，最初是由Flerchinger和Saxton[4]在1989年建立，并用于模拟土壤冻融过程中热量、水量、溶质通量的耦合迁移运动，模型的基本原理如图 1所示。本研究采用SHAW模型模拟冬季实验场内积雪覆盖条件下的土壤垂向剖面水热动态迁移过程，并探究输入参数对土壤温度以及冻结深度输出结果的影响。模型要求输入的信息包括：积雪、土壤温度和土壤含水量的初始状态；每日或每时的气象资料(气温、风速、湿度、降水和太阳辐射)；模拟地点的位置概况(坡度、坡向、纬度和表面糙率参数)以及描述积雪、土壤剖面物理特征的参数[14]。本研究选用2013年12月份的每日气象资料作为驱动，同时输入初始时刻的积雪特性参数(雪厚、雪温、雪密度)与初始、结束时刻的土壤温度的实测值。模型将土壤剖面划分为6个土层(0 ~ 20、20 ~ 40、40 ~ 60、60 ~ 100、100 ~ 140、140 ~ 180 cm)，每层由一个节点表示，并根据前期土壤颗粒分析结果建立两层系统，即0 ~ 20、20 ~ 40、40 ~ 60 cm为L1，60 ~ 100、100 ~ 140、140 ~ 180 cm为L2。将土壤特性参数分别视为一致输入模型。模拟时间为当年12月2日(第336天)至12日(第346天)，输出时间步长为24 h。
|
图 1 SHAW模型物理系统描述图 Fig. 1 Description diagram of SHAW model physical system |
本研究采用基于参数抽样的标准秩逐步回归方法进行不确定性分析，其步骤分析如下：①将输入数视为随机变量，确定各随机变量的概率分布和取值范围；②根据假设的分布与预估的取值区间，采用LHS方法生成若干组随即输入参数；③将生成的参数输入模型运行若干次；④对输出目标进行不确定性分析；⑤将输入参数、输出目标分别作为自变量和因变量进行标准秩逐步回归，并定量表示各参数敏感性强弱。
1.4.1 参数LHS抽样LHS方法作为定额抽样法的延伸，相比传统抽样方法，其优点在于能够确保采样点覆盖整个输入随机变量的分布区间[15]。在LHS方法中，记输入变量为k，抽样次数为n。对于每个输入变量X1，X2，…，Xk，按照假设的概率密度函数将其取值区间等概率地划分为n个互不重叠的间隔区间，然后在每个间隔内按其概率分布随机抽样。再将X1的n个抽样取值与X2的n个抽样取值随机配对，以此类推直到第k个变量，如此便得到n个k维输入函数组值[16]。
1.4.2 标准秩逐步回归标准秩逐步回归法将输入与输出数据转化为相关秩系数，即通过对输入和输出数据分别按照数值大小排序，再对排序后的数据进行标准化回归。该方法解决了复杂模型输入与输出之间的非线性关系，从而提高了回归估计的精度[17]。计算得到的指标包括标准秩回归系数(Standardized rank regression coefficients，SRRC)和可决系数R2。SRRC2表示参数对不确定性的影响大小，其值越大表明参数越敏感；R2表示回归的可行性，一般认为当R2 > 0.7时回归是可行的[18]。若各输入参数间相互独立，二者有如下关系：
| $ {R^{\text{2}}} = \sum\limits_{i = 1}^k {{\text{SRR}}{{\text{C}}^2}} $ | (1) |
SHAW模型的输入参数可分为4类：初始及边界条件，直接影响表层热迁移的参数，土壤层及残渣层的热传导参数，土壤特性参数[19]。根据上述参数分类，结合实验地气候特点与数值模拟结果分析，选择1个外部初始条件积雪厚度和4个土壤水力特性参数，即Cambell孔径大小分布指数、空气进入势、饱和导水率、饱和体积含水量作为模型输入参数。
通过查询参考文献[12，20]和气象统计资料，表 2列出了这5个输入参数的取值区间，并假定各参数符合均匀分布且相互独立。根据建立的土壤剖面两层系统，并将第一层系统命名为L1，第二层系统命名为L2，以区分两层系统的土壤特性参数。SHAW模型依据雪厚将积雪分为若干层，本文将初始积雪分为3层：①若积雪覆盖深度不足2 cm，则分为一层；②若积雪层厚度介于2 ~ 4 cm，则视为双层，同时将第二层厚度定为2 cm，其余划归于第一层；③若积雪层厚度介于4 ~ 6 cm，则视为三层，底层和中层定为2 cm，剩余作为表层厚度；④如果积雪层厚度大于6 cm，亦视为三层，表层定为2 cm，剩余均分给中层和底层，但中层不能超过18 cm，多余的雪层加在底层上。
|
|
表 2 LHS抽样分析模型输入参数取值范围 Table 2 Value ranges of input parameters of LHS sampling analysis model |
LHS的抽样次数至少为2k (k为因素个数)次[15]。根据试算，本研究采用LHS抽样方法生成300组输入参数，每组参数都输入SHAW模型运行一次，最后得到300组含有特定目标的输出结果，利用此300组输出数据定量研究5个参数对模型模拟结果不确定性的影响。
2 结果与分析将自然积雪边界参数输入模型，得到2013年11月29日至2014年4月15日的土壤冻结融化深度、温度的模拟结果，其中实测值的观测时间为2013年11月1日至实际融通时间。
2.1 土壤冻融过程和温度模拟 2.1.1 土壤冻融过程模拟由图 2可知，在模拟时间段内，自然积雪覆盖下土壤冻融模拟结果与实测值较吻合，最大冻结深度与实测值相差5 cm。土壤融通的实际时间为2014年4月28日，由模拟曲线可知，在模拟时间段内土壤均未融通，且土壤开始融化时间出现一定滞后，这可能与模型初始输入的积雪分层情况、模型模拟过程中忽略了不同深度雪密度的改变有关。
|
图 2 模拟与实测冻融深度 Fig. 2 Simulated and measured freeze-thaw depths |
不同深度土层温度的模拟值与实测值对比曲线如图 3所示。从图中可以看出，自然积雪覆盖下6个土层的温度模拟值变化趋势与实测值基本一致，能够反映土壤冻融期温度的变化规律。
|
图 3 自然积雪覆盖下不同深度土壤温度模拟值与实测值对比 Fig. 3 Comparison of simulated and measured soil temperatures at different depths under natural snow cover |
模拟值与实测值的均方根误差(root mean square error，RMSE)、最大绝对误差、平均绝对误差计算结果如表 3。从模拟的精度来看，大多数模拟值都小于实测值，RMSE除20 ~ 40、40 ~ 60 cm剖面分别为2.24、2.17 ℃外，其余剖面均小于1.6 ℃；且平均绝对误差均小于2 ℃，平均绝对误差最大值出现在20 ~ 40 cm土壤剖面，其值为1.89 ℃。
|
|
表 3 不同土层深度土壤温度模拟值与实测值的误差 Table 3 The errors between simulated and measured soil temperatures at different depths |
以模型模拟的第341天及2013年12月7日的土壤冻结深度、土壤温度为输出目标，针对参数组合产生的300组输出数据进行不确定性与敏感性分析。因各土层R2不同，为比较不同参数敏感型的相对大小，由SRRC2和R2的上述关系，定义不确定性贡献率为参数的SRRC2所占R2比例。
2.2.1 以土壤冻结深度为输出目标图 4为2013年12月7日的冻结深度的输出分布。300次输出的平均值为7.35 cm，由图可见，输出分布集中在2.7 ~ 5.6 cm，整体分布较均匀，其变差系数为2.59，这表明参数不确定性对土壤冻结深度的输出有较大影响。
|
图 4 土壤冻结深度输出分布图 Fig. 4 Output distribution map of soil freezing depth |
标准秩逐步回归及不确定性贡献率计算结果表明R2为0.812 9，大于0.7，表明回归结果是可行的。在模拟时段内冻结深度较浅，由图 5可见，L2土壤特性参数不确定性总贡献率较小(< 6%)，故在分析时忽略L2的不确定性贡献影响。以土壤冻结深度为输出目标，参数不确定性贡献率排序如下：初始积雪厚度(54.89%) > 空气进入势(15.97%) > 饱和含水率(13.83%) > Cambell孔径大小分布指数(8.94%) > 饱和导水率(< 1%)。
|
(图中h表示为初始积雪厚度；b表示为Cambell孔径分布指数；ψe表示为空气进入势；Ks表示为饱和导水率；θs表示为饱和含水率；L1表示为0 ~ 20、20 ~ 40、40 ~ 60 cm的不同深度土层系统；L2表示为60 ~ 100、100 ~ 140、140 ~ 180 cm的不同深度土层系统，下同) 图 5 以土壤冻结深度为输出目标参数不确定性贡献率 Fig. 5 Contribution rate of uncertainty of output target parameter with soil freezing depth |
可知，初始积雪厚度对土壤冻结深度的输出起主导影响，空气进入势与饱和含水率的不确定性贡献率接近，饱和导水率不确定性贡献率仅为0.7%，表明饱和导水率不是敏感参数。姚闯等人[20]在黄河源区积雪对冻土水热过程模拟研究中表明，积雪在土壤冻结期，具有一定的保温作用，而当土壤上积雪深度减少，保温作用随之减弱，土壤热量输送至大气。李杨等人[21]在北疆地区积雪与冻土变化的特征分析研究中，当积雪深度在0 ~ 15 cm时，最大冻土值上升较快；15 ~ 40 cm时最大冻土值缓慢下降，显示最大冻结深度随积雪深度逐渐增加。付强等人[22]在对不同积雪覆盖条件下冻融土壤水分运动规律研究中指出由于积雪的低热导性，积雪覆盖对土壤有良好的缓冲和保温作用。同时积雪能够减缓土壤的冻结速度，进而影响土壤水分的迁移过程。
2.2.2 以土壤温度为输出目标如图 6反映了6个土层温度的输出分布，总体看，土壤温度的输出波动远小于冻结深度的输出。0 ~ 20、140 ~ 180 cm深度土层温度输出最集中，这表明参数不确定性对土壤上下边界温度影响极小。在0 ~ 20 cm深度土壤温度的300次输出中出现了38次离群值，这可能是由于上边界积雪厚度的输入波动较大造成的。40 ~ 60、60 ~ 100、100 ~ 140 cm土壤温度输出分布亦较为集中，其四分位距分别为0.3、0.7、0.4 ℃，这表明参数不确定性对此3个深度土层温度的输出有一定影响但其影响不大。谢婷等人[23]在对不同深度土壤温度模拟的研究中指出，SHAW模型模拟土壤热量取决于周围环境土壤水热状况，随着深度增加，周围环境温湿度变化幅度越小，SHAW模型模拟土壤温度数值越精确。郭志强等人[24]在季节性冻融土壤水热耦合运移模拟的研究中发现地表由于受降水、蒸发、日照等影响，表层温度的模拟效果较差，而40和80 cm土壤温度实测值与模拟值吻合，模拟结果可靠。
|
图 6 不同深度土壤温度输出分布 Fig. 6 Output distribution of soil temperatures at different depths |
表 4为标准秩回归的结果。由于140 ~ 180 cm土层温度的输出值无变化，无法对其进行大小排序，根据标准秩逐步回归的定义，其回归不可行。由表可知，只有40 ~ 60、60 ~ 100 cm土层R2大于0.7，这表明回归效果较差，输入参数与输出之间关系不显著，可见选取的输入参数对土壤温度的输出较为不敏感。Flerchinger和Hardegree[25]在对土壤温度输出效果优化的研究中指出，若输入的土壤质地、容重等影响土壤热特性的参数在合理区间内，绝大多数参数对土壤温度的输出都不敏感。付强等人[26]在对积雪覆盖下土壤热状况及其对气象因素的响应研究中发现，随着积雪深度增加和密度增大，即覆盖物对地表与环境之间阻碍程度增加，环境温度、大气辐射对土壤温度关联度有所下降，土壤温度与大气温度的相关性减弱。
|
|
表 4 以土壤温度为输出目标标准秩逐步回归结果(SRRC2) Table 4 Stepwise regression results of standard rank with soil temperature as output target |
1) 土壤在冻融过程中表现出单向冻结双向融化，在自然积雪覆盖条件下，SHAW模型的模拟结果能够较好地反映出这一规律。土壤开始融化时间的模拟值均出现一定滞后。SHAW模型的土壤温度模拟值与实测值具有较一致的变化趋势，这表明SHAW模型的数值模拟结果能够反映出自然积雪覆盖下土壤冻融过程中的温度变化规律。
2) 在SHAW模型模拟积雪覆盖条件下冻融土壤热迁移过程中，参数不确定性的影响是客观存在的。其中，初始积雪厚度对模拟目标土壤冻结深度的输出影响起主导作用，空气进入势与饱和含水率的不确定性贡献率接近，饱和导水率对土壤冻结深度的输出敏感性较弱。
3) SHAW模型可用于积雪覆盖下土壤热过程动态变化的影响，基于LHS抽样的标准秩逐步回归法是对模型进行不确定性分析和参数敏感性分析的有效方法。
| [1] |
隽英华, 田路路, 刘艳, 等. 冻融作用对农田黑土有机氮组分的调控效应[J]. 土壤, 2020, 52(2): 262-270 (  0) |
| [2] |
朱龙祥, 范昊明, 马仁明. 冻融循环和土壤含水率对棕壤崩解特性的影响[J/OL]. 土壤学报, 2022: 1–14. (2022-01-22). https://kns.cnki.net/kcms/detail/32.1119.P.20220120.1355.004.html.
( 0) |
| [3] |
彭新华, 王云强, 贾小旭, 等. 新时代中国土壤物理学主要领域进展与展望[J]. 土壤学报, 2020, 57(5): 1071-1087 ( 0) |
| [4] |
Flerchinger G N, Saxton K E. Simultaneous heat and water model of a freezing snow-residue-soil system I. theory and development[J]. Transactions of the ASAE, 1989, 32(2): 565-571 DOI:10.13031/2013.31040 ( 0) |
| [5] |
Flerchinger G N, Kustas W P, Weltz M A. Simulating surface energy fluxes and radiometric surface temperatures for two arid vegetation communities using the SHAW model[J]. Journal of Applied Meteorology, 1998, 37(5): 449-460 DOI:10.1175/1520-0450(1998)037<0449:SSEFAR>2.0.CO;2 ( 0) |
| [6] |
DeGaetano A T, Cameron M D, Wilks D S. Physical simulation of maximum seasonal soil freezing depth in the United States using routine weather observations[J]. Journal of Applied Meteorology, 2001, 40(3): 546-555 DOI:10.1175/1520-0450(2001)040<0546:PSOMSS>2.0.CO;2 ( 0) |
| [7] |
Lu X H, Li R P, Shi H B, et al. Successive simulations of soil water-heat-salt transport in one whole year of agriculture after different mulching treatments and autumn irrigation[J]. Geoderma, 2019, 344: 99-107 DOI:10.1016/j.geoderma.2019.03.006 ( 0) |
| [8] |
Chen J F, Gao X G, Zheng X Q, et al. Simulation of soil freezing and thawing for different groundwater table depths[J]. Vadose Zone Journal, 2019, 18(1): 1-14 ( 0) |
| [9] |
杨军军, 高小红, 李其江, 等. 湟水流域SWAT模型构建及参数不确定性分析[J]. 水土保持研究, 2013, 20(1): 82–88, 93 ( 0) |
| [10] |
余红, 沈珍瑶. 三峡水库大宁河流域非点源污染参数的不确定性分析[J]. 中国环境科学, 2007, 27(4): 554-558 DOI:10.3321/j.issn:1000-6923.2007.04.025 ( 0) |
| [11] |
刘艳丽, 周惠成, 张建云. 不确定性分析方法在水库防洪风险分析中的应用研究[J]. 水力发电学报, 2010, 29(6): 47-53 ( 0) |
| [12] |
Flerchinger G N, Caldwell T G, Cho J, et al. Simultaneous heat and water (SHAW) model: Model use, calibration, and validation[J]. Transactions of the ASABE, 2012, 55(4): 1395-1411 DOI:10.13031/2013.42250 ( 0) |
| [13] |
McKay M D, Beckman R J, Conover W J. Comparison of three methods for selecting values of input variables in the analysis of output from a computer code[J]. Technometrics, 1979, 21(2): 239-245 ( 0) |
| [14] |
Flerchinger G N. The simultaneous heat and water (SHAW) model: User's manual. technical report NWRC-2000-10[R]. USDA-ARS, Northwest Watershed Research Center, Boise.
( 0) |
| [15] |
Helton J C, Davis F J. Latin hypercube sampling and the propagation of uncertainty in analyses of complex systems[J]. Reliability Engineering & System Safety, 2003, 81(1): 23-69 ( 0) |
| [16] |
Pan F, Zhu J T, Ye M, et al. Sensitivity analysis of unsaturated flow and contaminant transport with correlated parameters[J]. Journal of Hydrology, 2011, 397(3/4): 238-249 ( 0) |
| [17] |
Saltelli A, Ratto M, Tarantola S, et al. Sensitivity analysis practices: Strategies for model-based inference[J]. Reliability Engineering & System Safety, 2006, 91(10/11): 1109-1125 ( 0) |
| [18] |
Rerchinger G N. Sensitivity of soil freezing simulated by the shaw model[J]. Transactions of the ASAE, 1991, 34(6): 2381-2389 ( 0) |
| [19] |
Mertens J, Madsen H, Kristensen M, et al. Sensitivity of soil parameters in unsaturated zone modelling and the relation between effective, laboratory and in situ estimates[J]. Hydrological Processes, 2005, 19(8): 1611-1633 ( 0) |
| [20] |
姚闯, 吕世华, 李照国, 等. 黄河源区积雪对冻土水热过程影响的数值模拟[J]. 高原气象, 2020, 39(6): 1167-1180 ( 0) |
| [21] |
李杨, 李建刚, 刘艳, 等. 北疆地区积雪与冻土变化的特征分析[J]. 水土保持研究, 2015, 22(5): 342-348 ( 0) |
| [22] |
付强, 蒋睿奇, 王子龙, 等. 不同积雪覆盖条件下冻融土壤水分运动规律研究[J]. 农业机械学报, 2015, 46(10): 152-159 ( 0) |
| [23] |
谢婷, 马育军, 杨晨. 基于SHAW模型的青海湖流域土壤温度模拟[J]. 土壤, 2022, 54(1): 161-168 ( 0) |
| [24] |
郭志强, 彭道黎, 徐明, 等. 季节性冻融土壤水热耦合运移模拟[J]. 土壤学报, 2014, 51(4): 816-823 ( 0) |
| [25] |
Flerchinger G N, Hardegree S P. Modelling near-surface soil temperature and moisture for germination response predictions of post-wildfire seedbeds[J]. Journal of Arid Environments, 2004, 59(2): 369-385 ( 0) |
| [26] |
付强, 侯仁杰, 王子龙, 等. 积雪覆盖下土壤热状况及其对气象因素的响应研究[J]. 农业机械学报, 2015, 46(7): 154-161 ( 0) |
2. State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan 430072, China