第11章 数值方法详解

第10章从类层级和算法架构的角度描述了OGGM的动力学模型。本章深入到数值细节——SIA方程的离散化、两种时间步进方案的稳定性分析、交错网格的设计原理、以及冰速计算和通量合并的数学基础。理解这些数值方法对于调试模拟异常、解释模型行为以及实现自定义动力学方案至关重要。

第一遍可以怎么读

如果你的目标是先完成OGGM应用，可以先读11.1、11.2和11.5，理解SIA、求解器选择和时间步长。矩阵求解、交错网格和扩散率推导是解释数值异常或开发模型时才需要深入掌握的内容。

11.1 浅冰近似（SIA）在OGGM中的应用

11.1.1 完整的SIA方程

浅冰近似是冰川动力学的零阶近似，假设冰厚远小于水平范围，且侧向应力梯度可忽略。完整的SIA深度平均速度表达式为： 

# 深度平均冰流速度（SIA）
u_bar = -2 * A * (rho * g)^n / (n + 2)
      * |grad(S)|^(n - 1)
      * grad(S)
      * H^(n + 1)

其中:
  A   = glen_a  — Glen流动律参数 (Pa^-n s^-1)
  rho = 900    — 冰密度 (kg m^-3)
  g   = 9.81   — 重力加速度 (m s^-2)
  n   = 3      — Glen指数
  S   = bed_h + H  — 冰表面高程 (m)
  H   = 冰厚 (m)

Glen流动律是这个方程组的物理基础： 

# Glen流动律（Norton-Hoff定律）
epsilon_dot = A * tau^n

其中:
  epsilon_dot — 有效应变率 (s^-1)
  tau         — 有效剪应力 (Pa)
  n           — 应力指数，标准冰川学值 = 3

n = 3 意味着冰的黏度随应力的三次方变化——三倍的应力导致27倍的应变率。这种强烈的非线性是冰川动力学的核心特征：流速对厚度和表面坡度高度敏感。

11.1.2 SIA的有效性范围

SIA的数学推导基于小纵横比假设：冰厚H远小于水平尺度L，即H/L（纵横比）小到可以忽略纵向应力梯度。这在下列条件下成立：

在冰架（漂浮冰）、冰瀑（陡峭地形）或接地线附近，SIA的近似失效，需要更复杂的"高阶"模型。但OGGM的设计目标就是以SIA为物理基础，用经验修正弥补其不足。

11.2 显式方案与半隐式方案的对比

11.2.1 FluxBasedModel：向前欧拉

FluxBasedModel使用显式的向前欧拉时间离散： 

# 显式向前欧拉：H^(k+1) = H^(k) + dt * [MB - div(q)]
# 其中 div(q) = (q[i+1/2] - q[i-1/2]) / dx
# q 在交错格点上计算（使用当前时刻的H和S）

显式方案的两大限制：

  1. 时间步长限制 —— CFL条件 u * dt / dx < cfl_max。对于流速100 m/yr、网格200 m、CFL=0.02的情况，dt约0.05年。高流速或细网格需要更多子步。
  2. 无条件扩散不稳定 —— 当扩散项占主导时（厚冰 + 陡坡），即使CFL数满足名义值，向前欧拉法仍可能表现出高频空间振荡。这是显式扩散格式的固有问题。

稳定性分析（von Neumann）表明，对于纯对流（fs=0且低n），向前欧拉的放大因子为|G| = |1 - i*u*dt/dx*sin(k*dx)|，稳定条件为|G| <= 1，直接导出CFL约束。

11.2.2 SemiImplicitModel：类Crank-Nicolson格式

SemiImplicitModel将SIA重新表述为扩散问题： 

# SIA扩散形式：dH/dt = MB + div(D * grad(S))
# 其中 D = gamma * H^(n+2) * |grad(S)|^(n-1) 是扩散系数
# gamma = 2*A*(rho*g)^n / (n+2)

# 半隐式离散：
# H^(k+1) = H^(k) + dt*MB + dt*div(D^(k)*grad(S^(k+1)))
# 在新表面S^(k+1)的梯度上使用隐式处理

CFL条件从对流（u*dt/dx）变为扩散（D*dt/dx^2），容许的时间步长大幅增加： 

# FluxBased： dt <= 0.02 * dx / u     —— 对 dx 线性缩放
# SemiImplicit：dt <= 0.50 * dx^2 / D  —— 对 dx 平方缩放

# 当 dx 较大时（>200 m），半隐式的优势十分明显
中级
何时选择哪个方案

选择FluxBased：有支流系统（tributary system）时需要合并多条流线；需要显式通量诊断或崩解相关处理；低分辨率但大量的冰川集合模拟中计算成本可接受。注意OGGM v1.6.3配置文件默认evolution_model为SemiImplicit。
选择SemiImplicit：单条流线的精细研究；需要高空间分辨率（dx < 100 m）以避免过多的子步；需要严格数值稳定性（如含有急剧坡度变化的冰床横截面（bed cross-section））。
常见权衡：每步计算成本——FluxBased约O(nx)，SemiImplicit约O(nx)（三对角求解为线性），但FluxBased需要更多的步数。

11.3 交错网格（staggered grid）离散化

11.3.1 设计原理

OGGM的通量计算采用Arakawa C网格（交错网格）——厚度和表面高程定义在网格中心，冰通量定义在网格面（相邻格点之间）： 

# 网格布局示意：
#   cell:     [ i-1 ]    [  i  ]    [ i+1 ]
#   center:    H[i-1]     H[i]       H[i+1]
#   face:  q[i-1/2]   q[i+1/2]  q[i+3/2]

# 交错量计算：
slope_stag[i] = (surface_h[i] - surface_h[i+1]) / dx
thick_stag[i] = (thick[i] + thick[i+1]) / 2

# 通量在面上：
flux_stag[i] = u(thick_stag[i], slope_stag[i]) * section_stag[i]

# 厚度更新在中心：
new_thick[i] = thick[i] + dt * (flux_stag[i-1] - flux_stag[i]) / dx

11.3.2 为什么交错网格避免了棋盘格模式

非交错网格（所有量在网格中心）的一个已知问题是棋盘格模式——厚度场的2-dx波长振荡不会产生通量，因此不会被时间步进衰减。考虑中心梯度： 

# 非交错梯度（中心差分）：dsdx[i] = (S[i+1] - S[i-1]) / (2*dx)
# 振荡模式 S = [..., +a, -a, +a, -a, ...]
#  => dsdx = 0 （梯度为零，尽管表面有剧烈波动）
#  => 通量 = 0 （没有动力驱散这些波动）

# 交错梯度：dsdx[i] = (S[i] - S[i+1]) / dx
# 振荡模式 => dsdx = +/- 2a/dx （非零梯度）
#  => 通量将驱散这些波动

交错网格（staggered grid）在相邻格点间而非越过格点计算梯度，确保始终存在物理耦合，从而自然抑制数值网格尺度的振荡。

11.4 带状矩阵求解器（SemiImplicitModel）

11.4.1 三对角系统

半隐式离散化得到一个三对角线性系统： 

# 三对角系统：对于 i = 0, 1, ..., nx-1
d0[i] * S_new[i] + dm[i] * S_new[i-1] + dp[i] * S_new[i+1] = RHS[i]

其中:
  d0[i] = 1 + dt * (D_front[i] + D_back[i]) / dx^2
  dm[i] = -dt * D_back[i] / dx^2
  dp[i] = -dt * D_front[i] / dx^2
  RHS[i] = H_old[i] + dt * MB[i] + dt * b_corr[i]

# 使用scipy的带状矩阵求解器：
from scipy.linalg import solve_banded
S_new = solve_banded((1, 1), ab_matrix, RHS)

solve_banded(l_and_u, ab, b) 的(1,1)表示主对角线上下各有一条子对角线（三对角），ab矩阵以压缩格式存储非零对角线。nx通常为50-500，矩阵求解的计算成本O(nx)几乎可以忽略。

11.4.2 消融区边界条件

在冰川末端（消融区最后一个有冰的格点），施加零厚度边界条件：H[nx-1] = 0。对应地，三对角系统的最后一行被修正为简单的S_new[nx-1] = bed_h[nx-1]（表面等于冰床，无冰）。这确保了消融区边界处质量守恒的正确处理。

11.5 自适应时间步长

11.5.1 三层嵌套时间循环

OGGM动力学模拟的时间步进由三层嵌套控制： 

# 第1层：年度步（yearly step）
while yr < y1:
    mb_updated = False
    for month in range(12):           # 第2层：月度步（monthly step）
        if feedback_mode == 'monthly':
            update_mb_model()
        remaining = 1/12               # 目标时间 = 一个月（~0.083年）
        while remaining > 0:            # 第3层：CFL子步（sub-step）
            dt_cfl = cfl * dx / max_u
            dt = min(dt_cfl, remaining)
            step(dt)
            remaining -= dt

最小子步时长由 cfl_min_dt 限制（默认60秒）。当冰川流速极高（如 > 10000 m/yr，在巨型入海型冰川（tidewater glaciers）中可达此值）且网格很细时，这个下限防止了无限细分的退化情况。

11.5.2 计算效率的优化

在典型生产模拟中（dx ~ 100-200 m, u ~ 50-300 m/yr），CFL子步为每年10-50步。月度循环的开销可忽略（每年仅12次MB更新和循环管理）。年度步主要用于触发MB缓存更新和诊断检查，频率为每年一次。

11.6 质量守恒

11.6.1 MassConservationChecker包装器（第2031行）

MassConservationChecker 是对任何FlowlineModel的一个包装器，在每个时间步前后记录总质量并报告任何不守恒情况： 

class MassConservationChecker:
    def step(self, dt):
        mass_before = self.model.volume_m3
        self.model.step(dt)
        mass_after = self.model.volume_m3
        # 预期质量变化 = MB积分 + 通量边界条件
        # 实际质量变化 = mass_after - mass_before
        # 差异 = 当前步的数值质量泄漏

这对于开发和调试新动力学方案至关重要——任何非零的质量差异都表示数值方案中的实现错误。

11.6.2 通量闸门（flux gate）边界条件

在冰分（glacier divide，冰川最高点，u=0），冰通量为零——没有上游物质来源。这个物理约束在数值上通过通量闸门（flux gate）边界条件实现：q[0] = 0（流线起点的通量强制为零）。这防止了"凭空生成"的冰从上游进入模型域。

11.6.3 体积追踪与崩解桶（calving bucket）机制

质量守恒需要通过calving_bucket（崩解桶）模式处理崩解——崩解以离散的网格单元大小移除冰，但每次只有部分网格单元需要移除。桶机制（bucket mechanism）累积分数质量直到达到完整网格单元的体积，然后一次性移除——这种方法保持了体积追踪的精确性而不引入分数网格单元。

11.7 高程反馈机制

mb_elev_feedback 是OGGM在高程-物质平衡耦合上的核心控制参数。当冰表面因动力学而改变高度时，物质平衡必须更新以反映新的表面高程： 

# 反馈周期对数值行为的影响：
'always':  每个子步更新MB — 物理上最精确，但严重增加计算成本
'monthly': 每12个子步更新MB — 接近连续反馈
'annual':  每12*N_substep步更新MB — 默认
'never':   从不更新MB — 固定气候实验（物理上与"平衡态"不同）

高程反馈中的关键权衡：更频繁的反馈增加精确度但也增加了计算成本（尤其在MB模型基于月度气候数据时，每次更新可能触发对多变量格网数据（gridded data）的查询）。对于多数冰川，每年1-10米的高程变化意味着年度反馈的误差在约0.1-0.5 m w.e. MB偏差以内，这在气候强迫的不确定性中是可以接受的。

进阶
高程反馈与模型稳定性

高程反馈本身是一种正反馈：表面下降 → 更多消融（气温更高） → 更多下降。在显式方案中，这种正反馈可能与CFL不稳定耦合，触发"runaway thinning"。SemiImplicitModel的隐式处理能部分缓解这一问题，但合理的mb_elev_feedback频率（如'annual'而非'monthly'）仍然是一个重要的稳定性约束。对于入海型冰川（tidewater glacier）的快速消融场景，可能需要更保守的CFL设置。

11.8 支流通量合并（tributary flux merging）

11.8.1 高斯加权的动机

支流贡献到主流线的冰通量通过高斯加权核分布在多个格点上。核心问题是：支流在一个精确的几何点（交汇处）汇入主流线，但冰通量不应集中在一个网格点上——那样会造成急剧的通量梯度（数值不稳定）。高斯分布将通量"涂抹"在交汇点附近的几个格点上： 

# 支流通量分布权重
w[j] ~ exp(-(x[j] - x_junction)^2 / (2 * sigma^2))
# sigma 取决于冰川宽度：wider -> larger sigma -> more smoothing
# 核大小：5, 7, 或 9 个格点

核大小的自适应选择确保了较小的山谷冰川（支流交汇处下游宽度较窄）的支流贡献不会被过度涂抹，而较大的入海型冰川（tidewater glaciers）/冰盖的宽阔接收区能足够平滑地吸收支流输入。

11.9 冰速计算

11.9.1 深度平均速度与表面速度

SIA给出的速度是深度平均的，但表面速度（卫星InSAR和特征跟踪的直接观测量）是表面速度。两者之间的关系为： 

# 深度平均速度 → 表面速度
u_surface = u_bar * (n + 2) / (n + 1)

# 对于 n=3: u_surface = u_bar * 5/4 = 1.25 * u_bar
# 物理原因：速度剖面是垂向的(n+1)/(n+2)指数分布

# 速度在交错面上的平均值（用于诊断输出）：
u_at_center[i] = (u_stag[i-1] + u_stag[i]) / 2

11.9.2 速度诊断与观测校准

在诊断输出中，OGGM将交错速度平均到网格中心，然后乘以surface_velocity_factor（来自cfg.PARAMS）进行经验调整。这个因子的默认值为1.0，但在某些标定到卫星观测的项目中被调至0.5-2.0不等——反映了SIA速度理论值与观测值之间的系统性偏差。

11.10 SIA扩散率推导

将SIA重新表述为扩散形式是半隐式方案的基础： 

# 通量形式：q = -gamma * H^(n+2) * |dS/dx|^(n-1) * dS/dx
# 其中 gamma = 2*A*(rho*g)^n / (n+2)

# 连续性方程：dH/dt = MB - dq/dx
#            = MB + d/dx[D * dS/dx]
# 其中扩散系数 D = gamma * H^(n+2) * |dS/dx|^(n-1)

# 对 n=3: D = gamma * H^5 * (dS/dx)^2
# 扩散率对厚度的强敏感性（H^5！）是冰川响应非线性的根源

扩散系数D中H5的强依赖意味着当冰厚加倍（2x）时扩散率增加32倍，这解释了为什么厚冰川对气候变化更敏感（响应更快），而薄冰川的变化更缓慢。

11.11 通量散度作为残差

dh/dt - climatic_mb 即通量散度残差是一个核心诊断量。在平衡态中，残差应为零——MB在所有位置的积累都被冰流动力学精确平衡。在非平衡态中，残差表示局部正在发生质量净增（整理，div为正）或净损失（消融，div为负）。这个残差对于理解冰川不同部位的质量变化来源非常重要。 

# 诊断残差（用于模型验证）：
residual = dh/dt - climatic_mb
# 如果 residual = 0: 冰川处于平衡（冰流完美匹配气候）
# 如果 residual > 0: 动力学敛聚（局部冰流快于气候可提供的）
# 如果 residual < 0: 动力学扩散（局部冰流不足以带走所有积累）

11.12 关键参数速查

参数默认值含义
n3Glen流动律应力指数
glen_a2.4e-24Glen流动律参数 (Pa-n s-1)
cfl (FluxBased)0.02CFL数（对流稳定性）
cfl (SemiImplicit)0.50CFL数（扩散稳定性）
cfl_min_dt60 秒最小允许子步时长
(n+2)/(n+1)1.25表面速度/深度平均速度比率
D = gamma*H^5*(S')^2-SIA扩散系数（n=3时的形式）

11.13 总结

本章深入探讨了OGGM动力学代码背后的数值方法：SIA方程的两类离散化方案（显式通量形式和半隐式扩散形式）、CFL条件导出的稳定约束、交错网格的棋盘格模式抑制原理、以及支流通量的高斯合并技术。理解这些数值基础对于有效使用OGGM、扩展其模型系统、以及诊断模拟中的非预期行为都至关重要。

半隐式方案通过将SIA线性化并求解三对角系统获得巨大的稳定性优势，但其单流线限制限制了其在多条支流冰川上的应用。显式方案更灵活但步长更小——这是冰川模拟中经典的"精度-效率"张力的具体体现。