2024, Vol. 61
表1(Table 1)
2. 中国科学院大学, 北京 100049
2. University of Chinese Academy of Sciences, Beijing 100049, China
土壤有机碳（Soil organic carbon，SOC）作为土壤肥力形成的基础，不但影响土壤质量和粮食产量，而且在全球碳循环中扮演重要角色[1-3]。因此，准确把握SOC时空演变规律对于土壤资源高效可持续利用，以及应对气候变化和碳中和决策等均具有重要意义。
过程模型是SOC时空演变建模的常用方法。经典SOC过程模型（例如RothC、Century）已广泛应用于区域及国家尺度SOC时空演变建模，同时，经典SOC过程模型结构也被嵌入陆面过程模型（例如Community Land Model，CLM；Joint UK Land Environment Simulator，JULES），以实现全球尺度碳循环预测建模[4]。SOC过程模型大多在田块等较小的尺度上建立[5-6]，其内部参数的默认值则是基于模型开发环境条件下的长期定位试验数据来确定[7]。然而，由于气候、土壤、管理措施等的空间异质性，过程模型内部参数的默认值并非总是最优的。某些环境条件下，直接采用过程模型参数的默认值来建模SOC动态可能导致模拟结果产生极大的不确定性[8-9]，因此，在基于过程模型的SOC建模中，对过程模型内部参数进行优化校正已经成为降低模型预测不确定性、提升预测精度的重要手段。例如，Gurung等[10]基于欧洲及北美的19个长期定位试验获取的491个SOC观测数据，采用贝叶斯校正方法对DayCent模型（Century模型的日步长版本）中影响SOC过程的17个敏感参数进行优化，结果表明，SOC预测结果的不确定性由参数优化前的211%降低至参数优化后的33%。Hararuk等[11]采用马尔科夫链蒙特卡洛（Markov Chain Monte Carlo，MCMC）方法对CLM模型内部影响SOC过程的20个参数进行优化研究发现，参数优化后CLM模型预测的全球SOC能够解释观测值变异性的41%，较采用模型默认参数值预测精度提高了52%。Xu等[12]基于CLM-CASA（Community Land Model with Carnegie-Ames-Stanford Approach biogeochemistry sub-model，CLM-CASA）模型预测全球SOC的精度分析则表明，默认参数值下CLM-CASA模型预测结果仅能解释SOC观测值变异性的33%，而采用SBO（Surrogate-based Optimization，SBO）方法优化模型中影响SOC过程的20个内部参数后，模型解释的SOC变异性则提升至41%。
微生物是土壤碳循环的重要驱动者。一方面，微生物通过分解SOC获得自身生长所需要的养分和能量[13]；另一方面，微生物死亡后，其残留物是SOC的重要组成部分[14]。经典SOC过程模型考虑了微生物在SOC周转过程中的作用，例如，RothC和Century模型分别采用微生物生物量库（Microbial Biomass，BIO）和活性SOC库（包括活的微生物和微生物产物）表征微生物及其残体对SOC库的贡献。但在经典SOC过程模型（例如RothC、Century）中，微生物作为SOC分解者的作用仅被隐含在不同SOC亚库的分解速率常数中，未明确表达。而近些年新出现并应用于SOC建模的微生物模型，例如MIMICS（Microbial-Mineral Carbon Stabilization model，MIMICS）、MEND（Microbial ENzyme-mediated Decomposition，MEND）等，则直接将土壤碳周转和微生物生物量及生理机能耦合，从而明确反映了微生物在SOC分解和稳定化中的作用[15]。然而，微生物模型内部参数的默认值主要基于其开发环境条件下的短期田间试验以及实验室培养试验确定[9]，且大多为常值参数。但常值参数可能难以表征区域、国家、全球等大尺度上土壤碳循环动态显著的空间异质性[16-17]，从而对SOC时空演变预测精度产生深刻影响。针对这一问题，本研究以优化过程模型内部敏感参数有助于提高SOC预测精度为科学假设前提，选择有三期丰富实测数据的江苏省南部（苏南）为研究区，以明确表达微生物分解作用的微生物模型MIMICS为对象，以不同参数优化方法优化校正模型内部敏感参数为切入点，通过MCMC方法优化对MIMICS输出影响较大的模型内部参数，探究参数优化方法对MIMICS模型预测苏南农田表层（0~20 cm）SOC时空演变的影响，为提高微生物模型建模区域尺度SOC动态的可靠性提供理论和方法依据。
1 材料与方法 1.1 研究区概况研究区为江苏省南部地区（30°44′—32°16′N，118°29′—121°19′E），简称苏南地区（图 1）。该区地处长三角核心区，北依苏中，东临上海，南靠浙江，西接安徽，是我国经济社会发达、现代化程度较高的地区之一。该区属北亚热带季风气候，年平均气温16 ℃，年均降水量1 000 mm，气候温暖，雨量充足。地势西高东低，太湖平原区地势平坦，水系发达。土壤类型以水稻土为主，占研究区土壤总面积的83.9%；其次是黄棕壤和潮土，面积占比分别为8.1%和5.9%；其他土壤类型（如石灰（岩）土、沼泽土、红壤、紫色土等）则分布面积较小，不足3.0%。
|
图 1 研究地区地形及土壤采样点分布 Fig. 1 The topography and soil sampling sites in the study area |
MIMICS模型由美国国家大气研究中心Wieder等[18]开发，在SOC时空演变对气候变化响应等研究中得到了广泛应用[19-20]。该模型通过土壤碳周转和微生物生物量及生理机能的耦合，明确反映了SOC形成过程中，凋落物质量、微生物生理机能以及微生物副产物的物理化学保护机制三者之间的关系。MIMICS模型基于多分室建模思想，将SOC库划分为3个亚库，包括受物理保护的SOC亚库（SOCp）、可利用的SOC亚库（SOCa）和受化学保护的SOC亚库（SOCc）。该模型根据凋落物的质量，即木质素含量与氮含量比值（lignin/N）的线性函数（fmet，式（1）），将植被输入的碳分配到代谢性碳库（LITm）和结构性碳库（LITs）。根据生态学中的生物进化对策，MIMICS模型将土壤微生物分为两种功能类型：基于富养生长策略的微生物（r-策略，MICr）和基于贫养生长策略的微生物（k-策略，MICk），其中，MICr具有较高的增长率和周转率，倾向于消耗更多的LITm，而MICk则具有相对较低的增长率和周转率，在消耗LITs和SOCc时更具竞争力[19-21]。
| fmet=0.85−0.013×(lignin/N) | (1) |
| dCsdt=Is−MIC×Vmax×CsKm+Cs | (2) |
MIMICS模型以小时为模拟时间步长，采用具有温度敏感性的Michaelis-Menten动力学方程计算土壤微生物对LIT和SOC的分解速率（式（2））。式（2）中
（1）土壤样点数据。本研究收集了研究区1980年，2000年和2015年三期农田表层（0~20 cm）土壤样点数据（图 1）。其中，1980年土壤样点数量为399个，SOC含量、pH、机械组成等土壤属性数据来自研究区第二次全国土壤普查期间各县级土壤普查办公室编制的县级土种志；2000年土壤样点数量为413个，来自973项目“土壤质量演变规律与可持续利用”，主要根据二次土壤普查的样点位置进行布设；2015年土壤样点数量为407个，来自本研究的野外调查采样，主要根据2000年的土壤样点位置进行布设。2000年、2015年SOC、pH及机械组成的测定方法与1980年相同，即，SOC含量采用重铬酸钾氧化-外加热法测定，pH采用电位法（土：水=1︰2.5）测定，土壤机械组成采用吸管法测定[22]。土壤样点的SOC密度（kg∙m–2）采用以下公式计算：
| SOC密度= OC × BD × Depth × (1−δ)×100 | (3) |
式中，OC表示SOC含量（g∙kg–1），δ表示大于2 mm砾石含量（体积百分比含量）；Depth表示土层厚度（cm），本研究中为20 cm；BD表示土壤容重（g∙cm–3），采用Song等[23]建立的农田土壤容重土壤转换函数（pedotransfer function，PTF）进行估算。
（2）MIMICS模型输入数据。本研究以500 m × 500 m栅格为MIMICS模型的基本模拟单元。MIMICS模型输入数据主要包括土壤温度（Soil temperature，ST）、黏粒含量和土壤碳输入。其中，土壤温度数据来源于欧洲中期天气预报中心（European Centre For Medium-Range Weather Forecasts，ECMWF，https://apps.ecmwf.int/datasets/）[24]，空间分辨率为0.125°；黏粒含量来自本研究土壤样点黏粒含量数据；土壤碳输入基于遥感反演NPP（Net Primary Production）进行估算，来源于MUSES（Multi-resolution Satellite remote Sensing）植被生产力数据[25]，空间分辨率为5 km。
（3）SOC密度基线及敏感参数空间预测辅助变量。1980年SOC密度基线空间分布及MIMICS模型敏感参数空间分布预测的辅助变量主要包括：气候：年平均气温（Mean annual temperature，MAT）和年平均降水量（Mean annual precipitation，MAP），来源于CHELSA（Climatologies at high resolution for the earth’s land surface areas，CHELSA，https://chelsa-climate.org/）[26]，空间分辨率为1 km；土壤温度和土壤湿度（Soil moisture，SM）数据来源于ECMWF，空间分辨率为0.125°。
地形：高程（Elevation）、坡度（Slope）、坡向（Aspect）、地形位置指数（Topographic Position Index，TPI）、地形崎岖指数（Terrain Ruggedness Index，TRI）和地表粗糙度（Roughness）等地形因子，基于江苏省1︰10 000地形图生成的30 m分辨率DEM计算。
土壤类型及属性：土壤类型来源于江苏省1︰20万土壤图；土壤黏粒含量和pH来自本研究土壤样点数据。
植被和土地利用：归一化植被指数（Normalized Difference Vegetation Index，NDVI）和土地利用数据来源于中国科学院资源环境科学与数据中心（https://www.resdc.cn/），空间分辨率为1 km，其中，NDVI为1998─2005年的均值。
气候、地形及NDVI等连续型数据采用双线性内插重采样至500 m× 500 m栅格，土壤黏粒含量和pH分别采用块克里格方法插值到500 m × 500 m栅格。土壤类型及土地利用类型则根据面积优势原则分别聚合到500 m × 500 m栅格，其中，土壤类型和土地利用类型分别转换为0-1二值化变量进行哑变量编码，作为预测辅助变量。
1.4 敏感性分析方法MIMICS模型中，直接影响SOC分解过程的模型参数共30个。考虑到并非所有的参数值变动都会对研究区SOC模拟结果产生显著影响，本研究中采用Sobol全局敏感性分析方法对模型内部影响SOC过程的敏感参数进行识别。其中，MIMICS模型参数的先验不确定性区间采用Shi等[27]的方法进行设置，即，将U[a/3，3a]设置为参数的先验分布（其中，U表示均匀分布；a表示参数的默认值）。Sobol敏感性分析中，基于模型参数先验分布10 000次随机抽样计算的总敏感性指数（Total sensitivity indices，TSI）作为MIMICS模型敏感参数识别的依据。
1.5 参数优化方法基于研究区1980年样点的SOC密度观测数据，本研究采用MCMC方法对MIMICS模型敏感参数的后验分布进行推断。根据贝叶斯理论，给定SOC密度观测数据集D的条件下，MIMICS模型敏感参数集θ后验概率密度p（θ|D）的贝叶斯推断过程可表示为：
| \begin{array}{l} p\left( {\theta |D} \right) = \frac{{p\left( {D|{\rm{\mathsf{θ} }}} \right)p\left( \theta \right)}}{{p\left( D \right)}} = \hfill \\ \frac{{p\left( {D|{\rm{\mathsf{θ} }}} \right)p\left( \theta \right)}}{{\int {p\left( \theta \right)p\left( {D|{\rm{\mathsf{θ} }}} \right){\text{d}}\theta } }} \propto p\left( {D|\theta } \right)p\left( \theta \right) \hfill \\ \end{array} | (4) |
式中，p（θ）为模型参数的先验概率密度，p（D | θ）为给定参数θ条件下模型预测结果数据的似然性。式（4）表示模型参数的后验概率密度依比例收敛于模型参数的先验概率密度与数据似然性概率密度的乘积。
本研究采用自适应M-H（Metropolis-Hastings）算法实现MCMC。其中，模型预测的失配函数（misfit function）J（θ），以及建议参数值θnew被接受的概率p（θnew）分别采用以下公式计算：
| J\left( \theta \right) = \sum\nolimits_{i = 1}^n {\frac{{{{\left( {{D_i} - {M_i}} \right)}^2}}}{{2\delta _i^2}}} | (5) |
| p\left( {{\theta _{{\text{new}}}}} \right) = min\left\{ {1,exp\left[ { - \left( {J\left( {{\theta _{{\text{new}}}}} \right) - J\left( {{\theta _{{\text{old}}}}} \right)} \right)} \right]} \right\} | (6) |
式中，
本研究通过2个步骤实现模型参数后验分布的抽样[28]。
第1步：基于MCMC估算模型参数集的初始协方差矩阵C0。此步骤中，假设模型参数θ的先验分布p（θ）服从最小值为θmin、最大值为θmax的均匀分布U（θmin，θmax），采用以下公式产生建议的参数值θnew：
| {\theta _{{\text{new}}}} = {\theta _{{\text{old}}}} + r\frac{{{\theta _{\max }} - {\theta _{\min }}}}{S} | (7) |
式中，θold表示前一次迭代中接受的参数值；r为[–0.5，0.5]范围内的均匀分布随机数；S为步长控制系数，本研究中取值为5，即，新参数值的建议过程中，参数值变化步长为相应参数极差的1/10。迭代运行m次，基于第m/2到第m次迭代中被接受的建议参数值来估算模型参数集θ的初始协方差矩阵C0。
第2步：基于MCMC推断模型参数后验分布。此步骤中，假设模型参数服从多元正态分布。这里新产生的建议参数值依赖于k-1次接受的参数值
| {C_t} = \left\{ {\begin{array}{*{20}{l}} {{C_0}} \\ {sdCov({\theta _0}, \cdots \cdots ,{\theta _{t - 1}})} \end{array}} \right. | (8) |
式中的协方差矩阵C0是基于第1步中所接受的参数集估算产生，sd为常数，其值等于2.38sqrt（d），d为待优化参数总个数，本研究取值为9。本研究采用Gelman-Rubin（G-R）统计量来检验马尔科夫链是否收敛[17，27]。如果马尔科夫链收敛，链间方差与单条链方差接近，则G-R统计量的值接近于1。
本研究中1980年SOC密度观测数据用于MIMICS模型的参数优化，共设计了批处理（One Batch）和点对点（Site-by-Site）两种参数优化方案。批处理参数优化方案中，将1980年所有SOC密度观测值看作一个整体来共同约束参数值（即，不考虑模型参数的空间差异性），优化过程中采用3条马尔科夫链并行，第1步每条链迭代10 000次，第2步每条链迭代100 000次。第2步生成的马尔科夫链丢弃前50%样本后再计算每个参数的G-R统计量，以及每个参数的后验分布。以3条链后验分布中值的平均值作为模型参数的最优值，驱动MIMICS模型的模拟。该方案是针对整个研究区优化出一套最优参数集。
点对点参数优化方案中，对于1980年的每个采样点位置，分别采用其SOC密度观测值对样点位置的参数值进行约束（即，不同空间位置的模型参数值不同）。该方案中，对于1980年的每个采样点位置，分别采用3条马尔科夫链并行进行参数优化，其中，第1步每条链迭代10 000次，第2步每条链迭代50 000次。将第2步生成的马尔科夫链丢弃前50%样本后再计算每个采样点位置上每个参数的G-R统计量，以及每个参数的后验分布。以每个样点位置MIMICS模型参数的最优值（3条链后验分布中值的平均值）为目标变量，以MAP、MAT、ST、SM、NDVI、NPP、土壤黏粒、pH、地形因子、土壤类型和土地利用等22个协同变量为预测变量，采用随机森林（Random Forest，RF）预测模型敏感参数最优值的空间分布，用于驱动MIMICS模型的模拟。
1.6 模型精度评价本研究基于2000年和2015年两期SOC密度观测数据，对MIMICS模型模拟精度进行独立验证。采用决定系数（R2）、均方根误差（Root Mean Square Error，RMSE）和平均绝对误差（Mean Absolute Error，MAE）3个指标进行精度评价。
本研究中，MIMICS模型和MCMC算法采用Fortran语言编程实现，通过R语言调用MIMICS模型和MCMC算法，实现SOC时空动态模拟和参数优化。模型参数敏感性分析采用R语言中的sensitivity包[10]实现，敏感参数最优值空间分布的RF预测采用R语言的caret包实现[29]。
2 结果 2.1 SOC密度的描述性统计特征SOC密度的描述性统计分析结果表明（表 1），1980年、2000年和2015年苏南农田表层SOC密度（以下分别用SOC1980、SOC2000和SOC2015表示）中值分别为3.36 kg∙m–2、4.25 kg∙m–2和3.81 kg∙m–2。SOC密度变化表现为先增加后降低的趋势，其中，1980─2000年SOC密度净增幅为0.89 kg∙m–2，而2000─2015年净降幅为0.44 kg∙m–2。1980─2015年的35年间，研究区SOC密度净增幅为0.45 kg∙m–2，表层土壤年均固碳速率为129 kg∙hm–2，低于华东地区农田表层土壤1980─2011年间的年均固碳速率268 kg∙hm–2[30]。SOC1980、C2000和SOC2015的偏度系数分别为0.55、0.19和-0.09，偏度系数均 < |1|，表明极端值对3个年份SOC密度变异性的影响均不大。然而，SOC2015的变异系数（34%）和极差（7.15 kg∙m–2）均显著高于SOC1980（26%和4.82 kg∙m–2）和SOC2000（23%和6.70 kg∙m–2），表明相较前两个采样年份而言，SOC2015的变异性更强，因此，其可能受人为活动的影响更为强烈。
|
表 1 表层（0~20cm）土壤有机碳密度的描述性统计 Table 1 Descriptive statistics of topsoil(0~20 cm)SOC density(kg∙m–2)in different sampling date |
敏感性分析的目的在于量化模型参数变动对模型输出结果变动的影响程度，从而从众多的模型参数中识别出对模拟结果影响较大的敏感参数。过程模型敏感参数的识别不但有利于区域模拟中模型参数化过程的简化，而且可为模型参数优化提供更为明晰的优化对象。从基于Sobol方法的MIMICS模型参数全局敏感性分析结果来看（图 2），最敏感的参数为Vint，其TSI为0.85，随后依次分别为kint（TSI = 0.76）、Vslop（TSI = 0.37）、ak（TSI = 0.31）、aV（TSI = 0.16）和CUE_m_r（TSI = 0.06）等。参数fi_met的敏感性最弱，TSI仅为0.01，表明在本研究区的气候、土壤和C输入条件下，fi_met的变动对MIMICS模拟的SOC密度结果影响极小。此外，fi_met、KO_k及KO_r的敏感性亦较弱。本研究采用0.035作为敏感性指数的阈值（一般≤0.05[31]），以区分敏感和不敏感参数。共识别到TSI > 0.035的敏感性参数9个，包括与米氏方程中最大分解速率Vmax计算相关的Vslop、Vint和aV共3个，与半饱和浓度Km计算相关的Kslope_a、Kint、ak、Kmod_a_r_coeff共4个，以及CUE_m_k和CUE_m_r。
|
图 2 MIMICS模型参数全局敏感性分析指标 Fig. 2 Ranked MIMICS parameters based on the total sensitivity indices |
采用MCMC算法对识别的9个敏感参数进行优化，从模型参数的马尔科夫链收敛性诊断结果来看（表 2），批处理优化和点对点优化方法中，基于3条并行马尔科夫链计算的G-R统计量介于1.00~ 1.06之间，表明两种参数优化方法的马尔科夫链均收敛于平稳分布[32]。从模型参数优化的后验分布中值统计结果来看，Vslope、Vint、av、ak、Kint、Kslope_a和Kmod_a_r_coeff均大于模型参数默认值，而与CUE估算密切相关的CUE_m_k和CUE_m_r则显著低于模型的默认值（表 2）。批处理方法和点对点方法优化的Vslope、Vint、ak、Kint和Kmod_a_r_coeff的后验分布中值分别为0.13和0.11、12.48和12.73、18.22和15.73、5.95和4.01、0.53和0.42（表 2），其中Vint、aK和Kint与Shi等[27]采用批处理方法在全球尺度上的优化结果相比（Vslope = 0.14、Vint = 8.09、ak = 16.97、Kint = 5.38、Kmod_a_r_coeff = 0.42）具有明显差异，表明MIMICS模型在区域尺度应用时，对模型参数进行优化校正是十分必要的。此外，批处理方法进行参数优化时，CUE_m_k的后验分布中值为0.03，而点对点方法优化的CUE_m_k后验分布中值为0.27，为批处理方法优化结果的9倍，说明CUE_m_k的空间异质性对研究区表层SOC密度动态建模的影响极大。因此，如果整个研究区均采用批处理方法优化的单一参数值，则可能导致模拟结果误差增大。
|
|
表 2 马尔科夫链的收敛性诊断及MIMICS模型参数的后验分布中值统计 Table 2 The Gelman-Rubin convergence diagnostic and the median of model parameter posterior distribution |
以每个样点位置MIMICS模型参数的最优值为目标变量，以气候、地形、土壤类型和土地利用等22个协同变量为预测变量，采用RF预测模型敏感参数最优值的空间分布，用于驱动MIMICS模型的模拟。RF模型十折交叉验证结果表明，不同敏感参数空间预测的R2介于0.77~0.89之间，敏感参数空间分布预测的精度较高。
2.3 不同参数优化方法下MIMICS模型预测的SOC时空演变特征不同参数优化方法下MIMICS模型建模苏南农田表层SOC密度动态变化结果表明（图 3），默认参数与优化参数预测的SOC密度动态差异较大，但批处理和点对点两种参数优化方法下，SOC密度动态变化规律较为相似。具体而言，默认参数下，MIMICS模型预测的苏南农田表层SOC密度从1980年的3.39 kg∙m–2增加到2015年的4.33 kg∙m–2，SOC变化表现为逐年持续增加趋势。其中，2000年SOC密度预测值为4.10 kg∙m–2，较2000年观测值（4.25 kg∙m–2）低3.6%，而2015年SOC密度的预测值为4.33 kg∙m–2，较2015年观测值（3.81 kg∙m–2）高13.7%。批处理优化敏感参数的条件下，SOC密度先从1980年的3.39 kg∙m–2增加至2002年的4.26 kg∙m–2，再减少至2015年的3.68 kg∙m–2，总体呈先增加后降低的趋势。其中，2000年和2015年SOC密度的预测值分别为4.22 kg∙m–2和3.68 kg∙m–2，较观测值（4.25 kg∙m–2和3.81 kg∙m–2）分别低0.7%和3.3%。点对点优化模型敏感参数的条件下，研究区SOC密度先从1980年的3.39 kg∙m–2增加至1999年的4.18 kg∙m–2，再减少至2015年的3.87 kg∙m–2，亦呈现先增加后降低的趋势。其中，2000年SOC密度预测值为4.17 kg∙m–2，较观测值低1.7%，而2015年SOC密度预测值为3.87 kg∙m，–2较观测值高1.7%。
|
图 3 1980─2015年间苏南农田表层SOC密度变化趋势 Fig. 3 The dynamics of soil organic carbon density between 1980 and 2015 derived from default parameters, one batch, and site-by-site |
从农田表层SOC密度变化速率的空间分布来看（图5），不同参数优化方法下其空间分布模式具有明显的差异。MIMICS模型默认参数条件下，1980─2015年间，研究区西部、中部及北部沿江地区SOC密度显著增加，且大部分地区SOC增加速率高于0.02 kg∙m–2∙a–1；研究区东南部SOC密度降低，丢碳速率超过0.01 kg∙m–2∙a–1（图 4a）。批处理方法优化MIMICS模型敏感参数条件下，1980─2015年间苏南农田SOC密度总体呈增加趋势，其中，研究区西部SOC密度显著增加（图 4b），但增加速率明显低于默认参数的估计结果（图 4a），研究区中部SOC密度略有增加，但增速较低，而研究区东部增速更低，且有大面积的SOC密度降低区域（图 4b）。而点对点参数优化方法下，SOC密度变化速率的空间分异模式更为明显，其中，西部SOC密度增加的区域范围明显高于批处理参数优化方法的预测结果，而东部SOC密度降低的区域范围则较批处理参数优化方法的预测结果稍小（图 4b、图 4c）。
|
图 4 默认参数（a），批处理（b）及点对点（c）方法下1980─2015年农田表层SOC密度变化速率空间分布 Fig. 4 Spatial distributions of the rate of change in topsoil SOC density between 1980 and 2015 derived from the default parameters(a), one-batch optimized parameters(b), and site-by-site optimized parameters(c) |
从MIMICS模型预测精度的独立验证结果来看（表 3），采用模型默认参数值时的模拟误差最大，其中，2000年和2015年表层SOC密度预测RMSE分别为0.99 kg∙m–2和1.43 kg ∙m–2，仅能解释2000年、2015年表层SOC密度观测值变异性的31%和3%。当采用批处理方法优化MIMICS模型敏感参数后，2000年和2015年表层SOC密度预测的RMSE分别降低至0.81 kg∙m–2和1.25 kg∙m–2，较采用默认参数值的预测误差分别降低了18.2%和12.6%。同时，批处理优化敏感参数方法下，MIMICS模型解释的2000年、2015年表层SOC密度观测值变异性分别提高到37%和7%。当采用点对点优化MIMICS模型敏感参数，并利用RF模型估算的最优参数值空间分布数据驱动MIMICS模型时，2000年和2015年表层SOC密度预测的RMSE分别进一步降低至0.77 kg∙m–2和1.22 kg∙m–2，与采用默认参数值模拟相比，RMSE分别降低了22.2%（2000年）和14.7%（2015年），MIMICS模型解释的表层SOC密度观测值变异性则分别提高到43%（2000年）和13%（2015年），表明优化模型敏感参数并表征其空间异质性能有效提高MIMICS模型预测表层SOC密度时空演变的精度。
|
|
表 3 不同参数优化方法下MIMICS模型预测精度 Table 3 Prediction accuracy of MIMICS model with different parameter optimization methods |
土壤及气候条件等显著异于过程模型开发环境时，直接采用模型参数默认值来建模SOC时空演变可能导致模拟结果产生极大的不确定性[11，33]。因此，基于观测数据对模型参数进行本地化校正，是降低过程模型模拟SOC时空演变不确定性的重要手段。从本研究的结果来看，采用默认参数值驱动MIMICS模型，仅能够再现研究区农田前20年（1980─2000年）的表层SOC密度增加趋势，但不能反映随后15年（2000─2015年）的表层SOC密度降低趋势。而采用批处理和点对点两种方法优化敏感参数后，MIMICS模型能够反映与观测数据较为一致的SOC密度变化趋势，即1980─2015年间研究区农田表层SOC密度先增加后减少（图 3，表 1）。
此外，采用批处理方法优化MIMICS模型敏感参数后，2000年和2015年表层SOC密度预测RMSE分别较采用默认参数值方法降低18.2%和12.6%，而MIMICS模型解释的表层SOC密度变异性分别提高到37%和7%，表明采用批处理优化模型敏感参数有利于提高模型的时空预测能力，并显著降低由于模型参数偏差所导致的预测结果不确定性，这与Hararuk等[11]的研究结果一致。微生物模型内部参数大多为常值参数，但常值参数可能难以表征区域、国家、全球等大尺度上土壤碳循环动态的显著空间异质性[16-17]，从而对SOC时空演变预测精度产生深刻影响。比如，Ge等[34]研究发现，DALEC（Data Assimilation Linked Ecosystem Carbon，DALEC）模型中表征碳库周转过程的参数因气候带不同而具有明显差异。本研究结果表明，当采用点对点方法优化MIMICS模型敏感参数，并利用RF模型估算最优参数值的空间分布，进而利用其驱动MIMICS模型时，可使研究区2000年和2015年表层SOC密度预测RMSE分别降低22.2%和14.7%，其预测精度较采用批处理方法优化参数的预测精度更高（表 3）。由此可见，除了对模型敏感参数进行优化外，充分考虑模型参数优化值的空间异质性也是进一步降低MIMICS模型预测SOC时空演变预测误差的重要途径。
3.2 微生物模型模拟SOC时空动态的局限性尽管采用不同的参数优化方法能够在一定程度上降低研究区农田表层SOC密度时空动态的预测误差，然而，微生物模型建模区域尺度农田SOC密度时空演变依然具有一定的局限性。
首先，尽管点对点优化参数方法下MIMICS模型的预测精度最高，但2015年SOC密度的预测RMSE仍然较2000年SOC密度的预测RMSE高58%。表明，参数优化的MIMICS模型对强烈人为活动引起的高度变异SOC模拟精度依然偏低。这可能与MIMICS模型不能完全反映人为活动对SOC的影响过程有关。比如，Xie等[35]研究表明，人为活动导致1980─2000年苏南农田土壤平均酸化速率为每年0.028个pH单位，而2000─2015年则高达每年0.035个pH单位。土壤酸化过程影响土壤微生物活性，进而影响SOC分解速率，但MIMICS模型未考虑土壤pH变化的影响。同时，Xie等[36]研究发现，秸秆还田、化肥施用、以及土地利用变化是影响苏南农田表层SOC变化的主要人为因子。秸秆还田与合理施用化肥有利于增加农田土壤碳输入，促进农田土壤固碳[15，30]，而土地利用变化则通过改变土壤理化属性从而影响作物生长发育及土壤微生物活性等，进而影响土壤碳输入及SOC分解[37-38]。1980─ 2000年间苏南耕地面积减少了5.6%，城镇建设和农村居民点占地面积则增加了4.5%，而2000─2015年间苏南耕地面积减少了9.1%，城镇建设和农村居民点占地面积则增加了8.7%，然而，目前的微生物模型仅能以碳输入变化来反映土地利用变化的间接影响，还难以反映土地利用变化对SOC累积和分解过程的直接影响，这可能是导致2015年SOC预测精度偏低的主要原因之一。
其次，苏南不同区域1980年的初始SOC含量、施肥量及秸秆还田量、土地利用变化均存在一定的区域差异性[36]，因此，不同区域的SOC变化趋势亦存在明显的差异性。1980年研究区东部（江阴、锡山、惠山、滨湖、梁溪、新吴、张家港、常熟、太仓、昆山、虎丘、吴中、相城、姑苏和吴江）和西部（丹徒、丹阳、金坛、溧阳、武进、天宁、钟楼、新北和宜兴）的农田表层平均SOC密度分别为3.60 kg∙m–2和3.01 kg∙m–2，尽管东部地区SOC密度较西部稍高，但仍处于相对较低水平。因此，土壤碳输入的增加有助于促进SOC的累积，例如，2000年时，研究区东部农田表层平均SOC密度增加到4.55 kg∙m–2，而研究区西部则增加到3.96 kg∙m–2。2000年后研究区秸秆还田政策大面积推广[30]，与此同时，土地利用也发生了更为剧烈的变化，主要表现为农田面积加速减少、城镇建设用地面积持续增加，尤其以研究区东部更为明显[39]。然而，2000年之后土地利用剧烈变化导致的土壤碳输入量停滞不前，已不足以维持研究区较高的SOC水平，从而导致SOC降低。例如，2000─2015年，研究区东部SOC密度平均降幅为0.72 kg∙m–2，而西部地区SOC密度降幅稍小，为0.18 kg∙ m–2。尽管点对点参数优化方法通过MIMICS模型敏感参数的空间异质性表达，从而提供了关于研究区农田表层SOC动态区域差异性的更多局部细节信息，但MIMICS模型还难以全面表达土壤碳输入、农业管理措施及土地利用变化交互作用对SOC动态的影响过程，因此，进一步改进微生物模型的结构，将是未来提高微生物模型预测SOC时空动态精度所面临的重要挑战。
此外，MIMICS微生物模型采用具有温度敏感性的米氏动力学方程描述土壤碳的分解，这对于探讨SOC对于温度变化的响应非常有益，但目前在区域尺度上可免费获取的长时间序列土壤温度数据的空间分辨率仍然相对较低。比如，本研究中土壤温度数据空间分辨率为0.125°，在反映研究区土壤温度的局部空间变异性方面还具有一定的局限性，这也可能会对本研究中的SOC时空演变模拟结果产生一定的影响。
4 结论默认参数、批处理以及点对点优化过程模型参数对MIMICS微生物模型建模苏南地区农田表层SOC时空动态的对比分析表明：（1）采用默认参数值的MIMICS模型不能完全反映观测到的SOC变化总体趋势，但敏感参数优化后的MIMICS模型能够再现苏南农田表层SOC密度先增加后减少的总体变化趋势；（2）与采用默认参数值相比，采用批处理方法优化MIMICS模型敏感参数后，2000年和2015年SOC密度预测RMSE分别降低18.2%和12.6%，进一步考虑模型敏感参数的空间异质性，采用点对点方法优化MIMICS模型敏感参数后，2000年和2015年SOC密度预测RMSE则较默认参数法分别降低22.2%和14.7%，其预测精度最高；（3）尽管点对点参数优化方法下MIMICS模型的预测精度最高，但其2015年SOC密度的预测精度依然偏低，因此，进一步改进微生物模型的结构、提高模型输入数据的精度及分辨率，将是未来提高微生物模型预测SOC时空动态精度所面临的重要挑战。
| [1] |
Lal R. Societal value of soil carbon[J]. Journal of Soil and Water Conservation, 2014, 69(6): 186A-192A. DOI:10.2489/jswc.69.6.186A
(  1) |
| [2] |
Lal R. Soil carbon sequestration impacts on global climate change and food security[J]. Science, 2004, 304(5677): 1623-1627. DOI:10.1126/science.1097396
( 0) |
| [3] |
Lal R. Soil health and carbon management[J]. Food and Energy Security, 2016, 5(4): 212-222. DOI:10.1002/fes3.96
( 1) |
| [4] |
Vereecken H, Schnepf A, Hopmans J W, et al. Modeling soil processes: Review, key challenges, and new perspectives[J]. Vadose Zone Journal, 2016, 15(5): 1-57.
( 1) |
| [5] |
Jenkinson D S. The turnover of organic carbon and nitrogen in soil[J]. Philosophical Transactions of the Royal, 1990, B329: 361-368.
( 1) |
| [6] |
Li C S, Zhuang Y H, Frolking S, et al. Modeling soil organic carbon change in croplands of China[J]. Ecological Applications, 2003, 13(2): 327-336. DOI:10.1890/1051-0761(2003)013[0327:MSOCCI]2.0.CO;2
( 1) |
| [7] |
Luo Y Q, Wu L H, Andrews J A, et al. Elevated CO2 differentiates ecosystem carbon processes: Deconvolution analysis of duke forest face data[J]. Ecological Monographs, 2001, 71(3): 357-376.
( 1) |
| [8] |
Bonan G B, Doney S C. Climate, ecosystems, and planetary futures: The challenge to predict life in Earth system models[J]. Science, 2018, 359(6375): eaam8328. DOI:10.1126/science.aam8328
( 1) |
| [9] |
Luo Y Q, Ahlström A, Allison S D, et al. Toward more realistic projections of soil carbon dynamics by Earth system models[J]. Global Biogeochemical Cycles, 2016, 30(1): 40-56. DOI:10.1002/2015GB005239
( 2) |
| [10] |
Gurung R B, Ogle S M, Breidt F J, et al. Bayesian calibration of the DayCent ecosystem model to simulate soil organic carbon dynamics and reduce model uncertainty[J]. Geoderma, 2020, 376: 114529. DOI:10.1016/j.geoderma.2020.114529
( 2) |
| [11] |
Hararuk O, Xia J Y, Luo Y Q. Evaluation and improvement of a global land model against soil carbon data using a Bayesian Markov chain Monte Carlo method[J]. Journal of Geophysical Research: Biogeosciences, 2014, 119(3): 403-417. DOI:10.1002/2013JG002535
( 4) |
| [12] |
Xu H Y, Zhang T, Luo Y Q, et al. Parameter calibration in global land carbon models using surrogate-based optimization[J]. Geoscientific Model Development Discussions, 2017, 1-28.
( 1) |
| [13] |
Domeignoz-Horta L A, Pold G, Liu X J A, et al. Microbial diversity drives carbon use efficiency in a model soil[J]. Nature Communications, 2020, 11: 3684. DOI:10.1038/s41467-020-17502-z
( 1) |
| [14] |
Wang B R, An S S, Liang C, et al. Microbial necromass as the source of soil organic carbon in global ecosystems[J]. Soil Biology and Biochemistry, 2021, 162: 108422. DOI:10.1016/j.soilbio.2021.108422
( 1) |
| [15] |
Zhang X, Zhao Y C, Xie E Z, et al. Spatio-temporal change of soil organic carbon, progress and prospects (In Chinese)[J]. Journal of Agro-Environment Science, 2020, 39(4): 673-679. [张秀, 赵永存, 谢恩泽, 等. 土壤有机碳时空变化研究进展与展望[J]. 农业环境科学学报, 2020, 39(4): 673-679.]
( 3) |
| [16] |
Luo Y Q, White L W, Canadell J G, et al. Sustainability of terrestrial carbon sequestration: A case study in Duke Forest with inversion approach[J]. Global Biogeochemical Cycles, 2003, 17(1): 1021.
( 2) |
| [17] |
Tao F, Zhou Z H, Huang Y Y, et al. Deep learning optimizes data-driven representation of soil organic carbon in earth system model over the conterminous United States[J]. Frontiers in Big Data, 2020, 3: 17. DOI:10.3389/fdata.2020.00017
( 5) |
| [18] |
Wieder W R, Grandy A S, Kallenbach C M, et al. Representing life in the earth system with soil microbial functional traits in the mimics model[J]. Geoscientific Model Development, 2015, 8(6): 1789-1808. DOI:10.5194/gmd-8-1789-2015
( 1) |
| [19] |
Wieder W R, Grandy A S, Kallenbach C M, et al. Integrating microbial physiology and physio-chemical principles in soils with the microbial-mineral carbon stabilization(mimics)model[J]. Biogeosciences, 2014, 11(14): 3899-3917. DOI:10.5194/bg-11-3899-2014
( 2) |
| [20] |
Wieder W R, Hartman M D, Sulman B N, et al. Carbon cycle confidence and uncertainty: Exploring variation among soil biogeochemical models[J]. Global Change Biology, 2018, 24(4): 1563-1579. DOI:10.1111/gcb.13979
( 1) |
| [21] |
Zhang H C, Goll D S, Wang Y P, et al. Microbial dynamics and soil physicochemical properties explain large-scale variations in soil organic carbon[J]. Global Change Biology, 2020, 26(4): 2668-2685. DOI:10.1111/gcb.14994
( 1) |
| [22] |
Xie E Z, Zhao Y C, Lu F Y, et al. Comparison analysis of methods for prediction of spatial distribution of soil organic matter contents in farmlands south Jiangsu, China (In Chinese)[J]. Acta Pedologica Sinica, 2018, 55(5): 1051-1061. [谢恩泽, 赵永存, 陆访仪, 等. 不同方法预测苏南农田土壤有机质空间分布对比研究[J]. 土壤学报, 2018, 55(5): 1051-1061.]
( 1) |
| [23] |
Song G H, Li L Q, Pan G X, et al. Topsoil organic carbon storage of China and its loss by cultivation[J]. Biogeochemistry, 2005, 74(1): 47-62. DOI:10.1007/s10533-004-2222-3
( 1) |
| [24] |
Su Z, Rosnay P de, Wen J, et al. Evaluation of ECMWF's soil moisture analyses using observations on the Tibetan Plateau[J]. Journal of Geophysical Research: Atmospheres, 2013, 118(11): 5304-5318. DOI:10.1002/jgrd.50468
( 1) |
| [25] |
Wang J M, Sun R, Zhang H L, et al. New global MuSyQ GPP/NPP remote sensing products from 1981 to 2018[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2021, 14: 5596-5612. DOI:10.1109/JSTARS.2021.3076075
( 1) |
| [26] |
Karger D N, Conrad O, Böhner J, et al. Climatologies at high resolution for the earth's land surface areas[J]. Scientific Data, 2017, 4: 170122. DOI:10.1038/sdata.2017.122
( 1) |
| [27] |
Shi Z, Crowell S, Luo Y Q, et al. Model structures amplify uncertainty in predicted soil carbon responses to climate change[J]. Nature Communications, 2018, 9: 2171. DOI:10.1038/s41467-018-04526-9
( 3) |
| [28] |
Xu T, White L, Hui D F, et al. Probabilistic inversion of a terrestrial ecosystem model: Analysis of uncertainty in parameter estimation and model prediction[J]. Global Biogeochemical Cycles, 2006, 20(2): GB2007.
( 1) |
| [29] |
Taghizadeh-Mehrjardi R, Schmidt K, Amirian-Chakan A, et al. Improving the spatial prediction of soil organic carbon content in two contrasting climatic regions by stacking machine learning models and rescanning covariate space[J]. Remote Sensing, 2020, 12: 1-26.
( 1) |
| [30] |
Zhao Y C, Wang M Y, Hu S J, et al. Economics- and policy-driven organic carbon input enhancement dominates soil organic carbon accumulation in Chinese croplands[J]. Proceedings of the National Academy of Sciences of the United States of America, 2018, 115(16): 4045-4050.
( 3) |
| [31] |
Liu X Y, Zhao Y C, Shi X Z, et al. Sensitivity and uncertainty analysis of CENTURY-modeled SOC dynamics in upland soils under different climate-soil-management conditions: A case study in China[J]. Journal of Soils and Sediments, 2017, 17(1): 85-96. DOI:10.1007/s11368-016-1516-0
( 1) |
| [32] |
Brooks S P, Gelman A. General methods for monitoring convergence of iterative simulations[J]. Journal of Computational and Graphical Statistics, 1998, 7(4): 434-455.
( 1) |
| [33] |
Bloom A A, Exbrayat J F, van der Velde I R, et al. The decadal state of the terrestrial carbon cycle: Global retrievals of terrestrial carbon allocation, pools, and residence times[J]. Proceedings of the National Academy of Sciences of the United States of America, 2016, 113(5): 1285-1290.
( 1) |
| [34] |
Ge R, He H L, Ren X L, et al. Underestimated ecosystem carbon turnover time and sequestration under the steady state assumption: A perspective from long-term data assimilation[J]. Global Change Biology, 2019, 25(3): 938-953. DOI:10.1111/gcb.14547
( 1) |
| [35] |
Xie E Z, Zhao Y C, Li H D, et al. Spatio-temporal changes of cropland soil pH in a rapidly industrializing region in the Yangtze River Delta of China, 1980—2015[J]. Agriculture, Ecosystems & Environment, 2019, 272: 95-104.
( 1) |
| [36] |
Xie E Z, Zhang X, Lu F Y, et al. Spatiotemporal changes in cropland soil organic carbon in a rapidly urbanizing area of southeastern China from 1980 to 2015[J]. Land Degradation & Development, 2022, 33(9): 1323-1336.
( 2) |
| [37] |
Fan M M, Lal R, Zhang H, et al. Variability and determinants of soil organic matter under different land uses and soil types in Eastern China[J]. Soil and Tillage Research, 2020, 198: 104544. DOI:10.1016/j.still.2019.104544
( 1) |
| [38] |
Luo Y L, Li Q Q, Shen J, et al. Effects of agricultural land use change on organic carbon and its labile fractions in the soil profile in an urban agricultural area[J]. Land Degradation & Development, 2019, 30(15): 1875-1885.
( 1) |
| [39] |
Yirsaw E, Wu W, Temesgen H, et al. Socioeconomic drivers of spatio-temporal land use/land cover changes in a rapidly urbanizing area of China, the Su-Xi-Chang region[J]. Applied Ecology and Environmental Research, 2017, 15(4): 809-827.
( 1) |