用Matlab实现Ghil-Sellers能量平衡模型：从双稳态到气候突变模拟

我最近几个月一直在折腾一个看起来有点“复古”的模型——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)其中xsin(lat)是纬度的正弦这样处理的好处是极区网格点自然加密。T是纬向平均的地表温度C(x)是单位面积热容D是扩散系数代表大气和海洋的水平热量输送。Q是太阳常数s(x)是太阳辐射的纬度分布函数α(x,T)是行星反照率。最后一项ABT是向外长波辐射的线性参数化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(N1, 1); for j 2:N-1 F(j1) -p.D * (1 - x(j)^2) * (T(j1) - T(j)) / dx; end for j 2:N-1 dTdt(j) (F(j) - F(j1)) / (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); % 太阳辐射纬度分布几个参数需要特意解释。C5e7对应的是“深度混合层海洋”的热容量半个世纪时间尺度上影响温度演变速度。D0.2代表每单位温度梯度能输送多少热量这个值决定了极地和赤道的温差。Tc和delta共同控制冰-反照率反馈的强度是后面调参的重头戏。3. 跑通第一组模拟温度剖面、双稳态与太阳常数扫描3.1 默认参数下的平衡温度分布第一次跑通程序我做了个简单的测试给一个均匀的初始温度场比如T15°C然后积分300年看温度分布会收敛到什么状态。初始三条不同温度的曲线积分到最后全都收敛到同一个平衡剖面。这个剖面形状非常“教科书”——赤道温度接近295K约22°C极地温度大约245K约-28°C赤道和极地温差约50K。这个温差比真实地球略大一点原因是我这个版本没有考虑海洋感热输运的细节扩散系数D偏小。用Matlab画出来非常直观横轴从南极到北极三条彩色虚线从初始状态一路演变到最终的黑色粗线中间能看到缓慢的“传热”过程。这个图放到论文里就是一张很好的“模型验证”图。3.2 连续变化太阳常数捕捉“雪球—无冰”突变跑通基本状态后正戏来了扫描太阳常数Q。我从Q300扫描到Q420步长2 W/m²。对每个Q值先用上一个Q的平衡态作为初始条件然后积分到收敛记录最终的温度。正向扫一遍再从Q420反向扫回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.05C5e7D0.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模型是学习气候模式数值模拟最好的“第一辆车”。它结构简单却包含非线性、分岔、刚性问题、空间离散化这些核心概念每一个都值得反复琢磨。把这个模型吃透再去看那些大规模气候模式的文档你至少能看懂它们在做什么以及为什么会“跑飞”。