写在前面

最近在计算总湿稳定度，涉及到散度项、垂直积分等计算。这里简单记录一下。

1、GMS-总湿稳定度定义

为了描述对流和相关环流排放大气柱中水汽的效率，Neelin和Held [1987]提出了“总湿稳定性（GMS）”的最初概念，其定义为每单位垂直质量通量的净垂直MSE输出。

对于波动的不稳定发展，根据水汽模态理论，波动对流引起的环流必须进一步湿润大气，而这需要GMS的负值，或者如果与外部强迫的影响（包括表面热通量和辐射加热）相结合，则需要有效GMS的负值，以维持不稳定模态。

前人的一些研究中发现：GMS的垂直分量占总GMS的主要部分，并且与赤道波动相速度表现出更强的相关性（相关系数大约 0.9）。一些猜测认为：
行星尺度波动的对流区域静态稳定性降低使得斜压模式存在。斜压垂直结构允许对流层低层和高层之间的能量传递，导致较低的GMS。由于GMS较低，有效等效深度较低将反过来降低波动的相速度。

GMS最初的定义为：

$${\Gamma }_T=-\frac{T_R[\mathbf{v}\cdot \nabla s+\omega (\partial s/\partial p)]}{L[\nabla \cdot (r\mathbf{v})]},$$

并且进一步可以分解为水平项和垂直项：

$${\Gamma }_H=-\frac{T_R[\mathbf{v}\cdot \nabla s]}{L[\nabla \cdot (r\mathbf{v})]}$$

$$

\\Gamma_{V}=-\frac{T_{R}[\omega(\partial s/\partial p)]}{L[\mathbf{\nabla}\cdot(r\mathbf{v})]}
$$

好的扯远了，在计算GMS过程中就涉及到了垂直积分,这里主要讲一下相关注意点。之前老有人讨论计算垂直积分的问题：

• 在使用积分函数时，需要注意xarray.integrate()与numpy.trapz()的结果是等同的，但是xarray默认是根据level来自动运行的，如果level是hPa，那么积分后的数据仍然是hPa。

计算垂直积分需要注意两点：

• 1、先确保气压是的顺序，举个例子，我这边一般习惯确保level的排列是从100hPa-1000hPa；
• 2、转换气压层单位为pa，为国际基本单位，方便一些量级的统一。

为什么积分的顺序很重要呢，因为这涉及到后面你要不要乘上负号。因为在p坐标系下，利用静力平衡需要：

$$dp=-\rho gdz,\rho dz=-\frac{dp}{g}$$

这里以梯形积分计算整层水汽积分为例,如果你的level的顺序是从100hPa到1000hPa的话，很简单计算就直接：

$$\text{Column-integrated specific humidity}=\frac{1}{g}{\int }_{p_t}^{p_s}q(p)dp$$

• q 是比湿
• $p_s$是地表气压
• $p_t$是顶层气压

以下是xarray和numpy两个计算示例，得到的结果是相同的，再次强调下面数据中level是从100hPa-1000hPa排列的：

基于numpy函数计算

def numpy_intergrated(sph_interp):
    pressure_levels = sph_interp.level * 100  # 将 hPa 转换为 Pa
    # 重力加速度
    g = 9.81  # m/s²
    # 使用梯形积分法计算整层水汽积分
    sph_trap = np.trapz(sph_interp, pressure_levels, axis=1) / g
    print(sph_trap.min(),sph_trap.max())
    return sph_trap
sph_trap_numpy = numpy_intergrated(sph_interp) 

3.44126774828972 62.35031826524259

基于xarray函数计算

def xarray_integrate(sph):
    g = 9.81  # m/s²
    # 计算垂直积分 (Column-integrated specific humidity)
    column_integrated_sph = (sph).integrate(coord='level') / g *100
    print(column_integrated_sph.min(), column_integrated_sph.max())
    return column_integrated_sph
sph_integrate_xarray = xarray_integrate(sph_interp)    


<xarray.DataArray 'level' (level: 27)> Size: 216B
array([ 10000.,  12500.,  15000.,  17500.,  20000.,  22500.,  25000.,
        30000.,  35000.,  40000.,  45000.,  50000.,  55000.,  60000.,
        65000.,  70000.,  75000.,  77500.,  80000.,  82500.,  85000.,
        87500.,  90000.,  92500.,  95000.,  97500., 100000.])
Coordinates:
  * level    (level) float64 216B 100.0 125.0 150.0 175.0 ... 950.0 975.0 1e+03
<xarray.DataArray 'q' ()> Size: 8B
array(3.44126775) <xarray.DataArray 'q' ()> Size: 8B
array(62.35031827)

如果高度层是从1000-100hpa的降序排列的话，可以在进行积分时进行取反计算:

svdv_int  = np.trapz(s_vdv[:,::-1,:,:], plev[::-1], axis=1)*units('m**2/s**3')*units('Pa') / g

上图中，左图是反转高度的结果（即100-1000hPa），右图是不反转的结果（即1000-100hPa，明显的发现只是符号的正负差异。

或者对于数据的level进行升序排列，这可以通过xarray.sortby('level')实现

另外一个小技巧是，在进行计算时带上单位，保证量级范围是正确的。积分的过程是需要乘上高度的单位的，一般是Pa或者hPa

Benedict, J. J., E. D. Maloney, A. H. Sobel, and D. M. W. Frierson, 2014: Gross moist stability and MJO simulation skill in three full-physics GCMs. Journal of the Atmospheric Sciences, 71, 3327–3349, https://doi.org/10.1175/JAS-D-13-0240.1.

其他

$$\nabla \cdot (s\mathbf{v})=\frac{\partial }{\partial x}(su)+\frac{\partial }{\partial y}(sv)$$

进一步的分解为：

$$\mathbf{v}\cdot \nabla \mathbf{s}=u\frac{\partial \mathbf{s}}{\partial x}+v\frac{\partial \mathbf{s}}{\partial y}$$

$$s(\nabla \cdot \mathbf{v})=s(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y})$$