直接进入正题。近几年做气候模拟相关的工作时我经常被问到“有没有一个既简单又能抓住核心物理的模型可以练手”大多数时候我的回答都是先从一维能量平衡模型EBM开始。尤其是Ghil-Sellers模型它是我见过性价比最高的一类入门框架——代码量不大物理直觉完整还能复现出真实气候系统里“稳定态切换”和“滞后回线”这种教科书级的现象。这篇文章就围绕如何用Matlab把Ghil-Sellers能量平衡模型从零搭起来从方程推导、离散化到参数扫描把每一步细节和踩过的坑都摊开讲清楚。Ghil-Sellers模型本质上解决的是一个看似矛盾的问题地球从太阳接收的能量和向外太空辐射的能量大体平衡但为什么赤道那么热、极地那么冷这中间的温差由谁维持、由谁调节如果把全球平均温度当成唯一变量你只能看到“总收支”要想看到冰线位置、经向热量输送、冰-反照率反馈这些真正决定气候状态的结构就必须把纬度维度加进来。这个模型的价值就在于此它用一个关于纬度的扩散型能量方程把辐射、反照率、热输送三件事耦合在一起然后用Matlab几分钟就能跑出结果。它的适用范围很广对气候初学者它是理解辐射强迫和反馈机制的直观工具对做研究的它可以用来分析分岔点、临界过渡和气候敏感性问题对Matlab练习者它又是一个非常适合练手的有真实物理意义的偏微分方程求解案例。我下面从模型本身的物理设计开始讲然后再进入Matlab实现和实验分析。1. 为什么是能量平衡模型从零维到一维的关键一步很多资料常用零维模型来解释地球温度把地球当成一个均匀小球吸收的太阳辐射等于放出的长波辐射解一个代数方程得到约-18℃的全球平均温度再考虑温室效应修正到15℃左右。这个练习对理解辐射平衡很重要但它完全丢掉了气候的空间结构。Ghil-Sellers模型把“空间”这个维度加了回来不过只保留纬度方向经度方向做了平均这种设计跟大气环流的纬向平均特征密切相关。1.1 赤道到极地一个被温度差驱动的热机从观测上看赤道地区的净辐射收入大约是极地的两到三倍但赤道气温并没有因此无限升高极地也没有无限冷却原因就在于大气和海洋会把热量从低纬往高纬输送。这种输送在方程里可以用类似热扩散的项来近似温度梯度越大输送越强方向永远是从暖处到冷处。Ghil-Sellers模型用了一个带纬度的扩散算子来描述这个过程这是整个模型从“零维”走向“一维”的最本质改变。当初Sellers在1969年提出这个模型时核心目的就是为了研究“如果太阳辐射发生微小变化冰线会怎么移动”这类问题。它比零维模型多出来的东西恰恰是气候系统里最关键的现象冰-反照率反馈。大家可以这样理解当极地温度降低冰雪覆盖范围扩大地表反照率升高反射掉更多阳光地球吸收更少能量温度进一步降低冰继续扩大。这是一个正反馈过程。如果没有经向温度梯度和冰线这个反馈根本建立不起来。1.2 为什么选择Ghil-Sellers而不是Budyko模型在能量平衡模型这个家族里最经典的除了Ghil-Sellers还有Budyko模型。两者的差异主要在于经向热输送的参数化方式Budyko用了线性化的输送项即某纬度的温度直接与全球平均温度之差成正比Ghil-Sellers用的是扩散项更接近真实热传导的图像。扩散形式的优点是可以自然地在两极边界形成光滑的温度分布而且数学上更容易做稳定性分析也能更自然地匹配一些现代气候模型的子网格参数化。从Matlab实现角度来看扩散项处理起来更直接——空间二阶导数的差分格式是标准操作边界条件也清晰。这也意味着后续如果你想扩展比如加入云量反馈、季节循环或者把全球平均温度耦合进来扩散框架的可扩展性都更好。我实际做下来Budyko模型半小时就能跑通但Ghil-Sellers模型能玩出的花样更多复杂度的提升又完全在可控范围内所以更推荐作为练手项目。2. 模型核心方程与参数化方案这个模型的主控方程写出来其实非常紧凑大概就是这个样子C * dT/dt Q/4 * S(phi) * (1 - alpha(phi, T)) - (A B*T) 1/cos(phi) * d/dphi [D * cos(phi) * dT/dphi]这个式子里每一部分都有明确的物理含义理解它比背公式更重要。我一项一项拆开讲然后给出各个参数的取值和依据。2.1 入射太阳辐射项为什么是全球平均的四分之一方程第一项 Q/4 * S(phi) * (1 - alpha) 是系统吸收的太阳辐射。Q是太阳常数地球轨道处垂直于太阳光线的辐照度约为1368 W/m²。地球对阳光的等效截面是πR²但表面积是4πR²所以分摊到单位表面积上的平均入射能是Q/4约342 W/m²。这个“除以4”的几何因子是能量平衡模型最常被忽略但又最基础的地方。S(phi)是一个纬度分布函数表示不同纬度真实接收到太阳辐射相对于全球平均的比值。赤道地区接近1.2左右极地可能只有0.6左右。为了数学上的方便通常把S(phi)归一化为全球平均等于1这样整个方程对全球积分时入射项恰好就是Q/4乘以(1-平均反照率)。具体形式可以用勒让德多项式展开S(x) 1 - 0.477 * (3*x^2 - 1) / 2, x sin(phi)这个近似形式在Sellers和许多后续文献里都出现过它很好地抓住了“赤道略高、极地偏低”的分布特征。注意直接取S(phi) 1 - 0.477*(3*sin(phi)^2 - 1)/2的时候在赤道得到约1.238在极点得到约0.762这跟实际辐射的年平均分布很接近。2.2 长波辐射项把复杂的辐射传输线性化方程里的 A B*T 代表地球向太空放出的长波辐射。真实的出射长波辐射是温度的四次方关系斯蒂芬-玻尔兹曼定律但考虑到大气温室气体、云层和温度变化范围有限在气候模拟里经常把它在当地温度附近做线性化。比如取A 204 W/m², B 2.17 W/(m²·K)这样当T288 K约15℃时出射长波辐射约等于204 2.17288 204 625 829 W/m²。但是等等这个数看起来比入射的342 W/m²大很多其实这里T用的是开尔文温度而不是摄氏温度用Klvin算下来确实会得到800多但要注意这是“黑体辐射等效”的简化并不代表净出射。有些资料会写成A B(T - 273.15)这样A大约取210B取2左右。我更推荐用开尔文做绝对计算但需要在写代码时统一好单位避免出现“用摄氏度套开尔文系数”这种低级错误。2.3 反照率参数化冰-反照率反馈的核心反照率alpha在这里不是常数而是随温度和纬度变化。当某格点的温度低于结冰阈值时认为地表被冰雪覆盖反照率高典型值取0.62温度高于阈值时不结冰反照率低典型值取0.30。我建议在程序里用平滑过渡函数而不是硬阶跃比如alpha alpha0 0.5 * (alpha1 - alpha0) * (1 tanh((T - Tc)/delta))其中alpha0是无冰反照率0.3alpha1是冰雪反照率0.62Tc是冰点阈值我习惯取-10℃即263.15 K这代表的是实际冰雪能够持续存在的地表温度条件而不是纯水结冰温度delta是过渡宽度取1到2 K就够了。用平滑函数的理由是硬阶跃反照率会带来数值上的不连续在时间积分时容易在冰线附近产生振荡而且后续做稳态追踪和分岔分析时不连续的方程组会让求解器非常难受。tanh平滑化后的函数既保留了强非线性又让各种数值方法都能稳定工作。2.4 经向热扩散项大气和海洋的联合输送方程右边最后一项是经向热输送。D是扩散系数量纲是W/(m²·K)表示单位温度梯度下能输送多少热通量。为了体现球面几何扩散算子写成这个形式1/cos(phi) * d/dphi [D * cos(phi) * dT/dphi]在球坐标里cos(phi)的作用是度量因子纬度越高单位纬距对应的实际地表距离越小所以同样的温度梯度对应的热通量收敛/发散也不同。D的取值通常在0.3到0.6之间我一般取0.44 W/(m²·K)这是对现代气候条件下总经向热通量大气感热潜热海洋输送的一个综合估计。还有一个更巧妙的做法是把坐标从纬度phi换成x sin(phi)这样cos(phi)*d/dphi就变成d/dx扩散项简化为d/dx [D * (1-x²) * dT/dx]这个形式在做数值离散时非常清爽x在[-1,1]之间规则分布极点就是端点边界条件直接设为x±1处无通量。我强烈推荐用x坐标能少写不少代码也方便做谱方法扩展。3. Matlab实现从离散化到出图现在进入核心实操。我先把整体程序架构说一下参数定义、空间网格初始化、反照率函数、时间积分循环、结果可视化。整个脚本不需要依赖任何工具箱纯Matlab原生语法我经常在R2019b到R2023b的各个版本上跑都没问题。3.1 空间网格与扩散项差分格式使用x sin(phi)坐标空间网格可以这样建N 101; % 纬度网格点数取奇数包含赤道对称 x linspace(-1, 1, N); % 坐标x从南极到北极 dx x(2) - x(1);温度场T是一个N×1的列向量。扩散项需要中心差分在x坐标系里扩散通量F D * (1 - x²) * dT/dx散度就是dF/dx。我建议在交错网格上处理也就是在整数格点上存温度在半格点上算通量这样边界条件处理起来最干净。for i 2:N-1 % 半格点处的系数 xm 0.5 * (x(i-1) x(i)); xp 0.5 * (x(i) x(i1)); Dm D * (1 - xm^2); Dp D * (1 - xp^2); % 二阶中心差分 dTdxm (T(i) - T(i-1)) / dx; dTdxp (T(i1) - T(i)) / dx; diffusion(i) (Dp * dTdxp - Dm * dTdxm) / dx; end % x-1和x1处无通量边界diffusion(1)0, diffusion(N)0这里注意x坐标在-1和1时(1-x²)恰好为0所以即使dT/dx不为0通量也自动为0这就是极点的自然边界条件。这个性质非常好实现上几乎不用额外处理。3.2 时间积分与时间步长约束显式欧拉是最直观的方案但必须满足CFL条件。扩散项的稳定性条件是dt dx² / (2 * max(D/C * (1-x²)))这里D/C的量纲是1/s如果C取1e7 J/(m²·K)D取0.44 W/(m²·K)D/C 4.4e-8 s^-1。对于N101dx约0.02那么dt必须小于约1e4秒。这个限制非常严苛如果积分目标是100年就需要约30万步虽然每一步都是简单的向量运算但循环写不好就会很慢。我在实际中更推荐把时间方程写成一个函数句柄然后用Matlab自带的ode45求解。这个方程本质上是常微分方程组——对空间离散后每个格点的温度随时间变化都满足ODE所以ode45既能自动控制步长又能保证数值稳定性代码还更简洁[t, Tmat] ode45((t, T) ebm_rhs(T, params), [0, tspan_years*365*86400], T0);注意t的单位是秒因为各个常数项都是基于国际单位制算的。如果你觉得秒太大不好看可以自己在方程两边同时除以一年的秒数把时间基准改成“年”但要注意C也会同步缩放。我个人不建议改单位直接用秒输出的时间轴再除以3.156e7转成年更不容易出错。3.3 反照率平滑与冰线识别把反照率作为T的函数单独写一个子函数用tanh平滑。在算完稳态温度场之后如果需要找冰线即温度等于Tc的纬度可以在x坐标上找跨过Tc的点然后线性插值Tc 263.15; % 10°C in Kelvin crossings find(diff(sign(T - Tc)) ~ 0); for k 1:length(crossings) i crossings(k); x_ice interp1(T(i:i1), x(i:i1), Tc); lat_ice asin(x_ice) * 180 / pi; end这里输出的是绝对值由于模型在南北半球对称理论上应该得到两个对称的交叉点。如果你只取其中一个记得用abs或者只处理北半球。3.4 完整代码示例一套直接能跑的脚本下面我给出一个我自己调试过、可以直接运行的核心脚本重点是整体结构和主循环。为了让读者看清逻辑我稍微简化了一点输出部分但方程的物理部分都是完整的。% Ghil-Sellers 一维能量平衡模型 % 坐标取 x sin(lat)空间二阶中心差分时间用ode45 clear; clc; close all; % ------- 参数定义 ------- params.Q 1368; % 太阳常数 W/m^2 params.Acoef 204; % 长波辐射线性化系数 W/m^2 params.Bcoef 2.17; % 长波辐射线性化系数 W/(m^2*K) params.D 0.44; % 经向扩散系数 W/(m^2*K) params.C 1e7; % 热容量 J/(m^2*K) params.alpha0 0.30; % 无冰反照率 params.alpha1 0.62; % 冰雪反照率 params.Tc 263.15; % 冰线温度阈值 K params.delta 2.0; % 反照率过渡宽度 K params.S2 -0.477; % 太阳辐射分布勒让德系数 % ------- 空间网格 ------- N 101; x linspace(-1, 1, N); dx x(2) - x(1); params.x x; params.dx dx; % 入射辐射分布函数 S(x) 1 S2 * 0.5*(3*x.^2 - 1) params.Sdist 1 params.S2 * 0.5 * (3*x.^2 - 1); % ------- 初始温度取纬度线性分布 ------- T0 300 - 40*x.^2; % 赤道300K极地260K左右的初始猜值 % 如果对多个平衡态感兴趣可以换不同初值再试 % ------- 时间积分 ------- tspan [0, 100*365*86400]; % 模拟100年单位秒 [t, Tmat] ode45((t,T) ebm_rhs(t,T,params), tspan, T0); % 取最后时刻作为稳态或者取最后50年的平均 T_steady Tmat(end,:); % ------- 画图 ------- lat asin(x) * 180/pi; figure; plot(lat, T_steady - 273.15, linewidth, 2); xlabel(纬度 (deg)); ylabel(温度 (^{\circ}C)); grid on; title(Ghil-Sellers 模型稳态温度分布); % ------- 模型右端项核心方程 ------- function dTdt ebm_rhs(t, T, p) x p.x; dx p.dx; Nnum length(x); alpha p.alpha0 0.5*(p.alpha1 - p.alpha0)*(1 tanh((T - p.Tc)/p.delta)); absorbed p.Q/4 * p.Sdist .* (1 - alpha); olr p.Acoef p.Bcoef * T; % 扩散项d/dx [D*(1-x^2)*dT/dx] diffu zeros(size(T)); for i 2:Nnum-1 xm 0.5*(x(i-1)x(i)); xp 0.5*(x(i)x(i1)); Dm p.D*(1-xm^2); Dp p.D*(1-xp^2); diffu(i) (Dp*((T(i1)-T(i))/dx) - Dm*((T(i)-T(i-1))/dx))/dx; end dTdt (absorbed - olr diffu) / p.C; end好到这里一个能出稳态温度分布的“最小模型”已经成型了。上面这段代码如果复制到Matlab里应该能直接得到一条从赤道高温到极地低温的平滑曲线。我建议各位跑通之后立刻开始做下面的实验那才是这个模型真正好玩的地方。4. 数值实验与临界行为分析代码能跑只是第一步。Ghil-Sellers模型最迷人的部分是它能展示气候系统的不唯一性和突变。下面几个实验建议按顺序做你会发现每一步都在逼近真实的气候现象。4.1 基准实验暖态温度分布首先跑默认参数下的稳态。你会发现赤道温度大约在25到30℃之间极地温度在-30℃左右全球平均温度大约在10到15℃之间和真实地球气候比较接近。这个状态通常被称为“暖态”。可以尝试修改D值当D变大热输送增强赤道到极地的温差缩小当D变小温差拉大。这个实验直接让你感受到经向热输送对“全球温度格局”的控制力。我自己在第一次跑的时候犯过一个错初始温度给的是全球均匀300 K结果扩散项一开始几乎为零温度靠辐射项调整等稳定下来之后发现赤道温度还行但极地温度异常高。后来我才意识到在强非线性反照率反馈下初始场会影响最终落入哪个平衡态均匀初始场很容易把系统带到不真实的“冰雪线在极地”状态。所以初始条件建议用带温度梯度的线性或二次分布让模型更快进入合理区间。4.2 太阳常数扫描冰线的突变与滞后接下来做最关键的实验把太阳常数Q从0.9倍基准值逐渐增加到1.1倍再从1.1倍逐渐降回0.9倍每步都用上一步的稳态作为初值继续积分画冰线纬度随Q的变化曲线。这就是参数连续追踪结果出来你会看到典型的“滞后回线”。当Q下降时冰线逐渐向赤道方向移动温度越来越冷但在某个临界点附近系统会突然跳到“雪球态”——冰线直接冲到赤道全球被冰雪覆盖反过来从高Q的雪球态缓慢增大Q时系统并不会在原临界点跳回来而是等到Q高到一定程度才突然“解冻”。这两个临界点之间的区域系统存在两个稳定状态这就是所谓的“双稳态”和“分岔行为”。别小看这条滞后回线。它意味着即使外强迫完全回到原来水平系统也不会自动恢复原状。你可以在Matlab里复现它然后跟真实气候研究里讨论的“雪球地球”事件做对照——那是新元古代一次全球冰冻事件。当然真实过程要复杂得多但模型给出的机制性解释方向是和现代研究一致的。4.3 扩散系数对气候敏感性的影响如果你把D从0.3调到0.6会发现一个反直觉的现象更强的经向热输送让极地温度升高冰线向极地退缩全球平均温度也升高。原因在于更强的热输送让高纬地区变暖减少冰雪面积降低反照率地球整体吸收更多太阳辐射。这个正反馈被热输送所放大所以气候对太阳辐射变化的敏感程度会随D变化。这个实验值得做定量的分析定义气候敏感度为“全球平均温度变化量 / 太阳常数变化百分比”观察不同D下的变化。你可以用Matlab的polyfit做一个线性拟合会得到一个清晰的趋势。新手容易忽略的点是扩散系数D不是随便取的常数它代表了整个大气和海洋系统的热量输送效率不同D值对应的气候敏感度差异非常大这也是一部分气候模型“气候敏感度”不确定性的来源之一。5. 常见问题与调试实录最后把我在实际操作里遇到的坑集中整理一下。这些问题不踩一遍很难注意到但踩过之后对模型的理解会深一层。5.1 温度出现NaN或系统崩溃最常见的原因是时间步长过大导致显式格式发散如果你用ode45还出现NaN多半是方程本身有物理不合理的地方。我遇到过一次非常隐蔽的错误在反照率函数里传了摄氏温标但阈值却用了开尔文导致热带地区也判定为“冰覆盖”整个模型瞬间变成“雪球态”温度一路崩到负值。这类问题没有捷径只能把所有物理量写在同一张参数表里随时检查单位。建议做一个简单的“单位自查清单”入射项W/m²、长波辐射W/m²、扩散项W/m²、每项除以热容量之后量纲为K/s这样量纲平衡一眼就能验证。5.2 冰线“卡”在网格上位置不准确如果你用硬阶跃反照率冰线的位置会严重依赖网格间距而且随着时间步进会在相邻网格点之间跳来跳去。用tanh平滑过渡可以解决大部分问题但要注意delta不能太大否则冰线位置失真我常用的值是1到2 K既能保证数值稳定又不会把冰线“抹”得太模糊。另外尽量用插值法识别冰线而不是简单取最冷的“结冰网格”这样冰线在参数扫描时能连续变化画出来的滞后回线才平滑。5.3 参数扫描太慢怎么办在做滞后回线的时候如果每一步都用ode45从零积分到100年几百个参数点会非常耗时。一个技巧是每步都以相邻参数的稳态作为初值同时把时间跨度缩短到20到30年因为初值离稳态很近系统会很快收敛。另一个技巧是真正做分岔分析时干脆不用时间积分直接对稳态方程用fsolve求解。稳态方程就是把dT/dt置零后得到一个非线性方程组fsolve配合解析雅可比或者让你用自动微分可以秒级收敛。我个人的习惯是先用时间积分跑一次全局演化理解系统行为再用fsolve做精细的临界点定位。5.4 对称性破缺问题很多模型在南北半球对称的初始条件下得到的结果是严格对称的。但在某些参数区域微小扰动可能让南北极冰线不对称。这是真实系统里也可能发生的事比如北极和南极的冰覆盖确实不对称但在简化模型里出现通常意味着网格或参数数值误差放大了。如果想验证是否正常可以非对称扰动后观察冰线是否出现持续非对称稳态如果不确定就先加一个微小的白噪声扰动看系统是否会跳回对称解。我用的方法是给初始温度加一个幅度0.1 K的正弦扰动正常情况下模型会很快恢复对称。最后再分享一个实操建议给模型加一个简单的“诊断输出”函数每次积分一个时间周期就打印一次全球平均温度、冰线纬度和总能量收支残差。这样调试初期你会节省大量时间也能直观地感受到系统趋近稳态的松驰过程。这个模型的可扩展方向很多比如加云量参数化、季节强迫、以及用实测数据做订正都是很好的后续练习。我做下来最深的体会是这类看似“玩具”的模型其实是连接物理直觉和数值实现之间最短的桥。