第14章二维SIA模块

sia2d.py 约449行，实现了OGGM的二维SIA网格模型——与第10章讨论的一维流线模型不同，这个模块在规则网格上直接求解二维SIA方程。它适用于理想化实验、小冰川的详细模拟、以及1D流线与2D网格模型之间的方法对比。代码量虽小，但包含了从上游加权差分到CFL时间步进的完整二维冰流求解器。

先明确定位

OGGM的常规区域和全球模拟主要依赖一维流线工作流。二维SIA模块更适合理想化实验、教学和数值方法比较；除非已经验证输入、边界条件和适用性，不建议把它作为常规区域研究的默认工具。

14.1 概述与设计理念

1D流线模型（FlowlineModel及其子类）通过将冰川简化为沿中心线的带状几何来降低计算维度。这种方法在计算效率和物理精度之间取得了平衡，但以牺牲二维细节为代价——例如横谷宽度变化、侧向冰流效应、和不规则流域几何。2D SIA模块填补了这一空白，提供了在完整的二维地形网格上直接求解SIA方程的能力。

关键区别：

特性	1D流线模型（第10章）	2D SIA模块（本章）
几何表示	沿中心线的一维带状 + 冰床横截面形状（bed cross-section geometry）参数	二维规则网格（dx x dy）
侧向约束	通过宽度或梯形参数隐式表示	通过实际地形显式表示
支流系统（tributary system）	支持——支流汇入主流线	不支持——所有通量在网格中直接计算
崩解参数化	全功能k-崩解定律	未实现（简化设计）
计算成本	O(nx) per step	O(nx*ny) per step
适用场景	生产批量模拟	理想化实验、小冰川、方法研究

14.2 冰边界滤波：filter_ice_border（第10行）

在2D网格中，冰川边缘的处理是一个微妙的问题。冰覆盖区域与无冰区域之间的边界处，SIA梯度可能产生非物理的冰流。OGGM通过一个简单的域边界滤波来处理：

def filter_ice_border(ice_thick):
    # 将域边缘的一个像素宽的边界上的冰厚置为零
    ice_thick[0, :] = 0
    ice_thick[-1, :] = 0
    ice_thick[:, 0] = 0
    ice_thick[:, -1] = 0
    return ice_thick

这个边界条件在物理上等价于"冰川不能扩展到计算域之外"——在计算域边缘，冰被强制为零，起到了吸收边界的作用。默认的ice_thick_filter就是filter_ice_border，但用户可以替换为自定义的滤波器（例如更复杂的缓冲带或渐变过滤）。

中级

为什么不使用更复杂的边界条件？

完整的SIA模拟需要处理流动边界（flux boundary）——本质上是一个具有移动边界的自由边界问题。但在1D流线模型中通过跟踪每个格点的冰/无冰状态来处理（通过厚度截断），2D模块采用更简单的域级零厚度边界。对于从小冰期以来一直退缩的冰川，计算域可以简单地设置为比当前冰川范围大一些，确保在整个模拟期间冰川始终在域内。

14.3 Model2D基类（第19行）

14.3.1 构造函数

class Model2D:
    def __init__(self, bed_topo, init_ice_thick, dx=None, dy=None,
                 mb_model=None, glen_a=None, ice_thick_filter=None):

        # 2D 地形网格
        self.bed_topo = bed_topo          # shape (ny, nx) — 冰床高程 (m)
        self.init_ice_thick = init_ice_thick # shape (ny, nx) — 初始冰厚 (m)
        self.dx = dx or 100             # 格点间距 x 方向 (m)
        self.dy = dy or dx              # 格点间距 y 方向 (m)，默认与dx相同

        # MB模型、Glen参数、高程反馈（同FlowlineModel）
        self.mb_model = mb_model
        self.glen_a = glen_a or cfg.PARAMS['glen_a']
        self.mb_elev_feedback = 'annual'

        # 冰厚滤波器（默认：filter_ice_border）
        self.ice_thick_filter = ice_thick_filter or filter_ice_border

14.3.2 核心属性

属性	类型	描述
`yr`	float	当前模型年份
`area_m2`	float	冰川面积 (m²) = sum(ice_thick > 0) * dx * dy
`volume_m3`	float	冰川体积 (m³) = sum(ice_thick) * dx * dy
`surface_h`	2D array	冰表面高程 = bed_topo + ice_thick（仅在有冰处）
`area_km2`	float	冰川面积 (km²)
`volume_km3`	float	冰川体积 (km³)

14.3.3 MB缓存与时间步进

Model2D 共享与FlowlineModel相同的MB缓存策略和run_until/run_until_equilibrium/run_until_and_store接口。这种接口一致性使得在1D和2D模型之间切换时上层调用代码无需改动。

14.4 Upstream2D模型（第247行）

14.4.1 上游加权差分格式

Upstream2D 是OGGM中实际的2D SIA求解器。其名称源于使用的上游加权差分——在计算扩散系数时，对冰厚的插值偏向于上游方向（较高表面的方向）：

class Upstream2D(Model2D):
    # gamma = 2 * A * (rho * g)^n / (n + 2)    （第299行）
    # 这是预计算的SIA常数

    # CFL = 0.124 （比Hindmarsh理论略大）
    # 理论最大值：1 / (2*(n+1)) = 1/8 = 0.125
    # OGGM使用 0.124 —— 略低于理论极限以确保稳定性

    # max_dt = 31 天  （避免过于粗大的月度时间步）

为什么CFL取0.124？对于二维非线性扩散方程，理论分析（Hindmarsh 2001）给出了稳定性极限1/(2*(n+1))。对于n=3，这是1/8 = 0.125。OGGM选择了0.124（略低于0.125），为数值舍入误差和边界效应留出了安全余地。

14.4.2 numpy ix_ 索引技巧

为了高效地执行交错网格操作，Upstream2D使用numpy.ix_来避免显式的Python循环：

# 使用 numpy.ix_ 创建开放网格索引
# 这允许向量化切片操作，避免效率低下的 Python for 循环
iy = np.arange(0, ny)
ix = np.arange(0, nx)
iy_dn, ix_dn = np.ix_(iy[1:-1], ix[1:-1])

# 这使得可以对内部域一次性进行向量化差分运算
# 而非逐行逐列循环

14.5 diffusion_upstream_2d：核心计算（第332行）

这是sia2d.py中最重要的函数——实现了二维SIA扩散-对流方程的完整空间离散：

def diffusion_upstream_2d(H, bed, dx, dy, gamma):
    # ============================================================
    # 步骤1：上游加权冰厚
    # ============================================================
    # H_centered = 0.5 * (H[i+1] + H[i])    （中心平均）
    # H_upstream = H[i+1] if S[i+1] > S[i] else H[i]  （上游选择）
    H_up = upstream_weighted(H, bed, direction)

    # ============================================================
    # 步骤2：表面梯度范数
    # ============================================================
    dsdx = (S[i+1, j] - S[i, j]) / dx       # x方向梯度（交错面）
    dsdy_avg = (S[i, j+1] - S[i, j-1]) / (2*dy)  # y方向均值（非交错）

    # 梯度范数（用于计算扩散率中的非线性因子）
    grad_S_norm = sqrt(dsdx^2 + dsdy_avg^2)^((n-1)/2)

    # ============================================================
    # 步骤3：扩散系数
    # ============================================================
    # D = gamma * H^(n+1) * H_upstream * grad_S_norm
    # 注意这里使用了上游加权的 H_upstream（而非 H^(n+2)）
    # 当 n=3: D = gamma * H^4 * H_upstream * (grad_S)^(2)

    # ============================================================
    # 步骤4：通量散度
    # ============================================================
    # div_x = (D_up*ds_up - D_dn*ds_dn) / dx^2
    # div_y = (D_left*ds_left - D_right*ds_right) / dy^2
    div_q = div_x + div_y

    # ============================================================
    # 步骤5：CFL时间步长
    # ============================================================
    dt = cfl * min(dx^2, dy^2) / max(D)

    return div_q, dt

进阶

上游加权 vs 中心加权的物理含义

上游加权（在较高表面方向上取冰厚）的物理基础是：冰流方向是向下梯度的，来自上游（较高表面）的冰通量携带了上游的厚度信号。在陡峭地形或厚度剧烈变化的位置，中心平均会引入来自"错误方向"的虚假厚度信息——例如下游无冰区的零厚度被平均进来，低估了实际流过来的冰量。上游加权避免了这个问题，保证了扩散系数的物理一致性。

14.6 Step方法（第435行）

Upstream2D的步进方法是一个简单的向前欧拉更新：

def step(self, dt):
    # 1. 获取当前MB（取决于高程反馈模式）
    mb = self.get_mb(self.surface_h, self.yr)

    # 2. 计算通量散度和CFL时间步长
    div_q, dt_cfl = diffusion_upstream_2d(
        self.ice_thick, self.bed_topo,
        self.dx, self.dy, self.gamma)

    # 3. 使用实际时间步长（用户指定的 dt 或 CFL 限制的 dt_cfl，取较小者）
    dt_actual = min(dt, dt_cfl)

    # 4. 向前欧拉更新：dH/dt = MB + div_q
    new_surface_h = self.surface_h + (mb + div_q) * dt_actual

    # 5. 从表面高程和冰床高程反算冰厚
    new_ice_thick = new_surface_h - self.bed_topo

    # 6. 截断：冰厚不能为负
    new_ice_thick = np.clip(new_ice_thick, 0, None)

    # 7. 应用冰厚滤波器（默认：零边界条件）
    self.ice_thick = self.ice_thick_filter(new_ice_thick)

与1D FluxBasedModel的两步法（先算通量、再更新厚度）相比，2D的步进法将所有操作合并在一个函数中——通量散度和厚度更新同时完成。这是因为2D的网格结构使得交错网格的存访模式更简单：在每个时间步中，所有2D数组一次性计算，无需在交错面和中心之间维护多个网格表示。

14.7 与1D流线模型的差异总结

14.7.1 几何表示

维度	1D流线	2D SIA
侧向约束	冰床横截面形状（bed cross-section geometry）参数（f_d, trapezoid_lambdas）隐式	实际二维冰床横截面（bed cross-section）地形显式
冰厚分布	沿中心线一维变化	在二维平面内自由变化
宽度变化	预定义的宽度函数 w(x)	由冰厚和地形共同决定（自然涌现）
几何假设	冰床横截面形状（bed cross-section geometry）在各位置独立假设	使用完整的DEM冰床横截面（bed cross-section）

14.7.2 数值方案

特性	1D FluxBased	1D SemiImplicit	2D Upstream2D
离散化	交错通量（面上）	中心扩散（点内）	交错扩散（面上）
CFL限制	u*dt/dx < 0.02	D*dt/dx^2 < 0.5	D*dt/dx^2 < 0.124
矩阵求解	否（显式）	是（三对角）	否（显式）
支流支持	是	否	不适用（2D网格无支流概念）
崩解支持	是（k-定律）	是（有限制）	否

14.7.3 适用性考量

何时使用2D模型：

理想化实验——如在规则沟谷或人工地形上测试SIA的物理行为。
小冰川（< 5 km²）——计算成本可控，且1D的冰床横截面形状（bed cross-section geometry）假设对小冰川不太适用。
方法研究——比较1D和2D表示对同一冰川的影响。
教学用途——2D SIA方程的实现比1D更直观（直接对应数学公式）。

何时使用1D模型：

生产批量模拟——数万条冰川的模拟在计算上仅在1D表示下可行。
有支流系统（tributary system）的冰川——2D模块不支持支流合并。
需要崩解功能——入海型冰川需要崩解参数化。
大规模空间分析——1D模型的输出更易于汇总和统计。

14.8 gis.py集成

2D SIA模块与gis.py中的GIS工具配合使用。典型的使用流程为：

# 1. 从GIS数据准备二维网格
from oggm import gis
gdir = gis.define_glacier_region(rgi_id, ...)

# 2. 提取DEM和冰川掩码
dem = gdir.get_dem()              # 二维冰床高程网格
mask = gdir.get_glacier_mask()    # 冰川覆盖标记

# 3. 初始化冰厚（可选：使用简单常数或反演结果插值）
init_ice_thick = np.where(mask, 50, 0)  # 50 m常数厚度

# 4. 创建2D模型并运行
from oggm.core import sia2d
model = sia2d.Upstream2D(dem, init_ice_thick, dx=50,
                          mb_model=mb_model)
model.run_until_equilibrium()

14.9 关键参数速查

参数	默认值	位置	含义
`glen_a`	2.4e-24	cfg.PARAMS	Glen流动律参数
`cfl`	0.124	Upstream2D硬编码	CFL数（略低于理论极限0.125）
`max_dt`	31天	Upstream2D硬编码	最大单步时间步长
`gamma`	2A(rho*g)^n/(n+2)	第299行计算	SIA预计算常数
`mb_elev_feedback`	'annual'	Model2D默认	MB高程反馈频率

14.10 总结

sia2d.py 在紧凑的449行中实现了一个完整的二维SIA冰流求解器，为OGGM的1D流线表示提供了补充方案。上游加权差分、numpy ix_向量化、和略低于理论极限的CFL数共同确保了数值的稳定性和效率。虽然没有支流系统（tributary system）和崩解参数化，但其简化的结构和与Model2D/FlowlineModel一致的接口使得方法对比和教学实验可以无缝进行。

对于大多数批量生产模拟，1D流线模型（第10章）仍然是最佳选择，但2D模块为研究者提供了一个宝贵的"基准真相"——可以在少量冰川上运行以验证1D模型的适用性和假设限制。