FLEXPART 模型笔记

jokervTv 2020/11/30

写下这篇文章时，FLEXPART 停留在 version 8.2

拉格朗日模型与欧拉模型相比，一个显著的优势是没有数值扩散。此外，在欧拉模型中，从点源释放的示踪剂会立即在网格中混合，而拉格朗日模型与计算网格无关，并且理论上具有无限小的分辨率。

TO hard to me

边界层参数的物理参数化
粒子传输与扩散
羽流轨迹
Removal processes

# 边界层参数的物理参数化

如果 surface sensible heat fluxes 和 surface stresses 可以直接给出

total surface stress 由下式计算

$\tau = \sqrt{\tau_1^2 + \tau_2^2}$

其中 $\tau_1$ 和 $\tau_2$ 是东西和南北方向的 surface stress

摩擦速度是将在子过程 scalev.f 中计算

$u_* = \sqrt{\frac{\tau}{\rho}}$

其中 $\rho$ 为空气密度。

如果 surface sensible heat fluxes 和 surface stresses 不可以直接给出，则通过将 Berkowicz 和 Prahm(1982)(pbl_profile.f) 的方法应用于第二个模型高度： 10m（风）和 2m（温度）的风和温度数据 (以前是第一个模型)；由于 ECMWF 的第一个模型级别现在接近 10m，因此改用第二个级别：

$\displaystyle u_{*}=\frac{\kappa \Delta u}{\ln \frac{z_{l}}{10}-\Psi_{m}\left(\frac{z_{l}}{L}\right)+\Psi_{m}\left(\frac{10}{L}\right)} \\ \displaystyle \Theta_{*}=\frac{\kappa \Delta \Theta}{0.74\left[\ln \frac{z_{l}}{2}-\Psi_{h}\left(\frac{z_{l}}{L}\right)+\Psi_{h}\left(\frac{2}{L}\right)\right]} \\ \displaystyle L=\frac{\bar{T} u_{*}^{2}}{g \kappa \Theta_{*}}$

其中 $\kappa$ 是 von Karman 常数 (0.4)， $z_l$ 是第二个模型水平的高度， $\Delta u$ 是第二个模型高度与 10m 高度的水平风速之差， $\Delta \Theta$ 是在第二个模型水平的风速之差在 2m 模型水平，在 2m 处， $\Psi_{m}$ 和 $\Psi_{h}$ 是动量和热的稳定性校正函数 (Businger 等人，1971； Beljaars 和 Holtslag，1991)， $g$ 是重力引力， $\Theta_*$ 是温度标度， $T$ 是平均表层温度（在第一个模型级别取为 $T$ )。

然后通过以下公式计算热通量

$\left(\overline{w^{\prime} \Theta_{v}^{\prime}}\right)_{0}=-\rho c_{p} u_{*} \Theta_{*}$

其中 $\rho c_{p}$ 是空气的定压比热容。

ABL 高度是根据 Vogelezang 和 Holtslag(1996) 使用临界理查森数概念 (richardson.f)。

$R i_{l}=\frac{\left(g / \Theta_{v 1}\right)\left(\Theta_{v l}-\Theta_{v 1}\right)\left(z_{l}-z_{1}\right)}{\left(u_{l}-u_{1}\right)^{2}+\left(v_{l}-v_{1}\right)^{2}+100 u_{*}^{2}}$

ABL 高度 $h_{mix}$ 设置为 Richardson 数超过临界值 0.25 的首个模型级别 $l$ 的高度。 $\Theta_{v1}$ 和 $\Theta_{vl}$ 是虚位温， $z_1$ 和 $z_l$ 是其高度； $(u_1, v_1)$ 和 $(u_l, v_l)$ 是第 1 和第 l 个模型高度的风分量。通过将 $\Theta_{v1}$ 替换为下式，可以改进对流情况：

$\Theta_{v 1}^{\prime}=\Theta_{v 1}+8.5 \frac{\left(\overline{w^{\prime} \Theta_{v}^{\prime}}\right)_{0}}{w_{*} c_{p}}$

其中

$w_{*}=\left[\frac{\left(\overline{w^{\prime} \Theta_{v}^{\prime}}\right)_{0} g h_{m i x}}{\Theta_{v 1} c_{p}}\right]^{1 / 3}$

是对流速度标度。等式右侧的第二项表示过多的上升热量。由于 $w_*$ 事先未知，因此 $h_{mix}$ 和 $w _*$ 是互相迭代计算的。

::: note 标度还是尺度，我好混乱～～ :::

在 ECMWF 模型无法有效求解的尺度上，ABL 高度的时空变化在确定有效混合示踪剂层的厚度方面起着重要作用。由于复杂地形，土地利用或土壤湿度的变化，ABL 高度的空间变化也存在类似的观点（Hubbe 等，1997）。在一个复杂地形的表面上移动的示踪剂云的厚度将由最大值而不是平均 ABL 高度确定。

在 FLEXPART 中，使用某种程度的参数设置来避免在示踪剂云厚度和表面示踪剂浓度方面产生明显偏差。为了解决地形引起的空间变化，我们使用 "envelope" ABL 高度

$H_{e n v}=h_{m i x}+\min \left[\sigma_{Z}, c \frac{V}{N}\right]$

此处， $\sigma_Z$ 为 ECMWF 模型子网格地形的标准偏差，c 为常数（此处为 2.0），V 为 $h_{mix}$ 高度的风速，N 为 Brunt-Vaisala 频率。 Under convective conditions, the envelope ABL height is, thus, the diagnosed ABL height plus the subgrid topography (assuming that the ABL height over the hill tops effectively determines the dilution of a tracer cloud located in a convective ABL). Under stable conditions, air tends to flow around topographic obstacles rather than above it, but some lifting is possible due to the available kinetic energy. $\frac{V}{N}$ is the local Froude number (i.e., the ratio of inertial to buoyant forces) times the length scale of the sub-grid topographic obstacle. The factor $c\frac{V}{N}$ , thus, limits the effect of the subgrid topography under stable conditions, with c = 2 being a subjective scaling factor. $H_{env}$ rather than $h_{mix}$ is used for all subsequent calculations. In addition, $H_{env}$ is not interpolated to the particle position, but the maximum $H_{env}$ of the grid points surrounding a particle’s position in space and time is used.

# 粒子传输与扩散

# 粒子轨迹计算

FLEXPART 通常使用简单的“零加速度”方案（精确到一阶）来积分轨迹方程（Stohl，1998 年）

$X(t + \Delta t) = X(t) + v(X,t) \Delta t$

$t$ 是时间， $\Delta t$ 是时间增量， $X$ 是位置向量，并且 $v = \bar{v} + v_t + v_m$ 由网格尺度合成的风矢量 v，湍流风波动 $v_t$ 和中尺度风波动 $v_m$ 组成的风矢量。

从 FLEXPART 5.0 版开始，只要有可能，都对 Petterssen（1940）方案进行一次迭代（精确到二阶），从而提高了数值精度，但仅适用于网格尺度风。通过“零加速度”方案获得的位置进行的校正。在三种情况下，它无法应用。第一种情况，Petterssen 方案需要二次风场，该风可能在内存中保留的两个风场的时间间隔之外。第二种情况，粒子越过嵌套域的边界。第三种情况，则如果 ctl> 0（请参见下文），并在 ABL 中。

粒子传输和湍流扩散由子例程 advance.f 处理，在该子例程中调用了将风和其他数据插值到粒子位置的过程的调用，并求解了 Langevin 方程（请参见下文）。极点是纬度/经度网格上的奇点。因此，纬度（switchnorth，switchsouth）的极向的水平风（变量 uu，vv）被转换为极地立体投影（变量 uupol，vvpol），在该投影上计算了粒子对流。由于 uupol，vvpol 也存储在纬度/经度网格上，因此无需进行其他插值。

# Langevin 方程

TODO : 好像很难得样子

# 确定时间步长

FLEXPART 可以在两种不同的模式下使用。计算速度较快的一个（文件 COMMAND 中的 ctl<0）不能使计算时间步长适应拉格朗日时间尺度 $\tau_{L_i}$ （其中 i 是三个风分量之一），并且 FLEXPART 使用一个同步时间间隔（lsynctime，已指定）的恒定时间步长在 COMMAND 文件中，通常为 900 秒）。通常，在此模式下自相关非常低，并且湍流描述得不好。然而，对于大规模应用，FLEXPART 可以很好地使用此选项（Stohl 等，1998）。如果要更准确地描述湍流，则时间步长必须受 $\tau_L$ 限制。由于垂直风最重要，因此仅使用 $\tau_{L_w}$ 。用户必须在文件 COMMAND 中指定两个常量 ctl 和 ifine。第一个根据以下步骤确定时间步长 $\Delta t_i$

$\Delta t_{i}=\frac{1}{c_{t l}} \min \left(\tau_{L_{w}}, \frac{h}{2 w}, \frac{0.5}{\partial \sigma_{w} / \partial z}\right)$

$\Delta t_{i}$ 的最小值为 1 秒。 $\Delta t_{i}$ 用于求解水平湍流风分量的 Langevin 方程。为了求解垂直风分量的 Langevin 方程，使用了较短的时间步长 $\Delta t_{w} = \dfrac{\Delta t_{i}}{ifine}$ 。但是，请注意，在小于 $\Delta t_{i}$ 的时标上，水平和垂直风分量之间没有交互作用。这种策略（为 ctl 和 ifine 提供足够大的值）可确保粒子在非常不均匀的湍流中也保持垂直均匀混合，同时将计算成本保持在最低水平。

# 风波动的参数化

Hanna（1982）针对σvi 和τLi 提出了一种基于边界层参数 h，L，w ∗，z0 和 u ∗的参数化方案，即 ABL 高度，Monin-Obukhov 长度，对流速度标度，粗糙度长度和摩擦速度。在子程序 hanna.f，hanna1.f，hanna_short.f 中使用了 Ryall 等人的修改。（1997）对于σw，因为 Hanna 的方案并不总是在整个对流 ABL 中产生σw 的平滑剖面。在下文中，下标 u 和 v 分别表示顺风分量和侧风分量（在子程序 windalign.f 中转换为网格坐标），而 w 表示湍流速度的垂直分量； f 是科里奥利参数。所使用的最小τLu，τLv 和τLw 分别为 10 s，10 s 和 30 s，以避免对靠近表面的粒子的过多计算时间

# 中尺度风波动

中尺度运动既不能由 ECMWF 数据解析，也不能由湍流参数化覆盖。由于中尺度运动可以显着加速分散羽流的生长，因此至少必须以近似方式考虑该未解决的光谱间隔（Gupta 等，1997）。为此，我们使用与 Maryon（1998）类似的方法，即针对中尺度风速波动（以 Maryon 的术语“ meandering”）求解独立的 Langevin 方程。假设风的方差在栅格尺度上提供了一些关于其子栅格方差的信息，则将用于中尺度 Langevin 方程的风速标准偏差设置为 turbmesoscale（在文件 includepar 中设置）乘以围绕栅格的栅格点的标准偏差。粒子的位置。假设网格点之间的线性插值可以恢复一半的子网格变异性，则相应的时间尺度取为可获得风场的间隔的一半，这并非不可能（Stohl 等，1995）。这种经验方法没有描述实际的中尺度现象，但是它类似于用于评估轨迹精度的集成方法（Kahl，1996； Baumann 和 Stohl，1997； Stohl，1998）。

一个重要的传输机制是对流云的上升气流。它们与云层中的下降气流并补偿无云环境中的沉降共同发生。这些对流输运在垂直方向上是网格尺度，但在水平方向上是子网格尺度，并且不能用 ECMWF 垂直速度表示。为了表示颗粒扩散模型中的对流传输，必须在整个垂直列中重新分布颗粒。对于 FLEXPART，我们选择了 Emanuel 和 Zivkovi ˇ Rothman（1999）的对流参数化方案，因为它依赖于网格规模的温度和湿度场并计算位移矩阵，从而为粒子重新分布提供了必要的质量通量信息。使用文件 COMMAND 中的 lconvection 可以打开对流参数化。它的计算时间与垂直模型级别数的平方成比例，并且使用当前的 60 级 ECMWF 数据可能最多占 FLEXPART 的计算时间的 70％。对流是在子例程 convmix.f，calcmatrix.f，convect43c.f 和 redist.f 中计算的。它被称为每个 FLEXPART lsynctime 时间步长（通常为 900 s），其中包含时间插值的温度和 ECMWF 数据中的特定湿度曲线。请注意，对流方案中使用的是原始 ECMWF 模型级别，而不是笛卡尔坐标。出于效率考虑，在对流方案中每个网格列调用一次之前，根据粒子的水平网格位置（sort2.f）对其进行排序。

# 粒子分裂

在从大气中的点源扩散的初始阶段，粒子通常会形成致密的云团。相对较少的粒子足以正确模拟此初始阶段。但是，一段时间后，粒子云失真，粒子散布到更大的区域。现在需要更多的粒子。FLEXPART 允许用户指定时间常数Δts（命令 COMMAND）。在Δts，2Δts，4Δts，8Δts 等的传播时间（子例程 timemanager.f）之后，粒子被分为两部分（每个粒子都接受原始粒子质量的一半）。

# 羽流轨迹

在最近的一篇论文中，Stohl 等人。（2002 年）提出了一种使用聚类分析来浓缩复杂的大型 FLEXPART 输出的方法（Dorling 等，1992）。其背后的思想是在每个输出时间对源自释放点的所有粒子的位置进行聚类，并仅写出聚类的粒子位置以及其他信息（例如，ABL 和平流层中的粒子比例）。这样创建的信息几乎与传统轨迹一样紧凑，但会引起湍流和对流。可以通过在文件 COMMAND 中将 iout 设置为 4 或 5 来激活此选项。可以使用 includepar 文件中的参数 ncluster 设置群集数。由子例程 plumetraj.f 处理聚类并产生输出。

jokervTv 爬坑历程

Choose mode