Impact of Observational Time Window on Coupled Data Assimilation : Simulation with a Simple Climate Model

Climate signals are the results of interactions of multiple time scale media such as the atmosphere and ocean in the coupled earth system. Coupled data assimilation (CDA) pursues balanced and coherent climate analysis and prediction 15 initialization by incorporating observations from multiple media into a coupled model. In practice, an observational time window (OTW) is usually used to collect measured data for an assimilation cycle to increase observational samples that are sequentially assimilated with their original error scales. Given different time scales of characteristic variability in different media, what are the optimal OTWs for the coupled media so that climate signals can be most accurately recovered by CDA? With a simple coupled model that simulates typical scale interactions in the climate system and “twin” CDA experiments, we 20 address this issue here. Results show that in each coupled medium, an optimal OTW can provide maximal observational information that best fits characteristic variability of the medium during the data blending process. Maintaining correct scale interactions, the resulted CDA improves the analysis of climate signals greatly. This simple model results provide a guideline when the real observations are assimilated into a coupled general circulation model for improving climate analysis and prediction initialization by accurately recovering important characteristic variability such as sub-diurnal in the atmosphere 25 and diurnal in the ocean.


Introduction
Currently, the interactions between the earth climate system's major components, such as the atmosphere, ocean, land, and sea ice, have been reasonably simulated by coupled climate models, which can also give the evaluation of climate changes (Randall et al., 2007).However, because of the uncertainties and errors in models (e.g., parameterization is only an approximation to sub-grid processes and the dynamical core is imperfect), models always tend to produce different climate features and variability from the real world (e.g., Delworth et al., 2006;Collins et al., 2006;Zhang et al., 2014).Due to the significant importance of preserving the balance and coherence of different model components (or media) during the coupled model initialization, data assimilation for state estimation and prediction initialization should be performed within a coupled climate model framework (e.g., Chen et al., 1995;Zhang et al., 2007;Chen, 2010;Han et al., 2013).The characteristic variability timescales of different media within the coupled frameworks are usu-Y.Zhao et al.: Impact of an observational time window ally different.When the observed data included in one or more components of the coupled system framework are assimilated, the observational information will be able to be transferred among different media through the coupled dynamics so that all media gain consistent and coherent adjustments.Such an assimilation procedure is called coupled data assimilation (CDA), which can sustain the nature of multiple timescale interactions during climate estimation and prediction initialization (e.g., Zhang et al., 2007;Sugiura et al., 2008;Singleton, 2011), thus producing better climate analysis and prediction initialization and therefore improving the coupled models' predictability (Yang et al., 2013).Zhang et al. (2007) developed the first CDA system in a fully coupled general circulation model, version 2 of the Geophysical Fluid Dynamics Laboratory Coupled Model (GFDL CM2).The National Centres for Environmental Prediction (NCEP) also started using coupled models to generate first-guess forecasts for their Climate Forecast System Reanalysis (CFSR, Saha et al., 2010).Despite the enormous benefits and demand for CDA, it remains both theoretically and technically challenging to implement strong CDA in fully coupled models, including the estimation of the coupled model error covariance matrix and the huge computational costs (e.g., Han et al., 2013;Lu et al., 2015;Liu et al., 2016).
During the coupled data assimilation process, an observational time window (OTW) is usually used to collect measured data in each medium for an assimilation cycle (e.g., Pires et al., 1996;Hunt et al., 2004;Houtekamer and Mitchell, 2005;Laroche et al., 2007) to increase observational samples.As in Hunt et al. (2004), we expand the EnKF to include a time window in which the observations are treated as the exact assimilation times, even though their times are different in the window.That is, we just assume that all the collected data sample the "truth" variation at the assimilation time and will be sequentially assimilated with their original error scales.Thus the OTW is applied in a three-dimensional data assimilation fashion rather than a four-dimensional one.Apparently, while a large OTW provides more observational samples at the assimilation time, the assimilation process blends more data from different times and may distort the variability being retrieved.Given the fact that climate signals are the results of interactions of multiple timescale media, correct variability retrieved for each medium so that correct scale interaction is maintained in CDA is particularly important for climate analysis and prediction initialization.In this study we attempt to answer the following two questions.(1) What is the impact of varying OTWs for each coupled component within the coupled model framework on the quality of CDA? (2) Based on this impact, does an optimal OTW exist so that assimilation fitting has maximum observational information but minimum variability distortion?
With a simple conceptual coupled climate model and a sequential implementation of the ensemble Kalman filter, this study first analyses the characteristic variability timescale of each coupled medium and identifies the corresponding optimal OTW.Then the impact of an optimal OTW on the quality of CDA and its linkage with the corresponding timescale of characteristic variability are investigated.The simple coupled model consists of three typical components, including the synoptic atmosphere (Lorenz, 1963) and the seasonalinterannual slab upper ocean (Zhang et al., 2012) coupling with the decadal deep ocean (Zhang, 2011a, b).Although the simple conceptual coupled model does not share the similar complex physics with a coupled general circulation model (CGCM), it does reasonably simulate the typical interactions between multiple timescale components in the coupled climate system (see Zhang et al., 2013).The simple coupled model helps us understand the essence of the problem by revealing the relationship between the optimal OTWs and corresponding timescales of characteristic variability as well as their impact on CDA.The low-cost nature of the simple model also provides convenience for a large number of CDA experiments with different OTWs in optimal OTW detection.The ensemble Kalman filter (e.g., Evensen, 1994Evensen, , 2007;;Whitaker and Hamill, 2002;Anderson, 2001Anderson, , 2003) ) used in this study is the ensemble adjustment Kalman filter (EAKF, e.g., Anderson, 2001Anderson, , 2003;;Zhang and Anderson, 2003).Using the EAKF with the simple coupled model, we first establish a twin experiment framework.Within such a framework, the degree to which the state estimation based on a certain OTW recovers the truth is an assessment of the influence of the OTW on the quality of CDA.In such a way, the optimal OTW of each medium is detected and the impact of optimal OTWs on CDA is evaluated.We also discuss the influence of model bias on an optimal OTW through biased twin experiment setting.
This paper is organized as follows.Section 2 briefly describes the simple conceptual coupled model, the ensemble adjustment Kalman filter, as well as the twin experiment framework including perfect and biased settings.Also with a simplest case, we first show the influence of OTWs on assimilation quality and its linkage with the timescale of characteristic variability in this section.Then Sect. 3 presents results on detection of the optimal OTWs for different media and the impact of optimal OTWs on CDA.The influence of realistic assimilation scenarios on optimal OTWs is discussed in Sect. 4. Finally, a summary and discussions are given in Sect. 5.

The model
Due to the complicated physical processes and huge computational cost involved, it is inconvenient to use a CGCM to investigate the impact of the different OTWs on the analysis of climate signals so as to detect each coupled medium's optimal OTW.Instead, here we employ a simple coupled "climate" model developed by Zhang (2011a).This sim-ple model is based on Lorenz's three-variable chaotic model (Lorenz, 1963) that couples with a slab upper ocean (Zhang et al., 2012) and a simple pycnocline predictive model (Gnanadesikan, 1999).Although very simple with low computational cost, in terms of multi-scale interaction inducing low-frequency climate signals, this model shares a fundamental character with a CGCM, and it is very suitable for addressing the problem that is concerned here.And for the readers' convenience, here we simply review some key aspects of this conceptual coupled model.With all quantities being given in non-dimensional units, the governing equations are where X 1 , X 2 , and X 3 represent the atmospheric model states, while ω and η denote those for the upper and deep ocean, respectively.A dot above the variable denotes the time tendency.The atmosphere model states are the highfrequency variables, while the slab oceanic variable ω is of a lower frequency.To sustain the chaotic nature of the atmosphere in reality, the standard values of the parameters included in the atmospheric component (σ , k, and b) are set as 9.95, 28, and 8/3, respectively.In the equation of ω, the parameters O d and O m denote the damping coefficient and heat capacity of the upper slab ocean, respectively.Due to the lower frequency of ω than that of the model states in the atmospheric components, the timescale of the upper slab ocean variable must be much slower than that of the atmospheric model states.Thus the damping rate parameter (O d ) should be much smaller than the heat capacity, namely, O d O m .Here following Lorenz's idea (Lorenz, 1963), the atmospheric timescale is defined as the typical time by which the atmosphere goes through an attractive lob as 1 non-dimensional time unit (TU) ∼ O(1).We set the parameters (O m , O d ) as (10,1), which show that the slab oceanic variable's timescale is ∼ O(10), i.e., 10 times that of the atmospheric model states.While the S m + S s cos 2π t/S pd represents the external forcing, the parameter S pd denoting the model seasonal cycle is set as 10 to make sure that the period of the external forcing is comparable with the upper slab ocean variables' timescale.In this simple coupled model, the seasonal cycle is set as 10 TUs and thus a model year (decade) equals 10 (100) TUs.The parameters S s and S m , denoting the magnitudes of the external forcing's seasonal cycle and annual mean, are insensitive to the coupled model and set as (1,10).The coefficients C 1 and C 2 in the equations of X 2 and ω are used to implement the coupling between the fast atmosphere model states and the upper slab oceanic variable and are set as (0.1,1), where C 1 denotes the upper slab oceanic forcing on the atmosphere while C 2 denotes the atmosphere forcing on the ocean.In addition, C 3 and C 4 represent the deep oceanic forcing and the nonlinear interaction between the upper and deep ocean.
In order to make sure that the atmospheric forcing plays a dominant role in the upper slab ocean, the magnitudes of C 3 and C 4 should be lower than that of C 2 and both set as 0.01.As in Zhang (2011a), the deep ocean model state variable η, denoting the anomaly of pycnocline depth in the deep ocean, is derived from the two-term balance model of the zonal-time mean pycnocline (Gnanadesikan, 1999).
Within the equation of η, the parameter is kept constant and the ratio of and O d denotes the deep ocean variable's timescale.The timescale of the deep ocean variable is longer than that of the slab ocean, defined by the relative magnitude of to O d ( is set as 100).Similar to the equation of ω, the coefficients C 5 and C 6 denote the linear slab oceanic forcing and the nonlinear interaction between the upper and deep ocean.Also, to guarantee that the linear interaction is dominant and the nonlinear interaction is weaker than that in the deep ocean model, C 5 and C 6 are set as (1, 0.001).In summary, in this study the standard values of the parameters included in this simple coupled model ) are set as (9.95, 28, 8/3, 0.1, 1, 1, 10, 10, 1, 10, 100, 0.01, 0.01, 1, 0.001; e.g., Zhang, 2011a, b;Zhang et al., 2012;Han et al., 2013Han et al., , 2014)).
Following the study of Han et al. (2014), the fourth-order Runge-Kutta time-differencing scheme is used in this paper to resolve this simple coupled model, and the time step equals 0.01 TU (1 TU = 100 time steps).
Zhang (2011b) illustrated that, given the model parameters described above, the constructed simple coupled model can effectively simulate a fundamental feature of the realworld climate system in which different timescales interact with each other to develop climate signals.That is, the synoptic to decadal timescale signals are produced by the interactions between the transient atmosphere attractor, the slow slab ocean, and the even slower deep ocean (see Zhang, 2011a;Han et al., 2014).Again, although the simple coupled model does not have complex physics and cannot consider the issue of impact of localization and imbalance as in a CGCM, it can help us investigate the fundamental issue we want to address here more directly and clearly.

Ensemble coupled data assimilation
Following Zhang (2011a), during the state estimation, the error statistics evaluated from ensemble model integrations, such as the error covariance between model states, will be used in an ensemble filter to extract observational information to adjust the model states (e.g., Evensen, 1994Evensen, , 2007;;Anderson, 2001;Hamill et al., 2000;Zhang, 2011a, b;Zhang et al., 2012;Han et al., 2014).In this study, a derivative of www.nonlin-processes-geophys.net/24/681/2017/ Nonlin.Processes Geophys., 24, 681-694, 2017 the Kalman filter (Kalman, 1960;Kalman and Bucy, 1961) called the ensemble adjustment Kalman filter (EAKF, Anderson, 2001Anderson, , 2003;;Zhang and Anderson, 2003;Zhang et al., 2007), which is a sequential implementation of the ensemble Kalman filter under an "adjustment" idea, is used to implement the CDA scheme.The assumption of independence of observational error allows the EAKF to sequentially assimilate observations into corresponding model states (Zhang and Anderson, 2003;Zhang et al., 2007).While the sequential implementation provides much computation convenience for data assimilation, the EAKF maintains as much of the nonlinearity of background flows as possible (Anderson, 2001(Anderson, , 2003;;Zhang and Anderson, 2003).
Based on the two-step implementation of the EAKF scheme (Anderson, 2001(Anderson, , 2003)), the observational increment at an observation location is first computed.The observation is denoted as Y at time t (simply Y instead of Y t ) which has the observation value Y o and standard deviation σ o y (assumed to be Gaussian).Firstly, the reshaping of the model ensemble at the observation location, Y , is formulated as where i represents the ensemble index and k denotes the observation index.σ o k,k and σ p k,k are the standard deviation of observation error and its prior estimated ensemble standard deviation, respectively, while r k is the corresponding ratio.If r k > 1 , the ensemble spread is largely reduced by the observation; otherwise, the ensemble remains close to the prior.The shift of the ensemble mean induced by the observation is computed by We can see that if the prior estimated ensemble standard deviation is greater than that of the observation error, the ensemble mean shifts toward the observation value; otherwise, the ensemble mean remains close to the prior model ensemble mean Y p .Then the observational increment induced by the observation value Y o for the ith ensemble member at the kth observation location is computed as Once we get the observational increment at the observation location, then a least square fit is used to distribute the increment over the relevant grid points impacted by the observation using the covariance between the grid index j and the observation k, c p j,k , using where Z represents a certain state variable at the grid point j .The term Z i,j is the contribution of the kth observation to the ith ensemble member of the model state estimated at grid point j .When an observation is available, Eq. ( 5) will be applied to implement CDA for state estimation in a straightforward manner (Zhang et al., 2007;Zhang, 2011a).Although many sophisticated inflation algorithms (e.g., Anderson, 2007Anderson, , 2009;;Li et al., 2009;Miyoshi, 2011) exist for atmosphere data assimilation, the inflation scheme for a coupled model is a new subject due to the multipletimescale nature of the system.Furthermore, trial-and-error experiments show that the usual form of inflation (e.g., only inflate the atmosphere model states or inflate all the model states equally) will lead to the analysis becoming unstable.Thus, in this paper, for simplicity and computational convenience as well as convenience for comparison, no inflation is used in our assimilation experiments, just as in Han et al. (2014).

Perfect and biased twin experiment setups
In this study, a perfect twin experiment framework and a biased twin experiment framework are designed, respectively.In both perfect and biased twin experiments, a "truth" model using the standard parameter values listed in Sect.2.1 is used to generate the "true" solution of the model states and produce the observations sampling the "truth".Starting from the initial condition (0, 1, 0, 0, 0), the "truth" model is firstly integrated forward for 10 000 TUs (i.e., 1000 model years) for sufficient spinup and then integrated forward for another 10 000 TUs to generate the "truth" model states.The observations are produced by sampling the "truth" solution of the model states at an observational interval and superimposing with a white noise simulating the observational errors.As schematically shown in Fig. 1, all the observational intervals used in this study are assumed to be 1 time step (0.01 TU).Although in the real climate system, the oceanic observations are usually available less frequently than those in the atmosphere (that is, the oceanic observation interval is larger than that we set here), for this proof-of-concept study we will set the time interval of the oceanic observations as small as possible.The standard deviations of the observational errors are model also uses the standard parameter values, but starts from different initial conditions using the Gaussian white noises with the same standard deviation as observational errors (2 for X 1 , X 2 , and X 3 , 0.5 for ω, and 0.06 for η) added to the model states at different times during the spinup run to form the ensemble initial conditions for each ensemble filtering data assimilation experiment.Each assimilation experiment is integrated for 10 000 TUs and only the data obtained in the last 5000 TUs are used to conduct error statistics for evaluation.We choose the model states between 9000 and 10 000 TUs during the spinup at an interval of 50 TUs being perturbed to form 20 cases of ensemble initial conditions for each assimilation experiment analyzed in Sects.3 and 4. In this way, we attempt to minimize the dependence of the results of optimal OTWs on ensemble initial states.Then each assimilation experiment will be repeated 20 times starting from these 20 independent ensemble initial conditions and we will analyze the mean value and uncertainty evaluated from these 20 cases.
Then we use the biased experiment setting to simulate the real-world scenario.The biased twin experiment framework is similar to the perfect one except that the assimilation model in the biased twin experiment framework has a systematic discrepancy from the observations.Thus, in the biased twin experiment framework, the parameters included in the assimilation model will have 10 % errors relative to the standard values.The errors in the parameters will be the only model error source.
Figure 1 also illustrates the assimilation update intervals (the assimilation intervals are 5 time steps for atmosphere, 20 time steps for the slab ocean in all assimilation experiments, and 100 time steps for the deep ocean when using the η observations) as well as the length of the OTW, which will be used throughout the study.In addition, the coupling strength between the atmosphere and ocean may have influences on the characteristic variability timescale of each coupled medium, such as on the optimal OTW.We discuss this issue by changing the values of coupling coefficients C 1 and C 2 .In this simple model case, the model stability is sensitive to the coupling coefficient C 1 (Zhang et al., 2012), and changing C 1 only influences the chaotic component, so here we just change C 2 to investigate the impact of the coupling coefficient between the atmosphere and upper ocean on the optimal OTWs.As in Zhang and Anderson (2003), an ensemble size of 20 is applied in all assimilation experiments in this study.

Influence of the OTW on the accuracy of CDA
In order to exhibit the influence of the OTW on the quality of climate analysis, we show three simple assimilation experiments (the time series of ω s absolute errors) in Fig. 2: (1) standard CDA (green) (assimilating the observation rights at the analysis times without any atmospheric/oceanic OTW); (2) OCN-OTW(5) CDA (red) (assimilating all 11 ocean observations in an oceanic OTW with a half-width of 5, defined as the length of the OTW hereafter, but no atmospheric OTW is considered), and (3) OCN-OTW(100) CDA (blue) (assimilating all 201 ocean observations in an oceanic OTW with a length of 100, but no atmospheric OTW is considered).The three assimilation experiments above are all conducted using the perfect model setting (all the parameters use their standard values) and the univariate adjustment scheme.The atmospheric and oceanic update intervals are 0.05 and 0.2 TU, respectively.While the standard CDA does not use the atmospheric and oceanic OTWs and only assimilates the observations right at the analysis time, the OCN-OTW CDA incorporates all the valid observations collected in the oceanic OTW.All three assimilation experiments above do not use an atmospheric OTW.
From Fig. 2, we can see that a small OCN-OTW (a total of 11 observations in the oceanic OTW) can make a much better ocean analysis than the standard CDA (comparing the red line with the green line).We can understand that this is because an OTW can provide more observational information, thus enhancing the observational constraint so as to improve the accuracy of climate analysis.However, comparing the blue line to the green/red line, it is clear that an overly large OTW degrades the quality of the ocean analysis.The results of these simple assimilation experiments tell us that, if an ap- propriate OTW is used, we can gain optimal climate analysis.How can we determine such an optimal OTW? Next, starting from analyzing the characteristic variability of each coupled medium, we will discuss the methodology of how to determine an optimal OTW for each medium in a coupled climate system.

The timescale of characteristic variability and an optimal OTW
The key to improving the accuracy of climate analysis in CDA is by accurately recovering the characteristic variability of different media in the coupled system.Thus we can assume that the length of an optimal OTW for each medium will have some relationship with the corresponding characteristic variability timescale.Then, we should first analyze the timescale of characteristic variability in each medium.Figure 3 presents the power spectrum of X 2 ω and η based on the model states with a 4800 TU length (totally 480 000 data) after the spinup described in Sect.2.3.From Fig. 3, we learned that in this simple model, the characteristic variability timescales of atmosphere (X 2 ), upper ocean (ω), and deep ocean (η) are about 1-2 TUs (1-2 model months), 50-100 TUs (5-10 model years), and 500 TUs (5 model decades), respectively.That is, the characteristic variability timescale of the slab ocean is much larger than that of the atmosphere, but smaller than that of the deep ocean.
An optimal OTW aims to provide maximal observational information that best samples the characteristic variability of that medium during the data blending process.Thus the length of the optimal OTW should be smaller than the corresponding characteristic variability timescale, which means that the optimal OTW in the atmosphere must be much smaller than 1 TU (100 time steps), and in the ocean, the optimal OTW must be much smaller than 50 TUs (5000 time steps).If we take observations for η, the optimal OTW for η must be much smaller than 500 TUs (50 000 time steps).From Fig. 3, we also see that the characteristic variability timescales of different coupled media are a little larger than the corresponding ones set in Eq. ( 1).This is owing to the strong nonlinearity and smoothness of the fourthorder Runge-Kutta time-differencing scheme that prolongs the characteristic variability timescales of the simple coupled model.But they do not change the essence of the problem we address in this study.Given different timescales of characteristic variability in different media, in the following section we will further detect the optimal OTWs based on the corresponding characteristic variability timescales and examine their impact on the quality of climate analysis in CDA.

Detection of the optimal observational time window
In this section, with the perfect model framework described in Sect.2.3, we first conduct a series of CDA experiments with different ATM-OTWs and different OCN-OTWs to detect the optimal OTW for each medium.The assimilation scheme is the simple univariate adjustment scheme serving as a proof-of-concept study.To eliminate the dependency of results on initial states, each experiment is repeated 20 times starting from the 20 independent initial conditions described in Sect.2.2.Then the mean value and the spread of 20 cases of RMSEs are plotted in Fig. 4.
Figure 4a shows that the optimal ATM-OTW is 1; i.e., the optimal ATM-OTW includes only three atmosphere observations, with which the assimilation produces the lowest RMSE of the atmosphere and the smallest spread.(In this study each assimilation experiment will be repeated 20 times starting from 20 different independent initial ensemble conditions.Here the spread just represents the standard deviation of these 20 cases' results.Thus it will be smallest when using the optimal OTW.)In these experiments for detecting the optimal ATM-OTW, the ocean assimilation is kept as the standard setting (i.e., no OTW, 0.2 TU update interval).Then we keep the ATM-OTW as 1 and change the length of OCN-OTW to produce Fig. 4b.
From Fig. 4b, we can see that the optimal OCN-OTW is about 10 (i.e., each OTW includes a total of 21 observations), with which the lowest ω-RMSE and the smallest RMSE spread are produced.Compared to the case of standard CDA (denoted as CDA_NOTW), the uses of optimal ATM-OTW and OCN-OTW make the RMSEs of X 1,2,3 and ω significantly reduced.When the RMSE of ω has a distinguishable sensitive variation with respect to OCN-OTWs, the RMSE of X does not show such a sensitivity to the optimal Nonlin.Processes Geophys., 24, 681-694, 2017  Variations of root mean square errors (RMSEs) of (a) "atmospheric" states X 1,2,3 (namely, the average of X 1 , X 2 , and X 3 RMSEs) in the space of ATM-OTW length when the "oceanic" state (ω) only uses a single observation at the assimilation time; (b) "upper ocean" state (ω) in the space of OCN-OTW length when the ATM-OTW is fixed at 1 as shown in panel (a) (1 for the ATM-OTW, i.e., three observations in each window; see the caption of Fig. 1) but the OCN-OTW (for ω) is varying and (c) "deep ocean" state (η) in the space of η-OTW length when the "deep ocean" observations are assumed to be valid and the ATM-OTW and OCN-OTW are fixed as 1 and 10, respectively.The experiments are conducted in a perfect model setting with a simple univariate adjustment scheme.The red lines are the 20-case mean, each using different initial conditions taken from different periods in the control integration (see the description in Sect.2.2), and the blue lines represent the upper/lower bounds (mean ± standard deviations) of the RMSEs.An OTW with the length of 0 represents only assimilation of the observation at the assimilation time (i.e., with no OTW, dashed-black lines).The RMSE values of the control case (no observational constraint, called CTL) are marked in the parentheses.
OCN-OTW.(Because we just choose the optimal OCN-OTW from the figure of ω-RMSE.Thus in this study the variation of X-RMSEs in the OCN-OTW space is not shown.)This means that in this simple system, due to the strong nonlinearity and chaotic nature of the "atmosphere", the improved accuracy for ω from optimal observational constraint is not sufficient to impact the "atmosphere" (this point will be expanded in Sect.4.3).Similarly to the characteristic variability timescale of the slab ocean vs. that of the "atmosphere" shown by Fig. 3, the optimal OCN-OTW is much larger than that of ATM-OTW.
To further understand the relationship between the optimal OTW and characteristic variability timescale, we also examine the η-RMSEs in the space of η-OTWs.The assimilation interval of the pycnocline depth is set as 1 TU (100 time steps), which is much larger than that of the slab ocean.When we change the η-OTW, the optimal ATM-OTW and OCN-OTW detected from Fig. 4a and b are used.As shown in Fig. 4c, the optimal η-OTW is about 100 (i.e., a total of 201 observations), which is much larger than that of OCN-OTW and smaller than the characteristic variability timescale of the deep ocean pycnocline depth.With the optimal η-OTW, the RMSE of η is reduced by about 77.4 % from the level of the CDA_NOTW.
We also check the variation of the 20-case mean ensemble spread in the space of OTWs as shown in Fig. 5.The mean and standard deviation of the ensemble spreads of X 1,2,3 and ω (the uncertainty of the state estimation in each assimilation experiment is shown as the blue shadow in Fig. 5) gradually decrease when the ATM-OTW and OCN-OTW become larger.When the OTWs are set too large (here the ATM-OTW and OCN-OTW are greater than 20 and 250, respectively), the ensemble spreads of X 1,2,3 and ω decrease dramatically.This is owing to the fact that when we increase the length of the OTW, more observations will be included in the OTW and then assimilated into the corresponding model states, which can function as a smoother.The longer the lengths of OTWs are, the stronger the smoother will be.Also, under this circumstance, the overly strong smoother will distort the characteristic variability of the model states, which explains the blue line of Fig. 2. From Figs. 4 and 5, we can see that the mean of the ensemble spread is significantly smaller than that of the corresponding RMSE.It is owing to the fact that no inflation scheme is applied in this study.And the statistics for evaluation are conducted from the data obtained in the last 5000 TUs.Thus after the first 5000 TUs' assimilation in each assimilation experiment, the ensemble spreads of model states have been greatly reduced due to no inflation.Then the mean ensemble spread is significantly smaller than the mean RMSE.
To understand the essence of optimal OTWs, we show the auto-correlation for each model state and mark the time correlation coefficients at the timescales of optimal OTWs for X 2 (a), ω (b), and η (c) detected from Figs. 4 to 6.The result is the mean of 20 cases.In each case, the number of data are 10 000 steps (100 TUs), which are chosen from the period of 5000 to 9000 TUs in the truth run after spinup.From Fig. 6 we can see that all auto-correlations at the optimal OTW length are located at around 0.995.This means that the observations included in an optimal OTW are extremely highly correlated with the model state at the analysis time.This can be understood since in this sequential assimilation scheme all the observations included in an OTW are assumed to be sampled at the analysis time so that the difference among them must be in a negligible range.Under such a circumstance, the optimal OTWs provide maximal observational information that best fits the characteristic variability and minimizes the analysis error.

Influences of realistic assimilation scenarios on optimal OTWs
In this section, we first show the impact of the multi-variate adjustment scheme on the optimal OTWs in a perfect model setting.Then we discuss the influence of model bias through a biased model framework.We will also investigate the impact of coupling strength on the optimal OTWs.

Influence of multi-variate adjustment on optimal OTWs
While the experiments with the univariate adjustment scheme provide us with a direct understanding of the influence of the OTWs on CDA, we want to check whether or not it also applies to the multi-variate adjustment scheme.So we repeat the experiments described in Sect. 3 but with the multi-variate adjustment scheme.The results are shown in Fig. 7.Here the multi-variate adjustment scheme is only limited to the atmospheric observations (i.e., only the crosscovariances among X 1 , X 2 , and X 3 are used; as indicated in Han et al., 2013, the multi-variate adjustment scheme using the coupling cross-covariance between different coupled media involves complex scale interactions and may complicate In each case, an independent section (each has 10 000 data of the state -100 TUs with the interval of 0.01 TU) is used to evaluate the lag correlation coefficient.The 20 independent sections are taken from the model states apart each 200 TUs between 5000 and 9000 TUs integrations after the spinup of 10 000 TUs from the initial condition (0, 1, 0, 0, 0). the investigation of the problem we are addressing here).The results shown in Fig. 7 are similar to that in Fig. 4, suggesting the multi-variate adjustment scheme has little influence on the optimal OTWs, since it does not change the characteristic variability timescales (especially in this simple model).
The perfect experiment framework provides a direct guideline for the relationship between the optimal OTW and the corresponding characteristic variability timescale.However, in reality, the numerical model has errors and is biased with the observation.It is as necessary to investigate the influence of model bias on optimal OTWs as on the quality of CDA.

Influence of model bias on optimal OTWs
With the biased model experiment framework described in Sect.2.3, we repeat all the experiments above for detection of the optimal OTWs.The results are shown in Fig. 8. Compared to the results in the perfect model setting, the results in the biased model setting have two differences.First, the optimal ATM-OTW and OCN-OTW are larger than their counterparts in the perfect model setting, becoming 3 and 20 (that is, the total observations are 7 and 41, respectively).Second, the RMSE curves in the space of OTWs show more concavity and sensitive variation.This is more distinguishable in the curve of ω-RMSEs in the OCN-OTW space.All these phenomena can be explained by the influence of model bias on the assimilation quality.On the one hand, due to the existence of model bias, the assimilation not only needs observations to fit the observed variability, but also needs observations to reduce the mean discrepancy between the model and observations.This requires stronger observational constraints.An optimal OTW that makes the smallest RMSE of model states must include more observed data.On the other hand, the forecast ensemble in a biased model underestimates the forecast error, which results in the EAKF underweighting the observations.Therefore the optimal OTWs are larger than those in the perfect experiment case in that the observations included in the optimal OTWs will be assimilated for multiple times, which results in an improvement of filter performance.The test experiment for the optimal η-OTW is also consistent with this point (in Fig. 8c): the optimal η-OTW in the biased model setting is larger than that in the perfect model setting.Then we also investigate the influence of OTWs on the quality of CDA with the multi-variate adjustment scheme in the biased experiment framework (not shown here).The results are the same as the perfect model setting case; i.e., the multi-variate adjustment scheme does not change the optimal OTWs.
Comparing the results from two experiment frameworks, we can see that regardless of perfect or biased model setting used in the assimilation experiments, the optimal OTW must be associated with the corresponding characteristic variability timescale in the medium.It is clear that while using observations in an OTW increases observational information, an overly large OTW can distort the characteristic variability of coupled media during the information blending process.Therefore choosing an optimal OTW that is much smaller than the medium's characteristic variability timescale is very important.The simple model results suggest that the length of an optimal OTW is about 1-5 % of the medium characteristic timescale, with which characteristic variability of the medium can be retrieved most accurately.
In this study, the OTW validates the observations in a time window to the analysis time and all the observations included in the OTW are sequentially assimilated with their original error scales.Another general approach is to assimilate the average of the observations included in the OTW, but the observational errors decrease as 1/

√
N of their original error scales (N represents the number of observations included in the OTW).From the comparison of these two methods (not shown), we can see that the results obtained by them are almost the same.From the perspective of the calculation pro-cess of the EAKF method, owing to no inflation scheme being used, after many assimilation steps the ensemble spreads of the model states have been greatly reduced and are significantly smaller than the corresponding observational error scales.And the prior ensemble member will be very close to the prior ensemble mean.Thus the analysis adjustments obtained by these two methods will be almost the same.It is worth mentioning that although the resulting RMSEs obtained by these two assimilation schemes will be different when using the suitable inflation schemes, the lengths of the optimal OTWs are still the same and the essence of this study is still firm and does not change.
Also, among the above assimilation experiments in this study, we have not considered the temporal offset induced by the difference between the time of observations in the OTW and the analysis time.Here we can use the de-correlation coefficients to weight the observations included in the OTW and avoid overweighting them.The comparison of these two assimilation approaches (non-weighted and weighted) has been conducted (the results are not shown).From the comparison we learn that the lengths of the optimal OTWs obtained by these two assimilation schemes are similar, except that the RMSEs in the weighted observation experiment will be lower than that in the non-weighted one when using longer OTWs (when the length of ATM-OTW is greater than 4 and/or that of the OCN-OTW is larger than 50).This is owing to the high correlation between the observation included in the optimal OTWs and model states at the analysis time (exceeding 0.995).Thus the influence of the temporal offset can be ignored and the results obtained by these two schemes will almost be the same when using the shorter OTWs.When we use the longer ones, the correlation will decrease and the influence of the temporal offsets will be obvious in that the results of the weighted observation experiment will be better.For the CDA systems in the CGCMs, owing to the complex physics and dynamics, the influence of the time offsets will be obvious and the weights of the observations will be very necessary.But from this simple model case, we can see that regardless of using the weighted observations, the relationship between the characteristic variability timescales and the optimal OTWs will be robust, and the essence of this study is established.

Influence of coupling strength on optimal OTWs
Changing the coupling strength (controlled by the coupling coefficients C 1 and C 2 in this case) between the atmosphere and upper ocean may have some influence on the characteristic variability timescales of coupled media, as on the optimal OTWs.Test experiments show that changing the coupling coefficient C 1 has little influence on the characteristic variability timescales of X 1,2,3 and ω.This is because the characteristic timescale of X is determined by the chaotic nature of the Lorenz equations, not by the oceanic forcing associated with the coupling coefficient C 1 .Therefore, here we just change C 2 to investigate the coupling coefficient between the atmosphere and upper ocean on the optimal OTW of ω.Setting the values of C 2 as 1.5, 1.25, 1.0, 0.8, 0.5, and 0.1 and keeping C 1 as 0.1, we repeat all the biased CDA experiments with the multi-variate adjustment scheme.The results are shown in Fig. 9, which presents the power spectrum of X 2 and ω (panels a and b) of the six cases above based on the model states between 5000 and 9800 TUs, as well as the time series of model states between 5000 and 5100 TUs (panel c) after the spinup described in Sect.2.3.We can see that changing C 2 does not influence the characteristic variability timescale of the atmosphere, but strongly influences the variability of the slab ocean.From the equation of ω, the characteristic variability timescale of ω is determined by the combination of the atmospheric forcing and the periodic ex-ternal forcing.When C 2 is small, the forcing of atmosphere to ocean is weak, and then the periodic external forcing plays a dominant role in determining the characteristic variability timescale of the ocean component.
Then we examine the difference in the optimal OTW of ω in the six cases above, as shown in Fig. 10.The results show that changing C 2 does not have any influence on the optimal ATM-OTW (not shown).From panels a and b we can see that when C 2 is smaller, the optimal OCN-OTW is larger.This can be explained by the increasing role of the periodic external forcing in determining variability of the slab ocean, for which data assimilation needs more observational information to recover the periodic variation of ω, determined by the timescale defined by S pd (10 TU).When C 2 is larger than 1.0, changing it has little influence on the characteristic variability of the ω, as on the optimal OCN-OTW.
On the one hand, these experiments can further illustrate the idea that a close relationship between the length of the optimal OTW and the corresponding characteristic variability timescale exists.On the other hand, for a realistic CDA system, the coupling physics could be very complicated and affected by many factors.The results of this simple model give the insights that when determining the length of the optimal OTWs for a realistic CDA system, we can only consider such factors that have an obvious influence on the characteristic variability timescales.In this way, the process of determining the optimal OTWs in a realistic CDA system can be greatly simplified and make it possible to apply the method of using the optimal OTWs to the realistic CDA system.

Summary and discussions
With a simple conceptual climate model and the EAKF method, the impact of OTWs on the quality of CDA has been investigated in this study.This simple conceptual coupled model consists of a synoptic atmosphere (Lorenz, 1963) and seasonal-interannual slab upper ocean (Zhang et al., 2012) coupling with a decadal deep ocean (Zhang, 2011a, b), and reasonably simulates the typical interactions between multiple timescale components in the climate system.Determined from the characteristic variability timescale in each coupled medium, an optimal OTW provides maximal observational information to best fit the characteristic variability of the medium during the data blending process.With correct scale interactions within the coupled system, CDA can recover the climate signals most accurately by incorporating all observations in the optimal OTWs into the coupled model, although in an idealized and simple model circumstance, the conclusion addressing the best fitting characteristic variability in each medium with the optimal OTW is comprehensive and therefore provides a guideline for improving climate analysis and prediction initialization when real observations are assimilated into a CGCM.For example, as learned from the simple model results, we may consider improving the quality of climate analysis and prediction initialization by accurately recovering some important characteristic variability in the atmosphere (sub-diurnal variations, for instance) and ocean (diurnal cycle in the tropical oceans, for instance).
However, the current work can only serve as a proof-ofconcept study.Although CDA with the optimal OTWs has shown promising improvement in this simple model, serious challenges still exist for detecting optimal OTWs in the real world with a CGCM for improving climate analysis and prediction.First, the characteristic variability timescales in different media of the real world are complex, and great challenges remain to identify the characteristic variability of the different component models and the real atmosphere and upper and deep ocean, which need to be further studied.Also, in a real ocean model, the upper and deep ocean is inseparable, which bring some troubles in using different OTWs for different parts of the same ocean model.Second, due to model biases, characteristic variability in a CGCM may be different from the real world.The combination of variability of the real world and that of the model may further complicate the problem.Therefore, model bias and its influence on model variability need to be thoroughly analyzed before an optimal OTW is determined.Thirdly, the coupling physics between different coupled components are very complicated and are impacted by many factors for a realistic CDA system.Even though we only consider the factors which will obviously impact the characteristic variability timescales when determining the length of OTWs for different coupled components, it remains a heavy workload.In addition, in this study we assume that all observations in the OTWs have equal weights to contribute to the observational constraint.In the real observation case, the observation far away from the assimilation time should have less contribution to the state estimation at the assimilation time.How to take the time correlation into account in a sequential algorithm needs to be studied before implementing optimal OTWs in the assimilation with CGCM and real observations.

Figure 1 .
Figure1.The schematic for the assimilation interval, the length of the OTW, as well as the observational interval in terms of the model integration time step.Here L represents the time steps at one side of the OTW.For example, OCN-OTW (L) in the content stands for an ocean observational time window with total observations of 2L + 1.

Figure 2 .
Figure2.Time series of the absolute errors of the slab ocean variable (ω) in three assimilation experiments based on the model states between 9100 and 9110 TUs assimilation results in the perfect model experiment framework with the univariate adjustment scheme.Green -CDA control with the standard update intervals of 0.05 TU for X 1,2,3 and 0.2 TU for ω; red -CDA with an ocean observational time window (OCN-OTW) of five time steps (OCN-OTW (5)); blue -CDA with OCN-OTW (100).
Figure4.Variations of root mean square errors (RMSEs) of (a) "atmospheric" states X 1,2,3 (namely, the average of X 1 , X 2 , and X 3 RMSEs) in the space of ATM-OTW length when the "oceanic" state (ω) only uses a single observation at the assimilation time; (b) "upper ocean" state (ω) in the space of OCN-OTW length when the ATM-OTW is fixed at 1 as shown in panel (a) (1 for the ATM-OTW, i.e., three observations in each window; see the caption of Fig.1) but the OCN-OTW (for ω) is varying and (c) "deep ocean" state (η) in the space of η-OTW length when the "deep ocean" observations are assumed to be valid and the ATM-OTW and OCN-OTW are fixed as 1 and 10, respectively.The experiments are conducted in a perfect model setting with a simple univariate adjustment scheme.The red lines are the 20-case mean, each using different initial conditions taken from different periods in the control integration (see the description in Sect.2.2), and the blue lines represent the upper/lower bounds (mean ± standard deviations) of the RMSEs.An OTW with the length of 0 represents only assimilation of the observation at the assimilation time (i.e., with no OTW, dashed-black lines).The RMSE values of the control case (no observational constraint, called CTL) are marked in the parentheses.

Figure 5 .
Figure 5. Same as panels (a) and (b) in Fig. 4 but for the variation of ensemble spreads of the model states.In panel (b) the optimal ATM-OTW is also set as 1.The area between the lower and upper bounds (blue) represents the range evaluated from the 20 cases.And the blue shadow below the ensemble spreads represents the range of the uncertainty of state estimation in each assimilation experiment.

Figure 6 .
Figure6.The auto-correlation coefficient of (a) X 2 (b) ω, and (c) η in the space of lag times are marked by corresponding time correlation coefficients at the timescale (L) of optimal OTWs as detected by Fig.4for different media (the black-dashed lines).What are shown are the means of 20 cases.In each case, an independent section (each has 10 000 data of the state -100 TUs with the interval of 0.01 TU) is used to evaluate the lag correlation coefficient.The 20 independent sections are taken from the model states apart each 200 TUs between 5000 and 9000 TUs integrations after the spinup of 10 000 TUs from the initial condition (0, 1, 0, 0, 0).

Figure 7 .
Figure7.Same as Fig.4but using a multi-variate adjustment scheme.In panel (b) the optimal ATM-OTW is also set as 1.

Figure 8 .
Figure8.Same as Fig.4but using the biased model setting.In panel (b) the optimal ATM-OTW is set as 3.And in panel (c) the optimal ATM-OTW and OCN-OTW are kept as 3 and 20, when the "deep ocean" observations are assumed to be valid.

Figure 9 .Figure 10 .
Figure9.The power spectrum of (a) X 2 and (b) ω based on the model states between 5000 and 9800 TUs integrations after the spinup which integrates for 10 000 TUs from the initial condition (0, 1, 0, 0, 0) with different coupling strengths (C 2 is set as 1.5, 1.25, 1.0, 0.8, 0.5, and 0.1, while C 1 remains as 0.1).Panel (c) shows the time series of the model state ω between 5000 and 5100 TUs integrations corresponding to the six cases.