问题引入¶
上一讲用零维能量平衡模型描述地球温度，发现它只有一个稳定平衡态——前工业温度 $T = 14$°C。在前工业时期，系统总会回到那个唯一的平衡态。
然而，地质证据告诉我们，地球并非一直如此温和。大约7亿年前，地球经历过两次雪球地球事件——从赤道到两极，整个星球被冰雪覆盖。如果模型只有一个平衡态，这是不可能的。
更深入的问题是：
- 一个系统如何能拥有多个平衡态？
- 系统最终处于哪个状态，与什么因素有关？
- 能否在上一讲的线性模型基础上进行修改，来解释这一现象？
本讲将先从一个简化的数学模型出发，理解迟滞（hysteresis）现象的本质，再将这一概念应用到气候模式中——只需将反照率从常数改为温度的函数，就能让"雪球地球"和"温暖地球"共存。
准备工作¶
In [1]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from scipy.optimize import brentq
from ipywidgets import interact, FloatSlider
plt.rcParams['font.family'] = ['Noto Sans', 'Noto Sans CJK SC']
数学概念：sign函数¶
sign(x) 在 x=0 时返回 0，在 x 正/负时返回 ±1。
In [2]:
fig, ax = plt.subplots(figsize=(6, 3))
#xs = np.arange(-5, 5.01, 0.1)
xs = np.linspace(-5, 5, 101)
ax.scatter(xs, np.sign(xs), s=10, c='blue', alpha=0.7)
ax.set_xlabel('x')
ax.set_ylabel('sign(x)')
ax.set_title('sign函数在0处不连续')
ax.axhline(0, color='gray', ls='--', alpha=0.5)
plt.tight_layout()
plt.show()
print(f'sign(Inf) = {np.sign(np.inf)}')
print(f'sign(-Inf) = {np.sign(-np.inf)}')
sign(Inf) = 1.0 sign(-Inf) = -1.0
数学：多重平衡态¶
函数 $f(y,a) = \mathrm{sign}(y) + a - y$ 可以写为
$$f(y,a)= \begin{cases} -1+a-y & \text{if } y<0 \quad (\text{根在 } a-1 \text{, 若 } a<1 ) \\ 1+a-y & \text{if } y>0 \quad (\text{根在 } a+1 \text{, 若 } a>-1 ) \end{cases}$$
（忽略 $y=0$）。注意当 $-1<a<1$ 时有两个根。
In [3]:
def f(y, a):
return np.sign(y) + a - y
def f_ode(t, y, a):
return f(y, a)
函数 $z=f(y,a)$ 的图像由两个平行半平面组成。下图左侧是常数 $a$ 的截面，右侧是与 $z=0$ 的交线。
In [4]:
def plot_f(a=0.0):
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
eps = 1e-10
y_neg = np.linspace(-9, -eps, 500)
y_pos = np.linspace(eps, 9, 500)
ax1.plot(y_neg, [f(y, a) for y in y_neg], lw=3, color='blue')
ax1.plot(y_pos, [f(y, a) for y in y_pos], lw=3, color='blue')
ax1.axhline(0, ls='--', color='pink')
if a < 1:
root1 = -1 + a
ax1.scatter([root1], [0], s=100, c='red', zorder=5)
ax1.annotate(f'{root1:.1f}', (root1, -4), color='red', ha='center')
if a > -1:
root2 = 1 + a
ax1.scatter([root2], [0], s=100, c='blue', zorder=5)
ax1.annotate(f'{root2:.1f}', (root2, -4), color='blue', ha='center')
ax1.axvline(0, ls='--', color='pink', alpha=0.5)
ax1.set_xlim(-9, 9)
ax1.set_ylim(-5, 10)
ax1.set_xlabel('y')
ax1.set_ylabel('sign(y) - y + a')
ax1.set_title(f'f(y, a={a:.1f}) vs y')
a_vals_neg = np.linspace(-5, 1, 200)
a_vals_pos = np.linspace(-1, 5, 200)
ax2.plot(a_vals_neg, a_vals_neg - 1, lw=2, c='blue', label='y < 0 分支')
ax2.plot(a_vals_pos, a_vals_pos + 1, lw=2, c='blue', label='y > 0 分支')
if a < 1:
ax2.scatter([a], [-1+a], s=100, c='red', zorder=5)
ax2.axhline(-1+a, c='gray', ls='--', alpha=0.5)
if a > -1:
ax2.scatter([a], [1+a], s=100, c='blue', zorder=5)
ax2.axhline(1+a, c='gray', ls='--', alpha=0.5)
ax2.axvline(a, c='gray', ls='--', alpha=0.5)
ax2.set_xlim(-5, 5)
ax2.set_ylim(-5, 5)
ax2.set_xlabel('a')
ax2.set_ylabel('y')
ax2.set_title('f(y,a)=0 的解')
ax2.annotate('当 -1<a<1 时有两个根', (0, -4.7), color='blue', ha='center')
plt.tight_layout()
plt.show()
plot_f()
In [5]:
interact(plot_f, a=FloatSlider(min=-5, max=5, step=0.1, value=0.0, description='a'));
interactive(children=(FloatSlider(value=0.0, description='a', max=5.0, min=-5.0), Output()), _dom_classes=('wi…
y' = f(y,a) 的解¶
初始条件 $y(0)=y_0$
In [6]:
def plot_ode(y0=2.0, a=0.0):
sol = solve_ivp(f_ode, (0, 10), [y0], args=(a,), max_step=0.01, dense_output=True)
fig, ax = plt.subplots(figsize=(10, 5))
ax.set_facecolor('black')
t_plot = np.linspace(0, 10, 1000)
y_plot = sol.sol(t_plot)[0]
ax.plot(t_plot, y_plot, lw=3, c='red')
xs_arr = []
ys_arr = []
for x in np.linspace(0, 10, 25):
for y in np.linspace(-7, 7, 25):
v = np.array([1, f(y, a)])
v = v / (100 * 0.4)
xs_arr.extend([x - v[0], x + v[0], np.nan])
ys_arr.extend([y - v[1], y + v[1], np.nan])
ax.plot(xs_arr, ys_arr, alpha=0.5, c='yellow', lw=0.8)
if a < 1:
ax.axhline(-1 + a, c='white', ls='--')
if a > -1:
ax.axhline(1 + a, c='white', ls='--')
ax.axhline(0, c='pink', ls='--')
ax.set_ylim(-7, 7)
ax.set_xlabel('t')
ax.set_ylabel('y')
ax.set_title("y' = f(y,a) 的解")
plt.tight_layout()
plt.show()
plot_ode()
In [7]:
interact(plot_ode,
y0=FloatSlider(min=-6, max=6, step=0.1, value=2.0, description='y₀'),
a=FloatSlider(min=-5, max=5, step=0.1, value=0.0, description='a'));
interactive(children=(FloatSlider(value=2.0, description='y₀', max=6.0, min=-6.0), FloatSlider(value=0.0, desc…
迟滞：先增大再减小 a¶
将 $a$ 从 -4 增加到 4（步长 0.25），然后再从 4 减小到 -4。每次改变 $a$ 后，让系统演化 10 个时间单位（足以达到该 $a$ 值的平衡态），并观察 $y$ 值。
可以看到当 $-1<a<1$ 时，根据到达方式的不同，系统可能处于“负”平衡态或“正”平衡态。
In [8]:
dt = 10.0
t = 0.0
y0 = -5.0
a_val = -4.0
sol = solve_ivp(f_ode, (t, t + dt), [y0], args=(a_val,), max_step=0.1, dense_output=True)
fig, ax = plt.subplots(figsize=(14, 6))
ax.plot(np.linspace(t, t + dt, 100), sol.sol(np.linspace(t, t + dt, 100))[0])
a_annotations_forward = []
for i in range(64):
# 先增大， 再减小
a_val += 0.25 if i < 32 else -0.25
t += dt
y0 = sol.sol(t)[0]
sol = solve_ivp(f_ode, (t, t + dt), [y0], args=(a_val,), max_step=0.1, dense_output=True)
t_arr = np.linspace(t, t + dt, 100)
y_arr = sol.sol(t_arr)[0]
if -1 <= a_val <= 1:
ax.fill_between(t_arr, y_arr, -60, color='gray', alpha=0.3)
ax.plot(t_arr, y_arr, c='blue' if i < 32 else 'red')
if (i + 1) % 4 == 0:
ax.annotate(f'{a_val:.1f}', (t, -5.5), color='blue' if i < 32 else 'red', ha='center')
ax.set_ylim(-6, 6)
ax.set_ylabel('y')
ax.set_xlabel('t')
ax.set_title('迟滞：先增大a（蓝色）再减小a（红色）')
plt.show()
系统状态对其历史的依赖性，即上文观察到的现象，被称为迟滞（hysteresis，源自希腊语 ὑστέρησις，意为“迟滞”）。
应用：雪球地球与冰-反照率反馈¶
回顾上一讲气候内容¶
回顾上一讲的内容，零维能量平衡方程为
$$C \frac{dT}{dt} = \frac{(1 - \alpha)S}{4} - (A - BT) + a \ln \left( \frac{[\mathrm{CO}_2]}{[\mathrm{CO}_2]_{\mathrm{PI}}} \right)$$
现在忽略CO₂的变化，因此模型简化为
$$C \frac{dT}{dt} = \frac{(1 - \alpha)S}{4} - (A - BT)$$
这个常微分方程的动力学非常简单，因为它是线性的。线性ODE只允许一个稳定解。
本讲将展示如何通过一个小的修改使模型变为非线性的，从而完全改变其动力学，使“雪球地球”和相对温暖的前工业气候的同时存在得以解释。
1.1) 冰-反照率反馈¶
在第20讲中，使用了常数值 $\alpha = 0.3$ 作为地球的行星反照率，对于相对于现在的较小气候变异（如当代与前工业气候之间的差异），这是合理的。然而，在大的变化下，这种近似就不太可靠了。
海洋是深色的、吸热的，$\alpha_{\mathrm{ocean}} \approx 0.05$，而冰雪是明亮的、反射性的：$\alpha_{\mathrm{ice,snow}} \approx 0.5$ 到 $0.9$。因此，如果大部分海洋表面结冰，地球反照率会急剧上升，导致更多阳光被反射回太空，进而导致更多冷却和更多海洋结冰，等等。这种非线性正反馈效应被称为冰-反照率反馈。
可以通过允许反照率依赖于温度来粗略地表示冰-反照率反馈：
$$\alpha(T) = \begin{cases} \alpha_i & \text{if } T \leq -10\text{°C} & \text{(完全冻结)} \\ \alpha_i + (\alpha_0 - \alpha_i)\frac{T + 10}{20} & \text{if } -10\text{°C} \leq T \leq 10\text{°C} & \text{(部分冻结)} \\ \alpha_0 & \text{if } T \geq 10\text{°C} & \text{(无冰)} \end{cases}$$
1.2) 将冰-反照率反馈加入简单气候模型¶
In [9]:
# 能量平衡模型 (EBM) 参数
S_const = 1368.0 # 太阳辐照 [W/m^2]
alpha_0 = 0.3 # 默认反照率（行星反射率）
B_const = -1.3 # 气候反馈参数 [W/m^2/°C]
T0_const = 14.0 # 前工业温度 [°C]
A_const = S_const * (1.0 - alpha_0) / 4.0 + B_const * T0_const # [W/m^2]
a_const = 5.0 # CO2强迫系数 [W/m^2]
CO2_PI = 280.0 # 前工业CO2浓度 [ppm]
C_const = 51.0 # 大气和上层海洋热容量 [J/m^2/°C]
alpha_i = 0.5 # 冰面反照率
DT_threshold = 10.0 # 温度阈值 [°C]
def albedo(T):
if T < -DT_threshold:
return alpha_i
elif T < DT_threshold:
return alpha_i + (alpha_0 - alpha_i) * (T + DT_threshold) / (2 * DT_threshold)
else:
return alpha_0
def absorbed_solar_radiation(T, S=S_const):
return S * (1.0 - albedo(T)) / 4.0
def outgoing_thermal_radiation(T, A=A_const, B=B_const):
return A - B * T
def tendency(T, S=S_const):
return (1.0 / C_const) * (
absorbed_solar_radiation(T, S) - outgoing_thermal_radiation(T)
)
def ebm_ode(t, state, S):
T = state[0]
return [tendency(T, S)]
def run_ebm(T0, dt, end_year, S=S_const):
ts = [0.0]
Ts = [T0]
t = 0.0
T = T0
while t < end_year:
dT = tendency(T, S) * dt
T = T + dT
t = t + dt
Ts.append(T)
ts.append(t)
return np.array(ts), np.array(Ts)
print(f'前工业平衡温度: T0 = {T0_const} °C')
print(f'常数 A = {A_const:.2f} W/m^2')
print(f'冰面反照率 α_i = {alpha_i}')
print(f'水面反照率 α_0 = {alpha_0}')
前工业平衡温度: T0 = 14.0 °C 常数 A = 221.20 W/m^2 冰面反照率 α_i = 0.5 水面反照率 α_0 = 0.3
接下来绘制反照率随温度变化的函数：
In [10]:
T_example = np.arange(-20, 20.5, 1.0)
alpha_example = [albedo(T) for T in T_example]
fig, ax = plt.subplots(figsize=(6, 4))
ax.plot(T_example, alpha_example, lw=3, color='black', label='α(T)')
ax.axvspan(-20, -10, alpha=0.15, color='lightblue')
ax.axvspan(-10, 10, alpha=0.08, color='gray')
ax.axvspan(10, 20, alpha=0.10, color='red')
ax.annotate('完全冻结', (-15.5, 0.252), color='darkblue')
ax.annotate('部分冻结', (-2.5, 0.252), color='darkgrey')
ax.annotate('无冰', (13, 0.252), color='darkred')
ax.set_ylabel('反照率 α (行星反射率)')
ax.set_xlabel('温度 [°C]')
ax.set_ylim(0.2, 0.6)
plt.tight_layout()
plt.show()
In [11]:
fig, ax = plt.subplots(figsize=(8, 6))
ax.set_xlabel('年')
ax.set_ylabel('温度 [°C]')
ax.set_xlim(0, 205)
ax.set_ylim(-60, 30)
ax.axhspan(-60, -10, alpha=0.2, color='lightblue')
ax.annotate('完全冻结', (120, -20), color='darkblue')
ax.axhspan(10, 30, alpha=0.06, color='red')
ax.annotate('无冰', (120, 25), color='darkred')
for T0_sample in np.arange(-60, 30.5, 5):
ts, Ts = run_ebm(T0_sample, 1.0, 200)
ax.plot(ts, Ts)
# 不稳定平衡附近
T_un = 7.5472
for dT in [-0.02, -0.01, 0.0, 0.01, 0.02]:
ts, Ts = run_ebm(T_un + dT, 1.0, 200)
ax.plot(ts, Ts, ls='--')
ax.scatter([200], [T0_const], s=100, c='orange', marker='o', zorder=5,
label='前工业气候（稳定“暖”分支）')
ax.scatter([200], [-38.3], s=100, c='aqua', marker='o', zorder=5,
label='另一种宇宙前工业气候（稳定“冷”分支）')
ax.scatter([200], [T_un], s=100, c='lightgrey', marker='D', edgecolors='white',
zorder=5, label='不可能的另一种气候（不稳定分支）')
ax.legend(loc='lower right')
plt.tight_layout()
plt.show()
下面通过滑块改变初始温度 $T_0$，观察不同初始条件下的演化行为：
In [12]:
def plot_ebm_interactive(T0=24.0):
ts, Ts = run_ebm(T0, 1.0, 200)
fig, ax = plt.subplots(figsize=(8, 6))
ax.set_xlabel('年')
ax.set_ylabel('温度 [°C]')
ax.set_xlim(0, 205)
ax.set_ylim(-60, 30)
ax.axhspan(-60, -10, alpha=0.2, color='lightblue')
ax.annotate('完全冻结', (120, -20), color='darkblue')
ax.axhspan(10, 30, alpha=0.06, color='red')
ax.annotate('无冰', (120, 25), color='darkred')
# 绘制所有样本轨迹（灰色）
for T0_sample in np.arange(-60, 30.5, 5):
ts_s, Ts_s = run_ebm(T0_sample, 1.0, 200)
ax.plot(ts_s, Ts_s, color='gray', alpha=0.3)
# 不稳定平衡附近
T_un = 7.5472
for dT in [-0.02, -0.01, 0.0, 0.01, 0.02]:
ts_u, Ts_u = run_ebm(T_un + dT, 1.0, 200)
ax.plot(ts_u, Ts_u, ls='--', color='gray', alpha=0.3)
# 高亮当前选择的轨迹
ax.plot(ts, Ts, lw=3, color='green', zorder=4, label=f'T₀ = {T0:.1f}°C')
ax.scatter([200], [T0_const], s=100, c='orange', marker='o', zorder=5,
label='前工业气候（稳定暖分支）')
ax.scatter([200], [-38.3], s=100, c='aqua', marker='o', zorder=5,
label='另一种气候（稳定冷分支）')
ax.legend(loc='lower right')
plt.tight_layout()
plt.show()
interact(plot_ebm_interactive,
T0=FloatSlider(min=-60, max=30, step=0.5, value=24.0, description='T₀ [°C]'));
interactive(children=(FloatSlider(value=24.0, description='T₀ [°C]', max=30.0, min=-60.0, step=0.5), Output())…
可以看到当 $T_0 \gtrsim 7.55$ °C 时，所有曲线似乎都收敛到第20讲中看到的 $T = 14$°C 的平衡态（或不动点）。从该值以下出发的曲线升温，而从该值以上出发的曲线降温。然而，当 $T_0 \lesssim 7.55$ °C 时，温度收敛到一个更冷的平衡态，约为 $T = -40$°C。这就是雪球地球平衡态。这两个状态被称为稳定平衡态，因为即使状态暂时略微偏离平衡，它最终也会收敛回平衡。
那么当 $T_0 \approx 7.55$ °C 时会发生什么？在那个精确温度附近确实存在一个平衡态。然而，如果温度从该精确值偏离哪怕百分之一度，温度最终会收敛到其他两个平衡态之一。因此，这个中间平衡态被称为不稳定平衡态，因为任何无穷小的偏离都会使其偏离到另一个状态。
2.2) 辐射稳定性分析¶
可以通过应用动力系统理论的概念来理解为什么模型有两个稳定平衡态和一个不稳定平衡态。
在固定CO₂浓度的情况下，能量平衡模型微分方程可以表示为：
$$C\frac{dT}{dt} = \mathrm{ASR}(T) - \mathrm{OTR}(T),$$
现在吸收太阳辐射（ASR）也是温度依赖的，因为反照率 $\alpha(T)$ 是温度的函数。
通过绘制右侧趋势项作为状态变量 $T$ 的函数，可以得到系统的稳定性图，判断行星会升温（$C\frac{dT}{dt} > 0$）还是降温（$C\frac{dT}{dt} < 0$）。
In [13]:
T_samples = np.arange(-60, 30.5, 1.0)
OTR = [outgoing_thermal_radiation(T) for T in T_samples]
ASR = [absorbed_solar_radiation(T) for T in T_samples]
imbalance = np.array(ASR) - np.array(OTR)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
# 左图: ASR 和 OTR
ax1.plot(T_samples, OTR, label='向外热辐射 (OTR)', color='blue', lw=2)
ax1.plot(T_samples, ASR, label='吸收太阳辐射 (ASR)', color='orange', lw=2)
ax1.set_ylabel('能量通量 [W/m²]')
ax1.set_xlabel('温度 [°C]')
ax1.legend(loc='upper left')
# 右图: 辐射不平衡
ax2.fill_between([-60, 30], [-100, -100], [0, 0], color='blue', alpha=0.08)
ax2.fill_between([-60, 30], [0, 0], [100, 100], color='red', alpha=0.08)
ax2.annotate('降温', (-58, -40), color='darkblue')
ax2.annotate('升温', (-58, 38), color='darkred')
ax2.plot(T_samples, imbalance, label='辐射不平衡 (ASR - OTR)', color='black', lw=2)
ax2.scatter([T_un], [0], s=80, c='lightgrey', marker='D', edgecolors='darkgrey', zorder=5)
ax2.scatter([T0_const], [0], s=80, c='orange', edgecolors='black', zorder=5)
ax2.scatter([-38.3], [0], s=80, c='aqua', edgecolors='black', zorder=5)
ax2.set_ylim(-50, 45)
ax2.set_ylabel('能量通量 [W/m²]')
ax2.set_xlabel('温度 [°C]')
ax2.legend(loc='upper left')
plt.tight_layout()
plt.show()
3.2) 太阳变亮是否融化了雪球？¶
如果从新元古代气候出发，仅仅将太阳辐照增加到今天的值 $S = 1368$ W/m²，能让行星升温到前工业温度 $T=14$°C吗？
下面绘制太阳辐照的分岔图，观察温度如何随 $S$ 变化：
In [14]:
Smin = 1200
Smax = 1800
Smax_limited = 1650
# 计算分岔图：先增大S，再减小S
Svec_up = np.arange(Smin, Smax + 1, 1.0)
Svec_down = Svec_up[::-1]
Svec = np.concatenate([Svec_up, Svec_down])
Tvec = np.zeros_like(Svec)
T_restart = -100.0
for i, S_val in enumerate(Svec):
ts, Ts = run_ebm(T_restart, 5.0, 400.0, S=S_val)
T_restart = Ts[-1]
Tvec[i] = T_restart
n_half = len(Svec_up)
Svec_warming = Svec[:n_half]
Tvec_warming = Tvec[:n_half]
Svec_cooling = Svec[n_half:]
Tvec_cooling = Tvec[n_half:]
print(f'分岔图数据计算完成')
print(f'S范围: {Smin} - {Smax} W/m²')
print(f'前工业 S = {S_const} W/m²')
分岔图数据计算完成 S范围: 1200 - 1800 W/m² 前工业 S = 1368.0 W/m²
In [15]:
# 不稳定分支
def T_unstable_branch(S_val, A, B, ai, a0):
return (A/B - S_val/(4*B) * (1 - 0.5*(ai + a0))) / (1 + S_val*(ai - a0)/(80*B))
S_unstable = np.arange(Smin, Smax + 1, 2.0)
T_unstable = [T_unstable_branch(S, A_const, B_const, 0.5, 0.3) for S in S_unstable]
# 只保留反照率在过渡区域内的部分
mask = [0.3 < albedo(T) < 0.5 for T in T_unstable]
S_unstable_plot = np.where(mask, S_unstable, np.nan)
T_unstable_plot = np.where(mask, T_unstable, np.nan)
def plot_bifurcation(S=1368.0):
fig, ax = plt.subplots(figsize=(10, 7))
ax.axhspan(-60, -10, alpha=0.2, color='lightblue')
ax.annotate('完全冻结', (Smin + 12, -19), color='darkblue')
ax.axhspan(10, 80, alpha=0.06, color='red')
ax.annotate('无冰', (Smin + 12, 15), color='darkred')
ax.plot(Svec_warming, Tvec_warming, lw=3, color='blue', alpha=0.6, label='冷分支（S增大）')
ax.plot(Svec_cooling, Tvec_cooling, lw=3, color='red', alpha=0.6, label='暖分支（S减小）')
ax.plot(S_unstable_plot, T_unstable_plot, lw=3, color='gray', ls='--', alpha=0.5, label='不稳定分支')
ax.scatter([S_const], [T0_const], s=120, c='orange', marker='o', edgecolors='black',
zorder=5, label='前工业气候')
ax.scatter([S_const], [-38.3], s=120, c='aqua', marker='o', edgecolors='black',
zorder=5, label='另一种前工业气候')
Sneo = S_const * 0.93
ax.scatter([Sneo], [-48], s=120, c='lightblue', marker='o', edgecolors='black',
zorder=5, label='新元古代（7亿年前）')
# 标记当前S值
ax.axvline(S, color='green', lw=2, ls='-', alpha=0.8, label=f'S = {S:.0f} W/m²')
ax.set_xlim(Smin, Smax_limited)
ax.set_ylim(-55, 75)
ax.set_xlabel('太阳辐照 S [W/m²]')
ax.set_ylabel('全球温度 T [°C]')
ax.set_title('地球太阳辐照分岔图')
ax.legend(loc='upper left')
plt.tight_layout()
plt.show()
plot_bifurcation()
下面通过滑块改变太阳辐照 $S$，在分岔图上标记对应位置，并观察该辐照下的气候状态：
In [16]:
interact(plot_bifurcation,
S=FloatSlider(min=1200, max=1650, step=2, value=1368, description='S [W/m²]'));
interactive(children=(FloatSlider(value=1368.0, description='S [W/m²]', max=1650.0, min=1200.0, step=2.0), Out…
对于太阳辐照 $S$ 在 1200 W/m² 到 1650 W/m² 之间的值，温度始终保持在 -10°C 以下，行星保持完全冻结。如果扩展太阳辐照的上限，使太阳变得足够亮以开始融化冰呢？
In [17]:
def plot_bifurcation_extended(S=1368.0):
fig, ax = plt.subplots(figsize=(10, 7))
ax.axhspan(-60, -10, alpha=0.2, color='lightblue')
ax.annotate('完全冻结', (Smin + 12, -19), color='darkblue')
ax.axhspan(10, 80, alpha=0.06, color='red')
ax.annotate('无冰', (Smin + 12, 15), color='darkred')
ax.plot(Svec_warming, Tvec_warming, lw=3, color='blue', alpha=0.6, label='冷分支（S增大）')
ax.plot(Svec_cooling, Tvec_cooling, lw=3, color='red', alpha=0.6, label='暖分支（S减小）')
ax.plot(S_unstable_plot, T_unstable_plot, lw=3, color='gray', ls='--', alpha=0.5, label='不稳定分支')
ax.scatter([S_const], [T0_const], s=120, c='orange', marker='o', edgecolors='black',
zorder=5, label='前工业气候')
ax.scatter([S_const], [-38.3], s=120, c='aqua', marker='o', edgecolors='black',
zorder=5, label='另一种前工业气候')
Sneo = S_const * 0.93
ax.scatter([Sneo], [-48], s=120, c='lightblue', marker='o', edgecolors='black',
zorder=5, label='新元古代（7亿年前）')
ax.axvline(S_const, color='yellow', alpha=0.3, lw=8, label='前工业/当代辐照')
# 标记当前S值
ax.axvline(S, color='green', lw=2, ls='-', alpha=0.8, label=f'S = {S:.0f} W/m²')
ax.set_xlim(Smin, Smax)
ax.set_ylim(-55, 75)
ax.set_xlabel('太阳辐照 S [W/m²]')
ax.set_ylabel('全球温度 T [°C]')
ax.set_title('地球太阳辐照分岔图（扩展范围）')
ax.legend(loc='upper left')
plt.tight_layout()
plt.show()
plot_bifurcation_extended()
扩展太阳辐照范围后，拖动滑块观察雪球地球融化的临界点：
In [18]:
interact(plot_bifurcation_extended,
S=FloatSlider(min=1200, max=1800, step=2, value=1368, description='S [W/m²]'));
interactive(children=(FloatSlider(value=1368.0, description='S [W/m²]', max=1800.0, min=1200.0, step=2.0), Out…
在这个模型中，温度变化相当平滑，除非温度（从低温）升至 -10°C 以上或（从高温）降至 10°C 以下，此时冰-反照率正反馈发挥作用并导致突然气候转变。
虽然这只是一个简单的假设模型，但这种突然气候转变在古气候记录和更现实的气候模型模拟中经常出现。
这个模拟表明，不应认为气候的稳定性是理所当然的，将当前的气候推离其历史稳定范围可能会引发灾难性的突然气候转变。
4) 走向现实的气候建模¶
在这个简单模型中，前工业气候 $T=14$°C 如此温暖，以至于地球上任何地方都没有冰。实际上，唯一两个有效的稳定气候是无冰或全球冰覆盖。
那么地球的前工业气候——相对稳定了数千年——在两极是如何有大量冰盖的呢？
“水星球”——简单的三维海洋世界气候模型
“水星球”是一个假设完全被单一全球海洋覆盖的行星的三维全球气候模拟。虽然这与地球（27%被陆地覆盖）非常不同，但“水星球”表现出与地球许多相同的特征，比上面的零维气候模型更现实。
水星球模拟展示了第三个平衡态：大部分液态海洋但两极有冰盖，加上零维模型中发现的两个平衡态。
本讲小结¶
本讲从计算思维的角度，展示了"非线性如何带来质变"：
1. 抽象（Abstraction）
首先用一个极简的模型 $f(y,a) = \mathrm{sign}(y) + a - y$ 抽象出迟滞现象的核心机制——函数的不连续性导致多个根共存。在攻克真实的气候问题之前，先用最简单的数学对象理解本质，这是一种重要的计算思维策略：先简单，再复杂。
2. 渐进式建模（Incremental Modeling）
在上一讲线性气候模型的基础上，仅做了一个修改——将常数反照率 $\alpha$ 改为温度依赖的分段函数 $\alpha(T)$。模型从一个稳定解变为三个平衡态（两个稳定、一个不稳定）。这体现了建模中的关键洞察：线性是局部的近似，非线性效应可能主导全局行为。
3. 计算实验（Computational Experiment）
运用了多种计算手段探索非线性系统的丰富动力学：
- 参数扫描：遍历初始温度 $T_0$，发现双稳态与不稳定分界点
- 方向场可视化：用矢量场直观展示 ODE 解的走向
- 分岔图：通过太阳辐照 $S$ 的来回扫描，揭示迟滞环和突变点
- 交互式探索：滑块实时调参，建立对系统行为的直觉