HTML

A general form of PV (P) defined by Ertel (1942) is:
$P = \alpha {{{\xi}} _{\rm{a}}} \cdot \nabla \theta $ , where α is the specific volume of the air parcel,${{{\xi}} _{\rm{a}}}$ is the threedimensional absolute vorticity vector, and θ is the potential temperature. The PV equation proposed by Ertel (1942) is in the form$$ \dfrac{{{\rm{d}}P}}{{{\rm{d}}t}} = \alpha \left({{{{\xi}} _{\rm{a}}} \cdot \nabla \dot \theta + \nabla \times {{F}} \cdot \nabla \theta } \right), $$ (1) where
$\dot \theta $ is the diabatic heating$\left({\dot \theta = \dfrac{{{\rm{d}}\theta }}{{{\rm{d}}t}}} \right)$ ;${{F}}$ is the threedimensional frictional acceleration, where the zonal and meridional components are expressed as${F_x} = \alpha \dfrac{{\partial {\tau _x}}}{{\partial z}}$ and${F_y} = \alpha \dfrac{{\partial {\tau _y}}}{{\partial z}}$ , respectively; τ_{x} and τ_{y} are the zonal and meridional surface wind stresses, respectively; and their calculation formulas are${\tau _x} = \rho {C_{\rm{D}}}\left {{V}} \rightu$ and τ_{y} = ρC_{D}Vv, where ρ is the density of the air parcel, C_{D} is a nondimensional drag coefficient assumed to be 2 × 10^{–3}, and$\left {{V}} \right = \sqrt {{u^2} + {v^2}} $ is wind speed at 10 m. In the calculation, the height of the friction layer is taken as 1000 m, so Δz is equal to 1000 m in calculating the momentum friction.By introducing the notion of PVD (W), we obtain
$$ W = \rho P = {{{\xi}} _{\rm a}} \cdot \nabla \theta . $$ (2) According to Haynes and McIntyre (1987), W is the total PV per unit volume and is referred to as the PVD in this study (Schneider et al., 2003). Substituting Eqs. (2) into (1) and applying the continuity equation, we have
$$ \begin{aligned} \dfrac{{{\rm d}P}}{{{\rm d}t}} =& \alpha \dfrac{{{\rm d}W}}{{{\rm d}t}} + W\dfrac{{{\rm d}\alpha }}{{{\rm d}t}} = \alpha \left({\dfrac{{{\rm d}W}}{{{\rm d}t}} + W\nabla \cdot {V}} \right)\\ = & \alpha \left({{{\xi }_{\rm{a}}} \cdot \nabla \dot \theta + \nabla \times {F} \cdot \nabla \theta } \right). \end{aligned}$$ (3) By eliminating α from both sides of Eq. (3), the PVD equation can be obtained, after arranging terms, as:
$$ \dfrac{{{\rm d}W}}{{{\rm d}t}} =  W\nabla \cdot {V} + {{\xi }_{\rm{a}}} \cdot \nabla \dot \theta + \nabla \times {F} \cdot \nabla \theta . $$ (4) The above equation shows that in the Cartesian coordinate system, the rate of change of PVD following the air motion is given by the sum of the three terms on the right, i.e., the divergence term, the diabatic term, and the frictional term. Compared with the PV equation [Eq. (1)], which states that PV is conserved assuming a frictionless and adiabatic flow, PVD is not conserved and the concentration or dilution of PVD by the divergence field will cause the change of PVD.
Further simplifying Eq. (4) and applying the vector operation relationship
$\nabla \cdot \left({\varPhi {A}} \right) = \varPhi \nabla \cdot {A} + {A} \cdot \nabla \varPhi $ , the local change equation of the PVD can be obtained$$ \dfrac{{\partial W}}{{\partial t}} =  \nabla \cdot \left({W{V}  \theta \nabla \times {F}  {{\xi }_{\rm a}}\dot \theta } \right). $$ (5) The above equation indicates that the evolution of the local PVD is determined by the three terms on the right, namely, the divergence term of the PVD flux, frictional term, and diabatic heating term.
Since the expression of PV (P) is different under different coordinate systems, the expression of PVD (W) and its equation are also different in different coordinate systems. Because the isentropic coordinate system and the isobaric coordinate system are two commonly used coordinate systems in PV analysis, the following will focus on the forms of the PVD (W) and its equations in these two coordinate systems.
In the isentropic coordinate system, PV can be expressed as
$P = \dfrac{{{{\xi }_{\rm{a}}} \cdot \nabla \theta }}{{{\rm{\rho }}_\theta }}$ , where${{\rm{\rho }} _\theta } =  \dfrac{1}{g}\dfrac{{{\rm d}p}}{{{\rm d}\theta }}$ is the “density” in (x, y, θ) space (Holton, 2004). Introducing the PVD (W = ρ_{θ}P), then the following relationship exists in the isentropic coordinate system:$W = {{\xi }_{\rm{a}}} \cdot \nabla \theta = {{\xi }_{\rm{a}}} \cdot {k} = f + {\zeta _\theta }$ , indicating that the PVD is numerically equal to the isentropic absolute vorticity in the isentropic coordinate system (Rossby, 1940; Hoskins et al., 1985; Haynes and McIntyre, 1987). Similarly, by using this relationship and the continuous equation in the isentropic coordinate system, the equation of the PVD in the isentropic coordinate system can be obtained,$$ \dfrac{{\partial W}}{{\partial t}} =  {\nabla _\theta } \cdot \left[ {\left({f + {\zeta _\theta }} \right){V}  {F} \times {k} + \dot \theta \dfrac{{\partial {V}}}{{\partial \theta }} \times {k}} \right], $$ (6) where in the isentropic coordinate system, the vertical component (
$  \dfrac{\partial }{{\partial \theta }}\left[ {\left({f + {{\rm{\zeta }}_\theta }} \right)\dot \theta } \right]$ ) of the divergence term of PVD flux [$  \nabla \cdot \left({W{V}} \right)$ ] and the vertical component ($\dfrac{\partial }{{\partial \theta }}\left[ {\left({f + {{\rm{\zeta }}_\theta }} \right)\dot \theta } \right]$ ) of the diabatic heating term [$\nabla \cdot \left({{{\xi }_{\rm a}}\dot \theta } \right)$ ] are just offset, and the friction term$\nabla \times {F} \cdot \nabla \theta = \nabla \times {F} \cdot {k} = $ $ {\nabla _{\rm{\theta }}} \cdot \left({{F} \times {k}} \right)$ . Therefore, the evolution of the local PVD follows a horizontal motion on an isentropic surface. Equation (6) shows that in the isentropic coordinate system, the local variation of the PVD is related to the horizontal convergence or divergence of the PVD flux, the frictional effect, and the diabatic heating along the isentropic surface.The isobaric coordinate system is widely used in synoptic analysis. An infinitesimal control volume in isobaric coordinate with crosssectional area δA and vertical extent δp has a mass:
${\rm{\delta }}M = \rho_z{\rm{\delta }} z{\rm{\delta }}A = \left({  \dfrac{{\delta p}}{g}} \right){\rm{\delta }}A = $ ${\rho _{p}}\left({  {\rm{\delta }}p} \right){\rm{\delta }}A $ . Therefore, the “density” in (x, y, p) space is ρ_{p} = –ρ_{z}$\dfrac{{{\rm d}z}}{{{\rm d}p}} $ = g^{–1}, and the expressions of the PV (P) and PVD (W) in the isobaric coordinate system are$$ \left\{ {\begin{aligned} & {{P_p} =  g{{\xi }_{{\rm a}p}} \cdot {\nabla _p}\theta =  g\left[ {{k} \times \dfrac{{\partial {{V}_{\rm h}}}}{{\partial p}} \cdot {\nabla _{p{\rm h}}}\theta + \left({f + {\zeta _p}} \right)\dfrac{{\partial \theta }}{{\partial p}}} \right]}\\ & {{W_p} =  {{\xi }_{{\rm a}p}} \cdot {\nabla _p}\theta =  \left[ {{k} \times \dfrac{{\partial {{V}_{\rm h}}}}{{\partial p}} \cdot {\nabla _{p{\rm h}}}\theta + \left({f + {\zeta _p}} \right)\dfrac{{\partial \theta }}{{\partial p}}} \right]} \end{aligned}} \right.. $$ By comparing PV and PVD in the isobaric coordinate system, it can be found that, unlike in the isentropic coordinates, PV and PVD in the isobaric coordinate system are the same except for a factor g. PVD retains the dynamic and thermal properties of PV. Similarly, the local variation equation of PVD in the isobaric coordinate system can be obtained as follows,
$$ \dfrac{{\partial W}}{{\partial t}} =  {V} \cdot {\nabla _p}W  {{\xi }_{{\rm a}p}} \cdot {\nabla _p}\dot \theta  {\nabla _p}\theta \cdot {\nabla _p} \times {F}. $$ (8) The above equation shows that in the isobaric coordinate system, the local variation of the PVD is related to the PVD advection, diabatic heating, and frictional effects. In the case of an adiabatic and frictionless flow, Eq. (8) is reduced to
$\dfrac{{\partial W}}{{\partial t}} + { V} \cdot \nabla W = \dfrac{{dW}}{{dt}} = 0$ , indicating that the PVD is conserved following the flow motion. Like that PV is conserved following a horizontal motion on an isentropic surface with the absence of the diabatic heating and friction, PVD is conserved following a horizontal motion on an isobaric surface when the vertical velocity equals zero (ω = dp/dt = 0), which means that PVDanalysis in isobaric coordinate can be used to reveal the dynamic and thermal features of the atmospheric circulation following the quasihorizontal motions. However, when considering a highimpact synoptic process, with vertical velocity developing rapidly, it is necessary to diagnose the effect of the air ascending motion. 
Integrating Eq. (6) over an isentropic surface θ_{s} and employing Gauss’s theorem for the integration of the divergence of the PVD flux lead to
$$ \begin{aligned} & \dfrac{\partial }{{\partial t}}\int_S {\sigma P{\rm{d}}s} = \dfrac{\partial }{{\partial t}}\int_S {W{\rm{d}}s} \\ & =  \oint_\varGamma {\left[ {\left({f + {\zeta _\theta }} \right){V}  {F} \times {k} + \dot \theta \dfrac{{\partial {V}}}{{\partial \theta }} \times {k}} \right] \cdot {n}{\rm{d}}l}, \end{aligned} $$ (9) where σ = ρ_{θ}, S is the total area of the isentropic surface θ_{s} for the integration, Γ is the boundary of the θ_{s} surface intersecting with the ground surface, dl is the unit increment of the boundary Γ, which is positive when the linear integration takes an anticlockwise path, and n is a unit vector perpendicular to dl and oriented to its right side. For an isentropic layer with unit thickness (dθ = 1) and centered at the θ_{s} surface, the above equation defines the relationship between PV generation and PVD, i.e., the mass weighted PV increment of the isentropic layer is related to the distributions of PVD flux, friction, and diabatic heating along the boundary of the isentropic layer. For the isentropic surfaces in the overworld and the middleworld, since the isentropic surface does not intersect with the ground, there is no boundary (Γ→0) and the integral result must be zero. That is, the PV is conserved on these isentropic surfaces. However, the integral result in the underworld is related to the magnitude of the PVD transport, friction, and diabatic heating at the boundary. As a result, there will be PV generation at the boundary of the isentropic surface. In addition, Eq. (9) also shows that the change of the PVD on the boundary of the isentropic surface can intrude into the isentropic surface by the convergence or divergence of the PVD flux, and then influence the interior circulation.
By dividing the boundary Γ into N equal segments (
$\varGamma = \sum\nolimits_{n = 1}^N {{\varGamma _n}} $ ), Eq. (9) is transformed into$$ \dfrac{\partial }{{\partial t}}\int_S {W{\rm{d}}s} = \sum\nolimits_{n = 1}^N {\int_{{\varGamma _n}} {A{\rm{d}}l} } , $$ (10) where
$A =  \left[( {f + {\zeta _\theta }){V}  {F} \times {k} + \dot \theta \dfrac{{\partial {V}}}{{\partial \theta }} \times {k}} \right]\cdot{n}$ represents the normal component of the PVD flux along the isentropic surface boundary. Under the framework of the geostrophic balance,$\dfrac{{\partial {V}}}{{\partial \theta }} =  \dfrac{R}{{fp}}\dfrac{{\partial p}}{{\partial \theta }}{k} \times \nabla T$ . Therefore, the contribution of surface heating to A can be expressed as$$ {A_{{Q}}} =  \left({\dot \theta \dfrac{{\partial {V}}}{{\partial \theta }} \times {k}} \right) \cdot {n} =  \dfrac{g}{{f{\theta _z}}}\dfrac{{\dot \theta }}{\theta }{\theta _n}, $$ (11) where θ_{z} is the static stability, and θ_{n} is the surface potential temperature gradient, indicating the baroclinicity. Equations (10) and (11) then show that in the areas where the surface baroclinicity is stronger, there is more contribution of surface diabatic heating to the surface PV generation.
Suppose that at one moment, the lowest isentropic surface intersecting with the ground is θ_{l}, and the top isentropic surface is θ_{M}, then the total surface PV generation at that moment can be obtained by spatial integration of Eq. (10) from θ_{l} to θ_{M}
$$ \begin{aligned} & {{G}_{\rm{total}}}\left({{\theta }_{1}}, {{\theta }_{\rm M}} \right)=\int_{{{\theta }_{1}}}^{{{\theta }_{\rm M}}}{\left[ \dfrac{\partial }{\partial t}\int_{S}{W{\rm{d}}s} \right]{\rm{d}}\theta } \\ & = \sum\nolimits_{n=1}^{N}{\int_{{{\theta }_{1}}}^{{{\theta }_{\rm M}}}{\left[ \int_{{{\varGamma }_{n}}}{A{\rm{d}}l} \right]{\rm{d}}\theta }}=\dfrac{\partial }{\partial t}\iiint_{V}{P{\rm{d}}m}, \end{aligned}$$ (12) where dm = ρ_{θ}dxdydθ. Since the isentropic surfaces intersecting with the ground always change with time, G_{total} is a function of time. It can be seen from Eqs. (10)–(12) that the larger the value of A, the larger the value of G_{total}, indicating that the total amount of surface PV generation is related to the PVD flux, and frictional and diabatic heating effects at the isentropic surface boundary. In addition, the total amount of surface PV generation is also related to the distribution of the surface potential temperature (A_{Q}) in the integral region. The denser the surface potential temperature in the integral region (the larger the θ_{n}), the larger the value of A, and the more the surface PV generation.

The data used in this study are 3h averaged data in January and February 1980–2017 from MERRA2 (the second ModernEra Retrospective Analysis for Research and Applications) reanalysis dataset (Gelaro et al., 2017), provided by NASA Global Modeling and Assimilation Office. Three collections of the dataset are used here, including the basic assimilated meteorological fields on both pressure levels (“inst3_3d_asm_Np” also known as “M2I3NPASM”) and model levels (“tavg3_3d_asm_Nv” also known as “M2T3NVASM”), as well as the singlelevel atmospheric state variables (“tavg1_2d_slv_Nx” also known as “M2T1NXSLV”).
According to Sun and Zhao (2010), after 18 January, the convective systems over the TP became extremely active and frequently moved eastward. In order to highlight the role of surface PVD forcing over the TP, the extreme winter freezing rain/snow storm disaster from 18 January to 2 February 2008 is selected for this study and the time period is defined as the “storm period” in the following context. By comparing the difference of various fields during the storm period with the climatology from January to February over 1980–2017, the impact of surface PVD forcing over the TP on the extreme winter freezing rain/snow storm disaster is explored.
According to the surface PV generation theory introduced in Section 2.2, the isentropic coordinate system has a great advantage for PV and PVD analyses. However, it is difficult to calculate the surface PV and PVD in the isentropic coordinate system, which is mainly related to the complex distribution of the isentropic surfaces in the boundary layer. For example, due to the strong heating and cooling effects on the TP surface, the distribution of isentropic surfaces near the TP surface is usually tilted; the turbulent mixing in the boundary layer sometimes forms a constant flux layer near the ground surface, and the isentropic surface is almost perpendicular to the ground surface; in some areas, there is an unstable hydrostatic stratification (
$\dfrac{\partial \theta }{\partial z} < 0$ ), which disables the PVD diagnosis. In addition, since the isentropic surfaces intersect with the ground suface, when calculating the forcing terms on the right hand side of Eq. (6) at an isentropic surface boundary, it is necessary to use both surface data and isentropic surface data. Since most of the reanalysis data only have a few typical levels of isentropic surfaces, plenty of interpolations are needed, which may introduce inestimable errors. The same problem exists with the Cartesian coordinate system and the isobaric coordinate system. In order to avoid the introduction of errors, the forms of PV, PVD, and PVD equation are transformed to the σ–p coordinate system, so that each term in the PVD tendency equation can be directly calculated by using the data at model levels from MERRA2. After obtaining the PVD tendency at each model level, the PVD tendency in other coordinate systems can be obtained by coordinate transformation. The forms of PV, PVD, and PVD equation in the σ–p coordinate system are$$ \left\{ \begin{matrix} P=\dfrac{g}{{{p}_{\rm s}}}{{{\xi}}_{{\rm a}\sigma }}\cdot {{\nabla }_{\sigma }}\theta =\dfrac{g}{{{p}_{\rm s}}}\left[ {k}\times \dfrac{\partial {{{V}}_{\rm h}}}{\partial \sigma }\cdot {{\nabla }_{\sigma {\rm h}}}\theta +\left(f+{{\zeta }_{\sigma }} \right)\dfrac{\partial \theta }{\partial \sigma } \right] \\ \!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\! W={{\xi}_{{\rm a} \sigma }}\cdot {{\nabla }_{\sigma }}\theta =\left[ {k}\times \dfrac{\partial {{{V}}_{\rm h}}}{\partial \sigma }\cdot {{\nabla }_{\sigma {\rm h}}}\theta +\left(f+{{\zeta }_{\sigma }} \right)\dfrac{\partial \theta }{\partial \sigma } \right] \\ \!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\! \!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\! \!\!\!\!\!\!\! \!\!\!\!\!\!\! \!\!\!\!\!\!\! \dfrac{\partial W}{\partial t}={{\nabla }_{\sigma }}\cdot \left(W{V}+\theta \nabla \times {F}+{{\xi}_{{\rm a}\sigma }}\dot{\theta } \right) \\ \end{matrix} \right., $$ (13) where p_{s} is the surface pressure, and
${\nabla _{\rm{\sigma }}}$ is the threedimensional gradient operator in the σ–p coordinate system:${\nabla _{\rm{\sigma }}} = \dfrac{\partial }{{\partial x}}{i} + \dfrac{\partial }{{\partial y}}{j} + \dfrac{\partial }{{\partial \sigma }}{k}$ .