# Comparison of the Atmospheric 200 hPa Jet’s Analyses between Proper Orthogonal Decomposition and Advanced Dynamic Mode Decomposition Method

Authors: Mei Gao; Xiao-Qun Cao; Bai-Nian Liu; Zi-Hang Han; Shi-Cheng Hou; Guo-Gui Yang
Journal: Advances in Meteorology
 (2020)

---

## Abstract

In this paper, a frequently employed technique named the sparsity-promoting dynamic mode decomposition (SPDMD) is proposed to analyze the velocity fields of atmospheric motion. The dynamic mode decomposition method (DMD) is an effective technique to extract dynamic information from flow fields that is generated from direct experiment measurements or numerical simulation and has been broadly employed to study the dynamics of the flow, to achieve a reduced-order model (ROM) of the complex high dimensional flow field, and even to predict the evolution of the flow in a short time in the future. However, for standard DMD, it is hard to determine which modes are the most physically relevant, unlike the proper orthogonal decomposition (POD) method which ranks the decomposed modes according to their energy content. The advanced modal decomposition method SPDMD is a variant of the standard DMD, which is capable of determining the modes that can be used to achieve a high-quality approximation of the given field. It is novel to introduce the SPDMD to analyze the atmospheric flow field. In this study, SPDMD is applied to extract essential dynamic information from the 200 hPa jet flow, and the decomposed results are compared with the POD method. To further demonstrate the extraction effect of POD/SPDMD methods on the 200 hPa jet flow characteristics, the POD/SPDMD reduced-order models are constructed, respectively. The results show that both modal decomposition methods successfully extract the underlying coherent structures from the 200 hPa jet flow. And the DMD method provides additional information on the modal properties, such as temporal frequency and growth rate of each mode which can be used to identify the stability of the modes. It is also found that a fewer order of modes determined by the SPDMD method can capture nearly the same dynamic information of the jet flow as the POD method. Furthermore, from the quantitative comparisons between the POD and SPDMD reduced-order models, the latter provides a higher precision than the former, especially when the number of modes is small.


---

## Body

## 1. Introduction

With the rapid development of computer technology and continuous improvement of global observing system’s performance in recent decades, many new reanalysis datasets (also known as meteorological big data) from advanced operational and research centers are built, upgraded, and opened to the meteorology community. Nowadays, it is still a great challenge to gain maximum valuable information and find new motion mechanisms from reanalysis datasets using a new method or powerful algorithm, although many scientists have made much progress. It is very important to accurately describe and understand the structural changes and instability mechanisms of the complex atmospheric motion in both operational numerical weather predictions and scientific researches [1]. However, in the instability analysis methods developed, when analyzing the stability of the atmospheric fluid motion, it is necessary to greatly simplify the Navier–Stokes equations of atmospheric motion or to calculate the inverse and eigenvalues of large matrices; it is well known that the process of calculations is complex and costly [2]. While reanalysis datasets are available, they can help us to get rid of this current dilemma. So, the objective of the study can be changed from the Navier–Stokes equations to high-resolution meteorological fields such as temperature, velocity, vorticity, etc., because of their accuracy and near realness. To gain a deeper understanding of the complex atmospheric fluid motion, modal decomposition methods can be utilized; among several methods, the Proper Orthogonal Decomposition (POD) and the Dynamic Mode Decomposition (DMD) have been widely used to study the dynamic mechanism of the flows in different applications [3, 4], except for the analysis of atmospheric reanalysis datasets.The POD was introduced to the fluid dynamics community first by Lumley [5] as a mathematical technique to extract energetically and dynamically important features of the fluid flows that cannot be easily understood through direct inspection. The method itself also known as empirical eigenfunction decomposition (EOF) is applied to meteorological data [6]. In POD, the flow field is decomposed into spatial orthogonal basis functions (or modes), where the modes are arranged in terms of their energy content and truncated as required [7]. The POD method is easy to implement, and there is an explicit ranking criterion for the decomposed modes; thus, such technique has been extensively applied in numerous studies. For instance, several authors used the POD to extract coherent structures from the turbulent flows and oceanography [8–10]. Cazemier et al. [11] analyzed the flow in a cavity based on the POD. Perrin et al. [12] used POD to obtain phase averaged turbulence properties for flow past a cylinder. The orthogonal property of POD modes is useful for reconstructing the flow field which also has been shown in [13, 14]. However, the POD method is at its heart a second-order statistical method; value phase information is ignored, and the modes do not contain the dynamics encoded in the datasets. These problems can be overcome by the DMD method, which was recently developed by Schmid and Sesterhenn [15] as the basis of a Koopman spectral analysis of a dynamic system promoted by Mezić [16] and Rowley et al. [17] and has provided a compact and instructive manner for understanding the dynamic information of flow.In DMD, the flow field is decomposed into a set of modes whose dynamics are governed by their corresponding eigenvalues, and the isolated eigenvalues define the temporal frequency and growth/decay rate of each mode [18]. DMD makes almost no assumptions about the underlying system; it is an equation-free and data-driven method and isolates dynamic mode with a particular frequency unlike the POD mode with mixed frequencies. Much of the recent work involving DMD has focused on its application in fluid mechanics. Schmid and Sesterhenn [15] applied the DMD method on both numerical and experimental data, in which the plane channel flow and a jet passing two cylinders are investigated. Rowley et al. [17] performed a DMD analysis on the jet in the crossflow, and two dominant shedding frequencies are found. Muld et al. [19] applied POD and DMD to extract the most dominant flow structures of a simulated flow in the wake of a high-speed train model. Other works combined the two decomposition methods shown by Schmid et al. [20] and Goza and Colonius [21]. In many situations, it can be difficult to determine which modes are the most dynamic relevant to the flow field in the standard DMD, as there is no explicit criterion to rank eigenvalue importance, unlike the POD modes are arranged according to their energy content. For example, spatial nonorthogonality of the DMD modes may introduce a poor quality of approximation of experimentally or numerically generated snapshots when only a subset of modes with the largest amplitude is retained [23]. And there are different rank criteria of DMD modes that have been produced [23, 24] and some improved variants of standard DMD have been proposed including optimized DMD [25] and gradient-based DMD [26] to solve the problem. To strike a balance between the quality of approximation (in the least-squares sense) and the number of modes that are used to approximate the given fields, in this paper, we focus on the advanced modal decomposition method SPDMD, first proposed by Jovanovic et al. [27]. SPDMD method attempts to keep a balance between the accuracy of the approximation and the number of modes used to reconstruct the flow fields, which is achieved by using a least-squares deviation norm (between the matrix of snapshots and the reconstructed data sequence using the DMD modes) augmented with a l1-norm regularization term to specify the emphasis on the sparsity of the reduced-order model. The application of the SPDMD method to the analysis of the atmospheric motion flow field, especially the analysis of the weather system of a stratospheric jet, has not been reported publicly. In this paper, the SPDMD method is proposed to implement the instability analysis of the stratospheric 200 hPa jet using reanalysis datasets, which can gain many dynamical stable/unstable modes and to determine the modes that have the most profound influence on the jet flow. And the SPDMD reduced-order model is constructed based on the selected dynamical modes. Thus, the analysis of the whole flow field shows how to characterize the evolutions of several modes.The remainder of this article is organized as follows. In Section2, we recalled the principles governing the POD, DMD, and SPDMD methods and introduced the implementation of the POD and SPDMD algorithms for the stratospheric 200 hPa jet. The results and discussions of POD and SPDMD mode decomposition and reconstructions for the stratospheric 200 hPa jet are presented in Section 3. The paper ends with the conclusions in Section 4.


## 2. Materials and Methods

### 2.1. Proper Orthogonal Decomposition

Herein, a short review of the POD method is provided, whose details can be found in [28]. Considering the following dynamical system,(1)dxdt=fx,t;η,where xt∈ℝm is a vector representing the state of our dynamical system at time t, η contains parameters of the system and f⋅ represents the dynamics. Generally, the dynamic system is often nonlinear, and the dimension of state x is typically quite large m≫1. In practice, we adopt a corresponding discrete-time system representation induced from the continuous-time dynamics of equation (1), as we are ultimately interested in analyzing the discrete-time data:(2)xi+1=Fxi,where map F advances the states forward every Δt in time according to the dynamical system, and the subscript i denotes the time index i∈ℕ, so that xi=xiΔt.In applications of POD to fluid flow, we start with the vector fieldxt with its time-averaged value x¯t removed, then collect a sequence of snapshots from the fluctuating component of the vector field x′t, and stack them into a data matrix whose columns represent the individual data samples:(3)X′=x1′x2′x3′…xn′∈ℝm×n,in which the snapshots are collected at times ti from i=1,2,3,…,n for a total of n measurement times.The POD method aims to seek the optimal set of basis functions that can best represent the given flow field data, which can be obtained by solving the following eigenvalue problem:(4)Cφj=λjφj,where C is the covariance matrix of vector xi′, which expresses the difference between any two snapshots:(5)C=∑inx′tix′Tti=X′X′T∈ℝm×m.The solution of equation (4) gives the orthogonal eigenfunctions φj called the POD modes, and the corresponding eigenvalues λj are real and describe the energy content of the coherent flow structure represented by the eigenfunctions φj. The flow field snapshots can be represented as the linear combination of the POD modes φjwith their corresponding expansion coefficientsaj, where aj can be obtained by projecting the snapshot onto the POD modes,(6)aj=ΦTxj′,j=1,2,3,…,n,in which Φ=φ1,φ2,…,φn.Generally, only a subset of the POD modes retained to express the flow; the energy of these modes can be presented byEPOD, defined as [29](7)EPOD=∑j=1rλj∑j=1nλj,where the total energy is normalized such that ∑j=1nλj=1, and r is the number of the retained POD modes. In practice application, the truncated POD modes with energies that sum to 99% of the total energy are a plausible truncation.With the determination of the relative important POD modes, the flow field can be represented in a reduced dimensionality:(8)xj′≈Φaj,j=1,2,3,…,r.For practical applications, the decomposition can be carried out using the classical method or the snapshot method. In matrix form, the POD can be realized by the standard singular value decomposition (SVD) which is a robust approach to determine the POD modes [30]. In this paper, we use the SVD to implement the POD decomposition; the specific details are introduced below.


### 2.2. Computation of the Proper Orthogonal Decomposition

#### 2.2.1. Preprocessing

In the field of atmospheric/marine science, whether the forecast field data generated by the high-resolution numerical forecasting model or the analysis field produced by the global variational data assimilation system, they are distributed in the regular space grid at regular intervals. And they only describe the state of atmospheric motion at a certain moment, that is, strictly speaking, an instantaneous state field (actually, the state changes little before and after the moment) which is equivalent to store using the camera after taking pictures of the atmospheric motion state. Therefore, it can be named as the atmospheric motion snapshot analogy to the concept of flow snapshot in fluid mechanics. Thus, in this paper, we use the modal decomposition methods to analyze the snapshots of atmospheric motion over some time.Here, we collect data from a two-dimensional stratospheric jet flow velocity fieldu,v and construct snapshot matrix X in preparation to perform the modal decomposition analysis. The matrix X is constructed by the velocity components u,v of each measurement point of the M instantaneous stratospheric jet flow fields; each column of X contains all the velocity data at each time instance; hence, the matrix X is defined as follows:(9)X=u11u12⋯u1N⋮⋮⋮⋮uM1uM2⋯uMNv11v12⋯v1N⋮⋮⋮⋮vM1vM2⋯vMN,where the apex 1,2,3,…,N indicates the time number and the term 1,2,3,…,M is the number of measurement points.


#### 2.2.2. Implementation of POD

The key idea of the POD is to reduce the data dimensions by projecting high-dimensional data to a low-dimensional space while reserving the variation in the original variables as much as possible. It can be seen as an application of singular value decomposition (SVD) in a finite-dimensional space, and a brief review of the algorithm is shown as follows.(1)
Consider the given snapshot matrixX∈ℝM×N. Firstly, compute the mean flow X¯=1/N∑i=1Nxi and the mean-subtracted snapshot matrix X′=X−X¯.(2)
The snapshots matrixX′ can be decomposed with SVD as(10)X′=ΦΣΨ∗,whereΦ∈ℝM×N, Σ∈ℝN×NandΨ∈ℝN×N, matrix Σ holds the singular values σjj=1,2,…,N along its diagonal, that define how much “energy” contained in this mode, and the asterisk denotes the conjugate transpose. Besides, the left and right singular vectors of X are the eigenvectors of XXT and XTX, respectively. Moreover, the eigenvalues of XXTXTX and singular values of X are related by λj=σj2. Thus, the SVD can be directly performed on X to determine the POD modes Φ and corresponding eigenvalue and provide time-evolving information of the flow field contained in the matrix Ψ.(3)
In practical applications, one usually retainsrPOD several modes that capture 99% of the snapshots energy to achieve a high quality of model reduction, such that(11)EPOD=∑j=1rPODσj2∑j=1Nσj2≈99%.(4)
With the determination of the most important POD modes, the jet flow field can be represented only using the truncated modes, as follows:(12)X≈X¯+∑j=1rPODDjσj.To obtain the matrixDj, we multiply the corresponding columns of Φ and Ψ with an outer product. In the MATLAB notation, the computation of Dj is the following: Dj=Φ:,jV:,jT.


### 2.3. Sparsity-Promoting DMD

In this section, we first give a brief overview of the Koopman operator theory, which is the basis for DMD analysis, and then describe the basic procedure for the standard DMD and sparsity-promoting DMD method; the latter one leads to an optimal selection of the original DMD modes.First, we also consider the discrete-time dynamical system mentioned in Section2.1, in equation (2) xi+1=Fxi.The Koopman operatorK proposed by Koopman [18] in 1931 is an infinite-dimensional linear operator that operates on any scalar-valued function g in the following manner:(13)Kgx=gFx.The Koopman operator steps forward in time an observable. Note thatK is a linear operator; it is supposed that the operator K can be decomposed into eigenvalues μj and corresponding eigenfunctions ϕj, which satisfy(14)Kϕj=μjϕj,j=1,2,3,….According to Rowley et al. [17], we consider a vector-valued function g that lies within the span of eigenfunctions of Koopman operator, and g can be expanded in terms of these eigenfunctions as(15)gx=∑j=1∞ϕjxυj.Thus, the functiongx was expanded as a linear combination of the eigenfunctions ϕjof K, where υj are the vector coefficients which are referred to as Koopman modes of the map F, corresponding to the observable g.Under the action of Koopman operator, the evolution of vector functiongx can be expressed as follows:(16)gxi=Kigx0=∑j=1∞μjiϕjx0υj.Thus, the eigenvalueμj characterizes the temporal behavior of the Koopman modes υj, while the phase of μj determines its frequency, and the magnitude determines the growth rate. These are very important information for the stability analyses of the modes, as when the growth rate is positive, the corresponding mode is unstable, while when the growth rate is negative, the corresponding mode is stable, and the mode is neutrally stable when the growth rate is just zero. Note that, as described in [16], for a system evolving on an attractor, the Koopman eigenvalues always lie on the unit circle.In DMD architecture, consider finite measurements at timesti of the dynamic system in equation (2):(17)yi=gxi,i=1,2,3,…,N.In practice application, the state can be directly measured; thus,yi=xi and the DMD method [15] assumes that there is a locally best-fit linear operator A to step forward in time the snapshots(18)xi+1=Axi.It is easy to prove that the eigenvalues of the operatorA are also ones of the Koopman operator, and the Koopman modes coincide with the eigenvectors of A; thus, DMD method can be used for estimating an approximation of the decomposition in equation (16).Arrange theN snapshots into two data matrices, as follows:(19a)X1N−1=x1x2x3…xN−1,(19b)X2N=x2x3…xN.The locally linear approximation in equation (18) may be written in terms of these data matrices as(20)X2N≈AX1N−1.The matrixA can be given as(21)A=X2NX1N†,where † represents the Moore–Penrose pseudoinverse. As already highlighted, the dimensions of the snapshots xi are large m≫1; thus, the matrix A is difficult to represent or analyze directly.For this reason, DMD circumvents the eigendecomposition ofA by considering a rank-reduced representation in terms of a POD-projected matrixA˜ based on singular value decomposition (SVD) of X1N−1:(22)X1N−1=UΣV∗,where U∈CM×r, Σ∈Cr×r, and V∈CN×r. The left singular vectors U are the POD modes, matrix Σ is a diagonal matrix with the singular values σ1,σ2,…,σr that define how much “energy” contained in this mode，and r is the rank of X1N−1.Thus, the matrixA can be obtained as(23)A=X2NVΣ−1U∗.By projecting the high-dimensional matrixA onto POD, modes can obtain the representation of A˜ as(24)A˜=U∗AU=U∗X2NVΣ−1.The matrixA˜ defines a low-dimensional linear model of the dynamical system on POD coordinates, in which(25)x˜i+1=A˜x˜i,and the high-dimensional state xi can be reconstructed as xi=Ux˜i.The matrixA˜ is a similar transformation of A, so that the matrix A˜ contains the main eigenvalues of A; compute the eigendecomposition of A˜ as(26)A˜W=WΛ,where columns of W are eigenvectors and Λ is a diagonal matrix containing the corresponding eigenvalues μj.Therefore, the snapshot at a timeti can be estimated as(27)xi=UWΛi−1W−1U∗x1=∑j=1rμji−1αjυj,in which υj=UW are the DMD modes and αj is the amplitude of the jth DMD mode that represents the jth modal contribution of the initial condition defined as α=W−1x˜1=W−1U∗x1.Rewrite equation (27) in matrix form as(28)x1x2…xN−1≈υ1υ2…υrα10⋯00α2⋯0⋮⋮⋱⋮00⋯αr︸Dα=diagα1μ1⋯μ1N−11μ2⋯μ2N−1⋮⋮⋱⋮1μr⋯μrN−1︸Vand,in which the temporal evolution information of the DMD modes is contained in the Vandermonde matrix Vand.As there are no explicit ways for DMD to select modes, unlike the POD modes which are selected based on the rank of energy from the singular values. The sparsity-promoting DMD method is proposed by Jovanovic et al. which induces a sparse structure of the DMD modes by augmenting thel1-norm of the vector α through the above equation. That is, it is expressed as the following optimization problem:(29)minαX1N−1−υDαVandF2+γ∑jrαj,in which γ is a positive regularization parameter that reflects the emphasis on the sparsity of the amplitudes vector α, larger values of γ emphasize on the number of nonzero elements in the vector α, and ⋅F2 is the Frobenius norm. And the values of γ are selected in such a way that γmin induces a dense vector α (with 191 nonzero elements), and γmax induces α with a single nonzero element.After the sparse structure ofα is fixed, the problem becomes to determine only the nonzero amplitudes as the solution from the optimization problem in equation (29) with the constraint ETα=0, where the matrix E contains the sparsity structure of α. As α=0α20α4T, the matrix E is given by E=10000010T.A powerful and well-suited algorithm named Alternating Direction Method of Multipliers (ADMM) for solving the optimization problems of equation (29) can be found in [27].


### 2.4. Computation of the Sparsity-Promoting DMD

The sparsity-promoting DMD is a variant of the standard DMD method; here, the main steps of the SPDMD algorithm are summarized as follows:(1)
First, a matrixX1N−1 is formed with the first N−1 columns and matrix X2N contains the last N−1 columns of the snapshot matrix X shown in equation (9) And then perform the standard DMD for the matrices to obtain DMD solution A˜ and eigenvector matrix W and corresponding eigenvalue matrix Λ (and also U, V, and Σ).(2)
The SPDMD aims to solve the optimization problem in equation (29). We apply the alternating direction method of multipliers (ADMM) to solve this convex problem. ADMM is a state-of-the-art algorithm for solving large-scale and distributed optimization problems [31]. And there are four steps for the application of the ADMM.Step 1. Replace the vector of amplitude α in the l1-regularization with a new variable β; the optimization problem in equation (29) is transformed into the following constrained optimization problem:(30)minimizeα,βJα+γhβsubject toα−β=0,whereJα=minαX1N−1−υDαVandF2, and hβ=∑jrβj1. β is an auxiliary variable; for any feasible α and β, the optimization problems in equations (29) and (30) are equivalent.Step 2. To solve the problem, the augmented Lagrangian is introduced to transform the constraints in equation (30) to the objective function, as follows:(31)Lρα,β,λ=Jα+γhβ+12θ∗α−β+α−β∗θ+ρα−β22,hereθ is the vector of Lagrange multipliers and ρ is a positive parameter that introduces a quadratic penalty on the deviation between α and β, and ⋅2 is the Euclidean norm of a given vector.Step 3. ADMM algorithm is an iterative algorithm that combines the Dual Ascent and Multiplier method. It includes a α-minimization step, a β-minimization step, and a Lagrange multiplier update step; the details are shown as follows:(32a)αk+1=argminαLρα,βk,θk,(32b)βk+1=argminβLραk+1,β,θk,(32c)θk+1=λk+ραk+1−βk+1,in which theα-minimization step amounts to solving an unconstrained regularized quadratic program and that the β-minimization step amounts to an opportune use of the soft thresholding operator.α-Minimization step: completion of squares concerning α in the augmented Lagrangian Lρ can be used to show that the α-minimization step (33c) is equivalent to(33)minαJα+ρ2α−u2k2,uk=βk−1ρθk,in whichk is the step number of iterations. The solution is determined [27] by(34)αk+1=P+ρ2I−1q+ρ2uk,in whichP=W∗W∘VandV∗¯, q=diagVandVΣ∗W¯. And I is the unit matrix, an overline signifies the complex-conjugate of a vector or matrix, and ∘ is the elementwise multiplication of two matrices in this case [32].β-Minimization step: the completion of squares for β in the augmented Lagrangian Lρ can be used to show that the β-minimization step (32b) is equivalent to(35)minβγhβ+ρ2β−v2k2,vk=αk+1+1ρθk,and the solution is determined by(36)βjk+1=Sκvjk,κ=γρ,whereSκ⋅ denotes the soft thresholding operator:(37)Sκvjk=vjk−κ,vjk>κ,0,vjk∈−κ,κ,vjk−κ,vjk>−κ.Step 4. Start with an initial point α0,β0 and iterate until the required error is εprim and εdual then stop:(38)αk+1−β2k+1≤εprim,βk+1−β2k≤εdual.


## 2.1. Proper Orthogonal Decomposition

Herein, a short review of the POD method is provided, whose details can be found in [28]. Considering the following dynamical system,(1)dxdt=fx,t;η,where xt∈ℝm is a vector representing the state of our dynamical system at time t, η contains parameters of the system and f⋅ represents the dynamics. Generally, the dynamic system is often nonlinear, and the dimension of state x is typically quite large m≫1. In practice, we adopt a corresponding discrete-time system representation induced from the continuous-time dynamics of equation (1), as we are ultimately interested in analyzing the discrete-time data:(2)xi+1=Fxi,where map F advances the states forward every Δt in time according to the dynamical system, and the subscript i denotes the time index i∈ℕ, so that xi=xiΔt.In applications of POD to fluid flow, we start with the vector fieldxt with its time-averaged value x¯t removed, then collect a sequence of snapshots from the fluctuating component of the vector field x′t, and stack them into a data matrix whose columns represent the individual data samples:(3)X′=x1′x2′x3′…xn′∈ℝm×n,in which the snapshots are collected at times ti from i=1,2,3,…,n for a total of n measurement times.The POD method aims to seek the optimal set of basis functions that can best represent the given flow field data, which can be obtained by solving the following eigenvalue problem:(4)Cφj=λjφj,where C is the covariance matrix of vector xi′, which expresses the difference between any two snapshots:(5)C=∑inx′tix′Tti=X′X′T∈ℝm×m.The solution of equation (4) gives the orthogonal eigenfunctions φj called the POD modes, and the corresponding eigenvalues λj are real and describe the energy content of the coherent flow structure represented by the eigenfunctions φj. The flow field snapshots can be represented as the linear combination of the POD modes φjwith their corresponding expansion coefficientsaj, where aj can be obtained by projecting the snapshot onto the POD modes,(6)aj=ΦTxj′,j=1,2,3,…,n,in which Φ=φ1,φ2,…,φn.Generally, only a subset of the POD modes retained to express the flow; the energy of these modes can be presented byEPOD, defined as [29](7)EPOD=∑j=1rλj∑j=1nλj,where the total energy is normalized such that ∑j=1nλj=1, and r is the number of the retained POD modes. In practice application, the truncated POD modes with energies that sum to 99% of the total energy are a plausible truncation.With the determination of the relative important POD modes, the flow field can be represented in a reduced dimensionality:(8)xj′≈Φaj,j=1,2,3,…,r.For practical applications, the decomposition can be carried out using the classical method or the snapshot method. In matrix form, the POD can be realized by the standard singular value decomposition (SVD) which is a robust approach to determine the POD modes [30]. In this paper, we use the SVD to implement the POD decomposition; the specific details are introduced below.


## 2.2. Computation of the Proper Orthogonal Decomposition

### 2.2.1. Preprocessing

In the field of atmospheric/marine science, whether the forecast field data generated by the high-resolution numerical forecasting model or the analysis field produced by the global variational data assimilation system, they are distributed in the regular space grid at regular intervals. And they only describe the state of atmospheric motion at a certain moment, that is, strictly speaking, an instantaneous state field (actually, the state changes little before and after the moment) which is equivalent to store using the camera after taking pictures of the atmospheric motion state. Therefore, it can be named as the atmospheric motion snapshot analogy to the concept of flow snapshot in fluid mechanics. Thus, in this paper, we use the modal decomposition methods to analyze the snapshots of atmospheric motion over some time.Here, we collect data from a two-dimensional stratospheric jet flow velocity fieldu,v and construct snapshot matrix X in preparation to perform the modal decomposition analysis. The matrix X is constructed by the velocity components u,v of each measurement point of the M instantaneous stratospheric jet flow fields; each column of X contains all the velocity data at each time instance; hence, the matrix X is defined as follows:(9)X=u11u12⋯u1N⋮⋮⋮⋮uM1uM2⋯uMNv11v12⋯v1N⋮⋮⋮⋮vM1vM2⋯vMN,where the apex 1,2,3,…,N indicates the time number and the term 1,2,3,…,M is the number of measurement points.


### 2.2.2. Implementation of POD

The key idea of the POD is to reduce the data dimensions by projecting high-dimensional data to a low-dimensional space while reserving the variation in the original variables as much as possible. It can be seen as an application of singular value decomposition (SVD) in a finite-dimensional space, and a brief review of the algorithm is shown as follows.(1)
Consider the given snapshot matrixX∈ℝM×N. Firstly, compute the mean flow X¯=1/N∑i=1Nxi and the mean-subtracted snapshot matrix X′=X−X¯.(2)
The snapshots matrixX′ can be decomposed with SVD as(10)X′=ΦΣΨ∗,whereΦ∈ℝM×N, Σ∈ℝN×NandΨ∈ℝN×N, matrix Σ holds the singular values σjj=1,2,…,N along its diagonal, that define how much “energy” contained in this mode, and the asterisk denotes the conjugate transpose. Besides, the left and right singular vectors of X are the eigenvectors of XXT and XTX, respectively. Moreover, the eigenvalues of XXTXTX and singular values of X are related by λj=σj2. Thus, the SVD can be directly performed on X to determine the POD modes Φ and corresponding eigenvalue and provide time-evolving information of the flow field contained in the matrix Ψ.(3)
In practical applications, one usually retainsrPOD several modes that capture 99% of the snapshots energy to achieve a high quality of model reduction, such that(11)EPOD=∑j=1rPODσj2∑j=1Nσj2≈99%.(4)
With the determination of the most important POD modes, the jet flow field can be represented only using the truncated modes, as follows:(12)X≈X¯+∑j=1rPODDjσj.To obtain the matrixDj, we multiply the corresponding columns of Φ and Ψ with an outer product. In the MATLAB notation, the computation of Dj is the following: Dj=Φ:,jV:,jT.


## 2.2.1. Preprocessing

In the field of atmospheric/marine science, whether the forecast field data generated by the high-resolution numerical forecasting model or the analysis field produced by the global variational data assimilation system, they are distributed in the regular space grid at regular intervals. And they only describe the state of atmospheric motion at a certain moment, that is, strictly speaking, an instantaneous state field (actually, the state changes little before and after the moment) which is equivalent to store using the camera after taking pictures of the atmospheric motion state. Therefore, it can be named as the atmospheric motion snapshot analogy to the concept of flow snapshot in fluid mechanics. Thus, in this paper, we use the modal decomposition methods to analyze the snapshots of atmospheric motion over some time.Here, we collect data from a two-dimensional stratospheric jet flow velocity fieldu,v and construct snapshot matrix X in preparation to perform the modal decomposition analysis. The matrix X is constructed by the velocity components u,v of each measurement point of the M instantaneous stratospheric jet flow fields; each column of X contains all the velocity data at each time instance; hence, the matrix X is defined as follows:(9)X=u11u12⋯u1N⋮⋮⋮⋮uM1uM2⋯uMNv11v12⋯v1N⋮⋮⋮⋮vM1vM2⋯vMN,where the apex 1,2,3,…,N indicates the time number and the term 1,2,3,…,M is the number of measurement points.


## 2.2.2. Implementation of POD

The key idea of the POD is to reduce the data dimensions by projecting high-dimensional data to a low-dimensional space while reserving the variation in the original variables as much as possible. It can be seen as an application of singular value decomposition (SVD) in a finite-dimensional space, and a brief review of the algorithm is shown as follows.(1)
Consider the given snapshot matrixX∈ℝM×N. Firstly, compute the mean flow X¯=1/N∑i=1Nxi and the mean-subtracted snapshot matrix X′=X−X¯.(2)
The snapshots matrixX′ can be decomposed with SVD as(10)X′=ΦΣΨ∗,whereΦ∈ℝM×N, Σ∈ℝN×NandΨ∈ℝN×N, matrix Σ holds the singular values σjj=1,2,…,N along its diagonal, that define how much “energy” contained in this mode, and the asterisk denotes the conjugate transpose. Besides, the left and right singular vectors of X are the eigenvectors of XXT and XTX, respectively. Moreover, the eigenvalues of XXTXTX and singular values of X are related by λj=σj2. Thus, the SVD can be directly performed on X to determine the POD modes Φ and corresponding eigenvalue and provide time-evolving information of the flow field contained in the matrix Ψ.(3)
In practical applications, one usually retainsrPOD several modes that capture 99% of the snapshots energy to achieve a high quality of model reduction, such that(11)EPOD=∑j=1rPODσj2∑j=1Nσj2≈99%.(4)
With the determination of the most important POD modes, the jet flow field can be represented only using the truncated modes, as follows:(12)X≈X¯+∑j=1rPODDjσj.To obtain the matrixDj, we multiply the corresponding columns of Φ and Ψ with an outer product. In the MATLAB notation, the computation of Dj is the following: Dj=Φ:,jV:,jT.


## 2.3. Sparsity-Promoting DMD

In this section, we first give a brief overview of the Koopman operator theory, which is the basis for DMD analysis, and then describe the basic procedure for the standard DMD and sparsity-promoting DMD method; the latter one leads to an optimal selection of the original DMD modes.First, we also consider the discrete-time dynamical system mentioned in Section2.1, in equation (2) xi+1=Fxi.The Koopman operatorK proposed by Koopman [18] in 1931 is an infinite-dimensional linear operator that operates on any scalar-valued function g in the following manner:(13)Kgx=gFx.The Koopman operator steps forward in time an observable. Note thatK is a linear operator; it is supposed that the operator K can be decomposed into eigenvalues μj and corresponding eigenfunctions ϕj, which satisfy(14)Kϕj=μjϕj,j=1,2,3,….According to Rowley et al. [17], we consider a vector-valued function g that lies within the span of eigenfunctions of Koopman operator, and g can be expanded in terms of these eigenfunctions as(15)gx=∑j=1∞ϕjxυj.Thus, the functiongx was expanded as a linear combination of the eigenfunctions ϕjof K, where υj are the vector coefficients which are referred to as Koopman modes of the map F, corresponding to the observable g.Under the action of Koopman operator, the evolution of vector functiongx can be expressed as follows:(16)gxi=Kigx0=∑j=1∞μjiϕjx0υj.Thus, the eigenvalueμj characterizes the temporal behavior of the Koopman modes υj, while the phase of μj determines its frequency, and the magnitude determines the growth rate. These are very important information for the stability analyses of the modes, as when the growth rate is positive, the corresponding mode is unstable, while when the growth rate is negative, the corresponding mode is stable, and the mode is neutrally stable when the growth rate is just zero. Note that, as described in [16], for a system evolving on an attractor, the Koopman eigenvalues always lie on the unit circle.In DMD architecture, consider finite measurements at timesti of the dynamic system in equation (2):(17)yi=gxi,i=1,2,3,…,N.In practice application, the state can be directly measured; thus,yi=xi and the DMD method [15] assumes that there is a locally best-fit linear operator A to step forward in time the snapshots(18)xi+1=Axi.It is easy to prove that the eigenvalues of the operatorA are also ones of the Koopman operator, and the Koopman modes coincide with the eigenvectors of A; thus, DMD method can be used for estimating an approximation of the decomposition in equation (16).Arrange theN snapshots into two data matrices, as follows:(19a)X1N−1=x1x2x3…xN−1,(19b)X2N=x2x3…xN.The locally linear approximation in equation (18) may be written in terms of these data matrices as(20)X2N≈AX1N−1.The matrixA can be given as(21)A=X2NX1N†,where † represents the Moore–Penrose pseudoinverse. As already highlighted, the dimensions of the snapshots xi are large m≫1; thus, the matrix A is difficult to represent or analyze directly.For this reason, DMD circumvents the eigendecomposition ofA by considering a rank-reduced representation in terms of a POD-projected matrixA˜ based on singular value decomposition (SVD) of X1N−1:(22)X1N−1=UΣV∗,where U∈CM×r, Σ∈Cr×r, and V∈CN×r. The left singular vectors U are the POD modes, matrix Σ is a diagonal matrix with the singular values σ1,σ2,…,σr that define how much “energy” contained in this mode，and r is the rank of X1N−1.Thus, the matrixA can be obtained as(23)A=X2NVΣ−1U∗.By projecting the high-dimensional matrixA onto POD, modes can obtain the representation of A˜ as(24)A˜=U∗AU=U∗X2NVΣ−1.The matrixA˜ defines a low-dimensional linear model of the dynamical system on POD coordinates, in which(25)x˜i+1=A˜x˜i,and the high-dimensional state xi can be reconstructed as xi=Ux˜i.The matrixA˜ is a similar transformation of A, so that the matrix A˜ contains the main eigenvalues of A; compute the eigendecomposition of A˜ as(26)A˜W=WΛ,where columns of W are eigenvectors and Λ is a diagonal matrix containing the corresponding eigenvalues μj.Therefore, the snapshot at a timeti can be estimated as(27)xi=UWΛi−1W−1U∗x1=∑j=1rμji−1αjυj,in which υj=UW are the DMD modes and αj is the amplitude of the jth DMD mode that represents the jth modal contribution of the initial condition defined as α=W−1x˜1=W−1U∗x1.Rewrite equation (27) in matrix form as(28)x1x2…xN−1≈υ1υ2…υrα10⋯00α2⋯0⋮⋮⋱⋮00⋯αr︸Dα=diagα1μ1⋯μ1N−11μ2⋯μ2N−1⋮⋮⋱⋮1μr⋯μrN−1︸Vand,in which the temporal evolution information of the DMD modes is contained in the Vandermonde matrix Vand.As there are no explicit ways for DMD to select modes, unlike the POD modes which are selected based on the rank of energy from the singular values. The sparsity-promoting DMD method is proposed by Jovanovic et al. which induces a sparse structure of the DMD modes by augmenting thel1-norm of the vector α through the above equation. That is, it is expressed as the following optimization problem:(29)minαX1N−1−υDαVandF2+γ∑jrαj,in which γ is a positive regularization parameter that reflects the emphasis on the sparsity of the amplitudes vector α, larger values of γ emphasize on the number of nonzero elements in the vector α, and ⋅F2 is the Frobenius norm. And the values of γ are selected in such a way that γmin induces a dense vector α (with 191 nonzero elements), and γmax induces α with a single nonzero element.After the sparse structure ofα is fixed, the problem becomes to determine only the nonzero amplitudes as the solution from the optimization problem in equation (29) with the constraint ETα=0, where the matrix E contains the sparsity structure of α. As α=0α20α4T, the matrix E is given by E=10000010T.A powerful and well-suited algorithm named Alternating Direction Method of Multipliers (ADMM) for solving the optimization problems of equation (29) can be found in [27].


## 2.4. Computation of the Sparsity-Promoting DMD

The sparsity-promoting DMD is a variant of the standard DMD method; here, the main steps of the SPDMD algorithm are summarized as follows:(1)
First, a matrixX1N−1 is formed with the first N−1 columns and matrix X2N contains the last N−1 columns of the snapshot matrix X shown in equation (9) And then perform the standard DMD for the matrices to obtain DMD solution A˜ and eigenvector matrix W and corresponding eigenvalue matrix Λ (and also U, V, and Σ).(2)
The SPDMD aims to solve the optimization problem in equation (29). We apply the alternating direction method of multipliers (ADMM) to solve this convex problem. ADMM is a state-of-the-art algorithm for solving large-scale and distributed optimization problems [31]. And there are four steps for the application of the ADMM.Step 1. Replace the vector of amplitude α in the l1-regularization with a new variable β; the optimization problem in equation (29) is transformed into the following constrained optimization problem:(30)minimizeα,βJα+γhβsubject toα−β=0,whereJα=minαX1N−1−υDαVandF2, and hβ=∑jrβj1. β is an auxiliary variable; for any feasible α and β, the optimization problems in equations (29) and (30) are equivalent.Step 2. To solve the problem, the augmented Lagrangian is introduced to transform the constraints in equation (30) to the objective function, as follows:(31)Lρα,β,λ=Jα+γhβ+12θ∗α−β+α−β∗θ+ρα−β22,hereθ is the vector of Lagrange multipliers and ρ is a positive parameter that introduces a quadratic penalty on the deviation between α and β, and ⋅2 is the Euclidean norm of a given vector.Step 3. ADMM algorithm is an iterative algorithm that combines the Dual Ascent and Multiplier method. It includes a α-minimization step, a β-minimization step, and a Lagrange multiplier update step; the details are shown as follows:(32a)αk+1=argminαLρα,βk,θk,(32b)βk+1=argminβLραk+1,β,θk,(32c)θk+1=λk+ραk+1−βk+1,in which theα-minimization step amounts to solving an unconstrained regularized quadratic program and that the β-minimization step amounts to an opportune use of the soft thresholding operator.α-Minimization step: completion of squares concerning α in the augmented Lagrangian Lρ can be used to show that the α-minimization step (33c) is equivalent to(33)minαJα+ρ2α−u2k2,uk=βk−1ρθk,in whichk is the step number of iterations. The solution is determined [27] by(34)αk+1=P+ρ2I−1q+ρ2uk,in whichP=W∗W∘VandV∗¯, q=diagVandVΣ∗W¯. And I is the unit matrix, an overline signifies the complex-conjugate of a vector or matrix, and ∘ is the elementwise multiplication of two matrices in this case [32].β-Minimization step: the completion of squares for β in the augmented Lagrangian Lρ can be used to show that the β-minimization step (32b) is equivalent to(35)minβγhβ+ρ2β−v2k2,vk=αk+1+1ρθk,and the solution is determined by(36)βjk+1=Sκvjk,κ=γρ,whereSκ⋅ denotes the soft thresholding operator:(37)Sκvjk=vjk−κ,vjk>κ,0,vjk∈−κ,κ,vjk−κ,vjk>−κ.Step 4. Start with an initial point α0,β0 and iterate until the required error is εprim and εdual then stop:(38)αk+1−β2k+1≤εprim,βk+1−β2k≤εdual.


## 3. Results and Discussion

The data used in this paper is the ERA5 daily reanalysis dataset, the horizontal resolution of the data is 0.25° × 0.25°, and the time resolution is 1 hour [33, 34]. The SPDMD method is firstly applied to decompose a stratospheric 200 hPa jet dataset that includes 192 snapshots from 00:00 on September 24, 2009, to 24:00 on the 31st, and the analysis area is in the pan-Asian region (the range of latitude and longitude: 20°∼70°N, 40°∼130°E). Simultaneously, the decomposed results of the SPDMD method are compared with the ones from the POD method.


### 3.1. Analysis of the Stratospheric 200 hPa Jet

Before applying the POD and SPDMD to the reanalysis dataset, it seems instructive to treat the raw data to make ourselves familiar with the features of its output modes. For this reason, four snapshots of the original stratospheric 200 hPa jet are depicted to give a first impression of the flow field. In Figure1, the arrows show the vector field and the background color contours visualize the vertical vorticity component. As Figure 1 depicted, the spatial distribution of the flow field at four different times all present a fluctuating flow structure, dominated by zonal western wind, and this feature extends uniformly in the east-west direction along 20°～40°N. At high latitude, there exist some large-scale vortex structures with different shapes at different times. In this paper, the SPDMD method is introduced into the analysis of atmospheric jet aiming to improve the understanding of the dynamic characteristics buried in the complex atmospheric motion.Figure 1
Four snapshots of the stratospheric 200 hPa jet at time. (a)t = 1. (b) t = 24. (c) t = 48. (d) t = 60.
(a)(b)(c)(d)The POD and SPDMD methods are applied to decompose the snapshots of the stratospheric 200 hPa jet. Figure2(a) shows the percentage of energy in each POD and SPDMD mode, and the corresponding cumulative energy distribution is shown in Figure 2(b). As shown in Figure 2(a), most of the velocity fluctuation energy in POD is contained in the first few modes and shows a rather weak decay after the first few modes initial drop. As the first POD mode contains the most velocity fluctuation energy, reaching 44.15%, and the fluctuation energies of the second and third modes are rapidly reduced to 18.67% and 14.18%, respectively, and the energies of the fourth and fifth modes are both around 5.0%. In SPDMD, the energy is more evenly distributed [35], and as the number of modes increases, the energy of the mode decreases rate is slower than that in POD. The first two modes of energy are about 9.34% and 4.38%, respectively, just a little more than the other modes. As the number of modes increases, the energy of the mode decreases gradually and then remains at a low energy level. Correspondingly, as depicted in Figure 2(b), the convergence speed of the percent of the cumulative energy considerably slows down as more modes are included. In practice, to construct a reduced-order model for the given flow field, POD reserves around 26 number of modes that capture nearly 99% of the total energy, while nearly 151 number of modes in SPDMD are needed to reconstruct the same energy. Thus, considering the energy content across the modes, POD is preferable. Nevertheless, POD modes ranked in the order of energy content and not in the order of the dynamical importance as SPDMD; this may ignore some important modes that are dynamically relevant but contain low energy [22].Figure 2
(a) Percentage of the energy of each mode (b) and cumulative energy distribution of the 200 hPa jet.
(a)(b)The first eight POD spatial modes of the 200 hPa jet are shown in Figure3 (mode is normalized). As can be seen in Figure 3, the arrows show the vector field, and the background color contours visualize the vertical vorticity component. As Figure 3(a) depicted, the spatial distribution of the first POD mode is dominated by the large cyclone and anticyclonic circulation along the longitude direction. It is the most important coherent structure underlying the jet flow captured by the POD method, as the first POD mode contains the most fluctuation energy of the jet flow. The second and third POD modes have similar magnitude velocity fluctuation energy, and there is a strong resemblance of the spatial distribution between the two modes, with the nearly same sizes and locations of the vortices’ structures. The fifth and sixth modes that capture only about 8% of the total fluctuation energy highlight small-scale structures in space compared with the first four POD modes, and the spatial structures of the seventh and eighth modes are dominated by more vortices with smaller-scale structures as they have fewer fluctuation energy. The POD method can decompose the 200 hPa jet flow into different modes with different spatial scales; the first few modes capture the most velocity fluctuation energy with large-scale features, while the following modes with low-level energy content are dominated by small-scale characteristics. Different modes represent the motion components of different spatial scales of the jet flow, and the superposition of these modes can represent the main characteristics of the jet flow. It is worth noting that there is no strict correspondence between the spatial scales of the POD modes and the large-, medium-, and small-scales in atmospheric science; the latter corresponds to specific weather systems of different spatial scales, while the different modes represent the components of the flow field at different spatial scales.Figure 3
The first eight POD modes extracted from the 200 hPa jet (a–h correspond to the 1–8th modes; wind velocity visualized by vector plots and the background color contours visualize the vertical vorticity component).
(a)(b)(c)(d)(e)(f)(g)(h)Figure4(a) shows the frequency dependence of the absolute value of the amplitudes of DMD modes. As shown in Figure 4(a), a pronounced peak at the frequency is zero can be observed and the absolute value of amplitudes decreases as the increase of their respective frequencies on the whole. It is hard to identify subset modes that have an important impact on the dynamic of the jet flow field just based on the amplitudes of the DMD modes. As shown in Figure 4(b), most of the modes concentrated near-zero growth rate indicating the corresponding modes are stable, while there exist modes with large amplitudes originating from eigenvalues with strongly damped growth rates.Figure 4
The absolute value of the DMD amplitudesαi distribution versus (a) the frequency and (b) the real part of the corresponding DMD eigenvalues μi for the 200 hPa jet. (c) A zoomed version of (b).
(a)(b)(c)In SPDMD, as our emphasis on sparsity increases, a smaller number of modes with nonzero amplitudes are obtained, and the fidelity of the reduced-order model decreases, as quantified in Figure5 using the performance loss norm defined in [27] as follows:(39)Eloss=100×X1N−1−υDαVandFX1FN−1.Figure 5
Performance loss with the number of retained DMD modes resulting from SPDMD for the 200 hPa jet.As illustrated in Figure5, the SPDMD algorithm succeeds in isolating few modes, which have a significant influence on one performance loss as the level of sparsity increases, and the performance loss quickly decays with the level of sparsity. With nearly 11 DMD modes, the SPDMD algorithm can achieve a reasonable compromise between the quality of approximation and the size of the reduced-order model for the 200 hPa jet snapshots which introduces performance degradation just about 10.2%.Figure6(a) shows the full eigenvalues resulting from the 200 hPa jet flow field snapshots in which the subset eigenvalues selected by the SPDMD algorithm marked with the red-filled circles. The magnitude of μj relative to the unit circle indicates the stability of the corresponding DMD spatial mode υj, where modes whose eigenvalues lie exactly on the unit circle are neutrally stable, in the unit circle are stable, falling outside the unit circle are unstable [36]. As shown in Figure 6(a), most of the eigenvalues tend to distribute on the vicinity of the unit circle in the complex pane of μi=1 indicating that the corresponding modes are stable; few eigenvalues locate outside the unit circle corresponding to unstable modes. These stable modes signify a corresponding weather process that is steady or gradually decays in time; in contrast, the unstable modes with positive growth rates correspond to the structure of atmospheric flow field that will gradually strengthen in the future, which play a leading role in the development of weather systems, while their specific impact needs further research. The eigenvalues falling within the unit circle are corresponding to strongly damped modes that may be associated with large amplitudes as shown in Figure 4(b); however, these modes only influence early stages’ evolution of the jet flow, while the SPDMD algorithm focuses on a subset of DMD modes that have a profound influence on the full time of the jet flow. Figure 6(b) shows the magnified view of the relevant part of the unit pane, in which the numbers 1–5 are five modes (their complex conjugates are not shown) selected by the SPDMD algorithm that has the most important impact for the underlying dynamics of the 200 hPa jet flow data sequence. The first mode with zero frequency signifies the mean flow; other modes with low frequencies are near the unit circle.Figure 6
Distribution of eigenvalues of 200 hPa jet: open circles represent full DMD eigenvalues; filled circles that are numbered with 1–5 are the modes isolated by the SPDMD: (a) full spectrum, (b) zoom on the isolated modes.
(a)(b)We make a change of units by taking the logarithm ofμj and normalize it by Δt, ωj=logμj/2πΔt so that now the imaginary component ofωj is the frequency of oscillation in units of cycles per second (Hz) [36]. Table 1 lists the growth rates and frequencies of the five modes shown in Figure 6(b). The first mode signifies the main structure of the jet flow with zero frequency and a small growth rate. The modes 2 to 5 are conjugate, and the second mode is unstable with a positive growth rate, which plays a leading role in the development of the flow field. The rest of the modes are stable modes, with gradually increasing frequencies as the number increased.Table 1
Growth rate and frequency of the five SPDMD modes of the 200 hPa jet flow field.
DMD modeGrowth rateFrequency18.4299e − 07022.3828e − 074.1590e − 063−2.0152e − 061.7497e − 054−1.3161e − 063.2941e − 055−1.6680e − 065.3176e − 05The real part of the DMD modes identified in Figure6(b) is shown in Figure 7 (mode is normalized), in which the arrows show the vector field and the background color contours visualize the vertical vorticity component, and the imaginary part (not shown) is shifted 90°. As shown in Figure 7, the five DMD modes reveal the multiscale structure features of the 200 hPa jet flow. The first DMD mode with zero frequency shows the structure reminiscent of the mean flow or the temporally averaged flow field as depicted in Figure 7(a), the spatial structure of the first DMD mode is dominated by zonal western wind, and there exists large-scale vortex structure at high latitude. The spatial structure distribution of the second DMD mode characterized by localized vortices in the background of the zonal western wind, and the second mode is unstable; thus, the structure of this mode will strengthen in the future, which has an important impact on the evolution of the jet flow. It would be worth noticing that the direction of rotation of the vortex in the model does not represent that in the real flow field, as the spatial structure of the model is only the motion law extracted by the DMD method from the 200 hPa jet flow [37]. The spatial structure distribution of the third DMD mode is very similar to the fourth POD mode (see Figure 3(d)); there are some common features between the two graphs at nearly the same locations, and the same relationship also exists between the fourth DMD and the eighth POD mode, which indicates that the two decomposition methods can extract the dominant coherent structures from the jet flow, while the DMD method has clear advantages over the POD, as DMD separates the dominant modes with a single frequency and the SPDMD method selects modes based on their dynamical importance for the whole time evolution of the jet flow while the POD modes are arranged just by their energy content. In the atmospheric/marine science, it is easier to interpret the physical meaning contained in each mode which has a single frequency, and researchers focus more on the dynamical information underlying the flow field to identify the most important features that have an important impact on the development of the atmospheric/marine motion process. As shown in Figure 7(c), the third DMD mode has smaller vortex structures than the first two DMD modes which indicates a motion with higher frequency. The fifth mode spatial distribution is messier, and the scale of the vortex is increasingly smaller consistent with its higher associated frequency as depicted in Figure 7(d). Therefore, applying the SPDMD method to atmospheric/marine motion analysis is helpful to further understand the multiscale characteristics of atmospheric motion. The SPDMD method can decompose the 200 hPa jet flow into different modes associated with a particular frequency and growth rate which can be used to identify the stability of the model. The low-frequency modes represent the large-scale motion structures, while the high-frequency modes highlight the small-scale features. Similar to the POD method, the multiscale features captured by the different DMD modes represent different spatial scales of the motion components for the jet flow, and a linear combination of these modes can represent the main characteristics of the jet flow [18]. Compared with the POD spatial modes depicted in Figure 3, it can be seen that the SPDMD method selects a smaller number of modes that capture nearly the same dynamic structures of the 200 hPa jet flow as the POD method.Figure 7
The five DMD modes identified by the SPDMD from the 200 hPa jet flow (a–e correspond to the 1–5th modes, velocity is visualized by vector plots, and the background color contours visualize the vertical vorticity component).
(a)(b)(c)(d)(e)


### 3.2. Reconstruction of the 200 hPa Jet

To further demonstrate the extraction effect of POD/SPDMD methods on the 200 hPa jet flow characteristics, we consider constructing the POD/SPDMD reduced-order models, respectively. In POD, to reconstruct 99% of the total fluctuation energy, the POD modes are truncated atrPOD=26. While in SPDMD, there are only ten modes (modes shown in Figure 6, DM1, and DM2,3,4,5 with their complex conjugates) attempted. The initial transient flow field at the sampling period T=60h,124h is depicted in Figure 8.Figure 8
Initial velocity field of 200 hPa jet flow at (a)T = 60 h, (b) T = 124 h (velocity visualized by vector plots and the background color contours visualize the vorticity component).
(a)(b)The POD/SPDMD reconstruction flow fields atT=60h,124h are presented in Figures 9 and 10, respectively. The highly similar features between the graphics show that the POD/SPDMD reduced-order models successfully recover the original dynamic process to an impressive degree, which indicates that both approaches are effective in detecting and extracting the relevant dynamic features underlying the 200 hPa jet flow. To quantitatively evaluate the quality of the POD/SPDMD reduced-order models develop by equations (12) and (27), the average relative error is defined as follows:(40)error=1N∑i=1NXi−XROMFiXFi,where Xi represent a vector of the original snapshots matrix at the time i and XROMi is an approximation from the POD/SPDMD reduced-order models at the same time; the results are listed in Table 2.Figure 9
POD reconstruction of the 200 hPa jet flow field (26 POD modes) at (a)T = 60 h and (b) T = 124 h (velocity visualized by vector plots and the background color contours visualize the vorticity component).
(a)(b)Figure 10
SPDMD reconstruction of the 200 hPa jet flow field (11 SPDMD modes) at (a)T = 60 h and (b) T = 124 h (velocity visualized by vector plots and the background color contours visualize the vorticity component).
(a)(b)Table 2
The average relative error of the POD/SPDMD reduced-order model.
Number of modesError ofROMDMDError ofROMPOD11 modes0.09710.110017 modes0.06710.084325 modes0.05930.064047 modes0.04430.0383As shown in Table2, compared with the POD method, the results of the SPDMD method show that it has a smaller average relative error when the number of modes is small.


## 3.1. Analysis of the Stratospheric 200 hPa Jet

Before applying the POD and SPDMD to the reanalysis dataset, it seems instructive to treat the raw data to make ourselves familiar with the features of its output modes. For this reason, four snapshots of the original stratospheric 200 hPa jet are depicted to give a first impression of the flow field. In Figure1, the arrows show the vector field and the background color contours visualize the vertical vorticity component. As Figure 1 depicted, the spatial distribution of the flow field at four different times all present a fluctuating flow structure, dominated by zonal western wind, and this feature extends uniformly in the east-west direction along 20°～40°N. At high latitude, there exist some large-scale vortex structures with different shapes at different times. In this paper, the SPDMD method is introduced into the analysis of atmospheric jet aiming to improve the understanding of the dynamic characteristics buried in the complex atmospheric motion.Figure 1
Four snapshots of the stratospheric 200 hPa jet at time. (a)t = 1. (b) t = 24. (c) t = 48. (d) t = 60.
(a)(b)(c)(d)The POD and SPDMD methods are applied to decompose the snapshots of the stratospheric 200 hPa jet. Figure2(a) shows the percentage of energy in each POD and SPDMD mode, and the corresponding cumulative energy distribution is shown in Figure 2(b). As shown in Figure 2(a), most of the velocity fluctuation energy in POD is contained in the first few modes and shows a rather weak decay after the first few modes initial drop. As the first POD mode contains the most velocity fluctuation energy, reaching 44.15%, and the fluctuation energies of the second and third modes are rapidly reduced to 18.67% and 14.18%, respectively, and the energies of the fourth and fifth modes are both around 5.0%. In SPDMD, the energy is more evenly distributed [35], and as the number of modes increases, the energy of the mode decreases rate is slower than that in POD. The first two modes of energy are about 9.34% and 4.38%, respectively, just a little more than the other modes. As the number of modes increases, the energy of the mode decreases gradually and then remains at a low energy level. Correspondingly, as depicted in Figure 2(b), the convergence speed of the percent of the cumulative energy considerably slows down as more modes are included. In practice, to construct a reduced-order model for the given flow field, POD reserves around 26 number of modes that capture nearly 99% of the total energy, while nearly 151 number of modes in SPDMD are needed to reconstruct the same energy. Thus, considering the energy content across the modes, POD is preferable. Nevertheless, POD modes ranked in the order of energy content and not in the order of the dynamical importance as SPDMD; this may ignore some important modes that are dynamically relevant but contain low energy [22].Figure 2
(a) Percentage of the energy of each mode (b) and cumulative energy distribution of the 200 hPa jet.
(a)(b)The first eight POD spatial modes of the 200 hPa jet are shown in Figure3 (mode is normalized). As can be seen in Figure 3, the arrows show the vector field, and the background color contours visualize the vertical vorticity component. As Figure 3(a) depicted, the spatial distribution of the first POD mode is dominated by the large cyclone and anticyclonic circulation along the longitude direction. It is the most important coherent structure underlying the jet flow captured by the POD method, as the first POD mode contains the most fluctuation energy of the jet flow. The second and third POD modes have similar magnitude velocity fluctuation energy, and there is a strong resemblance of the spatial distribution between the two modes, with the nearly same sizes and locations of the vortices’ structures. The fifth and sixth modes that capture only about 8% of the total fluctuation energy highlight small-scale structures in space compared with the first four POD modes, and the spatial structures of the seventh and eighth modes are dominated by more vortices with smaller-scale structures as they have fewer fluctuation energy. The POD method can decompose the 200 hPa jet flow into different modes with different spatial scales; the first few modes capture the most velocity fluctuation energy with large-scale features, while the following modes with low-level energy content are dominated by small-scale characteristics. Different modes represent the motion components of different spatial scales of the jet flow, and the superposition of these modes can represent the main characteristics of the jet flow. It is worth noting that there is no strict correspondence between the spatial scales of the POD modes and the large-, medium-, and small-scales in atmospheric science; the latter corresponds to specific weather systems of different spatial scales, while the different modes represent the components of the flow field at different spatial scales.Figure 3
The first eight POD modes extracted from the 200 hPa jet (a–h correspond to the 1–8th modes; wind velocity visualized by vector plots and the background color contours visualize the vertical vorticity component).
(a)(b)(c)(d)(e)(f)(g)(h)Figure4(a) shows the frequency dependence of the absolute value of the amplitudes of DMD modes. As shown in Figure 4(a), a pronounced peak at the frequency is zero can be observed and the absolute value of amplitudes decreases as the increase of their respective frequencies on the whole. It is hard to identify subset modes that have an important impact on the dynamic of the jet flow field just based on the amplitudes of the DMD modes. As shown in Figure 4(b), most of the modes concentrated near-zero growth rate indicating the corresponding modes are stable, while there exist modes with large amplitudes originating from eigenvalues with strongly damped growth rates.Figure 4
The absolute value of the DMD amplitudesαi distribution versus (a) the frequency and (b) the real part of the corresponding DMD eigenvalues μi for the 200 hPa jet. (c) A zoomed version of (b).
(a)(b)(c)In SPDMD, as our emphasis on sparsity increases, a smaller number of modes with nonzero amplitudes are obtained, and the fidelity of the reduced-order model decreases, as quantified in Figure5 using the performance loss norm defined in [27] as follows:(39)Eloss=100×X1N−1−υDαVandFX1FN−1.Figure 5
Performance loss with the number of retained DMD modes resulting from SPDMD for the 200 hPa jet.As illustrated in Figure5, the SPDMD algorithm succeeds in isolating few modes, which have a significant influence on one performance loss as the level of sparsity increases, and the performance loss quickly decays with the level of sparsity. With nearly 11 DMD modes, the SPDMD algorithm can achieve a reasonable compromise between the quality of approximation and the size of the reduced-order model for the 200 hPa jet snapshots which introduces performance degradation just about 10.2%.Figure6(a) shows the full eigenvalues resulting from the 200 hPa jet flow field snapshots in which the subset eigenvalues selected by the SPDMD algorithm marked with the red-filled circles. The magnitude of μj relative to the unit circle indicates the stability of the corresponding DMD spatial mode υj, where modes whose eigenvalues lie exactly on the unit circle are neutrally stable, in the unit circle are stable, falling outside the unit circle are unstable [36]. As shown in Figure 6(a), most of the eigenvalues tend to distribute on the vicinity of the unit circle in the complex pane of μi=1 indicating that the corresponding modes are stable; few eigenvalues locate outside the unit circle corresponding to unstable modes. These stable modes signify a corresponding weather process that is steady or gradually decays in time; in contrast, the unstable modes with positive growth rates correspond to the structure of atmospheric flow field that will gradually strengthen in the future, which play a leading role in the development of weather systems, while their specific impact needs further research. The eigenvalues falling within the unit circle are corresponding to strongly damped modes that may be associated with large amplitudes as shown in Figure 4(b); however, these modes only influence early stages’ evolution of the jet flow, while the SPDMD algorithm focuses on a subset of DMD modes that have a profound influence on the full time of the jet flow. Figure 6(b) shows the magnified view of the relevant part of the unit pane, in which the numbers 1–5 are five modes (their complex conjugates are not shown) selected by the SPDMD algorithm that has the most important impact for the underlying dynamics of the 200 hPa jet flow data sequence. The first mode with zero frequency signifies the mean flow; other modes with low frequencies are near the unit circle.Figure 6
Distribution of eigenvalues of 200 hPa jet: open circles represent full DMD eigenvalues; filled circles that are numbered with 1–5 are the modes isolated by the SPDMD: (a) full spectrum, (b) zoom on the isolated modes.
(a)(b)We make a change of units by taking the logarithm ofμj and normalize it by Δt, ωj=logμj/2πΔt so that now the imaginary component ofωj is the frequency of oscillation in units of cycles per second (Hz) [36]. Table 1 lists the growth rates and frequencies of the five modes shown in Figure 6(b). The first mode signifies the main structure of the jet flow with zero frequency and a small growth rate. The modes 2 to 5 are conjugate, and the second mode is unstable with a positive growth rate, which plays a leading role in the development of the flow field. The rest of the modes are stable modes, with gradually increasing frequencies as the number increased.Table 1
Growth rate and frequency of the five SPDMD modes of the 200 hPa jet flow field.
DMD modeGrowth rateFrequency18.4299e − 07022.3828e − 074.1590e − 063−2.0152e − 061.7497e − 054−1.3161e − 063.2941e − 055−1.6680e − 065.3176e − 05The real part of the DMD modes identified in Figure6(b) is shown in Figure 7 (mode is normalized), in which the arrows show the vector field and the background color contours visualize the vertical vorticity component, and the imaginary part (not shown) is shifted 90°. As shown in Figure 7, the five DMD modes reveal the multiscale structure features of the 200 hPa jet flow. The first DMD mode with zero frequency shows the structure reminiscent of the mean flow or the temporally averaged flow field as depicted in Figure 7(a), the spatial structure of the first DMD mode is dominated by zonal western wind, and there exists large-scale vortex structure at high latitude. The spatial structure distribution of the second DMD mode characterized by localized vortices in the background of the zonal western wind, and the second mode is unstable; thus, the structure of this mode will strengthen in the future, which has an important impact on the evolution of the jet flow. It would be worth noticing that the direction of rotation of the vortex in the model does not represent that in the real flow field, as the spatial structure of the model is only the motion law extracted by the DMD method from the 200 hPa jet flow [37]. The spatial structure distribution of the third DMD mode is very similar to the fourth POD mode (see Figure 3(d)); there are some common features between the two graphs at nearly the same locations, and the same relationship also exists between the fourth DMD and the eighth POD mode, which indicates that the two decomposition methods can extract the dominant coherent structures from the jet flow, while the DMD method has clear advantages over the POD, as DMD separates the dominant modes with a single frequency and the SPDMD method selects modes based on their dynamical importance for the whole time evolution of the jet flow while the POD modes are arranged just by their energy content. In the atmospheric/marine science, it is easier to interpret the physical meaning contained in each mode which has a single frequency, and researchers focus more on the dynamical information underlying the flow field to identify the most important features that have an important impact on the development of the atmospheric/marine motion process. As shown in Figure 7(c), the third DMD mode has smaller vortex structures than the first two DMD modes which indicates a motion with higher frequency. The fifth mode spatial distribution is messier, and the scale of the vortex is increasingly smaller consistent with its higher associated frequency as depicted in Figure 7(d). Therefore, applying the SPDMD method to atmospheric/marine motion analysis is helpful to further understand the multiscale characteristics of atmospheric motion. The SPDMD method can decompose the 200 hPa jet flow into different modes associated with a particular frequency and growth rate which can be used to identify the stability of the model. The low-frequency modes represent the large-scale motion structures, while the high-frequency modes highlight the small-scale features. Similar to the POD method, the multiscale features captured by the different DMD modes represent different spatial scales of the motion components for the jet flow, and a linear combination of these modes can represent the main characteristics of the jet flow [18]. Compared with the POD spatial modes depicted in Figure 3, it can be seen that the SPDMD method selects a smaller number of modes that capture nearly the same dynamic structures of the 200 hPa jet flow as the POD method.Figure 7
The five DMD modes identified by the SPDMD from the 200 hPa jet flow (a–e correspond to the 1–5th modes, velocity is visualized by vector plots, and the background color contours visualize the vertical vorticity component).
(a)(b)(c)(d)(e)


## 3.2. Reconstruction of the 200 hPa Jet

To further demonstrate the extraction effect of POD/SPDMD methods on the 200 hPa jet flow characteristics, we consider constructing the POD/SPDMD reduced-order models, respectively. In POD, to reconstruct 99% of the total fluctuation energy, the POD modes are truncated atrPOD=26. While in SPDMD, there are only ten modes (modes shown in Figure 6, DM1, and DM2,3,4,5 with their complex conjugates) attempted. The initial transient flow field at the sampling period T=60h,124h is depicted in Figure 8.Figure 8
Initial velocity field of 200 hPa jet flow at (a)T = 60 h, (b) T = 124 h (velocity visualized by vector plots and the background color contours visualize the vorticity component).
(a)(b)The POD/SPDMD reconstruction flow fields atT=60h,124h are presented in Figures 9 and 10, respectively. The highly similar features between the graphics show that the POD/SPDMD reduced-order models successfully recover the original dynamic process to an impressive degree, which indicates that both approaches are effective in detecting and extracting the relevant dynamic features underlying the 200 hPa jet flow. To quantitatively evaluate the quality of the POD/SPDMD reduced-order models develop by equations (12) and (27), the average relative error is defined as follows:(40)error=1N∑i=1NXi−XROMFiXFi,where Xi represent a vector of the original snapshots matrix at the time i and XROMi is an approximation from the POD/SPDMD reduced-order models at the same time; the results are listed in Table 2.Figure 9
POD reconstruction of the 200 hPa jet flow field (26 POD modes) at (a)T = 60 h and (b) T = 124 h (velocity visualized by vector plots and the background color contours visualize the vorticity component).
(a)(b)Figure 10
SPDMD reconstruction of the 200 hPa jet flow field (11 SPDMD modes) at (a)T = 60 h and (b) T = 124 h (velocity visualized by vector plots and the background color contours visualize the vorticity component).
(a)(b)Table 2
The average relative error of the POD/SPDMD reduced-order model.
Number of modesError ofROMDMDError ofROMPOD11 modes0.09710.110017 modes0.06710.084325 modes0.05930.064047 modes0.04430.0383As shown in Table2, compared with the POD method, the results of the SPDMD method show that it has a smaller average relative error when the number of modes is small.


## 4. Conclusion

In this paper, an advanced modal decomposition method named sparsity-promoting DMD (SPDMD) has been introduced into atmospheric motion flow-field analysis. And its performance and advantages are demonstrated on sequence snapshots of the stratospheric 200 hPa jet flow firstly and compared with the POD method. The SPDMD method can successfully identify the most important dynamic modes associated with their respective frequencies and growth rates, and the characteristics of each model, including the energy distribution and stability. The different dynamic modes with different spatial scales, which represent the motion components of different spatial scales of the jet flow, and the superposition of these modes can represent the main characteristics of the jet flow. Compared with the traditional modal decomposition method POD, fewer sparsity-promoting DMD modes can capture nearly the same coherent structures and provide a high-quality approximation of the original jet flow field.



---

*Source: 8383657-2020-09-17.xml*