四维变分同化在非常规卫星观测资料的直接同化方面具有明显优势, 所以在最近十几年中, 越来越多受到关注。中国气象科学研究院从2005年开始针对完全可压的非静力GRAPES (Global/Regional Assimilation and PrEdiction System) 模式开发可投入业务使用的四维变分同化系统 (GRAPES 4DVAR)。目前, 国际上只有少数几个类似的四维变分同化系统。GRAPES是中国气象局新一代全球/区域多尺度通用同化与数值预报系统^[1-4]。GRAPES 4DVAR以GRAPES三维变分同化系统^[5-10]为基础, 通过GRAPES切线性模式及其伴随模式的积分将不同时刻的观测与初始时刻分析场联系在一起, 从而实现四维变分同化^[11]。切线性模式是非线性模式的一阶线性近似, 它被用来描述扰动态在基本态附近的演变。切线性模式同时计算基本态和扰动态随时间演变, 其中基本态的计算和非线性模式完全一致, 扰动态由切线性方案利用基本态计算得到。

2005年底, 中国气象局数值预报创新基地与国防大学计算学院合作完成了只包括动力框架部分的GRAPES切线性模式的开发, 随后开始研究物理过程方案的线性化问题。物理过程的基本特点是强非线性和不连续性, 一阶线性近似误差可能比较明显, 所以在设计物理过程的切线性方案时至少会面临两个问题:①强非线性和不连续性对一阶线性近似有什么样的影响? ②当一阶线性近似的误差超过可以接受的范围时能否找到有效的应对方法?

本文以GRAPES模式垂直扩散方案线性化为切入点, 研究简单物理过程线性化问题, 目的是将切线性垂直扩散方案合理地接入GRAPES切线性模式, 为GRAPES 4DVAR系统的改进打下基础。

GRAPES模式采用MRF非局地边界层方案计算边界层和自由大气中的湍流扩散倾向。本文只研究自由大气垂直扩散方案的线性化问题, 下一步工作将在此基础上研究边界层部分的内容。

湍流扩散方程为, 其中C代表模式物理量风速u, v, 位温θ和比湿q, K是湍流扩散系数。

在MRF方案中, 自由大气中的垂直扩散系数用下面的公式计算:K=l ² f , 其中l是摩擦长度; f是稳定度函数, 它根据理查孙数R_i计算得到; 是垂直风切变。

垂直扩散的下边界条件由地表动量通量、地表热量通量和地表水汽通量提供, 它们分别是:

(1)

线性垂直扩散方案的改进涉及到地表动量通量, 计算方案为

(2)

式 (2) 中, u^*是摩擦速度, 是全风速。

有关MRF非局地边界层方案的具体细节参见文献[12]。

在GRAPES模式垂直扩散方案源代码的基础上进行逐句线性化可以得到切线性垂直扩散方案。为了检验切线性垂直扩散方案的正确性, 对增加切线性垂直扩散方案前后的GRAPES切线性模式进行了批量的正确性测试^①。表 1和表 2给出了2005年8月27日的个例测试结果。增加切线性垂直扩散方案前后, GRAPES切线性模式积分12 h的正确性测试精度都能够达到6位的精度。但是, 切线性模式能够通过正确性测试只能说明切线性模式在原则上没有问题, 当扰动态的量级足够小时它是非线性模式的一阶线性近似。在实际四维变分同化试验中, 扰动态 (即同化分析增量) 的量级显然是比较大的, 有时甚至与背景场量级相当, 此时一阶线性近似误差是否可以接受还是一个问题。

①张林, 朱宗申.GRAPES四维变分同化系统中线性大尺度凝结方案的开发与试验.待发表.

表 1

表 1 只包括动力框架的切线性模式的正确性测试结果 Table 1 Test results of tangent-linear model without physics

表 1 只包括动力框架的切线性模式的正确性测试结果 Table 1 Test results of tangent-linear model without physics

表 2

表 2 包括动力框架和垂直扩散方案的切线性模式的正确性测试结果 Table 2 Test results of tangent-linear model with vertical diffusion

表 2 包括动力框架和垂直扩散方案的切线性模式的正确性测试结果 Table 2 Test results of tangent-linear model with vertical diffusion

切线性模式描述扰动态在基本态附近的演变, 所以可以根据切线性模式中的扰动和非线性模式中扰动的偏差来评估当扰动态的量级比较大的时候切线性模式的计算精度。如果用X代表非线性模式的初始场, 用δX代表初始场的扰动, 分别用X和X+δX积分非线性模式得到两次积分结果M_N(X) 和M_N(X+δX), 它们之间的差异M_N(X+δX)-M_N(X) 就是不同时刻非线性模式中的扰动, 称之为非线性扰动。用X作为基本态的初始场, 用δX作为扰动态的初始场, 积分切线性模式就可以得到不同时刻切线性模式中的扰动M_T(δX), 称之为切线性扰动。用切线性扰动和非线性扰动的平均偏差s作为指标来评估切线性模式一阶线性近似的误差, s的表达式如下:

(3)

式 (3) 中, M_N表示非线性模式的积分, M_T表示切线性模式的积分。

以2005年8月7—27日00:00(世界时, 下同) 国家气象中心业务中期预报系统T213的分析场作为初始场进行了21个个例的批量试验。GRAPES非线性模式和切线性模式在水平方向上都采用等经纬度网格, 水平格距是1.0°, 格点数是37×31, 在垂直方向上分为15层。非线性模式需要积分两次, 两次积分的不同在于初始场之间相差一个扰动场。选择实际探空资料四维变分同化试验的分析增量场作为非线性模式初始场的扰动, 同时它们也是切线性模式扰动态的初值。变分分析增量场的选择不影响本文的结论, 所以本文没有介绍同化试验的具体方案。为了比较切线性扰动和非线性扰动的差异, 将非线性模式和切线性模式同样积分12 h。

本文的研究对象是垂直扩散方案, 所以在积分非线性GRAPES模式时只包括动力框架和垂直扩散方案, 将这样计算得到的非线性扰动作为衡量标准, 如果切线性扰动与它之间的偏差越小就认为切线性模式计算结果的质量越高。图 1给出了切线性模式有无垂直扩散方案时切线性扰动和非线性扰动的平均偏差, 其中无量纲气压、风场u、风场v、比湿q都是GRAPES模式的预报变量, 前3个变量的平均偏差是21个个例所有垂直层次平均, 比湿q则是21个个例最低5层的平均。当切线性模式只包括动力框架时, 切线性扰动和非线性扰动的偏差是较大的, 增加切线性垂直扩散方案后切线性模式计算结果的质量有明显的改进, 尤其是比湿场, 但是在个别时刻 (例如07:00和09:00), 切线性扰动气压场和风场与非线性扰动的偏差甚至超过没有切线性垂直扩散方案时的结果 (图 1中的虚线)。在21个个例中, 只有8月9日和8月27日的计算结果会出现异常, 而且异常时刻也只局限在07:00和09:00, 所以认为直接在非线性垂直扩散方案源代码基础上线性化, 不做任何修改而得到的切线性垂直扩散方案在少数情况下是不合适的。

图 1

图 1. 切线性模式与只有动力框架和垂直扩散方案非线性模式平均偏差 (虚线:切线性模式只有动力框架; 实线:切线性模式包括动力框架和垂直扩散方案) (a) 无量纲气压, (b) u风场, (c) v风场, (d) 比湿q Fig 1. Mean departure between the tangent-linear model and the nonlinear model including vertical diffusion (dashed line denotes the results are obtained with a dry tangent-linear model; solid line denotes a tangent-linear model including vertical diffusion) (a) non-dimensional pressure, (b) u-wind, (c) v-wind, (d) specific humidity

增加切线性垂直扩散方案后切线性扰动与非线性扰动的偏差反而增大, 这说明垂直扩散方案一阶线性近似的误差已经超出可以接受的范围。在这种情况下, 首先, 要找到造成一阶线性近似误差增大的原因, 在此基础上对切线性方案进行修正。本文认为非线性方案是实际大气较正确的描述, 所以应该尽量不修改非线性垂直扩散方案。对切线性垂直扩散方案进行修正的原则是:①不能明显降低切线性模式计算结果的质量; ②修正后的切线性方案存在相应的伴随方案, 不然就无法在四维变分同化系统中应用。

引起一阶线性近似误差增大的原因可能来自于垂直扩散方案的不连续性或者强非线性。物理过程不连续性是物理过程中的开关造成的。在GRAPES模式垂直扩散方案中只存在一个开关, 是对层结不稳定性的判断, 它决定理查孙数R_i用不同的公式进行计算。为了研究开关的影响, 对GRAPES模式的非线性垂直扩散方案进行了简化, 不对层结不稳定性作判断, 用同样的公式计算理查孙数R_i, 这样就避免了开关的出现。但是, 开关是否出现对GRAPES切线性模式的计算结果没有明显影响?切线性垂直扩散方案没有开关时, 切线性模式的扰动气压场和扰动风场也同样存在异常, 这说明垂直扩散方案的不连续性不是造成一阶线性近似误差增大的根本原因。

Mahfouf^[13]指出, 在设计切线性垂直扩散方案时忽略垂直扩散系数的扰动可以避免误差的异常增长, 但是根据试验结果, 这一经验在GRAPES切线性垂直扩散方案中并不适用, 忽略垂直扩散系数的扰动没有避免异常出现, 反而降低了切线性模式的计算精度。仔细分析切线性模式计算结果后发现, 异常首先出现在最低层的扰动风场上。在MRF边界层方案的数值求解中, 地表动量通量-和-提供下边界条件, 它们的线性近似误差在极少数格点上非常大, 从而导致了最低层扰动风场的异常, 随后周围气压场由于散度场的异常而调整。因为只在极少数格点上出现这样的情况, 随着模式的积分调整, 扰动气压场和扰动风场的异常只能维持很短的时间。图 2给出了异常格点上切线性模式最低层扰动变量的演变情况, 它来自于2005年8月27日的个例。08:00, 地表动量通量-的扰动开始发生异常, 引起最低层扰动u风场的异常。09:00, 最低层扰动u风场甚至达到-161.9 m/s。最低层扰动v风场的异常情况与u风场的类似。如果在切线性垂直扩散方案中忽略可能出问题的地表动量通量扰动, 也就是不计算扰动地表动量通量-和-, 那么异常不再出现, 最低层扰动风场的变化与非线性模式中的情况在量级和趋势上是一致的 (图 3)。图 4给出了修改切线性垂直扩散方案前后切线性扰动和非线性扰动的平均偏差, 其中虚线代表修改切线性垂直扩散方案前的结果, 与图 1中的实线是一样的。修改切线性垂直扩散方案后, 切线性模式的计算结果不再出现异常, 而且对原来没有出现异常的时刻影响也很小, 这正是希望得到的结果。

图 2

图 2. 修改切线性垂直扩散方案前异常单格点上最低层扰动变量随时间演变 (a) 地表动量通量, (b) u风场 Fig 2. Temporal evolution of perturbed variable at an abnormal grid from the original linearied vertical diffusion scheme (a) momentum flux at the surface, (b) u-wind

图 3

图 3. 修改切线性垂直扩散方案后单格点最低层扰动u风场与非线性扰动u风场随时间演变 (实线:切线性模式的结果; 虚线:非线性模式的结果; 格点位置同图 2) Fig 3. Temporal evolution of perturbed u-wind at the same location as Fig.2 after the modification of linearied vertical diffusion scheme (solid line denotes the results from the tangent-linear model, dashed line denotes those from the nonlinear model)

图 4

图 4. 切线性模式和非线性模式的平均偏差 (a) 无量纲气压, (b) u风场 (均只有动力框架和垂直扩散方案; 虚线:切线垂直扩散方案修改前; 实线:切线垂直扩散方案修改后) Fig 4. Mean departure between the tangent-linear model and the nonlinear model including vertical diffusion (a) non-dimensional pressure, (b) u-wind (dashed line denotes tangent-linear results are obtained with the original linearied vertical diffusion scheme, solid line denotes the modified scheme)

实际上, 切线性扰动应该和包括全物理过程的非线性模式中的扰动保持一致, 因为包括全物理过程的非线性模式才是实际大气正确的描述。当非线性模式包括全物理过程 (辐射过程、边界层过程、对流参数化过程等) 时, 只包括动力框架的切线性模式中的扰动和非线性扰动的偏差明显增大 (表 3~6中方案1的结果与图 1的结果对比)。总的来说, 垂直扩散方案对模式计算结果的影响较小, 但是本文修改后的切线性垂直扩散方案使得切线性扰动和非线性扰动的平均偏差在所有变量和不同时刻上都一致降低 (只有表 6中2 h的结果升高非常小的0.2%), 其中气压场和风场的平均偏差降低了1%左右, 湿度场则能够达到3.5%左右。

表 3

表 3 切线性模式扰动无量纲气压场与包括全物理过程非线性模式扰动场平均偏差 (单位:10-4) Table 3 Mean departure of the evolution of perturbed non-dimensional pressure between the tangent-linear model and the nonlinear model with full physics (unit:10-4)

表 3 切线性模式扰动无量纲气压场与包括全物理过程非线性模式扰动场平均偏差 (单位:10-4) Table 3 Mean departure of the evolution of perturbed non-dimensional pressure between the tangent-linear model and the nonlinear model with full physics (unit:10-4)

表 4

表 4 切线性模式扰动u风场与包括全物理过程非线性模式扰动场平均偏差 (单位:m/s) Table 4 Mean departure of the evolution of perturbed u-wind between the tangent-linear model and the nonlinear model with full physics (unit:m/s)

表 4 切线性模式扰动u风场与包括全物理过程非线性模式扰动场平均偏差 (单位:m/s) Table 4 Mean departure of the evolution of perturbed u-wind between the tangent-linear model and the nonlinear model with full physics (unit:m/s)

表 5

表 5 切线性模式扰动v风场与包括全物理过程非线性模式扰动场平均偏差 (单位:m/s) Table 5 Mean departure of the evolution of perturbed v-wind between the tangent-linear model and the nonlinear model with full physics (unit:m/s)

表 5 切线性模式扰动v风场与包括全物理过程非线性模式扰动场平均偏差 (单位:m/s) Table 5 Mean departure of the evolution of perturbed v-wind between the tangent-linear model and the nonlinear model with full physics (unit:m/s)

表 6

表 6 切线性模式扰动比湿场与包括全物理过程非线性模式扰动场平均偏差 (单位:10-4kg/kg) Table 6 Mean departure of the evolution of perturbed specific humidity between the tangent-linear model and the nonlinear model with full physics (unit:10-4kg/kg)

表 6 切线性模式扰动比湿场与包括全物理过程非线性模式扰动场平均偏差 (单位:10-4kg/kg) Table 6 Mean departure of the evolution of perturbed specific humidity between the tangent-linear model and the nonlinear model with full physics (unit:10-4kg/kg)

针对GRAPES四维变分同化系统的升级, 本文研究了GRAPES模式垂直扩散方案的线性化问题。通过21个个例的批量试验, 发现在非线性垂直扩散方案源代码的基础上逐句线性化而得到的切线性垂直扩散方案虽然能够通过正确性测试, 但是当初始扰动场的量级较大时, 它可能会降低切线性模式计算结果的质量, 具体表现在扰动气压场和扰动风场在个别格点上会发生异常, 一阶线性近似的误差明显增大。

相对于强非线性来说, 垂直扩散方案中开关不连续性对切线性模式计算结果影响较小。扰动气压场和扰动风场异常主要和地表动量通量的强非线性有关。针对这个问题, 解决方案是忽略可能出问题的地表动量通量扰动, 也就是将扰动地表动量通量取为零值。从批量试验结果来看, 该方案能够避免最低层扰动风场可能出现的异常, 同时没有降低正常情况下的计算精度。在其他物理过程方案的开发过程中, 该修正方案具有一定的借鉴意义。

增加切线性垂直扩散方案后, GRAPES切线性模式计算的所有扰动态在不同时刻质量得到了一致的改进, 但改进程度较小, 其中扰动气压场和扰动风场与非线性扰动偏差降低了1%左右, 相对来说, 比湿场的改进较明显, 能够达到3.5%, 这些结果说明切线性垂直扩散方案在GRAPES切线性模式中起到了积极作用。