我最近几个月一直在折腾一个看起来有点“复古”的模型——Ghil-Sellers能量平衡模型。说它复古，是因为这模型比现在动辄几十万行代码的大气环流模式（GCM）老了半个世纪，可它模拟出来的结果，却让人惊讶地“现代”：只要调一个参数，就能看着地球从冰封状态突然“融化”，或者同一个太阳常数下出现多个稳定温度状态。这篇文章就聊聊我是怎么用Matlab把它从论文公式变成一整套能跑、能出图、能玩参数扫描的仿真程序，以及在这个过程里踩过的一堆坑。如果你是学气候、地球物理或者数值计算的学生，或者手里正好有个类似的项目要做，这篇应该对你有用。
1. Ghil-Sellers模型到底在算什么事？
1.1 为什么叫“能量平衡模型”，它和大气环流模型差在哪
地球气候系统的数值模拟，从方法论上大概分两条路线。一条是直接求解流体力学方程，把大气、海洋、陆面过程全部耦合起来，这就是我们常说的GCM，跑一次实验动辄几千核并行几个月。另一条就是今天要聊的“能量平衡模型”（Energy Balance Model, EBM），它不关心风怎么吹、云怎么飘，只盯着一个变量——地表温度T，用一维或零维的能量收支方程描述气候状态。
Ghil-Sellers模型就是EBM家族里非常经典的一个版本。它的名字里有两个人：Sellers在1969年提出了一维能量平衡模型，考虑了冰-反照率反馈；Ghil后来把这个模型做了数学化处理，用非线性动力学方法分析，发现这个模型在太阳常数变化时会出现多个平衡态和滞后现象。换句话说，Sellers提供了物理骨架，Ghil把它变成了一个可以严格讨论分岔、稳定性和敏感度的动力系统。
这模型最大的价值，不是精度，而是“可解释性”。GCM里的温度变化是成千上万个变量相互作用的结果，很多时候你很难说清楚“到底是谁导致的”；但EBM里每个项都能拆开看，温度变化就是辐射收支和热输运之间的拔河，物理图像非常清楚。
1.2 核心方程与物理参数：从偏微分方程到可解形式
我实现的一维Ghil-Sellers模型，基本形式是这样：
C(x) ∂T(x,t)/∂t = ∂/∂x [ D(1-x²) ∂T/∂x ] + Q s(x) [1 - α(x,T)] - (A + B T)
其中x=sin(lat)，是纬度的正弦，这样处理的好处是极区网格点自然加密。T是纬向平均的地表温度，C(x)是单位面积热容，D是扩散系数，代表大气和海洋的水平热量输送。Q是太阳常数，s(x)是太阳辐射的纬度分布函数，α(x,T)是行星反照率。最后一项A+BT是向外长波辐射的线性参数化，B就是这个模型的“辐射反馈强度”。
这个方程本质上就是一个带非线性源项的热扩散方程。非线性来自反照率α对T的依赖：温度低，冰雪覆盖面积大，反射率高，吸收的太阳辐射少，温度更低，这就是冰-反照率正反馈。
我用的Sellers参数版本里，反照率函数是这样的：
alpha = 0.62 - 0.42 * tanh((T - T_c) / delta);T_c是冰覆盖的临界温度，delta控制过渡带的宽度。用tanh而不是原论文里的阶跃函数，是为了让系统可微分，后续求雅可比矩阵、用隐式求解器都方便。
1.3 冰-反照率反馈：让模型“活”起来的非线性项
把反照率做了温度依赖之后，整个系统的行为就完全不一样了。你可以把冰-反照率反馈理解成一个“开关”：
- 当全球温度较高时，反照率低，系统吸收更多辐射，温度维持在高位；
- 当全球温度较低时，反照率高，系统吸收更少辐射，温度进一步降低；
- 但在中间过渡带上，同一个小扰动可能把系统推向完全不同方向。
Ghil最关键的发现，就是在这个非线性反馈下，同一个太阳常数Q可能对应多个稳定平衡解。有的解对应“无冰地球”，全球温度很高；有的解对应“雪球地球”，全球大部分被冰覆盖。这就像同一个水杯放在桌上，既可以稳稳站着，也可以倒扣着，都算是稳定状态，但能不能从一个状态到另一个状态，取决于你推它的力度。
这个特性，就是整个模型最值得模拟、也最能体现“气候临界点”概念的地方。
2. 用Matlab把模型搬进电脑：离散化、求解器与代码骨架
2.1 空间离散化：把连续纬度切成一串网格点
要把偏微分方程变成Matlab能算的东西，第一步是在空间上离散化。我选了等间距的x网格，从-0.999到0.999，取了41个点。为什么不是恰好从-1到1？因为x=±1是方程里扩散项的奇异点，(1-x²)在边界上为0，直接落在边界会导致雅可比矩阵退化，取靠近边界的点可以避免这个问题。
扩散项的处理我用了中心差分，但注意系数D(1-x²)不是常数，所以离散形式要写成通量守恒的形式：
function dTdt = ebm_rhs(t, T, p) x = p.x; N = length(x); dx = x(2) - x(1); dTdt = zeros(N, 1); % 扩散通量: F = -D*(1-x^2)*dT/dx F = zeros(N+1, 1); for j = 2:N-1 F(j+1) = -p.D * (1 - x(j)^2) * (T(j+1) - T(j)) / dx; end for j = 2:N-1 dTdt(j) = (F(j) - F(j+1)) / (p.C * dx); dTdt(j) = dTdt(j) + (p.Q * p.s(x(j)) * (1 - alpha(T(j), p)) - (p.A + p.B * T(j))) / p.C; end end其实我这里为了可读性，用循环写通量，并没有做向量化优化。41个网格点规模很小，这样写完全够快。如果你要上几百个点，再改用矩阵运算。
2.2 时间积分：为什么我最后选了自适应ode15s而不是自己写欧拉
刚开始我图省事，用了显式欧拉。结果发现一个问题：扩散项的时间步长限制极其严苛。D的量级大约是0.2 W/m²K，在41个网格点上，显式格式的稳定条件大约是Δt < 0.05年。这个限制并不是因为物理过程需要这么小步长，纯粹是数值稳定性逼着你走这么小。跑几十年模拟，显式欧拉需要几千步，每一步还有截断误差累积，到后面温度场会越来越“毛糙”。
后来我换成了ode15s，这是一个变步长的隐式求解器，专门处理刚性问题。隐式格式的好处是没有显式稳定性限制，步长可以自适应。我跑一次100年的模拟，ode15s大概只需要几十步，而且每步的误差有控制，算出来的曲线明显更光滑。
提示：如果你用的是Matlab R2020以后的版本，也可以用ode23t或ode23tb，它们对轻度刚性问题的性能表现也不错。但我最终留下ode15s，原因只有一个——它最稳，参数怎么乱改都不至于直接发散。
2.3 基础代码结构与参数表
下面是我整理的一份参数结构体定义，基本照搬Sellers原始论文量级，再按常用EBM文献做了一点微调：
p.x = linspace(-0.999, 0.999, 41)'; p.N = length(p.x); p.D = 0.2; % 扩散系数 W/(m^2 K) p.C = 5.0e7; % 热容 J/(m^2 K)，约等于50米海洋混合层 p.Q = 340; % 全球平均入射太阳辐射 W/m^2 p.A = 190; % OLR拟合常数 W/m^2 p.B = 2.0; % OLR温度反馈系数 W/(m^2 K) p.Tc = -5; % 反照率过渡带中心温度 °C p.delta = 5; % 反照率过渡带宽度 °C p.s = @(x) 1 - 0.241 * (3 * x.^2 - 1); % 太阳辐射纬度分布几个参数需要特意解释。C=5e7对应的是“深度混合层海洋”的热容量，半个世纪时间尺度上影响温度演变速度。D=0.2代表每单位温度梯度能输送多少热量，这个值决定了极地和赤道的温差。Tc和delta共同控制冰-反照率反馈的强度，是后面调参的重头戏。
3. 跑通第一组模拟：温度剖面、双稳态与太阳常数扫描
3.1 默认参数下的平衡温度分布
第一次跑通程序，我做了个简单的测试：给一个均匀的初始温度场，比如T=15°C，然后积分300年，看温度分布会收敛到什么状态。
初始三条不同温度的曲线，积分到最后全都收敛到同一个平衡剖面。这个剖面形状非常“教科书”——赤道温度接近295K（约22°C），极地温度大约245K（约-28°C），赤道和极地温差约50K。这个温差比真实地球略大一点，原因是我这个版本没有考虑海洋感热输运的细节，扩散系数D偏小。
用Matlab画出来非常直观，横轴从南极到北极，三条彩色虚线从初始状态一路演变到最终的黑色粗线，中间能看到缓慢的“传热”过程。这个图放到论文里就是一张很好的“模型验证”图。
3.2 连续变化太阳常数，捕捉“雪球—无冰”突变
跑通基本状态后，正戏来了：扫描太阳常数Q。我从Q=300扫描到Q=420，步长2 W/m²。对每个Q值，先用上一个Q的平衡态作为初始条件，然后积分到收敛，记录最终的温度。正向扫一遍，再从Q=420反向扫回300，把两条曲线叠在一起。
结果出来的那一刻，确实很震撼：在Q大约340~360区间内，正向扫描终点的全球平均温度是280K左右，反向扫描却停留在250K左右。同一个太阳常数，两个完全不同的气候态，中间隔着一条“不可跨越”的缝隙。这背后就是前面提到的冰-反照率反馈在临界点附近的“爆发”。
这个滞后回线（hysteresis loop）就是Ghil模型最经典的输出。你把这张图画出来，基本就等于复现了上世纪70年代那篇著名论文的核心结果。
3.3 从分岔结果反推气候敏感度
扫描完Q之后，我顺手算了一下气候敏感度，就是平衡态温度对太阳常数变化的响应速率。在高温分支上，dT/dQ大约0.4 K/(W/m²)，对应2×CO₂增温（约3.7W/m²强迫）的增温约1.5K。在低温分支上，dT/dQ明显更大，因为冰雪覆盖区大，反照率反馈更强，增温可以达到3K以上。
注意：这里的气温敏感度是“冰雪动力敏感度”，比实际气候系统的敏感度要小。真实气候还有水汽反馈、云反馈、碳循环反馈，EBM里都没有。所以在解释结果时，要特别注明这是“仅辐射-反照率反馈下的敏感度”，不能直接外推到真实气候。
4. 实操中的坑：数值发散、初值依赖、边界条件与调参技巧
4.1 显式扩散的稳定性条件，以及我踩过的CFL坑
这是我在这个项目里踩的第一个坑，也是新手最容易忽略的。
之前用显式欧拉，把扩散系数调到0.35后，程序直接“爆炸”。温度场在相邻网格点间出现“锯齿”状摆动，振幅越来越大，最后NaN。这就是扩散方程典型的违反CFL条件的表现。
扩散方程的显式格式稳定性条件大概是：
Δt ≤ (Δx)² · C / (2D)
我41个点的网格，x方向跨度2，Δx = 0.05，C=5e7，D=0.2，这个条件算出来大约是Δt ≤ 0.003年，也就是大约1天。你想想，我要积分300年，这就是显式欧拉最大的问题。
所以后来我直接放弃显式格式，改用ode15s。如果你一定要用显式格式，那就老老实实按这个公式推算步长，别指望“慢慢调”。
4.2 反照率台阶函数导致求解器“卡死”的处理
Sellers原始论文里，反照率对温度的分段函数带有一个发射率跃变，T在-5°C到0°C之间时，反照率像台阶一样跳变。这种不光滑的函数放进ode15s里，会导致求解器频繁缩小步长来解析突变点，严重的时候200个时间步都跑不完。
我后来用tanh函数把反照率过渡带从“突变”改成了“平滑过渡”，把delta设为5°C。这一改，积分速度提升了十倍不止，而且结果几乎没有变化。Ghil在1976年的论文里也提到过类似的处理，他用了一个连续可微的简化函数来分析系统的不动点。
提示：如果你只想复现“突变型”反照率的行为，也可以在反照率函数里用if-else写分段逻辑，但最好设置一个很小的过渡区间，给求解器留一点“呼吸空间”。完全不可微的函数的后果，就是ode15s会把大量时间花在检测事件点上。
4.3 边界条件怎么设：对称边界还是周期边界
一开始我在x=±1处用了“绝热边界条件”（零通量），也就是假设极地没有跨极热量输运。这听起来挺合理，但实际运行发现，最高纬度的温度会比相邻网格点低很多，形成一个“边界冰帽”畸变。
后来改成“对称边界条件”，假设南极和北极处的温度梯度为0，且两侧对称。具体做法是虚构两个网格点：x(-1)点上的温度等于x(1)点上的温度，然后把边界通量设成0。这样做之后，极地温度分布就平滑多了。
在这里推荐一个更稳定的做法：在离散的时候直接让方程系数在边界处趋近于0。所以我在代码里把边界网格点取在x=±0.999而不是±1，这样边界通量F自然就是0，不需要额外处理。
4.4 参数调优的几条经验
参数调优的核心原则：先固定一切参数，只改变一个量，观察系统行为变化。这个项目的关键参数优先级排序如下：
- Q（太阳常数）：影响系统是单稳态还是双稳态，是全局性参数。
- Tc（冰覆盖临界温度）：决定反照率反馈的触发位置，对温度分布形态影响很大。
- D（扩散系数）：决定极地与赤道的温差，温差过小或过大，平衡曲线形态会很怪。
- delta（过渡带宽度）：影响分岔回线的“拐角”锐度，但不会改变临界Q的大致位置。
我实测中发现，Tc从-5°C调到0°C，滞后回线的宽度会明显增加；delta从5°C调到15°C，回线会变窄甚至消失。这说明反照率过渡带的平滑度，直接决定了系统是否还保留双稳态。如果你发现模拟跑不出双稳态，先检查delta是不是调太大了。
5. 这个模型还能往哪走：从教学玩具到科研前哨
5.1 加季节循环、云反馈与海洋热输运
这个模型虽然是“玩具”，但它的框架完全支持继续加东西。我给模型加过两个扩展：
第一个是季节循环。把太阳辐射s(x)从固定分布改成随时间变化的函数，加入自转轴倾角引起的季节变化，模型就能输出“最热月”和“最冷月”的温度分布。这个扩展对研究冰盖季节性融化很有意义，代码改动也不大，只要把s(x)变成s(x, t)就行。
第二个是云反馈。在OLR项后面加一个云辐射强迫项，让云的覆盖比例随着温度变化。比如温度升高时云增多，反射太阳辐射，削弱增温趋势。这种负反馈对双稳态回线的形状有直接影响，做起来也不复杂。
不过要提醒一下：每加一个自由度，你就得多面对一个“参数不确定”的问题。EBM的优势就在于参数少、机理清楚，加太多东西反而会失去这个优势。
5.2 和GCM/LSMs差异对比
如果你以后要往更复杂的模型方向走，可以把Ghil-Sellers模型当作“验证单元”来用。
比如，你可以在GCM的输出里提取全球平均温度剖面，再放到EBM里跑一遍，看看EBM能否再现GCM的气候敏感度。如果差得远，说明某个反馈过程在EBM里没有被正确参数化，这就给了你一个定位问题的线索。
我见过一些做古气候模拟的研究组，把EBM当作“探路兵”：先跑几百个参数组合找关键区间，再用GCM做精细模拟。2000行不到的代码，但作用并不比几十万行的GCM小。
5.3 给Matlab学习者的建议
最后，说点给准备复现这个模型的Matlab初学者的建议：
- 先跑通最基本的零维模型，也就是去掉空间扩散项，只算全球平均温度。这个只有一行微分方程，最容易验证思路。
- 再扩展到一维。扩散项用通量守恒形式写，不要太早优化性能。
- 调参时做“扫描图”，把结果保存成结构体或表格，用subplot把温度剖面和分岔图放在同一张图里看。
- 尽量用函数句柄把参数和右端项封装起来，后续做参数扫描会省很多事。
我一直觉得，Ghil-Sellers模型是学习气候模式数值模拟最好的“第一辆车”。它结构简单，却包含非线性、分岔、刚性问题、空间离散化这些核心概念，每一个都值得反复琢磨。把这个模型吃透，再去看那些大规模气候模式的文档，你至少能看懂它们在做什么，以及为什么会“跑飞”。
` 编排与节奏控制)