做陆地生态系统碳循环模拟的人,手里大概都绕不开 Biome-BGC 这个名字。我最初接触它是在做森林碳汇评估的时候,当时需要估算一个省份的人工林净生态系统交换量,翻了一圈文献发现,要么是统计模型太粗糙,要么是遥感反演缺乏机制支撑,最后几乎所有同行都指向了同一个工具——Biome-BGC。这个上世纪九十年代就成型的生物地球化学过程模型,到今天依然是站点尺度到区域尺度碳氮水循环模拟的标配之一。它不花哨,没有漂亮的图形界面,甚至输入文件格式还停留在固定列宽的纯文本,但就是这种"老派"工具,在生态学、林学、气候变化响应的研究里,承担了大量"机制解释"的硬活。
这篇内容我打算把 Biome-BGC 从原理到实操完整梳理一遍。适合刚接触过程模型、想用它跑通第一个模拟案例的研究生,也适合已经跑过模型但总在数据准备和参数校准上反复踩坑的同行。我会把模型的核心逻辑拆开讲清楚,把输入数据的每个坑都标出来,最后附上我实际运行中的踩坑记录,都是文档里不会写的东西。
1. 这个模型到底在算什么账
Biome-BGC 的全称是 Biome BioGeochemical Cycles,出自 Montana 大学 Running 团队。它是一个以日为步长的过程模型,核心任务是模拟陆地生态系统的碳、氮、水三个生物地球化学循环,以及三者之间的耦合关系。和 Miami 模型、CASA 这类光能利用率模型最大的区别在于,它不是靠遥感植被指数直接推算生产力,而是从叶片尺度出发,模拟光合作用、呼吸作用、碳分配、凋落物分解、土壤有机质周转这一整套机制过程。
1.1 碳账本、水账本和氮账本
模型里的碳循环主线是这样的:每天根据气象条件、叶面积指数和叶片氮含量,用 Farquhar 生化模型算出总初级生产力 GPP,然后减去自养呼吸得到净初级生产力 NPP,再按照根系、茎、叶、粗根等不同库的分配比例把碳分到各个器官。凋落物和枯死组织进入土壤后,经过异养呼吸释放回大气,剩下的就是净生态系统交换量 NEE。
水循环的驱动逻辑是 Penman-Monteith 方程。模型通过气孔导度和叶面积指数计算冠层蒸腾,加上土壤蒸发和冠层截留蒸发,得到总蒸散量。土壤水分通过一个多层土壤水分的收支平衡来模拟,水分亏缺会反馈到气孔导度和光合作用上,形成水-碳耦合。这一点在干旱半干旱区尤其重要,我记得在黄土高原做模拟的时候,如果不考虑水分胁迫对光合的抑制,模拟出的 GPP 会偏高一倍以上。
氮循环是很多人容易忽略的一块。Biome-BGC 里氮不是陪跑,它直接限制光合作用的速率:叶片氮含量决定 Rubisco 酶活性,进而决定羧化速率上限。同时,植物生长分配需要氮,凋落物分解需要氮,土壤有机质周转受碳氮比调节。简单说,碳积累的速度被氮供应卡着脖子,这就是所谓的氮限制机制。
1.2 时间尺度与空间尺度
模型默认以日为步长运行,但气象输入可以是日值,输出也是日值。你要想跑出年尺度的净生态系统生产力,一般会连续运行几十年,再逐年汇总。空间上,Biome-BGC 最初是站点模型,后来通过地理分区参数化的方式扩展到区域尺度,也就是把研究区按植被类型和气候区划分成若干单元,每个单元独立跑一遍模拟,最后汇总。这种方式在数据条件有限时很实用,但要注意它不是真正的网格模型,没有空间交互过程,每个像元之间是独立模拟的。
运行模式上有个非常关键的概念叫"平衡态运行"。模型初始化时,土壤碳氮库需要被"预热"到接近稳定的状态,否则初值会对模拟结果产生长期影响。通常做法是先让模型空转几千上万年,直到土壤碳氮库达到动态平衡,再在此基础上接上真实气象数据做模拟。这一步往往决定整个模拟的成败,后面我会专门讲。
2. 碳氮水循环的齿轮组:模型核心机制拆解
很多人拿到模型就开始填参数,对内部机制一知半解,最后出了问题也不知道是哪里错了。我建议至少把下面几个核心模块的逻辑先吃透,调试起来会顺手很多。
2.1 光合作用:Farquhar模型怎么在代码里落地
Biome-BGC 的光合模块用的是 Farquhar 生化模型,它把叶片光合作用分解为两个可能限制步骤:一个是 Rubisco 酶活性限制的羧化速率,一个是 RuBP 再生能力限制的电子传递速率。实际光合速率取两者之间的最小值。在这里,叶片氮含量是核心参数,因为 Rubisco 酶就占叶片总氮的很大比例。
气象变量里,辐射决定了光合有效辐射 PAR 的供给,温度决定酶的活性范围,水汽压亏缺通过气孔导度影响 CO2 向内扩散。这三者的交互在模型里表现得非常敏感。我做过一个实验,把夏季高温期的 VPD 数据人工抬高 0.5 kPa,模拟出的整个生长季 GPP 就下降了 18%。所以气象数据的质量对光合模块的准确性影响极大。
2.2 呼吸与分配:维持性呼吸和生长性呼吸分开算
呼吸在模型里分两类。维持性呼吸和生物量成正比,和温度呈指数关系,Q10 通常取 2.0;生长性呼吸则按新建组织的构建成本计算,一般占 NPP 的一小部分比例。要特别注意,叶片的维持性呼吸在冬季可能大于光合积累,这时净碳交换会变为负值,这在北方森林的模拟中很常见。
碳分配系数控制光合产物往叶、茎、细根、粗根的流向。模型按物候阶段动态调整分配比例,比如落叶林在展叶期优先分配叶片,生长季后期转向茎和根。很多初学者把这些系数当作固定值填,结果模拟出的叶面积指数季节曲线完全对不上实测。实际上,分配系数应该和实测的器官生物量数据反复校对。
2.3 凋落物与土壤碳周转
植物枯死后,碳和氮进入凋落物库。模型把凋落物按结构性和代谢性两类区分,前者木质素含量高、分解慢,后者易分解。土壤有机质分成多个周转速率不同的库,分解速率受温度和水分的修正。这些库之间的碳氮流动关系用一个连锁矩阵控制,参数多,但稳定性很好。
这里有个容易忽略的点:模型的分解过程不区分土壤层次,所有土壤碳都放在一个混合库里。这意味着它没法解释深层土壤有机碳和表层有机碳的周转差异。如果你的研究恰好关注土壤剖面分布,Biome-BGC 就不合适,得换 Century 或者使用像 CLM 这类更复杂的模型。
2.4 氮循环的闸门效应
氮进入生态系统的主要途径是矿化、大气沉降和生物固氮,输出途径是淋溶和气体损失。模型里矿物氮库的丰缺直接反馈到植物吸收和微生物分解上。高碳氮比的凋落物在分解初期会固定土壤矿物氮,造成植物氮供应不足,这个机制在模型里被忠实复现了。正因为有这层耦合,模拟长期碳汇时才会出现"氮限制饱和"现象——大气 CO2 升高刺激的光合增益,会逐渐被氮供应不足抵消。这也是为什么很多用 Biome-BGC 跑未来情景的研究,得到的碳汇增长远小于纯碳模型预测。
3. 数据准备:驱动文件质量决定模拟结果上限
模型再精巧,最后输出的好坏很大程度上取决于输入数据。我见过太多人把大量时间花在调参上,结果气象数据本身就有严重的系统误差,越调越发散。这一节我把输入数据的准备要点逐条列清楚。
3.1 气象驱动数据:七个变量一个都不能少
Biome-BGC 每天需要七类气象数据:日最高温(°C)、日最低温(°C)、日均温(°C)、日降水量(cm)、水汽压亏缺(Pa)、太阳辐射(MJ/m²/day)、日长(秒)。其中水汽压亏缺通常用最高温、最低温和相对湿度计算,辐射也可以从日照时数估算,但如果站点数据缺这个,宁可先做插值再算,也别直接从再分析数据集随便拉一条就填进去。
时间序列的连续性比长度更重要。模型对有缺口的数据非常敏感,尤其是生长季中间断了一个月,那一年模拟出的 NPP 可能就废了。我的习惯是先用 ERA5 或 CRU 再分析数据补全缺测,再用站点实测做偏差校正。校正方法不复杂,按月的平均值和标准差做线性缩放即可。
3.2 CO2 浓度与氮沉降
CO2 驱动数据用全球平均的 Mauna Loa 观测序列就基本够用,但要注意和模拟时段对应,不要出现 2020 年模拟却用 1980 年的 CO2 浓度这种低级错误。氮沉降数据站点尺度可以用观测值,区域模拟通常用区域大气化学模型产品,分辨率可能较粗,空间插值时要注意梯度合理性。
3.3 初始化文件:平衡态是你最忠实的起点
初始化文件包含各碳库、氮库的初始值以及年龄、土壤质地等属性。建初始化文件最稳妥的方法不是手填,而是让模型先跑一个"spin-up":用驱动数据的前十年反复循环运行(或者用重复的多年平均气候),直到土壤碳氮库达到动态平衡,再把此时的输出状态写成新的初始化文件。这个过程在文档里叫 initialize,实际操作时建议至少循环 2000 年以上,不然高纬度冷湿地区的土壤碳库根本稳定不下来。
我习惯用这样的判断标准:最后 100 年的土壤碳总量变化不超过 0.1%,才算平衡态达标。如果差得多,就加大预热年份,别急着往下跑。
4. 参数化:生态生理参数是模型的灵魂
如果说气象数据决定了模拟的上限,生态生理参数就决定了模型能不能真正代表你研究的那片林子或草地。Biome-BGC 每个植被功能型(PFT)都对应一套生态生理参数,这些参数才是模型的灵魂。
4.1 核心参数清单
我用一个表格列出我每次都要重点检查的参数,这些参数对模拟结果的敏感度极高:
| 参数 | 含义 | 对输出的影响 |
|---|---|---|
| 比叶面积 SLA | 叶片面积与干重比(m²/kg) | 直接控制 LAI 和光截获,影响 GPP |
| 叶碳氮比 | 叶片碳与氮的比值 | 决定光合酶氮含量,限制光合速率 |
| 最大气孔导度 | 单位叶面积最大气孔导度 | 影响蒸腾和碳同化耦合 |
| 展叶日/落叶日 | 物候起止日期 | 决定生长季长度,强烈影响年 NPP |
| 细根周转率 | 细根更新速率 | 改变碳分配和土壤碳输入 |
| 分配系数 | 各器官碳分配比例 | 直接决定各库增长速率 |
| 木质素含量 | 叶片和凋落物木质素比例 | 影响分解速率和氮固定 |
4.2 参数从哪来
最理想的是直接测自己研究样地的叶片氮含量、比叶面积、生物量分配。但多数人做不到全量实测,那就优先参考 White 等人 2000 年那篇经典文献给出的全球各生物群系参数表,这是 Biome-BGC 参数化引用率最高的一套数据。另外实测数据能测多少算多少,比如比叶面积和叶片碳氮比这两项,野外采集叶片后回实验室一两个星期就能出数据,性价比极高。
如果连实测都做不了,至少要做一次敏感性分析,把影响最大的参数识别出来,用文献或遥感数据加以约束。我在黄土高原做过一个简单的一次一个变量(OAT)敏感性测试,结果发现 SLA、展叶日期和最大气孔导度排前三。那就重点约束这三个,其余参数用文献默认值也能接受。
4.3 用遥感数据同化物候参数
物候参数(展叶日、落叶日)靠人工设定容易跑偏,尤其空间范围一扩大,每个像元都用同一日期显然不合理。我的替代方案是用 MODIS NDVI 时间序列反演每个像元的生长季开始和结束日期,反算成模型需要的物候参数,再代入模型。这个方法比手工设定好一截,跑出来的 LAI 季节曲线的峰型、相位都能和遥感 LAI 对得比较好。
5. 站点尺度运行完整步骤:从零开始跑通一个案例
我以一个假设的温带落叶阔叶林站点为例,从下载源码到解读输出,把完整操作过程过一遍。Biome-BGC 的代码可以从官方渠道拿到,常见的是 Fortran 版本,我用的是 4.2 系列,Linux 环境编译运行。
5.1 编译环境准备
首先确保系统有 Fortran 编译器(gfortran 即可)。源码包里有 Makefile,直接make就能编出可执行文件。如果遇到编译报错,大多是缺库或者 Makefile 里的编译器路径不对,改一下即可。Windows 用户建议用 WSL 或者直接在服务器上跑,省去很多环境折腾的时间。
5.2 目录结构和文件清单
建立工作目录后,至少需要四类文件:
- EPC 文件:存放生态生理参数的文本文件,一行对应一个参数。
- init 文件:存放初始碳氮库和站点属性。
- met 文件:气象驱动数据,格式为逗号分隔或固定列宽,含日序、年份、七类气象变量。
- 运行配置文件:指定模拟年份、输出变量、输出格式等信息。
运行配置文件是入口,它告诉模型读哪个 EPC、哪个 met 文件,模拟从哪年到哪年,以及每个输出变量要不要写。我习惯把输出设置为日尺度,后面再自己聚合到月、年,灵活度更高。
一个常用的运行命令示例:
./biome-bgc run_conf.txt5.3 输出文件解读
运行完成后会生成多个变量文件的日值序列。比较关键的输出变量包括:
- GPP、NPP、NEE(碳通量)
- 叶面积指数 LAI
- 蒸散量 ET
- 土壤含水量
- 各碳库、氮库的大小
注意 NEE 的符号约定:模型输出为正值代表生态系统吸收碳(碳汇),负值代表释放碳(碳源)。不同版本符号约定可能不同,一定要看头文件注释,我吃过这个亏,正负号整反了差点把一个区域碳汇算成碳源。
6. 校准验证:如何让模型输出经得起审稿人推敲
跑通模型只是第一步,真正费时间的是校准和验证。审稿人不会因为你用了个机制模型就觉得结论可信,他们要看评价指标。
6.1 用涡度相关数据做基准
FLUXNET 的涡度相关通量塔数据是校准碳通量模拟的金标准。拿到站点数据后,先把实测 NEE 插补成连续日值,再和模型输出的日尺度和月尺度 NEE 对比。评价指标我习惯用三个:相关系数 R²、均方根误差 RMSE、Nash-Sutcliffe 效率系数 NSE。
我自己的经验是,日尺度能跑到 R²>0.6、NSE>0.5 就算不错了,月尺度应该能上到 R²>0.8。如果日尺度表现很差,先检查物候参数和水分胁迫逻辑,这两处出问题的概率最大。
6.2 参数敏感性分析与自动校准
手动调参太累且主观,我现在都建议先用敏感性分析筛出高敏感参数,再做自动校准。常用的工具是 PEST 或基于贝叶斯的算法。原理不复杂:给参数一个合理的先验范围,迭代求解使模拟值和观测值的加权残差最小。只校准前五到八个最敏感的参数,参数太多会过拟合,结果换个站点就崩。
校准一定要分批。先校准物候(展叶日、落叶日),再校准 SLA 和气孔导度,最后调分配系数和周转率。一次全上,容易出现参数补偿现象,比如把 SLA 调高掩饰了气孔导度偏低的问题,离开了这个站点就全乱了。
6.3 验证中的常见陷阱
验证数据和时间序列要有独立性,别用校准时段的数据来验证,至少要留出两年的数据做非校准验证。另外要注意时间尺度匹配,站点通量塔的足迹范围和模型模拟的站点范围不完全一致,地形复杂时误差很大。还有,通量塔观测值本身就有不确定性和能量不闭合问题,不必追求完全重合,重点看季节变化趋势和年际动态是否一致。
7. 应用场景组合拳:碳汇评估、气候变化响应、生态恢复
模型校准好之后才能谈应用。Biome-BGC 的典型应用方向主要集中在以下几类。
7.1 国家或区域碳汇评估
用站点校准后的参数外推到同类植被覆盖区,结合网格化气象数据和土壤质地数据,可以估算区域尺度的植被碳汇。这类结果常用于碳收支清单、森林资源清查的补充。但外推时要格外小心 PFT 分类的粗放问题,比如把落叶松林和常绿针叶林混为一类,模拟出的碳汇差异能差 30%。
7.2 气候变化情景模拟
BIome-BGC 对 CO2 浓度上升、温度升高的响应机制明确,很适合做未来气候情景模拟。做法是用 CMIP6 等气候模式输出的未来气象数据驱动模型,对比基准期和未来期的 NPP、碳汇变化。注意要用偏差校正后的气候输出,直接使用气候模式的原始气象预报,辐射和降水的系统偏差会让模型输出彻底失真。
7.3 造林与生态修复效益评价
在退耕还林、天然林保护等生态工程的评估中,Biome-BGC 可以模拟不同造林树种在给定立地条件下的碳积累潜力。这时候最关键的是选对 PFT 参数,还要考虑幼龄林与成熟林的碳分配差异。模型默认的分配参数偏成熟林,模拟幼龄林前期碳汇会严重偏高,需要针对年龄结构修正分配系数。
8. 踩坑记录:文档里不会写的那些事
最后分享几个我实际运行中踩过的坑,每一个都花了不小代价才爬出来。
8.1 初始化年份的 "伪平衡"
spin-up 跑了三千年,以为平衡了,结果一接真实气象数据,前十年土壤碳持续下降。原因是种子数据里土壤碳初始值和该地区实际承载能力差太远,预热循环的气象序列又恰好降雨偏多,掩盖了问题。现在我的做法是预热完成后,把最后一天的状态文件保存下来,再用偏差较大的另一套气象数据跑 500 年,观察土壤碳是否还能回来。能回来才是真平衡。
8.2 气象文件的单位陷阱
模型要求降水单位是厘米,辐射单位是 MJ/m²/day,VPD 单位是 Pa。ERA5 下载的是米制单位,如果不转换,降水会差 10 倍,模拟出的 LAI 直接就崩了。每次换新数据源,先花十分钟画一下年均气象数据的空间分布图,确认量级合理再进模型。
8.3 版本之间的输出差异
Biome-BGC 4.1 和 4.2 在氮沉降处理上有差异,输出变量命名也不完全一致。我在一个项目里用到 4.1 的气象预处理脚本,结果读 4.2 的头文件,对不上列数,调试了一整天。项目中途尽量不要换版本,非换不可的话,所有输入输出格式和正负号约定都要重新过一遍。
8.4 大面积网格模拟的存储与计算策略
区域模拟经常会遇到几十万个像元,每个像元独立跑一遍几十年的模拟,计算量不容小觑。我的经验是先按 PFT 和气候带分桶,合并相同参数和气象序列的像元,再并行计算。输出结果只保留年尺度变量,日尺度文件硬盘撑不住。另外,检查点(checkpoint)机制非常有用,每十年写一个中间状态,程序意外中断时不用从头再来。
结尾:一点个人体会
做模型的人常常被问"你的模拟准不准"。说实话,这类过程模型的准确性从来都不是和实测数据逐点相减那么简单的。Biome-BGC 的价值在于它提供了一个机制框架,让你能拆解每个因子对碳氮水循环的相对贡献。比起"准不准",我更愿意问"这个模型帮助你理解了什么"。在实际操作中,我最珍惜的是每次模型失败的过程,它往往比成功结果更能暴露你对生态系统机制理解的盲区。如果你正在用这个模型碰壁,别急着抄别人的参数,回到站点数据本身,看看是气象驱动错了,还是参数没表达出这片地真正的性格。模型只是一种思考方式,而数据才是你和真实生态系统之间最短的那座桥。