1. 边界条件为什么是区域温室气体模拟的“地基”区域气候和空气质量模式圈出一块计算域但大气本身是没有边界的。CO2、CH4这类长寿命温室气体随气团流动从边界外源源不断进入你圈定的区域。模拟域内部排放清单、VPRM植被通量算得再准也只是“内部账”。外部本底浓度一旦给错模拟结果会像被污染的观测记录一样整段整段地偏向错误方向。我帮团队搭建WRF-GHG区域模拟系统时一开始把注意力全放在排放清单和生物圈通量上。直到一次测试我故意用了一套旧年份的全球背景浓度结果模拟出来的CO2时间序列整体比观测高8~10ppm而且这个偏差随主导风向不停地从边界往里渗透。从那之后我彻底明白了一个道理背景场处理特别是边界条件才是区域温室气体模拟最容易被低估的环节。本系列前面几篇已经讲了气象背景场、排放清单这类输入的准备这篇专门聊WRF-GHG-Prepy工具中CAMS-Inversion数据做边界条件时的函数解析。CAMS-Inversion是哥白尼大气监测中心Copernicus Atmosphere Monitoring Service发布的全球反演产品它把地面站、卫星等观测数据同化进全球化学传输模型对温室气体通量进行优化约束再输出一套全球网格化浓度。相比普通全球分析场反演产品背后的通量更接近真实地气交换状态对区域模拟来说这是一套物理上更自洽的“外部大环境”。1.1 区域模拟的“外部记忆”边界条件在模式方程中的角色从数值计算角度讲区域模式解的是局部区域内的输运方程。求解偏微分方程时如果没有边界上的约束条件方程组解不唯一模式根本没法收敛到一个有物理意义的答案。通俗地说你在大气里圈出一块“自留地”但周围的气体每时每刻都在往这块地里流动。边界条件描述的就是“自留地”边缘处浓度怎么变化模式内部再通过平流、湍流混合把这些边界信息传输到全域。WRF-GHG以及与之耦合的VPRM生物圈模块会在模拟域的四个侧边界和顶部动态更新CO2、CH4等示踪物浓度。气象场边界通常用ERA5、GFS等再分析数据逐小时更新。温室气体边界则往往依赖全球大气成分产品比如CAMS-Inversion。这里要特别提醒一句温室气体边界和气象边界是两套完全独立的输入不要混为一谈。很多人把气象边界条件配好就以为完事了结果模式跑出来的CO2边界上每6小时才更新一次而且浓度数值完全不合理整个模拟就废了。在实际运行中边界条件更新频率、浓度数值的合理性、垂直方向上的插值精度会直接影响模拟区域内部的浓度分布。尤其对CO2这种本底浓度已经接近420ppm的组分来说边界浓度偏差5ppm其影响往往比排放清单某些分项误差还大。如果你做的实验是区域本底变化分析、碳源汇反演、或者城市CO2扩散模拟边界条件几乎就是决定性因素之一。1.2 CAMS-Inversion和CAMS分析场有什么区别很多刚接触CAMS数据的朋友会在数据列表里看到两类产品一个是CAMS全球分析场一个是CAMS全球反演产品名字相近用途却大不相同。分析场是连续同化模拟更多反映的是全球模型对大气状态的即时估计优点是有近实时的逐小时场适合做快速预报背景。反演产品则使用传输模型模拟浓度并用观测进行通量反演优化输出结果通常滞后发布时间较长但通量和浓度的可信度更高。从做区域温室气体模拟的角度两套数据都能充当边界条件但开展通量相关研究时我更倾向于用CAMS-Inversion。理由很简单你区域内部的VPRM计算的是光合与呼吸通量排放清单给出的是人为源汇这套内部通量体系需要和区域外的大尺度通量反演结果保持一致。如果外部边界用的是未经通量优化的全球分析浓度内外两套体系在边界面附近容易出现“接缝”导致模拟场在边界处出现明显不连续。下表是我自己常用的产品对比思路对比项CAMS分析场CAMS-Inversion反演产品数据滞后时间近实时一般需数日到数周通量是否优化否是通过观测同化反演浓度连续性较高较高且与优化通量自洽适用场景业务预报、快速模拟碳通量研究、精细区域模拟是否适合做GHG边界可以但需谨慎推荐物理一致性更好我个人的经验是做多天到数月尺度的模拟实验优先选择CAMS-Inversion的浓度场做边界。如果实验要求近实时模拟、来不及等反演产品发布再用分析场过渡但后续对比结果时要留个心眼边界本底的不确定性会传导到内部浓度场。1.3 预处理为什么要做“函数解析”WRF-GHG-Prepy不是直接把CAMS数据拖过来就能跑。CAMS-Inversion发布格式是为了全球模型和再分析体系服务的不会迁就WRF的有限区域网格。CAMS输出的全球规则经纬网格、气压/海拔层次、每日或每3小时时间频率与WRF要求的顺序、嵌套域、η坐标、逐小时边界更新这些约定完全对不上。“函数解析”这个词说白了就是把CAMS数据翻译成WRF边界条件格式的整套逻辑。Prepy工具里一般会有一组函数负责文件路径解析、变量读取、子区域裁剪、水平插值、垂直插值、时间插值最后写入wrfbdy文件。之所以要强调函数化而不是写死一个脚本是因为CAMS产品本身也在迭代变量名从co2换成co2vmr网格分辨率从0.5度变到0.25度时间频率从月均变成日均如果预处理脚本硬编码一个版本升级就要改一遍。合理的设计是对外暴露少量参数配置对内把数据读取、插值、写文件解耦成独立函数每个函数只解决一个问题。这也是我在这篇里要重点做函数解析的原因。理解了函数背后的输入输出、插值逻辑和单位换算你才能在自己的实验中灵活适配不同版本的CAMS数据。否则遇到一次变量名变化就卡壳一次莫名其妙地制造加班。基辛格有句话叫“如果你控制了石油你就控制了所有国家”在WFR预处理里“如果你理解了数据解析函数你就控制了所有格式坑”。2. CAMS-Inversion数据本身长什么样解析前必须看清CAMS-Inversion产品的底层格式并不可怕无非是netCDF或grib。真正麻烦的是不同版本里的命名、单位、维度顺序各不相同。我见过最坑的情况是下载回来的文件里同时包含CO2浓度和通量两类变量名称高度相似稍不注意就拿错变量当边界。做函数解析之前第一步永远是“摸底”——把数据文件的维度、变量、属性彻底摸清。2.1 网格、变量命名与文件格式常见版本如CAMS-GLOBAL-INV-EGG4数据覆盖全球水平分辨率从0.25°到0.5°不等格式通常为netCDF4。用ncdump -h查看文件头能看到类似latitude、longitude、level、time这样的维度以及co2vmr、ch4vmr这类浓度变量。有的版本还会带co2_flux、ch4_flux通量变量负责描述地表与大气之间的碳交换。一定要注意边界条件需要的是大气浓度场不是地表通量。通量单位一般是mol m-2 s-1或kg m-2 s-1浓度单位则多为1e-9、1e-6或具体摩尔分数。两者物理意义差异巨大绝对不能互换。用Python检查时xarray是最趁手的工具import xarray as xr ds xr.open_dataset(cams_egg4_co2_20230101.nc) print(ds.dims) print(ds.data_vars) print(ds[co2vmr].attrs)我每次拿到新数据都会先跑这样三行确认三件事维度名是不是time, lev, lat, lon变量名是否和配置里映射一致单位属性写的是什么。这三件事搞清楚了后面插值和写入才不会翻车。2.2 时间维和垂直维CAMS-Inversion数据的时间频率可能是日均、月均或3小时一次具体要看产品版本。做区域WRF模拟时气象边界条件通常是逐小时温室气体边界也建议至少逐小时更新这样就需要在时间维上做插值。日均数据插值到小时原则上是用前后两天的日平均值做线性过渡虽然无法还原真实的日变化细节但比整个模拟周期只用一个固定浓度要好得多。垂直维度也要格外注意。CAMS的垂直层次一般是模型的大气层从地表到高层可能是等气压面或等η面。WRF模式本身用的是地形追随坐标所以预处理时要把CAMS层次上的浓度插值到WRF格点对应的气压高度。低层大气梯度大边界层内CO2浓度随高度变化明显插值算法建议在气压坐标下做对数线性插值不要直接用线性插值否则近地面浓度容易失真。还有一个常被忽略的坑CAMS模式顶不一定覆盖到WRF模式顶。当CAMS最高层低于WRF顶层时高层缺数据怎么办我一般用最顶层的浓度值进行近等比外推或者干脆用气候态背景浓度填充。对于温室气体这种垂直混合比较充分的高层大气外推带来的误差远小于边界条件完全缺失带来的问题。2.3 单位摩尔分数、质量混合比与通量单位单位换算是函数解析里最容易出错的一环。先列出一张我常用的转换参考表变量类型惯用单位进入WRF边界文件时的建议CO2摩尔分数1e-6即ppm按ppm直接保留CH4摩尔分数1e-9即ppb看WRF-GHG要求若要求ppm则除以1000CO2通量kg m-2 s-1不用于边界条件不要写入CH4通量mol m-2 s-1不用于边界条件不要写入不同CAMS版本的浓度变量单位可能直接写在units属性里也可能写在文档里。有的旧版本CO2变量名是co2实际内容是摩尔分数单位却是1e-6有些版本CH4是ch4vmr单位却是1e-9。最稳妥的做法是在读取函数中自动判断单位属性统一换算成你WRF-GHG配置里要求的单位。我踩过一次特别愚蠢的坑某版本CH4变量实际上是ppb我按ppm读取结果边界浓度放大了1000倍模拟区域内部CH4浓度直接飙到几十ppm一眼假数据。当时为了排查这个错误浪费了整整半天。后来我学乖了在函数解析第一步就强制打印min、max、mean三个统计量看到数值明显超出合理范围就立即停止。3. 核心函数解析Prepy里的读取、插值和写边界到了这篇的关键部分。以下函数名和调用逻辑基于我当时使用的Prepy工具常见实现方式不同版本可能有差异但核心思想是通用的。理解这些函数背后的数据流比死记API更有价值。3.1 load_cams_bc把CAMS文件读进xarray函数的第一步是打开CAMS文件抽取你需要的浓度变量并做最基础的单位与范围检查。我会把它实现成一个独立的读取函数而不是在插值步骤里顺手读目的是保持单一职责。import xarray as xr def load_cams_bc(cams_file, var_nameco2vmr): ds xr.open_dataset(cams_file) if var_name not in ds.data_vars: raise KeyError(f{var_name} 不存在请用 ncdump -h 确认变量名) da ds[var_name] print(shape:, da.shape) print(units:, da.attrs.get(units, N/A)) print(min/max/mean:, float(da.min()), float(da.max()), float(da.mean())) return da这个函数的输出是整个预处理流程里最重要的“体检报告”。如果min或max远超合理范围你根本不应该继续往下跑而是先回去查数据文件本身。我只是把这一步骤写在函数里但实际运行时我建议直接单独跑这个函数确认数据健康再全流程执行。3.2 regrid_cams_to_wrf水平与垂直一致化CAMS是全球规则网格WRF是区域网格两者几何关系完全不同。水平插值我用的是双线性插值通过scipy.interpolate.RegularGridInterpolator实现。实现时先判断WRF网格点是否落在CAMS数据范围内避免出现大面积缺测。import numpy as np from scipy.interpolate import RegularGridInterpolator def cams_horizontal_interp(cams_da, lons_wrf, lats_wrf): # 假设 cams_da 的维度为 (lev, lat, lon) lat cams_da.lat.values lon cams_da.lon.values interp RegularGridInterpolator( (lat, lon), cams_da.values, methodlinear, bounds_errorFalse, fill_valuenp.nan ) pts np.column_stack([lats_wrf.ravel(), lons_wrf.ravel()]) out interp(pts).reshape(lats_wrf.shape) return out垂直插值就要麻烦一点。CAMS的气压层次通常与WRF的η层不一致标准做法是把WRF网格点的气压算出来然后沿气压维进行插值。高层大气的温室气体垂直梯度较小低层梯度大因此我推荐沿气压做对数插值这样可以更好保留近地面浓度变化的细节。代码实现大致是from scipy.interpolate import interp1d def cams_vertical_interp(cams_on_wrf_horiz, pressure_cams, pressure_wrf): # cams_on_wrf_horiz 形状: (lev, n_y, n_x) # pressure_cams: CAMS 各层气压pressure_wrf: WRF 各层气压 n_lev, n_y, n_x cams_on_wrf_horiz.shape out np.full((len(pressure_wrf), n_y, n_x), np.nan) for i in range(n_y): for j in range(n_x): f interp1d( np.log(pressure_cams), cams_on_wrf_horiz[:, i, j], kindlinear, bounds_errorFalse, fill_valueextrapolate ) out[:, i, j] f(np.log(pressure_wrf)) return out这个双重循环在数据量很大时会慢所以在实际项目里我会先裁剪水平区域范围只保留WRF域外扩几十公里范围的CAMS网格点大幅减少插值计算量。你如果跑的模拟域特别大还可以用dask做分块并行但一般情况下不需要这么复杂。3.3 interp_time_and_write_bc时间插值与写入wrfbdyCAMS-Inversion如果是日平均数据而模拟要求逐小时边界就需要时间插值。最朴素的做法是对每个格点沿时间维做线性插值from scipy.interpolate import interp1d def cams_time_interp(cams_array, time_in, time_out): f interp1d( time_in, cams_array, axis0, kindlinear, bounds_errorFalse, fill_valueextrapolate ) return f(time_out)真正写入wrfbdy_d01文件时不同WRF-GHG版本对变量名有要求常见的是co2_b、co2_p这类前缀组合。_b代表基础值_p代表扰动值。Prepy函数一般会把插值后的边界浓度拆成这两部分按WRF时间维顺序依次写入对应格点。这里有一个特别容易踩的坑WRF边界文件里的时间维长度通常是模拟时次数加一因为边界需要在起始时刻前多一个侧边值供差分使用。写之前一定要确认数组形状完全匹配否则netCDF4写入时会报维度不匹配或者更糟糕地静默截断。边界条件并不需要把整个三维场全部写入wrfbdy只需要提取南北东西四个侧边以及四个角点附近的一圈格点。Prepy的写函数一般会根据边界标记位从全场中抽取这些点这样生成的文件体积非常小读取速度也快。函数调用顺序整体上是读文件→水平插值→垂直插值→时间插值→抽取边界带→写入wrfbdy。链路上每一步的中间结果都可以保存成日志里的统计量方便排查。3.4 函数调用顺序与整体流程把上面几个函数串起来整体流程大致是这样读配置得到CAMS文件路径、模拟时间段、WRF网格模板文件路径对每个CAMS文件调用load_cams_bc确认数据健康读取WRF网格的经纬度和气压场调用cams_horizontal_interp调用cams_vertical_interp把浓度插值到WRF的η层对模拟时间序列调用cams_time_interp生成逐小时三维浓度抽取边界带按WRF边界文件格式写入wrfbdy_d01Prepy工具通常会在配置文件中暴露interp_method、variable_map、time_frequency这些选项让你不需要改代码就能适配新数据。理解每个函数的职责后遇到问题你就能很快定位是读取阶段、插值阶段还是写入阶段出了错而不是在几千行日志里盲目找。4. 实操过程从CAMS-Inversion下载到边界条件生成光讲原理容易飘上一段踩在实地上。我以一个为期10天的模拟实验为例演示从下载数据到生成边界条件的完整过程。声明一下以下是基于我自己项目目录结构和当时工具使用习惯整理的一套流程你拿到Prepy仓库后可以对照调整。4.1 下载与文件整理去CAMS官方数据门户检索CAMS global inversion产品选好变量CO2、CH4和起止时间。下载时建议选择netCDF格式方便xarray处理。文件下载回来后统一放在一个目录下project/ ├── data/ │ └── cams_inversion/ │ ├── cams_egg4_co2_20230101.nc │ ├── cams_egg4_co2_20230102.nc │ └── ch4/ │ └── cams_egg4_ch4_20230101.nc ├── wrf/ │ ├── wrfinput_d01 │ └── wrfbdy_d01 └── scripts/ └── prepy_bc_config.yaml把数据文件命名统一成cams_egg4_species_YYYYMMDD.nc后面函数解析模式匹配会省很多事。不同物种分开目录存储也是避免混淆的一种好习惯。4.2 配置Prepy参数配置文件我习惯写成YAML核心字段是边界数据路径、物种映射、时间范围和WRF网格文件路径。类似于下面这样bc: dataset: cams_inversion source_dir: /home/user/project/data/cams_inversion start_time: 2023-01-01 00:00:00 end_time: 2023-01-10 00:00:00 species: [co2, ch4] variable_map: co2: co2vmr ch4: ch4vmr interp_method: linear wrf_file: /home/user/project/wrf/wrfinput_d01 output_bdy: /home/user/project/wrf/wrfbdy_d01variable_map的作用我之前提过就是把CAMS文件中的变量名与Prepy内部使用的通用物种名对应起来。CAMS版本一旦变化通常只需要改这里而不需要改代码。4.3 运行并验证配置好后执行Prepy的主脚本。不同仓库入口不同一般用法类似python wrf_ghg_prepy.py --mode bc --config prepy_bc_config.yaml运行过程中日志里如果出现Loading CAMS file ... OK、Horizontal interpolation ... OK、Vertical interpolation ... OK这类信息说明链路基本通畅。跑完以后不要急着直接提交模式运行任务先做验证。最基本的验证是用ncdump看生成的wrfbdy文件ncdump -h wrfbdy_d01 | grep -i co2重点看两件事变量名是否出现co2_b和co2_p时间维长度是否符合预期。再用ncks抽取边界格点数据看一眼数值ncks -v co2_p -d south_north,0,2 -d west_east,0,2 wrfbdy_d01 | head -40正常情况下CO2边界浓度应该在400~430ppm的范围内波动CH4应该在1800~2000ppb左右。如果数值差了数量级十有八九是单位换算环节出错回看2.3节。5. 常见问题排查与避坑实录这部分是实打实的排障记录。CAMS-Inversion数据预处理说难不难但坑确实不少我把遇到过的典型问题和解决办法整理如下。5.1 变量名找不到或维度名混乱CAMS旧版和新版变量的命名风格不一致。你调用load_cams_bc时如果报KeyError不要怀疑代码先回去看文件头。用xarray打印ds.data_vars看看实际叫co2、co2vmr还是co2_conc。维度名也是一样有的版本叫level有的叫lev有的叫generic。遇到这种情况最简单的处理是在读取函数里加一个重命名操作da da.rename({level: lev})我建议把维度名统一成time, lev, lat, lon四个后续插值函数全按这个约定写这样可以避免在每个函数里重复判断。5.2 插值后出现NaN或负数NaN的来源通常有三种。第一CAMS文件本身部分区域有掩膜比如海洋某些变量的低层格点无数据第二WRF网格范围超出了CAMS数据经纬度边界双线性插值时返回了NaN第三垂直插值时WRF气压超出了CAMS气压范围。排查方法是用numpy.isnan统计NaN比例再画个空间分布图看哪些区域有问题。如果只是超出边界把bounds_errorFalse改成fill_valueNone可以让插值用最近边缘值填充或者把WRF网格裁剪到CAMS覆盖范围内。如果CAMS文件本身有缺失建议直接用周边格点的平均或时间邻近日填充必要时做一次空间平滑。负数浓度问题主要出在时间插值或垂直插值在边界处的过冲。CO2浓度出现负值没有物理意义对CH4也一样。稳妥的处理是在插值输出后进行截断np.clip(out, 0, None, outout)这种方法虽然粗糙但至少能保证写入边界文件的数值不会让模式崩溃。更高阶的做法是改用一次单调插值或者先对数变换再插值再还原。5.3 边界条件与内部通量不一致导致的怪异结果用CAMS-Inversion做边界时一个绕不开的问题是外部反演通量和内部VPRM/排放清单之间存在时空尺度差异。CAMS日平均浓度变化平滑而VPRM响应太阳辐射、温度和湿度有明显日变化两者在边界附近容易“打架”。你可能会看到模拟区域边缘出现带状浓度异常或者浓度场随风向出现锯齿。缓解方法有几个。一是提高边界更新频率CAMS虽然是日均但你可以用前后两天的日均值做小时插值避免边界上出现阶梯跳跃。二是在模拟区域边缘设置缓冲区比如把研究关注区放在区域内部边界外推区内不进行分析。三是在做结果分析时弃掉边界附近若干网格点的数值只分析内部区域。我自己实际跑下来的一个教训是初始条件和边界条件最好来自同一个数据源。一次模拟中边界浓度用的是CAMS平均值大约414ppm而初始场来自另一套全球再分析数据区域中心初始值只有404ppm。模式起跑后前6小时整个区域浓度都朝边界值调整产生了一个明显的“假趋势”。后来我把初始场和边界场统一用CAMS问题迎刃而解。5.4 针对Prepy工具的具体排查手段如果Prepy脚本本身报错我的排查顺序是先看WPS/real生成的wrfbdy模板是否存在、变量命名是否完整再检查配置文件中variable_map是否和数据文件一致接着用一个小脚本单独测试读取和插值两步跳过量大的写文件环节最后才跑完整流程。调试时尽量在配置里开启debug日志。Prepy这类工具往往会把数据文件名、插值范围、数组shape打印出来这些信息比任何报错信息都值钱。如果工具没有debug开关就在关键函数里手动加几行print输出数组min/max和shape跑完再删掉。我特别推荐一个小技巧每次处理新的CAMS版本数据先用小脚本把第一个文件读进来打印dimensions、variables、units、min/max/mean确认全部合理后再跑全流程。这一步看起来很基础但真的能省掉后面十倍的时间。预处理工具出问题多数时候不是算法难而是数据描述文件本身的信息没吃透。把数据脾性摸清楚后Prepy里那些解析函数天然就会老实听话。