观潮 · 模式预报与气象 AI 观察

Suêtes:把非静力模式动力核写成可微的(JAX 实现)

2026-10-04

> Suêtes 是一个用 JAX 从头写的非静力有限区域大气动力核。它的特别之处不在于算得更快,而在于整条模式链路都能反向传播——一次反向计算就能得到预报量对地形高度、初始场、边界强迫、物理参数和数值离散本身的敏感度。

缺的是哪一块

现在能拿到梯度的气象模式,都是全球尺度的。NeuralGCM 那一类(可微谱动力核 + 神经网络次网格参数化,在线训练)已经证明梯度有用,但全球尺度跑不动对流允许分辨率。

区域模式这边正好反过来:WRF、MC2/GEM、AROME-France、COSMO、ICON-CLM 都能跑到 4 km 以下,但没有一个是可微的。作者原话:目前不存在可微的有限区域动力核,上面这些例子全是全球尺度的。

区域尺度做不了在线优化,根源在两处:

缺口 · 后果

没有区域可微动力核 · 学习到的组件无法在模式积分过程上做梯度优化

现有可微核都在 O(100) km 量级 · 训练 rollout 的成本随反向传播携带的自由度增长,全球尺度下 1 km 分辨率算不动

静力近似是第二道墙。现有可微大气核大多建立在静力近似上,而静力假设把垂直动量方程消掉了——对流爆发时的垂直加速度算不出来。要在 1 km 量级做对流允许模拟,这道墙绕不过去。

作者点明这是三个区别点中的第三个:梯度要穿过此前没在这个场景下被微分过的时间积分器——一整套两时间层半隐式半拉格朗日方案(继承自加拿大 MC2 模式,耦合声重波系统用 GMRES 预条件反演),加上带声学子步进的 Split-Explicit Runge-Kutta(沿袭 WRF)。

名字的来历

Suêtes 指的是加拿大布雷顿角高地西缘的下坡风暴(southeasterly downslope windstorms)。这种风背后的驱动山脊只有几十公里宽,在全球模式被平滑过的地形里直接消失了。

论文的 3.5 节拿 2025 年 2 月 14 日的一次真实 Wreckhouse 风暴做算例,域中心就放在 Wreckhouse,地形来自 GEBCO 高程。这个名字顺带说明了作者选区域模式而不是全球模式的动机:不是分辨率洁癖,是有些强迫本来就装不进全球模式的网格。

方程与离散

积分的是全可压缩 Euler 方程组,不是静力近似。离散用交错 Arakawa C 网格(Arakawa et al., 1977),垂直方向是一般地形跟随坐标,侧边界用 Davies 松弛(1976)向驱动数据过渡,标量通量形式输送(Lin and Rood, 1996),物理过程打包成一个紧凑的可微物理包。

热力学常数全文统一:p₀ = 10⁵ Pa,R_d = 287 J kg⁻¹ K⁻¹,c_p = 1004 J kg⁻¹ K⁻¹,c_vd = 717 J kg⁻¹ K⁻¹,密度在每个算例里都由方程 (6) 诊断得到。

作者明确交代可微的环节有哪些:地形跟随几何、空间离散、时间积分、初始条件、侧边界处理、内含的物理参数化。侧边界那一段是关键——很多区域模式在侧边界插值上做文章,那里断掉可微性,整条链就白搭。

模式几何与计算网格
模式几何与计算网格

论文图 1:Suêtes 的空间离散。有限区域用斜轴 stereographic 投影定义计算域,物理域的度量由地形跟随坐标变换给出。水平方向的计算坐标与物理域之间的这一层变换是后续所有梯度能回传的前提。

两套动力核,运行时切换

动力核 · 构成 · 优势区间 · 可微方式

**半隐式半拉格朗日(SISL)** · 轨迹法平流 + GMRES 求解耦合压力–速度系统 · 约 10 km 及更粗——这个尺度上显式时间步受稳定性而非精度限制,长时间步才有意义 · 用转置线性求解微分隐式更新,不存储正向 Krylov 迭代序列

**Split-Explicit Runge-Kutta** · 严格欧拉、声学子步进 · 对流允许(公里级)——精度本来就要求较短平流时间步,避开半拉格朗日插值有助于保留短波长结构 · 用检查点限制中间 RK 步与声学子步状态存储,反向时重算

SISL 的设计灵感来自加拿大 MC2 模式,而这套方案如今用在 ECMWF 的 IFS、加拿大的 GEM、英国气象局的 UM。作者特意点这个渊源,是为了说明它不是玩具方案。

作者强调这两个核不是可互换的实现,也没有普适的性能排序。它们是为互补区间设计的,用哪一个取决于分辨率、可达到的时间步、求解器收敛性、rollout 长度和可用加速器内存。

算例清单

论文第 3 节的六个算例,每个都同时做两件事:正向验证 + 一个不同的可微能力演示。Table 1 给出全部配置。

算例 · 网格 · Δt / 子步 · 终止时刻 · 演示什么

上升热气泡 · 200×3×200,Δ=50/50 m · 1 s / 12 · 1000 s · 浮力驱动非静力流;两核互比 + 收敛阶

Schär 山岳波(正向) · 200×3×50,Δ=500/400 m · 2 s / 6 · 7200 s · 两核与基准比

Schär 山岳波(优化地形) · 300×3×50,Δ=500/400 m · 4 s / 6 · 1800 s · 对地形高度求梯度并优化地形

示踪剂源反演 · 50×50×20,Δ=500/250 m · 2 s / 6 · 1200 s · 对初始示踪剂分布求梯度,反推源位置与强度

三维飑线(正向) · 480×320×96,Δ=312.5/250 m · 0.375 s / 2 · 7200 s · 对流允许尺度基准;Kessler 微物理

三维飑线(敏感度) · 150×100×24,Δ=1000/1000 m · 0.5 s / 2 · 6750 s · 对初始干位温、水汽求 adjoint 敏感度

Wreckhouse 飓风 · 200×200×40,Δ=2000/350 m · 5 s / 4 · 9 h · ERA5 驱动,用 adjoint 构造热力扰动

NEUVE 坐标 · 32×12×16,Δ=500/1000 m · 0.8 s / 20 · 400 s · 对离散方案本身求梯度

Schär 山岳波和飑线各占两列,因为它们先正向积分、再在另一套配置下微分。

验证一:两套核差多少

上升热气泡。 背景虚位温恒为 300 K,扰动是余弦平方型,Δθ = 2 K,半径 R = 1500 m,中心在 (0, 2 km)。三个均质网格格子的准二维形式,跑在完整三维网格上,y 方向 Davies 松弛关闭(孤立理想算例没有外加边界强迫)。

两核的正异常质心相差 5.3 m,垂直展宽相差 3.0 m,积分正温度异常相差约 0.5%。整个扰动场的余弦相似度 0.9953,相对 L₂ 差异 9.69%。差异集中在热边界附近的窄带,反映的是界面位置和幅值的小偏移,峰值幅度比 1.142;主体位置、展宽、异常含量和涡对几何都吻合得不错。左右对称性误差低于 2.5×10⁻⁷ K,相对质量变化低于 5.4×10⁻⁷。

上升热气泡算例:两套核的结果与差异
上升热气泡算例:两套核的结果与差异

论文图 2:t = 1000 s 时准二维上升热气泡在中心 x–z 平面上的虚位温扰动。(a) SISL 核,(b) Split-Explicit 核,(c) 两者的有符号差异。两套核都给出了预期的蘑菇状顶盖和两侧对称反向旋涡对;(c) 中差异集中在热边界附近的窄带,主体结构几乎重合。

山岳波。 地形是解析式:h(x) = h₀ exp[−(x/a)²] cos²(πx/λ),h₀ = 250 m,a = 5 km,λ = 4 km。背景风 U = 10 m s⁻¹,N = 0.01 s⁻¹。

这里有个容易被忽略的处理:全球初始条件会在下坡侧激发出一个启动瞬变,污染稳定驻波。作者把比较区限定在 |x| ≤ 20 km——既排除启动瞬变,又远离约 |x| = 46 km 开始的侧边界松弛带。

在这个区间,两套解的相对 L₂ 范数差 3.12%,余弦相似度 0.999575。垂直廓线差 0.89%,余弦相似度 0.999973。

山岳波算例:两套核的垂直速度与差异
山岳波算例:两套核的垂直速度与差异

论文图 3:t = 7200 s 时准二维 Schär 山岳波算例在中心 x–z 平面上的垂直速度。(a) SISL 核,(b) Split-Explicit 核,(c) 两者差异。显示区域限定在 |x| ≤ 20 km,用来隔离已经建立的地形强迫波、排除下风向启动瞬变。

梯度实验一:用 20 次迭代设计一座山

目标不是模拟山岳波,而是反过来问山该怎么长——这直接对应山岳波、地形降水和严重下坡风的预报问题。

地形用 8 个高斯径向基函数(高斯 RBF)线性组合表示,标准差固定 σ = 2.5 km,中心 μᵢ 在 x = −15 km 到 10 km 之间均匀分布。关键约束是严格的质量预算:振幅之和被限制为 A_tot = 1500 m,否则优化器会直接堆一座任意高的山。每个振幅通过归一化 sigmoid 从无界优化变量映射得到,保证非负且总和固定。

目标函数是目标区域 x ∈ [15, 25] km、z ∈ [3.2, 8.0] km 内垂直速度的平方和,在终止时刻算。

优化设置:20 次 Adam 迭代,梯度裁剪到全局范数 1.0,学习率从 0.5 余弦衰减到 0.005。

结果:优化出的地形在下游目标框内给出的能量比随机质量分布高出一个到两个数量级。 从一个相对平坦、分散的初始猜测出发,优化器系统性地把可用质量集中成一道陡峭山脊,峰值在 x = 10 km 附近。

机制上,最优地形用它有限的物质让向上传播的重力波相位对齐,使它们在下风向目标框内建设性干涉、峰值叠加。

优化地形与三个随机地形的垂直速度
优化地形与三个随机地形的垂直速度

论文图 4:t = 30 min 时的垂直速度。左上为优化地形,其余三幅为在相同总质量预算下抽取的随机地形。黑色虚线矩形标出目标框。随机地形在框内几乎没有能量,而优化地形在框内形成集中的上升区。

梯度实验二:从 6 个塔的浓度反推污染源

地形是三维两峰:主峰 H₀ = 1.8 km 位于 (−6, −2) km,次峰 H₁ = 1.0 km 位于 (2, 4) km。背景风既切变又转向,西风分量从地面 5 m s⁻¹ 增至模式顶 15 m s⁻¹,南风分量 v_max = 4 m s⁻¹。

被动示踪剂初始化为三维高斯羽流,固定展宽 σ_x = σ_y = 2 km、σ_z = 0.8 km,控制参数是 p = (x_s, y_s, z_s, A) 四个量——源的三维位置和幅值。

观测是"合成"的:用真值积分一次,在 6 个监测塔最近网格点上记录浓度。这 6 个塔的位置是 (−2,−2,1)、(−2,2,1)、(0,0,1.5)、(2,−1,1.5)、(2,1,2)、(3,0,2) km。目标函数是模拟与目标浓度历史之间的均方失配。真实释放历史完全未知,反演要对稀疏且不规则的观测做正则化。

优化:Adam 100 步,学习率从 0.25 余弦衰减到其 2%,每次更新把 x_s、y_s 裁剪到 [−20, 20] km,z_s ≥ 0.1 km,A ≥ 0。

初始猜测故意放错:p⁰ = (−4 km, 3 km, 3.5 km, 2)——水平位置错、比真实高度高 1.5 km、幅值只有真值的五分之一。

真实示踪剂羽流在三维地形上的演化
真实示踪剂羽流在三维地形上的演化

论文图 5:真实示踪剂羽流在三维复杂地形上的时间演化,从 t = 5 min 到 t = 20 min。上排是 z ≈ 2 km 的水平 x–y 切片,下排是沿 y = −1.5 km 平面的垂直 x–z 切片。绿色十字标出真实源位置。羽流被地形抬升和拉伸,这是仅凭塔上浓度难以直接反推源位置的典型情形。

优化器在水平面和垂直面内的轨迹
优化器在水平面和垂直面内的轨迹

论文图 6:优化器轨迹。左为水平 xy 平面,右为垂直 xz 平面。轨迹(紫)从初始猜测(橙)出发,收敛到真实源(绿色十字)附近,灰色等高线和剖面轮廓表示地形。轨迹从错误位置单调移动到真实源,说明在这个算例里目标函数不存在明显的多解陷阱。

梯度实验三:飑线的降水对环境场有多敏感

这是最有气象使用价值的一块。论点是:降水的敏感度应该直接当成风暴机制的诊断量用——问"这场风暴的降水实际在响应环境的哪些特征",而不是靠一堆扰动集合去反推。

正向算例是标准 Weisman-Klemp 风暴:θ₀ = 300 K,θ_tr = 343 K,T_tr = 213 K,对流层顶 12 km 以上近似等温;q_v0 = 14 g kg⁻¹,饱和混合比用 Tetens 公式算。初始热泡 Δθ = 3 K,中心 (0, 0, 2 km),水平半径 10 km、垂直半径 2 km。背景风在最低 2.5 km 内增加 12 m s⁻¹(u_s = −12 m s⁻¹,du/dz = 4.8×10⁻³ s⁻¹)。

积分到 7200 s:云盾发展出宽阔的砧,雨水仍集中在主要对流核下方的局地柱体里。地面累积形成沿风暴移动方向拉长的多峰带,最大 141.1 mm。

Weisman-Klemp 三维飑线的演化
Weisman-Klemp 三维飑线的演化

论文图 7:三维 Weisman-Klemp 风暴的演化。上图为 t = 1800 s 和 7200 s 的云水与雨水等值面,以及干位温扰动 θ_d′ = θ_d − θ̄_d(z) 的水平与垂直切片(θ̄_d(z) 是初始水平均匀的环境廓线);云水画在 q_c = 10⁻⁵ kg kg⁻¹,雨水画在 q_r = 10⁻⁴ kg kg⁻¹。下图为 t = 7200 s 的地面累积降水。

敏感度实验另起一个更粗的积分(150×100×24,Δ = 1 km),从静止 spin-up 到 t = 6600 s 得到可比的成熟风暴。目标函数是随后 150 s 窗口内、固定风暴覆盖区内的平均雨率。覆盖区 Ω_s 在未扰动轨迹上诊断一次(取超过窗口累积最大值 10% 的内部列,外扩四格halo),之后微分过程中保持不变。

成熟风暴降水的 adjoint 敏感度
成熟风暴降水的 adjoint 敏感度

论文图 8:成熟风暴降水的 adjoint 敏感度。上图为目标窗口累积降水,以及对边界层干位温和水汽扰动的线性雨率响应。黑色等值线表示固定的风暴覆盖区,虚线标出垂直剖面位置。下图为穿过最强降水柱的敏感度剖面,橙色和黑色等值线分别对应 q_c = 0.1 g kg⁻¹ 和 w = 5 m s⁻¹。可微动力核在这里给出的不是集合统计,而是一次反向传播就能拿到的逐点梯度场。

梯度实验四:真实个例上把风速"抬" 2.4%

这是论文里最接近业务应用的一块:不再用理想算例,而是拿 ERA5 实际驱动的真实风暴,用 adjoint 找出一个能加强预报的风。

设置:2025 年 2 月 14 日的逐小时 ERA5 场初始化并作侧边界强迫,域中心在 Wreckhouse,地形取自 GEBCO 高程、轻微平滑后并入侧边界海绵层。不开启任何次网格物理参数化——响应完全来自已解析的动力过程,因此干净。前 6 小时作为动力调整的基线态,基线预报里存在一个紧贴陡崖的低空急流,伴随地形上明显偏移的等熵面。

在 t = 6 h 施加一个加性虚位温扰动 Δθ_v,Exner 压力保持固定,密度用 θ_v + Δθ_v 重新诊断,保证热力一致;不直接扰动压力或风。目标诊断量是 Wreckhouse 中心 3×3 网格点水平风速的瞬时空间平均,取最低地形跟随质量层(k = 0)——该层在中心点约高于当地地形 174 m、高于海平面 219 m。几何高度随地形水平变化,不做向地面或固定离地高度的插值。

这里有一个作者主动交代的技术细节,值得记:求梯度用的控制模型比正向评估用了更强的稳定化(c_h = c_div = 0.20 对 0.05),以压制一小时反向轨迹上的网格尺度噪声;算出来的扰动再拿到弱稳定化配置里评估。这是一个刻意的交叉检验——问的是"从稳定化模型得到的敏感度,在标准设置下还能不能兑现预测的响应"。

梯度 $g = ∇_{\Deltaθ_v} J$ 经过一个 6 km × 6 km × 1.05 km 的三维平滑算子并屏蔽侧边界松弛区,再归一化施加上限 1 K 的扰动。

结果:

量 · 控制配置 · 评估配置

3×3 格点平均风速 · 164.9 km h⁻¹ · 未扰动 166.5 → 扰动后 **170.5 km h⁻¹**

中心格点风速 · — · 165.6 → 169.7 km h⁻¹

放大 4.0 km h⁻¹,约 2.4%。这个幅度看着小,但对照扰动的量级就很清楚:扰动点值上限 1 K,全 domain 均方根幅值只有 0.015 K——几乎全域均匀加热的反面。

一维单侧 Taylor 检验(沿归一化方向施加 0.01 K 扰动)给出的相对残差 2.92%,确认控制模型的梯度能预测非线性响应的符号和局部幅值。

adjoint 找到的是一个上游海洋气流里的紧致有符号热偶极,而不是空间均匀增温。垂直剖面上,扰动与基线表现为相邻的加强与减弱区,目标点落在加强区内。

Wreckhouse 下坡风事件:adjoint 指向的修改
Wreckhouse 下坡风事件:adjoint 指向的修改

论文图 9:adjoint 指向的 Wreckhouse 下坡风事件修改。左上为 Wreckhouse 中心格点最低地形跟随质量层(约当地地形上方 174 m)的基线与不利(扰动后)风速,点线和虚线分别标出 t = 6 h 引入热力控制与 t = 7 h 的验证时刻。右上为 t = 6 h 最低质量层上的虚位温扰动,附地形等值线、基线低层风与下方面板所用剖面,星号为 Wreckhouse 区域。左下为 t = 7 h 穿过陡崖的基线风速与虚位温等值线;右下为同一时刻的不利减基线风速响应,并叠加基线等熵面作参照。

梯度实验五:优化模式自己的离散方案

这一节最能说明"可微"这件事的分量——控制变量是垂直坐标变换本身。

原理是一个自洽的诊断量:静力平衡的大气静止状态是连续非静力方程的精确解,但在地形跟随坐标里,水平压力梯度力包含两个相互抵消的贡献,垂直压力梯度与浮力项也要保持平衡。它们的离散近似在斜坐标面上不必精确抵消,残余加速度就产生虚假运动。从静止初态出发的动能,因此提供了一个自监督的离散平衡误差度量。

实验设计有一个精妙处:Split-Explicit 核把完整热力状态分解成一个固定的静力参考廓线加扰动,作者故意把参考的层结设成测试大气的两倍(N_ref = 0.01 s⁻¹ 对测试的 0.02 s⁻¹)。这样初始化时位温和 Exner 压力扰动都非零,尽管完整大气是平衡静止的。理由作者写得很清楚:如果用测试大气本身作参考,扰动会消失,实验就只检验模型参考态里内建的平衡。两倍对比提供了一个困难但稳定的离散良平衡性检验,而且反映了正常预报情形——一个固定的一维参考廓线不可能匹配每一根水平变化、不断演变的 atmospheric column。

NEUVE 坐标的构造:把地形衰减函数 B_φ 建成一个严格正的层密度 ρ_φ 的归一化积分,单调性、B_φ(0) = 1、B_φ(1) = 0 都由构造保证,对任意参数值都成立。ρ_φ 的正性通过对数参数化施加(对数无约束,展开在钳位开均匀三次 B-spline 基 N_{k,3} 上)。基函数二阶连续可微且局部支撑,两个积分用固定的 12 点 Gauss-Legendre 求积计算,替代原 NEUVE 表述里的离散累积积分。

样条系数水平变化,一个全局共享廓线加上地形驱动的局地修正:

目标函数是内部区域的时均、体积加权虚假比动能,加一个屏障项防止层变得任意薄(满足最小层厚约束时屏障为零,对更薄的层二次惩罚)。

训练:10 个独立生成的三维地形,随机位置、取向和各向异性,横跨多种水平尺度。先用一个便宜的几何二分法把全局廓线设为"训练集最小层厚仍高于 175 m 的最强均匀衰减",这一步不用 SLEVE 参数、也不用动力学损失信息。然后 Adam 联合优化 60 个 epoch,梯度裁剪到全局范数 1.0,每次更新在十个地形上取平均,学习率分别是 10⁻³ 和 3×10⁻³ 余弦衰减。

对照是在同一组十地形训练集上调的单组分解析 SLEVE 坐标,两个可调参数(衰减尺度 λ 和指数 h),受同样的最小层厚约束。作者主动区分了这和 Schär et al. (2002) 更一般的双组分 SLEVE 构造——后者先把地形分解成大、小尺度分量再各给一个衰减廓线。这里用的单分量形式提供了一个定义明确的双参数解析对照物,作用在未分解的同一地形场上。调参方式是:几何搜索先为十个候选衰减尺度各自确定最强可行指数,再只对可行边界候选做动力学评估。所以最终 SLEVE 配置是训练集调优的解析基线,不是默认参数选择。

三种坐标在三个未见过地形上的界面分布
三种坐标在三个未见过地形上的界面分布

论文图 10:Gal-Chen、训练集调优 SLEVE 与 NEUVE 三种坐标在三个未见过地形上的对比。上排把最低 6 km 内的坐标界面叠在一起,这一带是地形引起的变形集中的地方;下排给出 NEUVE − SLEVE 的界面位移。两类地形衰减坐标都比 Gal-Chen 拉平更快,且二者的大尺度几何很接近——NEUVE 的差别体现在随地形结构局地重分布的变形上。

累积虚假动能对比
累积虚假动能对比

论文图 11:三个代表性未见过地形上,相对训练集调优 SLEVE 的累积虚假比动能。上排对数坐标面板显示 Gal-Chen 大得多的响应(深蓝,虚线远高于基准),下排面板分辨 SLEVE 与 NEUVE 的差别(浅蓝量级)。SLEVE 定义为参考基准,低于 1 表示优于 SLEVE。NEUVE 在整个积分过程中都在基准以下,不是在最后时刻靠一个瞬态赢下来的。

结果数字:

指标 · 数值

三个展示地形上积分末的累积 NEUVE/SLEVE 比 · 0.75 / 0.72 / 0.79

十个未见过地形上 Gal-Chen 平均虚假动能 · 7.88×10⁻¹ m² s⁻²

训练集调优 SLEVE · 1.23×10⁻¹ m² s⁻²

NEUVE · **9.30×10⁻² m² s⁻²**

相对 Gal-Chen 的逐地形平均降幅 · 86.99%

相对调优 SLEVE 的降幅 · 23.88%

相对 SLEVE 改善的配对 bootstrap 95% 置信区间 · 21.95%–25.92%

十个留出地形上 NEUVE 优于 SLEVE 的个数 · 10 / 10

论文图 11:三个代表性未见过地形上,相对训练集调优 SLEVE 的累积虚假比动能。上排对数坐标面板显示 Gal-Chen 大得多的响应(深蓝,虚线远高于基准),下排面板分辨 SLEVE 与 NEUVE 的差别(浅蓝量级)。SLEVE 定义为参考基准,低于 1 表示优于 SLEVE。NEUVE 在整个积分过程中都在基准以下,不是在最后时刻靠一个瞬态赢下来的。

两个关键反直觉点:

NEUVE 的优势不来自大尺度几何形态的根本差异。 两者界面很接近,但地形相关的界面位移可达约 200 m,重新分配了变形和层压缩。作者的解读是:最小化离散压力梯度误差,不只取决于高层快速拉平,还取决于整根柱子里坡度、曲率和层压缩分布的相当微妙的变化。 学习到的曲面在响应地形高度、坡度、曲率变化的同时保持光滑。

改进不靠"允许更薄的层"。 两个优化坐标在训练集上都受 150 m 最小层厚约束。在更极端的留出地形上,最小层厚降到调优 SLEVE 的 81.0 m 和 NEUVE 的 78.7 m——NEUVE 的层更薄一点,但所有映射保持单调、所有积分保持稳定,且 NEUVE 在整个留出集上普遍保持与 SLEVE 可比的层厚,同时高层坡度和曲率更小。

作者还做了分层结稳健性检验:把三个固定坐标(不重训)用于 N_bv ∈ {0.01, 0.0125, 0.015, 0.0175, 0.02} s⁻¹ 的静力静止态,保留 N_ref = 0.01 s⁻¹。

也就是说,学到的优势对测试族里所有偏离参考的静力平衡都成立,不是只在训练时用的两倍对比下才成立。

附录 A.1:正向和反向的收敛性

这一节回答一个质疑:可微的梯度到底准不准,还是只是"能算出来"。

正向自收敛。 上升热气泡在 10 km × 10 km 的笛卡尔 x–z 域、y 方向三个网格格子上积分到 t_end = 24 s。网格 Δx = Δz = 250、125、62.5、31.25 m,大时间步按比例从 1 s 细化到 0.125 s,Split-Explicit 全程 12 个声学子步。细网格解先限制到次粗网格再在内部子域上算差,排除邻接 x 和 z 边界的 1.5 km 条带以减少数值边界处理的影响;排除宽度以物理单位固定,因此各分辨率比较的是同一块物理区域——报告的收敛率刻画的是内部收敛。

结论:两套独立构造的核在速度、Exner 压力和虚位温上都趋近二阶收敛。

跨核收敛。 同一套四次积分在 t = 24 s 上比较 SISL 与 Split-Explicit,用各公共网格上的内部 L₂ 差异。差异随细化系统下降,最细一级的收敛率在所有预报量上落在 1.79–1.99。作者给这句结果加了解释:因为两套核对轨迹、声学模和压力耦合的处理方式根本不同,它们彼此的收敛可以作为额外的核健康度检查——说明两种离散在逼近同一个被解析的解,而不只是各自内部收敛。

离散伴随的验证。 第二个上升热泡实验,固定 12×3×8 网格,Δx = Δz = 375 m,终止时刻 t = 2 s。标量目标函数是体积积分的上升气流动能加虚位温扰动方差,离散梯度对初始 u、w、Exner 压力扰动和格点中心虚位温四个场求。两核都用中心、无阻尼的加密族,Δt ∈ {0.5, 0.25, 0.125, 0.0625} s。

检查 · 结果

一阶方向 Taylor 余项 R₁(扣除 adjoint 方向导数后) · 二次衰减——直接验证每个离散反向实现

Split-Explicit 四个梯度场的时间自收敛率 · 1.998–2.023

SISL 四个梯度场的时间自收敛率 · 1.911–2.008

30 s 非线性热实验,独立生成的虚位温敏感度场余弦相似度 · 0.99955

这一节的含义值得单独记住:伴随场是独立生成的,与正向解的收敛性不是同一回事。 一个模式在预报上误差很小,不代表它的反演梯度也可靠。

附录 A.2:CPU-GPU 扩展,以及作者为什么不肯排名

作者先立好界限:这是一次有限硬件实验,用来展示两套 JAX 核都能在加速器上执行,不是建立优化性能基线,也不是给算法排名。所以表里没有别的模式对照。

硬件:Intel Core i9-14900K(32 逻辑核)与 NVIDIA GeForce RTX 4090(24 GB)。干燥中性 10×10×5 km 域,除一个 z = 1.5 km 中心、半径 1 km 的 2 K 球形热泡外初始静止。物理、侧边界强迫、水平扩散和散度阻尼全部关闭,4 km 以上保留一个弱 w 海绵。网格按 N³ 加密,Δx = Δy = 10 km/N,Δz = 5 km/N,外步长正比于 Δx。Split-Explicit 用 8 个声学子步,SISL 用十倍长的外步长、GMRES 相对容差 10⁻⁴、最多十次迭代。全部计时用 32 位算术,编译与预热不计入,重复评估同步执行。

负载 · CPU · GPU · 加速比

Split-Explicit,96³,四步正向 · 1.11 s · 17.9 ms · ≈ 60×

Split-Explicit,96³,值+梯度 · 3.90 s · 62.3 ms · ≈ 60×

SISL,96³,四步正向 · 54.4 s · 1.42 s · 38×

SISL,96³,值+梯度 · 87.9 s · 2.13 s · 41×

吞吐上,Split-Explicit 正向在 48³ 和 64³ 达到 193 和 230 百万 cell-steps 每秒,值+梯度在 48³ 是 71 百万,96³ 上保持在 57–67 百万。SISL 的 GPU 吞吐在 64³ 达到 4.08 百万(正向)和 2.36 百万(值+梯度)cell-steps 每秒。

最有信息量的是反向开销的差别:

作者立刻加了一句防止过度解读的话:这不代表 SISL 普遍更快——它的正向前向步本身贵得多,而且更大的时间步推进的是不同长度的物理时间。准确的表述是:Split-Explicit 核高效映射到规则加速器核函数,SISL 则用一个更昂贵的隐式步换取更长的时间步和相对温和的反向增量成本。

64³ 上 SISL 的加速比看起来更高(约 111 和 94),但作者明确说不强调这些峰值比,因为 CPU 扩展性和小问题开销在加密序列上是变化的。绝对耗时对编译器与 JAX 版本、CPU 线程设置、加速器代际、精度、Krylov 收敛性和问题规模都敏感。这些比值只解释为"这套实现在这些代表性负载、在所述设备上取得的加速",不是与其他大气模式的比较、不是硬件无关的吞吐主张、也不是精度归一化的算法定量排名。

作者承认的限制

实验都是理想化形式。 论文自己写明:第 3 节的算例代表既有科学问题中理想化的形态——

物理包是紧凑的,不是完整的。 Suêtes 当前只包含基本湿过程和湍流的紧凑表示,一个完整的大气模式还需要边界层、辐射传输、微物理和地表通量的可微参数化。作者列了三条下一步,都属于扩展系统本身而非应用:

作者还说明现有算例没有开启次网格物理参数化(Wreckhouse 个例明确写"响应完全来自已解析的动力"),并且这些应用本身"超出本工作范围"。

一句话判断

这不是"又一个 JAX 天气模式",而是把梯度做成了模式的输出。真正值得单独记住的是 3.6 节那条结论:可微动力核可以让模式优化自己的离散方案,而 NEUVE 对调优 SLEVE 的 23.88% 改善在 10/10 个留出地形上都成立、对分层结稳健、且不靠允许更薄的层——这是能拿去和业务模式讨论的硬结果。做 LUOS / 四维同化、反问题或地形敏感度诊断的人可以顺着这条线看。

在观潮中打开(含配图与交互)→