Mitigation of coupled model biases induced by dynamical core misfitting through parameter optimization : simulation with a simple pycnocline prediction model

Imperfect dynamical core is an important source of model biases that adversely impact on the model simulation and predictability of a coupled system. With a simple pycnocline prediction model, in this study, we show the mitigation of model biases through parameter optimization when the assimilation model consists of a “biased” timedifferencing. Here, the “biased” time-differencing is defined by a different time-differencing scheme from the “truth” model that is used to produce “observations”, which generates different mean values, climatology and variability of the assimilation model from the “truth” model. A series of assimilation experiments is performed to explore the impact of parameter optimization on model bias mitigation and climate estimation, as well as the role of different media parameter estimations. While the stochastic “physics” implemented by perturbing parameters can enhance the ensemble spread significantly and improve the representation of the model ensemble, signal-enhanced parameter estimation is able to mitigate the model biases on mean values and climatology, thus further improving the accuracy of estimated climate states, especially for the low-frequency signals. In addition, in a multiple timescale coupled system, parameters pertinent to low-frequency components have more impact on climate signals. Results also suggest that deep ocean observations may be indispensable for improving the accuracy of climate estimation, especially for low-frequency signals.


Introduction
Imperfect dynamical core, empirical physical schemes and improper parameter values are several sources of couple model bias (Zhang et al., 2012).Simulated climate by a coupled model often tends to drift away from the real world due to the existence of model bias (Collins et al., 2006;Delworth et al., 2006;Smith et al., 2007).However, it is quite difficult to improve the model simulation and forecast capability through using observations to correct the dynamical core and physics that are "built-in" in the coupled model.One expects that the parameter optimization can partly compensate for the deficiencies of both numerics and physics of a coupled model and improve the model performance to some degree.
To constrain model biases and improve the quality of climate estimation and prediction, Zhang et al. (2012) designed a coupled data assimilation scheme with what these authors called "enhancive" parameter correction (DAEPC) based on an ensemble Kalman filter with the adjustment idea (Anderson, 2001).With the DAEPC algorithm, Zhang (2011a) investigated the impact of observation-optimized model parameters on decadal predictions with a simple pycnocline prediction model.Then in a follow-up study (Zhang, 2011b), the author also investigated the impacts of coupled model initial shocks and state-parameter optimization on climate predictions using this simple model.Results show that model parameter optimization with observations can effectively mitigate the model bias, thus constraining the model drift in long timescale predictions.Wu et al. (2012) further introduced a geographic dependent parameter optimization (GPO) scheme to increase the signal-to-noise ratio of the background error covariance in parameter estimation, and examined the impact of this new scheme on climate estimation and prediction using an intermediate coupled model within a perfect model framework (Wu et al., 2013).Recently, Zhang et al. (2013b) investigated the impact of parameter estimation on climate estimation and prediction in an intermediate coupled model with biased physics within a biased twin experiment framework, which indicates that the adverse impact of biased physical schemes in a coupled model on climate estimation and prediction can be compensated partly by optimizing the most sensitive parameters employed in the physical schemes.The impact of estimated parameters on the behavior of model simulation has also been examined (Zhang et al., 2013a), with results showing that biased climate simulated by "biased" physics in that intermediate coupled model can be well corrected through parameter estimation.
While coupled model parameter estimation has shown a great potential to improve the quality of climate estimation and prediction as well as model simulation, the impact of imperfect dynamical cores such as imperfect numerical schemes has not been examined yet.To address the question, based on the DAEPC algorithm (Zhang et al., 2012), we study how to mitigate coupled model bias induced by imperfect time-differencing schemes through parameter optimization.Here we use the simple pycnocline prediction model described by Zhang (2011a) to investigate this issue within a biased twin experiment framework.Under such circumstances, one model simulation that uses a leap-frog time-differencing scheme with a Robert-Asselin time filter (Robert, 1969;Asselin, 1972) is treated as a "truth" that is sampled to produce "observations."Then the "observations" are assimilated into the assimilation model that uses the fourth-order Runge-Kutta time-differencing scheme.The degree to which the assimilation result recovers the truth is an assessment of the impact of parameter optimization on the climate estimation with a "biased" time-differencing.
The paper is organized as follows.After describing the simple pycnocline prediction model and the method of ensemble coupled data assimilation for parameter estimation, two different time-differencing schemes are introduced and the setting of the biased twin experiment framework is discussed in Sect. 2. Sections 3 and 4 investigate the impact of parameter optimization on climate estimation and the impact of parameter estimation in different media on model bias mitigation, respectively.Summary and discussions are given in Sect. 5.

The model
To address the fundamental issues raised in Sect. 1 clearly, we employ a simple decadal prediction model developed by Zhang (2011a).The model consists of a conceptual atmosphere-ocean coupled model that couples three Lorenz chaotic atmosphere variables, x 1 , x 2 , and x 3 (Lorenz, 1963), to a slab-ocean variable w and a simple pycnocline predictive model (Gnanadesikan, 1999).The governing equations with all quantities being given in non-dimensional units are where the five model variables represent the atmosphere (x 1 , x 2 , and x 3 ) and the ocean (w for the slab ocean, η for the deep ocean pycnocline).A dot above the variable denotes time tendency.For the equation of w, O m is the heat capacity of the ocean, and O d denotes the damping coefficient of the slab ocean variable w.An important feature of w is that it must have a much slower timescale than the atmosphere, which needs a much larger heat capacity than the damping rate, that is O m O d .For example, the values of (10, 1) for (O m , O d ) define the oceanic timescale as ∼ O(10), 10 times the atmospheric timescale ∼ O(1).The parameters S m and S s define the magnitudes of the annual mean and seasonal cycle of the external forcings.S pd is chosen as 10 so that the period of the forcing is comparable with the oceanic timescale, defining the timescale of the model seasonal cycle.The coupling between the fast atmosphere and the slow ocean is realized by choosing the values of the coupling coefficients c 1 and c 2 , with c 1 representing the upper oceanic forcing on the atmosphere, and c 2 representing the atmospheric forcing on the upper ocean.In addition, c 3 and c 4 denote the linear forcing of the deep ocean and the nonlinear interaction of the upper and deep oceans.In the pycnocline model, η represents the anomaly of ocean pycnocline depth, and its tendency equation is derived from the two-term balance model of the zonal-time mean pycnocline (Gnanadesikan, 1999). is a constant of proportionality.The ratio of and O d determines the timescale of variations of η, for example, a value of 100 for defining 10 "seasonal" cycles of w (a model decade) as the typical timescale of η variability.To simulate the effects of the nonlinear advection in the upper and deep oceans, the nonlinear terms are introduced into w and η equations.c 5 and c 6 represent the linear forcing of the upper ocean and the nonlinear interaction of the upper and deep oceans.Following Zhang (2011a), the values of 15 model parameters (σ , κ, b, c 1 , c 2 , O m , O d , S m , S s , S pd , , c 3 , c 4 , c 5 , c 6 ) are set as (9.95, 28, 8/3, 10 −1 , 1, 1, 10, 10, 1, 10, 100, 10 −2 , 10 −2 , 1, 10 −3 ), in which the parameters of σ , κ, and b in the Lorenz atmosphere keep their standard values to sustain the chaotic nature of the "atmosphere."Zhang (2011b) illustrated that, given the model parameters prescribed above, the built simple coupled model can effectively simulate a fundamental feature of the real world climate system in which different timescales interact with each other to develop climate signals.For example, the transient atmospheric attractor, the slow upper ocean and the even-slower deep ocean interact to produce synoptic decadal timescale signals.

Ensemble coupled data assimilation for parameter estimation
The coupled data assimilation scheme with "enhancive" parameter correction (DAEPC) (Zhang et al., 2012) mentioned above is employed to perform the model state and parameter optimization, which is a modification of the standard data assimilation with adaptive parameter estimation (e.g., Kulhavy, 1993;Tao, 2003).Some details of the DAEPC algorithm are given below to make it easy to follow.Based on a two-step ensemble adjustment Kalman filter (EAKF; Anderson, 2003;Zhang and Anderson, 2003), the observational increment for the ith ensemble member produced by the kth observation, y o i,k , is computed firstly following Zhang et al. (2007) as where the first two terms on the right-hand side represent the shift of ensemble mean and the third term is the adjustment of ensemble spread.y i,k is the ith prior ensemble member of the kth observation.y k is the model estimate ensemble for observation y o k .An overbar represents the ensemble mean.r(y k , y o k ) is the ratio of the model ensemble standard deviation and the observational error standard deviation, that is, σ y k /σ y o k .In the second step of EAKF, the observational increment is projected onto the corresponding model variables using a uniform linear regression formula as where x i,k is the contribution of the kth observation to the ith ensemble member of each model variable x. cov(x, y k ) denotes the error covariance between the prior ensemble of x and the model-estimated observation ensemble of y k .The observational increment is also projected onto the parameters being optimized using the uniform linear regression formula as where β i,k is the contribution of the kth observation to the ith ensemble member of the parameter being optimized, called β. cov(β, y k ) denotes the error covariance between the prior ensemble of β and the model-estimated observation ensemble of y k , and is calculated as where N is the number of the ensemble member.β i is the ith ensemble member of each parameter being optimized.An overbar represents the ensemble mean.σ β is the prior standard deviation of the parameter being optimized.
Unlike the model state variables, model parameters do not have any dynamically supported internal variability in general.Therefore, the successfulness of parameter estimation entirely depends on the accuracy of the state-parameter covariance in Eq. ( 4).Parameter estimation is activated after state estimation reaches quasi-equilibrium where the uncertainty of model states is sufficiently constrained by observations so that the state-parameter covariance is signal dominant.Otherwise, the parameters being optimized are likely to be deteriorated by the noised state-parameter covariance in Eq. ( 4).
In addition, the inflation scheme for the DAEPC algorithm follows Zhang et al. (2012), which is formulated as where β l and βl represent the prior and the inflated ensemble of the lth parameter.σ l,t and σ l,0 are the prior spreads of β l at time t and the initial time.α 0 is a constant tuned by a trial-and-error procedure.σ l is the sensitivity of the model state with regard to β l .The overbar represents the ensemble mean.Equation ( 6) indicates that if the prior spread of β l is less than α 0 /σ l times the initial spread, it will be enlarged to this amount.

Two different time-differencing schemes
Here, we introduce two time-differencing schemes.The first one is the leap-frog (LF) time-differencing with a Robert-Asselin time filter, which has the form where ϕ represents state variables (x 1 , x 2 , x 3 , w, and η) in Eq. (1).t is the time interval.F is the right term of state variables in Eq. ( 1).An overbar denotes a time-filtering value.The time-filtering coefficient is set as γ = 0.25 in this study, following Zhang (2011b).
The second one is the fourth-order Runge-Kutta (RK4) time-differencing scheme, which can be described as where k 0 -k 3 represent four time levels.

Model bias induced by different time-differencing schemes
Starting from the initial conditions (x 1 , x 2 , x 3 , w, η) = (0, 1, 0, 0, 0), the model is run for 10 4 non-dimensional time units (TUs, 1 TU = 100 time steps given t = 0.01) with the LF time-differencing scheme and the RK4 time-differencing scheme respectively, described in Sect.2.3.Figure 1a, b,  and c show the time series of x 2 in the first 100 TUs, w in the first 10 3 TUs, and η in 10 4 TUs obtained from the LF timedifferencing scheme (see the red line in Fig. 1) and from the RK4 time-differencing scheme (see the black line in Fig. 1), respectively.From Fig. 1a, the two lines of x 2 are almost coincident in the first 5 TUs and separate gradually then to follow different paths, which indicates that the difference originating in these two time-differencing schemes can generate a dramatic effect due to the strong nonlinear nature of the climate system.Due to the different time-differencing schemes and the coupling, the low-frequency signals are also affected significantly (see w in Fig. 1b and η in Fig. 1c, respectively).
Especially the mean value of η derived from the RK4 timedifferencing scheme (the black line in Fig. 1c) is larger than that derived from the LF time-differencing scheme (the red line in Fig. 1c) after a sufficient spin-up of η.Therefore, for the high-frequency signal such as the strong nonlinear atmosphere, the different time-differencing schemes can result in a difference in phase, while for a low-frequency signal such as the deep ocean, the different time-differencing schemes can result in a difference not only in phase, but also in amplitude.The plots in Fig. 2 also reveal the same fact, in which the variation of x 2 in the space of x 3 (Fig. 2a) and the variation of η in the space of w (Fig. 2b) are shown for those derived from both the LF (red circle) and the RK4 (black circle) time-differencing schemes.We can see from Fig. 2a that both projections on the x 2 -x 3 plane from the two schemes lie in the similar equilibrated positions with two attractor lobs even though they are not coincident with each other in details.However, projections on the η-w plane from the two time-differencing schemes have different positions after reaching equilibrium (Fig. 2b). Figure 3 shows the power spectra of w (Fig. 3a) and η (Fig. 3b) based on the model results between 10 3 and 10 4 TUs derived from the LF time-differencing scheme (red line) and from the RK4 time-  differencing scheme (black line), respectively.According to the governing equation of w, the seasonal cycle is defined by 10 TUs, and therefore a model year (decade) is 10 (100) TUs.
It can be seen from Fig. 3 that the power spectra from the The above discussions show that both the time mean and the variability of the model simulation are affected strongly by the type of time-stepping scheme used.We expect to constrain the model drift caused by biased time-differencing through the DAEPC algorithm that tunes model parameters optimally according to the observational information.Next, we will design a biased twin experiment to investigate the impact of parameter optimization on the mitigation of coupled model bias induced by imperfect time-differencing schemes.

Biased twin experiment setup
A biased twin experiment is designed, which assumes that the imperfect time-differencing scheme is the only source of model biases.Here we define a "truth" model in which the LF time-differencing scheme is used.Starting from the initial conditions described in Sect.2.4, after the truth model is spun up for 10 3 TUs, a time series of true states is generated over a period of 10 4 TUs."Observations" are produced by adding a Gaussian white noise to the true values every 5 time steps for x 1,2,3 and every 20 time steps for w, respectively.These observational frequencies simulate the feature of the real climate observing system in which the atmospheric observations are available more frequently than the ocean.At the same time, these observational frequencies are considered as the assimilation frequencies that the observations are assimilated into the model.To simulate the lack of deep ocean measurements in the real observing system, no observation is available for η.According to Zhang (2011a), the standard deviations of the observational errors are 2 for x 1,2,3 and 0.5 for w, respectively, which are derived from a long-time model free run.
The assimilation model uses the RK4 time-differencing scheme.A Gaussian white noise with the same standard deviation as the observational error is added to the atmospheric variable of x 2 at the end of spin-up (with the same integration period as the truth model) to form the ensemble initial conditions from which the ensemble filtering data assimilation starts.The total data assimilation period is 10 4 TUs, and parameter optimization is started after 10 3 TUs when state estimation reaches its quasi-equilibrium.
Starting from the ensemble initial conditions created above, three assimilation experiments are conducted.First, the experiment of state estimation only (SEO) that only the model states are estimated by assimilating observations into the model.Second, SEO with perturbed parameters, denoted as SEO_PP, is performed, namely, the 9 parameters (σ , κ, b, c 1 , c 2 , c 3 , c 4 , c 5 , and c 6 ) are perturbed at the same time when the ensemble initial conditions are formed by adding a Gaussian white noise with the standard deviation being 5 % of the default values.Third, parameter estimation is performed to optimize all these 9 parameters using the DAEPC algorithm, denoted as PP_PO.It should be noted that the other 6 parameters (O m , O d , S m , S s , S pd , and ) act as the part of the dynamic core that determines the timescale and period of the outer forcing of this coupled system.Therefore, they will not alter once they are determined using the default values in this study.In addition, a control run without any observational constraint, called CTRL, serves as a reference for the evaluation of the assimilation results.

Impact of parameter optimization on climate estimation
Figure 4 shows the time series of root mean square errors (RMSEs) of x 2 in the last 100 TUs (Fig. 4a), w in the last 10 3 TUs (Fig. 4b), and η in 10 4 TUs (Fig. 4c) obtained from CTRL (pink line), SEO (black line), SEO_PP (red line), and PP_PO (blue line).Here, all the RMSE time series are computed from the difference between the ensemble mean of the model run and the truth at each time integral step over a 1-TU time window.From the black line in Fig. 4a, we can see that the estimate of x 2 in SEO does not show obvious improvement compared with that in CTRL (pink line).Meanwhile, SEO has little contribution to the estimate of the slow-varying variables such as w or η (black line vs. pink line in Fig. 4b and c).The time mean RMSEs of x 2 , w, and η during the entire data assimilation period in SEO are actually reduced, with respect to CTRL, by 2, 9, and 12 %, respectively.The improvements are not much and can hardly be observed in the time series of RMSEs.This means that within the framework of the specific dynamical core biased twin experiment setting with this simple model, the dynamical core misfitting of the assimilation model with respect to the truth model can function as a significant obstacle for the traditional SEO.In what follows, we will see the significant improvement with parameter optimization on assimilation quality, which addresses the potential role of parameter optimization in the mitigation of dynamical core model bias.Compared with SEO, estimates of all three variables (x 2 , w, and η) in SEO_PP demonstrate significant improve-ment (red lines in Fig. 4a, b, and c), evidenced by one-order RMSE magnitude reduction with a small oscillation.This suggests that the perturbed physical parameters can improve the quality of state estimation greatly, since the perturbation of parameters can increase the spread of model states so as to increase the representation of model ensemble for the model uncertainty.The RMSE of x 2 in PP_PO is reduced with the same magnitude as in SEO_PP (blue line vs. red line in the small panel in Fig. 4a).That is to say, the improvement in state estimation via parameter optimization is similar to that via the initial perturbation of parameters for the high-frequency atmospheric variables.However, this is not the case for the low-frequency oceanic states.The RMSE of w in PP_PO (blue line in the small panel in Fig. 4b) is reduced by about 66 % compared to that in SEO_PP (red line).So, the RMSE is further reduced by observation-optimized model parameters via PP_PO relative to SEO_PP.For the deep ocean, η in PP_PO is further improved significantly from SEO_PP.From Fig. 4c, it can be seen that the blue line and the red line overlap completely before 10 3 TUs during the early state estimation only period.When the process of parameter optimization starts after 10 3 TUs, the RMSE of η in PP_PO decreases from about 0.68 to about 0.2, while the RMSE in SEO_PP remains at a level of 0.68.This means that observational information can be effectively extracted by the process of parameter optimization for both the upper ocean and the deep ocean.
Figure 5 shows the time series of the ensemble means of κ (Fig. 5a) and c 2 (Fig. 5b) between 900 and 1200 TUs in PP_PO.We can see that the ensemble means of the parameters remain at their default values before 10 3 TUs in which parameter optimization is not activated yet, but oscillate after 10 3 TUs to compensate for the bias of state estimation induced by the biased time-differencing through the background error covariance between parameters and observations.The oscillation frequency of κ (see Fig. 5a) in the atmospheric control equation is higher than that of c 2 (see Fig. 5b) in the oceanic equation, showing that the oscillation frequencies of parameters are affected strongly by the frequencies of state variables.As noted above, model parameters do not have any dynamically supported internal variability, and their variation depends completely on the variation of the state-parameter covariance.
Figure 6 shows the power spectra of w (Fig. 6a) and η (Fig. 6b) derived from the truth (red line), SEO (black line), and PP_PO (blue line), respectively.The significant characteristic timescales of w and η in PP_PO are almost the same as those in the truth (blue line vs. red line), while the significant characteristic timescales of w and η in SEO are distinctly different from the truth (black line vs. red line).This means that the characteristics of climate variability in the upper and deep oceans can be reconstructed accurately by ensemble coupled parameter optimization when biased timedifferencing is the main source of model bias.

Impact of parameter estimation in different media on model bias mitigation
In this section, we investigate the impact of different media parameter estimations on model bias mitigation.A series of assimilation experiments, named Case1, Case2, Case3 till Case11, as listed in Table 1, is designed to accomplish this objective.
Figure 7 presents the time series of RMSEs of x 2 in the last 100 TUs (Fig. 7a), w in the last 100 TUs (Fig. 7b), and η in 10 4 TUs (Fig. 7c) in Case1 (black line), Case2 (green line), Case3 (pink line), and Case4 (cyan line), respectively.For x 2 in the transient atmosphere, Case2 is the worst (green line in Fig. 7a) among all the four cases presented.The reason is that c 1 and c 2 being optimized in Case2 are the coupling parameters between the atmosphere variable x 2 and the upper ocean variable w.Compared with the other three cases, optimizing these two coupling parameters alone cannot efficiently enhance the signal-to-noise ratio of the background error covariance between the parameter and the model state.For the upper ocean, it is noticed from Fig. 7b that the RM-SEs of w in both Case1 and Case4 remain similarly larger values than the other two cases.This can be understood as follows.Case1 optimizes parameters in the atmosphere (σ , κ, and b) and Case4 optimizes parameters in the deep ocean (c 5 and c 6 ), respectively.Neither one has a direct effect on the upper ocean variable w, and therefore neither one is able to affect w significantly.
For η, the deep ocean component, compared to SEO_PP (red line in Fig. 7c), although the RMSEs can be reduced  in Case1 (black line in Fig. 7c) and Case2 (green line in Fig. 7c) to some degrees, both of which remain unstable with a large oscillation.That is because observational information is used only to adjust the parameters pertinent to the highfrequency atmospheric variables instead of adjusting those pertinent to the low-frequency oceanic variables.In contrast, Case3 has the best estimate of η (pink line in Fig. 7c), with an RMSE reduction over 75 % compared to the case in which all parameters are optimized in PP_PO (blue line in Fig. 7c).This means that some parameters being optimized introduce noises under given conditions.Case3 optimizes c 3 and c 4 , which represent the linear and the nonlinear interactions between the upper and deep oceans.The observational information of w can be retrieved efficiently to improve the The setups of the four experiments are described in Table 1.Note that the time series of RMSEs in SEO_PP (red line; only in c) and PP_PO (blue line) are also plotted as references.
estimation of w and η significantly through the parameter optimization and the coupling between the upper ocean and pycnocline depth.It is noticed that the performance of Case4 (cyan line in Fig. 7c) is the worst among all four cases, which has an even larger RMSE than SEO_PP (red line in Fig. 7c).Parameters c 5 and c 6 being optimized in Case4 are in the η equation.Due to the absence of the observation of η, c 5 and c 6 are estimated only indirectly through the observation of w.This suggests that only optimizing parameters pertinent to the deep ocean component may lead to the noisedominating background error covariance when other parameters retain their default values.To explore this issue further, we performed two additional experiments, Case5 and Case6 (Table 1).Case5 optimizes c 5 and c 6 together with c 3 and c 4 .Case6 optimizes all parameters but excluding c 5 and c 6 .It is already learned from the discussions before that the case with only c 3 and c 4 being optimized (i.e., Case3) has the best performance.However, the additional inclusion of c 5 and c 6 in Case5 deteriorates the results of η (see green line in Fig. 8), yielding larger RMSEs than the case in which all parameters are optimized in PP_PO (blue line in Fig. 8).In contrast, the RMSEs of η in Case6 (black line in Fig. 8) have a similar mitigation to that in Case3, indicating that the exclusion of c 5 and c 6 increases the signal-to-noise ratio during parameter optimization.
The poor performance of additional optimization of c 5 and c 6 may be associated with the insufficient ensemble representation and the lack of the η observation.To verify the first reason, we increase the ensemble size of SEO_PP, PP_PO, and Case4 to 200, and denote them as Case7, Case8, and Case9 (see Table 1).The time-averaged RMSE of η in Case7 is reduced by more than 30 % from SEO_PP, which can be seen from the green and red lines in Fig. 9.This means that the capability of ensemble spread is enhanced with a large ensemble size.However, the RMSEs of η in Case8 are not smaller than ones in PP_PO (black line vs. blue line in Fig. 9), which indicates that increasing the ensemble size does not improve the η estimate due to the existence of noises during signalenhanced parameter optimization, even if the ensemble representation can be improved by increasing the ensemble size.The RMSEs of η in Case9 are slightly larger than those in Case7 (cyan line vs. green line in Fig. 9).Thus, the noises induced by optimizing c 5 and c 6 cannot be removed successfully through increasing the ensemble size.However, it is clear that the performance of Case4 is indeed improved greatly when the ensemble representation is enhanced in Case9 (comparing both cyan lines in Figs.7c and 9).
Regarding the second reason, to demonstrate the effect of the assimilation of η observations, we performed two additional experiments, Case10 and Case11 (see Table 1).Case10 and Case11 are the counterpart of SEO_PP and Case4 but with the assimilation of η observations that have the same sampling frequency as w.Compared to SEO_PP, the RMSEs  of η in Case10 (green line in Fig. 10) are reduced significantly, but still a little larger than those in PP_PO (blue line in Fig. 10).This suggests that the signal-enhanced parameter optimization can improve the η estimation even in the absence of observations of the deep ocean.The benefit of assimilating deep ocean observations is significant when c 5 and c 6 are optimized (Case11).The RMSEs of η in Case11 (black line in Fig. 10) remain well consistent with those in Case3, indicating that the deep ocean observations may be indispensable when parameters pertinent to low-frequency components are optimized through signal-enhanced parameter estimation.Otherwise, noises of the background error covariance will be enlarged to deteriorate the quality of climate estimation, especially for the low-frequency signals.

Summary and discussions
A biased twin experiment framework is designed to study the mitigation of coupled model bias induced by imperfect time-differencing through parameter optimization.In this framework, the assimilation model includes a "biased" time-differencing scheme from the truth model that is used to produce "observations".The leap-frog time-differencing scheme with a Robert-Asselin time filter serves as the truth run from which "observations" are drawn.The fourth-order Runge-Kutta time-differencing scheme is used for the assimilation runs.A series of assimilation experiments is performed to examine the impact of parameter optimization on climate estimation and model bias mitigation, as well as the role of different media parameter estimations.Results show that the initial perturbations of parameters can enhance the ensemble spread greatly and improve the representation of the model ensemble to the model error that consists of biases arising from different differencing schemes.Furthermore, parameter estimation enhances the accuracy of climate estimation, especially for low-frequency signals.In addition, in a multiple timescale coupled system, parameters pertinent to low-frequency components have more impact on climate signals.
Although the parameter optimization shows promise with the simple coupled model to mitigate the model bias arising from dynamical core misfitting, further research is required to understand the impact of such model biases in coupled general circulation models (CGCMs) on climate estimation and prediction.It should be kept in mind that the concept model used in this study contains a rather large set of parameters compared to the dimension of the model attractor.This will not be the case for the state-of-the-art CGCMs in general.There might also not be enough data to constrain both the model states and additional parameters under the modern atmospheric and oceanic observing networks.Thus potential benefits of parameter estimation vs stochastic physics, achieved in this study, need to be investigated further in the realistic scenario.Nevertheless, imperfect numerical schemes are usually used in CGCMs.Therefore, besides the development of more robust numerical schemes, how to optimize model parameters needs to be examined further to compensate for the deficiencies of the currently used numerical schemes.

Fig. 1 .
Fig. 1.Time series of (a) x 2 in the first 100 TUs, (b) w in the first 10 3 TUs, and (c) η in 10 4 TUs derived from both the LF (red line) and the RK4 (black line) time-differencing schemes.

Fig. 2 .
Fig. 2. Variation of (a) x 2 in the space of x 3 and (b) η in the space of w derived from both the LF (red circle) and the RK4 (black circle) time-differencing schemes.

Fig. 3 .
Fig. 3. Power spectra of (a) w and (b) η based on the model data between 10 3 and 10 4 TUs derived from both the LF (red line) and the RK4 (black line) time-differencing schemes.

Fig. 4 .
Fig. 4. Time series of RMSEs of (a) x 2 in the last 100 TUs, (b) w in the last 10 3 TUs, and (c) η in 10 4 TUs in CTRL (pink line), SEO (black line), SEO_PP (red line), and PP_PO (blue line).Small panels in (a) and (b) show the detailed evolutions of RMSEs of x 2 and w in SEO_PP (red line) and PP_PO (blue line) through enlarging the scale of the vertical coordinate.

Fig. 7 .
Fig. 7. Time series of RMSEs of (a) x 2 in the last 100 TUs, (b) w in the last 100 TUs, and (c) η in 10 4 TUs in Case1 (black line), Case2 (green line), Case3 (pink line), and Case4 (cyan line), respectively.The setups of the four experiments are described in Table1.Note that the time series of RMSEs in SEO_PP (red line; only in c) and PP_PO (blue line) are also plotted as references.

Fig. 8 .
Fig. 8. Time series of RMSEs of η in 10 4 TUs in Case5 (green line) and Case6 (black line).The setups of the two experiments are described in Table 1.Note that the time series of RMSEs of η in SEO_PP (red line), PP_PO (blue line), and Case3 (pink line) are also plotted as references.

Fig. 9 .
Fig. 9. Time series of RMSEs of η in 10 4 TUs in Case7 (green line), Case8 (black line), and Case9 (cyan line).The setups of the three experiments are described in Table1.Note that the time series of RMSEs of η in SEO_PP (red line) and PP_PO (blue line) are also plotted as references.

Fig. 10 .
Fig. 10.Time series of RMSEs of η in 10 4 TUs in Case10 (green line) and Case11 (black line).The setups of the two experiments are described in Table 1.Note that the time series of RMSEs of η in SEO_PP (red line), PP_PO (blue line), and Case3 (pink line) are also plotted as references.