利用 NetCDF4 数据计算海峡断面通量

本文整理一次 PlioMIP3 相关数据处理中的断面通量计算流程。核心问题是：如何从 POP 模型输出的 NetCDF4 数据中，截取指定海峡断面，并计算体积输运、海洋热输运以及沿深度累计的物理量。

这篇文章主要参考李振学长的代码 ITF.ipynb、以及NCL 原脚本，并通过对比脚本中验证结果。总结了“数据维度、物理公式、单位换算、NetCDF 输出和误差来源”的工作流程。

背景：为什么要算海峡断面通量

PlioMIP3 的气候模拟实验关注古地理边界条件、温室气体强迫和地球系统反馈对上新世气候的影响。在海洋环流问题中，海峡开闭或通道深度变化会改变不同海盆之间的水体交换，例如 ITF（Indonesian Throughflow）和 Panama 相关断面。

对这些断面的诊断通常至少需要两个量：

• 体积输运（Volume Transport）：某一断面单位时间通过的水体体积，常用单位是 Sv，其中 $1\mathrm{Sv}={10}^6{\mathrm{m}}^3{\mathrm{s}}^{-1}$。
• 海洋热输运（Ocean Heat Transport, OHT）：水体携带的热量通量，常用单位是 PW，其中 $1\mathrm{PW}={10}^{15}\mathrm{W}$。

在当前工程里，主要比较三套实现：

文件	作用
ITF.ipynb	单断面诊断和绘图，用于理解计算过程
step2_monthly_ocn_variables_ITF.py	Python 批处理脚本，用 NetCDF4 生成逐月和气候态输出
NCL脚本/step2_monthly_ocn_varibles_ITF_by_LiZhen_kimi*.ncl	原始 NCL 版本，是 Python 重写的参照

数据结构：从 NetCDF 变量到断面

计算使用的关键输入有三类 NetCDF 文件：

文件类型	主要变量	含义
*-ocn-UVEL-monthly-*.nc	UVEL	纬向流速，单位通常为 cm/s
*-ocn-TEMP-monthly-*.nc	TEMP	海水位温，单位为 ${}^{\circ }\mathrm{C}$
*.pop.gx3v7_gridinfo.nc	DYU	POP 网格经向宽度，原始单位为 cm

脚本中读取的速度和温度一般具有类似维度：

(time, depth, lat, lon)

计算海峡断面时，不是对整个三维海洋场积分，而是在固定经度上截取一条沿纬向展开的断面。因此每个断面需要三类定位信息：

• 固定经度：用于找到对应的经度列；
• 起止纬度：用于确定断面覆盖的纬向网格点；
• 网格索引辅助量：用于从 POP 的曲线网格变量中取出 ULONG、ULAT 和 DYU。

当前脚本中定义了两个断面：

channel_names = ["ITF_W", "Panama"]

dx_indices = np.array([52, 88])
dy_indices = np.array([43, 67])

其中 ITF_W 对应印尼贯穿流西侧断面，Panama 对应巴拿马海道附近断面。实际截取时，脚本使用最近邻方法找到目标经纬度最接近的模型网格点。

def find_nearest(array, value):
    return int(np.argmin(np.abs(array - value)))

这个步骤很简单，但很重要：模型网格点并不一定正好落在目标经纬度上，所以断面定位本质上是“目标地理位置到模型网格”的映射。

1. array - value：让数组中的每一个元素都减去目标值。
2. np.abs(...)：计算上述差值的绝对值。这时候，差值的绝对值越小，就代表原本那个元素距离目标值越近。
3. np.argmin(...)：寻找这个绝对值数组中的最小值，并返回该最小值所在的索引（位置）。
4. int(...)：将 NumPy 返回的索引类型（通常是 np.int64）强制转换为 Python 的原生整数类型 int。

符号和物理量

设深度层索引为 $k$，断面纬向网格点索引为 $j$，时间索引为 $t$。后续计算使用的主要符号如下：

符号	含义	单位
$\mathrm{\Delta }{H}_k$	第 $k$ 层层厚	m
$\mathrm{\Delta }{y}_j$	断面网格点经向宽度	cm $\to$ m
$u_{t,k,j}$	纬向流速 UVEL	cm/s $\to$ m/s
$T_{t,k,j}$	海水位温 TEMP	${}^{\circ }\mathrm{C}$
$\rho$	海水密度	kg m${}^{-3}$
$c_p$	海水定压比热	J kg${}^{-1}$ K${}^{-1}$

在 step2_monthly_ocn_variables_ITF.py 中，物理常数当前设为：

rho = 1025.0
cp  = 3990.0

需要注意， .py / NCL 与 notebook 存在 cp=4178 和 cp=3990 的差异。热输运与 $c_p$ 成正比，因此不同取值会直接带来系统性比例差异。做跨脚本对比时必须先统一这个常数。

第一步：截取断面

对固定经度断面，脚本先在 ULONG 中找经度索引，再在 ULAT 中找纬度起止索引。由于 Python 切片不包含末端，因此截取纬度范围时需要 +1：

u_sec = UVEL[:, :, idy_lat1:idy_lat2+1, idx_lon]
t_sec = TEMP[:, :, idy_lat1:idy_lat2+1, idx_lon]

截取后，断面变量维度变为：

(time, depth, section_lat)

也就是说，每个月、每个深度层，都有一条沿断面方向的速度和温度序列。

第二步：计算断面面积元

体积通量需要面积元。对于第 $k$ 层、第 $j$ 个断面网格点，面积元为：

$$\mathrm{d}{A}_{k,j}=\mathrm{\Delta }{H}_k\cdot \mathrm{\Delta }{y}_j$$

其中 $\mathrm{\Delta }{H}_k$ 是层厚，$\mathrm{\Delta }{y}_j$ 是网格宽度。脚本中 DYU 原始单位是 cm，因此需要先除以 100 转为 m：

dy = ds_grid.variables["DYU"][:, dy_idx] / 100.0
dwidth = dy[idy_lat1:idy_lat2+1]
ds_layer = dH[:, np.newaxis] * dwidth[np.newaxis, :]

得到的 ds_layer 维度为：

(depth, section_lat)

它描述了每个深度层、每个断面网格点对应的面积。

第三步：计算体积输运

逐层体积通量是面积元与法向流速的积分：

$${\mathrm{\Phi }}_{t,k}=\sum _j\mathrm{d}{A}_{k,j}{u}_{t,k,j}\times {10}^{-2}[{\mathrm{m}}^3{\mathrm{s}}^{-1}]$$

这里的 ${10}^{-2}$ 来自速度单位换算：UVEL 是 cm/s，而面积元已经是 m${}^2$，所以速度必须转成 m/s。

Python 脚本中的实现是：

volume_flux = np.sum(
    ds_layer[np.newaxis, :, :] * u_sec,
    axis=2
) * 1.0e-2

volume[:, :, ch] = volume_flux / 1.0e6

其中：

• ds_layer[np.newaxis, :, :] 把面积元扩展到时间维；
• axis=2 表示沿断面纬向积分；
• * 1.0e-2 将流速从 cm/s 转为 m/s；
• / 1.0e6 将 m${}^3$/s 转为 Sv。

最终得到的 VOLUME 维度为：

(time, depth, channel)

第四步：计算断面平均温度

热输运需要温度。当前脚本采用断面简单平均温度：

$${\overline{T}}_{t,k}=\frac{1}{N_j}\sum _j{T}_{t,k,j}$$

对应代码为：

temp_avg = np.mean(temp_sec, axis=2)

这里的 axis=2 仍然是断面纬向。得到的 TEMPAVG 维度是：

(time, depth, channel)

严格来说，温度平均也可以设计为面积权重或流量权重平均。当前脚本使用简单平均，主要是为了与 notebook 中的计算口径保持一致。后续如果要做更严格的热输运诊断，需要明确温度代表的是断面平均、流量加权平均，还是逐网格点 $uT$ 直接积分。

第五步：计算海洋热输运

在当前口径下，逐层海洋热输运写作：

$${\mathrm{OHT}}_{t,k}=\frac{\rho c_p{\mathrm{\Phi }}_{t,k}{\overline{T}}_{t,k}}{{10}^{15}}[\mathrm{PW}]$$

对应代码为：

oht[:, :, ch] = rho * cp * volume_flux * temp_avg / 1.0e15

这一步把体积通量、温度、密度和比热容组合起来，得到单位时间通过断面的热量。由于 $\mathrm{\Phi }$ 已经是 m${}^3$/s，乘以 $\rho$ 后得到 kg/s，再乘 $c_{p}T$ 得到 W，最后除以 ${10}^{15}$ 得到 PW。

第六步：沿深度累计

很多时候我们不仅关心每一层的通量，还关心从表层累积到某一深度的通量剖面。这个操作就是沿深度方向做前缀和：

$${\mathrm{cVOLUME}}_{t,k}=\sum _{k^{\prime }=0}^k{\mathrm{VOLUME}}_{t,k^{\prime }}$$$${\mathrm{cOHT}}_{t,k}=\sum _{k^{\prime }=0}^k{\mathrm{OHT}}_{t,k^{\prime }}$$

Python 中直接用 np.cumsum：

cvolume = np.cumsum(volume, axis=1)
coht    = np.cumsum(oht,    axis=1)

这里 axis=1 是深度维。

第七步：写入 NetCDF4 文件

批处理脚本会为每个实验、每个时间块生成中间逐月文件，随后拼接成完整的 monthly 文件，再计算 climatology 文件。

输出变量包括：

变量	含义	单位	维度
VOLUME	逐层体积输运	Sv	time/month, depth, channel
TEMPAVG	断面平均温度	degC	time/month, depth, channel
OHT	逐层海洋热输运	PW	time/month, depth, channel
cvolume	累计体积输运	Sv	time/month, depth, channel
coht	累计海洋热输运	PW	time/month, depth, channel

相比 NCL 原脚本，Python 版本有一个工程上的改进：不再依赖外部 ncrcat 命令拼接 NetCDF，而是在 Python 内部使用 np.concatenate 合并时间维。

all_volume = np.concatenate(all_volume, axis=0)
all_oht    = np.concatenate(all_oht,    axis=0)

这样脚本对运行环境的依赖更少，也更适合在 Python 项目中维护。

第八步：多年同月平均与气候态

气候态文件的目标是得到 12 个月的多年同月平均。假设总时间长度可以整除 12，脚本先把时间维重排为：

(year, 12, depth, channel)

再对年份维求平均：

def clim(arr):
    return arr.reshape(nyears, 12, ndH_file, nch_file).mean(axis=0)

volume_clim  = clim(VOLUME_all)
tempavg_clim = clim(TEMPAVG_all)

一个容易被忽略的问题是：气候态 OHT 到底应该怎样平均？

有两种口径：

$$\overline{U\cdot T}$$

和

$$\bar{U}\cdot \bar{T}$$

二者一般不相等，差异来自协方差项。NCL 原脚本采用的是先逐月算 OHT，再对年份维平均，也就是 $\overline{U\cdot T}$。Python 版本为了与 notebook 口径一致，在 climatology 阶段使用平均后的体积输运和平均温度重新计算：

oht_clim = rho * cp * (volume_clim * 1.0e6) * tempavg_clim / 1.0e15

这相当于 $\bar{U}\cdot \bar{T}$。因此对比 NCL 与 Python 的气候态 OHT 时，不能只看数值差异，还要先确认平均口径是否一致。

与 NCL 对比时发现的单位问题

对比脚本 compare_NCL_Li.py 暴露出一个关键问题：NCL 原脚本中体积输运少乘了一个 1.0e-2，导致 VOLUME 结果整体偏大 100 倍。

NCL 中的相关逻辑是：

dwidth = dy(idy_lat1:idy_lat2) / 100.0
ds_layer = dH(m) * dwidth
volume_flux = sum(ds_layer * u_section(t, m, :))
VOLUME(t, m, ch) = dble2flt(volume_flux) / 1.0e6

这里 dwidth 已经从 cm 转成 m，但 u_section 仍然是 cm/s。面积元是 m${}^2$，流速却是 cm/s，因此 volume_flux 实际上比 m${}^3$/s 大 100 倍。

Python 版本补上了这一步：

volume_flux = np.sum(ds_layer[np.newaxis, :, :] * u_sec, axis=2) * 1.0e-2

对 PIdefault 的检验结果显示，notebook 与 step2 在同一原始数据和相同常数设置下完全一致：

cp=3990, rho=1025
channel=ITF_W, months=120
notebook_total_transport_sv=-14.5356266077
step2_total_transport_sv=-14.5356266077
notebook_total_heat_pw=-0.508581136882
step2_total_heat_mean_u_times_mean_t_pw=-0.508581136882
step2_total_heat_mean_of_monthly_pw=-0.510048246123

其中 mean_u_times_mean_t 与 notebook 一致，而 mean_of_monthly 略有差异，正是因为平均顺序不同。

一条完整计算流程

整体流程可以概括为：

读取 UVEL / TEMP / grid NetCDF 文件
        |
        v
根据目标经纬度定位断面网格点
        |
        v
截取 (time, depth, section_lat) 的 UVEL 和 TEMP
        |
        v
用 dH 和 DYU 构造面积元 dA
        |
        v
VOLUME = sum(dA * UVEL * 1e-2) / 1e6
        |
        v
TEMPAVG = mean(TEMP, section_lat)
        |
        v
OHT = rho * cp * volume_flux * TEMPAVG / 1e15
        |
        v
沿 depth 做 cumsum 得到累计剖面
        |
        v
拼接所有时间块，输出 monthly.nc
        |
        v
reshape(year, 12, depth, channel)，输出 climatology.nc

这个流程中最容易出错的地方不是公式本身，而是单位和维度：

• DYU 是 cm，需要转成 m；
• UVEL 是 cm/s，需要转成 m/s；
• Python 切片不包含末端，纬度范围要 +1；
• time 维必须能正确重排为 (year, 12, ...)；
• 气候态 OHT 的平均口径必须在结果对比前说清楚。

小结

利用 NetCDF4 计算海峡断面通量，本质上是把模型网格上的三维变量投影到一个二维断面，再沿断面方向做积分。对体积输运来说，关键是面积元和速度单位；对热输运来说，关键是体积通量、温度代表性和物理常数。

从这次 PlioMIP3 脚本对比可以得到几个实践经验：

1. 物理公式要和数据单位一起检查，不能只看公式形式。
2. 批处理脚本应尽量保存中间诊断量，例如 TEMPAVG，这样后续可以重算不同口径的 OHT。
3. notebook 适合验证单断面过程，生产脚本适合批量生成 NetCDF，两者最好用同一组原始数据做交叉检验。
4. NCL 到 Python 的重写不只是语法迁移，还包括单位检查、向量化、文件拼接方式和实验参数管理的重构。
5. 对气候态热输运，$\overline{U\cdot T}$ 和 $\bar{U}\cdot \bar{T}$ 是两个不同诊断量，不能混用。

这也是海洋模式后处理里最典型的一类问题：计算看似只是数组乘加，但每个数组背后都有网格、单位、维度和物理含义。只有把这几件事同时对齐，最后得到的 Sv 和 PW 才真正有解释价值。