## 1. Introduction

The flow-dependent background error covariance is an advantage of the ensemble Kalman filter (EnKF; Evensen 2007) over variational analysis methods. However, a commonly used finite model ensemble often underestimates prior variance and introduces spurious correlation between a state variable and remote observations. While various variance inflation schemes (e.g., Anderson and Anderson 1999; F. Zhang et al. 2004; Anderson 2008; Houtekamer et al. 2009; Miyoshi 2011) have been developed to address the first issue, a localization technique is usually used to ease the second problem. The original localization scheme (e.g., Houtekamer and Mitchell 1998; Anderson 2001; Hamill et al. 2001; Szunyogh et al. 2008) was the fixed localization approach, which is usually realized by a Schur product (an element-by-element multiplication) of the ensemble-evaluated error covariance with an analytic localization operator. All fixed localization models need to determine an impact radius that defines a maximum impact range of observations and significantly influences the analysis quality of the EnKFs. The optimal impact radius depends on the ensemble size, the observing system, and the numerical model itself (Houtekamer and Mitchell 1998; Mitchell et al. 2002). However, it is time consuming to tune the optimal impact radius, especially for a general circulation model (GCM).

Given the above disadvantages of fixed localization models, many adaptive localization algorithms (e.g., Anderson 2007; Bishop and Hodyss 2007, 2009a,b, 2011; Bishop et al. 2011; Jun et al. 2011; Anderson and Lei 2013; Lei and Anderson 2014) have been proposed. Several adaptive schemes (e.g., Anderson 2007) entail trivial additional cost to a standard filter and are routinely used in GCMs and other very large models. A more detailed review of these localization methods can be found in Berre and Desroziers (2010). Because of the low computational cost and feasible implementation, the fixed method is commonly used in actual data assimilation. However, the assimilation quality of EnKF with a fixed localization strongly depends on the impact radius. An overly small (large) impact radius may cause the loss (contamination) of longwave information (analysis states). Accordingly, some multiscale analysis methods have been presented to extract the multiscales of the observations. Miyoshi and Kondo (2013) proposed a dual-localization method, which was demonstrated to be superior to a single-localization method under twin experiments. Wu et al. (2014, hereafter W14) utilized a multiscale analysis (MSA) technique to extract the multiscale information from the observational residuals (the difference between the observations and the interpolated analysis ensemble means produced by the EnKF without inflation) caused by extreme (that is overly small or overly large) impact radii and then corrected the analysis ensemble mean of the EnKF without inflation. Their results demonstrated that the EnKF-MSA method can greatly improve the assimilation quality of the EnKF without inflation for extreme impact radii and has less dependence on the impact radius than the EnKF without inflation. However, the computational cost of the MSA will be huge for high-resolution numerical models and intensive observing systems, and the EnKF-MSA method made things worse than the EnKF without inflation for moderate impact radii.

In this study, we present an adaptive approach to make up for the above demerits. At each analysis step, after the EnKF without inflation is completed, if the root-mean-square error (RMSE) between the interpolated analysis ensemble mean and observations exceeds a threshold, a multigrid analysis (MGA) is adaptively activated to retrieve the multiscale information from the observational residual; then the analysis field of the MGA is added to the analysis ensemble mean of the EnKF without inflation to form the final analysis ensemble mean. With a global barotropic spectral model and biased twin experiments, as well as an idealized observing system, the performance of the proposed algorithm is investigated. Although localization can also be applied to the background error covariance in the observation space (e.g., Houtekamer and Mitchell 1998; Ott et al. 2004; Hunt et al. 2007; Greybush et al. 2011), we only focus on the background error covariance localization in the state space in this study. Note that the EnKF-MSA and the adaptive EnKF-MGA are different from the running in place and quasi-outer-loop (RIP/QOL) method of Yang et al. (2012) in some aspects. The RIP/QOL method proposed activation criteria to determine when the second loop is activated to use the same observations more than once. In contrast, the MSA and the MGA are applied only to the observational residuals rather than the original observations. There is some similarity between the adaptive EnKF-MGA and the RIP/QOL: that is, both methods use criteria to decide whether to enter the next loop; in contrast, the EnKF-MSA triggers the MSA at each analysis step without any judgment.

The remainder of this paper is organized as follows. Section 2 briefly describes the global barotropic spectral model and the data assimilation algorithms used. Section 3 introduces biased twin experiments. Detailed performance of the proposed method is investigated in section 4. Summary and discussion are given in section 5.

## 2. Methodology

### a. The model

*f*represent the relative vorticity and planetary vorticity (i.e., Coriolis parameter), respectively, and

*H*is the depth of the atmospheric layer.

*f*denotes the Coriolis parameter,

*y*represents the northward meridional distance from equator, and

A rhomboidal 21 truncation is applied to the transformation between spectral coefficients and grid values. For the time stepping of the model, the state variables are spectral coefficients. The integration step size is half an hour. A leapfrog time step is used to integrate the model, and a Robert–Asselin time filter (Robert 1969; Asselin 1972) with the filter coefficient *γ* is applied to damp spurious computational modes. For the data assimilation, the state variable is the atmospheric streamfunction at the 64 (longitude) × 54 (latitude) Gaussian grid points.

### b. The EAKF algorithm in Anderson (2001)

*y*

^{o}, the implementation of EAKF can be compressed to the following two steps. First, compute the observational increment as follows:where

*i*indexes ensemble member;

*y*

^{o}; and

*r*and

*y*

^{o}, respectively. The

*i*th prior ensemble member of

*y*

^{o},

*i*th prior ensemble member of state variable

*x*. Second, use the following linear regression to project the observational increment to the model grids:where

*j*indexes state variable, while

*x*

_{j}and

*y*

^{o}, and Δ

*x*represents the state increment.

*y*

^{o}and

*x*

_{j}only appears in the linear regression formula Eq. (4). Therefore, when the covariance localization is introduced into the EAKF, the local support correlation is imposed in the numerator of the coefficient of linear regression as follows:where

*ρ*

_{j,y}represents the localization function between

*y*

^{o}and

*x*

_{j}.

*b*denotes the physical distance between

*y*

^{o}and

*x*

_{j}, and

*a*represents the half-width of the GC function. In this study, we use

*a*to denote the impact radius of observations for the EnKFs. Table 1 lists the key notations used in this study.

Glossary of notations in this study.

Because of the existence of model error, model nonlinearity, sampling error, and so on, the prior variance of the model state is underestimated, which may cause the filter to diverge. Thus, inflation in some form (normally explicit) is used universally in large applications. In this study, we use “the standard EnKF” to represent the EnKF with inflation.

### c. Multiscale analysis

**y**

^{o},

*K*× 1,

*K*×

*M*, and

*M*× 1, respectively.

*l*th scale level in the MSA iswhere

*L*

_{MSA}is the number of scale levels, and the superscript “T” indicates the transpose of the matrix;

**x**

_{MSA}, the linear projection operator from the observation space to the state space, the observation error covariance matrix in the state space, the smoothing operator, and the observational innovation vector for the

*l*th scale level in the MSA. The dimensions of the above five matrices are

*M*× 1,

*M*×

*K*,

*M*×

*M*,

*M*×

*M*, and

*K*× 1, respectively.

*l*th scale level. For each scale level, the dimension of

*K*×

*M*. The mapping operator

**y**

^{res}in Eq. (7)],

**y**

^{res}. Thus,

*l*− 1)th scale level. The observational error covariance matrix

**x**produced by the MSA is

Keeping the analysis grid as the model grid for all scale levels, the MSA extracts multiscale information from the observational residual by varying the correlation factor, which is involved in *i*, *j*)th element of *L*_{ij}^{(l)}, denotes the weight of the *j*th observation on the *i*th state variable [see Eq. (27) in W14] for the *l*th scale level. The numerator of *L*_{ij}^{(l)} is the GC localization function *a*^{(l)} represents the impact radius (i.e., the GC half-width) of observations for the *l*th scale level in the MSA; *b*_{ij} represents the physical distance between the *j*th observation and the *i*th state variable; and the denominator of *L*_{ij}^{(l)} is a normalization factor.

*l*th scale level, the computational cost of the MSA consists of the calculations of

*MK*and

*MK*, respectively. For each iteration of the L-BFGS algorithm, the computational costs of the cost function [Eq. (8)] and the gradient [Eq. (10)] are

*M*

^{2}+

*MK*+ 2

*M*+ 2 and

*M*

^{2}+

*MK*, respectively. Thus, the computational cost of the MSA for the

*l*th scale level iswhere

*T*

_{MSA}represents the number of iterations of the L-BFGS algorithm in the MSA. The total computational cost of the MSA is

From the above analysis, we can see that the computational cost of the MSA is mainly caused by the large dimensions of the analysis grid for each scale level and by the mapping from the observation space to the state space. It can be expected that the computational cost of the MSA will be very large for high-resolution and/or high-dimension numerical models and intensive observing systems.

### d. The EnKF-MSA method

W14 presented the EnKF-MSA method, the implementation of which at each analysis step can be expressed as the following five steps: First, use the EnKF without inflation to assimilate all available observations with *a*-GC half-width (i.e., the GC half-width in the EAKF, see *a* in Table 1). Second, compute the observational residual using Eq. (7). Third, use the MSA to extract multiscale information from the observational residual with a prescribed configuration of *a*^{(l)}. Fourth, add the analysis field of the MSA to the analysis ensemble mean produced by the EnKF without inflation to obtain the final ensemble mean. Last, add the new ensemble mean to the analysis ensemble perturbations produced by the EnKF without inflation to generate the final ensemble members. More detail on the EnKF-MSA method can be found in W14.

### e. Limitations of the EnKF-MSA method

Since the EnKF-MSA activates the MSA at each analysis step, the limit on the high computational cost of the MSA for high-resolution and/or high-dimension numerical models and intensive observing systems is also shared by the EnKF-MSA.

In addition, as we stated in the introduction, although the EnKF-MSA method can greatly improve the performance of the EnKF without inflation for extreme (that is, overly small or overly large) impact radii, it made things worse than the EnKF without inflation for moderate impact radii (see Fig. 1 in W14).

**x**

^{t}represents the truth of model state. The first term on the RHS of Eq. (13) represents the unknown true observational error, while the second term indicates the unknown true analysis error that is projected to the observation space. In other words, the above two terms represent the noise and the signal, respectively. Under specific configurations of ensemble size, model error, and observing system, the analysis error may be smaller than the observational error for moderate impact radii, which is what happened in W14. Under this circumstance, the observational residual is dominated by the noise. It is difficult for the MSA to extract the multiscale information of the signal, because the EnKF-MSA method activates the MSA at each analysis step without checking the signal-to-noise ratio (SNR) of the observational residual.

**y**

^{res}as follows:where

*k*indexes the observations;

**h**

_{k}represents the

*k*th row vector of

**y**

^{res}so as to determine whether the multiscale information should be extracted, the performance of the EnKF-MSA method may be further improved for moderate impact radii.

In summary, there are mainly two limitations of the EnKF-MSA method. One is the high computational cost of the MSA. The other is the inferior performance to the EnKF without inflation for moderate impact radii.

### f. Multigrid analysis

To ease the first limitation of the EnKF-MSA, an economical algorithm that can realize similar functions of the MSA is needed. Inspired by this idea, we introduce an MGA method that was initially suggested for solving differential equations (Briggs et al. 2000) and later introduced into the data assimilation community (e.g., Li et al. 2008, 2010; Xie et al. 2011).

*l*th grid in the MGA is formulated aswhere

*L*

_{MGA}is the number of levels of analysis grids. The expressions

**x**

_{MGA}, the mapping operator from the state space to the observation space, the observational error covariance matrix, and the observational innovation vector and the smoothing matrix for the

*l*th grid in the MGA. The dimensions of the above five matrices are, respectively,

*m*

_{l}× 1,

*K*×

*m*

_{l},

*K*×

*K*,

*K*× 1, and

*m*

_{l}×

*m*

_{l}, where

*m*

_{l}represents the dimension of the

*l*th grid. Note that, to retrieve multiple signals from longwave to shortwave,

*m*

_{l}monotonically increases from

*m*

_{1}to

*l*increases from 1 to

*L*

_{MGA}. Here,

*m*

_{l}is set to (2

^{l−1}+ 1) × (2

^{l−1}+ 1) = 4

^{l−1}+ 2

^{l}+ 1: that is, the resolution of the

*l*th analysis grid is approximately twice that of the (

*l*− 1)th grid.

*l*th analysis grid to the observation space;

*l*− 1)th grid level.

*i*,

*j*)th element of

*L*

_{MGA}th grid) with the operator

**x**produced by the MGA iswhere

*l*th grid level, the computational cost of the MGA consists of the calculation of

*Km*

_{l−1}. For each iteration of the L-BFGS algorithm, the computational costs of the cost function [Eq. (15)] and the gradient [Eq. (18)] are

*l*th grid level iswhere

*T*

_{MGA}represents the number of iterations of the L-BFGS algorithm in the MGA. The total computational cost of the MGA is

Here, we compare the computational cost of the MGA to that of the MSA. For small *l*, *m*_{l} is much smaller than *M*. If *T*_{MGA} equals *T*_{MSA}, *L*_{MGA} usually has a similar magnitude to that of *L*_{MSA}, the computational cost of the MGA is also much smaller than that of the MSA. Thus, the MGA is expected to be much more computationally effective than the MSA.

Since the MGA has similar functions as the MSA, a direct way to ease the first limitation of the EnKF-MSA is to replace the MSA by the MGA so as to put forward the EnKF-MGA method.

### g. An adaptive EnKF-MGA method

To further avoid the second limitation of the EnKF-MSA method, we present an adaptive EnKF-MGA method based on the above EnKF-MGA method. At each analysis step, the adaptive approach includes the following sequential steps:

- Sequentially adjust the ensemble members using the EnKF without inflation with
*a*-GC half-width. - Project the analysis ensemble mean produced by the EnKF without inflation to the observation positions to get the EnKF-estimated posterior observations. Then the observational residual is computed by Eq. (7).
- Compute the RMSE between the original observations and the EnKF-estimated posterior observations. If the RMSE is larger than a critical value (denoted as
*θ*), we go through.- 3.1) Retrieve the multiscale information from the observational residual using the MGA method.
- 3.2) Add the analysis field generated by the MGA to the ensemble mean produced by the EnKF without inflation to obtain the final analysis ensemble mean.
- 3.3) Add the final analysis ensemble mean to the ensemble perturbations produced by the EnKF without inflation to generate the final ensemble members.

If the RMSE is equal to or smaller than *θ*, the analysis of the EnKF without inflation will be taken as that of the adaptive EnKF-MGA method: that is, the adaptive method reduces to the EnKF without inflation.

Note that no inflation is applied to the first EnKF step in the adaptive EnKF-MGA. There are two motivations: one is to facilitate the comparison between the EnKF-MSA and the adaptive EnKF-MGA; the other is to investigate whether the adaptive EnKF-MGA without inflation can produce comparable (or better) assimilation results to (than) the standard EnKF, which is meaningful for the future application of the adaptive EnKF-MGA in GCMs. In addition, similar to the implementation of the EnKF-MSA, the adaptive EnKF-MGA only updates the ensemble mean produced by the first EnKF step but keeps on using the same ensemble perturbations from the first EnKF step. The unadjusted ensemble perturbations or covariance will introduce an inconsistency that is common to the RIP/QOL method and will be discussed in section 5 in detail.

Since the threshold *θ* serves as the on–off switch of the MGA, it is a key parameter of the adaptive EnKF-MGA method. Here, we use the following test to determine the value of *θ* given a significance level *α*:

- Null hypothesis: All elements of
**y**^{res}are dominated by the observational error. - Alternative hypothesis: At least one element of
**y**^{res}is not dominated by the observational error.

The null hypothesis indicates that each element of **y**^{res} (denoted as *r*, respectively. With a commonly made assumption that observational errors are uncorrelated, the sum of the square of the normalized observational errors [i.e., *χ*^{2} probability distribution function with *K* degrees of freedom. If the real observed value of the statistic *α*. Here, *α*) upper fractile of the *χ*^{2}(*K*) distribution.

^{res}to denote the RMSE between

*θ*is accordingly set to the RHS of the above inequality. If the RMSE

^{res}is greater than

*θ*, the null hypothesis will be denied with the significance level

*α*and the MGA will be activated to extract the multiscale information from

**y**

^{res}. Obviously, the unit of

*θ*is the same as that of the observation (i.e., m

^{2}s

^{−1}in this study). Note that, because there are no Gaussian and uncorrelation assumptions of the analysis errors in the observation space, the analysis error is not involved in the determination of

*θ*.

According to Eq. (23), *θ* is a function of the significance level *α*, the number of observations *K*, and the standard deviation of observational error *r*. Since *r* is usually assumed to be a constant (10^{6} m^{2} s^{−1} in this study), *θ* mainly depends on *K* and *α*. Figure 1 displays the variation of *θ* with respect to *K* for *α =* 0.01 (red), 0.05 (black), and 0.10 (blue). Obviously, *θ* varies inversely with both *K* and *α*. The sensitivity mostly occurs for small values of *K*. When *K* exceeds a critical value, *θ* gradually saturates to a constant for each value of *α*. This justifies that *θ* can be set to a constant value when *K* is sufficiently large (e.g., 10^{7} in practice).

## 3. Biased twin experiments setup

Biased twin experiments are designed to investigate the performances of the standard EnKF, the EnKF-MSA, and the adaptive EnKF-MGA methods. The model error is assumed to arise from the uncertainty of the time filter coefficient *γ*. We assume *γ* = 0.01 in the truth model and *γ* = 0.02 in the assimilation model. The sole difference produces an apparent bias between the assimilation model and the truth model.

Starting from the streamfunction at 1200 UTC 1 January 1991 derived from the European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis 500-hPa *u* and *υ* data, both the truth model and the assimilation model are spun up for 30 model days. Then the truth model is further integrated for another 200 model days to sample “observations.” The observational interval is set to 6 h (12 time steps). A Gaussian noise with the standard deviation of 10^{6} m^{2} s^{−1} is imposed to the “truth” streamfunction to simulate the observational error. Here, we set the standard deviation of the observational error to be 4% of the global mean (2.5 × 10^{7} m^{2} s^{−1}) of the natural variability (i.e., climatological standard deviation) of the streamfunction. For the observing system, we assume observations are randomly and uniformly distributed in area A (0°–180°, 0°–90°N) and in area B (180°–360°, 0°–90°N) and the Southern Hemisphere (denoted as area C). To reflect the spatial distribution of the observing system, we assume the spatial sampling density of the observations in area A is twice that in area B and three times that in area C. In area A, the number of observations is 864, which equals the number of model grids in area A. Figure 2 displays the observing system (dots) and the model grids (pluses). The values of *K* and *M* in this study are 1872 and 3456, respectively.

Initial model ensembles of atmospheric streamfunction for the assimilation are generated by adding a Gaussian white noise with the mean being zero and the same standard deviation of observational error to the biased initial condition generated by the assimilation model. Note that, because of the coarse resolution of the barotropic model, the Gaussian noise is independently added to each model grid without any coherent length scale. The ensemble size is set to a typical value of 20. Additionally, because the leapfrog scheme is used to integrate the model, a two-time-level adjustment method (S. Zhang et al. 2004) is applied to the data assimilation: namely, the observations at time *t* are used to adjust the model states at times *t* and *t* − 1.

Four types of experiment are conducted to evaluate the performances of the three assimilation algorithms. The first is the ensemble free run (without observational constraints), serving as the reference experiment (CTL). The second is the standard EnKF implemented by the EAKF. The third is the EnKF-MSA method, and the last one is the adaptive EnKF-MGA method. Note that the EnKF-MSA and adaptive EnKF-MGA methods do not employ inflation in this study. In addition, because the EnKF-MSA and the adaptive EnKF-MGA are based on the EnKF without inflation, we also conduct the experiment of the EnKF without inflation so as to compare the performance of the adaptive EnKF-MGA to that of the EnKF without inflation, although an EnKF without inflation simply is not an interesting case for comparison in reality. All assimilation schemes are applied to global observational residuals and activated at the initial time. The three data assimilation algorithms use the same observing system and ensemble initial conditions. The number of iterations of the optimization algorithm for both the MGA and the MSA is set to 10 (i.e., *T*_{MSA} = *T*_{MGA} = 10). For the adaptive EnKF-MGA method, to determine the threshold *θ* we set the significance level *α* to the typical value 0.01. With the above observing system, *θ* can be computed as 1.04 × 10^{6} m^{2} s^{−1}. Discarding the assimilation results in the first 100 days as the spinup, the results of the last 100 days are used to conduct error statistics and analysis. Note that, because of the overly large error of

To investigate the dependence of the adaptive EnKF-MGA method on *a*, 16 values of *a* starting from 250 to 4000 km with an even increment of 250 km are used. According to the description of the EnKF-MSA method in W14, the configuration of *L*_{MSA} equals 5.

*υ*represents the inflation factor, and

*i*th ensemble member for the

*j*th model state. The bar stands for the ensemble mean. For each assimilation experiment of the EAKF, the best value of

*υ*is determined by minimizing the time-mean RMSE through varying

*υ*from 1.00 to 3.00 with an even increment of 0.01. The time-mean RMSE is defined aswhere

*S*represents the number of analysis steps,

*s*indexes the analysis step, and RMSE

_{s}denotes the RMSE at the

*s*th analysis step;

*im*and

*jm*are the dimensions of zonal and meridional model grids (which are 64 and 54), respectively. The ensemble mean of the streamfunction

*t*” denote the prior and the truth, respectively. Note that, because of the coarse resolution of the model, we can use this simple method to determine the best

*υ*. In GCM applications, some adaptive methods (e.g., Anderson 2008) are usually used. The black curve in Fig. 3 shows the variation of

*υ*for the standard EnKF. For

*a*= 250 km,

*υ*is 1.0 (i.e., no need to apply the variance inflation). As the value of

*a*increases,

*υ*gradually increases. Note that, when

*a*equals 4000 km, the EnKF crashes for all values of

*υ*varying from 1 to 3. For values of

*υ*larger than 3, it also fails. This may be caused by the spatially and temporally uniform inflation factor.

## 4. Results

### a. Determination of the number of grid levels

There is a factor that may complicate the application of the MGA. It is the value of *L*_{MGA}. Thus, a question arises as follows: How does one determine *L*_{MGA}?

To answer this question, we take *a* = 250 km as an example to investigate the dependence of the adaptive EnKF-MGA method on *L*_{MGA}. Since the dimension of the first analysis grid level in the MGA is usually set to 2 × 2, the value of *L*_{MGA} depends on the selection of *m*_{l} = (2^{l−1} + 1) × (2^{l−1} + 1). We perform 11 experiments of the adaptive EnKF-MGA with *L*_{MGA} varying from 0 to 10. When *L*_{MGA} = 0, the adaptive EnKF-MGA reduces to the EnKF without inflation, which is also the standard EnKF because the optimal *υ* = 1.0 (the black curve in Fig. 3).

^{6}m

^{2}s

^{−1}) of the prior ensemble mean of

*L*

_{MGA}. The vertical line at each

*L*

_{MGA}value represents the ±

*ζ*bound of the time-mean RMSE, where

*ζ*represents the standard deviation of the RMSE:The optimal value of

*L*

_{MGA}is around 9. However, the RMSE seems to have virtually saturated at

*L*

_{MGA}= 7, which has an RMSE close to the optimal value but at less cost. The standard deviations of the RMSE for

*L*

_{MGA}= 7, 8, 9, and 10 are smaller than those for

*L*

_{MGA}= 0, 1, 2, 3, and 4, demonstrating that the adaptive EnKF-MGA is more stable for the former

*L*

_{MGA}values. This can also be justified by Fig. 4b, which plots the RMSE (10

^{6}m

^{2}s

^{−1}) time series of the prior ensemble mean of

*L*

_{MGA}= 7, the dimension of the finest analysis grid is 65 × 65, which is close to the dimensions of the model grid (64 × 54). Because area A in the observing network has the largest sampling density, which is close to the resolution of the model grid, further refining the finest analysis grid of the MGA cannot provide additional useful information. In contrast, too few levels of the MGA cannot sufficiently reduce the analysis error. Results of other experiments with different

*a*values also demonstrate that 7 is a nearly optimal value of

*L*

_{MGA}(not shown). In the following experiments of the adaptive EnKF-MGA method,

*L*

_{MGA}is set to 7. Note that the observation sampling density may be high in some areas in reality. In such cases, according to the above selection criterion,

*L*

_{MGA}will be set to a large value, causing the resolution of the finest grid to be much higher than that of the model grids. When the analysis of the MGA on the fine grid is projected onto the coarser model grids, the small-scale (smaller than the model grid scale) information will be lost rapidly. In addition, overly large

*L*

_{MGA}will also significantly increase the computational cost. Thus,

*L*

_{MGA}should be determined by both observation sampling density and the resolution of the model grids in actual data assimilation. Generally speaking,

*L*

_{MGA}should be set to a value corresponding to the finest grid, the resolution of which is close to the smaller value of the highest model resolution and the largest observation sampling density.

### b. Performance of the MGA

We investigate the performance of the MGA from the following four aspects in this section.

First, we examine whether the MGA can further reduce the error produced by the EnKF without inflation at each analysis step. As in W14, we use the same three error statistics, which are as follows: 1) the RMSE of the analysis ensemble mean of _{EnKF-MGA}; see Eq. (34) in W14]; 2) the RMSE of the analysis ensemble mean of _{EnKF}, see Eq. (35) in W14]; and 3) the difference between RMSE_{EnKF-MGA} and RMSE_{EnKF} (i.e., RMSE_{EnKF-MGA} − RMSE_{EnKF}). A negative difference means that the MGA is valid.

Figure 5 shows the time series of the RMSE differences (10^{6} m^{2} s^{−1}) for *a* = 250 km (black) and *a* = 2500 km (blue), where the dashed line means no difference. The MGA can rapidly correct the model state during the spinup period, which is about 10 model days. After the spinup, the MGA can maintain the error of the model state at an equilibrium, which does not mean that the MGA is invalid. If the MGA is removed at a certain model step, the error of the model state will rebound to the level of the EnKF without inflation. To verify the correctness of the inference, we perform a test experiment of the adaptive EnKF-MGA, which artificially shuts off the MGA on the 40th model day for *a* = 250 km. Figure 6 plots the time series of the RMSE (10^{6} m^{2} s^{−1}) of the prior ensemble mean of

Second, to check whether the MGA can extract the multiscale information from the observational residual, we take the results of the first data assimilation cycle as an instance. We also adopt the following two quantities used in W14: 1) the “true” difference *a* = 250 km. Apparently, ^{6} m^{2} s^{−1}) of the observational error, causing the pattern in the observational residual to look like

Third, we check the time series of RMSE^{res} [Eq. (22)], because this parameter acts as the switch to invoke the MGA. Figure 9 shows the time series of RMSE^{res} (10^{6} m^{2} s^{−1}) for *a* = 250 km (Fig. 9a), *a* = 1500 km (Fig. 9b), and *a* = 2500 km (Fig. 9c), where the dashed line represents the critical value of 1.04. For *a* = 250 and 2500 km, the MGA is frequently activated. However, for *a* = 1500 km, the MGA is rarely called, justifying that the SNR of the observational residual is low. Thus, the assimilation quality of the adaptive EnKF-MGA method is nearly the same as that of the EnKF without inflation. Note that the activation of the MGA is completely adaptive rather than manual.

Last, because of the existence of the observational error in the observational residual, the analysis quality of the MGA strongly depends on the SNR of the observational residual. Here, taking the results of the adaptive method with *a* = 250 and 1500 km at 1800 UTC on model day 103 as examples, we analyze the SNR of the observational residuals. Figure 10 plots the spatial distributions of the observational residual (10^{6} m^{2} s^{−1}; Fig. 10a), the observational error (10^{6} m^{2} s^{−1}; Fig. 10b), the analysis error (10^{6} m^{2} s^{−1}; Fig. 10c) (in the observation space) of the first EnKF step in the adaptive EnKF-MGA, and the SNR for *a* = 250 km (Fig. 10d). Comparison between Figs. 10a–c illustrates that the observational residual is dominated by the analysis error rather than by the observational error, especially in area B and area C, where observations are sparsely distributed. The RMSE^{res} in this case is 1.11 × 10^{6} m^{2} s^{−1}, which is larger than the threshold *θ* of 1.04 × 10^{6} m^{2} s^{−1}. Thus, the MGA will be adaptively triggered. Figure 10d verifies that the SNRs are greater than 1.0 (the black contour in Fig. 10d) in most places, especially in area B and area C. Under this circumstance, the MGA will further refine the quality of the analysis ensemble mean of the first EnKF step (not shown). Figure 11 shows the same results as Fig. 10, but for *a* = 1500 km. Comparison between Figs. 11a–c reveals that the observational residual is dominated by the observational error, causing the SNRs in most places to be smaller than 1.0. The RMSE^{res} here is 0.98 × 10^{6} m^{2} s^{−1}, which automatically switches off the MGA. Although the above results are snapshots, the same conclusions can be drawn from other cases (not shown). Thus, the threshold is an effective criterion to measure whether the observational residual is signal dominant or not.

### c. Comparison with the standard EnKF

In this section, we compare the performance of the adaptive EnKF-MGA with that of the standard EnKF method from the following three aspects.

First, we investigate the dependence of the time-mean RMSE of the prior ensemble mean of the streamfunction on *a*. The black and red curves in Fig. 12 show the results of the standard EnKF and the adaptive EnKF-MGA, respectively. For comparison, we also plot the results (the dashed curve) of the EnKF without inflation. Note that the overly large assimilation errors of the EnKF without inflation for *a* larger than 2250 km are not shown in Fig. 12 and that the standard EnKF only works for *a* smaller than 4000 km. The vertical line at each *a* represents the ±*ζ* bounds of the time-mean RMSE, where *ζ* represents the standard deviation of the RMSE. If the error produced by method *A* falls out of the ±*ζ* bounds of the error produced by method *B*, we define that *A* is significantly different from *B*.

For small *a* (say 250 km), there is no overlap between the ±*ζ* bound of the time-mean RMSE for the standard EnKF and that for the adaptive EnKF-MGA, demonstrating that the assimilation quality of the adaptive EnKF-MGA is significantly better than that of the standard EnKF.

For moderate *a* (between 500 and 1750 km), the adaptive EnKF-MGA reduces to the EnKF without inflation because the MGA is rarely called. The ±*ζ* bound of the time-mean RMSE for the adaptive method is nearly overlapped by that for the standard EnKF, justifying that the difference between the adaptive EnKF-MGA and the EnKF without inflation is not significant.

For relatively large *a* (between 2000 and 3500 km), error bars show that the adaptive EnKF-MGA is much better than the EnKF without inflation but worse than the standard EnKF. According to Fig. 9c, the MGA is frequently called by the adaptive method for these *a* values: that is, the analysis ensemble mean of the first EnKF step in the adaptive method is often adjusted by the MGA. Under this circumstance, the analysis ensemble perturbation of the first EnKF step in the adaptive method should also be consistently adjusted. Thus, there are two possible reasons that cause the inferior performance of the adaptive EnKF-MGA compared to that of the standard EnKF: one is the lack of variance inflation, and the other is the preservation of the ensemble perturbations. We will discuss which one is the dominant factor in section 5.

For overly large *a* (say 3750 km), the RMSE produced by the adaptive method falls out of the ±*ζ* bound produced by the standard EnKF, proving that there is significant difference between the adaptive method and the standard EnKF, and the former method is better than the latter method. Thus, the adaptive EnKF-MGA can generally weaken the dependence of the assimilation quality on *a.* Besides, the adaptive EnKF-MGA works for a broader range of *a* than the standard EnKF. For example, the adaptive EnKF-MGA produces a small assimilation error for *a* = 4000 km, while the standard EnKF does not work in this case. We even performed an extreme experiment of the adaptive EnKF-MGA, which applied the MGA to an unlocalized EnKF: that is, all observations were allowed to impact global model grids during the first EnKF step in the adaptive EnKF-MGA. Results showed that the EnKF system can run without localization when it is used with the MGA, though the error is very large (i.e., the time-mean RMSE = 2.4 × 10^{6} m^{2} s^{−1}).

Second, taking three typical values of *a* (i.e., 250, 1500, and 3750 km) as examples, we conduct error analysis for the standard EnKF and the adaptive EnKF-MGA data assimilation schemes. Figure 13 shows the RMSE time series of the prior ensemble mean of the streamfunction for *a* = 250 km (Fig. 13a), *a* = 1500 km (Fig. 13b), and *a* = 3750 km (Fig. 13c), where the black, red, and green curves represent the results of the standard EnKF, the adaptive EnKF-MGA and the EnKF-MSA methods, respectively. For overly small or overly large values of *a*, the error produced by the adaptive EnKF-MGA is smaller than that produced by the standard EnKF, and the spinup period of the data assimilation is shortened. Although the adaptive EnKF-MGA method is worse than the standard EnKF when *a* is moderate, the difference looks small. Figure 14 gives the spinup periods of the standard EnKF (black) and the adaptive EnKF-MGA (red) methods for different *a* values. If we roughly average the spinup periods for small (250 and 500 km) and large (larger than 2000 km) values of *a*, the adaptive EnKF-MGA can shorten the spinup period of the standard EnKF by 53%.

*a*= 250 km (Figs. 15a,b),

*a*= 1500 km (Figs. 15d,e), and

*a*= 3750 km (Figs. 15g,h) as examples. Here, RMSE is defined aswhere (

*i*,

*j*) indexes the model grid. For

*a*= 250 km and

*a*= 3750 km, the adaptive EnKF-MGA method can greatly reduce the errors in area B and area C compared to the standard EnKF. This indicates that the adaptive EnKF-MGA is superior to the standard EnKF in sparsely observed areas for extreme

*a*values. For

*a*= 1500 km, the adaptive EnKF-MGA enhances the global errors.

In summary, the adaptive EnKF-MGA method has four advantages over the standard EnKF in terms of assimilation quality: smaller assimilation errors and shorter spinup periods for extreme *a* values, weaker dependence on *a*, and a broader application range of *a*. These merits are shared by the EnKF-MSA method. It should be mentioned that the EnKF-MSA was compared to the EnKF without inflation in W14. Unlike for the EnKF-MSA method, the adaptive EnKF-MGA can reduce to the EnKF without inflation, which is better than the EnKF-MSA (the green curve in Fig. 12) for moderate *a* values.

### d. Comparison with the EnKF-MSA

In this section, we compare the performance of the adaptive EnKF-MGA to that of the EnKF-MSA on the assimilation quality based on Fig. 12, in which the green curve shows the results of the EnKF-MSA. For moderate *a* values, the RMSE produced by the adaptive method (the red curve in Fig. 12) is smaller than the −*ζ* bound produced by the EnKF-MSA, demonstrating that the adaptive EnKF-MGA is better than the EnKF-MSA that triggers the MSA at each analysis. Since the SNR of the observational residual is low (see Fig. 11d for the *a* = 1500-km case), it is hard for the MSA to effectively extract the analysis error (in the observation space) in the observational residual. Thus, the EnKF-MSA worsens the assimilation quality of the EnKF without inflation. This can also be demonstrated by the results of the EnKF-MSA in Fig. 13b (the green curve) and Fig. 15f for *a* = 1500 km.

For small *a* values (250 and 500 km), the *ζ* bound produced by the EnKF-MSA is smaller than the −*ζ* bound produced by the adaptive EnKF-MGA, demonstrating the EnKF-MSA is better than the adaptive EnKF-MGA. Taking *a* = 250 km as an instance, the green curve in Fig. 13a shows the RMSE time series, and Fig. 15c shows the spatial distribution of the RMSE of the prior ensemble mean of the streamfunction for the EnKF-MSA. Comparison between these results and those (see the red curve in Fig. 13a and Fig. 15b) produced by the adaptive EnKF-MGA also supports the above conclusion.

For large *a* values (over 2000 km), because the ±*ζ* bound produced by the adaptive EnKF-MGA (the EnKF-MSA) covers the RMSE produced by the EnKF-MSA (the adaptive EnKF-MGA), the superiority of the EnKF-MSA over the adaptive EnKF-MGA is not significant. Taking *a* = 3750 km for example, the green curve in Fig. 13c shows the RMSE time series, and Fig. 15i shows the spatial distribution of the RMSE of the prior ensemble mean of the streamfunction for the EnKF-MSA. Comparison between these results and those (see the red curve in Fig. 13c and Fig. 15h) produced by the adaptive EnKF-MGA also justifies that there is no significant difference between the RMSE produced by the EnKF-MSA and that produced by the adaptive EnKF-MGA.

Thus, the adaptive EnKF-MGA has an incremental improvement over the EnKF-MSA for moderate *a* values. In contrast, the EnKF-MSA is superior to the adaptive EnKF-MGA for small *a* values. However, it should be noted that this superiority is based on specific selections of *a*^{(l)} and *L*_{MSA} in the MSA, which are tuned, rather than adaptively determined.

### e. Computational cost

In this section, we simply use the wall-clock time (in minutes) of an experiment to measure the computational cost. Figure 16 shows the computational costs of the standard EnKF (black), the adaptive EnKF-MGA (red), the MGA (red dashed), the EnKF-MSA (green), and the MSA (green dashed) algorithms. All experiments are sequentially conducted on the same computational platform. Here, the computational cost of the MGA (MSA) is computed as the difference between the counterparts of the adaptive EnKF-MGA (the EnKF-MSA) and the standard EnKF.

The adaptive EnKF-MGA requires only slightly more wall-clock time compared to the standard EnKF, especially for moderate impact radii, which only trigger the MGA at a handful of times. If we calculate the average value of the computational costs for all impact radii, the computational cost of the MGA (2.1 min) is about 6% of that of the standard EnKF (37.0 min). Thus, relative to the computational expense of the standard EnKF, the cost of the MGA is almost negligible. In contrast, the EnKF-MSA increases the computational cost of the standard EnKF by 76% (from 37.0 to 65.2 min). The MGA can reduce the computational cost of the MSA by 93% (from 28.2 to 2.1 min). Results here are consistent with the analysis of the computational cost of the MSA in section 2c and of the MGA in section 2f. Thus, the MGA is much more computationally effective than the MSA.

## 5. Summary and discussion

The fixed covariance localization is often used to suppress the spurious long-distance correlation evaluated by a finite ensemble in EnKF. Although fixed covariance localization can greatly improve analysis quality, it has significant limitations, including insufficient longwave information with a small cutoff distance, contaminated analysis states with a large cutoff distance, and nearly unavailable optimal cutoff distance. W14 presented an EnKF and MSA hybrid method to make up for the above demerits. However, the compensatory method has some disadvantages. For example, the computational cost of the MSA is huge for high-resolution and/or high-dimension numerical models and dense observing systems; and the assimilation quality of the compensatory method is worse than the EnKF without inflation for moderate impact radii. To avoid these issues, we presented an adaptive EnKF and MGA hybrid method in this paper. At each analysis step, after the EnKF without inflation is completed, a threshold, which depends on the observing system, is used to decide whether to activate the MGA or not. Similar to the MSA algorithm, the MGA technique is used to extract multiscale information from observational residuals (the differences between observations and interpolated analysis ensemble means produced by the first EnKF step) by gradually refining the analysis grid. The analysis of the MGA is added to the analysis ensemble mean of the first EnKF step to form the final analysis ensemble mean. Then the analysis ensemble perturbations produced by the first EnKF step are added to the final analysis ensemble mean to form the final analysis ensemble members. Within a biased twin experiment framework consisting of a global barotropic spectral model and an idealized observing system, the performance of the proposed method was examined.

There are four advantages of the adaptive EnKF-MGA over the standard EnKF. First, the adaptive EnKF-MGA can produce smaller assimilation errors for overly large or overly small *a* values. Second, the adaptive EnKF-MGA can shorten the spinup period of the standard EnKF by 53% for both small and large *a* values. Third, the adaptive EnKF-MGA has less dependence on *a* than the standard EnKF. Last, the adaptive EnKF-MGA works for a broader range of *a* than the standard EnKF. For the twin experiments in this study, the EnKF system can run without localization when it is used with the MGA, though the error is very large. At present, although several adaptive localization schemes (e.g., Anderson 2007) have been presented, the fixed localization scheme is still widely used in some operational centers that establish an ensemble data assimilation system. Because of the flow-dependent nature of model states, it is difficult for them to set a fixed impact radius working for all times and spaces. Under this circumstance, the MGA method appears to be a potential cost-effective means of adjusting the output of an EnKF.

There are two demerits of the adaptive EnKF-MGA over the standard EnKF. On the one hand, for moderate *a* values, because of the absence of the variance inflation and the low SNR of the observational residual, the adaptive EnKF-MGA automatically reduces to the EnKF without inflation, which is worse than the standard EnKF. This demerit is expected to be corrected through introducing variance inflation. On the other hand, for relatively large *a* values, the adaptive EnKF-MGA is much better than the EnKF without inflation but worse than the standard EnKF. As we analyzed in section 4c, there are two possible reasons causing the above disadvantages of the adaptive EnKF-MGA over the standard EnKF. One is the absence of variance inflation in the first EnKF step of the adaptive EnKF-MGA. The other is not dealing with the ensemble perturbation. To figure out which is the dominant causation, we also apply the multiplicative variance inflation to the adaptive EnKF-MGA to find the optimal inflation factor (the red curve in Fig. 3) for each *a* value. Note that the variance inflation is only applied to the analysis step whose prior analysis step does not trigger the MGA. The inflation factors for the adaptive EnKF-MGA experiments are much smaller than those for the standard EnKF (the black curve in Fig. 3). The reason is that, once the analysis ensemble mean of the first EnKF step is improved by the MGA, the ensemble spread of the analysis ensemble should be reduced. Thus, the inflation needed for the standard EnKF is reduced. Without consistently adjusting the ensemble perturbations, the optimal inflation factor is also reduced. The red dashed curve in Fig. 17 shows the dependence of the assimilation RMSE on *a* for the adaptive EnKF-MGA with inflation. For comparison, results of the adaptive EnKF-MGA (the red curve) and the standard EnKF (the black curve) are also presented. The red and black curves here are the same as those in Fig. 12. Comparison between the red curve and the red dashed curve indicates that the inflation can improve the performance of the adaptive EnKF-MGA to approximate the standard EnKF for moderate *a* values, which rarely call the MGA. The RMSE produced by the adaptive EnKF-MGA with inflation is smaller than the −*ζ* bound produced by the adaptive EnKF-MGA, demonstrating that inflation can significantly improve the assimilation quality of the adaptive EnKF-MGA for moderate *a* values. For relatively large *a* values, there is almost no difference between the RMSE and its bound produced by the adaptive EnKF-MGA and those produced by the adaptive EnKF-MGA with inflation, justifying that the inflation is invalid when the MGA is frequently called. This demonstrates that the lack of consistent adjustment of the ensemble perturbations is the main cause of the inferior performance of the adaptive EnKF-MGA. This issue should be addressed in future work. In addition, in the current version of the adaptive EnKF-MGA, the computational cost of the MGA is negligible relative to that of the standard EnKF.

Comparison between the adaptive EnKF-MGA and the EnKF-MSA shows that the MGA can realize the basic functions of the MSA and that the adaptive EnKF-MGA has three superiorities over the EnKF-MSA. First and foremost, the MGA can reduce the computational cost of the MSA by 93%. Second, on the assimilation quality, the adaptive EnKF-MGA has an incremental improvement over the EnKF-MSA for moderate *a* values. That is, the adaptive EnKF-MGA automatically reduces to the EnKF without inflation, which is better than the EnKF-MSA for moderate *a* values. Last, on the selection of the assimilation parameter, the determination of the analysis grids of the MGA is straightforward and can be adaptive, while the *a*^{(l)} and *L*_{MSA} in the MSA have to be manually determined. When the values of *a*^{(l)} and *L*_{MSA} are properly selected, the assimilation errors produced by the EnKF-MSA may be smaller than those produced by the adaptive EnKF-MGA, which is like what the EnKF-MSA did in this study.

Some other challenges exist before the adaptive EnKF-MGA can be applied for actual weather and climate analysis. First, the MGA mainly serves as a mathematical tool. As we know, the coarser the grid, the larger the representative error of the observations. This effect should be included in ^{(l)}. Also, the shape of the smoothing term could be further improved to have more physical meanings. Second, the model used in this study is fairly simple compared to realistic models that have highly nonlinear processes like condensation. This raises issues that may affect the performance of the adaptive EnKF-MGA. Simply shifting the whole ensemble according to the adaptive EnKF-MGA may introduce undesirable effects in more complex models. Can the proposed method beat the standard EnKF with a reasonable but not optimal length scale in GCMs? These problems would not be visible in the system studied here. Thus, the MGA should be applied to assimilate actual instrumental measurements into a sophisticated atmospheric, oceanic, or atmosphere–ocean coupled GCM to further verify the validity of the adaptive EnKF-MGA method. Last, for a large ensemble size, the optimal impact radius would increase, which reduces the loss of longwave information. Thus, the performance of the proposed method should also be investigated using large ensemble sizes for extreme impact radii.

## Acknowledgments

This research is cosponsored by grants from the National Basic Research Program (2013CB430304), the National Natural Science Foundation (41306006, 41376015, 41376013, 41176003, and 41206178), the National High-Tech R&D Program (2013AA09A505), and Global Change and Air–Sea Interaction (GASI-01-01-12) of China.

## REFERENCES

Anderson, J. L., 2001: An ensemble adjustment Kalman filter for data assimilation.

,*Mon. Wea. Rev.***129**, 2884–2903, doi:10.1175/1520-0493(2001)129<2884:AEAKFF>2.0.CO;2.Anderson, J. L., 2007: Exploring the need for localization in ensemble data assimilation using a hierarchical ensemble filter.

,*Physica D***230**, 99–111, doi:10.1016/j.physd.2006.02.011.Anderson, J. L., 2008: Spatially and temporally varying adaptive covariance inflation for ensemble filters.

,*Tellus***61A**, 72–83, doi:10.1111/j.1600-0870.2008.00361.x.Anderson, J. L., , and S. L. Anderson, 1999: A Monte Carlo implementation of the nonlinear filtering problem to produce ensemble assimilations and forecasts.

,*Mon. Wea. Rev.***127**, 2741–2758, doi:10.1175/1520-0493(1999)127<2741:AMCIOT>2.0.CO;2.Anderson, J. L., , and L. L. Lei, 2013: Empirical localization of observation impact in ensemble Kalman filters.

,*Mon. Wea. Rev.***141**, 4140–4153, doi:10.1175/MWR-D-12-00330.1.Asselin, R., 1972: Frequency filter for time integrations.

,*Mon. Wea. Rev.***100**, 487–490, doi:10.1175/1520-0493(1972)100<0487:FFFTI>2.3.CO;2.Berre, L., , and G. Desroziers, 2010: Filtering of background error variance and correlations by local spatial averaging: A review.

,*Mon. Wea. Rev.***138**, 3693–3720, doi:10.1175/2010MWR3111.1.Bishop, C. H., , and D. Hodyss, 2007: Flow adaptive moderation of spurious ensemble correlations and its use in ensemble-based data assimilation.

,*Quart. J. Roy. Meteor. Soc.***133**, 2029–2044, doi:10.1002/qj.169.Bishop, C. H., , and D. Hodyss, 2009a: Ensemble covariances adaptively localized with ECO-RAP. Part 1: Tests on simple error models.

,*Tellus***61A**, 84–96, doi:10.1111/j.1600-0870.2008.00371.x.Bishop, C. H., , and D. Hodyss, 2009b: Ensemble covariances adaptively localized with ECO-RAP. Part 2: A strategy for the atmosphere.

,*Tellus***61A**, 97–111, doi:10.1111/j.1600-0870.2008.00372.x.Bishop, C. H., , and D. Hodyss, 2011: Adaptive ensemble covariance localization in ensemble 4D-VAR state estimation.

,*Mon. Wea. Rev.***139**, 1241–1255, doi:10.1175/2010MWR3403.1.Bishop, C. H., , D. Hodyss, , P. Steinle, , H. Sims, , A. M. Clayton, , A. C. Lorenc, , D. M. Barker, , and M. Buehner, 2011: Efficient ensemble covariance localization in variational data assimilation.

,*Mon. Wea. Rev.***139**, 573–580, doi:10.1175/2010MWR3405.1.Briggs, W. L., , V. E. Henson, , and S. F. McCormick, 2000:

2nd ed. Society for Industrial and Applied Mathematics, 193 pp.*A Multigrid Tutorial.*Evensen, G., 2007:

*Data Assimilation: The Ensemble Kalman Filter*. Springer, 187 pp.Gaspari, G., , and S. E. Cohn, 1999: Construction of correlation functions in two and three dimensions.

,*Quart. J. Roy. Meteor. Soc.***125**, 723–757, doi:10.1002/qj.49712555417.Greybush, S. J., , E. Kalnay, , T. Miyoshi, , K. Ide, , and B. R. Hunt, 2011: Balance and ensemble Kalman filter localization techniques.

,*Mon. Wea. Rev.***139**, 511–522, doi:10.1175/2010MWR3328.1.Haltiner, G. J., , and R. T. Williams, 1980:

. 2nd ed. Wiley, 477 pp.*Numerical Prediction and Dynamic Meteorology*Hamill, T. M., , J. S. Whitaker, , and C. Snyder, 2001: Distance-dependent filtering of background error covariance estimates in an ensemble Kalman filter.

,*Mon. Wea. Rev.***129**, 2776–2790, doi:10.1175/1520-0493(2001)129<2776:DDFOBE>2.0.CO;2.Houtekamer, P. L., , and H. L. Mitchell, 1998: Data assimilation using an ensemble Kalman filter technique.

,*Mon. Wea. Rev.***126**, 796–811, doi:10.1175/1520-0493(1998)126<0796:DAUAEK>2.0.CO;2.Houtekamer, P. L., , H. L. Mitchell, , and X. Deng, 2009: Model error representation in an operational ensemble Kalman filter.

,*Mon. Wea. Rev.***137**, 2126–2143, doi:10.1175/2008MWR2737.1.Hunt, B. R., , E. J. Kostelich, , and I. Szunyogh, 2007: Efficient data assimilation for spatiotemporal chaos: A local ensemble transform Kalman filter.

,*Physica D***230**, 112–126, doi:10.1016/j.physd.2006.11.008.Jun, M., , I. Szunyogh, , M. G. Genton, , F. Zhang, , and C. H. Bishop, 2011: A statistical investigation of the sensitivity of ensemble-based Kalman filters to covariance filtering.

,*Mon. Wea. Rev.***139**, 3036–3051, doi:10.1175/2011MWR3577.1.Lei, L. L., , and J. L. Anderson, 2014: Empirical localization of observations for serial ensemble Kalman filter data assimilation in an atmospheric general circulation model.

,*Mon. Wea. Rev.***142**, 1835–1851, doi:10.1175/MWR-D-13-00288.1.Li, W., , Y. Xie, , Z. He, , G. Han, , K. Liu, , J. Ma, , and D. Li, 2008: Application of the multigrid data assimilation scheme to the China Seas’ temperature forecast.

,*J. Atmos. Oceanic Technol.***25**, 2106–2116, doi:10.1175/2008JTECHO510.1.Li, W., , Y. Xie, , S.-M. Deng, , and Q. Wang, 2010: Application of the multigrid method to the two-dimensional Doppler radar radial velocity data assimilation.

,*J. Atmos. Oceanic Technol.***27**, 319–332, doi:10.1175/2009JTECHA1271.1.Liu, D. C., , and J. Nocedal, 1989: On the limited memory BFGS method for large scale optimization.

,*Math. Program.***45**, 503–528, doi:10.1007/BF01589116.Mitchell, H. L., , P. L. Houtekamer, , and G. Pellerin, 2002: Ensemble size, balance, and model-error representation in an ensemble Kalman filter.

,*Mon. Wea. Rev.***130**, 2791–2808, doi:10.1175/1520-0493(2002)130<2791:ESBAME>2.0.CO;2.Miyoshi, T., 2011: The Gaussian approach to adaptive covariance inflation and its implementation with the local ensemble transform Kalman filter.

,*Mon. Wea. Rev.***139**, 1519–1535, doi:10.1175/2010MWR3570.1.Miyoshi, T., , and K. Kondo, 2013: A multi-scale localization approach to an ensemble Kalman filter.

,*SOLA***9**, 170–173, doi:10.2151/sola.2013-038.Ott, E., and et al. , 2004: A local ensemble Kalman filter for atmospheric data assimilation.

,*Tellus***56A**, 415–428, doi:10.1111/j.1600-0870.2004.00076.x.Robert, A., 1969: The integration of a spectral model of the atmosphere by the implicit method.

*Proc. WMO/IUGG Symp. on Numerical Weather Prediction*, Tokyo, Japan, Japan Meteorological Society, 19–24.Szunyogh, I., , E. J. Kostelich, , G. Gyarmati, , E. Kalnay, , B. R. Hunt, , E. Ott, , E. Satterfield, , and J. A. Yorke, 2008: A local ensemble transform Kalman filter data assimilation system for the NCEP global model.

,*Tellus***60A**, 113–130, doi:10.1111/j.1600-0870.2007.00274.x.Wu, X., , W. Li, , G. Han, , S. Zhang, , and X. Wang, 2014: A compensatory approach of the fixed localization in EnKF.

,*Mon. Wea. Rev.***142**, 3713–3733, doi:10.1175/MWR-D-13-00369.1.Xie, Y., , S. Koch, , J. McGinley, , S. Albers, , P. E. Bieringer, , M. Wolfson, , and M. Chan, 2011: A space–time multi-scale analysis system: A sequential variational analysis approach.

,*Mon. Wea. Rev.***139**, 1224–1240, doi:10.1175/2010MWR3338.1.Yang, S.-C., , E. Kalnay, , and B. R. Hunt, 2012: Handling nonlinearity in an ensemble Kalman filter: Experiments with the three-variable Lorenz model.

,*Mon. Wea. Rev.***140**, 2628–2646, doi:10.1175/MWR-D-11-00313.1.Zhang, F., , C. Snyder, , and J. Sun, 2004: Impacts of initial estimate and observation availability on convective-scale data assimilation with ensemble Kalman filter.

,*Mon. Wea. Rev.***132**, 1238–1253, doi:10.1175/1520-0493(2004)132<1238:IOIEAO>2.0.CO;2.Zhang, S., , J. L. Anderson, , A. Rosati, , M. J. Harrison, , S. P. Khare, , and A. Wittenberg, 2004: Multiple time level adjustment for data assimilation.

,*Tellus***56A**, 2–15, doi:10.1111/j.1600-0870.2004.00040.x.