第一次听说 WRF-CMAQ 是在课题组一次组会上师兄顺手把一张模型流程图甩给我看上面密密麻麻排着 WRF、MCIP、ICON、BCON、CCTM 几个缩写。我当时的第一个反应是这不就是四五个模型串起来跑呗跑通一个总流程就行了吗结果真上手之后才发现这套系统里的每一个模块都有自己的脾气任何一环出了问题最后高时间分辨率的 O3、PM2.5 模拟结果都会变成一堆漂亮但没法用的数字。这篇就来把 WRF-CMAQ 这套空气质量数值模拟系统从头到尾讲清楚适合刚接手课题要做空气质量模拟、读不懂模型文档、或者准备搭建本地模拟平台的读者看完至少能搞明白这条流水线里谁在干什么、你跑 CCTM 之前缺了什么、以及结果离谱的时候该从哪里查起。1. 先看大图WRF-CMAQ 这条模拟流水线里谁在干什么1.1 WRF 的任务把“风、温、湿、雨”算对把天气场交给下游WRFWeather Research and Forecasting本质上是气象模型它在整条模拟链路里的角色不是算污染物而是提供空气质量模型赖以运行的“天气舞台”。风场决定了污染物往哪边输送温度场决定了化学反应速率常数水汽和云的过程决定了气态污染物会不会被降水清除、光解速率有多大行星边界层高度更是直接锁死了污染物的垂直稀释空间。换句话说WRF 输出文件的每一个变量几乎都会被 CMAQ 用到。我在实际跑的过程中发现新手最容易犯的错是把 WRF 当黑箱跑完 WRF 看一眼模拟的 2 米温度大致合理就直接扔给 MCIP 了。这种做法风险很高。因为 WRF 的气象场如果存在系统性偏差比如边界层高度白天明显偏低、夜间风速偏大那么后面 CCTM 算出来的 O3 峰值和 PM2.5 日均值都会跟着歪掉而且歪得非常有规律——白天臭氧偏低、夜间颗粒物偏高或者反过来。所以拿到 wrfout 之后最好先花半天时间和实际观测站点对比几个关键变量10 米风、2 米温度、2 米比湿、降水。日子一天天对不上没关系只要统计指标和观测量级在一个合理范围内就可以进入下一环如果连量级都对不上那问题多半在 WRF 本身得回去修参数化方案或驱动数据而不是硬着头皮把气象场灌进 CMAQ。这一步的省力程度直接决定了你后面查错花费的精力。1.2 CMAQ 的任务把“排放、化学反应、沉降”算清楚回答污染物从哪来到哪去CMAQCommunity Multiscale Air Quality Modeling System是美国环境保护署开发的第三代空气质量模型它和第一代、第二代模型的本质区别是“一个大气”理念用一个 CCTMCMAQ Chemistry-Transport Model核心模块同时处理 O3、PM2.5、酸沉降、能见度和有毒空气污染物。CCTM 拿到 MCIP 处理好的气象输入之后会逐时间步地解算平流输送、湍流扩散、气相化学、气溶胶动力学、液相化学、云化学、干湿沉降以及气溶胶的酸碱平衡。用一句话概括WRF 管“大气物理”CMAQ 管“大气化学”。二者之间靠离线文件传递数据。我不建议一上来就去啃 CMAQ 的源代码——你只要知道 CCTM 里面每一层的计算顺序大致是先做水平平流和垂直对流然后是扩散跟着处理气相化学EBI 或类似数值解法再处理云和液相过程最后算完气溶胶过程之后更新干沉降速度输出浓度场。这个顺序的重要性在于当你去看 debug log 发现某一步数值发散时你能大概猜出是哪一段出了问题。1.3 离线耦合与在线耦合新手先别急着追 WRF-Chem你可能听说过 WRF-Chem 是“在线耦合”模型说它把气象和化学放在同一个时间步里计算还考虑了气溶胶对辐射的反馈。这个方向听起来确实更先进但对刚入门的团队来说我建议先踏踏实实把 WRF-CMAQ 这套离线系统跑通。原因很简单离线耦合让你能把“气象误差”和“化学误差”分开排查。如果模拟结果差你可以分别检验 WRF 单独的气象表现和 CMAQ 的化学表现变量少好归因。而在线模型一旦结果不对气象和化学两个模块互相影响排查链路会非常长。WRF-CMAQ 的“离线”也只是指两个模型分步运行实际在物理机制上CMAQ 读取的气象场仍然是高时间分辨率的通常 WRF 输出频率设置为 1 小时或更密CMAQ 内部会自动做时间插值。所以时间分辨率上并不会明显吃亏。真正吃亏的是缺少气溶胶对辐射的在线反馈但多数城市尺度空气质量模拟场景中这个反馈不是核心矛盾核心矛盾往往还是排放清单的精度。1.4 全套工具链CMAQ 不只是 CCTM 一个可执行文件CMAQ 发布包默认有多个可执行程序我平时最常用的是这么几个模块全称作用MCIPMeteorology-Chemistry Interface Processor把 WRF 气象输出转为 CMAQ 可用格式ICONInitial Conditions生成初始浓度场BCONBoundary Conditions生成边界浓度场CCTMCMAQ Chemistry-Transport Model核心化学传输计算CCTM_ICON / CCTM_BCON内嵌脚本从 CCTM 输出结果生成嵌套域边界 / 初始场很多人刚开始被这一串缩写吓住其实你只需要按顺序理解WRF 输出先给 MCIPICON 和 BCON 分别给出模拟域的初始状态和边界状态CCTM 是最终的计算主体如果你想跑嵌套网格那么父域跑完之后CCTM 的 CONC 文件就会被拿来生成子域的 BCON。这一套就是“用户指南第 13 章”里最核心的流程。把这四个模块的分工刻在脑子里后面看脚本就不会再被绕晕。2. 启动 CCTM 前必须备齐的输入三件套气象、排放、边界场2.1 气象驱动与 WRF 输出文件的验收我通常把 CMAQ 模拟的输入拆成三件套气象场、排放清单、初始和边界条件。这三样东西缺一不可而且每一件都有各自的验收标准。气象场这边除了 WRF 本身的配置还有一个经常被忽略的点WRF 输出文件的时区默认是 UTC而 CMQA 的输出、站点观测、空气质量标准比较都基于本地时间。如果这里弄错你得到的结果会整体平移几个小时臭氧的午后峰值会被算到上午去排查起来非常隐蔽。我的习惯是在写后处理脚本的第一步就把时间轴统一换算成北京时间或目标时区的本地时间并且在下一次后处理里固定用同一套时间戳函数避免每次手动改。另外WRF 输出频率也值得注意。虽然 CMAQ 能在 MCIP 阶段做时间插值但如果 WRF 输出间隔太粗比如 3 小时某些短时强降水、边界层塌陷的过程会被抹平影响夜间 PM2.5 累积的模拟。我的通行配置是外层域输出频率 1 小时内层域 1 小时或 30 分钟如果磁盘够用尽量用 1 小时性价比最高。2.2 排放清单CMAQ 模拟结果准不准一半看这里如果说气象场决定模拟结果“跑向哪里”排放清单就决定模拟结果“起点有多高”。我用 CMAQ 跑了这么多年最深的体会是化学机制选错可以改配置重跑但排放清单给错了后面所有工作都相当于在错误的地基上盖楼。排放清单在我们这边的通行做法分三步选定清单源常规应用里用得最多的是 MEIC 这类区域清单全球尺度可以看 EDGAR。MEIC 的优点是中国区域的排放活动水平和时空分配都比较贴合实际网格分辨率可以根据模拟域选 0.25°或 0.1°直接覆盖到常见嵌套域。把清单数据重格到 CMAQ 网格这一步是最容易出 bug 的。清单原始的网格定义、投影、纬度顺序和 CMAQ 的不一定一致必须用统一的网格描述文件重投影。我建议你写一个重格脚本的单元测试专门输入一个单位脉冲排放检查它到底落到哪个网格就可以快速发现经度方向是否翻转了 180 度、纬度是否倒序。物种映射与时间分配清单里给的是 NOx、SO2、VOCs 等总排放量化学机理需要的是 NO、NO2、各类 VOC 物种和粒径分段的气溶胶前体物。CMAQ 官方发布包里有 profile 工具帮你做映射但不同清单的分子谱不一样你必须确认用的机制和清单谱是配套的比如 MEIC 与 CB6r3 的映射表有一些默认分配直接拿来用没问题但 VOCs 拆分比例会显著影响 O3 生成敏感性因此建议对照文献验证一次。关于清单还有一个经常被忽略的地方时间分配。很多清单给的是年总量你需要按月和按小时生成时间权重的“时间分配因子”。如果缺了这一步模拟结果会产生一种“每天 24 小时排放强度恒定”的假象这会让污染物日变化曲线变得非常钝尤其峰值期明显偏低。我吃过这个亏跑出来的 O3 午后峰值抬不上去找了两天原因最后发现是排放清单压根没做时间分配。2.3 ICON 与 BCON初始条件和边界条件的推荐做法初始条件ICON决定模拟开始时刻整个域内的大气浓度场边界条件BCON决定模拟期间从域外流入的污染物浓度。对区域模拟来说边界条件的影响比初始条件更长尾因为模拟域边缘的空气会随主导风持续吹进模拟区。新手起步阶段可以直接用 CMAQ 发布目录里的 profile 文件跑一个“清洁大气”的边界条件也就是把域外背景 AOTaerosol optical thickness和 O3 季节背景值等设为低值然后让它先模拟 3~5 天作为 spin-up。但如果你的研究区域有明显的境外污染长距离输送影响比如春季沙尘或者夏季台风前期的气团回流单靠 profile 会产生系统性低估。进阶做法是从全球模型里抽取边界场格局常用的是 GEOS-Chem 或 MOZART 输出也有人从 CAMS 全球模式直接下载边界条件重格为 CMAQ BCON 格式。这里需要注意的是边界条件的时变频率一般取 3 小时间隔域外浓度变化剧烈时间隔太大容易让污染物“一块一块”地往域内涌后处理图上看起来像噪声。我自己的偏好是如果磁盘和下载流量允许边界场尽量用 3 小时间隔。初始条件因为 spin-up 的作用影响会衰减得比较快。模拟域越小边界场的影响越小模拟域越大初始场影响越大尤其第一天偏大。因此建议把结果分析从模拟开始后的第二天起算避开初始场的调整期。3. 网格转换的关键抽屉MCIP 在前处理里到底做了什么3.1 MCIP 的数据流从 wrfout 到 CMAQ-ready 的 netCDFMCIP 全称 Meteorology-Chemistry Interface Processor是整个 CMAQ 工具链里最不显眼、却最容易被忽视的一环。它的任务是把 WRF 输出的气象场转换成一整套 CMAQ 需要的中间文件包括 METCRO3D三维温度、气压、水汽、风的分量、METDOT3D风场在网格交错点上的值、METCRO2D二维地表参数和气象量和 METBDY3D边界上的气象场。我第一次跑 MCIP 的时候以为这就是个“格式转换”后来打开输出文件比对了半天才发现它远不止格式转换它要重新计算 CMAQ 自己定义的坐标系里的变量、插值垂直速度、处理干沉降速度所需的植被和土地利用参数、计算云覆盖和降水过程所需的辐射变量。这些变量直接决定了 CCTM 里干沉降算得准不准、云光解算得准不准。MCIP 的配置里最要紧的选项是输出时间步长和网格定义。时间步长好理解一般直接继承 WRF 的输出步长网格定义则需要从 WRF 的 geo_em 文件或 wrfout 里读取网格信息这里很容易出现的一个错误是你自以为设定了一个 12 km 的模拟域结果跑出来一看实际网格间距和排放清单的网格间距对不上。这类问题通常在第一次生成地图叠加图层时爆发——站点落点明显地偏移到道路或行政区边界之外这时候第一反应不是去改画图脚本而应该回到 MCIP 和清单重格环节检查网格参数。3.2 嵌套网格与投影参数改一处就要重查所有下游CMAQ 支持多层嵌套域常见做法是从 36 km 的母域嵌套到 12 km 的中间域再到 4 km 或更细的城区域。嵌套的方向是单向的母域 CCTM 算完输出边界场给子域使用。实际操作时我一般把母域网格设置得足够大让边界尽量远离研究重点区域减少边界条件对兴趣区的影响。比如研究一个城市群最外层域至少要把整个大区域的污染物来源包进去而不仅仅圈住城市本身。投影参数是嵌套设置里最磨人的地方。WRF 和 CMAQ 用的都是 Lambert 投影这个没问题但投影参数两条标准纬线、中心经度、中心纬度必须严格一致。任何一处不一致MCIP 输出的网格定义值会以非常隐蔽的方式偏移。我的排查技巧是把 wrfout 里的网格角和中心点经纬度与 MCIP 输出文件里的坐标做一次算术比对如果偏差超过一个网格边长的一半就能立刻断定投影不匹配。3.3 垂直层数设置为什么 31 层和 35 层模拟 PM2.5 差很多垂直层数的选择很多新手觉得只是“精度问题”的次要因素实际上它对结果的改变比想象中大得多。WRF-CMAQ 的垂直坐标都是地形追随坐标CMAQ 默认的层结构需要和 WRF 的层结构保持一致或经过 MCIP 重映射。你可以在 MCIP 阶段选择保持相同的层数和层顶也可以选择不同的垂直结构。我的建议是至少用 35 层且低层要有足够的层数。原因很简单污染物排放和化学反应都集中在边界层内秋冬季节边界层高度经常只有 300~500 米。如果层数太少每一层的厚度会非常大污染物会被垂直平滑掉地面浓度容易系统性偏低。我们做过一个对比同一个模拟域把垂直层从 31 层加到 35 层低层加密PM2.5 的地面日均值整体抬升了大约 15%而 O3 的日最大 8 小时值变化不大。这说明颗粒物模拟对垂直分辨率更敏感。如果你手里的机器内存有限我不建议在水平分辨率上过度追求 1 km 之类的极端设置先保证垂直层数达标性价比更高。4. CCTM 运行与机制选择化学方案比算力更影响结果4.1 CB6r3、SAPRC07 与 AERO7 怎么选CCTM 的化学机制是整个模型里最有“主观选择”成分的地方。CMAQ 6.x 官方最常用的组合是 CB6r3气相化学碳键机制 AERO7气溶胶机制此外 SAPRC07 也是常用选择。CB6r3 是碳键机制把挥发性有机物按结构片段归类计算快对臭氧模拟稳健是目前大多数基准模拟的首选SAPRC07 的物种拆分更细对 VOCs 敏感分析的细致程度更高但输入清单需要更多物种映射计算量也更大。对你刚开始跑基准模拟的团队我建议直接用官方默认组合 CB6r3 AERO7不要图新鲜乱改。原因有两个一是文档和社区经验多出错容易搜到二是官方默认的 profile 和 namelist 就是围绕这套机制调试过的排错成本最低。等以后你做敏感性分析比如某类 VOCs 对臭氧生成的贡献再考虑切换到 SAPRC 或做机制对比。4.2 CCTM 配置若干项化学机制、气溶胶方案与输出频次CCTM 的 namelistCMAQ 里叫CMAQ_Control或类似配置文件里有几个选项直接影响模拟结果Control CTM_CHEM_MECH cb6r3_ae7_aq ! 气相化学 气溶胶 液相化学的组合 CTM_AEROSYM AERO7 CTM_ABFLUX .TRUE. ! 开启氨气等地表通量 CTM_CLD_PHOT .TRUE. ! 云光解反馈 CTM_WETDEP .TRUE. ! 开启湿沉降 CTM_DRYDEP .TRUE. ! 开启干沉降 /这些开关看起来很基础但每一项改动的物理含义都不同。比如CTM_ABFLUX决定是否从地表计算 NH3 通量如果关掉它农业区的 NH3 排放只能依赖排放清单模拟出来的硝酸铵浓度可能会偏低。云光解反馈也很重要如果开启多云天气下光解速率会改变直接影响 O3 生成。我建议你把各项开关的默认值先全部理解一遍尤其注意CTM_WETDEP和CTM_DRYDEP不要误关否则污染物的沉降清除都没有区域浓度会高到离谱。输出频次方面CCTM 默认可以按小时输出 CONC 文件我通常设置为 1 小时。如果做污染过程分析或轨迹分析可以额外输出更密集的数据但文件体积会成倍增长对磁盘和 I/O 压力不小。实际项目中先保证 1 小时输出绝大多数分析都够用。4.3 MPI 分区那点事NPCOL、NPROW 怎么定log 里看什么CCTM 用 MPI 做并行配置里有两个参数决定进程按二维网格如何分区NPCOL和NPROW。它们的乘积应等于你申请的总 CPU 核数。选择时最好让分区的形状尽量接近模拟域网格的形状比如模拟域是 100 × 80你分 4 × 5 还是 5 × 4负载均衡差别不大但如果是狭长区域分区不当会显著拖慢。我遇到过一种情况CCTM 跑到某个时步之后日志不再增长看起来像卡死实际上是因为某个进程负载特别高而其他进程在等它。解决办法是看CTM_LOG_前缀的日志文件里有没有某一行一直不更新同时用系统监控看各核心的使用率。如果发现个别核 100% 而其他核空闲说明负载不均可以调整NPCOL和NPROW或者检查是不是某一个嵌套域的边界处理导致进程间通信陡增。并行核数不是越多越好。我们实际测试网格规模较小时进程数超过一定阈值后加速比明显放缓反而通信开销变大。常规做法12 km 域、几十万网格规模配 24~48 核已经够用4 km 细网格再考虑上更高并行度。与其堆核数不如先把 I/O 的瓶颈解决掉比如把输出文件写到本地磁盘而不是网络盘。5. 我反复踩过的坑CCTM 跑通后结果离谱的五个排查起点5.1 O3 日变化峰谷错位先去查边界场和 spin-up模拟结果的 O3 如果出现“日出前有一个明显高值午后反而低”这类违背常识的现象不要先怀疑化学机制先查三件事本地时间转换、边界场浓度、spin-up 时长。时间戳问题前面说过是最常见的低级错误边界场给了特别高的地表 O3 背景或者边界层夜间存了大量 O3也会造成夜间高值。有时候你只是新跑了一个嵌套域但 BCON 是从太粗糙的父域转过来的父域在城市附近臭氧偏低那么子域夜间就会出现边界附近污染被“稀释”成高值或低值的不真实结构。spin-up 的时长我建议设置在 3~5 天尤其新启动 ICON 时。我见过有论文用 1 天 spin-up 也觉得够但在复杂下垫面区域化学物种的初始调整期比你想象中长。跑完前几天的结果观察它与后几天是否已经“忘记”初始场如果仍有明显差异就加长 spin-up。5.2 PM2.5 整体偏低排放清单和干沉降都要重新过一遍PM2.5 模拟值整体偏低最常见的原因依次是排放清单总量被低估、清单与机制映射不匹配、干沉降速度过于活跃、垂直层数不足。我的经验是先不要想复杂化学机制问题先跑一个“无湿沉降”测试来检查湿沉降是否在气溶胶上乱扣浓度再检查干沉降模块里气溶胶沉积速度是否因为地表类型数据错配而过快沉降。我们调试过一个案例地表类型文件版本装错导致城区被分成了裸土和草地气溶胶干沉降速度从 0.1 cm/s 跳升到 1 cm/sPM2.5 地面浓度掉了一半。如果排除了这些问题再去核对清单重格脚本里的单位换算。排放清单经常有“吨/年”和“克/月”两种写法CCTM 内部统一用mol/(s·m²)或ug/(m²·s)类型单位换错的 bug 在清单环节里极其常见。我建议你在第一次跑完全流程后挑一个站点把模拟和观测的 PM2.5 按照月尺度做一次平均偏差分析看看是不是全时段系统性偏低如果是那大概率是清单总量问题。5.3 时间轴对不上UTC 与本地时间转换的代价刚才提到过这问题但值得单列一条。CMAQ 所有 IO 时间戳都基于 UTC画图时如果不做转换你会看到臭氧峰值出现在晚上 20 点PM2.5 谷值出现在午夜全反过来。这个坑实在太容易踩了尤其是第一次把模拟结果和国控站点数据叠图时。我的习惯是写一个统一的时间戳函数在后处理脚本里直接调用并明确标注“输入 UTC输出北京时间”每次出图前自动检查时间轴范围是不是覆盖了完整的模拟时段。5.4 网格对不齐MCIP 与排放清单的分辨率必须一致网格错位的问题经常以“浓度图上出现一条条整齐的锯齿线”或者“重点城市旁边出现一个异常高值岛”的形式被发现。这时候你要做的第一件事是检查 MCIP 输出的网格定义和排放清单重格用的网格定义是否完全一致。我们曾有一版清单是从 36 km 重格到 12 km 时没有做面积守恒修正导致某些网格集中的排放量被放大模拟结果在城市周围出现一团团“小太阳”。面积守恒的意义在于当你从粗网格重采样到细网格时不仅要看格点交错还要保证每个格点的排放通量在空间上守恒否则总量就变了。5.5 嵌套域没开反馈D02 结果差得像另一个世界最后一个坑是关于嵌套域的。很多人以为只要跑完 D01 再跑 D02子域结果就一定是 D01 的“高分辨率版本”但 CCTM 的嵌套运行里D01 的结果通过 BCON 传递给 D02D02 的高分辨率结果并不会反过来修正 D01。如果你在 D02 里显著改进了排放清单而 D01 的清单没有同步更新那么在父域边界附近会看到一种奇怪的现象子域内部浓度平滑合理一到边界附近就跳变到另一个量级。这个现象不能简单归结为“边界效应正常”而应检查 D01 和 D02 的排放清单输入是否一致。另一个常见问题是D02 跑完后你没有重新生成新 BCON而是用旧 D01 的输出当边界导致污染物通过子域边界进来之后浓度断层。正确的做法是每一次 D01 重跑后都重新执行一次 BCON 生成脚本再启动 D02。写在最后一条最朴素的入坑建议如果你刚接触 WRF-CMAQ我给的建议是先别急着把模拟域设得又大又细更别急着换高级化学机制。用官方基准配置选一个较小但具有代表性的区域把 WRF、MCIP、ICON、BCON、CCTM 全流程跑通一遍做一个月的模拟拿站点观测做一次粗略对比先确保结果在量级上不离谱。这时候再去调排放清单、调垂直层数、加嵌套域——每一次只改一个变量你会发现问题的归因会清晰很多。我在这些年里踩过的坑十有八九不是 CCTM 算法本身的错而是前处理环节里某个网格、单位、时间轴或开关没有对齐。把这套检查做扎实你的 CMAQ 模拟会省心得多。