Formulations for Estimating Spatial Variations of Analysis Error Variance to Improve Multiscale and Multistep Variational Data Assimilation

When the coarse-resolution observations used in the first step of multiscale and multistep variational data assimilation become increasingly nonuniform and/or sparse, the error variance of the first-step analysis tends to have increasingly large spatial variations. However, the analysis error variance computed from the previously developed spectral formulations is constant and thus limited to represent only the spatially averaged error variance. To overcome this limitation, analytic formulations are constructed to efficiently estimate the spatial variation of analysis error variance and associated spatial variation in analysis error covariance. First, a suite of formulations is constructed to efficiently estimate the error variance reduction produced by analyzing the coarseresolution observations in oneand two-dimensional spaces with increased complexity and generality (from uniformly distributed observations with periodic extension to nonuniformly distributed observations without periodic extension). Then, three different formulations are constructed for using the estimated analysis error variance tomodify the analysis error covariance computed from the spectral formulations.The successively improved accuracies of these three formulations and their increasingly positive impacts on the two-step variational analysis (or multistep variational analysis in first two steps) are demonstrated by idealized experiments.


Introduction
Multiple Gaussians with different decorrelation length scales have been used at NCEP to model the background error covariance in variational data assimilation (Wu et al. [1], Purser et al. [2]), but mesoscale features are still poorly resolved in the analyzed incremental fields even in areas covered by remotely sensed high-resolution observations, such as those from operational weather radars (Liu et al. [3]).This problem is common for the widely adopted singlestep approach in operational variational data assimilation, especially when patchy high-resolution observations, such as those remotely sensed from radars and satellites, are assimilated together with coarse-resolution observations into a high-resolution model.To solve this problem, multiscale and multistep approaches were explored and proposed by several authors (Xie et al. [4], Gao et al. [5], Li et al. [6], and Xu et al. [7,8]).For a two-step approach (or the first two steps of a multistep approach) in which broadly distributed coarse-resolution observations are analyzed first and then locally distributed high-resolution observations are analyzed in the second step, an important issue is how to objectively estimate or efficiently compute the analysis error covariance for the analyzed field that is obtained in the first step and used to update the background field in the second step.To address this issue, spectral formulations were derived by Xu et al. [8] for estimating the analysis error covariance.As shown in Xu et al. [8], the analysis error covariance can be computed very efficiently from the spectral formulations with very (or fairly) good approximations for uniformly (or nonuniformly) distributed coarse-resolution observations and, by using the approximately computed analysis error covariance, the twostep analysis can outperform the single-step analysis under the same computational constraint (that mimics the operational situation).

Analysis Error Variance Formulations for
One-Dimensional Cases

Error Variance Reduction
Produced by a Single Observation.When observations are optimally analyzed in terms of the Bayesian estimation (see chapter 7 of Jazwinski [10]), the background state vector b is updated to the analysis state vector a with the following analysis increment: and the background error covariance matrix B is updated to the analysis error covariance matrix A according to where R is the observation error covariance matrix, d = y − h(b) is the innovation vector (observation minus background in the observation space), y is the observation vector, h() denotes the observation operator, and H is the linearized h().
For a single observation, say, at   in the one-dimensional space of , the inverse matrix (HBH T + R) −1 reduces to ( 2  +  2  ) −1 , so the th diagonal element of A in (1b) is simply given by where   =  2  /( 2  +  2  ),  2  (or  2  ) is the background (or observation) error variance,   () is the background error correlation function,   denotes the th point in the discretized analysis space   , and  is the number of grid points over the analysis domain.The length of the analysis domain is  = Δ, where Δ is the analysis grid spacing and  is assumed to be much larger than the background error decorrelation length scale .
Note that   () is a continuous function of , so (2) can be written into  2  () ≡  2  − Δ 2  () also as a continuous function of , where is the error variance reduction produced by analyzing a single observation at  =   .The error variance reduction in (3) decreases rapidly as | −   | increases, and it becomes much smaller than it peak value of    2  C 2  at  =   as | −   | increases to .This implies that the error variance reduction produced by analyzing  sparsely distributed coarse-resolution observations can be estimated by properly combining the error variance reduction computed by (3) for each coarse-resolution observation as a single observation.This idea is explored in the following three subsections for one-dimensional cases with successively increased complexity and generality: from uniformly distributed coarseresolution observations with periodic extension to nonuniformly distributed coarse-resolution observations without periodic extension.

Uniform Coarse-Resolution Observations with Periodic
Extension.Consider that there are  coarse-resolution observations uniformly distributed in the above analysis domain of length  with periodic extension, so their resolution is Δ co ≡ /.In this case, the error variance reduction produced by each observation can be considered as an additional reduction to the reduction produced by its neighboring observations, and this additional reduction is always smaller than the reduction produced by the same observation but treated as a single observation.This implies that the error variance reduction produced by analyzing the  coarse-resolution observations, denoted by Δ 2  (), is bounded above by ∑  Δ 2  (); that is, where ∑  denotes the summation over  for the  observations.The equality in (4) is for the limiting case of Δ co / → ∞ only.The inequality in (4) implies that the domainaveraged value of ∑  Δ 2  () is larger than the true averaged reduction estimated by Δ 2  ≡  2  −  2  , where  2  is the domain-averaged analysis error variance estimated by the spectral formulation for one-dimensional cases in Section 2.2 of Xu et al. [8].
The domain-averaged value of ∑  Δ 2  () can be computed by where ∫   denotes the integration over the analysis domain, ∑  denotes the summation over  for the  grid points, and  = Δ is used in the last step.By extending  2  ( −   ) with the analysis domain periodically, Δ 2  can be also estimated analytically as follows: where ∫  denotes the integration over the infinite space of , is used in the second to last step, and  1 ≡ ∫ 2  ()/ is used with Δ co ≡ / in the last step.For the double-Gaussian form of   () = 0.6 exp(− 2 /2 2 ) + 0.4 exp(−2 2 / 2 ) used in (5) of Xu et al. [8], we have  1 = (2) 1/2 (0.44/2 1/2 + 0.48/5 The analysis error variance,  2  (), is then estimated by As shown by the example in Figure 1 (in which  = 110.4km and  = 10 so Δ co = / = 11.04 km is close to  = 10km), the estimated  2  * () in ( 7) has nearly the same spatial variation as the benchmark  2  () that is computed precisely from (1b), although the amplitude of spatial variation of   monotonically toward its asymptotic upper limit of    2  = 20 m 2 s −2 ) as Δ co / decreases to 0.5 and then to Δ/ = 0.1 (or increases toward ∞), and this decrease (or increase) of the amplitude of spatial variation of  2  () with Δ co / is closely captured by the amplitude of spatial variation of the estimated  2   * () as a function of Δ co /.Using the estimated   * () in (7), the previously estimated analysis error covariance matrix, denoted by A  with its th element   ≡  2     (  −   ) obtained from the spectral formulations, can be modified into A  , A  , or A  with its th element given by or   ≡   + {  The formulation in (8a) is conventional, as in (2.1) of Purser et al. [2] or originally (11) of Rutherford [11], in which the covariance is modified by applying   * () separately to each entry (indexed by  and ) of   (  −   ) to retain the self-adjointness.The second equation in (8a) shows that the conventional approach can be viewed alternatively as   plus a correction term, the last term in (8a).Ideally, the correction term should completely offset the deviation of   from the true covariance, but the correction term in (8a) offsets only a part of the deviation.
For the case in Figure 1, the benchmark analysis error covariance matrix, denoted by A, is computed precisely from (1b) and is plotted in Figure 3, while the deviations of A  , A  , A  , and A  from the benchmark A are shown in Figures 4(a), 4(b), 4(c), and 4(d), respectively.As shown, the deviation becomes increasingly small when A  is modified successively to A  , A  , and A  .Note that the correction term in (8a) is   (  −  ) modulated by   * (  )  * (  )− 2   .This modulation has a chessboard structure, while the desired modulation revealed by the to-be-corrected deviation of A  in Figure 4(a) has a banded structure (along the direction of   +   = constant, perpendicular to the diagonal line).This explains why the correction term in (8a) offsets only a part of the deviation as revealed by the deviation of A  in Figure 4(b).On the other hand, the correction term in (8b) is modulated by  2  * (  /2 +   /2) −  2  .This modulation not only retains the self-adjointness but also has the desired banded structure, so the correction term in (8b) is an improvement over that in (8a), as shown by the deviation of A  in Figure 4(c) versus that of A  in Figure 4(b).However, as revealed by Figure 4(c), the deviation of A  still has two significant maxima (or minima) along each band on the two sides of the diagonal line of   =   , while the to-be-corrected deviation of A  in Figure 4(a) has a single maximum (or minimum) along each band.This implies that the function form of   (  −  ) is not sufficiently wide for the correction.As a further improvement, this function form is widened to   (  −   ) for the correction term in (8c), so the deviation of A  in Figure 4(d) is further reduced from that of A  in Figure 4(c).
When an estimated A is used to update the background error covariance in the second step for analyzing the highresolution observations in the nested domain, the accuracy of the second-step analysis depends not only, to a certain extent, on the number of iterations performed by the minimization algorithm but also on the accuracy of the estimated A over the nested domain plus its extended vicinities within the distance of 2  outside the nested domain.Here,   is the decorrelation length scale of   () defined by  2   ≡ [−  ()/ 2    ()]| =0 according to (4.3.10) of Daley [12], and   (=4.45 km for the case in Figures 1 and 3) can be easily computed as a by-product from the spectral formulation.Over this extended nested domain, the relative error (RE) of the estimated A  with respect to the benchmark A can be measured by where I  denotes the unit matrix in the subspace associated with the grid points in the extended nested domain and thus I  (A  − A)I  (or I  AI  ) is the submatrix of A  − A (or A) associated only with the grid points in the extended nested domain and ‖()‖  denotes the Frobenius norm of () defined by the square root of the sum of the squared absolute values of the elements of the matrix in () according to (2.2-4) of Golub and Van Loan [13].The REs of A  , A  , and A  can be measured by the same form of Frobenius norm ratio as that defined for A  in (9).The REs of A  , A  , A  , and A  are computed for the case in Figure 1 and listed in the first column of Table 1.As shown by the listed values, the RE becomes increasingly small when A  is modified successively to A  , A  , and A  , and this is consistent with and  3 plotted by color contours every 0.5 m 2 s −2 .Deviations of A  , A  , and A  from benchmark A are plotted by color contours every 0.2 m 2 s −2 in panels (b), (c), and (d), respectively.Here, A  is the previously estimated analysis error covariance matrix with its th element   ≡  2     (  −   ) obtained from the spectral formulation, while A  , A  , and A  are the newly modified estimates of A as shown in (8a), (8b), and (8c), respectively.also quantifies the successively reduced deviation shown in Figures 4(a)-4(d).

Nonuniform Coarse-Resolution Observations with Periodic Extension.
Consider that the  coarse-resolution observations are now nonuniformly distributed in the analysis domain of length  with periodic extension, so their averaged resolution can be defined by Δ co ≡ /.The spacing of a concerned coarse-resolution observation, say the th observation, from its right (or left) adjacent observation can be denoted by Δ co+ (or Δ co− ).Now we can consider the following two limiting cases.
First, we consider the case of Δ co+ → 0 with Δ co− = Δ co (or Δ co− → 0 with Δ co+ = Δ co ).In this case, the concerned th observation collapses onto the same point with its right (or left) adjacent observation, that is, the ( + 1)th [or ( − 1)th] observation.The two collapsed observations should be combined into one superobservation with a reduced error variance from  2  to  2  /2.The error variance reduction produced by this superobservation still can be estimated by (3) but with Table 1: Entire-domain averaged RMS errors (in ms −1 ) for the analysis increments obtained from SE, TEe, TEa, TEb, and TEc applied to the first set of innovations with periodic extension and consecutively increased , where  is the number of iterations.All the RMS errors are evaluated with respect to the benchmark analysis increment.The relative error (RE) of the estimated analysis error covariance for updating the background error covariance in the second step of the two-step analysis is listed with the experiment name in the first column for each two-step experiment.On the other hand, without super-Obbing, the error variance reduction produced by the two collapsed observations will be overestimated by (3) with

Experiment
By comparing (10b) with (10a), it is easy to see that this overestimation can be corrected if the error variance is inflated from  2  to  2  +  2  for each of the two collapsed observations.
Then, we consider the case of Δ co+ → 0 and Δ co− → 0. In this case, the concerned th observation collapses with its two adjacent observations, that is, the ( + 1)th and ( − 1)th observations.The three collapsed observations should be combined into one superobservation with a reduced error variance from  2  to  2  /3.The error variance reduction produced by this superobservation still can be estimated by (3) but with On the other hand, without super-Obbing, the error variance reduction produced by the three collapsed observations will be overestimated by (3) with By comparing (10d) with (10c), it is easy to see that this overestimation can be corrected if the error variance is inflated from  2  to  2  + 2 2  for each of the three collapsed observations.
The above results suggest that for the definition of This modification can improve the similarity of the spatial variation of ∑  Δ 2  () to that of the true error variance reduction, denoted by Δ 2  () ≡   6) for uniform coarse-resolution observations but with Δ co decreased to Δ omn (or increased to Δ omx ), where Δ omn (or Δ omx ) is the minimum (or maximum) spacing between two adjacent observations among all nonuniformly distributed coarse-resolution observations in the one-dimension analysis domain.By adjusting Δ 2 emx to Δ 2 mx and Δ 2 emn to Δ 2 mn , the error variance reduction can be estimated by where as in (7), except that Δ 2  () is computed by (12a) instead of (6).As shown by the example in Figure 5, the estimated  2  * () captures closely not only the maximum and minimum but also the spatial variation of the benchmark  2  () computed from (1b).Using this estimated   * (), the previously estimated A  from the spectral formulation can be modified into A  , A  , or A  with its th element given by the same formulation as shown in (8a), (8b), or (8c).For the case in Figure 5, the benchmark A is plotted in Figure 6, while the deviations of A  , A  , A  , and A  from the benchmark A are shown in Figures 7(a), 7(b), 7(c), and 7(d), respectively.As shown, the deviation becomes increasingly small when the estimated analysis error covariance matrix is modified successively to A  , A  , and A  .
As explained in Section 2.2, the accuracy of the secondstep analysis depends on the accuracy of the estimated A over the extended nested domain (i.e., the nested domain plus its extended vicinities within the distance of 2  on each side outside the nested domain), while the latter can be measured by the smallness of the RE of the estimated A with respect to the benchmark A, as defined for A  in (9).The REs of A  , A  , A  , and A  computed for the case in Figure 5 are listed in the first column of Table 2.As listed, the RE becomes increasingly small when A  is modified successively to A  , A  , and A  , which quantifies the successively reduced deviation shown in Figures 7(a)-7(d).

Nonuniform Coarse-Resolution Observations without
Periodic Extension.Consider that the  coarse-resolution observations are still nonuniformly distributed in the onedimensional analysis domain of length  but without periodic extension.In this case, their produced error variance reduction Δ 2  () still can be estimated by (12a) except for the following three modifications.
(i) The maximum (or minimum) of ∑  Δ 2  (), that is, Δ 2 emx (or Δ 2 emn ), should be found in the interior domain between the leftist and rightist observation points.
(ii) For the leftist (or rightist) observation that has only one adjacent observation to its right (or left) in the one-dimensional analysis domain, its error variance is still (iii) Note from (12a) that ∑  Δ 2  () → 0 and thus () → Δ 2 mn −Δ 2 emn as  moves outward far away from the leftist (or rightist) measurement point and thus far away from all the observations points.In this case, if Δ 2 mn − Δ 2 emn < 0 (as for the case in this section), then Δ 2  () estimated by () in (12a) may become unrealistically negative as  moves outward beyond the leftist (or rightist) measurement point, denoted by   .To avoid this problem, (12a) is modified into where  1 is a factor defined by It is easy to see from (12b) that for Δ 2 mn − Δ 2 emn < 0 and thus The analysis error variance is estimated by  2  () ≈  2   * () ≡  2  − Δ 2  () as in (7), except that Δ 2  () is computed by (12a) [or (12b)] for  within (or beyond) the interior domain.As shown by the example in Figure 8, the estimated  2  * () captures closely the spatial variation of the benchmark  2  () not only within but also beyond the interior domain.Using this estimated   * (), A  can be modified into A  , A  , or A  with its th element given by the same formulation as shown in (8a), (8b), or (8c).For the case in Figure 8, the benchmark A (not shown) has the same interior structure (for interior grid points  and ) as that for the case with periodic extension in Figure 6, but significant differences are seen in the following two aspects around the four corners (similar to those seen from Figures 7(a corners along the diagonal line (which is consistent with the increased analysis error variance toward the two ends of the analysis domain as shown in Figure 8 in comparison with that in Figure 5).(ii) The element value becomes virtually zero toward the two off-diagonal corners (because there is no periodic extension).The deviations of A  , A  , A  , and A  from the benchmark A are shown in Figures 9(a), 9(b), 9(c), and 9(d), respectively, for the case in Figure 8.As shown, the deviation becomes increasingly small when the estimated analysis error covariance matrix is modified successively to A  , A  , and A  .The REs of A  , A  , A  , and A  are listed in the first column of Table 3.As listed, the RE becomes increasingly small when A  is modified successively to A  , A  , and A  , which quantifies the successively reduced deviation shown in Figures 9(a)-9(d).

Analysis Error Variance Formulations for
Two-Dimensional Cases   in x, where Δ (or Δ) is the grid spacing in the  (or ) direction and Δ = Δ is assumed for simplicity.Since   (x) is a continuous function of x, the aforementioned formulation for the th diagonal element of A can be written into  2  (x) ≡  2  − Δ 2  (x) also as a continuous function of x, where is the error variance reduction produced by analyzing a single observation at x = x  .This reduction decreases rapidly and becomes much smaller than it peak value of    2   2  at x = x  as |x − x  | increases to  and beyond.

Uniform Coarse-Resolution Observations with Periodic
Extension.Consider that there are  coarse-resolution observations uniformly distributed in the above analysis domain of length   and width   with periodic extension along  and , so their resolution is Δ co ≡ (    ) 1/2 / 1/2 , where  =     ,   (or   ) denotes the number of observations uniformly distributed along the  (or ) direction in the two-dimensional analysis domain, and   /  =   /  is assumed (so Δ co =   /  =   /  ).In this case, as explained for the one-dimensional case in Section 2.2, the error variance reduction produced by each observation can be considered as an additional reduction to the reduction produced by its neighboring observations.This additional reduction is smaller than the reduction produced by a single observation, so the error variance reduction produced by analyzing the  coarse-resolution observations is bounded above by ∑  Δ 2  (x), which is similar to that for the onedimensional case in (4).For the same reason as explained for the one-dimensional case in (4), this implies that the domain-averaged value of ∑  Δ 2  (x) is larger than the true averaged reduction estimated by Δ 2  ≡  2  −  2  , where  2  is the domain-averaged analysis error variance estimated by the spectral formulation for two-dimensional cases in Section 2.3 of Xu et al. [8].
The domain-averaged value of ∑  Δ 2  (x) can be computed by where ∬  x denotes the integration over the twodimensional analysis domain, ∑  denotes the summation over  for the  grid points, and     =   Δ  Δ = ΔΔ is used in the last step.By extending  2  (x − x  ) with the analysis domain periodically in both the  and  directions, Δ 2   can be estimated analytically as follows: where ∬ x = ∬ denotes the integration over the entire space of x, is used in the second to last step, and  2 ≡ ∬x  (x)/ 2 is used with Δ 2 co ≡     / in the last step.For the double-Gaussian form of   (x) = 0.6 exp(−|x| 2 /2 2 ) + 0.4 exp(−2|x| 2 / 2 ) used in Section 4 of Xu et al. [8], we have  2 = 2(0.2+ 0.48/5).The derived value in (15b) is very close to the numerically computed value from (15a).
With the domain-averaged value adjusted from Δ 2  to Δ 2  , Δ 2  (x) can be estimated by the same formulation as in (6) except that  is replaced by x.The analysis error variance is then estimated by As shown by the example in Figure 10, the estimated  2  * (x) in ( 16) is very close to the benchmark  2  (x) computed precisely from (1b), and the deviation of  2  * (x) from the benchmark  2  (x) is within (−0.21, 0.35) m 2 s −2 .On the other hand, the constant analysis error variance ( 2  = 6.7 m 2 s −2 ) estimated by the spectral formulation deviates from the benchmark  2  (x) widely from −1.91 to 2.22 m 2 s −2 .
Using the estimated   * (x) in ( 16), the previously estimated analysis error covariance matrix, denoted by A  with its th element   ≡  2    (x  − x  ) obtained from the spectral formulation, can be modified into A  , A  , or A  with its th element given by the same formulation as in (8a), (8b), or (8c) except that   (or   ) is replaced by x  (or x  ).Again, as explained in Section 2.2 but for the two-dimensional case here, the accuracy of the second-step analysis depends on the accuracy of the estimated A over the extended nested domain, that is, the nested domain plus its extended vicinities within the distance of 2  outside the nested domain.Here,   is the decorrelation length scale of   (x) defined by  2   ≡ [−2  (x)/∇ 2   (x)]| x=0 according to (4.3.12) of Daley [12], and   (=4.52 km for the case in Figure 10) can be easily computed as a by-product from the spectral formulation.Over this extended nested domain, the relative error (RE) of each estimated A with respect to the benchmark A computed precisely from (1b) can be defined in the same way as that for A  in (9), except that the extended nested domain is two-dimensional here.The REs of A  , A  , A  , and A  computed for the case in Figure 10 are listed in the first column of Table 4.As listed, the RE becomes increasingly small when A  is modified successively to A  , A  , and A  .

Nonuniform Coarse-Resolution Observations with Periodic Extension.
Consider that the  coarse-resolution observations are now nonuniformly distributed in the analysis domain of length   and width   with periodic extension, so their averaged resolution can be defined by Δ co ≡ (    ) 1/2 / 1/2 .The spacing of a concerned coarseresolution observation, say the th observation, from its th adjacent observation (among the total 4 adjacent observations), can be denoted by Δ co .Now we can consider the limiting case of Δ co → 0 for  (≤4) adjacent observations with Δ co = Δ co for the remaining 4 −  (≥0) adjacent observations.In this case, the concerned th observation collapses onto the same point with its  adjacent observations.The  + 1 collapsed observations should be combined into one superobservation with a reduced error variance from  2  to  2  /( + 1).The error variance reduction produced by this superobservation still can be estimated by (14) but with On the other hand, without super-Obbing, the error variance reduction produced by the  + 1 collapsed observations will be overestimated by ( 14) with By comparing (17b) with (17a), it is easy to see that this overestimation can be corrected if the error variance is inflated from  2  to  2  =  2  +  2  for each of the ( + 1) collapsed observations.Based on the above analyses, when the error variance reduction produced by the concerned th observation is estimated by ( 14), the error variance should be adjusted for this observation unless Δ co = Δ co for  = 1, 2, 3, and 4. In particular,  2   can be adjusted to  2  =  2  +    2  with   given empirically by where ∑  denotes the summation over  for the four adjacent observations nearest to the concerned th observation.With   given by (18a), the adjusted  2  =  2  +    2  recovers not only the inflated observation error variance derived above for each limiting case [with Δ co → 0 for  = 1, 2, . . .,  (≤4) and Δ co = Δ co for the remaining 4 −  (≥0) observations] but also the original observation error variance  2   for uniformly distributed coarse-resolution observations.The above results suggest that   =    is the total number of adjacent observation points (nearest to x  ) used for estimating Δ omn (with  = 2) or Δ omx (or  = 4).By adjusting Δ 2 emx to Δ 2 mx and Δ 2 emn to Δ 2 mn , the error variance reduction can be estimated by where 16), except that Δ 2  (x) is computed by (19a).As shown by the example in Figure 11, the estimated  2  * (x) is fairly close to the benchmark  2  (x), and the deviation of  2  * (x) from the benchmark  2  (x) is within (−2.40, 4.20) m 2 s −2 .On the other hand, the constant analysis error variance ( 2  = 6.7 m 2 s −2 ) estimated by the spectral formulation deviates from the benchmark  2  (x) widely from −9.98 to 3.83 m 2 s −2 .Using this estimated   * (x), the previously estimated A  from the spectral formulation can be modified into A  , A  , or A  with its th element given by the same two-dimensional version of (8a), (8b), or (8c) as explained in Section 3.2.The REs of A  , A  , A  , and A  computed for the case in Figure 11 are listed in the first column of Table 5.As listed, the RE becomes increasingly small when A  is modified successively to A  , A  , and A  .

Nonuniform Coarse-Resolution Observations without Periodic Extension.
Consider that the  coarse-resolution observations are still nonuniformly distributed in the analysis domain of length   and width   but without periodic extension.In this case, their averaged resolution is still defined by Δ co ≡ (    /) 1/2 .To estimate their produced error variance reduction, we need to modify the formulations constructed in the previous subsection with the following preparations.First, we need to identify four near-corner observations among all the coarse-resolution observations.Each near-corner observation is defined as the one that nearest to one of the four corners of the analysis domain.Then, we need to identify   − 2 (or   − 2) near-boundary observations associated with each -boundary (or -boundary), where   (or   ) is estimated by the nearest integer to   /Δ co (or   /Δ co ).The total number of near-boundary observations is thus given by 2(  +   ) − 8. To identify these near-boundary observations, we need to divide the 2D domain uniformly along the -direction and -direction into     boxes, so there are 2(  +   ) − 8 boundary boxes (not including the four corner boxes).If a boundary box contains no coarse-resolution observation, then it is an empty box and should be substituted by its adjacent interior box (as a substituted boundary box).From each nonempty boundary box (including substituted boundary box), we can find one near-boundary observation that is nearest to the associated boundary.A closed loop of observation boundary can be constructed by piece-wise linear segments with every two neighboring near-boundary observation points connected by a linear segment and with each near-corner observation point connected by a linear segment to each of its two neighboring near-boundary observation points.
(ii) For each above defined near-boundary (or nearcorner) observation that has only three (or two) adjacent observations, its error variance is still adjusted from  2  to  2  +    2  but   is calculated by setting  2  (Δ co ) = 0 in (18a) for  = 4 (or  = 3 and 4).
(iii) Note from (19a) that ∑  Δ 2  (x) → 0 and thus Δ 2  (x) → Δ 2 mn − Δ 2 emn < 0 as x moves outward far away from all the observations points.In this case, if Δ 2 mn − Δ 2 emn < 0 (as for the case in this section), then Δ 2  (x) estimated by (19a) may become unrealistically negative as x moves outward beyond the above constructed observation boundary loop.To avoid this problem, (19a) is modified into where  2 is a factor defined by x  is the projection of x on the observation boundary loop and the projection from x is along the direction normal to the x-associated domain boundary (nearest to x).However, if x is closer to a corner observation point than the remaining part of the observation boundary loop, then x  is simply that near-corner observation point.It is easy to see from (19b) that, for Δ 2 mn − Δ 2 emn < 0 and thus The analysis error variance is estimated by  2  (x) ≈  2   * (x) ≡  2  − Δ 2  (x) as in (16), except that Δ 2  (x) is computed by (19a) [or (19b)] for x inside (or outside) the closed observation boundary loop.As shown by the example in Figure 12, the estimated  2  * (x) is fairly close to the benchmark  2  (x), and the deviation of  2  * (x) from the benchmark  2  (x) is within (−4.08, 5.54) m 2 s −2 .On the other hand, the constant analysis error variance ( 2  = 6.7 m 2 s −2 ) estimated by the spectral formulation deviates from the benchmark  2  (x) very widely from −16.1 to 3.82 m 2 s −2 .Using the estimated   * (x), the previously estimated A  from the spectral formulation can be modified into A  , A  , or A  with its th element given by the same two-dimensional version of (8a), (8b), or (8c) as explained in Section 3.2.The REs of A  , A  , A  , and A  computed for the case in Figure 12 are listed in the first column of Table 6.As listed, the RE becomes increasingly small when A  is modified successively to A  , A  , and A  .

Numerical Experiments for
One-Dimensional Cases 4.1.Experiment Design and Innovation Data.In this section, idealized one-dimensional experiments are designed and performed to examine to what extent the successively improved estimate of A in (8a), (8b), and (8c) can improve the two-step analysis.In particular, four types of two-step experiments, named TEe, TEa, TEb, and TEc, are designed for analyzing the high-resolution innovations in the second step with the background error covariance updated by A  , A  , A  , and A  , respectively, after the coarse-resolution innovations are analyzed in the first step.The TEe is similar to the first type of two-step experiment (named TEA) in Xu et al. [8], but the TEa, TEb, and TEc are new here.As in Xu et al. [8], a singlestep experiment, named SE, is also designed for analyzing all the innovations in a single step.In each of the above five types of experiments, the analysis increment is obtained by using the standard conjugate gradient descent algorithm to minimize the cost-function (formulated as in ( 7) of Xu et al. [8]) with the number of iterations limited to  = 20, 50, or 100 before the final convergence to mimic the computationally constrained situations in operational data assimilation.Three sets of simulated innovations are generated for the above five types of experiments.The first set consists of  (=10) uniformly distributed coarse-resolution innovations over the analysis domain (see Figure 1) with periodic extension and   (=74) high-resolution innovations in the nested domain of length /6 (similar to those shown by the purple × signs in Figure 1 of Xu et al. [8] but generated at the grid points not covered by the coarse-resolution innovations within the nested domain).The second (or third) set is the same as the first set except that the coarse-resolution innovations are nonuniformly distributed with (or without) periodic extension as shown in Figure 5 (or Figure 8).All the innovations are generated by simulated observation errors subtracting simulated background errors at observation locations.Observation errors are sampled from computer-generated uncorrelated Gaussian random numbers with   = 2.5 ms −1 for both coarse-resolution and high-resolution observations.Background errors are sampled from computer-generated spatially correlated Gaussian random fields with   = 5 ms −1 and   () modeled by the double-Gaussian form given in Section 2.2 (also see the caption of Figure 1).The coarseresolution innovations in the first, second, and third sets are thus generated in consistency with the three cases in Figures 1, 5, and 8, respectively.

Results from the First Set of Innovations.
The first set of innovations is used here to perform each of the five types of experiments with the number of iterations limited to  1, the TEe outperforms SE for  = 20, 50, and 100 but not for  increased to the final convergence.The improved performance of TEe over SE is similar to but less significant than that of TEA over SE in Table 1 of Xu et al. [8].The reduced improvement can be largely explained by the fact that the coarse-resolution innovations are generated here more sparsely and the deviation of A  from the benchmark A is thus increased (as seen from Figure 4 2 versus those from the SE.As shown, the TEe outperforms SE for  = 20 but not so for  = 50.The improvement of TEe over SE is similar to but much less significant than that of TEA over SE in Table 2 of Xu et al. [8].This reduced improvement can be largely explained by the fact that the coarse-resolution innovations are generated here not only more sparsely but also more nonuniformly than those in Section 3.3 of Xu et al. [8] and the deviation of A  from the benchmark A becomes much larger in Figure 7(a) here than that in Figure 7 3 versus those from the SE.As shown, the TEe outperforms SE for  = 20 but not so for  = 50.The improvement of TEe over SE is much less significant than that of TEA over SE in Table 3 of Xu et al. [8], and this reduced improvement can be explained by the same fact as stated for the previous case in Section 4.  In particular, the first set consists of  (= M x M y = 12 × 6) uniformly distributed coarse-resolution innovations over the periodic analysis domain (as shown in Figure 10) and   (=66) high-resolution innovations generated at the grid points not covered by the coarse-resolution innovations within the nested domain.The nested domain (  /6 = 20 km long and   /6 = 10 km wide) is the same as that shown in Figure 16 of Xu et al. [8].Again, all the innovations are generated by simulated observation errors subtracting simulated background errors at observation locations.Observation errors are sampled from computer-generated uncorrelated Gaussian random numbers with   = 2.5 ms −1 for both coarse-resolution and high-resolution observations.Background errors are sampled from computer-generated spatially correlated Gaussian random fields with   = 5 ms −1 and   (x) modeled by the double-Gaussian form given in Section 3.2 (also see the caption of Figure 10).The second (or third) set is the same as the first set except that the coarse-resolution innovations are nonuniformly distributed with (or without) periodic extension as shown in Figure 11 (or Figure 12).

Conclusions
In this paper, the two-step variational method developed in Xu et al. [8] for analyzing observations of different spatial resolutions is improved by considering and efficiently estimating the spatial variation of analysis error variance produced by analyzing coarse-resolution observations in the first step.The constant analysis error variance computed from the spectral formulations in Xu et al. [8] can represent the spatial averaged value of the true analysis error variance but it cannot capture the spatial variation in the true analysis error variance.As revealed by the examples presented in this paper (see Figures 1,2,5, and 8 for one-dimensional cases and Figures 10-12 for two-dimensional cases), the true analysis error variance tends to have increasingly large spatial variations when the coarse-resolution observations become increasingly nonuniform and/or sparse, and this is especially true and serious when the separation distances between neighboring coarse-resolution observations become close to or even locally larger than the background error decorrelation length scale.In this case, the spatial variation of analysis error variance and associated spatial variation in analysis error covariance need to be considered and estimated efficiently in order to further improve the two-step analysis.The analysis error variance can be viewed equivalently and conveniently as the background error variance minus the total error variance reduction produced by analyzing all the coarse-resolution observations.To efficiently estimate the latter, analytic formulations are constructed for three types of coarse-resolution observations in one-and two-dimensional spaces with successively increased complexity and generality.The main results and major findings are summarized below for each type of coarse-resolution observations: (i) The first type consists of uniformly distributed coarseresolution observations with periodic extension.For this simplest type, the total error variance reduction is estimated in two steps.First, the error variance reduction produced by analyzing each coarse-resolution observation as a single observation is equally weighted and combined into the total.Then, the combined total error variance reduction is adjusted by a constant to match to the domain-averaged total error variance reduction estimated by the spectral formulation [see (5a), (5b), (15a), and (15b)].The estimated analysis error variance (i.e., the background error variance minus the adjusted total error variance reduction) captures not only the domain-averaged value but also the spatial variation of the benchmark truth (see Figures 1, 2, and 10).
(ii) The second type consists of nonuniformly distributed coarse-resolution observations with periodic extension.For this more general type, the total error variance reduction is also estimated in two steps: The first step is similar to that for the first type but the combination into the total is weighted based on the averaged spacing of each concerned observation from its neighboring observations [see (11a), (11b), (18a), and (18b)].In the second step, the combined total error variance reduction is adjusted and scaled to match the maximum and minimum of the true total error variance reduction estimated from the spectral formulation for uniformly distributed coarse-resolution observations but with the observation resolutions set, respectively, to the minimum spacing and maximum spacing of the nonuniformly distributed coarse-resolution observations [see (12a) and (19a)].The estimated analysis error variance captures not only the maximum and minimum but also the spatial variation of the benchmark truth (see Figures 5 and 11).
(iii) The third type consists of nonuniformly distributed coarse-resolution observations without periodic extension.For this most general type, the total error variance reduction is estimated with the same two steps as for the second type, except that three modifications are made to improve the estimation near and at the domain boundaries [see (i)-(iii) in Sections 2.4 and 3.4].The analysis error variance finally estimated captures the spatial variation of the benchmark truth not only in the interior domain but also near and at the domain boundaries (see Figures 8 and 12).
The above estimated spatially varying analysis error variance is used to modify the analysis error covariance computed from the spectral formulations of Xu et al. [8] in three different forms [see (8a), (8b), and (8c)].The first is a conventional formulation in which the covariance is modulated by the spatially varying standard deviation separately via each entry of the covariance to retain the self-adjointness.This modulation has a chessboard structure

Figure 3 :
Figure 3: Structure of benchmark A plotted by color contours every 1 m 2 s −2 for the case in Figure 1.

Figure 4 :
Figure 4: (a) Deviation of A  from benchmark A in Figure3plotted by color contours every 0.5 m 2 s −2 .Deviations of A  , A  , and A  from benchmark A are plotted by color contours every 0.2 m 2 s −2 in panels (b), (c), and (d), respectively.Here, A  is the previously estimated analysis error covariance matrix with its th element   ≡ 2     (  −   ) obtained from the spectral formulation, while A  , A  , and A  are the newly modified estimates of A as shown in (8a), (8b), and (8c), respectively.

Figure 7 :
Figure 7: As in Figure 4 but for the case in Figure 5.

Figure 8 :
Figure 8: As in Figure 5 but without periodic extension.

Figure 9 :
Figure 9: As in Figure 7 but for the case in Figure 8.

Figure 11 :
Figure 11: As in Figure 10 but for the second set of innovations with nonuniformly distributed coarse-resolution observations, and the colored contours are plotted every 1 m 2 s −2 for the deviation of  2  * (x) from  2  (x) in panel (b).

Figure 12 :
Figure 12: As in Figure 11 but without periodic extension.
(a) in comparison with Figure 5(b) of Xu et al. [8]).The TEa outperforms TEe for  = 20 and 50 before  increased to 100 (which is very close to the final convergence at  = 116 for TEa).The improvement of TEa over TEe is consistent with and can be largely explained by the improved accuracy of A  [RE(A  ) = 0.156] over A  [RE(A  ) = 0.229].The TEb outperforms TEa for  = 20 and 50 (before the final convergence at  = 67).The improvement of TEb over TEa is consistent with the improved accuracy of A  [RE(A  ) = 0.101] over A  .The TEc outperforms TEb for each listed value of , and the improvement is consistent with the improved accuracy of A  [with RE(A  ) = 0.042] over A  .4.3.Results from the Second Set of Innovations.The second set of innovations is used here to perform each of the five types of experiments with the number of iterations limited to  = 20, 50, or 100 before the final convergence.The domain-averaged RMS errors of the analysis increments obtained from the four two-step experiments are shown in Table (b) of Xu et al. [8].The TEa outperforms TEe for  = 20 and 50 but still underperforms SE for  increased to 50 and beyond.The improvement of TEa over TEe is consistent with the improved accuracy of A  [RE(A  ) = 0.238] over A  [RE(A  ) = 0.355].The TEb outperforms TEa for each listed value of  and also outperforms SE for  up to 100.The improvement of TEb over TEa is consistent with the improved accuracy of A  [RE(A  ) = 0.197] over A  .The TEc outperforms TEb for each listed value of , and the improvement is consistent with the improved accuracy of A  [RE(A  ) = 0.148] over A  .4.4.Results from the Third Set of Innovations.The third set of innovations is used here to perform each of the five types of experiments with the number of iterations limited to  = 20, 50, or 100 before the final convergence.The domain-averaged RMS errors of the analysis increments obtained from the four two-step experiments are shown in Table 3. The TEa outperforms TEe for  = 20 and 50, and the improvement is consistent with the improved accuracy of A  [RE(A  ) = 0.238] over A  [RE(A  ) = 0.355].The TEb outperforms TEa for each listed value of , which is consistent with the improved accuracy of A  [RE(A  ) = 0.196] over A  .The TEc outperforms TEb for each listed value of , which is consistent with the improved accuracy of A  [RE(A  ) = 0.147] over A  .
).The analytically derived value in (5b) is very close to (slightly larger than) the numerically computed value from (5a).With the domain-averaged value of ∑  Δ 2 () adjusted from Δ 2  to Δ 2  , Δ 2  () can be estimated by Δ co decreased to Δ omn (or increased to Δ omx ), where Δ omn (or Δ omx ) is the minimum (or maximum) spacing of adjacent observations among all nonuniformly distributed coarse-resolution observations in the two-dimension analysis domain.Specifically, Δ omn (or Δ omx ) is estimated by min  (∑  |x  − x  |)/ with  = 2 and Δ omx is estimated by max  (∑  |x  − x  |)/ with  = 4, where x  denotes the th observation point, x  denotes the observation point that is th nearest to x  , min m (or max m ) denotes the minimum (or maximum) over index  for all the coarseresolution observation points in the two-dimension analysis domain, ∑  denotes the summation over  from 1 to , and

Table 2 :
As in Table 1 but for the second set of innovations with periodic extension.

Table 3 :
As in Table2but for the third set of innovations without periodic extension.50,or100before the final convergence.The accuracy of the analysis increment obtained from each experiment with each limited  is measured by its domain-averaged RMS error (called RMS error for short hereafter) with respect to the benchmark analysis increment computed precisely from (1a).Table1lists the RMS errors of the analysis increments obtained from the SE, TEe, TEa, TEb, and TEc with the number of iterations increased from  = 20 to 50, 100, and/or the final convergence.As shown in Table

Table 4 :
As in Table1but for the two-dimensional case in Figure10in which the first set of two-dimensional innovations is used with periodic extension.

Table 5 :
As in Table4but for the two-dimensional case in Figure11where the second set of innovations is used with periodic extension.
5.1.Experim11,and 12, and Innovation Data.In this section, idealized two-dimensional experiments are designed and named similarly to those in Section 4 except that simulated innovations are generated in three sets for the twodimensional cases inFigures 10,11,and 12, respectively.
[8]st Set of Innovations.The first set of innovations is used here to perform each of the five types of experiments with the number of iterations limited to  = 20, 50, or 100 before the final convergence.The domain-averaged RMS errors of the analysis increments obtained from the four two-step experiments are shown in Table4versus those from the SE.As shown, the TEe outperforms SE for each listed value of  before the final convergence, which is similar to the improved performance of TEA over SE shown in Table4of Xu et al.[8].The TEa outperforms TEe as  increases to 100 and beyond, which is consistent with the improved accuracy of A The second set of innovations is used here to perform each of the five types of experiments with the number of iterations limited to  = 20, 50, or 100 before the final convergence.The domainaveraged RMS errors of the analysis increments obtained from the four two-step experiments are shown in Table5versus those from the SE.As shown, the TEe outperforms SE for each listed value of  before the final convergence.The TEa outperforms TEe slightly, and the improved performance is consistent with the improved accuracy of A  [RE(A  ) = 0.274] over A  [RE(A  ) = 0.462].The TEb outperforms TEA for each listed value of , which is consistent with the improved accuracy of A  [RE(A  ) = 0.244] over A  .The TEc outperforms TEb for  > 20, and the improved performance is consistent with the improved accuracy of A  [RE(A  ) = 0.165] over A  .
[RE(A  ) = 0.181] over A  [RE(A  ) = 0.233].The TEb outperforms TEa as  increases to 50 and beyond, which is consistent with the improved accuracy of A  [RE(A  ) = 0.130] over A  .The TEc outperforms TEb for each listed value of , which is consistent with the improved accuracy of A  [RE(A  ) = 0.038] over A  .5.3.Results from the Second Set of Innovations.50, or 100 before the final convergence.The domain-averaged RMS errors of the analysis increments obtained from the four two-step experiments are shown in Table6versus those from the SE.As shown, the TEe outperforms SE for each listed value of  before the final convergence.The improved performance of TEe over SE is similar to but less significant

Table 6 :
[8]in Table5but for the two-dimensional case in Figure12where the third set of innovations is used without periodic extension.that of TEA over SE in Table5of Xu et al.[8], and the reason is mainly because the coarse-resolution innovations are generated more sparsely and nonuniformly than those in Section 4.3 of Xu et al.[8].The TEa outperforms TEe for  = 20 but not so as  increases to 50 and beyond, although A  has an improved accuracy [RE(A  ) = 0.305] over A  [RE(A  ) = 0.462].The TEb outperforms TEa for each listed value of , and the improved performance is consistent with the improved accuracy of A  [RE(A  ) = 0.258] over A  .The TEc outperforms TEb for each listed value of , which is consistent with the improved accuracy of A  [RE(A  ) = 0.240] over A  . than