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 Px=(x'1,
,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,
,lmax (lmax is the maximum iteration number), where Py=(y'1,y'2,
,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 (=t1,
,tS).
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 (Px); 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+,k1/2)(R+,k1/2) T=Rk. 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 Py,k=(y'k,1,
,y'k,N) and y'k,j=Lk(x b+x'j)-Lk(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 xa,iδ'=Px,ρβρlmax.
Step 4: Compute ||xa,iδ'$$ and then μi for the ith possible observation area Ωi and determine the sensitive areas in decreasing order of μi.