第10章冰流动力学：Flowline模型

flowline.py 是OGGM中规模最大的单一源文件，约4853行，包含冰川动力学的主要实现：从冰床横截面形状（bed cross-section geometry）参数化到SIA方程数值求解，再到多种演化模型的构建。本章关注它如何把第9章反演的冰厚作为初始条件、把第8章的物质平衡作为气候强迫，进而驱动冰川随时间演化。

初学者阅读路径

第一遍只需读懂三个对象：Flowline保存一维几何，MassBalanceModel提供气候强迫，FlowlineModel根据SIA推进冰川几何。冰床横截面（bed cross-section）类、支流通量合并（tributary flux merging）和求解器细节可以在能跑通教程6之后再回看。

模型假设边界

OGGM主工作流使用一维流线表达复杂山地冰川，并以SIA为核心动力学近似。它适合区域到全球尺度的大量冰川模拟，但不显式求解纵向应力、完整三维流场或复杂基底水文。解释个别陡峭冰瀑、跃动冰川或快速入海型冰川时，应把这些假设作为结果讨论的一部分。

10.1 Flowline基础类

10.1.1 Flowline类（第49行）

Flowline 是OGGM中冰川几何表示的核心数据结构，继承自 Centerline 并添加了冰厚相关的属性。每条Flowline代表冰川的一条分支——主流线或支流线——携带该分支沿线的所有物理状态（厚度、表面高程、宽度、体积、冰通量等）。

class Flowline(Centerline):
    # 核心几何属性（从Centerline继承）
    # line: (n,)-形 numpy数组 — 沿流线的累积距离 (m)
    # nx: int — 网格点数
    # dx: float — 网格间距 (m)

    # 冰川状态属性
    thick = None          # (nx,)-形数组 — 冰厚 (m)
    surface_h = None      # (nx,)-形数组 — 表面高程 (m)
    widths = None         # (nx,)-形数组 — 表面宽度 (m)
    widths_m = None       # (nx,)-形数组 — 中间宽度（用于通量计算）
    section = None        # (nx,)-形数组 — 横截面积 (m2)
    bed_h = None          # (nx,)-形数组 — 冰床高程 (m)

    # 诊断属性
    volume_m3 = None      # 冰川总体积
    length_m = None       # 冰川总长度
    flux = None           # (nx,)-形数组 — 冰通量 (m3 s-1)
    is_rectangular = [...]  # bool 数组 — 每点的冰床横截面（bed cross-section）类型

surface_h = bed_h + thick，这是基本的几何关系。但surface_h不保证处处大于bed_h——当冰厚为零时两者相等。 widths_m 用于交错网格上的通量计算，定义在相邻格点之间（沿流线方向的半步点）。

10.1.2 冰床横截面形状（bed cross-section geometry）类层级

冰床横截面（bed cross-section）的形状决定了冰流通过横截面的方式，这是SIA理论中形状因子（shape factor）f_d的来源。OGGM定义了四种冰床横截面形状（bed cross-section geometry）类，构成一个功能递进的层级：

ParabolicBedFlowline（第301行）

抛物线形冰床横截面（bed cross-section）是最经典的冰川沟谷表示，对应于Nye (1965)和Cuffey & Paterson (2010)的经典处理：

class ParabolicBedFlowline(Flowline):
    # 宽度由抛物线公式给出：w = sqrt(4 * h / bed_shape)
    # bed_shape 是抛物线曲率参数
    # 横截面积：section = (2/3) * w * h
    # 形状因子（shape factor） f_d = 0.5

抛物线的几何含义：对于给定的冰厚h，宽度w由抛物线曲率参数bed_shape决定。w和h呈平方根关系，意味着冰厚的增厚会导致宽度的渐进式扩展。横截面积为(2/3)*w*h，形状因子（shape factor）固定为0.5——这是半圆冰床横截面（bed cross-section）的理论值，SIA形状因子（shape factor）恰好等于抛物线面积与等效矩形面积之比。

RectangularBedFlowline（第347行）

class RectangularBedFlowline(Flowline):
    # section = w * h  （矩形面积）
    # 形状因子（shape factor） f_d = 1.0
    # 无侧向约束——全部应力由冰床承担

矩形冰床横截面（bed cross-section）假设侧壁垂直，并在一维SIA框架中不显式表示侧向阻力。这种假设适用于横向约束较弱的宽阔冰体；若用于狭窄山谷冰川，可能高估冰流速度。

TrapezoidalBedFlowline（第394行）

class TrapezoidalBedFlowline(Flowline):
    # 宽度随冰厚线性增长：w = w0 + lambda * h
    # lambda 是梯形张开参数（无量纲）
    # 横截面积：section = (w + w0) / 2 * h
    # 形状因子（shape factor）随 w0/w 比值变化，介于0.5到1.0之间

梯形是更灵活的冰床横截面形状（bed cross-section geometry）表示，通过参数lambda（来自配置或反演中的trapezoid_lambdas）调整侧壁张开程度。w0表示冰床底部宽度。

MixedBedFlowline（第451行）——默认方案

class MixedBedFlowline(Flowline):
    # 在每个网格点独立选择冰床横截面（bed cross-section）类型
    # is_trapezoid 布尔掩码区分应用哪种公式
    # 抛物线点：section = 2/3 * w * h, fd = 0.5
    # 梯形点：    section = (w+w0)/2 * h, fd = f(lambda)
    # 矩形点：    section = w * h, fd = 1.0

这是OGGM常用的冰床横截面形状（bed cross-section geometry）类。关键理念是：同一条冰川的不同位置可能需要不同的横截面表示。例如，在冰川上游宽阔盆地中梯形截面可能更合适，而在狭窄U形谷中抛物线截面可能更合适。is_trapezoid掩码允许逐格点切换公式，实现几何表示的灵活性。

10.2 FlowlineModel基类（第571行）

10.2.1 构造函数

class FlowlineModel:
    def __init__(self, flowlines, mb_model=None, glen_a=None,
                 fs=0, mb_elev_feedback='annual'):
        # flowlines: Flowline对象列表（主流线 + 支流线）
        # mb_model: MassBalanceModel实例（气候强迫）
        # glen_a: Glen流动律参数 (Pa^-3 s^-1)，默认2.4e-24
        # fs: 基底滑动因子，默认0（无滑动）
        # mb_elev_feedback: 'annual','monthly','always','never'

模型的核心物理常数在构造时预计算：_fd = 2 * glen_a / (n + 2)（形变流参数，第653行）。这个预计算避免了在每个时间步中重复评估多项式。参数fs（基底滑动）在构造后通过_fs属性被类似处理。

10.2.2 支流系统（tributary system）：reset_flowlines()

reset_flowlines() 构建 _tributary_indices——一个记录每条支流线与主流线连接位置和方式的内部数据结构。支流贡献的冰通量通过高斯核分布在主流线上：

# 支流通量在主流线上的贡献
# 分布核大小自适应：5, 7, 或 9 个格点
# 更宽的冰川 → 更大的核 → 更平滑的通量合并
for trib in tributary_flowlines:
    junction_pos = trib.line[-1]          # 支流汇入主流线的位置
    kernel_width = adapt_to_glacier_size  # 自适应核宽度
    # 高斯权重分布到附近的主流线格点

高斯核的使用避免了支流通量集中到一个网格点导致的数值震荡。核大小的自适应逻辑考虑了冰川宽度——更宽的冰川下游使用9点核，更窄的山谷冰川使用5点核。

10.2.3 MB缓存策略

mb_elev_feedback 参数控制何时重新计算物质平衡：

模式	行为	适用场景
`'always'`	每次需要MB时重新计算	物理精确但极慢
`'monthly'`	每个月重新计算一次（12次/年）	亚季节分辨率的模拟
`'annual'`	每年重新计算一次	常用折中：兼顾精度与性能
`'never'`	使用初始MB，永不更新	测试或教学目的

缓存的核心设计逻辑是：在时间步之间，表面高程变化通常很小（典型值为每年几厘米到几米），因此每月或每年更新一次MB不会显著损失精度。通过跳过不必要的MB计算（尤其是基于气候数据的月度MB模型），性能可提升1-2个数量级。

10.2.4 时间步进：run_until()与run_until_and_store()

run_until(y1)（基本方法）将模型从当前时间向前积分到目标年份y1。它采用自适应时间步长策略：

外层：按年份推进。每年开始时检查是否需要更新MB（取决于mb_elev_feedback模式）。
内层：在每一年内，以月度步长推进（12步）。在每个月度步中，CFL条件可能触发进一步的子步分解。

run_until_and_store(y1)（第922行）是run_until的扩展版本，将演化过程记录到xarray.Dataset中：

def run_until_and_store(self, y1):
    # 创建xarray Dataset，包含时间坐标
    ds = xr.Dataset(
        coords={'time': [self.yr, y1]},
    )
    # 在每一步存储诊断变量
    # 变量包括：volume, area, length, ela, calving_rate
    # 每个变量是 shape=(time, flowline_idx, nx) 的DataArray
    return ds

存储的xarray Dataset使得后处理、可视化和分析变得非常方便——时间序列的体积变化、末端的进退、ELA的移动等都可以直接从Dataset中提取。

10.2.5 平衡态求解：run_until_equilibrium()（第1513行）

这个方法迭代运行模型直到冰川达到平衡态：

def run_until_equilibrium(self, rate=1e-3):
    while True:
        V_before = self.volume_m3
        self.run_until(self.yr + 1)
        V_after = self.volume_m3
        if abs(V_after - V_before) / V_before < rate:
            break   # |dV/V| < 0.1% 时视为平衡

收敛判据通常基于体积的年际变化率。平衡态求解常用于：(1) 理想化预热（spinup）或平衡实验；(2) 敏感性实验；(3) 理论分析，即理解给定气候下冰川的平衡响应。第13章讨论的动态预热（dynamic spinup）并不简单等同于求平衡态，而是以观测面积或体积为目标进行近期历史初始化。

10.3 FluxBasedModel（第1568行）

10.3.1 算法架构

FluxBasedModel 是OGGM的重要显式流线动力学模型。OGGM v1.6.3的配置文件默认evolution_model为SemiImplicit，但FluxBasedModel仍常用于多流线、崩解和需要显式通量处理的实验。它采用交错网格通量形式，每个时间步包含两个循环：

循环1（通量计算，第1771-1866行）：

for i in range(nx):
    # 步骤1: 支流汇合 — 在交汇点附加下游表面高程
    append_tributary_surface_h()

    # 步骤2: 崩解限制器 — 将水面以下冰厚截断
    if is_tidewater and bed_below_water:
        thick = clip_to_water_level()

    # 步骤3: 交错表面梯度
    slope_stag = (surface_h[i] - surface_h[i+1]) / dx

    # 步骤4: 交错厚度（相邻格点平均）
    thick_stag = (thick[i] + thick[i+1]) / 2

    # 步骤5: SIA流速（形变 + 滑动）
    tau = rho * g * slope_stag   # 基底剪应力
    u_deform = thick_stag^(n+1) * fd * tau^n
    u_sliding = thick_stag^(n-1) * fs * tau^n
    u_total = u_deform + u_sliding

    # 步骤6: 冰通量 = 流速 × 交错截面积
    flux_stag[i] = u_total * section_stag[i]

# 步骤7: CFL条件 — 确定子时间步长
dt = cfl * dx / max_u   # cfl 默认为 0.02

循环2（质量重新分布，第1874-1971行）：

for i in range(nx):
    # 连续性方程：d(section)/dt = -div(flux) + mb
    new_section = section[i] \
        + (flux_stag[i] - flux_stag[i+1]) * dt / dx \
        + trib_flux_contribution * dt / dx \
        + mb[i] * dt

    # 截断：截面积不能为负
    new_section = max(0, new_section)

    # 从截面积反算新厚度（取决于冰床横截面形状（bed cross-section geometry））
    new_thick = invert_section_to_thick(new_section, bed_shape)

    # 通过高斯核向下游传递通量
    distribute_flux_downstream()

中级

CFL数0.02的含义

FluxBasedModel使用非常保守的CFL数0.02。对于典型冰川流速（100 m/yr）和网格间距（200 m），这意味着dt ~ 0.04年，每年需要约25个时间步。与SemiImplicitModel的CFL 0.5相比，FluxBasedModel每个物理年需要更多的数值步，但每个步的计算量远少（无需矩阵求解）。

10.3.2 表面速度因子

SIA方程给出的是深度平均速度，但观测和诊断通常使用表面速度。两者之间的转换由经典冰川学关系给出：

# 深度平均速度 → 表面速度的转换因子
_surf_vel_fac = (n + 2) / (n + 1)   # n=3 时 = 5/4 = 1.25

这个因子来源于SIA速度剖面的解析积分（Cuffey & Paterson 2010, Eq. 8.34）。对于n=3，表面速度比深度平均速度高25%。在OGGM中，这个因子还可以乘以cfg.PARAMS['surface_velocity_factor']进行进一步的经验调整。

10.4 备选模型方案

10.4.1 SemiImplicitModel（第2152行）

半隐式方案通过将SIA方程转换为扩散问题来获得更好的稳定性。核心是一种类似Crank-Nicolson的格式：

# 扩散系数：D = (fd*h^(n+2) + fs*h^n) * rho*g * (w0+w)/2 * |dsdx|^(n-1)
# 对于纯形变流（fs=0）且抛物线冰床横截面（bed cross-section）：
# D = fd * h^(n+2) * rho*g * w * |dsdx|^(n-1)

# 三对角系统：d0[i]*S[i] + dm[i]*S[i-1] + dp[i]*S[i+1] = RHS
# 通过 scipy.linalg.solve_banded 求解

CFL条件从对流约束变为扩散约束：dt <= cfl * dx^2 / max(D/w)，cfl约为0.5。这比FluxBasedModel的CFL 0.02大幅宽松，在空间网格较粗时半隐式模型可以省去很多子步。但代价是：(1) 需要构建和求解三对角矩阵；(2) 仅支持单条流线（无支流系统（tributary system））。

进阶

SemiImplicitModel的床高程修正项

半隐式扩散方案的一个独特特征是包含了床高程修正项：b_corr = -D * d(bed)/dx。当冰床表面存在急剧的高程变化时（如反向坡度或冰瀑），这个项防止数值不稳定，确保质量守恒在变坡处保持精确。FluxBasedModel中这个效应的处理是隐式的（通过交错梯度），但SemiImplicitModel要求显式的修正项。

10.4.2 KarthausModel（第2067行）

KarthausModel是一个简化的显式方案，设计用于教学和概念验证。它直接从SIA计算流速和通量，使用简单的向前欧拉时间步进。默认参数经过调优使模型行为直观易懂，但物理精度不如FluxBasedModel。

10.4.3 MassRedistributionCurveModel（第2810行）

这是一个基于经验的方法，使用Huss et al. (2010)的delta-h参数化。核心思想是：

# 冰川响应由经验的"质量重新分布曲线"参数化
# delta_h(z) = f(z, glacier_size_class)
# 曲线来自Huss等人对34条瑞士冰川的观测研究
mass_redistribution_curve_huss(...)   # 第3068行

# 三种前进方法（advance_method 0, 1, 2）
# 方法0：固定曲线，截面积被MB比例缩放
# 方法1：MB驱动冰川末端前进/后退
# 方法2：MB同时驱动厚度变化和末端变动

该方法完全跳过了SIA冰流计算，用统计规律代替了物理模拟。它的计算成本极低（每冰川每年只需几分之一秒），适用于数千到数万条冰川的集合模拟，但牺牲了物理精确性——无法表示冰流的动力学非线性反馈。

10.4.4 FileModel（第2658行）

FileModel是一个"鸭子类型"读取器，允许将之前存储的模拟输出重新加载为一个可查询的FlowlineModel。它不执行任何物理计算，而是提供与FlowlineModel相同的查询接口（volume, area, length等），方便后处理脚本的复用。

10.5 模型初始化与运行函数

10.5.1 init_present_time_glacier（第3188行）

这是模型初始化的主入口。它负责：

加载反演的冰厚和流线。
调用 decide_evolution_model（第3342行）选择适当的动力学模型类。
处理崩解参数（如果适用）。
设置初始表面高程和体积。
可选：应用动态预热（dynamic spinup）的温度偏差校准。

10.5.2 decide_evolution_model（第3342行）

模型选择逻辑根据配置参数自动决定使用哪个模型类：

def decide_evolution_model(gdir):
    # 检查 cfg.PARAMS['evolution_model']
    # 可选值：'FluxBased', 'SemiImplicit', 'Karthaus',
    #         'MassRedistributionCurve', 'File'
# v1.6.3 params.cfg 默认：'SemiImplicit'

10.5.3 顶层运行函数

OGGM提供了封装完整模拟流程的顶层函数：

flowline_model_run（第3362行）：单条冰川的完整前向模拟，从初始化到指定结束年份。
run_from_climate_data（第3786行）：使用GCM/RCM的气候数据驱动模拟。处理气候模型输出的网格化降水和温度场。
run_random_climate（第3546行）：使用随机气候生成器进行不确定性量化或统计实验。
run_constant_climate（第3670行）：使用固定气候（通常为参考期气候态）进行平衡态实验。
run_with_hydro（第3939行）：耦合冰川动力学与水文模型，输出径流时间序列。

10.6 停止准则（第4371行）

模拟的终止由以下任一条件触发：

准则	类型	描述
`zero_glacier`	消融完成	冰川体积降为0（某些区域可能保留多年冰）
`spec_mb`	超过气候范围	MB模型在其标定范围外被查询（如年MB超出+/-5 m w.e.）
`equilibrium`	达到平衡	体积变化率`\|dV/V\| < 阈值`（用于平衡态求解）

停止准则的检查集成在run_until的每个年度步中。一旦触发，模型返回一个状态码，供调用代码确定是继续其他步骤还是终止模拟。

10.7 关键参数速查

参数名	默认值	作用	文件位置
`glen_a`	2.4e-24	Glen流动律参数	第653行
`fs`	0.0	基底滑动因子	构造函数参数
`n`	3	Glen流动律指数	硬编码
`cfl`	0.02 (Flux) / 0.5 (SemiImp)	CFL数	模型类默认
`mb_elev_feedback`	'annual'	MB高程反馈频率	构造函数参数
`cfl_min_dt`	60秒	最小子步时长	第1832行
`surface_velocity_factor`	1.0	表面速度经验乘数	cfg.PARAMS
`evolution_model`	'SemiImplicit'	配置文件默认动力学模型	cfg.PARAMS

10.8 总结

flowline.py 是OGGM动力学核心的完整实现：从冰床横截面（bed cross-section）几何的四种形状参数化，到基于SIA的两种数值方案（显式FluxBasedModel和半隐式SemiImplicitModel），再到经验的MassRedistributionCurveModel。支流系统（tributary system）通过高斯核将多条流线合并为统一的冰流通量场。自适应时间步进和MB缓存策略在物理精度与计算效率之间取得了工程平衡。

FluxBasedModel凭借显式、直观的通量计算流程和支持多流线的能力，在许多OGGM应用中仍十分重要；SemiImplicitModel则在默认配置中承担一维动力学正演的角色。二者共同构成第9章冰厚反演与后续工作流执行之间的动力学连接环节。