问题引入¶
全球变暖是当今人类面临的最严峻挑战之一。每天都有关于极端天气、海平面上升的新闻报道。但你是否想过：
- 地球的温度是由什么决定的？
- 为什么工业革命以来全球温度持续升高？
- 如果不减少排放，未来会怎样？
要回答这些问题，离不开建立模型。物理学家和气候学家经过长期研究，发现可以用一个简洁的微分方程来描述地球温度的变化——这就是零维能量平衡模型。
本讲将从数学出发，逐步构建这个模型，用 Python 进行求解和可视化，最终与真实观测数据对比，预测未来气候走势。
准备工作¶
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy.integrate import solve_ivp
from ipywidgets import interact, FloatSlider, IntSlider
plt.rcParams['font.family'] = ['Noto Sans CJK SC', 'IBM Plex Sans']
知识回顾：常微分方程¶
记号：
$y=y(t)$ 是时间 $t$ 的函数
初始时刻 $t=0$，$y(0)=y_0$
记 $y'=dy/dt$。
导数是关于 y 的线性函数¶
$$y' = a - by$$
解：$y(t) = (- \frac{a}{b})(e^{-bt}-1) + y_0 e^{-bt}$
注释：当 $b \ne 0$ 时，令 $y'=0$ 求解 $y$ 可得平衡解 $y=a/b$。也可以通过在解中令 $t \to \infty$ 使指数项趋于零来获得。
导数是 关于 y 的线性函数加上关于 t 的强迫项¶
$$y' = a - by + f(t)$$
解：$y(t) = e^{-bt} \left( y_0 + \int_0^t e^{bu}(a+f(u)) du \right)$
用 Python 求解常微分方程¶
- 定义 $f(y, t)$
- 使用
scipy.integrate.solve_ivp，传入 $f$、初始值、时间跨度、参数 - 求解
- 绘图等
def forcing(c):
return lambda t: c * t
def ode_func(t, y, a, b, c):
return a - b * y + forcing(c)(t)
编程知识：闭包（Closure）
forcing(c)返回一个lambda函数lambda t: c * t。这个匿名函数捕获了外部变量c——即使forcing函数已经执行完毕，c的值仍然被内部函数"记住"。这种机制称为闭包。类似地，后面代码中的
lambda t: get_co2_concentration(t, 'high')捕获了字符串'high'，使得这个 lambda 在被solve_ivp调用时始终使用高排放情景。
交互演示：求解 ODE¶
@interact(
a=FloatSlider(min=0, max=10, step=0.01, value=0, description='a'),
b=FloatSlider(min=0, max=5, step=0.1, value=0, description='b'),
y0=FloatSlider(min=-5, max=15, step=0.1, value=2.0, description='y₀'),
c=FloatSlider(min=0, max=5, step=0.1, value=0.0, description='c'),
)
def plot_ode(a, b, y0, c):
sol = solve_ivp(
ode_func, [0, 10], [y0], args=(a, b, c),
dense_output=True, max_step=0.05
)
t_eval = np.linspace(0, 10, 500)
y_eval = sol.sol(t_eval)[0] # [0] 取第一个状态变量（本例为一维 ODE）
fig, ax = plt.subplots(figsize=(8, 5))
ax.set_facecolor('black')
ax.plot(t_eval, y_eval, color='red', lw=3)
ax.set_ylim(0, 10)
ax.set_ylabel('y')
# 方向场
xs_field = []
ys_field = []
lrx = np.linspace(0, 10, 30)
lry = np.linspace(0, 10, 30)
step = lrx[1] - lrx[0]
for x in lrx:
for y in lry:
v = np.array([1, a - b * y + forcing(c)(x)])
v = v / (20 * step)
xs_field.extend([x - v[0], x + v[0], np.nan]) # np.nan 断开线段，避免箭头首尾相连
ys_field.extend([y - v[1], y + v[1], np.nan])
ax.plot(xs_field, ys_field, alpha=0.7, color='yellow')
eq = 0 if b == 0 else a / b
ax.axhline(eq, color='white', ls='--')
ax.annotate('y₀', xy=(-0.5, y0), color='red', fontsize=12)
ax.set_title("y'(t) = a - by + forcing(c)(t) 的解")
plt.tight_layout()
plt.show()
interactive(children=(FloatSlider(value=0.0, description='a', max=10.0, step=0.01), FloatSlider(value=0.0, des…
背景：气候物理¶
最简单的气候模型可以概念化为：
$$\text{热量变化} = + \text{吸收的太阳辐射（来自太阳射线的能量）} - \text{向外热辐射（即向太空的黑体冷却）} + \text{人类引起的温室效应（截留的向外辐射）}$$
其中每一项都被解释为全球的平均值（因此是"零维"的）。
在地球距离太阳的轨道距离处，拦截地球的太阳射线的功率等于
S = 1368 # 太阳辐照度 [W/m²]（单位时间单位面积的能量）
一小部分
alpha = 0.3 # 反照率，即行星反射率 [无量纲]
在数学上只需写下一个微分方程，但在物理世界中需要识别物理变量。
在"烘烤地球"的例子中，将确定以下物理量：
- 工业革命开始年份：1850
- 1850 年平均温度：14.0 °C
- 太阳辐照度 $S = 1368$ W/m²：来自太阳的能量
- 反照率或行星反射率：$\alpha = 0.3$
- 大气和上层海洋热容量：$C = 51$ J/m²/°C
地球烘烤公式：
$$C \cdot \mathrm{temp}'(t) = S(1-\alpha)/4$$
print(f'S(1-α)/4 = {S * (1 - alpha) / 4:.4f}')
S(1-α)/4 = 239.4000
这些入射太阳辐射中的一小部分被反射回太空（通过白云、雪和冰等反射面），剩余比例 $(1-\alpha)$ 被吸收。
由于入射太阳射线在距离太阳这么远的地方大致平行，拦截它们的地球截面积只是一个面积 $\pi R^2$ 的圆盘。由于后续涉及的所有其他项都作用于球形地球的整个表面积 $4\pi R^2$，因此单位表面积吸收的太阳辐射（全球平均）减少了 4 倍。
因此，单位面积吸收的太阳辐射为
$$\text{吸收的太阳辐射} \equiv \frac{S(1-\alpha)}{4}$$
absorbed_solar_radiation = S * (1 - alpha) / 4 # [W/m²]
C = 51.0 # 大气和上层海洋热容量 [Wyr/m²/°C]
temp0 = 14.0 # 工业革命前温度 [°C]
def solar_only(t, temp, C, absorbed_solar):
return (1 / C) * absorbed_solar
sol1 = solve_ivp(
solar_only, [0, 170], [temp0],
args=(C, absorbed_solar_radiation),
dense_output=True, max_step=0.5
)
t_eval = np.linspace(0, 170, 500)
y_eval = sol1.sol(t_eval)[0]
fig, ax = plt.subplots(figsize=(8, 5))
ax.set_facecolor('black')
ax.plot(t_eval + 1850, y_eval, lw=2)
ax.axhline(temp0, color='white', ls='--')
ax.set_xlabel('年份（从 1850 年起）')
ax.set_ylabel('温度 °C')
ax.annotate(f'工业革命前温度 = {temp0}°C', xy=(1850 + 80, temp0 + 5), color='white')
ax.set_title('仅吸收太阳辐射')
plt.tight_layout()
plt.show()
热含量 $C \cdot \mathrm{temp}$ 由温度 $\mathrm{temp}$（以开尔文为单位）和气候系统的热容量决定。虽然大气温度才是关注的焦点（其热容量很小），但它与上层海洋的热量密切耦合，而上层海洋的热容量要大得多。
因此，热含量随时间的变化简单地由 $\frac{d(CT)}{dt}$ 给出。由于海水的热容量几乎不随温度变化，可以用温度随时间的变化来表示：
$$\text{热含量变化} = C \frac{d\mathrm{temp}}{dt}$$
向外热辐射¶
向外热辐射项 $\mathcal{G}(T)$（或"向太空的黑体冷却"）代表了抑制变暖的负反馈（如黑体辐射）和放大变暖的正反馈（如水汽反馈）的综合效应。
由于这些物理过程过于复杂，在此线性化模型，将入射和出射结合起来。
假设工业革命前的世界处于能量平衡状态，因此平衡温度就是工业革命前的温度。
于是假设
$$\mathrm{temp}'(t) = B(\mathrm{temp}(0) - \mathrm{temp}(t))$$
其中 $B$ 为某个值。$\mathrm{temp}(t)$ 前面的负号表示它恢复平衡。
所选取的值为
B = 1.3 # 气候反馈参数 [W/m²/°C]
定义健康地球的能量平衡模型，改变起始温度并绘图：
@interact(
start_temp=IntSlider(min=0, max=28, value=14, description='起始温度'),
)
def plot_energy_balance(start_temp):
def energy_balance(t, temp, C, B, temp0):
return (1 / C) * B * (temp0 - temp)
sol2 = solve_ivp(
energy_balance, [0, 170], [float(start_temp)],
args=(C, B, temp0),
dense_output=True, max_step=0.5
)
t_eval = np.linspace(0, 170, 500)
y_eval = sol2.sol(t_eval)[0]
fig, ax = plt.subplots(figsize=(8, 5))
ax.set_facecolor('black')
ax.plot(t_eval, y_eval, lw=2)
ax.axhline(temp0, color='white', ls='--')
ax.set_xlabel('从起始年份算起的年数')
ax.set_ylabel('温度 °C')
ax.set_ylim(0, 30)
ax.annotate(f'工业革命前温度 = {temp0}°C', xy=(80, temp0 - 1), color='white')
ax.set_title('能量平衡模型（健康的地球）')
plt.tight_layout()
plt.show()
interactive(children=(IntSlider(value=14, description='起始温度', max=28), Output()), _dom_classes=('widget-intera…
编程知识：嵌套函数定义
在上面的代码中，
energy_balance函数被定义在plot_energy_balance函数的内部。这样做的好处是：
- 内部函数可以直接访问外部函数的局部变量（如
start_temp），无需通过参数传递；- 内部函数只在被
@interact调用时才需要存在，不会污染全局命名空间。这种模式在交互式绘图（
@interact+solve_ivp）中非常常见。
温室效应¶
根据经验，温室效应已知是气态二氧化碳（CO₂）浓度的对数函数：
$$\text{人类引起的温室效应} = \text{强迫系数} \cdot \ln \left( \frac{[\text{CO}_2]}{[\text{CO}_2]_{\text{工业革命前}}} \right)$$
这如何随时间变化取决于人类行为！上述方程中并未对时间进行建模。
其中
forcing_coef = 5.0 # CO₂ 强迫系数 [W/m²]
CO2_PreIndust = 280.0 # 工业革命前 CO₂ 浓度 [百万分之一; ppm]
def greenhouse_effect(CO2):
return forcing_coef * np.log(CO2 / CO2_PreIndust)
CO2_present = 420.0
CO2_range = 280 * (2 ** np.linspace(-1, 3, 100))
fig, ax = plt.subplots(figsize=(6, 4))
ax.plot(CO2_range, greenhouse_effect(CO2_range), lw=2.5, color='black', label=None)
ax.plot(CO2_PreIndust, greenhouse_effect(CO2_PreIndust), 'o', ms=6, color='blue', label='工业革命前 (PI)')
ax.plot(CO2_present, greenhouse_effect(CO2_present), 'o', ms=6, color='red', label='当前 (2020)')
ax.set_xticks([280, 280*2, 280*4, 280*8])
ax.set_xticklabels(['280', '560', '1120', '2240'])
ax.legend(loc='lower right')
ax.set_ylabel('辐射强迫 [W/m²]')
ax.set_xlabel('CO₂ 浓度 [ppm]')
ax.set_title('温室效应 vs CO₂ 浓度')
plt.tight_layout()
plt.show()
def CO2(t):
return CO2_PreIndust * (1 + (t / 220) ** 3)
years = np.arange(1850, 2031)
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(years, CO2(years - 1850), lw=3) # 向量化运算：years 是数组，整个数组传入函数无需循环
ax.set_xlabel('年份')
ax.set_ylabel('CO₂ (ppm)')
ax.set_title('CO₂ 浓度模型（三次多项式拟合）')
plt.tight_layout()
plt.show()
def climate_model_with_co2(t, temp, C, B, temp0):
return (1 / C) * (B * (temp0 - temp) + greenhouse_effect(CO2(t)))
start_temp_val = 14.0
sol3 = solve_ivp(
climate_model_with_co2, [0, 170], [start_temp_val],
args=(C, B, temp0),
dense_output=True, max_step=0.5
)
绘制包含 CO₂ 排放的温度随时间变化：
t_eval = np.linspace(0, 170, 500)
y_eval = sol3.sol(t_eval)[0]
fig, ax = plt.subplots(figsize=(8, 5))
ax.set_facecolor('black')
ax.plot(t_eval + 1850, y_eval, lw=2)
ax.axhline(temp0, color='white', ls='--')
ax.set_xlabel('年份（从 1850 年起）')
ax.set_ylabel('温度 °C')
ax.set_ylim(10, 20)
ax.annotate(f'工业革命前温度 = {temp0}°C', xy=(1850 + 80, temp0 - 0.5), color='white')
ax.set_title('包含 CO₂ 的模型')
plt.tight_layout()
plt.show()
CO2_historical_data_url = "https://gml.noaa.gov/webdata/ccgg/trends/co2/co2_mm_mlo.csv"
CO2_historical_data_raw = pd.read_csv(
CO2_historical_data_url, comment='#', header=0 # comment='#' 跳过以 # 开头的注释行
)
CO2_historical_data_raw.head(10)
| year | month | decimal date | average | deseasonalized | ndays | sdev | unc | |
|---|---|---|---|---|---|---|---|---|
| 0 | 1958 | 3 | 1958.2027 | 315.71 | 314.44 | -1 | -9.99 | -0.99 |
| 1 | 1958 | 4 | 1958.2877 | 317.45 | 315.16 | -1 | -9.99 | -0.99 |
| 2 | 1958 | 5 | 1958.3699 | 317.51 | 314.69 | -1 | -9.99 | -0.99 |
| 3 | 1958 | 6 | 1958.4548 | 317.27 | 315.15 | -1 | -9.99 | -0.99 |
| 4 | 1958 | 7 | 1958.5370 | 315.87 | 315.20 | -1 | -9.99 | -0.99 |
| 5 | 1958 | 8 | 1958.6219 | 314.93 | 316.21 | -1 | -9.99 | -0.99 |
| 6 | 1958 | 9 | 1958.7068 | 313.21 | 316.11 | -1 | -9.99 | -0.99 |
| 7 | 1958 | 10 | 1958.7890 | 312.42 | 315.41 | -1 | -9.99 | -0.99 |
| 8 | 1958 | 11 | 1958.8740 | 313.33 | 315.21 | -1 | -9.99 | -0.99 |
| 9 | 1958 | 12 | 1958.9562 | 314.67 | 315.43 | -1 | -9.99 | -0.99 |
数据在 "average" 列中。
注意：有缺失数据（-99.99），需要过滤。
valid_mask = CO2_historical_data_raw['average'] > 0
CO2_historical_data = CO2_historical_data_raw[valid_mask].copy() # .copy() 创建独立副本，避免 SettingWithCopyWarning
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(
CO2_historical_data['decimal date'],
CO2_historical_data['average'],
label='莫纳罗亚 CO₂ 数据（基林曲线）'
)
ax.plot(years, CO2(years - 1850), lw=3, label='三次拟合')
ax.legend(loc='upper left')
ax.set_title('CO₂ 观测数据与拟合')
ax.set_xlabel('年份')
ax.set_ylabel('CO₂ (ppm)')
plt.tight_layout()
plt.show()
交互模型：调整气候参数¶
考虑使气候反馈参数 $B$ 和海洋热容参数 $C$ 可调，分别用 $BB$ 和 $CC$ 表示，并加载 NASA 温度观测数据进行对比。
@interact(
BB=FloatSlider(min=0, max=4, step=0.1, value=1.3, description='气候反馈 BB'),
CC=FloatSlider(min=10, max=200, step=0.1, value=51, description='海洋热容量 CC'),
)
def plot_with_params(BB, CC):
def model4(t, temp, BB, CC, temp0):
return (1 / CC) * (BB * (temp0 - temp) + greenhouse_effect(CO2(t)))
sol4 = solve_ivp(
model4, [0, 170], [14.0],
args=(BB, CC, temp0),
dense_output=True, max_step=0.5
)
years_plot = np.arange(1850, 2031)
y_model = sol4.sol(years_plot - 1850)[0]
# 下载 NASA 温度观测数据
T_url = "https://data.giss.nasa.gov/gistemp/graphs_v4/graph_data/Global_Mean_Estimates_based_on_Land_and_Ocean_Data/graph.txt"
try:
T_df = pd.read_csv(T_url, sep=r'\s+', skiprows=5, header=None)
T_years = T_df.iloc[:, 0].astype(float).values
T_anomaly = T_df.iloc[:, 1].astype(float).values + 14.15
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(years_plot, y_model, lw=2, label='模型预测温度')
ax.plot(T_years, T_anomaly, color='black', label='NASA 观测数据')
ax.legend(loc='upper left')
ax.set_xlabel('年份')
ax.set_ylabel('温度 °C')
ax.set_title(f'气候模型（BB={BB:.1f}, CC={CC:.1f}）')
except Exception as e:
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(years_plot, y_model, lw=2, label='模型预测温度')
ax.legend(loc='upper left')
ax.set_xlabel('年份')
ax.set_ylabel('温度 °C')
ax.set_title(f'气候模型（BB={BB:.1f}, CC={CC:.1f}）- 无法加载 NASA 数据')
plt.tight_layout()
plt.show()
interactive(children=(FloatSlider(value=1.3, description='气候反馈 BB', max=4.0), FloatSlider(value=51.0, descript…
未来气候预测：最好和最坏情景¶
考虑两种截然不同的假设未来：
- 低排放世界：排放量减少，使得到 2100 年 CO₂ 浓度保持在 500 ppm 以下（在气候界称为"RCP2.6"）
- 高排放世界：排放量持续增加，CO₂ 浓度飙升至 1200 ppm 以上（"RCP8.5"）
T0 = 14.0
C_rcp = 100.0
B_rcp = 0.3
alpha_rcp = 1.5
def get_co2_concentration(year, scenario='high'):
if year <= 2020:
return 280 * np.exp(0.0023 * (year - 1850))
else:
t_future = year - 2020
if scenario == 'high':
return 415 * np.exp(0.0133 * t_future)
else:
return 415 + (t_future * 0.5) - (t_future**2 * 0.001)
def temperature_ode(t, T, co2_func):
co2 = co2_func(t)
forcing = alpha_rcp * np.log(co2 / 280)
dT_dt = (1 / C_rcp) * (forcing - B_rcp * (T - T0))
return dT_dt
years_rcp = np.linspace(1850, 2100, 250)
co2_high = [get_co2_concentration(y, 'high') for y in years_rcp]
co2_low = [get_co2_concentration(y, 'low') for y in years_rcp]
sol_high = solve_ivp(temperature_ode, [1850, 2100], [T0], args=(lambda t: get_co2_concentration(t, 'high'),), t_eval=years_rcp)
sol_low = solve_ivp(temperature_ode, [1850, 2100], [T0], args=(lambda t: get_co2_concentration(t, 'low'),), t_eval=years_rcp)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 6))
ax1.plot(years_rcp, co2_low, color='#00aaff', lw=2.5, label='Low emissions')
ax1.plot(years_rcp, co2_high, color='#e57347', lw=2.5, label='High emissions')
ax1.scatter([2020], [415], color='black', s=80, zorder=5, label='You Are Here')
ax1.set_title("Cause...", fontsize=20)
ax1.set_ylabel("CO₂ concentration [ppm]", fontsize=14)
ax1.set_xlabel("year", fontsize=14)
ax1.legend()
ax1.grid(alpha=0.2)
ax2.plot(years_rcp, sol_low.y[0], color='#00aaff', lw=2.5)
ax2.plot(years_rcp, sol_high.y[0], color='#e57347', lw=2.5)
ax2.axhline(y=16, color='gray', linestyle='--', label='Paris Agreement threshold (2°C warming)')
ax2.set_title("... and effect", fontsize=20)
ax2.set_ylabel("temperature [°C]", fontsize=14)
ax2.set_xlabel("year", fontsize=14)
ax2.legend()
ax2.grid(alpha=0.2)
plt.tight_layout()
plt.show()
在低排放情景下，到 2100 年温度增幅保持在 $\Delta T = 2$ °C 以下，而在高排放情景下，温度比工业革命前水平上升超过 3.5 °C。
虽然人类 CO₂ 排放引起的温室效应是历史和未来预计变暖的主要强迫因素，但现代气候建模考虑了相当详尽的其他强迫因素（气溶胶、其他温室气体、臭氧、土地利用变化等）。下面的视频展示了最先进的气候模型在历史时期模拟中这些强迫因素的分解。
视频：是什么让地球变暖？¶
%%HTML
<video width="960" controls>
<source src="https://pacman.cs.tsinghua.edu.cn/~hanwentao/cpct/19/whats-really-warming-the-world.mp4" type="video/mp4">
</video>
本讲小结¶
本讲从计算思维的四个核心要素出发，构建了一个简单的气候模型：
1. 抽象（Abstraction）
复杂的地球气候系统被抽象为一个零维能量平衡模型：忽略经纬度、海拔等空间细节，只关注全球平均温度。三个物理过程（太阳辐射、热辐射、温室效应）被抽象为微分方程中的三个数学项，对应三种 ODE 形式（$y'=a$、$y'=a-by$、$y'=a-by+f(t)$）。
2. 分解（Decomposition）
采用逐步叠加的策略分解问题：先只考虑太阳辐射（温度无限上升），再加入热辐射（恢复平衡），最后加入温室效应（温度再次上升）。每一步都有清晰的物理意义和对应的数学形式，最终组合成完整模型。
3. 建模与仿真（Modeling & Simulation）
物理规律被翻译为 Python 函数（solar_only、energy_balance、climate_model_with_co2），使用 scipy.integrate.solve_ivp 进行数值求解，实现了用代码模拟真实世界。通过调整参数（$B$、$C$），用模型与 NASA 观测数据对比，验证了模型的合理性。
4. 数据驱动的分析（Data-Driven Analysis）
从 NOAA 在线获取莫纳罗亚火山的真实 CO₂ 观测数据，经过清洗（过滤缺失值）、可视化后与模型拟合对比。最后基于 RCP2.6 和 RCP8.5 两种情景进行未来预测，从"因"（CO₂ 浓度）到"果"（温度变化）给出了定量分析。