An Adjoint-Free CNOP-4DVar Hybrid Method for Identifying Sensitive Areas in Targeted Observations: Method Formulation and Preliminary Evaluation

The inclusion of a large amount of observation data into numerical weather prediction (NWP) models by data assimilation has become increasingly vital for accurately modeling high-impact weather events (HWEs) (Joly et al., 1997, Joly et al., 1999; Langland et al., 1999a, b, Aberson, 2003; Wu et al., 2007a; Aberson et al., 2011; Chou et al., 2011). On the other hand, acquiring sufficiently dense field observations to fully cover the vast (spatial) areas of HWEs is prohibitively expensive (Mu, 2013). To circumvent the dilemma, a targeted (or adaptive) observation strategy was proposed to obtain a limited number of observations distributed over sensitive areas. This strategy has been considered as a significant contribution to the field of weather forecasting. Recall that, to better predict an event at a future verification time t_v, the idea of targeted observations is to make additional observations at a proceeding future target time t_a (t_a<t_v) in sensitive areas. These additional observations are often acquired by a data assimilation system to provide a more reliable initial state and, as a result, prediction errors at the verification time t_v are substantially reduced (Majumdar, 2016). To realize this idea, identifying the sensitive (spatial) areas is the first and the most important step in targeted observations.

Existing methods for identifying sensitive areas in targeted observations can be divided into two groups. The first group of methods aims at capturing the most sensitive initial errors, which have the greatest influence on prediction accuracy. Areas containing a high concentration of the most sensitive initial errors are identified as sensitive areas by these methods. Examples of the first group of methods include the quasi-inverse linear method (Pu et al., 1997), the linear singular vector (LSV; Palmer et al., 1998), the sensitivity vector (Langland and Rohaly, 1996), ensemble sensitivity methods (Ancell and Hakim, 2007), and the conditional nonlinear optimal perturbation (CNOP) method (Mu et al., 2003). In particular, CNOP extends the LSV into a nonlinear regime by overcoming the limitations associated with the linearity assumption in the LSV and obtains promising targeted observation results (Mu et al., 2009; Zhou and Mu, 2011). The second group of methods aims at directly reducing prediction uncertainty by assimilating observations throughout all possible observation areas (e.g. Bishop and Toth, 1999; Bishop et al., 2001; Daescu and Navon, 2004; Hossen et al., 2012; Zhang et al., 2016). These methods compare the error reductions in all possible observation areas based on the extent of the decrease in the forecast error variance in each area, and regard the areas with the greatest assimilation corrections as sensitive areas. Among the methods in the second group, ensemble Kalman filter-based methods have been widely adopted due to their implementation simplicity and low computational cost. On the other hand, these methods were formulated under the linearity assumption, which limits their applicability to general nonlinear NWP models.

It is concluded that the first group of methods focuses only on initial errors; they do not consider the relationships between those errors and the observations as well as the impact of the data assimilation methods used. On the other hand, the second group of methods do not specifically take initial errors into consideration.

It seems a natural idea is to combine the two groups of methods to take the advantages of both types while avoiding their deficiencies. Indeed, the need for developing such methods has already been noticed (Mu, 2013). The primary objective of this paper is to attempt to meet the need and to fill the void of missing such methods. Specifically, we propose a CNOP-4DVar (four-dimensional variational) hybrid method to identify sensitive areas in targeted observations. In this hybrid method, we first compute the CNOP $x'_\delta\in\mathbb{R}^{m_x}$ (where m_x is the dimension of the vector x'_δ and the superscript '‘' denotes the perturbation) of the background state $x_b\in \mathbb{R}^{m_x}$ for the verification time t_v and utilize x _b+x'_δ as the most deviant "true" initial field to initialize the forecast model, and to obtain the most deviant "true" model integrations and the possible corresponding observations. We then adopt a 4DVar-based approach (Rabier et al., 2000) to determine the sensitive areas by utilizing the possible observations obtained in the first step. An important advantage of the proposed hybrid approach is that it can handle both linear and nonlinear models because both 4DVar and CNOP are based on nonlinear optimization formulations. Moreover, it is clear that this CNOP-4DVar method can organically integrate errors from all three (i.e., the initial, target and verification) stages into the overall system. Moreover, the proposed hybrid method also has a modular property, which makes the method easily expandable. For example, the 4DVar component in the method may be replaced by the ensemble transform Kalman filter (ETKF) method (Bishop et al., 2001), which then results in a new hybrid CNOP-ETKF method for identifying sensitive areas in targeted observations.

From the definition it is clear that the proposed hybrid CNOP-4DVar method comprises two nonlinear optimization sub-problems——namely, the CNOP and 4DVar problems. The default approach for solving these sub-problems is the adjoint-based approach. However, it is well-known that developing and maintaining the adjoint model is computationally very expensive. Moreover, the adjoint model is often difficult to obtain when the forecast model is highly nonlinear and/or when the model physics involve discontinuous parameters (Xu, 1996). It is well understood that requiring the use of an adjoint model has seriously limited the broader applicability of adjoint-based methods (Tian et al., 2008, Tian et al., 2010, Tian et al., 2011, Tian et al., 2016, Tian et al., 2018; Tian and Feng, 2015, 2017). To alleviate the dependence of using an adjoint model, several ensemble-based approaches have recently been developed for the CNOP and 4DVar problems. Among them, the nonlinear least-squares (NLS)-based ensemble method with a penalty strategy for the CNOP problem (NLS-CNOP; Tian et al., 2016; Tian and Feng, 2017) and the NLS-4DVar method for the 4DVar problems (Tian et al., 2011, Tian et al., 2018; Tian and Feng, 2015) have shown promising and competitive performance due to their implementation simplicity, high accuracy, and lack of dependence on tangent and adjoint models. Both methods will be employed to solve the CNOP and 4DVar sub-problems in the proposed hybrid CNOP-4DVar method.

The rest of this paper is organized as follows: In Section 2, we present our CNOP-4DVar hybrid formulation for identifying sensitive areas in targeted observations. Mathematically, it is a double constrained optimization problem. In section 3, we introduce ensemble-based NLS (adjoint-free) solvers/algorithms for solving the CNOP and 4DVar sub-problems in the CNOP-4DVar hybrid formulation; together they give rise to a hybrid CNOP-4DVar method for solving the problem of identifying sensitive areas. In section 4, we present a comprehensive performance evaluation for the proposed CNOP-4DVar method and compare it with the CNOP and CNOP-ETKF methods based on 10 000 observing system simulation experiments (OSSEs) on the shallow-water equation model. Finally, we close the paper with a summary and concluding remarks in section 5.

The targeted (or adaptive) observation strategy involves obtaining a limited number of observations distributed over sensitive areas, where a high concentration of the most sensitive initial errors occurs.

Its key idea is to make additional observations at a proceeding future target time t_a (t_a<t _v) in sensitive areas. To formulate our CNOP-4DVar hybrid formulation for the problem of identifying sensitive areas, we first apply this idea to compute the CNOP x'_δ of the background state x_b for the verification time t _v by maximizing the following cost function: \begin{equation} J({x}'_\delta)=\mathop{\max}\limits_{\|{x}'\|\le\delta}\|M_{t_0\to t_{\rm v}}({x}_{\rm b}+{x}')-M_{t_0\to t_{\rm v}}({x}_{\rm b}+{x}')\| , \ \ (1)\end{equation} where $M_{t_0\to t_v}$ is the nonlinear prediction model from t₀ to t_v, x' is an initial perturbation (IP) superposed on x _b, ||x'||≤δ is the constraint on the IP with ||·|| denoting a given norm, and δ is a positive constant (Mu et al., 2003). (Qin and Mu, 2011) described how to choose δ and to scale the multivariate variables appropriately in real numerical models. (Mu et al., 2003) demonstrated that the CNOP x'_δ grows faster than other types of perturbations (including those between the true and background states) due to the dynamic and physical mechanisms of HWEs.

Secondly, we define x_b+x'_δ as the most deviant "true" initial field and integrate the forecast model $x_{t_a}=M_{t_0\to t_a}(x_b+x'_{\delta})$ to produce the corresponding observations $y_{obs,k}^{\delta,i}\in \mathbb{R}^{m_y,k}$ for all possible observation areas Ω_i (i=1,…,n _s, where n _s is the total number of all possible observation areas Ω_i) at the targeted observation time (s) t _a.

Thirdly, we propose to assimilate the above possible observations y _obs,k^δ,i from the ith possible observation areas Ω_i by maximizing the following 4DVar objective functional: \begin{eqnarray} E({x}')&=&\frac{1}{2}({x}')^{\rm T}{B}^{-1}({x}')+\frac{1}{2}\sum_{k=0}^S[H_kM_{t_0\to t_{\it k}}({x}_{\rm b}+{x}')-{y}_{{\rm obs},k}^{\delta,i}]^{\rm T}\nonumber\\ &&\cdot{{R}}_k^{-1}[H_k M_{t_0\to t_k}({x}_{\rm b}+{x}')-{y}_{{\rm obs},k}^{\delta,i}] , \ \ (2)\end{eqnarray} where x' denotes the perturbation of the background field x _b at t₀, the superscript `T' stands for the matrix transpose operation, t_k denotes the target time, m_y,k is the dimension of the observational vector y_obs,k^δ,i, b is the background value, S+1 is the total number of targeted observation time points in the assimilation window, H_k is the observation operator, and matrices B and R_k are respectively the background and observation error covariance matrices.

Finally, the analysis increment x_a,i^δ' (where the superscript δ' denotes x_a,i^δ' corresponding to x'_δ) for the ith possible observation area Ω_i is obtained as the maximizer of E(x') and the areas with the highest concentrations of analysis increment indexes μ_i ($\mu_i\triangleq\frac{\|x_{a,i}^{\delta'}\|}{\begin{matrix}{max\|x_{a,i}^{\delta'}\|} \\ {\Omega_i}\end{matrix}}$ and 0≤ μ_i≤ 1) are designated as sensitive areas.

The above proposed CNOP-4DVar formulation for the problem of identifying sensitive area contains two optimization sub-problems, which must be solved numerically to implement the formulation. Solving these optimization problems could be very expensive if the adjoint-based approach is employed. To circumvent this difficulty, we propose to solve these optimization sub-problems by using an ensemble-based NLS approach developed by Tian and Feng (2015, 2017) and (Tian et al., 2016).

Specifically, we propose the following four-step algorithm to solve the problem of identifying sensitive areas for an ensemble-based forecast system based on the CNOP-4DVar formulation:

Step 1: Solve the CNOP problem using the NLS-CNOP approach. The CNOP problem (1) is first approximated by a sequence of unconstrained optimization problems (UOPs, parameterized by the small positive parameter ε) using a novel penalty strategy as follows (see Tian and Feng, 2017): \begin{align*} J_\varepsilon({x}'_\ast)=\min J_\varepsilon({x}') ,\tag*{(3)} \ \ (3)\end{align*} where \begin{align*} {y}'&=M_{t_0\to t_{\rm v}}({x}_{\rm b}+{x}')-M_{t_0\to t_{\rm v}}({x}_{\rm b}) ,\tag*{(4a)} \ \ (4a)\\ J_\varepsilon&=\frac{1}{2}\left(\frac{1}{\|{y}'\|^2}\cdot\frac{1}{\|{y}'\|^2}+f_\delta\right) ,\tag*{(4b)} \ \ (4b)\\ f_\delta:&=\varepsilon^{-1}(\|{x}'\|^2-\delta^2)^+:=\left\{ \begin{array}{l@{\quad}l} 0 & 0\le \|{x}'\|\le\delta\\ \varepsilon^{-1}\|{x}'\|^2 & \|{x}'\|>\delta \end{array} \right. ,\tag*{(4c)} \ \ (4c)\end{align*} in which $y'\in\mathbb{R}^m_y$ and $x'_\ast\in\mathbb{R}^m_x$ denote a solution to UOP (3). It can be shown that for a small enough ε the solution x'_* of UOP (3) is a close approximation to the solution x'_δ of the CNOP problem (1). We note that this procedure only guarantees finding one solution to CNOP (1), even if it has multiple solutions (Tian and Feng, 2017), and the computed solution is determined by the starting value of the iterative method used to solve the nonlinear problem.

Following the usual ensemble-based approach (e.g., Tian et al., 2016), we begin by preparing an ensemble of N linearly independent IPs (x'_j, j=1,…,N) and their corresponding prediction increments $(y'_j=M_{{t_0}\to t_\rm v}(x_\rm b+x'_j)-M_{{t_0}\to t_\rm v}(x_\rm b)$, j=1,…,N) at time t _v. Assume that any IP x' can be expressed by the linear combination of the sampling IPs (x'_j, j=1,…,N), i.e.: \begin{align*} {x}'={P}_x{v} ,\tag*{(5)} \end{align*} where P_x=(x'₁,…,x'_N) and $v[=(v_1,\cdots,v_N)^{T}]$ is the coefficient vector of the sampling IPs.

Substituting Eq. (5) into Eq. (3) expressing the cost function in terms of the new control variable v yields: \begin{align*} J_\varepsilon({v})=\frac{1}{2}Q^{\rm T}({v})Q({v}) ,\tag*{(6a)} \end{align*} where \begin{align*} Q({v})=\left( \begin{array}{c} \frac{1}{\|{y}'({P}_x{v})\|^2}\\[3mm] f_\delta^{1/2} \end{array} \right) ,\tag*{(6b)} \end{align*} and \begin{align*} f_\delta^{1/2}=:\left\{ \begin{array}{l@{\quad}l} 0 & 0\le\|{P}_x{v}\|\le\delta\\[1mm] \frac{1}{\sqrt\varepsilon}{P}_x{v} & \|{P}_x{v}\|>\delta \end{array} \right. .\tag*{(6c)} \end{align*}

Thus, we have transferred UOP (3) into the above NLS formulation, which can be solved by using the Gauss-Newton iterative method (Dennis and Schnabel, 1996). Before doing that, we first need to do some necessary modifications to guarantee its robustness. Let \(\upsilon\) be a small positive number defined as follows (see Tian and Feng, 2017 for further details): \begin{align*} &\left[2\frac{1}{\|{y}'^{l-1}\|^2}{P}_y^{\rm T}({y}'^{l-1})({y}'^{l-1})^{\rm T}{P}_y+\upsilon {I}\right]\Delta {v}^l={P}_y^{\rm T}{y}'^{l-1} ,\\ &0\le\|{P}_x{v}^{l-1}\|\le\delta ,\tag*{(7a)} \end{align*} and \begin{align*} &\left[4\frac{\varepsilon}{\|{y}'^{l-1}\|^8}{P}_y^{\rm T}({y}'^{l-1})({y}'^{l-1})^{\rm T}{P}_y+{P}_x^{\rm T}{P}_x+\upsilon {I}\right]\Delta {v}^l\\ =\,&2\frac{\varepsilon}{\|{y}'^{i}\|^6}{P}_y^{\rm T}{y}'^{l-1}-{P}_x^{\rm T}{P}_x{v}^{l-1},\quad \|{P}_x{v}^{s-1}\|>\delta ,\tag*{(7b)} \end{align*} for l=1,2,…,l_max (l_max is the maximum iteration number), where P_y=(y'₁,y'₂,…,y'_N) and I is an N× N identity matrix.

An efficient localization process is often used in this ensemble-based approach to filter out spurious correlations resulting from inadequate sampling or false long-range correlations (Houtekamer and Mitchell, 2001). After the localization process is done, the Gauss-Newton iteration (7) can be rewritten as: \begin{align*} \Delta {v}_\rho^l=({P}_y^\ast<e>{\rho}_y)^{\rm T}{y}'^{l-1},\quad 0\le\|{P}_{x,\rho}{v}_\rho^{l-1}\|\le\delta ,\tag*{(8a)} \end{align*} and \begin{align*} \Delta {v}_\rho^l= & 2\frac{\varepsilon}{\|{y}'^{l-1}\|^6}({P}_y^{\#}< e>{\rho}_y)^{T}{y}'^{l-1}-({P}_x^\ast < e>{\rho}_x)^{T}{P}_{x,\rho}{v}^{l-1} ,\\ &\|{P}_{x,\rho}{v}_\rho^{l-1}\|>\delta ,\tag*{(8b)} \end{align*} where \begin{align*} {P}_y^\ast &={P}_y\left[2\frac{1}{\|{y}'^{l-1}\|^2}{P}_y^{\rm T}({y}'^{l-1})({y}'^{l-1})^{\rm T}{P}_y+\upsilon{I}\right]^{-1} ,\\ {P}_y^{\#}&={P}_y\left[4\frac{\varepsilon}{\|{y}'^{l-1}\|^8}{P}_y^{\rm T}({y}'^{l-1})({y}'^{l-1})^{\rm T}{P}_y+{P}_x^{\rm T}{P}_x+\upsilon{I}\right]^{-1} ,\\ {P}_x^\ast&={P}_x\left[4\frac{\varepsilon}{\|y'^{l-1}\|^8}{P}_y^{\rm T}({y}'^{l-1})({y}'^{l-1})^{\rm T}{P}_y+{P}_x^{\rm T}{P}_x+\upsilon {I}\right]^{-1} , \end{align*} and \begin{align*} {P}_{x,\rho}=({P}_x<e>{\rho}_x) . \end{align*} Here, \(\rho_x\rho_x^\rm T=C\in\mathbb{R}^m_x\times m_x\) is the spatial correlation matrix (Tian et al., 2016; Zhang and Tian, 2018a). \(\rho_y\in \mathbb{R}^m_y\times r\) should be computed together with \(\rho_x\in \mathbb{R}^m_x\times r\), and r is the selected truncation mode number (Zhang and Tian, 2018a). Consequently, we obtain v_{ρlmax=vρlmax-1+∆vρlmax} and x'_{δ=Px,ρvρlmax}. The subscript `ρ' represents "localization".

Step 2: Integrate the forecast model $x_{t_\rm a}=M_{t_{0}\to t_\rm a}(x_\rm b+x'_\delta)$ to obtain the most deviant observations y _obs,k^δ,i at the target time(s) t _a (=t₁,…,t_S).

Step 3: Solve the 4DVar problem (2) with observations y _obs,k^δ,i using the NLS-4DVar. For ease of presentation, we first rewrite Eq. (2) as follows: \begin{align*} J({x}')=\,&\frac{1}{2}({x}')^{\rm T}{B}^{-1}({x}')+\frac{1}{2}\sum_{k=1}^S[L'_k({x}')-{y}_{{\rm obs},k}^{\delta,i'}]^{\rm T}\\[1mm] &\cdot{R}_k^{-1}[L'_k({x}')-{y}_{{\rm obs},k}^{\delta,i'}] ,\ \ (9) \end{align*} where \begin{align*} L'_k({x}')=L_k({x}_{\rm b}+{x}')-L_k({x}_{\rm b}) ,\quad {y}_{{\rm obs},k}^{\delta,i'}={y}_{{\rm obs},k}^{\delta,i}-L_k({x}_{\rm b}) , (9)\end{align*} and $L_k=H_k M_{t_0\to t_k}$.

Similarly, we assume that the analysis increment x' can be expressed as a linear combination of the sampling IPs (P_x); namely, \begin{align*} {x}'={P}_{x}{\beta} ,\tag*{(10)} \end{align*} where $\beta[=(\beta_1,\cdots,\beta_N)^\rm T]$ is the vector of coefficients of the sampling IPs. Substituting Eq. (10) and the ensemble-estimated $B\approx\frac{P_{x}P_{x}^{T}}{N-1}$ (Evensen, 1994) into Eq. (9) and expressing the cost function in terms of the new control variable β yields \begin{align*} J({\beta})=\,&\frac{1}{2}(N-1){\beta}^{\rm T}{\beta}+\frac{1}{2}\sum_{k=0}^S[L'_k({P}_x{\beta})-{y}_{{\rm obs},k}^{\delta,i'}]^{\rm T}\\ &{R}_{k,i}^{-1}[L'_k({P}_x{\beta})-{y}_{{\rm obs},k}^{\delta,i'}] .\tag*{(11)} \end{align*}

Equation (11) can be further rewritten as the following quadratic cost function: \begin{align*} J({\beta})=\frac{1}{2}Q({\beta})^{\rm T}Q({\beta}) ,\tag*{(12)} \end{align*} where \begin{align*} Q({\beta})=\left( \begin{array}{c} \sqrt{N-1}{\beta}\\ {R}_{+,0}^{-1/2}[L'_0({P}_{x,\rho}{\beta})-{y}_{{\rm obs},0}^{\delta,i'}]\\ \vdots\\ {R}_{+,S}^{-1/2}[L'_S({P}_{x,\rho}{\beta})-{}_{{\rm obs},S}^{\delta,i'}] \end{array} \right) ,\tag*{(13)} \end{align*} and (R_+,k^1/2)(R_+,k^1/2) ^T=R_k. The first-derivative (or Jacobian) matrix J _acQ(β) of Q(β) is given by \begin{align*} J_{\rm ac}Q({\beta})=\frac{\partial Q({\beta})}{\partial{\beta}}\approx\left( \begin{array}{c} \sqrt{N-1}{I}\\ {R}_{+,0}^{-1/2}{P}_{y,0}\\ \vdots\\[1.5mm] {R}_{+,S}^{-1/2}{P}_{y,S} \end{array} \right) ,\tag*{(14)} \end{align*} where P_y,k=(y'_k,1,…,y'_k,N) and y'_k,j=L_k(x _b+x'_j)-L_k(x _b). The nonlinear least-squares formulation [Eqs. (12)-(14)] is then solved by the Gauss-Newton iterative method [see (Tian and Feng, 2015) and (Tian et al., 2018) for further details], as follows: \begin{align*} {\beta}^l=\,&{\beta}^{l-1}-\left[(N-1){I}+\sum_{k=0}^S({P}_{y,k})^{\rm T}{R}_k^{-1}({P}_{y,k})\right]^{-1}\\ &\times\left[\sum_{k=0}^S({P}_{y,k})^{\rm T}{R}_k^{-1}(L'_k({P}_x{\beta}^{i-1})-{y}_{{\rm obs},k}^{\delta,i'})+(N-1){\beta}^{l-1}\right] ,\tag*{(15a)} \end{align*} and therefore, \begin{align*} &\left[(N-1){I}+\sum_{k=0}^S({P}_{y,k})^{\rm T}{R}_k^{-1}({P}_{y,k})\right]\Delta{\beta}^l\\ =\,&-\left[\sum_{k=0}^S({P}_{y,k})^{\rm T}{R}_k^{-1}(L'_k({P}_x{\beta}^{i-1})-{y}_{{\rm obs},k}^{\delta,i'})+(N\!-\!1){\beta}^{l-1}\right] .\tag*{(15b)} \end{align*}

Following (Tian and Feng, 2015), we obtain \begin{align*} \Delta{\beta}^l=\sum_{k=0}^S({P}_{y,k}^\ast)^{\rm T}L'_k({P}_x{\beta}^{l-1})\!+\!\sum_{k=0}^S({P}_{y,k}^{\#})^{\rm T}[{y}_{{\rm obs},k}^{\delta,i'}\!-\!L'_k({P}_x{\beta}^{l-1})] ,\tag*{(16)} \end{align*} where \begin{align*} {P}_{y,k}^\ast\!\!=\!\!-(N\!-\!1){P}_{y,k}\!\left[\!\sum_{k=0}^S\!{P}_{y,k}^{\rm T}{R}_k^{-1}\!{P}_{y,k}\!\!+\!(N\!\!-\!\!1){I}\!\right]^{-1}\!\left[\!\sum\limits_{k=0}^S{P}_{y,k}^{\rm T}{P}_{y,k}\!\right]^{-1}, \end{align*} and \begin{align*} {P}_{y,k}^{\#}={R}_k^{-1}{P}_{y,k}\left[\sum_{k=0}^S{P}_{y,k}^{\rm T}{R}_k^{-1}{P}_{y,k}+(N-1){I}\right]^{-1} . \end{align*}

The analysis increment for the ith possible observation area is thus given by \begin{align*} {x}_{a,i}^{\delta'}={P}_x{\beta}^{l_{\max}} .\tag*{(17)} \end{align*}

Like the ensemble-based approach to the CNOP problem shown in Step 1, we adapt an efficient localization scheme (Zhang and Tian, 2018a) that was originally proposed for the NLS-4DVar by (Tian et al., 2018) by transforming Eq. (16) into the following form: \begin{align*} \Delta{\beta}_\rho^l=\,&\sum_{k=0}^S({P}_{y,k}^\ast<e>\rho_y)^{\rm T}L'_k({P}_{x,\rho}{\beta}_\rho^{l-1})\\ &+\sum_{k=0}^S({P}_{y,k}^{\#} <e>\rho_y)^{\rm T}[{y}_{{\rm obs},k}^{i,'}-L'_k({P}_{x,\rho}{\beta}_\rho^{l-1})] ,\tag*{(18)} \end{align*} and thus x_a,i^δ'=P_x,ρβ_ρlmax.

Step 4: Compute ||x_a,i^δ'$$ and then μ_i for the ith possible observation area Ω_i and determine the sensitive areas in decreasing order of μ_i.

We utilize a two-dimensional (2D) shallow-water equation model (Qiu et al., 2007) to evaluate the proposed CNOP-4DVar hybrid method, particularly in comparison with the CNOP method and the CNOP-ETKF hybrid method. The 2D shallow-water equations are given in the f-plane by \begin{align*} \frac{\partial u}{\partial t}&=-u\frac{\partial u}{\partial x}-v\frac{\partial u}{\partial y}+fv-g\frac{\partial h}{\partial x} ,\tag*{(19a)}\\ \frac{\partial v}{\partial t}&=-u\frac{\partial u}{\partial x}-v\frac{\partial u}{\partial y}-fu-g\frac{\partial h}{\partial y} ,\tag*{(19b)}\\ \frac{\partial h}{\partial t}&=-u\frac{\partial (h-h_{\rm s})}{\partial x}-v\frac{\partial(h-h_{\rm s})}{\partial y}-(H+h-h_{\rm s})\left(\frac{\partial u}{\partial x}+\frac{\partial u}{\partial y}\right) ,\tag*{(19c)} \end{align*} where g is the acceleration of gravity, f=7.272× 10^-5 s^-1 is the Coriolis parameter, H=3000 m is the basic-state depth, $h_s=h_0 sin(4\pi x/L_x)[\sin(4\pi x/L_y)]^2$ is the terrain height, h₀=250 m and D=L_x=L_y=13200 km are the lengths of the two sides of the computational domain, respectively, and the uniform grid size is taken as d=∆_x=∆_y=300 km. The computational domain is partitioned into a square mesh with 45 grid points in each coordinate direction and the periodic boundary conditions are imposed at x=0, L_x and y=0, L_y. The spatial/local time derivatives are discretized by the second order central finite difference scheme and two-step backward difference scheme of (Matsuno, 1966), respectively. We choose the time step size ∆ t=360 s (or 6 min). The model state vector is represented by height h and horizontal velocity components u and v at the grid points.

Similar to (Tian and Xie, 2012), we prepare the background (initial) state x _b (Fig. 1) by integrating the shallow-water equation model with the imperfect h₀=0 m from the following initial conditions (ICs): \begin{align*} h&=360\left[\sin\left(\frac{\pi y}{D}\right)\right]^2+120\sin\left(\frac{2\pi x}{D}\right)\sin\left(\frac{2\pi y}{D}\right) ,\tag*{(20a)}\\ u&=-f^{-1}g\frac{\partial h}{\partial y},\ {\rm and}\ v=-f^{-1}g\frac{\partial h}{\partial x} ,\tag*{(20b)} \end{align*} at t=-60 h for 60 hours. We assume the spatially averaged root-mean-square errors (RMSEs) to be 23.4 m, 1.53 m and 2.58 m s^-1 for the h, u and v fields, respectively. The verification time is 12 h and the target times are 3, 6, and 9 h, respectively. Additionally, $\|\cdot\|\doteq\sqrt{\sum_{i, j=1}^{45}\left[\left(\frac{u_{i,j}}{u_\delta}\right)^2+\left(\frac{v_{i, j}}{v_\delta}\right)^2+\left(\frac{h_{i, j}}{h_\delta}\right)^2\right]}$ and $\delta=\sqrt{3}$ for the CNOP problem. At this stage, we utilized the NLS-CNOP method to compute the CNOP x'_δ of the background state x _b for the verification time t _v=12 h, which in turn determines the CNOP-type sensitive areas. Subsequently, the most deviant observations y _obs,k^δ,i can be determined at the target times (i.e., t _a=3, 6 and 9 h) by $x_{t_a}=M_{t_0\to t_a}(x_b+x'_\delta)$, where $M_{t_0\to t_a}$ is the shallow-water equation model with h₀=250. Thus, we apply the NLS-4DVar and ETKF methods to determine the sensitive areas among all possible observation areas Ω_i (i=1,…,n_s and n_s=45× 45) using the observations y _obs,k^δ,i obtained at the target times.

Figure 1.  Spatial distribution of the ICs at 0 h: (a) u (m s^-1); (b) v (m s^-1) and h (m s^-1).

The second part of our numerical experiments involves designing a large number of 4DVar-based OSSEs, by assimilating the synthetic observations in the identified sensitive areas. For this purpose, we first generate a basic "true" state x_t,0 by integrating the shallow-water equation model with the perfect h₀=250 m from the same ICs [see Eq. (20)] for 60 h. Subsequently, 10 000 groups of "true" initial states x_t,s (where the subscript 's' denotes the sth "true" initial state and s=1,2,…,10000) are generated by adding random constraint perturbations x'_t,s (||x'_t,s||≤δ) to the basic "true" states x_t,0 as x_t,s=x_t,0+x'_t,s. Model-produced "true" fields at 3, 6, 9 and 12 h, respectively, are therefore produced after 12-h integrations from the "true" initial x_t,s at the same starting time in the first part of our numerical experiments. The synthetic observations in the identified sensitive areas at the target times (t_k=3, 6, and 9 h) in the CNOP, CNOP-4DVar and CNOP-ETKF methods are generated by adding random noise to the model-produced "true" fields, respectively. The number of observations m_y,k is taken to be 225 for each target time (3, 6, and 9 h).

We evaluate the observations in the sensitive areas (consisting of 225 grid points) obtained by the CNOP, CNOP-4DVar and CNOP-ETKF methods using 10 000 4DVar-based OSSEs with the 10 000 "true" initial states x_t,s (s=1,…,10000) and their corresponding synthetic observations for a single observational time level (i.e., t₀=3, 6, or 9 h), and for all three observational time levels [i.e., t_k=3(k+1) h, k=0, 1, and 2]. The background x _b is the same as used in the first part of our numerical experiments. We then further compare the CNOP-4DVar observation strategy with the commonly used strategy involving the selection of sparse data points in the observation area.

In addition, we also use the NLS-4DVar approach for data assimilation in the second part of the numerical experiments. The standard parameters are an ensemble size N=90 and a covariance localization Schur radius of d_h,0=d_v,0=10× 30 km for all three ensemble-based (NLS-CNOP, NLS-4DVar, and ETKF) methods.

A comparison of the spatial distribution of the CNOP $u'_{\delta}[(u'_{δ},v'_{δ},h'_{δ})=x'_{δ}]$ with a verification time of 12h, and its corresponding model perturbations $u'_{t_k}[(u'_{t_k},v'_{t_k},h'_{t_k})=x'_{t_k} \text{and}\ x'_{t_k}=M_{{t_0}\to t_k}(x_b+x'_\delta)-M_{t_0\to t_k}(x_\rm b)]$ at the target times (t_k=3, 6, and 9 h) are shown in Fig. 2. The locations of areas with a high concentration of errors change with model integration. Therefore, we check whether the CNOP-determined sensitive areas (at the initial time) perfectly match with those at the target times. We observe that the sensitive areas identified by the CNOP-4DVar/ETKF methods are significantly different from those obtained by the CNOP method at all target times (Figs. 2 and 3). Furthermore, there is a large variability in the CNOP-4DVar/ETKF-determined sensitive areas at different target times (Figs. 3a-c and Figs. 3d and e). Moderate differences are also found between the CNOP-4DVar- and CNOP-ETKF-determined sensitive areas, even at the same target times. This is the consequence of using different optimization methods (i.e., NLS-4DVar vs ETKF) by the two methods to identify sensitive areas.

Figure 2.  Spatial distribution of (a) the CNOP u'_δ (m s^-1) and its corresponding model perturbations u'_tk (m s^-1) at (b) 3 h, (c) 6 h and (d) 9 h.

Figure 3.  Spatial distribution of the (a-c) CNOP-4DVar-determined and (d-f) CNOP-ETKF-determined analysis increment indexes at (a, d) 3 h, (b, e) 6 h and (c, f) 9 h.

We then quantitatively examine the performance of the three methods using normalized RMSEs, which are computed by the formula \begin{align*} r_{\rm sn}=\sqrt{\sum_{i,j=1}^{45}\left[\left(\frac{u_{i,j}^{a(b)}-u_{i,j}^t}{u_\delta}\right)^2+\left(\frac{v_{i,j}^{a(b)}-v_{i,j}^t}{v_\delta}\right)^2 +\left(\frac{h_{i,j}^{a(b)}-h_{i,j}^t}{h_\delta}\right)^2\right]} , \end{align*} where a, b and t indicate the assimilated, background and true solutions, respectively. We compare the background-derived RMSEs with the corresponding NLS-4DVar-derived RMSEs by assimilating the observations in the sensitive areas identified by the CNOP, CNOP-4DVar and CNOP-ETKF methods at a single targeted observation time (t₀=3, 6, or 9 h), respectively, and at all three target times (t_k=3, 6, and 9 h). The results show that all three methods are satisfactory in terms of the overall lower RMSEs compared with the background-derived RMSEs for the total 10 000 experiments. Generally, if the target time is closer to the verification time, the assimilation performance becomes better; and using multiple target times produces a substantially more accurate result than using a single target time. The proposed CNOP-4DVar hybrid method is the best performing method among the tested methods. It moderately surpasses the CNOP-ETKF method but substantially outperforms the CNOP method. Presumably, the sensitive areas with high concentrations of errors at the initial rather than target time of the CNOP method is credited for its inferior performance compared with the other two methods. The 4DVar nonlinear optimization of the CNOP-4DVar method is superior to the ETKF method adopted by the CNOP-ETKF approach, leading to lower RMSEs for the CNOP-4DVar method than for the CNOP-ETKF method, both at one and three target time(s) (Fig. 4).

Figure 4.  Comparisons between the background-derived RMSEs and the corresponding NLS-4DVar-derived RMSEs by assimilating the observations in the sensitive areas obtained by the (a-d) CNOP, (e-h) CNOP-4DVar, and (i-l) CNOP-ETKF methods at a single targeted observation time (t₀=3, 6, or 9 h), respectively, and at all three target times (t_k=3, 6, and 9 h).

The use of the maximum iteration number l_max for the NLS-4DVar solver in Step 3 of the proposed CNOP-4DVar method clearly influences to a large extent the computational accuracy and complexity. To determine the influence of the number of iterations on the CNOP-4DVar results, we evaluate the CNOP-4DVar-determined sensitive areas at target time t _a=9 h with l_max=1, 2, 3, and 4. As Fig. 5 shows, only small differences are observed when l_max≥ 2. We also compare RMSEs by assimilating the CNOP-4DVar-determined observations at t _a=9 h for l_max=2 and 3; Fig. 6 shows that their differences are mainly concentrated near the diagonal line y=x in both cases, which implies that they generally have a similar performance, and that the Gauss-Newton iterative method ensures that NLS-4DVar reaches its minimum after only one to two iterations, which is also consistent with the results of our previous studies (Tian and Feng, 2015; Tian et al., 2018).

Figure 5.  Spatial distribution of the CNOP-4DVar-determined analysis increment indexes at the single targeted observation time (t₀=9 h) for l_max=1, 2, 3, and 4, respectively.

Figure 6.  Comparisons between the NLS-4DVar-derived RMSEs by assimilating the observations in the sensitive areas obtained by the CNOP-4DVar method for l_max=2 and 3, respectively, at a single targeted observation time (t₀=9 h).

We also compare the proposed CNOP-4DVar method with the commonly used sparse observation strategy, by selecting 225 sparse grid points within the observation area. For the CNOP-4DVar observation strategy, we chose observation points (j=1,…,m_y,k), in decreasing order of ||x_a,i^δ,'||, with a grid size ≥ 3× 300 km. For the sparse observation strategy, observation points are chosen at sparsely selected grid points equally spaced every 3d=900 km in each coordinate direction. As expected, the proposed CNOP-4DVar hybrid method substantially outperforms the sparse strategy in all 10 000 experiments (Fig. 7). The results of these experiments show that the proposed hybrid method is indeed capable of achieving higher assimilation accuracy.

Figure 7.  Comparisons between the NLS-4DVar-derived RMSEs by assimilating the observations in the sensitive areas obtained by the CNOP-4DVar method and the commonly used strategy (CUS), respectively, at a single targeted observation time (t₀=9 h).

Existing methods for identifying sensitive areas in targeted observations can generally be divided into two groups, according to whether they aim to directly reduce initial errors or prediction errors; these two groups of methods are only weakly linked. In this paper we propose a hybrid method, called the CNOP-4DVar method, which links these two groups of methods and takes the advantages of both the CNOP and the 4DVar method. The proposed CNOP-4DVar hybrid method is capable of capturing the most sensitive IP, which causes the greatest perturbation growth at the time of verification, and its 4DVar capability allows it to identify sensitive areas by evaluating their assimilation effects for eliminating the most sensitive IP, which has already been captured by its CNOP capability. The adjoint-based approach is typically used to solve the CNOP and 4DVar nonlinear optimization problems, but it is computationally very expensive for large-scale problems. To avoid the limitations of the adjoint-based approach, we adopt two adjoint-free methods——namely, the NLS-CNOP and NLS-4DVar methods——to solve the CNOP and 4DVar sub-problems, respectively.

We carry out a comprehensive performance evaluation for the proposed CNOP-4DVar hybrid method by comparing its results with those obtained by the CNOP and CNOP-ETKF methods over 10 000 OSSEs based on the shallow-water equation. The results show that the proposed CNOP-4DVar method generally performs better than the CNOP-ETKF method, but substantially better than the CNOP method. Compared with the CNOP method, the CNOP-4DVar hybrid method has an additional advantage of being able to precisely locate sensitive areas at the target time. Our tests also show that the 4DVar-based nonlinear optimization used in the CNOP-4DVar hybrid method is superior to the ETKF method used in the CNOP-ETKF hybrid method.

Step 3 of the CNOP-4DVar hybrid method requires the use of the NLS-4DVar method with l_max≈ 2, throughout all possible observation areas. Admittedly, this requirement may be computationally costly, especially in complicated real-world nonlinear weather and climate models. This issue may be overcome and evaluated by the multigrid iterative scheme (Zhang and Tian, 2018b) and will be studied in a forthcoming work.